#this procedure generates all words of length m in L letters... p:=proc(L,m) option remember; local L1; if m = 1 then L else L1 := p(L,m-1); [op(map(add_one,L1,0)),op(map(add_one,L1,1)),op(map(add_one,L1,2))] fi; end; add_one:=(u,d)-> [op(u),d]; #This function returns all vertices in Vm in a big list vertices:= proc(m) local i,outlist; option remember; if m = 1 then RETURN([[0],[1],[2],[0,1],[0,2],[1,2]]) fi; if m=2 then [[0],[1],[2],[0,1],[0,2],[1,2],[0,0,1],[0,0,2],[0,1,2],[1,0,1],[1,0,2],[1,1,2],[2,0,1],[2,0,2],[2,1,2]]; else outlist:=[[0] , [1] , [2] , [0,1], [0,2], [1,2]]; for i from 2 to m do outlist:=[op(outlist),op(map(add_pair, p([[0],[1],[2]], (i-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (i-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (i-1)), [1,2]))]; od; outlist; #[[0] , [1] , [2] , [0,1], [0,2], [1,2], op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [1,2]))]; fi; end; #vertices_m creates just the vertices in Vm, excluding Vk, k 1 then [op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,1])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,2])), op(map(add_pair, p([[0],[1],[2]], (m-1)), [1,2]))]; fi; end; add_pair := (u,v) -> [op(u), op(v)]; #This function returns a table of the gasket out to level m... gasket:=proc(m) local i, vm, T; vm:= vertices(m); T:=table(); for i from 1 to nops(vm) do T[vm[i], 0] := vm[i] od; T; end; #This function makes stores the neighbors of each vertex out to level m in #the gasket table T. The neighbors are ordered clockwise; i.e., the 1st neighbor #of 0,1 in V3 is 1,0,0,1, etc... make_neighbors:=proc(T,m,vm) local i,j,depth,nb1,nb2,nb3,nb4,iverts,last_verts; #first, we have to deal with the special case of V1... if m=1 then print(`Yo, V1 is STOOOOOOPID!`); T[[0],1]:=[0,1]; T[[0],2]:=[0,2]; T[[1],1]:=[1,2]; T[[1],2]:=[0,1]; T[[2],1]:=[0,2]; T[[2],2]:=[1,2]; T[[0,1],3]:=[0,2]; T[[0,1],4]:=[0]; T[[0,1],1]:=[1]; T[[0,1],2]:=[1,2]; T[[1,2],3]:=[0,1]; T[[1,2],4]:=[1]; T[[1,2],2]:=[0,2]; T[[1,2],1]:=[2]; T[[0,2],1]:=[0]; T[[0,2],2]:=[0,1]; T[[0,2],3]:=[1,2]; T[[0,2],4]:=[2]; #hard coding reverse neighbors... T[T[[0],1],4] := [0]; T[T[[0],2],1] := [0]; T[T[[1],1],4] := [1]; T[T[[1],2],1] := [1]; T[T[[2],1],4] := [2]; T[T[[2],2],1] := [2]; #done with hard coding (yay!) fi; if m <>1 then #but if we are not in V1... #hard coding the neighbors for V0,V1 T[[0],1]:=[seq(0,i=1..m),1]; T[[0],2]:=[seq(0,i=1..m),2]; T[[1],1]:=[seq(1,i=1..m),2]; T[[1],2]:=[seq(1,i=1..m-1),0,1]; T[[2],1]:=[seq(2,i=1..m-1),0,2]; T[[2],2]:=[seq(2,i=1..m-1),1,2]; T[[0,1],3]:=[0,seq(1,i=1..m-2),1,2]; T[[0,1],4]:=[0,seq(1,i=1..m-2),0,1]; T[[0,1],1]:=[1,seq(0,i=1..m-2),0,1]; T[[0,1],2]:=[1,seq(0,i=1..m-2),0,2]; T[[1,2],3]:=[1,seq(2,i=1..m-2),0,2]; T[[1,2],4]:=[1,seq(2,i=1..m-2),1,2]; T[[1,2],2]:=[2,seq(1,i=1..m-2),0,1]; T[[1,2],1]:=[2,seq(1,i=1..m-2),1,2]; T[[0,2],1]:=[0,seq(2,i=1..m-2),0,2]; T[[0,2],2]:=[0,seq(2,i=1..m-2),1,2]; T[[0,2],3]:=[2,seq(0,i=1..m-2),0,1]; T[[0,2],4]:=[2,seq(0,i=1..m-2),0,2]; #hard coding reverse neighbors... T[T[[0],1],4] := [0]; T[T[[0],2],1] := [0]; T[T[[1],1],4] := [1]; T[T[[1],2],1] := [1]; T[T[[2],1],4] := [2]; T[T[[2],2],1] := [2]; T[T[[0,1],1],4] := [0,1]; T[T[[0,1],2],1] := [0,1]; T[T[[0,1],3],4] := [0,1]; T[T[[0,1],4],1] := [0,1]; T[T[[1,2],1],4] := [1,2]; T[T[[1,2],2],1] := [1,2]; T[T[[1,2],3],4] := [1,2]; T[T[[1,2],4],1] := [1,2]; T[T[[0,2],1],4] := [0,2]; T[T[[0,2],2],1] := [0,2]; T[T[[0,2],3],4] := [0,2]; T[T[[0,2],4],1] := [0,2]; #done with hard coding (yay!) # First, we compute the neighbors of all pts in Vk, 1= minlen and nops(L[1]) <= maxlen); end; #This function is called by make_points. It calls itself recursively 'n stuff. gen_points:=proc(T,m,add0,add1,add2) if m <> 0 then #So the general strategy is as follows: We can't assign values to the next level of #vertices without knowing the orientation of the triangle (because we need to know whether we need to assign #to a 0,1 point, 1,2 point, etc.). So what we do is call the function with 3 addresses #and make it so that the 0th address (add0) is one of the vertices of the outside triangle #(i.e., going from the V0 triangle, the triangles would be T1=([0],[0,1],[0,2]),T2=([1],[1,2],[0,1]) #and so forth). Thus, by looking at add1 of this triangle, we can determine the orientation of #the triangle by checking whether add1 ends in 0,1, 0,2 or 1,2. Then we assign the values #accordingly, after which we call gen_points recursively, making sure to keep the orientation #of the 3 addresses in the recursive calls. We use this same method to generate harmonic functions #and n-harmonic functions, as well as several other things. Notice that this recursive process #has it's base case as V1, therefore requiring 3 seperate calls to gen_points from make_points. #This is due to the idiosyncratic addressing of the vertices of V0. Oh well.... if add1[nops(add1)-1..nops(add1)] = [0,1] then T[[op(add1[1..nops(add1)-2]),0,0,1],5] := .5 * T[add0,5] + .5 * T[add1,5]; T[[op(add1[1..nops(add1)-2]),0,1,2],5] := .5 * T[add1,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),0,0,2],5] := .5 * T[add0,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),0,0,1],6] := .5 * T[add0,6] + .5 * T[add1,6]; T[[op(add1[1..nops(add1)-2]),0,1,2],6] := .5 * T[add1,6] + .5 * T[add2,6]; T[[op(add1[1..nops(add1)-2]),0,0,2],6] := .5 * T[add0,6] + .5 * T[add2,6]; gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2]); gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1]); gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2]); fi; if add1[nops(add1)-1..nops(add1)] = [1,2] then T[[op(add1[1..nops(add1)-2]),1,1,2],5] := .5 * T[add0,5] + .5 * T[add1,5]; T[[op(add1[1..nops(add1)-2]),1,0,2],5] := .5 * T[add1,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),1,0,1],5] := .5 * T[add0,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),1,1,2],6] := .5 * T[add0,6] + .5 * T[add1,6]; T[[op(add1[1..nops(add1)-2]),1,0,2],6] := .5 * T[add1,6] + .5 * T[add2,6]; T[[op(add1[1..nops(add1)-2]),1,0,1],6] := .5 * T[add0,6] + .5 * T[add2,6]; gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1]); gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2]); gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2]); fi; if add1[nops(add1)-1..nops(add1)] = [0,2] then T[[op(add1[1..nops(add1)-2]),2,0,2],5] := .5 * T[add0,5] + .5 * T[add1,5]; T[[op(add1[1..nops(add1)-2]),2,0,1],5] := .5 * T[add1,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),2,1,2],5] := .5 * T[add0,5] + .5 * T[add2,5]; T[[op(add1[1..nops(add1)-2]),2,0,2],6] := .5 * T[add0,6] + .5 * T[add1,6]; T[[op(add1[1..nops(add1)-2]),2,0,1],6] := .5 * T[add1,6] + .5 * T[add2,6]; T[[op(add1[1..nops(add1)-2]),2,1,2],6] := .5 * T[add0,6] + .5 * T[add2,6]; gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2]); gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2]); gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1]); fi; fi; end; #This function returns a list of x,y,z coordinates from a function #stored in T. This list can be sent directly into pointplot3d for plotting SG_plot:=proc(T,m,funcno) local i,plot_pts,func; #func:=select(the_points,[indices(T)],5,0,m+1); func:=vertices(m); plot_pts:=[seq([T[func[i],5],T[func[i],6],T[func[i],funcno]],i=1..nops(func))]; plot_pts; end; #This is used to animate a time series of basis values. #It takes the parameters in vals, which is a list of lists, #with each sublist corresponding to a particular spline function #representation. Funcno is just some place in which to store the #actual values of the function. The output is a list of plot #structures, which can then be sent directly to display, and if you #use insequence = true, it will animate for you... make_animation:=proc(T,m_basis,m_func,vals,funcno,harm_funcno) local S,vm,SG_plots,final_points,inds,i; SG_plots:=[]; S:=gasket(m_basis); vm:=vertices(m_basis); make_neighbors(S,m_basis,vm); make_points(T,m_func); for i from 1 to nops(vals) do print(`computing plot number `,i); make_pretty_picture(T,S,m_basis,m_func,vm,funcno,harm_funcno,vals[i]); final_points:=SG_plot(T,m_func,funcno); SG_plots:=[op(SG_plots),pointplot3d(final_points,color=BLUE,symbol=POINT)]; od; SG_plots; end; #This function will take a line of the sierpinski gasket, and #returns, in order, a list of the vertices along that line. It is #used in conjunction with SG_plot_line in order to create line #plots of a function along a line of the SG. The line is defined #by starting at startpt. From there, we go to first_neighbor. After #that, we continue following rest_neighbor, until the ending (ie, last #2 digits) change, thereby indicating we have gotten to the end of the line. #for example, the line from [0,2] to [1,2] would be gotten by #line_vertices(4,[0,2],3,1), since we would start by taking the 3rd neighbor, #and then continue along the 1st neighbor (the endings would all be 0,1) #until we got to 1,2, at which point, the ending of rest_neighbor would be #1,2, indicating that we have gotten to the end of the line. line_vertices:=proc(m_func,startpt,first_neighbor,rest_neighbor) local lineT,vm,outlist,ending,last_pt; lineT:=gasket(m_func); vm:=vertices(m_func); make_neighbors(lineT,m_func,vm); outlist:=[startpt,lineT[startpt,first_neighbor]]; ending:=last_digits(lineT[startpt,first_neighbor]); last_pt:=lineT[startpt,first_neighbor]; while last_digits(lineT[last_pt,rest_neighbor])=ending do outlist:=[op(outlist),lineT[last_pt,rest_neighbor]]; last_pt:=lineT[last_pt,rest_neighbor]; od; outlist:=[op(outlist),lineT[last_pt,rest_neighbor]]; outlist; end; #This function is just like make_animation, except that it makes 2d line plots #of the function along a line defined by startpt, first_neighbor, rest_neighbor. #You also need to input linelength to tell the program how long the line should #be. For instance, the line from [0,2] to [1,2] would require linelength to be #1/2. For more information about how startpt, first_neighbor, and rest_neighbor #define a line, look at line_vertices above. make_animation_line:=proc(T,m_basis,m_func,vals,funcno,harm_funcno,startpt,first_neighbor,rest_neighbor,linelength) local linepts,S,vm,SG_plots,final_points,inds,i; SG_plots:=[]; S:=gasket(m_basis); vm:=vertices(m_basis); make_neighbors(S,m_basis,vm); make_points(T,m_func); linepts:=line_vertices(m_func,startpt,first_neighbor,rest_neighbor); for i from 1 to nops(vals) do print(`computing plot number `,i); make_pretty_picture(T,S,m_basis,m_func,vm,funcno,harm_funcno,vals[i]); final_points:=SG_plot_line(T,m_func,funcno,linepts,linelength); SG_plots:=[op(SG_plots),plot(final_points,color=BLUE,symbol=POINT)]; od; SG_plots; end; #This function is just like SG_plot, except that it creates a 2d #line plot along a line which is defined by linepts, and is of length #linelength. SG_plot_line:=proc(T,m,funcno,linepts,linelength) local i,plot_pts,numpts; numpts:=nops(linepts); plot_pts:=[seq([(i-1)*linelength/(numpts-1),T[linepts[i],funcno]],i=1..numpts)]; plot_pts; end; #this is just a function to close up extraneous files... close_files:=proc(maxk) local S,vm,EM,IP,k,outfile,i,numops,ilist,jlist,elist,iplist,filename_EM,filename_IP, filename_EMi,filename_EMj,filename_EMx,filename_IPi,filename_IPj,filename_IPx,fEMx,fEMi,fEMj,fIPi,fIPj,fIPx,directory,inds; directory:="/tahoe/home/raj/"; for k from 1 to maxk do filename_EMi:=cat(directory,"EM",convert(k,string),"_i.dat"); filename_EMj:=cat(directory,"EM",convert(k,string),"_j.dat"); filename_EMx:=cat(directory,"EM",convert(k,string),"_x.dat"); filename_IPi:=cat(directory,"IP",convert(k,string),"_i.dat"); filename_IPj:=cat(directory,"IP",convert(k,string),"_j.dat"); filename_IPx:=cat(directory,"IP",convert(k,string),"_x.dat"); close(filename_EMi); close(filename_EMj); close(filename_EMx); close(filename_IPi); close(filename_IPj); close(filename_IPx); od; end; #This function generates EM and IP files which can be read as sparse matrices #into matlab. It generates files like EM1_i, EM2_j,IP4_x. The EM and the IP are #stand for energy and inner product (who woulda thought...) and the i, j, and x #stand for the ith row, jth column, and value at that coordinate. These #3 files can be read into matlab as 3 column vectors, which can then be used #with the sparse command to generate a sparse matrix. The number after the #EM and the IP is the level of m used with a biharmonic basis. This function #is really not very robust, and should be used with care. Note also that the #output files may have to have their unix permissions changed with chmod. Ask #a computer geek to help you out with that if you don't know what I'm talking #about... generate_EM_IP_files:=proc(maxk) local S,vm,EM,IP,k,outfile,i,numops,ilist,jlist,elist,iplist,filename_EM,filename_IP, filename_EMi,filename_EMj,filename_EMx,filename_IPi,filename_IPj,filename_IPx,fEMx,fEMi,fEMj,fIPi,fIPj,fIPx,directory,inds; for k from 1 to maxk do print(k); S:=gasket(k): eval(S): vm:=vertices(k): make_neighbors(S,k,vm): eval(S): print(`calculating matrices...`); EM:=eval(gen_energy_matrix(S,k,vm)); print(`calculated EM...`); IP:=eval(gen_iprod_matrix(S,k,vm)); print(`calculated IP...`); numops:=nops([indices(EM)]); #ilist:=[seq([indices(EM)][i][1],i=1..numops)]; #jlist:=[seq([indices(EM)][i][2],i=1..numops)]; #elist:=evalf([seq(EM[[indices(EM)][i][1],[indices(EM)][i][2]],i=1..numops)],20); #iplist:=evalf([seq(IP[[indices(IP)][i][1],[indices(IP)][i][2]],i=1..numops)],20); filename_EM:=cat("EM" , convert(k,string)); filename_IP:=cat("IP" , convert(k,string)); directory:="/tahoe/home/raj/"; print(`saving .m files...`); save EM,cat(directory,filename_EM,".m"); save IP,cat(directory,filename_IP,".m"); filename_EMi:=cat(directory,"EM",convert(k,string),"_i.dat"); filename_EMj:=cat(directory,"EM",convert(k,string),"_j.dat"); filename_EMx:=cat(directory,"EM",convert(k,string),"_x.dat"); filename_IPi:=cat(directory,"IP",convert(k,string),"_i.dat"); filename_IPj:=cat(directory,"IP",convert(k,string),"_j.dat"); filename_IPx:=cat(directory,"IP",convert(k,string),"_x.dat"); #close(filename_EMi); #close(filename_EMj); #close(filename_EMx); #close(filename_IPi); #close(filename_IPj); #close(filename_IPx); fEMi:=open(filename_EMi,WRITE); fEMj:=open(filename_EMj,WRITE); fEMx:=open(filename_EMx,WRITE); fIPi:=open(filename_IPi,WRITE); fIPj:=open(filename_IPj,WRITE); fIPx:=open(filename_IPx,WRITE); print(`copying to text files...`); inds:=[indices(EM)]; print(`(got indices...)`); for i from 1 to numops do # fprintf(fEMi,`%20.20g\n`,ilist[i]); # fprintf(fEMj,`%20.20g\n`,jlist[i]); # fprintf(fEMx,`%20.20g\n`,elist[i]); # fprintf(fIPi,`%20.20g\n`,ilist[i]); # fprintf(fIPj,`%20.20g\n`,jlist[i]); # fprintf(fIPx,`%20.20g\n`,iplist[i]); fprintf(fEMi,`%20.20g\n`,inds[i][1]); fprintf(fEMj,`%20.20g\n`,inds[i][2]); fprintf(fEMx,`%20.20g\n`,EM[inds[i][1],inds[i][2]]); fprintf(fIPi,`%20.20g\n`,inds[i][1]); fprintf(fIPj,`%20.20g\n`,inds[i][2]); fprintf(fIPx,`%20.20g\n`,IP[inds[i][1],inds[i][2]]); od; close(fEMi); close(fEMj); close(fEMx); close(fIPi); close(fIPj); close(fIPx); od; end; #generate_EM_IP_files(6); #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; 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; #This function returns a table R which contains the values of p and q as defined #by Mike Usher for determining the values of n-harmonic extensions. The values #a and b are values of inner products. The function basically implements the #formula for determining extensions based on an iterative method involving inner #products. See "Splines on Fractals" from more information (ie, don't ask me, #I just work here...). ps_n_qs:=proc(n) local i,l,j,R,m; R:=table(); R[p,0]:=2/5; R[q,0]:=1/5; R[a,0]:=2; R[b,0]:=-1; R[p,1]:=-3/25; R[q,1]:=-7/75; R[a,1]:=7/45; R[b,1]:=4/45; for i from 1 to n-1 do R[a,i+1]:=(2/15)*(1/(5^(2*i)-6*5^(i-1) + 1/5)) * sum( (((4*5^i-4/15)*R[p,i+1-j]+(5^i+13/15)*R[q,i+1-j])*(R[a,j]+2*R[b,j])), j=1..i); R[b,i+1]:=(2/15)*(1/(5^(2*i)-6*5^(i-1) + 1/5))*sum( (((3*5^i-13/15) * R[p,i+1-j] + (2*5^i-14/15) * R[q,i+1-j]) * (R[a,j] + 2*R[b,j])),j=1..i); R[p,i+1]:=(-1/5)*((9/5)*(R[a,i+1] + 2*R[b,i+1]) + sum( ((4*R[p,i+1-m] + R[q,i+1-m])*R[a,m] + (3 * R[p,i+1-m] + 2*R[q,i+1-m]) * R[b,m]),m=1..i)); R[q,i+1]:=(-1/5)*((7/5)*(R[a,i+1] + 2*R[b,i+1]) + sum( ((2*R[p,i+1-m] + 3*R[q,i+1-m])*R[a,m] + (4 * R[p,i+1-m] + R[q,i+1-m]) * R[b,m]),m=1..i)); od; R; end; #As long as I'm here, why don't we just calculate this thing out to 5 #so that we can do sex-harmonic functions (Yeah, baby!) R:=ps_n_qs(5); eval(R): #This function creates a n-harmonic function, using the values for harmonic extensions calculated #in R. It uses the recursive method of calculating the n-harmonic extension based on Mike #Usher's work. In order to understand how the recursion addressing works, I recommend you look #at the documentation of make_points and gen_points in SG.mws, where I fully explain how the #recursion works (in terms of the addressing). Basically, after that, this function #is quite straightforward; it simply calculates the next level based on the values of the previous #level, as well as the necessary laplacians. This procedure is called by the wrapper function #n_harm. nharm:=proc(T,R,n,m,add0,add1,add2,funcno) local i,j; if m <> 0 then #print(funcno); #print(add0,add1,add2); if add1[nops(add1)-1..nops(add1)] = [0,1] then #This loop is because we need to calculate up to n-1 laplacians of the function #at this level as well... for i from 0 to n-1 do #the sums here determine the value at the next level based on the values of #various laplacians and the such before. Again, see Mike Usher's work for a #derivation of this formula. T[[op(add1[1..nops(add1)-2]),0,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),0,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),0,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i); od; nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno); nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1],funcno); nharm(T,R,n,m-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 for i from 0 to n-1 do T[[op(add1[1..nops(add1)-2]),1,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),1,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),1,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i); od; nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno); nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],funcno); nharm(T,R,n,m-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 for i from 0 to n-1 do T[[op(add1[1..nops(add1)-2]),2,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),2,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i); T[[op(add1[1..nops(add1)-2]),2,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i); od; nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2],funcno); nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno); nharm(T,R,n,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1],funcno); fi; fi; end; #So this is basically a wrapper function around the function nharm. This function is #necessary because nharm won't work until V1, due to idiosyncracies in the labeling system. #This function basically computes the first n-harmonic extension, and then sends the rest out to #nharm to compute. n stands for n-harmonic, and R is the table of harmonic extension values #calculated using ps_and_qs. In order to make this function work, you have to put in the values #of the function at the boundaries in funcno, as well as the values of the laplacians of the #function in funcno+1, L^2 of the function in funcno+2... and so forth. This uniquely determines #an n-harmonic function (ie, you need to specify boundary values of L^(0->n-1)). Notice that #the algorithm developed required that the function as well as all of it's laplacians be #defined out to all levels. This has the added bonus that while calculating the function, you have #calculated the laplacian of it in funcno+1, L^2(f) in funcno+2... Fun, isn't it? n_harm:=proc(T,R,n,m,add0,add1,add2,funcno) local i,j; for i from 0 to n-1 do T[[0,1],funcno+n-1-i]:=sum((1/(5^(j*nops(add1))))*(R[p,j]*(T[add0,funcno+n-1-i+j]+T[add1,funcno+n-1-i+j])+R[q,j]*T[add2,funcno+n-1-i+j]),j=0..i); T[[1,2],funcno+n-1-i]:=sum((1/(5^(j*1)))*(R[p,j]*(T[add1,funcno+n-1-i+j]+T[add2,funcno+n-1-i+j])+R[q,j]*T[add0,funcno+n-1-i+j]),j=0..i); T[[0,2],funcno+n-1-i]:=sum((1/(5^(j*1)))*(R[p,j]*(T[add0,funcno+n-1-i+j]+T[add2,funcno+n-1-i+j])+R[q,j]*T[add1,funcno+n-1-i+j]),j=0..i); od; nharm(T,R,n,m-1,[0],[0,1],[0,2],funcno); nharm(T,R,n,m-1,[1],[1,2],[0,1],funcno); nharm(T,R,n,m-1,[2],[0,2],[1,2],funcno); end; #This function recursively generates a harmonic function. It is a much simpler algorithm than #the n-harmonic algorithm, so the addressing business is much more transparent here. Basically, #it implements the 2/5, 2/5, 1/5 rule for generating the harmonic extension, and recursively #generates values all the way down into the gasket. The recursion is a little bit complicated, #because each step of the way, you need to know the orientation of the triangle in order to know #what type of points are assigned what values. Basically, we used a rotational method of #determining what the orientation is: we made it so that the 0th address always was a vertex of the #previous triangle (ie, the 0th address of the first triangle is [1]...). This way, we always #which address to assign to, based on the last digits of address 1 (if address 1 is [1,2], then #we know the triangle is a 1 triangle. Thus, you know what addresses to assign what values. #Anyway, it works. If you want a somewhat more transparent example, see make_points and gen_points #in SG.mws. The recursion method used there is exactly the same as used here, and should be easily #generalizable to other fractals (ie, pentagasket...). harm:=proc(T,m,add0,add1,add2,funcno) if m <> 0 then #print(funcno); #print(add0,add1,add2); if add1[nops(add1)-1..nops(add1)] = [0,1] then T[[op(add1[1..nops(add1)-2]),0,0,1],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno]; T[[op(add1[1..nops(add1)-2]),0,1,2],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno]; T[[op(add1[1..nops(add1)-2]),0,0,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno]; harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno); harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1],funcno); harm(T,m-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 T[[op(add1[1..nops(add1)-2]),1,1,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno]; T[[op(add1[1..nops(add1)-2]),1,0,2],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno]; T[[op(add1[1..nops(add1)-2]),1,0,1],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno]; harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno); harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],funcno); harm(T,m-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 T[[op(add1[1..nops(add1)-2]),2,0,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno]; T[[op(add1[1..nops(add1)-2]),2,0,1],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno]; T[[op(add1[1..nops(add1)-2]),2,1,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno]; harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2],funcno); harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno); harm(T,m-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 is basically a wrapper for harm #It generates a harmonic function with val0,val1,val2 on the boundaries. #This function calculates the first harmonic extension, and then calls harm #on the 3 parts to generate the rest of the harmonic function. It can't do it directly due #to idiosyncracies in the nomenclature of the V0 points. harmonic:=proc(T,m,val0,val1,val2,funcno) T[[0],funcno]:=val0; T[[1],funcno]:=val1; T[[2],funcno]:=val2; T[[0,1],funcno]:=2/5 * val0 + 2/5 * val1 + 1/5 * val2; T[[0,2],funcno]:=2/5 * val0 + 2/5 * val2 + 1/5 * val1; T[[1,2],funcno]:=2/5 * val1 + 2/5 * val2 + 1/5 * val0; 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 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; #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; #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; #we are going to make a bunch of functions to be used as f's and q's... #5 -> x #6 -> y #7 -> constant = 1 #8 -> harm(5,0,0) #9 -> harm(.1,0,0) #10 -> harm(100,0,0) #11 -> (harm(5,0,0))^2 #12 -> quad-harm(-100,40,-130) #16 -> tri-harm(-100,40,-130) # function number 7 should be plugged in for funcno... make_f_and_q:=proc(T,m,funcno) make_points(T,m); harmonic(T,m,1,1,1,funcno); harmonic(T,m,5,0,0,funcno + 1); harmonic(T,m,.1,0,0,funcno + 2); harmonic(T,m,100,0,0,funcno + 3); mult(T,m,funcno,funcno,funcno + 4); T[[0],funcno+4]:=0: T[[0],funcno+5]:=0: T[[0],funcno+6]:=0: T[[0],funcno+7]:=-100: T[[1],funcno+4]:=0: T[[1],funcno+5]:=0: T[[1],funcno+6]:=0: T[[1],funcno+7]:=40: T[[2],funcno+4]:=0: T[[2],funcno+5]:=0: T[[2],funcno+6]:=0: T[[2],funcno+7]:=-130: n_harm(T,R,4,m,[0],[1],[2],funcno+5); end;