 |
Well, here is the Maple code for various functions used throughout
our programming:
>
#This function multiplies func1 and func2 pointwise and stores
#the result in outfunc...
mult:=proc(T,m,func1,func2,outfunc)
local inds,i;
inds:=vertices(m);
for i from 1 to nops(inds) do
T[inds[i],outfunc]:=T[inds[i],func1]*T[inds[i],func2];
od;
end;
>
#This function takes the sine of func1 (pointwise) and stores
#it in outfunc.
sinof:=proc(T,m,func1,outfunc)
local inds,i;
inds:=vertices(m);
for i from 1 to nops(inds) do
T[inds[i],outfunc]:=sin(T[inds[i],func1]);
od;
end;
>
#This procedure takes func1,func2,mult1, and mult2, and stores
# mult1*func1+mult2*func2 in outfunc (again, pointwise).
add_func_lin:=proc(T,m,func1,func2,mult1,mult2,outfunc)
local inds,i;
inds:=vertices(m);
for i from 1 to nops(inds) do
T[inds[i],outfunc]:=mult1*T[inds[i],func1]+mult2*T[inds[i],func2];
od;
end;
>
#So this function stores the laplacian of harm_funcno+(odd number up to 11)
#and stores it in harm_funcno+12+(odd number up to 11). It stores this by
#using the laplacians of the "old" basis functions, which are stored in
#harm_funcno+(odd number up to 11) (values are overridden later).
#See calc_biharms in splines.mws for full details on this process.
f_g_lap:=proc(T,m,harm_funcno)
local inds,i;
inds:=vertices(m);
for i from 1 to nops(inds) do
#Uses the 30-15-15 rule...
T[inds[i],harm_funcno + 12+1]:=T[inds[i],harm_funcno+1]-30*T[inds[i],harm_funcno+7] +15 * T[inds[i],harm_funcno + 9]+15 * T[inds[i],harm_funcno + 11];
T[inds[i],harm_funcno + 12+5]:=T[inds[i],harm_funcno+3]-30*T[inds[i],harm_funcno+9] +15 * T[inds[i],harm_funcno + 7]+15 * T[inds[i],harm_funcno + 11];
T[inds[i],harm_funcno + 12+9]:=T[inds[i],harm_funcno+5]-30*T[inds[i],harm_funcno+11] +15 * T[inds[i],harm_funcno + 9]+15 * T[inds[i],harm_funcno + 7];
T[inds[i],harm_funcno + 12+3]:=11*T[inds[i],harm_funcno + 7] -4 * T[inds[i],harm_funcno + 9]-4 * T[inds[i],harm_funcno + 11];
T[inds[i],harm_funcno + 12+7]:=11*T[inds[i],harm_funcno+9] -4 * T[inds[i],harm_funcno + 7]-4 * T[inds[i],harm_funcno + 11];
T[inds[i],harm_funcno + 12+11]:=11*T[inds[i],harm_funcno+11] -4 * T[inds[i],harm_funcno + 9]-4 * T[inds[i],harm_funcno + 7];
od;
end;
#This function creates the f and the g type biharmonic splines
#from the old basis splines by using the 30-15-15 rule and the 4-11-11 rule.
#It puts the f-type vector in out_funcno, and then puts
#the g-type vector in out_funcno+2.
f_g_bih:=proc(T,m,num,harm_funcno,out_funcno)
local inds,i;
inds:= vertices(m);
for i from 1 to nops(inds) do
T[inds[i],out_funcno]:=T[inds[i],harm_funcno+num]-30*T[inds[i],harm_funcno+6+num] + 15 * T[inds[i],harm_funcno+((8+num) mod 6) + 6] + 15 * T[inds[i],harm_funcno+((10+num) mod 6) + 6];
T[inds[i],out_funcno+2]:=11 * T[inds[i],harm_funcno+((6+num) mod 6) + 6 ] -4 * T[inds[i],harm_funcno+((8+num) mod 6) + 6]-4 * T[inds[i],harm_funcno+((10+num) mod 6) + 6];
od;
end;
>
#This function implements the trapzoid rule for integration (out to level m).
#It should be noted that simpson's rule is significantly better
#than the trapezoid rule, and simp_rule should be used whenever
#possible.
trap_rule:=proc(T,m,funcno)
local i, func, V0, trap,Vother,Vother_e;
V0:=vertices_m(0);
Vother_e:=vertices(m); Vother:=Vother_e[4..nops(Vother_e)];
trap := eval(3^(-m-1) * ( sum(2 * T[Vother[i],funcno], i = 1..nops(Vother) ) + sum(T[V0[i],funcno], i = 1..nops(V0)) ));
trap;
end;
>
#This function implements simpson's rule. This function gives
#much better results than trap_rule. It should be noted that this algorithm is
#used for all the numerical integrations performed when implementing the FEM.
#However, most of the time, instead of using simp_rule, the algorithm is
#implemented within the procedure itself as a matter of convenience and efficiency.
simp_rule:=proc(T,m,funcno)
local i, func, V0,Vm, simp,Vother_e,Vother;
#get V0, Vm, Veverythinginbetween...
V0:=vertices_m(0);
if m =1 then Vother:=[] else Vother_e:=vertices(m-1); Vother:=Vother_e[4..nops(Vother_e)]; fi;
Vm:=vertices_m(m);
#and then perform the sums...
simp := eval( 1/(6*3^m)) * ( 5* sum(T[Vm[i],funcno], i = 1..nops(Vm) ) + 2*sum(T[Vother[i],funcno], i = 1..nops(Vother)) + sum(T[V0[i],funcno],i=1..nops(V0)) );
simp;
end;
>
#This function takes the laplacian of funcno and stores it in funcno_out.
#Make_neighbors must be run on T before this function is run.
#Also, it should be noted that the laplacian doesn't exist on the
#boundaries of the gasket when using the pointwise formula used here;
#instead, I just set the values on the edges to 0 for ease of programming.
#It can be changed quite easily...
find_lap:=proc(T,m,funcno,funcno_out)
local i, j, Vm;
Vm:=vertices(m);
#assign the values on the boundaries...
for i from 1 to 3 do
T[Vm[i],funcno_out]:=0; od;
# T[Vm[i], funcno_out] := (3/2) * (5^(m)) * (add( T[T[Vm[i],j],funcno] - (2 * T[Vm[i],funcno]), j = 1..2)); od;
#and now do the rest...
for i from 4 to nops(Vm) do
T[Vm[i], funcno_out] := (3/2) * (5^(m)) * add( T[T[Vm[i],j],funcno] - T[Vm[i],funcno], j = 1..4); od;
end;
>
#This function returns the last 2 digits of address.
#It is used to determine orientation.
last_digits:=proc(address)
address[nops(address)-1..nops(address)]; end;
>
#This function returns the addresses of the boundaries of the 2
#triangles of which basis[1] is a common vertex. It requires that the
#spline basis table S has been make_neighbored. If which_triangle=1 then
#it returns the first triangle (in clockwise orientation), otherwise it returns
#the second triangle (in clockwise orientation). It returns the addresses in
#clockwise fashion, with the first vertex being the one furthest up (ie,
#largest y value in the euclidean sense). These addresses are often used in
#coord_transfer.
get_addresses:=proc(S,basis,m_basis,which_triangle)
local digits;
digits:=last_digits(basis[1]);
if digits=[0,1] then
if which_triangle=1 then RETURN([basis[1],S[basis[1],1],S[basis[1],2]]);
else RETURN([S[basis[1],4],basis[1],S[basis[1],3]]); fi;
fi;
if digits = [0,2] then
if which_triangle=1 then RETURN([S[basis[1],1],S[basis[1],2],basis[1]]);
else RETURN([basis[1],S[basis[1],3],S[basis[1],4]]); fi;
fi;
if digits = [1,2] then
if which_triangle=1 then RETURN([S[basis[1],2],basis[1],S[basis[1],1]]);
else RETURN([S[basis[1],3],S[basis[1],4],basis[1]]); fi;
fi;
end;
>
#This generates a function with values 0 on (V1-V0), and sets the boundaries to value,
#and then fills in the rest of the function harmonically, thus generating
#a piecewise harmonic function with full D3 symmetry.
sym_func:=proc(T,m,value,funcno)
T[[0],funcno]:=value;
T[[1],funcno]:=value;
T[[2],funcno]:=value;
T[[0,1],funcno]:=0;
T[[0,2],funcno]:=0;
T[[1,2],funcno]:=0;
harm(T,m-1,[2],[0,2], [1,2], funcno);
harm(T,m-1,[0],[0,1], [0,2], funcno);
harm(T,m-1,[1],[1,2], [0,1], funcno);
end;
>
#This procedure takes funcno, and (given IPinv = inverse(IP) where
#IP is the inner product matrix created by gen_iprod_matrix) expresses
#the function in the spline basis.
express_in_basis:=proc(T,IPinv,S,m_basis,m_func,funcno,harm_funcno)
local vec,vm;
vm:=vertices(m_basis);
vec:=gen_f_vector(T,S,m_basis,m_func,vm,funcno,harm_funcno);
evalm(IPinv&*vec);
end;
>
#This function takes a vector of length n_elts
#and writes it to a file. The subsequent
#column vector can be read into matlab quite easily.
write_vec:=proc(filename,vec,n_elts)
local fd,i;
fd:=fopen(filename,WRITE);
for i from 1 to n_elts do
fprintf(fd,"%g\n",vec[i]);
od;
fclose(fd);
end;
>
#This function simply takes funcno1 and copies it to funcno2
copy_function:=proc(T1,T2,m,funcno1,funcno2)
local i,vm;
vm:=vertices(m);
for i from 1 to nops(vm) do
T2[vm[i],funcno2]:=T1[vm[i],funcno1];
od;
end;
>
#This function will read (out of filename), and put the
#results in a list of lists. Each sublist will contain n_elements
#and corresponds to a row of the file.
#ie, if your file looks like:
#1 2
#3 4
#5 6
#then running read_basis_vals returns:
#[[1,2],[3,4],[5,6]]
#This is used for reading in time-series solutions to the wave
#and heat equations.
read_basis_vals:=proc(filename,n_elements)
local i,fd,vals,line,outvals;
fd:=fopen(filename,READ);
vals:=readdata(fd,n_elements,float);
fclose(fd);
vals;
end;
|
