 |
Here is the Maple code for working out the FEM using Bi-Harmonic splines:
>
with(plots):with(linalg):
>
#This procedure returns the matrix E(fi,fj) where fi, fj are spline basis functions.
#it stores values in sparse matrix form, and also takes advantage of the fact that the
#matrix is symmetric...
gen_energy_matrix:=proc(S,m_basis,vm)
local i,j,en_matrix,en;
en_matrix:=array(sparse,1..3^(m_basis+1),1..3^(m_basis+1));
for i from 1 to 3^(m_basis+1) do
for j from 1 to i do
en:=energy(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 the biharmonic basis element corresponding to number. The bases are numbered by the
#ordering provided by vm: the first 3 elements of vm are the boundaries, and there is no f type basis vectors
#at those points. The first bunch of basis vectors are basically f-type basis elements, and the rest are g-type.
get_basis:=proc(vm,number)
option remember;
if number <= nops(vm)-3 then [vm[number+3],f]
else [vm[number-nops(vm)+3],g] fi; end;
>
#This operates just the same way as gen_energy_matrix, except using inner products...
gen_iprod_matrix:=proc(S,m_basis,vm)
local i,j,i_matrix,ip;
i_matrix:=array(sparse,1..3^(m_basis+1),1..3^(m_basis+1));
for i from 1 to 3^(m_basis+1) do
for j from 1 to i do
ip:=inner_prod(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 the f vector, which is basically f integrates versus all the different basis functions.
#It calls int_vs_basis, which contains the meat of the functionality.
gen_f_vector:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno)
local i,f_vec;
f_vec:=array(1..3^(m_basis+1));
for i from 1 to 3^(m_basis+1) do
printf(`%d `,i);
if ((i mod 15)=1) then
printf(`\n`): fi:
f_vec[i]:=int_vs_basis(T,S,m_basis,m_func,get_basis(vm,i),funcno,harm_funcno);
od;
f_vec;
end;
>
>
#This function creates a actual function (which can then be plotted and the such) from the values
#of the spline projection stored in C. It creates the function simply by going through each basis
#element and adding in the appropriate amount of it. Sounds, simple, huh?
make_pretty_picture:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno,C)
local i,j,vm2,k,vmh,triangle1,triangle2,basis,add0,add1,add2,which_harm1,which_harm2,coord1,coord2;
#first, let's set the initial values of our function to 0.
vm2:=vertices(m_func);
for i from 1 to nops(vm2) do
T[vm2[i],funcno]:=0;
od;
vmh:=vertices(m_func-m_basis); #These are the vertices which we will need from our spline basis elements
for i from 1 to 3^(m_basis+1) do
basis:=get_basis(vm,i); #This is the basis element we will be working with.
printf(`%d `,i);
if ((i mod 15)=1) then
printf(`\n`): fi:
#So if the basis element is in V0, we know it is a g-type basis vector.
#the scheme here is to find the appropriate triangle (which will be some
#sequence of 0's or 1's or 2's) and also get the coordinates of the boundaries
#of this triangle (using the neighbors defined in S; like S[basis[1],2]).
#also, depending on the orientation of the basis, we need to pick the appropriate
#basis element (so for basis = [[0],g], we need a basis element which has values
#of zero everywhere, and has normal derivative 1 at [0], and normal derivative 0
#at the other boundaries. This is stored in harm_funcno+3). This explains the assignments
#of triangle1, add0, add1, add2, and which_harm1. The function coord_transfer is then used
#to take a coordinate from the basis function and transform it to the appropriate part
#of the function we are trying to generate (basically, by adding some more digits to
#the front). Then add in the value of the basis function * C[i] and voila!
if nops(basis[1]) = 1 then #this is if the basis is an element of V0...
if basis[1] = [0] then
triangle1:=[seq(0,k=1..m_basis)];
add0:=basis[1]; add1:=S[basis[1],1]; add2:=S[basis[1],2];which_harm1:=3;
fi;
if basis[1] = [1] then
triangle1:=[seq(1,k=1..m_basis)];
add0:=S[basis[1],2]; add1:=basis[1]; add2:=S[basis[1],1];which_harm1:=7;
fi;
if basis[1] = [2] then
triangle1:=[seq(2,k=1..m_basis)];
add0:=S[basis[1],1]; add1:=S[basis[1],2]; add2:=basis[1];which_harm1:=11;
fi;
for j from 1 to nops(vmh) do
coord1:=coord_transfer(vmh[j],triangle1,add0,add1,add2);
T[coord1,funcno]:=T[coord1,funcno]+C[i]*T[vmh[j],harm_funcno + which_harm1];
od;
else
#Here, the strategy is the same, except that we have 2 triangles to deal with: one on each side of the
#basis vertex. Notice that the values we feed into coord_transfer as the triangle coordinates
#are different depending on the orientation of the basis element (ie, what the last 2 digits are).
#Also, the basis functions we want to use change as well, again depending on the last 2 digits.
#All of this hinges around the fact that the neighbors are ordered in a clockwise fashion, so that the
#orientation of the basis element tells you which neighbors are the vertices of which triangle.
#first let's deal with the basis element being an f...
if basis[2]=f then
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funcno]:= T[coord_transfer(vmh[j],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funcno] + C[i] * (T[vmh[j],harm_funcno+1]);
T[coord_transfer(vmh[j],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno]:= T[coord_transfer(vmh[j],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno] + C[i] * (T[vmh[j],harm_funcno+5]);
od;
T[basis[1],funcno]:=T[basis[1],funcno]-C[i]; #we added in the value of the basis here twice, so subtract one away
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funcno]:= T[coord_transfer(vmh[j],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funcno] + C[i] * (T[vmh[j],harm_funcno+9]);
T[coord_transfer(vmh[j],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno]:= T[coord_transfer(vmh[j],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno] + C[i] * (T[vmh[j],harm_funcno+1]);
od;
T[basis[1],funcno]:=T[basis[1],funcno]-C[i];
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funcno]:= T[coord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funcno] + C[i] * (T[vmh[j],harm_funcno+5]);
T[coord_transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funcno]:= T[coord_transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funcno] + C[i] * (T[vmh[j],harm_funcno+9]);
od;
T[basis[1],funcno]:=T[basis[1],funcno]-C[i];
fi;
fi;
#now let us deal with the basis element being a g. The situation here is identical to the f situation, except
#that the basis functions change (increase by 2), and in triangle 2, we must subtract instead of add
#because the basis values are negative there...
if basis[2]=g then
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funcno]:= T[coord_transfer(vmh[j],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funcno] + C[i] * (T[vmh[j],harm_funcno+3]);
T[coord_transfer(vmh[j],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno]:= T[coord_transfer(vmh[j],triangle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno] - C[i] * (T[vmh[j],harm_funcno+7]);
od;
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funcno]:= T[coord_transfer(vmh[j],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funcno] + C[i] * (T[vmh[j],harm_funcno+11]);
T[coord_transfer(vmh[j],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno]:= T[coord_transfer(vmh[j],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno] - C[i] * (T[vmh[j],harm_funcno+3]);
od;
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);
for j from 1 to nops(vmh) do
T[coord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funcno]:= T[coord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funcno] + C[i] * (T[vmh[j],harm_funcno+7]);
T[coord_transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funcno]:= T[coord_transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funcno] - C[i] * (T[vmh[j],harm_funcno+11]);
od;
fi;
fi;
fi;
od;
end;
>
#This function generates the matrix for non-constant q (ie, integral(q*basisi*basisj) ).
#Again, the meat of the functionality is in int_vs_basis12, which performs the actual integration
gen_q_matrix:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno)
local i,j,q_matrix,q_val;
q_matrix:=array(sparse,1..3^(m_basis+1),1..3^(m_basis+1));
for i from 1 to 3^(m_basis+1) 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(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;
>
#This function creates the laplacian of the function defined by the spline
#projection C. Notice that the values at the Vm_basis are completely wrong
#because the laplacian is discontinuous at those points. See the posteriori
#error function in error_funcs.mws for a method of storing the discontinuous
#values of the laplacian at the basis vertices. Notice that the way this
#function works is just by calling make_pretty picture, but instead of sending
#the biharmonic basis elements, we send it the laplacians of the biharmonic
#basis elements (hence harm_funcno+12).
make_pretty_picture_lap:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno,C)
local i,vmfunc;
vmfunc:=vertices(m_func);
make_pretty_picture(T,S,m_basis,m_func,vm,funcno,harm_funcno+12,C);
for i from 1 to nops(vmfunc) do
T[vmfunc[i],funcno]:=T[vmfunc[i],funcno]*5^m_basis;
od;
end;
|
