 |
Here is the Maple code for the error functions:
>
#This procedure creates a function in out_func which is u-Pu
#It also creates a function in out_func + 1 which is (u-Pu)^2 (used for the L2_norm)
error_fns:=proc(T,m,func_T,func_A,out_func)
local i,vm;
vm:=vertices(m);
for i from 1 to nops(vm) do
T[vm[i],out_func] := T[vm[i],func_T] - T[vm[i],func_A];
T[vm[i],out_func+1]:= (T[vm[i],func_T] - T[vm[i],func_A])^2;
od;
end;
>
#This function takes the L2_norn of error_func, assuming that the square of error_func
#is in error_func+1, as it would be if generated by error_fns.
L2_norm:=proc(T,m,error_func) #!!!!NOTE, ASSUMES THE SQUARE OF THE FUNCTION IS IN ERROR_FUNC+1!!
local err,i;
err:=simp_rule(T,m,error_func+1);
sqrt(err);
end;
>
#This procedure finds the square root of the energy of func_in.
#Note that this function requires that make_neighbors has been run
#on T out to level m...
get_root_energy:=proc(T,m,func_in)
local err, Vm,V0,j, i;
err:=0;
Vm:=vertices(m);
V0:=vertices(0);
for i from 4 to nops(Vm) do
err := err + (5/3)^m * (add( (T[T[Vm[i],j],func_in] - T[Vm[i],func_in])^2, j = 1..4))/2;
od:
for i from 1 to 3 do
err := err + (5/3)^m * (add( (T[T[Vm[i],j],func_in] - T[Vm[i],func_in])^2, j = 1..2))/2;
od;
sqrt(err);
end;
>
#This procedure finds the maximum of a function (used to find max error...)
find_max:=proc(T,m,func_in)
local i,abs_diff,temp,Vm;
Vm:=vertices(m);
temp:=0;
for i from 1 to nops(Vm) do
abs_diff := abs(T[Vm[i],func_in]):
if abs_diff > temp then temp := abs_diff; fi; od:
temp;
end;
>
#So this function estimates the posteriori error with q nonconstant (if q is constant, just stick in a
#function with constant values). It calculates L2_norm(-Lp+qp-f). This function requires you send in C,
#C being the coefficients of the spline representation of p. This allows us to calculate a theoretical
#laplacian, since we know the exact values of the laplacian of the basis functions (they're biharmonic...).
#Also, l_funcno is just some random function number where you want the laplacian of the function stored.
#Note that the laplacian stored in l_funcno has incorrect values at the vertices of the spline basis, because
#the laplacian is discontinuous at those points. The solution to this problem is documented within this procedure.
posteriori_error_q_nonc:=proc(T,S,m_basis,m_func,C,p_funcno,q_funcno,f_funcno,harm_funcno,l_funcno)
local vm,Vm,fact,V0,Vother,Vother_e,i,basis,triangle1,triangle2,addfuncs,addresses1,addresses2,tris,total_int,Vmsum1,Vothersum1,V0sum1,j;
vm:=vertices(m_basis);
fact:=1/(6*3^m_func); #This is just the factor from simpson's rule.
#Here, we are just getting the vertices of V0, Vm, and (in Vother) everything in between.
V0:=vertices_m(0);
if m_basis = 1 then Vother:=[] else Vother_e:=vertices(m_func-m_basis-1);Vother:=Vother_e[4..nops(Vother_e)]; fi;
Vm:=vertices_m(m_func-m_basis);
#Here, we are getting all the triangles in our spline basis and storing them in tris
if m_basis=1 then tris:=[[0],[1],[2]] else tris:=p([0,1,2],m_basis) fi;
#Okay, so here's the general idea for dealing with the discontinuous laplacian:
#We will store the addresses of the vertices of the triangles in tris in S[tris[i],10-12],
#which is performed by get_triangle_adds. Then, in S[tris[i],7-9], we store the values of the
#laplacians at the boundaries of tris[i]. Thus, we have stored both the values of the laplacian
#at each spline basis vertex: once for each triangle. Then we will piecewise integrate our function
#by summing the integrals over all the tris[i], using the boundary values which correspond to that
#particular triangle. Pretty tricky, huh?
for i from 1 to nops(tris) do S[tris[i],7]:=0; S[tris[i],8]:=0; S[tris[i],9]:=0; od; #initialize values to 0...
get_triangle_adds(S,m_basis); #get's the addresses of the boundaries of tris and stores them as indicated
#Here, we generate the laplacian for all values EXCEPT for the spline vertices
print(`making laplacian...`);
make_pretty_picture_lap(T,S,m_basis,m_func,vm,l_funcno,harm_funcno,C);
print(`finished laplacian...`);
for i from 1 to 3^(m_basis+1) do #First, let's put the values of the laplacians at the spline basis vertices (DISCONTINUOUS) in the table S...
basis:=get_basis(vm,i);
if basis[1]<>[0] and basis[1]<>[1] and basis[1]<>[2] then #ie, not a boundary point
triangle1:=get_triangle(S,basis,m_basis,1);
triangle2:=get_triangle(S,basis,m_basis,2);
#Here, we determine which basis functions we need, which again
#depend entirely on the last two digits of the basis elt (ie, the orientation)
if last_digits(basis[1])=[0,1] then addfuncs:=[1,5,3,7] fi;
if last_digits(basis[1])=[1,2] then addfuncs:=[5,9,7,11] fi;
if last_digits(basis[1])=[0,2] then addfuncs:=[9,1,11,3] fi;
#And then just start adding the laplacian values in.
#Notice that the +12 after harm_funcno means that we are
#accessing the laplacians of the basis functions, rather than the
#basis functions themselves.
if basis[2]=f then
S[triangle1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],harm_funcno+12+addfuncs[1]]*5^m_basis;
S[triangle1,7+1]:=S[triangle1,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs[1]]*5^m_basis;
S[triangle1,7+2]:=S[triangle1,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs[1]]*5^m_basis;
S[triangle2,7+0]:=S[triangle2,7+0]+C[i]*T[[0],harm_funcno+12+addfuncs[2]]*5^m_basis;
S[triangle2,7+1]:=S[triangle2,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs[2]]*5^m_basis;
S[triangle2,7+2]:=S[triangle2,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs[2]]*5^m_basis;
else #DON'T FORGET THE NEGATIVE SIGN FOR G-TYPE BASIS ELTS!!!
S[triangle1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],harm_funcno+12+addfuncs[3]]*5^m_basis;
S[triangle1,7+1]:=S[triangle1,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs[3]]*5^m_basis;
S[triangle1,7+2]:=S[triangle1,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs[3]]*5^m_basis;
S[triangle2,7+0]:=S[triangle2,7+0]-C[i]*T[[0],harm_funcno+12+addfuncs[4]]*5^m_basis;
S[triangle2,7+1]:=S[triangle2,7+1]-C[i]*T[[1],harm_funcno+12+addfuncs[4]]*5^m_basis;
S[triangle2,7+2]:=S[triangle2,7+2]-C[i]*T[[2],harm_funcno+12+addfuncs[4]]*5^m_basis;
fi;
else #If we have a basis elt which is on the boundary
triangle1:=[seq(op(basis[1]),j=1..m_basis)]; #So the triangle we're dealing with is just a sequence of
#0's, 1's, or 2's, depending on the value of the boundary
if basis[1]=[0] then addfuncs:=3 fi;
if basis[1]=[1] then addfuncs:=7 fi;
if basis[1]=[2] then addfuncs:=11 fi;
#We know it's a g-type basis elt, and that the values will always be positive...
S[triangle1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],harm_funcno+12+addfuncs]*5^m_basis;
S[triangle1,7+1]:=S[triangle1,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs]*5^m_basis;
S[triangle1,7+2]:=S[triangle1,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs]*5^m_basis;
fi;
od;
print(`integrating difference fn... `);
#So now we're ready to integrate each triangle...
total_int:=0;
for j from 1 to nops(tris) do
printf(`%d `,j);
if ((j mod 15)=1) then
printf(`\n`): fi:
#So basically, we call func_val to get the value of the at the vertex defined by V0 (Vother,Vm) [i], when considered as
#a vertex of the subtriangle tris[j]. This basically encapsulates the functionality of coord_transfer. So we just
#add up -Lp+qp-f. Notice that in the V0sum1, instead of using -func_val(laplacian), we use -S[tris[j],7+i-1], which is
#the value of the laplacian at the boundary when considering the triangle tris[j]. Thus, we use the value of the
#laplacian appropriate to each triangle, thereby dealing with the multiple values of the laplacian.
V0sum1:=add( ( -S[tris[j],7+i-1]+func_val(T,S,tris[j],V0[i],q_funcno)*func_val(T,S,tris[j],V0[i],p_funcno)-func_val(T,S,tris[j],V0[i],f_funcno) )^2 ,i=1..nops(V0) );
#print(`sum1`);
if m_basis <> 1 then Vothersum1:=add( ( -func_val(T,S,tris[j],Vother[i],l_funcno)+func_val(T,S,tris[j],Vother[i],q_funcno)*func_val(T,S,tris[j],Vother[i],p_funcno)-func_val(T,S,tris[j],Vother[i],f_funcno) )^2,i=1..nops(Vother)) else Vothersum1:=0 fi;
#print(`sum2`);
Vmsum1:=add( ( -func_val(T,S,tris[j],Vm[i],l_funcno)+func_val(T,S,tris[j],Vm[i],q_funcno)*func_val(T,S,tris[j],Vm[i],p_funcno)-func_val(T,S,tris[j],Vm[i],f_funcno) )^2,i=1..nops(Vm));
#print(`sum3`);
total_int:=total_int + fact*(5*Vmsum1 + 2*Vothersum1 + V0sum1);
od;
sqrt(total_int); #Don't forget the square root!
end;
>
#This function takes the value of funcno at address add when add is considered an address on the subtriangle
#triangle. For instance, if triangle is 0 (ie, the top 1/3 of the SG), then func_val(T,S,triangle,[0,1],funcno)
#will return the value of funcno at 0,0,1, and if add=[1], then we get the value of funcno at 0,1. Notice
#that this requires that we have run get_triangle_adds, since coord_transfer requires the addresses of the vertices
#of the triangles to operate...
func_val:=proc(T,S,triangle,add,funcno)
T[coord_transfer(add,triangle,S[triangle,10],S[triangle,11],S[triangle,12]),funcno] end;
>
>
#So this function stores the addresses of the boundaries of the triangles
#out to m_basis in S[triaddress,10-12], which 10,11,12 corresponding to add0,add1,add2 respectively.
#It does V1 manually, and then calls get_triangle_addresses, which recursively does the rest.
get_triangle_adds:=proc(S,m_basis)
if m_basis=1 then
S[[0],10]:=[0];
S[[0],11]:=[0,1];
S[[0],12]:=[0,2];
S[[1],10]:=[0,1];
S[[1],11]:=[1];
S[[1],12]:=[1,2];
S[[2],10]:=[0,2];
S[[2],11]:=[1,2];
S[[2],12]:=[2];
else
get_triangle_addresses(S,[0],[0],[0,1],[0,2],m_basis-1);
get_triangle_addresses(S,[1],[0,1],[1],[1,2],m_basis-1);
get_triangle_addresses(S,[2],[0,2],[1,2],[2],m_basis-1);
fi;
end;
>
#This is the recursive call which get_triangle_adds calls to generate the triangle_addresses for m_basis>1...
get_triangle_addresses:=proc(S,triangle,add0,add1,add2,m_basis)
#This is our base case. We simply assign the values according to triangle
if m_basis=1 then
S[[op(triangle),0],10]:=add0;
S[[op(triangle),0],11]:=[op(triangle),0,1];
S[[op(triangle),0],12]:=[op(triangle),0,2];
S[[op(triangle),1],10]:=[op(triangle),0,1];
S[[op(triangle),1],11]:=add1;
S[[op(triangle),1],12]:=[op(triangle),1,2];
S[[op(triangle),2],10]:=[op(triangle),0,2];
S[[op(triangle),2],11]:=[op(triangle),1,2];
S[[op(triangle),2],12]:=add2;
else
#if we aren't in the base case, we simply call the procedure again 3 times,
#except that this time, we add on the appropriate extra digit onto triangle.
#Notice that nothing gets assigned until we reach the base case; the rest is
#just setting it up (primarily just figuring out which triangles to assign).
get_triangle_addresses(S,[op(triangle),0],add0,[op(triangle),0,1],[op(triangle),0,2],m_basis-1);
get_triangle_addresses(S,[op(triangle),1],[op(triangle),0,1],add1,[op(triangle),1,2],m_basis-1);
get_triangle_addresses(S,[op(triangle),2],[op(triangle),0,2],[op(triangle),1,2],add2,m_basis-1);
fi;
end;
>
#I don't know if this function is ever used, but it returns the triangle
#address of triangle of which add is a midpoint.
get_tri_digits:=proc(add)
if nops(add)=1 or nops(add) = 2 then [] else add[1..nops(add)-2] fi; end;
|
