 |
Here is the Maple code for the splines stuff on the gasket. Remember if you
need more info on splines you can download the "
Splines on Fractals" paper.
>
#This function returns the inner product of two basis vectors in the biharmonic basis.
#It returns theoretical values from the splines on fractals paper (scaled down by the appropriate
#factor of 3 because we are dealing with smaller triangles...). Notice there are several
#cases to deal with, and one must also be careful to use the fact that we have defined the
#normal derivative to be negative in the 2nd triangle...
inner_prod:=proc(S,m,basis_1,basis_2)
local out,temp,basis1,basis2;
basis1:=basis_1;#We do this so that we can change basis1 and basis2
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;
#first deal with special case of either basis1 or basis2 elt of v0:
if basis1[1]=[0] or basis1[1] = [1] or basis1[1] = [2] or basis2[1]=[0] or basis2[1] = [1] or basis2[1] = [2] then
#first make basis2 to be definitely in V0 so that we can keep track of it.
if nops(basis1[1])=1 then
temp:=basis2;
basis2:=basis1;
basis1:=temp;
fi;
if basis1[1]=basis2[1] then #I have included other cases here for reasons I'm not even sure of.
#but the only case which matters is g and g, because only g type basis
#elements can exist on the boundaries...
if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then out:=-47/1395 / (3^m) fi;
if basis1[2] = f and basis2[2] = f then out:=190/837 / (3^m) fi;
if basis1[2] = g and basis2[2] = g then out:=206/37665 / (3^m) fi;
else
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
#we will always make it inside of this if statement because the basis elts are always neighbors...
if basis1[2] = f and basis2[2] = f then out:=89/1674 / (3^m) fi; #This case never happens...
if basis1[2] = g and basis2[2] = g then #If we get within this case, then we have to figure out the signs
if S[basis1[1],1] = basis2[1] then out:=83/37665 / (3^m) fi; #what helps is the fact that the sign of
if S[basis1[1],4] = basis2[1] then out:=-83/37665 / (3^m) fi; #the boundary is always positive...
fi;
if (basis1[2] = f and basis2[2] = g) then out:=-61/5580 / (3^m) fi;
if (basis1[2] = g and basis2[2] = f) then #This case never happens, because basis2 can never be f
if S[basis1[1],1] = basis2[1] then out:=-61/5580 / (3^m) fi;
if S[basis1[1],4] = basis2[1] then out:=61/5580 / (3^m) fi;
fi;
fi;
fi;
else
#OTHERWISE, if neither are on the boundaries...
if basis1[1]=basis2[1] then #This is the case of the basis elts being in the same place...
if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then out:=0 fi; #by odd symmetry of g, this integral is 0
if basis1[2] = f and basis2[2] = f then out:=190/837 * 2 / (3^m) fi;
if basis1[2] = g and basis2[2] = g then out:=206/37665 * 2 / (3^m) fi;
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
#again, we will always enter this if statement because they are always neighbors...
if basis1[2] = f and basis2[2] = f then out:=89/1674 / (3^m) fi;
if basis1[2] = g and basis2[2] = g then #Here, we have to be careful of signs, using the fact that
#the 1 and 2 neighbors are in the 1 triangle, which is positive...
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=83/37665 / (3^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=-83/37665 / (3^m) fi;
fi;
if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then # out:=-83/37665 / (3^m) fi;
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=-83/37665 / (3^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=83/37665 / (3^m) fi;
fi;
fi;
if (basis1[2] = f and basis2[2] = g) then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=-61/5580 / (3^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=61/5580 / (3^m) fi;
fi;
if (basis1[2] = g and basis2[2] = f) then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=-61/5580 / (3^m) fi;
if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then out:=61/5580 / (3^m) fi;
fi;
fi;
fi;
out;
end;
>
>
#This function calculates the energy of basis1, basis2. It works in exactly the same way as the inner_prod function above,
#since the energy form is bilinear, just like the inner product. In fact, the only things which change are the theoretical
#values (again, found in the splines paper), and the scaling factor becomes 5/3...
energy:=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;
#first deal with special case of either basis1 or basis2 elt of v0:
if basis1[1]=[0] or basis1[1] = [1] or basis1[1] = [2] or basis2[1]=[0] or basis2[1] = [1] or basis2[1] = [2]
then
# print(basis1,basis2);
#first make basis1 the one not in V0...
if nops(basis1[1])=1 then
temp:=basis2;
basis2:=basis1;
basis1:=temp;
fi;
#print (basis1,basis2);
if basis1[1]=basis2[1] then
if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then out:=-7/18 * ((5/3)^m) fi;
if basis1[2] = f and basis2[2] = f then out:=19/6 * ((5/3)^m) fi;
if basis1[2] = g and basis2[2] = g then out:=5/27 * ((5/3)^m) fi;
else
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 basis1[2] = f and basis2[2] = f then out:=-19/12 * ((5/3)^m) fi;
if basis1[2] = g and basis2[2] = g then
if S[basis1[1],1] = basis2[1] then out:=-1/108 * ((5/3)^m) fi;
if S[basis1[1],4] = basis2[1] then out:=1/108 * ((5/3)^m) fi;
fi;
if (basis1[2] = f and basis2[2] = g) then out:=7/36 * ((5/3)^m) fi;
if (basis1[2] = g and basis2[2] = f) then
if S[basis1[1],1] = basis2[1] then out:=7/36 * ((5/3)^m) fi;
if S[basis1[1],4] = basis2[1] then out:=-7/36 * ((5/3)^m) fi;
fi;
fi;
fi;
# fi;
else
#OTHERWISE...
if basis1[1]=basis2[1] then
if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then out:=0 fi;
if basis1[2] = f and basis2[2] = f then out:=19/6 * 2 * ((5/3)^m) fi;
if basis1[2] = g and basis2[2] = g then out:=5/27 * 2 * ((5/3)^m) fi;
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
if basis1[2] = f and basis2[2] = f then out:=-19/12 * ((5/3)^m) fi;
# if basis1[2] = g and basis2[2] = g then
# if S[basis1[1],3] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=-1/108 * ((5/3)^m) fi;
# if S[basis1[1],1] = basis2[1] or S[basis1[1],4] = basis2[1] then out:=1/108 * ((5/3)^m) fi;
# fi;
if basis1[2] = g and basis2[2] = g then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=-1/108 * ((5/3)^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=1/108 * ((5/3)^m) fi;
fi;
if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then # out:=-83/37665 / (3^m) fi;
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=1/108 * ((5/3)^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=-1/108 * ((5/3)^m) fi;
fi;
fi;
if (basis1[2] = f and basis2[2] = g) then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=7/36 * ((5/3)^m) fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=-7/36 * ((5/3)^m) fi;
fi;
if (basis1[2] = g and basis2[2] = f) then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=7/36 * ((5/3)^m) fi;
if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then out:=-7/36 * ((5/3)^m) fi;
fi;
fi;
fi;
out;
end;
>
# After running this function, you will end up having:
# funcno, funcno+2,fucno+4 are (100),(010),(001) basis f's (old basis)
# funcno+6,funcno+8,funcno+10 are Lap(f)=(100),etc. for the old basis
# funcno+1,funcno+5,funcno+9 are f-type bases for the new basis
# funcno+3,funcno+7,funcno+11 are g-type bases for the new basis
# funcno+12+1,12+5,12+9,12+3...->12+11 are laplacians of 1,5,9,3,7,11
# (used for calculating laplacians of splines)
calc_biharms:=proc(T,R,m_harm,funcno)
#first, let's set up to calculate the old basis elements...
T[[0],funcno+0]:=1: T[[1],funcno+0]:=0: T[[2],funcno+0]:=0;
T[[0],funcno+1]:=0: T[[1],funcno+1]:=0: T[[2],funcno+1]:=0;
T[[0],funcno+2]:=0: T[[1],funcno+2]:=1: T[[2],funcno+2]:=0;
T[[0],funcno+3]:=0: T[[1],funcno+3]:=0: T[[2],funcno+3]:=0;
T[[0],funcno+4]:=0: T[[1],funcno+4]:=0: T[[2],funcno+4]:=1;
T[[0],funcno+5]:=0: T[[1],funcno+5]:=0: T[[2],funcno+5]:=0;
T[[0],funcno+6]:=0: T[[1],funcno+6]:=0: T[[2],funcno+6]:=0;
T[[0],funcno+7]:=1: T[[1],funcno+7]:=0: T[[2],funcno+7]:=0;
T[[0],funcno+8]:=0: T[[1],funcno+8]:=0: T[[2],funcno+8]:=0;
T[[0],funcno+9]:=0: T[[1],funcno+9]:=1: T[[2],funcno+9]:=0;
T[[0],funcno+10]:=0: T[[1],funcno+10]:=0: T[[2],funcno+10]:=0;
T[[0],funcno+11]:=0: T[[1],funcno+11]:=0: T[[2],funcno+11]:=1;
#Now, we'll actually calculate the bi-harmonic functions, leaving the laplacian of the function
#in funcno+1... (hence the need to space them out by 2).
n_harm(T,R,2,m_harm,[0],[1],[2],funcno);
n_harm(T,R,2,m_harm,[0],[1],[2],funcno+2);
n_harm(T,R,2,m_harm,[0],[1],[2],funcno+4);
n_harm(T,R,2,m_harm,[0],[1],[2],funcno+6);
n_harm(T,R,2,m_harm,[0],[1],[2],funcno+8);
n_harm(T,R,2,m_harm,[0],[1],[2],funcno+10);
#Here, we will apply the transforms to the new basis to the laplacians of the old basis, thereby
#creating the laplacians of the new basis functions...
f_g_lap(T,m_harm,funcno);
#and lastly, we will create the new basis elements, sticking their values in between the old basis elts (to save memory).
f_g_bih(T,m_harm,0,funcno,funcno+1);
f_g_bih(T,m_harm,2,funcno,funcno+5);
f_g_bih(T,m_harm,4,funcno,funcno+9);
end;
>
#This function is used almost everywhere. Basically, let's say you have an address
#which is given relative to a subtriangle, and then want to find out what the address
#of that point is relative to the total triangle. Then you would call this function,
#with input the relative address in triangle, and with add0,add1,add2 being the addresses
#of the 3 vertices of triangle (relative to the main triangle).
coord_transfer:=proc(input, triangle,add0,add1,add2)
if nops(input) = 1 then
if input = [0] then RETURN(add0); fi;
if input = [1] then RETURN(add1); fi;
if input = [2] then RETURN(add2); fi;
fi;
if nops(input) > 1 then
RETURN([op(triangle),op(input)]); fi;
end;
>
#Basically, given a basis element, it is often useful to determine the
#two neighboring triangles (often used in coord_transfer). Thus, this function
#takes a basis element as its input, and if which_tri=1, it returns the triangle
#address of the first triangle, and if which_tri=2, it returns the triangle address
#of the second triangle. Basically, if the basis element is an element of Vm_basis,
#then we can just cut off the last 2 digits and tack on a 0 or a 1 or a 2, depending
#on the orientation of the basis element (ie, the last 2 digits). Otherwise,
#we can use the fact that the neighbors must be in Vm_basis, so we can choose the
#appropriate neighbor and chop off some digits, and so forth...
get_triangle:=proc(S,basis,m_basis,which_tri)
if which_tri = 1 then
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),1]);
else RETURN([op(S[basis[1],1][1..nops(S[basis[1],1])-2]),0]); fi;
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),2]);
else RETURN([op(S[basis[1],1][1..nops(S[basis[1],1])-2]),1]); fi;
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),0]);
else RETURN([op(S[basis[1],1][1..nops(S[basis[1],1])-2]),2]); fi;
fi;
else
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),0]);
else RETURN([op(S[basis[1],4][1..nops(S[basis[1],4])-2]),1]); fi;
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),1]);
else RETURN([op(S[basis[1],4][1..nops(S[basis[1],4])-2]),2]); fi;
fi;
if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then
if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),2]);
else RETURN([op(S[basis[1],4][1..nops(S[basis[1],4])-2]),0]); fi;
fi;
fi;
end;
>
#This procedure takes a function and returns integral(funct*basis).
#the basic strategy is to determine which triangles basis is non zero on
#and then use coord_transfer to determine which values of funct to multiply
#the value of the appropriate basis function by and then integrate.
int_vs_basis:=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); #factor for simpson's rule...
#get the vertices (ie, V0, Vm, Veverything in between)...
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);
if nops(basis[1]) = 1 then #this is if the basis is an element of V0...
#Here, we basically know the basis element is a g type vector. Then, we use the fact that
#the triangle in question is just a sequence of 0's or 1's. Then, we multiply by the appropriate
#biharmonic basis function (ie, harm_funcno+3 or some such number), and sum...
if basis[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
V0sum1:=add(T[coord_transfer(V0[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funct]*( T[V0[i],harm_funcno+3] ),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+3]),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+3]),i=1..nops(Vm));
fi;
if basis[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*( T[V0[i],harm_funcno+7]),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+7]),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+7]),i=1..nops(Vm));
fi;
if basis[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),funct]*(T[V0[i],harm_funcno+11]),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+11] ),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+11]),i=1..nops(Vm));
fi;
RETURN(fact * (5*(Vmsum1) + 2*(Vothersum1) + (V0sum1))); #And then return the simpson's rule formula...
else #So we will come here if the basis element is not on the boundary (could be either f or g...)
#first let's deal with the basis element being an f...
if basis[2]=f then
#So now we have to get 2 different triangles (the two triangles associated with the vertex)
#Hence we have 2 calls to get_triangle. So we have 2 sums, one for triangle 1 and one for triangle2.
#We also need to keep track of the orientation (ie, the last 2 digits of the basis) in order to
#determine which basis functions we want to use.
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+1] ),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+1] ),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+1] ),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+5]),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+5]),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+5]),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+9]),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+9]),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+9]),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+1]),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+1]),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+1]),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+5]),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+5]),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+5]),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+9]),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+9]),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+9]),i=1..nops(Vm));
fi;
fi;
#now let us deal with the basis element being a g...
if basis[2]=g then
#The code is exactly the same as for the f-type vector, except that
#we use g-type basis vectors (which basically means adding 2 to the
#harm_funcno+k part).
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+3]),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+3]),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+3]),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+7]),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+7]),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+7]),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+11]),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+11]),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+11]),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+3]),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+3]),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+3]),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+7]),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+7]),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+7]),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+11]),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+11]),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+11]),i=1..nops(Vm));
fi;
fi;
#Return simpson's rule. Notice that the second sum must be subtracted if it's a g-type vector
#because the normal derivative is defined as negative in the 2nd triangle...
if basis[2] = f then RETURN(fact * (5*(Vmsum1+Vmsum2) + 2*(Vothersum1+Vothersum2) + (V0sum1+V0sum2))); fi;
if basis[2] = g then RETURN(fact * (5*(Vmsum1-Vmsum2) + 2*(Vothersum1-Vothersum2) + (V0sum1-V0sum2))); fi;
fi;
end;
>
#so this function returns integral(funct*basis1*basis2).
#it is used to calculate the q-matrix. This function basically is like
#a combination of the logic of inner_prod with the integration of int_vs_basis.
#several help functions are defined to help do this (it's a pretty complicated function).
int_vs_basis12:=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); #integration scaling factor...
basis1:=basis_1; #again, so that we can switch around basis1 and basis2, if need be.
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;
#first deal with special case of either basis1 or basis2 elt of v0:
if basis1[1]=[0] or basis1[1] = [1] or basis1[1] = [2] or basis2[1]=[0] or basis2[1] = [1] or basis2[1] = [2] then
#first make basis1 the one not in V0...
#ie, make basis2 the one which is in V0 for sure
if nops(basis1[1])=1 then
temp:=basis2;
basis2:=basis1;
basis1:=temp;
fi;
if basis1[1]=basis2[1] then #The only way we'll ever get here is if both elements are g-type basis vectors
if basis1[2] = g and basis2[2] = g then
if basis1[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
#The function integrate_basis12 basically does all the work; alls we have to do is tell it which functions...
out:=integrate_basis12(T,m_func-m_basis,triangle1,basis1[1],S[basis1[1],1],S[basis1[1],2],funct,harm_funcno+3,harm_funcno+3);
fi;
if basis1[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis1[1],2],basis1[1],S[basis1[1],1],funct,harm_funcno+7,harm_funcno+7);
fi;
if basis1[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis1[1],1],S[basis1[1],2],basis1[1],funct,harm_funcno+11,harm_funcno+11);
fi;
fi;
else
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
# Notice that the only way we'll get in here is if basis2 is a g-type vector; however, basis1 can still be either g or f...
if basis1[2] = g and basis2[2] = g then
if S[basis1[1],1] = basis2[1] then
if basis2[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],funct,harm_funcno+3,harm_funcno+11);
fi;
if basis2[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+7,harm_funcno+3);
fi;
if basis2[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcno+11,harm_funcno+7);
fi;
fi;
if S[basis1[1],4] = basis2[1] then
if basis2[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
out:=-integrate_basis12(T,m_func-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],funct,harm_funcno+3,harm_funcno+7);
fi;
if basis2[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
out:=-integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+7,harm_funcno+11);
fi;
if basis2[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
out:=-integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcno+11,harm_funcno+3);
fi;
fi;
fi;
if (basis1[2] = f and basis2[2] = g) then #This is the only remaining case, because f-type vectors can never be at the boundaries
if S[basis1[1],1] = basis2[1] then
if basis2[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],funct,harm_funcno+3,harm_funcno+9);
fi;
if basis2[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+7,harm_funcno+1);
fi;
if basis2[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcno+11,harm_funcno+5);
fi;
fi;
if S[basis1[1],4] = basis2[1] then
if basis2[1] = [0] then
triangle1:=[seq(0,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],funct,harm_funcno+3,harm_funcno+5);
fi;
if basis2[1] = [1] then
triangle1:=[seq(1,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+7,harm_funcno+9);
fi;
if basis2[1] = [2] then
triangle1:=[seq(2,j=1..m_basis)];
out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcno+11,harm_funcno+1);
fi;
fi;
fi;
fi;
fi;
# fi;
else #i.e, neither are in V0....
if basis1[1]=basis2[1] then #they're at the same vertex.
#here, the only thing that matters is the orientation to determine which basis functions to use...
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);
addresses2:=get_addresses(S,basis1,m_basis,2);
if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then #out:=0 fi;
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+1,harm_funcno+3)-
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+5,harm_funcno+7); 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+9,harm_funcno+11)-
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+1,harm_funcno+3); 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+5,harm_funcno+7)-
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+9,harm_funcno+11); fi;
fi;
if basis1[2] = f and basis2[2] = f then #out:=190/837 * 2 / (3^m) fi;
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+1,harm_funcno+1)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+5,harm_funcno+5); 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+9,harm_funcno+9)+
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]) = [1,2] then
out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+5,harm_funcno+5)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+9,harm_funcno+9); fi;
fi;
if basis1[2] = g and basis2[2] = g then #out:=206/37665 * 2 / (3^m) fi;
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+3,harm_funcno+3)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+7,harm_funcno+7); 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+11,harm_funcno+11)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+3,harm_funcno+3); 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+7,harm_funcno+7)+
integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+11,harm_funcno+11); fi;
fi;
fi; #END OF CASE WHERE BASIS1 and BASIS2 are in same place
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 they are neighbors, then we're in here. Notice that we just call
#which_funcs to determine which basis functions to use, as well as getting the negative sign to make the
#sum negative if we are in a negative normal derivative region...
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = 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
addresses1:=get_addresses(S,basis1,m_basis,1);
which_harms:=which_funcs(S,basis1,basis2);
out:=which_harms[3] * 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
triangle1:=get_triangle(S,basis1,m_basis,2);
addresses1:=get_addresses(S,basis1,m_basis,2);
which_harms:=which_funcs(S,basis1,basis2);
out:=which_harms[3] * 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;
fi;
out*fact;
end;
>
#This function basically integrates funcno*h_func1*h_func2 over triangle.
#it works by using coordinates relative to the subtriangle triangle, and
#then using coord_transfer to determine which address of funcno in the
#main triangle to multiply by. Used in both int_vs_basis12 and int_vs_basis12_h...
integrate_basis12:=proc(T,m,triangle,add0,add1,add2,funcno,h_func1,h_func2)
local fact,V0,Vm, Vother_e,Vother,V0sum1,Vothersum1,Vmsum1,i;
fact:=1/(6*3^m);
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);
V0sum1:=add(T[coord_transfer(V0[i],triangle,add0,add1,add2),funcno]*T[V0[i],h_func1]*T[V0[i],h_func2],i=1..nops(V0));
Vothersum1:=add(T[coord_transfer(Vother[i],triangle,add0,add1,add2),funcno] *T[Vother[i],h_func1]*T[Vother[i],h_func2],i=1..nops(Vother));
Vmsum1:=add(T[coord_transfer(Vm[i],triangle,add0,add1,add2),funcno]*T[Vm[i],h_func1]*T[Vm[i],h_func2],i=1..nops(Vm));
RETURN(fact*(5*Vmsum1 + 2 * Vothersum1 + V0sum1));
end;
>
#This function is called by int_vs_basis12.
#it basically assumes basis1 and basis2 are neighbors,
#and using that, determines which harmonic basis functions
#you would want to multiply against when you integrate.
#it uses the orientation of the vertex, as well as the
#f or g status. Also, it includes, as the last element of the
#output list, a multiplicative factor of 1 or -1, depending on whether
#you want the positive or the negative of the resultant integral. This
#sign depends on the g-vectors and which triangle you are in, relatively
#speaking...
which_funcs:=proc(S,basis1,basis2)
local out;
if last_digits(basis1[1]) = [0,1] then
if S[basis1[1],1] = basis2[1] then out:=[1,5] fi;
if S[basis1[1],2] = basis2[1] then out:=[1,9] fi;
if S[basis1[1],3] = basis2[1] then out:=[5,9] fi;
if S[basis1[1],4] = basis2[1] then out:=[5,1] fi;
fi;
if last_digits(basis1[1]) = [0,2] then
if S[basis1[1],1] = basis2[1] then out:=[9,1] fi;
if S[basis1[1],2] = basis2[1] then out:=[9,5] fi;
if S[basis1[1],3] = basis2[1] then out:=[1,5] fi;
if S[basis1[1],4] = basis2[1] then out:=[1,9] fi;
fi;
if last_digits(basis1[1]) = [1,2] then
if S[basis1[1],1] = basis2[1] then out:=[5,9] fi;
if S[basis1[1],2] = basis2[1] then out:=[5,1] fi;
if S[basis1[1],3] = basis2[1] then out:=[9,1] fi;
if S[basis1[1],4] = basis2[1] then out:=[9,5] fi;
fi;
#if they are both f's...
if basis1[2] = f and basis2[2] = f then out:=[op(out),1]; fi;
#if 1 is an f and the other is a g...
if basis1[2] = f and basis2[2] = g then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=[op(out),1];
else out:=[op(out),-1]; fi;
fi;
if basis1[2] = g and basis2[2] = f then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=[op(out),1];
else out:=[op(out),-1]; fi;
fi;
#if they're both g's (word!)...
if basis1[2] = g and basis2[2] = g then
if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=[op(out),1] fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=[op(out),-1] fi;
fi;
if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then
if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=[op(out),-1] fi;
if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=[op(out),1] fi;
fi;
fi;
if basis1[2] = g then out[1]:=out[1] + 2; fi;
if basis2[2] = g then out[2]:=out[2] + 2; fi;
out;
end;
|
