 |
Here is the Maple code for working out the FEM using Harmonic splines:
>
#This function creates the basis functions for the harmonic
#finite element method. The basis functions are quite simple:
#They are just 1 at one boundary, 0 at the other boundaries,
#and filled in harmonically the rest of the way. Here, I store
#(1,0,0) in funcno, (0,1,0) in funcno+1, and (0,0,1) in funcno+2,
#where the first number corresponds to the value at [0], the second num
#to the value at [1], and I think you can figure out the last case.
#It should be noted that funcno, funcno+1, and funcno+2 are all used up
#by this function.
calc_harms:=proc(T,m,funcno)
T[[0],funcno]:=1;
T[[1],funcno]:=0;
T[[2],funcno]:=0;
T[[0,1],funcno]:=2/5;
T[[1,2],funcno]:=1/5;
T[[0,2],funcno]:=2/5;
harm(T,m,[0],[0,1],[0,2],funcno);
harm(T,m,[1],[1,2],[0,1],funcno);
harm(T,m,[2],[0,2],[1,2],funcno);
T[[0],funcno+1]:=0;
T[[1],funcno+1]:=1;
T[[2],funcno+1]:=0;
T[[0,1],funcno+1]:=2/5;
T[[1,2],funcno+1]:=2/5;
T[[0,2],funcno+1]:=1/5;
harm(T,m,[0],[0,1],[0,2],funcno+1);
harm(T,m,[1],[1,2],[0,1],funcno+1);
harm(T,m,[2],[0,2],[1,2],funcno+1);
T[[0],funcno+2]:=0;
T[[1],funcno+2]:=0;
T[[2],funcno+2]:=1;
T[[0,1],funcno+2]:=1/5;
T[[1,2],funcno+2]:=2/5;
T[[0,2],funcno+2]:=2/5;
harm(T,m,[0],[0,1],[0,2],funcno+2);
harm(T,m,[1],[1,2],[0,1],funcno+2);
harm(T,m,[2],[0,2],[1,2],funcno+2);
end;
>
>
>
#This function computes the inner product of two harmonic basis elements.
#It uses theoretical values for inner products of harmonic functions. It
#then scales down the answer by a factor of 3^m, because the triangles are
#smaller by a factor of 3^m.
inner_prod_h:=proc(S,m,basis_1,basis_2)
local out,temp,basis1,basis2;
basis1:=basis_1;
basis2:=basis_2;
#if they aren't neighbors, we get nothing...
if basis1[1] <> basis2[1] and S[basis1[1],1] <> basis2[1] and S[basis1[1],2] <> basis2[1] and S[basis1[1],3] <> basis2[1] and S[basis1[1],4] <> basis2[1] then RETURN(0) fi;
#If the 2 basis elts are the same, then we get 7/45*2 (the two comes from the 2 neighboring triangles)
#the value 7/45 is theoretical
if basis1[1]=basis2[1] then
out:=7/45 * 2 / (3^m);
fi;
#If the 2 basis elts are neighbors, then we get 4/45 (theoretical value)
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[basis1[1],3] = basis2[1] or
S[basis1[1],4] = basis2[1]
then
out:=4/45 / (3^m);
fi;
out;
end;
>
#This function returns the energy of two basis functions.
#It works just like inner_prod_h, except that it uses a scaling factor of
#(5/3)^m, which is the energy scaling factor from kigami.
energy_h:=proc(S,m,basis_1,basis_2)
local out,temp,basis1,basis2;
basis1:=basis_1;
basis2:=basis_2;
#if they aren't neighbors, we get nothing...
if basis1[1] <> basis2[1] and S[basis1[1],1] <> basis2[1] and S[basis1[1],2] <> basis2[1] and S[basis1[1],3] <> basis2[1] and S[basis1[1],4] <> basis2[1] then RETURN(0) fi;
if basis1[1]=basis2[1] then
out:=2 * 2 * (5/3)^m;
fi;
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[basis1[1],3] = basis2[1] or
S[basis1[1],4] = basis2[1]
then
out:=-1 * (5/3)^m;
fi;
out;
end;
>
#This function takes a basis element and a function, and returns integral(basis*function).
#The nice thing is that the basis element is zero everywhere except for 2 triangles. Thus,
#we can speed computations by just integrating over those particular triangles.
#Note that the integration method used here is simpson's rule out to maximum depth in the SG.
int_vs_basis_h:=proc(T,S,m_basis,m_func,basis,funct,harm_funcno)
local V0,Vother,Vm,fact,triangle1,triangle2,V0sum1,V0sum2,Vothersum1,Vothersum2,Vmsum1,Vmsum2,i,j,Vother_e;
fact:= 1/(6*3^m_func); #This is the factor for simpson's rule...
#get the vertices for V0, Vm, and everything in between...
#Notice that here we are taking m to be m_func-m_basis, because what we want to do is
#consider vertices scaled down by m_basis in funct and multiply by vertices in the basis
V0:=vertices_m(0);
if m_func-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);
#all basis elements are of type f...
if basis[2]=f then
#So basically, we will first check the orientation of the basis element (i.e, the last two digits).
#Then, we get the triangles associated with the basis element. Then, given a point in Vm_func-m_basis,
#we can use coord_transfer to determine (using the triangle) which is the corresponding coordinate
#in Vm_func, and then multiply that by the corresponding value of the basis function, and sum.
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1] then
triangle1:=get_triangle(S,basis,m_basis,1);
triangle2:=get_triangle(S,basis,m_basis,2);
#let's calculate sum's for triangle 1, shall we?
V0sum1:=add(T[coord_transfer(V0[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funct] *(T[V0[i],harm_funcno]),i=1..nops(V0));
Vothersum1:=add(T[coord_transfer(Vother[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funct] *(T[Vother[i],harm_funcno]),i=1..nops(Vother));
Vmsum1:=add(T[coord_transfer(Vm[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funct] *(T[Vm[i],harm_funcno]),i=1..nops(Vm));
#and now on to triangle numero 2!
V0sum2:=add(T[coord_transfer(V0[i],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funct]*(T[V0[i],harm_funcno+1] ),i=1..nops(V0));
Vothersum2:=add(T[coord_transfer(Vother[i],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funct] *(T[Vother[i],harm_funcno+1] ),i=1..nops(Vother));
Vmsum2:=add(T[coord_transfer(Vm[i],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funct]*(T[Vm[i],harm_funcno+1]),i=1..nops(Vm));
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then
triangle1:=get_triangle(S,basis,m_basis,1);
triangle2:=get_triangle(S,basis,m_basis,2);
#let's calculate sum's for triangle 1, shall we?
V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funct]*(T[V0[i],harm_funcno+2] ),i=1..nops(V0));
Vothersum1:=add(T[coord_transfer(Vother[i],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funct] *(T[Vother[i],harm_funcno+2]),i=1..nops(Vother));
Vmsum1:=add(T[coord_transfer(Vm[i],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funct]*(T[Vm[i],harm_funcno+2]),i=1..nops(Vm));
#and now on to triangle numero 2!
V0sum2:=add(T[coord_transfer(V0[i],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funct]*(T[V0[i],harm_funcno] ),i=1..nops(V0));
Vothersum2:=add(T[coord_transfer(Vother[i],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funct] *(T[Vother[i],harm_funcno]),i=1..nops(Vother));
Vmsum2:=add(T[coord_transfer(Vm[i],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funct]*(T[Vm[i],harm_funcno] ),i=1..nops(Vm));
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then
triangle1:=get_triangle(S,basis,m_basis,1);
triangle2:=get_triangle(S,basis,m_basis,2);
#let's calculate sum's for triangle 1, shall we?
V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*(T[V0[i],harm_funcno+1]),i=1..nops(V0));
Vothersum1:=add(T[coord_transfer(Vother[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct] *(T[Vother[i],harm_funcno+1]),i=1..nops(Vother));
Vmsum1:=add(T[coord_transfer(Vm[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*(T[Vm[i],harm_funcno+1]),i=1..nops(Vm));
#and now on to triangle numero 2!
V0sum2:=add(T[coord_transfer(V0[i],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funct]*(T[V0[i],harm_funcno+2]),i=1..nops(V0));
Vothersum2:=add(T[coord_transfer(Vother[i],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funct] *(T[Vother[i],harm_funcno+2]),i=1..nops(Vother));
Vmsum2:=add(T[coord_transfer(Vm[i],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funct]*(T[Vm[i],harm_funcno+2]),i=1..nops(Vm));
fi;
fi;
#And then just return the formula for simpson's rule...
RETURN(fact * (5*(Vmsum1+Vmsum2) + 2*(Vothersum1+Vothersum2) + (V0sum1+V0sum2)));
end;
>
#Here, we generate the matrix of E(basisi,basisj) out to level m_basis
#The final matrix is a sparse, symmetric matrix, and we take that into
#account here.
gen_energy_matrix_h:=proc(S,m_basis,vm)
local i,j,en_matrix,n_verts,en;
n_verts:=nops(vm)-3; sum(3^i,i=1..m_basis);
en_matrix:=array(sparse,1..n_verts,1..n_verts);
for i from 1 to n_verts do
for j from 1 to n_verts do
en:=energy_h(S,m_basis,get_basis(vm,i),get_basis(vm,j));
if en <> 0 then en_matrix[i,j]:=en; en_matrix[j,i]:=en; fi;
od;
od;
en_matrix;
end;
>
#This function returns a basis element given any particular basis element number.
#Notice that this function works for both harmonic and biharmonic splines; if the
#number is greater than the number of f-type vectors, we return a g-type vector,
#but if we are using harmonic splines, we will never get a number that high, and
#so all we get are the f-type vectors.
get_basis:=proc(vm,number)
if number <= nops(vm)-3 then [vm[number+3],f]
else [vm[number-nops(vm)+3],g] fi; end;
>
#This is just like gen_energy_matrix_h, except that creates the inner product matrix.
gen_iprod_matrix_h:=proc(S,m_basis,vm)
local i,j,i_matrix,n_verts,ip;
n_verts:=nops(vm)-3; add(3^i,i=1..m_basis);
i_matrix:=array(sparse,1..n_verts,1..n_verts);
for i from 1 to n_verts do
for j from 1 to n_verts do
ip:=inner_prod_h(S,m_basis,get_basis(vm,i),get_basis(vm,j));
if ip <> 0 then i_matrix[i,j]:=ip; i_matrix[j,i]:=ip; fi;
od;
od;
i_matrix;
end;
>
#This procedure generates a vector of integral(funcno*basisi).
#It calls int_vs_basis_h, where most of the hard work is actually
#done.
gen_f_vector_h:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno)
local i,f_vec,n_verts;
n_verts:=add(3^i,i=1..m_basis);
f_vec:=array(1..n_verts);
for i from 1 to n_verts do
f_vec[i]:=int_vs_basis_h(T,S,m_basis,m_func,get_basis(vm,i),funcno,harm_funcno);
od;
f_vec;
end;
>
#This function is called by make_pretty_picture_h. It basically goes through not assigning values to
#any vertices until it reaches the end of m_basis, and then it just fills everything in harmonically,
#basically making a piecewise harmonic function with boundary values at the harmonic spline basis
#vertices.
make_h_answer:=proc(T,m_basis,m_func,add0,add1,add2,funcno);
if m_basis = 0 then
# print(add0,add1,add2);
harm(T,m_func,add0,add1,add2,funcno);
else
if add1[nops(add1)-1..nops(add1)] = [0,1] then
make_h_answer(T,m_basis-1,m_func-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno);
make_h_answer(T,m_basis-1,m_func-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1],funcno);
make_h_answer(T,m_basis-1,m_func-1,add2,[op(add1[1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [1,2] then
make_h_answer(T,m_basis-1,m_func-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno);
make_h_answer(T,m_basis-1,m_func-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],funcno);
make_h_answer(T,m_basis-1,m_func-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [0,2] then
make_h_answer(T,m_basis-1,m_func-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2],funcno);
make_h_answer(T,m_basis-1,m_func-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno);
make_h_answer(T,m_basis-1,m_func-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1],funcno);
fi;
fi;
end;
>
#This function basically copies the values of C (a projective representation of a function
#in the harmonic spline basis) into a function, and then calls make_h_answer to fill in the
#rest of the function. This procedure is far more efficient than the method used in
#make_pretty_picture because there, we have to add in each basis function manually.
#Here, we can use the fact that the harmonic functions are completely defined
#by knowing the values of the function on the boundary and then just using the harmonic
#algorithm to determine the rest.
make_pretty_picture_h:=proc(T,m_basis,m_func,C,funcno)
local vm,j;
vm:=vertices(m_basis);
#first, set the initial values...
T[[0],funcno]:=0 ;T[[1],funcno]:=0; T[[2],funcno]:=0;
for j from 1 to nops([indices(C)]) do
T[vm[j+3],funcno]:=C[j];
od;
#then call make_h_answer to do all the real work...
make_h_answer(T,m_basis-1,m_func-1,[0],[0,1],[0,2],funcno);
make_h_answer(T,m_basis-1,m_func-1,[1],[1,2],[0,1],funcno);
make_h_answer(T,m_basis-1,m_func-1,[2],[0,2],[1,2],funcno);
end;
>
#This function returns a list of which harmonic basis functions
#should be used for integrating versus two basis elements.
#It assumes that the basis elements are neighbors, because otherwise
#it's a fairly trivial question. It is called by int_vs_basis12_h.
#Notice that the values depend on orientation as well as neighbors...
which_funcs_h:=proc(S,basis1,basis2)
local out;
if last_digits(basis1[1]) = [0,1] then
if S[basis1[1],1] = basis2[1] then out:=[0,1] fi;
if S[basis1[1],2] = basis2[1] then out:=[0,2] fi;
if S[basis1[1],3] = basis2[1] then out:=[1,2] fi;
if S[basis1[1],4] = basis2[1] then out:=[1,0] fi;
fi;
if last_digits(basis1[1]) = [0,2] then
if S[basis1[1],1] = basis2[1] then out:=[2,0] fi;
if S[basis1[1],2] = basis2[1] then out:=[2,1] fi;
if S[basis1[1],3] = basis2[1] then out:=[0,1] fi;
if S[basis1[1],4] = basis2[1] then out:=[0,2] fi;
fi;
if last_digits(basis1[1]) = [1,2] then
if S[basis1[1],1] = basis2[1] then out:=[1,2] fi;
if S[basis1[1],2] = basis2[1] then out:=[1,0] fi;
if S[basis1[1],3] = basis2[1] then out:=[2,0] fi;
if S[basis1[1],4] = basis2[1] then out:=[2,1] fi;
fi;
out;
end;
>
#This function is used to generate the matrix integral(q,basisi,basisj).
#It basically has to go through the logic of inner_prod_h, except that instead
#of just returning a theoretical value, it returns the appropriate integral.
int_vs_basis12_h:=proc(T,S,m_basis,m_func,basis_1,basis_2,funct,harm_funcno)
local out,temp,basis1,basis2,fact,triangle1,triangle2,addresses1,addresses2,j,which_harms;
fact:=1/(3^m_basis); #scaling factor due to m_basis
basis1:=basis_1;
basis2:=basis_2;
# print(basis1,basis2);
#if they aren't neighbors, we get nothing...
if basis1[1] <> basis2[1] and S[basis1[1],1] <> basis2[1] and S[basis1[1],2] <> basis2[1] and S[basis1[1],3] <> basis2[1] and S[basis1[1],4] <> basis2[1] then RETURN(0) fi;
if basis1[1]=basis2[1] then
triangle1:=get_triangle(S,basis1,m_basis,1); #first, let's get the triangles, as well as the addresses of the vertices of those triangles
triangle2:=get_triangle(S,basis1,m_basis,2);
addresses1:=get_addresses(S,basis1,m_basis,1); #this returns the addresses of the coordinates of the triangles
addresses2:=get_addresses(S,basis1,m_basis,2);
#Notice that if basis1 and basis2 are in the same place, then we have to integrate over two triangles
#Another important thing is that the harmonic basis elements that you integrate against dependent solely
#upon the orientation of the basis elements, since they are in the same place.
#Notice that a lot of the work here is done in integrate_basis12, which is in splines.mws.
if last_digits(basis1[1]) = [0,1] then
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+0,harm_funcno+0)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+1,harm_funcno+1); fi;
if last_digits(basis1[1]) = [0,2] then
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+2,harm_funcno+2)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+0,harm_funcno+0); fi;
if last_digits(basis1[1]) = [1,2] then
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+1,harm_funcno+1)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+2,harm_funcno+2); fi;
fi; #END OF CASE WHERE BASIS1 and BASIS2 are in same place
#And so now what if the basis elements are neighbors
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then
#This case is if the overlapping triangle is the first triangle of basis1 (in the clockwise sense).
triangle1:=get_triangle(S,basis1,m_basis,1); #first, let's get the triangle...
addresses1:=get_addresses(S,basis1,m_basis,1); #and the corresponding address...
which_harms:=which_funcs_h(S,basis1,basis2); #This function does most of the work for determining which basis functions to use.
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct, harm_funcno+which_harms[1],harm_funcno+which_harms[2]);
else
#This case is if the overlapping triangle is the second triangle of basis1 (in the clockwise sense).
#Otherwise, it's just like before.
triangle1:=get_triangle(S,basis1,m_basis,2);
addresses1:=get_addresses(S,basis1,m_basis,2);
which_harms:=which_funcs_h(S,basis1,basis2);
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct, harm_funcno+which_harms[1],harm_funcno+which_harms[2]);
fi;
fi;
out*fact; #Don't forget the scaling factor!
end;
>
#This function generates the matrix integral(q*basisi,basisj) which is used for q-nonconstant
#in the finite element method. The matrix is sparse and symmetrical, and both of these properties are
#used here. Most of the work is done in int_vs_basis12_h.
gen_q_matrix_h:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno)
local i,j,q_matrix,q_val,n_verts;
n_verts:=nops(vm)-3;
q_matrix:=array(sparse,1..n_verts,1..n_verts);
for i from 1 to n_verts do
printf(`%d `,i);
if ((i mod 15)=1) then
printf(`\n`): fi:
for j from 1 to i do
q_val:=int_vs_basis12_h(T,S,m_basis,m_func,get_basis(vm,i),get_basis(vm,j),funcno,harm_funcno);
if q_val <> 0 then q_matrix[i,j]:=q_val; q_matrix[j,i]:=q_val;fi;
od;
od;
q_matrix;
end;
|
