{VERSION 5 0 "IBM INTEL LINUX" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 0 0 0 0 0 1 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Outpu t" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 } 1 3 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 232 "#this procedure generates all words of leng th m in L letters...\n\np:=proc(L,m)\noption remember;\nlocal L1;\n i f m = 1 then L\n else\n L1 := p(L,m-1);\n [op(map(add_one,L1,0) ),op(map(add_one,L1,1)),op(map(add_one,L1,2))]\n fi;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%\"pGf*6$%\"LG%\"mG6#%#L1G6#%)rememberG6\"@ %/9%\"\"\"9$C$>8$-F$6$F2,&F0F1F1!\"\"7%-%#opG6#-%$mapG6%%(add_oneGF5\" \"!-F<6#-F?6%FAF5F1-F<6#-F?6%FAF5\"\"#F-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "add_one:=(u,d)-> [op(u),d];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(add_oneGf*6$%\"uG%\"dG6\"6$%)operatorG%&arrowGF)7$-% #opG6#9$9%F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1237 "#This function returns all vertices in Vm in a big list\nvertices:= proc(m) \n\nlocal i,outlist;\noption remember;\n\nif m < 0 then RETURN([]) fi; \nif m = 0 then RETURN([[0],[1],[2]]) fi;\nif m = 1 then RETURN([[0],[ 1],[2],[0,1],[0,2],[1,2]]) fi;\nif m=2 then \n[[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 \n outlist:=[[0] , [1] , [2] , [0,1], [0,2], [1,2]]; \n for i from 2 to m do\n 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] ))];\n od;\n outlist; \n\n \n#[[0] , [1] , [2] , [0,1], [0,2], \+ [1,2], op(map(add_pair, p([[0],[1],[2]], (m-1)), [0,1])), op(map(add_p air, p([[0],[1],[2]], (m-1)), [0,2])), op(map(add_pair, p([[0],[ 1],[2]], (m-1)), [1,2]))];\nfi;\nend;\n\n#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;\nend ;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)verticesGf*6#%\"mG6$%\"iG%(out listG6#%)rememberG6\"C&@$29$\"\"!-%'RETURNG6#7\"@$/F1F2-F46#7%7#F27#\" \"\"7#\"\"#@$/F1F>-F46#7(F7$F2F@7$F>F@@%/F1F@71F7%F2F2F@7%F2F>F@7%F>F2F>7%F>F2F@7%F>F>F@7%F@F2F>7%F@F2F@7%F@F> F@C%>8%FE?(8$F@F>F1%%trueG>FW7&-%#opG6#FW-Fhn6#-%$mapG6%%)add_pairG-% \"pG6$F;,&FYF>F>!\"\"FF-Fhn6#-F]o6%F_oF`oFG-Fhn6#-F]o6%F_oF`oFHFWF-F-F -" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+vertices_mGf*6#%\"mG6\"6#%)rem emberGF(C%@$/9$\"\"!-%'RETURNG6#7%7#F/7#\"\"\"7#\"\"#@$/F.F6-F16#7%7$F /F67$F/F87$F6F8@$2F6F.7%-%#opG6#-%$mapG6%%)add_pairG-%\"pG6$F3,&F.F6F6 !\"\"F>-FE6#-FH6%FJFKF?-FE6#-FH6%FJFKF@F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "add_pair := (u,v) -> [op(u), op(v)];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%)add_pairGf*6$%\"uG%\"vG6\"6$%)operatorG%& arrowGF)7$-%#opG6#9$-F/6#9%F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 220 "#This function returns a table of the gasket out to \+ level m...\ngasket:=proc(m)\n local i, vm, T;\n\n vm:= vertices( m);\n T:=table();\n for i from 1 to nops(vm) do \n T[vm[i ], 0] := vm[i] od;\n\n\n T; \nend;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'gasketGf*6#%\"mG6%%\"iG%#vmG%\"TG6\"F,C&>8%-%)verticesG6#9$>8 &-%&tableGF,?(8$\"\"\"F:-%%nopsG6#F/%%trueG>&F56$&F/6#F9\"\"!FBF5F,F,F ," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4949 "#This function makes stores the neighbors of each vertex out to level m in\n#the gasket ta ble T. The neighbors are ordered clockwise; i.e., the 1st neighbor\n# of 0,1 in V3 is 1,0,0,1, etc...\nmake_neighbors:=proc(T,m,vm)\n \n l ocal i,j,depth,nb1,nb2,nb3,nb4,iverts,last_verts;\n\n#first, we have t o deal with the special case of V1...\n if m=1 then \n\nprint(`Yo, V 1 is STOOOOOOPID!`);\n T[[0],1]:=[0,1];\n T[[0],2]:=[0,2];\n \n T[ [1],1]:=[1,2];\n T[[1],2]:=[0,1]; \n\n T[[2],1]:=[0,2];\n T[[2],2] :=[1,2];\n\n T[[0,1],3]:=[0,2];\n T[[0,1],4]:=[0];\n T[[0,1],1]:=[1 ];\n T[[0,1],2]:=[1,2];\n\n T[[1,2],3]:=[0,1];\n T[[1,2],4]:=[1];\n T[[1,2],2]:=[0,2];\n T[[1,2],1]:=[2];\n\n T[[0,2],1]:=[0];\n T[[0 ,2],2]:=[0,1];\n T[[0,2],3]:=[1,2];\n T[[0,2],4]:=[2];\n\n\n#hard co ding reverse neighbors...\n T[T[[0],1],4] := [0]; \n T[T[[0],2],1] \+ := [0];\n\n T[T[[1],1],4] := [1];\n T[T[[1],2],1] := [1];\n\n T[T[[ 2],1],4] := [2];\n T[T[[2],2],1] := [2];\n\n#done with hard coding (y ay!)\n fi;\n\n if m <>1 then #but if we are not in V1...\n#hard codi ng the neighbors for V0,V1\n\n T[[0],1]:=[seq(0,i=1..m),1];\n T[[0], 2]:=[seq(0,i=1..m),2];\n \n T[[1],1]:=[seq(1,i=1..m),2];\n T[[1],2] :=[seq(1,i=1..m-1),0,1]; \n\n T[[2],1]:=[seq(2,i=1..m-1),0,2];\n T[ [2],2]:=[seq(2,i=1..m-1),1,2];\n\n T[[0,1],3]:=[0,seq(1,i=1..m-2),1,2 ];\n T[[0,1],4]:=[0,seq(1,i=1..m-2),0,1];\n T[[0,1],1]:=[1,seq(0,i=1 ..m-2),0,1];\n T[[0,1],2]:=[1,seq(0,i=1..m-2),0,2];\n\n T[[1,2],3]:= [1,seq(2,i=1..m-2),0,2];\n T[[1,2],4]:=[1,seq(2,i=1..m-2),1,2];\n T[ [1,2],2]:=[2,seq(1,i=1..m-2),0,1];\n T[[1,2],1]:=[2,seq(1,i=1..m-2),1 ,2];\n\n T[[0,2],1]:=[0,seq(2,i=1..m-2),0,2];\n T[[0,2],2]:=[0,seq(2 ,i=1..m-2),1,2];\n T[[0,2],3]:=[2,seq(0,i=1..m-2),0,1];\n T[[0,2],4] :=[2,seq(0,i=1..m-2),0,2];\n\n\n#hard coding reverse neighbors...\n T [T[[0],1],4] := [0]; \n T[T[[0],2],1] := [0];\n\n T[T[[1],1],4] := \+ [1];\n T[T[[1],2],1] := [1];\n\n T[T[[2],1],4] := [2];\n T[T[[2],2] ,1] := [2];\n\n \n T[T[[0,1],1],4] := [0,1];\n T[T[[0,1],2],1] := [ 0,1];\n T[T[[0,1],3],4] := [0,1];\n T[T[[0,1],4],1] := [0,1];\n\n T [T[[1,2],1],4] := [1,2];\n T[T[[1,2],2],1] := [1,2];\n T[T[[1,2],3], 4] := [1,2];\n T[T[[1,2],4],1] := [1,2];\n\n T[T[[0,2],1],4] := [0,2 ];\n T[T[[0,2],2],1] := [0,2];\n T[T[[0,2],3],4] := [0,2];\n T[T[[0 ,2],4],1] := [0,2];\n\n#done with hard coding (yay!)\n \n\n# First, w e compute the neighbors of all pts in Vk, 1%/make_neighborsGf*6%%\"TG%\"mG%#vm G6+%\"iG%\"jG%&depthG%$nb1G%$nb2G%$nb3G%$nb4G%'ivertsG%+last_vertsG6\" F4C%@$/9%\"\"\"C;-%&printG6#%7Yo,~V1~is~STOOOOOOPID!G>&9$6$7#\"\"!F97$ FDF9>&FA6$FC\"\"#7$FDFI>&FA6$7#F9F97$F9FI>&FA6$FNFIFE>&FA6$7#FIF9FJ>&F A6$FVFIFO>&FA6$FE\"\"$FJ>&FA6$FE\"\"%FC>&FA6$FEF9FN>&FA6$FEFIFO>&FA6$F OFgnFE>&FA6$FOF[oFN>&FA6$FOFIFJ>&FA6$FOF9FV>&FA6$FJF9FC>&FA6$FJFIFE>&F A6$FJFgnFO>&FA6$FJF[oFV>&FA6$F@F[oFC>&FA6$FGF9FC>&FA6$FLF[oFN>&FA6$FQF 9FN>&FA6$FTF[oFV>&FA6$FXF9FV@$0F8F9CM>F@7$-%$seqG6$FD/8$;F9F8F9>FG7$Fa rFI>FL7$-Fbr6$F9FdrFI>FQ7%-Fbr6$F9/Fer;F9,&F8F9F9!\"\"FDF9>FT7%-Fbr6$F IFasFDFI>FX7%FgsF9FI>Fen7&FD-Fbr6$F9/Fer;F9,&F8F9FIFdsF9FI>Fin7&FDF]tF DF9>F]o7&F9-Fbr6$FDF_tFDF9>F`o7&F9FftFDFI>Fco7&F9-Fbr6$FIF_tFDFI>Ffo7& F9F\\uF9FI>Fio7&FIF]tFDF9>F\\p7&FIF]tF9FI>F_p7&FDF\\uFDFI>Fbp7&FDF\\uF 9FI>Fep7&FIFftFDF9>Fhp7&FIFftFDFI>F[qFC>F^qFC>FaqFN>FdqFN>FgqFV>FjqFV> &FA6$F]oF[oFE>&FA6$F`oF9FE>&FA6$FenF[oFE>&FA6$FinF9FE>&FA6$F\\pF[oFO>& FA6$FioF9FO>&FA6$FcoF[oFO>&FA6$FfoF9FO>&FA6$F_pF[oFJ>&FA6$FbpF9FJ>&FA6 $FepF[oFJ>&FA6$FhpF9FJ?(FerFIF9Fcs%%trueGC(>8+-%0pick_vertices01G6$9&F er?(8%F9F9-%%nopsG6#FjxFgxC.>8'7&-%#opG6#&&Fjx6#F`y6#;F9,&FerF9F9Fds-F iy6#&F]o6#;F9,&F8F9FerFdsFDF9>8(7&Fhy-Fiy6#&F`oFdzFDFI>8)7&Fhy-Fiy6#&F enFdzF9FI>8*7&Fhy-Fiy6#&FinFdzFDF9>&FA6$F\\zF9Ffy>&FA6$F\\zFIFhz>&FA6$ F\\zFgnF^[l>&FA6$F\\zF[oFd[l>&FA6$FfyF[oF\\z>&FA6$FhzF9F\\z>&FA6$F^[lF [oF\\z>&FA6$Fd[lF9F\\z>Fjx-%0pick_vertices12GF]y?(F`yF9F9FayFgxC.>Ffy7 &Fhy-Fiy6#&F\\pFdzF9FI>Fhz7&Fhy-Fiy6#&FioFdzFDF9>F^[l7&Fhy-Fiy6#&FcoFd zFDFI>Fd[l7&Fhy-Fiy6#&FfoFdzF9FI>Fj[lFfy>F]\\lFhz>F`\\lF^[l>Fc\\lFd[l> Ff\\lF\\z>Fi\\lF\\z>F\\]lF\\z>F_]lF\\z>Fjx-%0pick_vertices02GF]y?(F`yF 9F9FayFgxC.>Ffy7&Fhy-Fiy6#&F_pFdzFDFI>Fhz7&Fhy-Fiy6#&FbpFdzF9FI>F^[l7& Fhy-Fiy6#&FepFdzFDF9>Fd[l7&Fhy-Fiy6#&FhpFdzFDFI>Fj[lFfy>F]\\lFhz>F`\\l F^[l>Fc\\lFd[l>Ff\\lF\\z>Fi\\lF\\z>F\\]lF\\z>F_]lF\\z>8,-F\\y6$F^yF8?( F`yF9F9-Fby6#FdalFgxC&>Fhz7%-Fiy6#&&FdalF]z6#FbsF9FI>F^[l7%F]blFDFI>&F A6$F`blFIFhz>&FA6$F`blFgnF^[l>Fdal-Fc]lFfal?(F`yF9F9FhalFgxC&>FhzFcbl> F^[l7%F]blFDF9>FeblFhz>FhblF^[l>Fdal-Fd_lFfal?(F`yF9F9FhalFgxC&>FhzF`c l>F^[lF\\bl>FeblFhz>FhblF^[lFAF4F4F4" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "#picks vertices out of vm of length i with ending [0 ,1]\npick_vertices01:=proc(vm,i)\n select(len_list_eq, vm,i,[0,1]);\n end;" }{TEXT -1 1 "\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%0pick_verti ces01Gf*6$%#vmG%\"iG6\"F)F)-%'selectG6&%,len_list_eqG9$9%7$\"\"!\"\"\" F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 93 "len_list_eq:=proc (l,i,ending)\n evalb(nops(l) = i+1 and l[nops(l)-1..nops(l)] = ending );\nend;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,len_list_eqGf*6%%\"lG% \"iG%'endingG6\"F*F*-%&evalbG6#3/-%%nopsG6#9$,&9%\"\"\"F6F6/&F36#;,&F0 F6F6!\"\"F09&F*F*F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "#lik e pick_vertices01...\npick_vertices12:=proc(vm,i)\n select(len_list_e q, vm,i,[1,2]);\nend;\n\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%0pick_v ertices12Gf*6$%#vmG%\"iG6\"F)F)-%'selectG6&%,len_list_eqG9$9%7$\"\"\" \"\"#F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "#like pick_v ertices01...\npick_vertices02:=proc(vm,i)\n select(len_list_eq, vm,i, [0,2]);\nend;\n\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%0pick_vertices0 2Gf*6$%#vmG%\"iG6\"F)F)-%'selectG6&%,len_list_eqG9$9%7$\"\"!\"\"#F)F)F )" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 675 "#This function stores the x and y coordinates of each vertex\n#in the 5 and 6 slots of T, r espectively\nmake_points:=proc(T,m)\n local i;\n\n #func:=select(the _points,[indices(T)],0,3);\n\n\n# first, we have to manually set the v alues of x and y out to V1...\n T[[0],5]:=0;\n T[[1],5]:=.5;\n T[[2 ],5]:=-.5;\n\n T[[0],6]:=sqrt(3)/2;\n T[[1],6]:=0;\n T[[2],6]:=0;\n \n T[[0,1],5]:=.25;\n T[[0,1],6]:=sqrt(3)/4;\n T[[1,2],5]:=0;\n T[ [1,2],6]:=0;\n T[[0,2],5]:=-.25;\n T[[0,2],6]:=sqrt(3)/4; \n\n#gen_p oints recursively creates the rest of the coordinates out to Vm... \n \ngen_points(T,m-1,[0],[0,1],[0,2]);\ngen_points(T,m-1,[1],[1,2],[0,1] );\ngen_points(T,m-1,[2],[0,2],[1,2]); \n\nend;" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%,make_pointsGf*6$%\"TG%\"mG6#%\"iG6\"F+C1>&9$6 $7#\"\"!\"\"&F2>&F/6$7#\"\"\"F3$F3!\"\">&F/6$7#\"\"#F3$!\"&F:>&F/6$F1 \"\"',$-%%sqrtG6#\"\"$#F8F?>&F/6$F7FEF2>&F/6$F>FEF2>&F/6$7$F2F8F3$\"#D !\"#>&F/6$FUFE,$FG#F8\"\"%>&F/6$7$F8F?F3F2>&F/6$F\\oFEF2>&F/6$7$F2F?F3 $!#DFX>&F/6$FcoFEFfn-%+gen_pointsG6'F/,&9%F8F8F:F1FUFco-Fjo6'F/F\\pF7F \\oFU-Fjo6'F/F\\pF>FcoF\\oF+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 119 "the_points:=proc(L,funcno,minlen,maxlen)\n evalb(L[ 2] = funcno and nops(L[1]) >= minlen and nops(L[1]) <= maxlen);\nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+the_pointsGf*6&%\"LG%'funcnoG%' minlenG%'maxlenG6\"F+F+-%&evalbG6#33/&9$6#\"\"#9%19&-%%nopsG6#&F36#\" \"\"1F99'F+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3930 "#This \+ function is called by make_points. It calls itself recursively 'n stu ff.\ngen_points:=proc(T,m,add0,add1,add2)\n if m <> 0 then\n\n#So th e general strategy is as follows: We can't assign values to the next \+ level of \n#vertices without knowing the orientation of the triangle ( because we need to know whether we need to assign\n#to a 0,1 point, 1, 2 point, etc.). So what we do is call the function with 3 addresses\n #and make it so that the 0th address (add0) is one of the vertices of \+ the outside triangle\n#(i.e., going from the V0 triangle, the triangle s would be T1=([0],[0,1],[0,2]),T2=([1],[1,2],[0,1])\n#and so forth). \+ Thus, by looking at add1 of this triangle, we can determine the orien tation of\n#the triangle by checking whether add1 ends in 0,1, 0,2 or \+ 1,2. Then we assign the values\n#accordingly, after which we call gen _points recursively, making sure to keep the orientation\n#of the 3 ad dresses in the recursive calls. We use this same method to generate h armonic functions\n#and n-harmonic functions, as well as several other things. Notice that this recursive process\n#has it's base case as V 1, therefore requiring 3 seperate calls to gen_points from make_points .\n#This is due to the idiosyncratic addressing of the vertices of V0. Oh well....\n\n if add1[nops(add1)-1..nops(add1)] = [0,1] then \n T[[op(add1[1..nops(add1)-2]),0,0,1],5] := .5 * T[add0,5] + .5 * T[add1,5]; \n T[[op(add1[1..nops(add1)-2]),0,1,2],5] \+ := .5 * T[add1,5] + .5 * T[add2,5];\n T[[op(add1[1..nops(add1 )-2]),0,0,2],5] := .5 * T[add0,5] + .5 * T[add2,5]; \n T[[op( add1[1..nops(add1)-2]),0,0,1],6] := .5 * T[add0,6] + .5 * T[add1,6]; \+ \n T[[op(add1[1..nops(add1)-2]),0,1,2],6] := .5 * T[add1,6] + .5 * T[add2,6];\n T[[op(add1[1..nops(add1)-2]),0,0,2],6] := \+ .5 * T[add0,6] + .5 * T[add2,6];\n gen_points(T,m-1,add0,[op( add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2]);\n \+ gen_points(T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add 1[1..nops(add1)-2]),0,0,1]);\n gen_points(T,m-1,add2,[op(add1 [1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2]);\n f i;\n if add1[nops(add1)-1..nops(add1)] = [1,2] then\n T[ [op(add1[1..nops(add1)-2]),1,1,2],5] := .5 * T[add0,5] + .5 * T[add1,5 ]; \n T[[op(add1[1..nops(add1)-2]),1,0,2],5] := .5 * T[add1, 5] + .5 * T[add2,5];\n T[[op(add1[1..nops(add1)-2]),1,0,1],5] := .5 * T[add0,5] + .5 * T[add2,5];\n T[[op(add1[1..nops(add 1)-2]),1,1,2],6] := .5 * T[add0,6] + .5 * T[add1,6]; \n T[[o p(add1[1..nops(add1)-2]),1,0,2],6] := .5 * T[add1,6] + .5 * T[add2,6]; \n T[[op(add1[1..nops(add1)-2]),1,0,1],6] := .5 * T[add0,6] + .5 * T[add2,6];\n gen_points(T,m-1,add0,[op(add1[1..nops(add 1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1]);\n gen_point s(T,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)- 2]),1,1,2]);\n gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2 ]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2]);\n fi;\n if add 1[nops(add1)-1..nops(add1)] = [0,2] then\n T[[op(add1[1..nops (add1)-2]),2,0,2],5] := .5 * T[add0,5] + .5 * T[add1,5]; \n \+ T[[op(add1[1..nops(add1)-2]),2,0,1],5] := .5 * T[add1,5] + .5 * T[add2 ,5];\n T[[op(add1[1..nops(add1)-2]),2,1,2],5] := .5 * T[add0, 5] + .5 * T[add2,5];\n T[[op(add1[1..nops(add1)-2]),2,0,2],6] := .5 * T[add0,6] + .5 * T[add1,6]; \n T[[op(add1[1..nops(a dd1)-2]),2,0,1],6] := .5 * T[add1,6] + .5 * T[add2,6];\n T[[o p(add1[1..nops(add1)-2]),2,1,2],6] := .5 * T[add0,6] + .5 * T[add2,6]; \n gen_points(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[o p(add1[1..nops(add1)-2]),2,1,2]);\n gen_points(T,m-1,add1,[op (add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2]);\n \+ gen_points(T,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(ad d1[1..nops(add1)-2]),2,0,1]);\n fi;\n fi; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+gen_pointsGf*6'%\"TG%\"mG%%add0G%%add1G%%add 2G6\"F,F,@$09%\"\"!C%@$/&9'6#;,&-%%nopsG6#F5\"\"\"F&9 $6$7&-%#opG6#&F56#;F<,&F9F<\"\"#F=F0F0F<\"\"&,&&FB6$9&FM$FMF=*&FRF<&FB 6$F5FMF&FB6$7&FEF0F&FB6$7&FEF0F0FLF M,&FOFR*&FRF&FB6$FD\"\"',&&FB6$FQFboFR*&FRF<&FB6$F5FboF&F B6$FYFbo,&FgoFR*&FRF<&FB6$FhnFboF&FB6$F\\oFbo,&FdoFR*&FRF&FB6$7&FEF&FB6$7&FEF&FB6$7&FEF&FB6$F cqFboFco>&FB6$FgqFboF\\p>&FB6$F[rFboFcp-F$6'FBFgpFQFcqF[r-F$6'FBFgpF5F gqFcq-F$6'FBFgpFhnF[rFgq@$/F47$F0FLC+>&FB6$7&FEFLF0FLFMFN>&FB6$7&FEFLF 0F&FB6$7&FEFLF&FB6$FbsFboFco>&FB6$FfsFboF\\p>&FB6$FjsFb oFcp-F$6'FBFgpFQFbsFjs-F$6'FBFgpF5FfsFbs-F$6'FBFgpFhnFjsFfsF,F,F," }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 373 "#This function returns a li st of x,y,z coordinates from a function\n#stored in T. This list can \+ be sent directly into pointplot3d for plotting\nSG_plot:=proc(T,m,func no)\n local i,plot_pts,func;\n\n #func:=select(the_points,[indices(T )],5,0,m+1);\n func:=vertices(m);\n \n plot_pts:=[seq([T[func[i],5] ,T[func[i],6],T[func[i],funcno]],i=1..nops(func))];\n \n plot_pts;\n end; \n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%(SG_plotGf*6%%\"TG%\"mG %'funcnoG6%%\"iG%)plot_ptsG%%funcG6\"F.C%>8&-%)verticesG6#9%>8%7#-%$se qG6$7%&9$6$&F16#8$\"\"&&F>6$F@\"\"'&F>6$F@9&/FB;\"\"\"-%%nopsG6#F1F7F. F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 932 "#This is used to a nimate a time series of basis values.\n#It takes the parameters in val s, which is a list of lists,\n#with each sublist corresponding to a pa rticular spline function\n#representation. Funcno is just some place \+ in which to store the\n#actual values of the function. The output is \+ a list of plot\n#structures, which can then be sent directly to displa y, and if you\n#use insequence = true, it will animate for you...\n\nm ake_animation:=proc(T,m_basis,m_func,vals,funcno,harm_funcno)\n\nlocal S,vm,SG_plots,final_points,inds,i;\n\nSG_plots:=[];\n\nS:=gasket(m_ba sis);\nvm:=vertices(m_basis);\nmake_neighbors(S,m_basis,vm);\nmake_poi nts(T,m_func);\n\nfor i from 1 to nops(vals) do\n\n print(`computing plot number `,i);\n make_pretty_picture(T,S,m_basis,m_func,vm,funcn o,harm_funcno,vals[i]);\n final_points:=SG_plot(T,m_func,funcno);\n \+ SG_plots:=[op(SG_plots),pointplot3d(final_points,color=BLUE,symbol=P OINT)];\nod;\n\nSG_plots;\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#> %/make_animationGf*6(%\"TG%(m_basisG%'m_funcG%%valsG%'funcnoG%,harm_fu ncnoG6(%\"SG%#vmG%)SG_plotsG%-final_pointsG%%indsG%\"iG6\"F4C)>8&7\">8 $-%'gasketG6#9%>8%-%)verticesGF=-%/make_neighborsG6%F:F>F@-%,make_poin tsG6$9$9&?(8)\"\"\"FM-%%nopsG6#9'%%trueGC&-%&printG6$%7computing~plot~ number~GFL-%4make_pretty_pictureG6*FIF:F>FJF@9(9)&FQ6#FL>8'-%(SG_plotG 6%FIFJFen>F77$-%#opG6#F7-%,pointplot3dG6%Fjn/%&colorG%%BLUEG/%'symbolG %&POINTGF7F4F4F4" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 1388 "#This function will take a line of the sierpinski gasket, and \n#returns, in order, a list of the vertices along that l ine. It is\n#used in conjunction with SG_plot_line in order to create line\n#plots of a function along a line of the SG. The line is defin ed\n#by starting at startpt. From there, we go to first_neighbor. Af ter\n#that, we continue following rest_neighbor, until the ending (ie, last\n#2 digits) change, thereby indicating we have gotten to the end of the line.\n#for example, the line from [0,2] to [1,2] would be got ten by\n#line_vertices(4,[0,2],3,1), since we would start by taking th e 3rd neighbor,\n#and then continue along the 1st neighbor (the ending s would all be 0,1)\n#until we got to 1,2, at which point, the ending \+ of rest_neighbor would be\n#1,2, indicating that we have gotten to the end of the line.\n\n\nline_vertices:=proc(m_func,startpt,first_neighb or,rest_neighbor)\nlocal lineT,vm,outlist,ending,last_pt;\n\nlineT:=ga sket(m_func);\nvm:=vertices(m_func);\nmake_neighbors(lineT,m_func,vm); \n\noutlist:=[startpt,lineT[startpt,first_neighbor]];\n\nending:=last_ digits(lineT[startpt,first_neighbor]);\nlast_pt:=lineT[startpt,first_n eighbor];\n\nwhile last_digits(lineT[last_pt,rest_neighbor])=ending do \n outlist:=[op(outlist),lineT[last_pt,rest_neighbor]];\n last_pt: =lineT[last_pt,rest_neighbor];\nod;\n\noutlist:=[op(outlist),lineT[las t_pt,rest_neighbor]];\n\noutlist;\nend;\n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.line_verticesGf*6&%'m_funcG%(startptG%/first_neighbo rG%.rest_neighborG6'%&lineTG%#vmG%(outlistG%'endingG%(last_ptG6\"F1C+> 8$-%'gasketG6#9$>8%-%)verticesGF7-%/make_neighborsG6%F4F8F:>8&7$9%&F46 $FC9&>8'-%,last_digitsG6#FD>8(FD?(F1\"\"\"FOF1/-FJ6#&F46$FM9'FHC$>FA7$ -%#opG6#FAFS>FMFS>FAFXFAF1F1F1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1099 "#This function is just like make_ animation, except that it makes 2d line plots\n#of the function along \+ a line defined by startpt, first_neighbor, rest_neighbor.\n#You also n eed to input linelength to tell the program how long the line should\n #be. For instance, the line from [0,2] to [1,2] would require linelen gth to be\n#1/2. For more information about how startpt, first_neighb or, and rest_neighbor\n#define a line, look at line_vertices above.\n \n\nmake_animation_line:=proc(T,m_basis,m_func,vals,funcno,harm_funcno ,startpt,first_neighbor,rest_neighbor,linelength)\n\nlocal linepts,S,v m,SG_plots,final_points,inds,i;\n\nSG_plots:=[];\n\nS:=gasket(m_basis) ;\nvm:=vertices(m_basis);\nmake_neighbors(S,m_basis,vm);\nmake_points( T,m_func);\n\nlinepts:=line_vertices(m_func,startpt,first_neighbor,res t_neighbor);\n\nfor i from 1 to nops(vals) do\n\n print(`computing p lot number `,i);\n make_pretty_picture(T,S,m_basis,m_func,vm,funcno, harm_funcno,vals[i]);\n final_points:=SG_plot_line(T,m_func,funcno,l inepts,linelength);\n SG_plots:=[op(SG_plots),plot(final_points,colo r=BLUE,symbol=POINT)];\nod;\n\nSG_plots;\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%4make_animation_lineGf*6,%\"TG%(m_basisG%'m_funcG%%va lsG%'funcnoG%,harm_funcnoG%(startptG%/first_neighborG%.rest_neighborG% +linelengthG6)%(lineptsG%\"SG%#vmG%)SG_plotsG%-final_pointsG%%indsG%\" iG6\"F9C*>8'7\">8%-%'gasketG6#9%>8&-%)verticesGFB-%/make_neighborsG6%F ?FCFE-%,make_pointsG6$9$9&>8$-%.line_verticesG6&FO9*9+9,?(8*\"\"\"FZ-% %nopsG6#9'%%trueGC&-%&printG6$%7computing~plot~number~GFY-%4make_prett y_pictureG6*FNF?FCFOFE9(9)&Fhn6#FY>8(-%-SG_plot_lineG6'FNFOFboFQ9->F<7 $-%#opG6#F<-%%plotG6%Fgo/%&colorG%%BLUEG/%'symbolG%&POINTGF " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 355 "#This function is just like SG_plot, except that i t creates a 2d\n#line plot along a line which is defined by linepts, a nd is of length\n#linelength.\n\nSG_plot_line:=proc(T,m,funcno,linepts ,linelength)\n local i,plot_pts,numpts;\n\n numpts:=nops(linepts);\n \n plot_pts:=[seq([(i-1)*linelength/(numpts-1),T[linepts[i],funcno]], i=1..numpts)];\n\n plot_pts;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-SG_plot_lineGf*6'%\"TG%\"mG%'funcnoG%(lineptsG%+linelengthG6% %\"iG%)plot_ptsG%'numptsG6\"F0C%>8&-%%nopsG6#9'>8%7#-%$seqG6$7$*(,&8$ \"\"\"FB!\"\"FB9(FB,&F3FBFBFCFC&9$6$&F76#FA9&/FA;FBF3F9F0F0F0" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 863 "#this is just a function to close up extraneous fi les...\n\nclose_files:=proc(maxk)\n\nlocal S,vm,EM,IP,k,outfile,i,numo ps,ilist,jlist,elist,iplist,filename_EM,filename_IP, filename_EMi,file name_EMj,filename_EMx,filename_IPi,filename_IPj,filename_IPx,fEMx,fEMi ,fEMj,fIPi,fIPj,fIPx,directory,inds;\n\ndirectory:=\"/tahoe/home/raj/ \";\n\nfor k from 1 to maxk do\nfilename_EMi:=cat(directory,\"EM\",con vert(k,string),\"_i.dat\");\nfilename_EMj:=cat(directory,\"EM\",conver t(k,string),\"_j.dat\");\nfilename_EMx:=cat(directory,\"EM\",convert(k ,string),\"_x.dat\");\n\nfilename_IPi:=cat(directory,\"IP\",convert(k, string),\"_i.dat\");\nfilename_IPj:=cat(directory,\"IP\",convert(k,str ing),\"_j.dat\");\nfilename_IPx:=cat(directory,\"IP\",convert(k,string ),\"_x.dat\");\n\nclose(filename_EMi);\nclose(filename_EMj);\nclose(fi lename_EMx);\nclose(filename_IPi);\nclose(filename_IPj);\nclose(filena me_IPx);\n\n\nod;\n\n\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%,cl ose_filesGf*6#%%maxkG6>%\"SG%#vmG%#EMG%#IPG%\"kG%(outfileG%\"iG%'numop sG%&ilistG%&jlistG%&elistG%'iplistG%,filename_EMG%,filename_IPG%-filen ame_EMiG%-filename_EMjG%-filename_EMxG%-filename_IPiG%-filename_IPjG%- filename_IPxG%%fEMxG%%fEMiG%%fEMjG%%fIPiG%%fIPjG%%fIPxG%*directoryG%%i ndsG6\"FEC$>8>Q1/tahoe/home/raj/FE?(8(\"\"\"FL9$%%trueGC.>82-%$catG6&F HQ#EMFE-%(convertG6$FK%'stringGQ'_i.datFE>83-FS6&FHFUFVQ'_j.datFE>84-F S6&FHFUFVQ'_x.datFE>85-FS6&FHQ#IPFEFVFZ>86-FS6&FHFcoFVFin>87-FS6&FHFco FVF^o-%&closeG6#FQ-F]p6#Ffn-F]p6#F[o-F]p6#F`o-F]p6#Feo-F]p6#FioFEFEFE " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3269 "#This function genera tes EM and IP files which can be read as sparse matrices\n#into matlab . It generates files like EM1_i, EM2_j,IP4_x. The EM and the IP are \n#stand for energy and inner product (who woulda thought...) and the \+ i, j, and x\n#stand for the ith row, jth column, and value at that coo rdinate. These \n#3 files can be read into matlab as 3 column vectors , which can then be used \n#with the sparse command to generate a spar se matrix. The number after the\n#EM and the IP is the level of m use d with a biharmonic basis. This function\n#is really not very robust, and should be used with care. Note also that the\n#output files may \+ have to have their unix permissions changed with chmod. Ask\n#a compu ter geek to help you out with that if you don't know what I'm talking \n#about... \n\n \ngenerate_EM_IP_files:=proc(maxk)\n\nlocal S,vm,EM,I P,k,outfile,i,numops,ilist,jlist,elist,iplist,filename_EM,filename_IP, filename_EMi,filename_EMj,filename_EMx,filename_IPi,filename_IPj,file name_IPx,fEMx,fEMi,fEMj,fIPi,fIPj,fIPx,directory,inds;\n\nfor k from 1 to maxk do\n\nprint(k);\nS:=gasket(k):\neval(S):\nvm:=vertices(k):\nm ake_neighbors(S,k,vm):\neval(S):\n\nprint(`calculating matrices...`); \nEM:=eval(gen_energy_matrix(S,k,vm));\nprint(`calculated EM...`);\nIP :=eval(gen_iprod_matrix(S,k,vm));\nprint(`calculated IP...`);\n\nnumop s:=nops([indices(EM)]);\n\n#ilist:=[seq([indices(EM)][i][1],i=1..numop s)];\n#jlist:=[seq([indices(EM)][i][2],i=1..numops)];\n\n#elist:=evalf ([seq(EM[[indices(EM)][i][1],[indices(EM)][i][2]],i=1..numops)],20);\n #iplist:=evalf([seq(IP[[indices(IP)][i][1],[indices(IP)][i][2]],i=1..n umops)],20);\n\nfilename_EM:=cat(\"EM\" , convert(k,string));\nfilenam e_IP:=cat(\"IP\" , convert(k,string));\n\ndirectory:=\"/tahoe/home/raj /\";\n\nprint(`saving .m files...`);\nsave EM,cat(directory,filename_E M,\".m\");\nsave IP,cat(directory,filename_IP,\".m\");\n\nfilename_EMi :=cat(directory,\"EM\",convert(k,string),\"_i.dat\");\nfilename_EMj:=c at(directory,\"EM\",convert(k,string),\"_j.dat\");\nfilename_EMx:=cat( directory,\"EM\",convert(k,string),\"_x.dat\");\n\nfilename_IPi:=cat(d irectory,\"IP\",convert(k,string),\"_i.dat\");\nfilename_IPj:=cat(dire ctory,\"IP\",convert(k,string),\"_j.dat\");\nfilename_IPx:=cat(directo ry,\"IP\",convert(k,string),\"_x.dat\");\n\n\n#close(filename_EMi);\n# close(filename_EMj);\n#close(filename_EMx);\n#close(filename_IPi);\n#c lose(filename_IPj);\n#close(filename_IPx);\n\nfEMi:=open(filename_EMi, WRITE);\nfEMj:=open(filename_EMj,WRITE);\nfEMx:=open(filename_EMx,WRIT E);\nfIPi:=open(filename_IPi,WRITE);\nfIPj:=open(filename_IPj,WRITE); \nfIPx:=open(filename_IPx,WRITE);\n\nprint(`copying to text files...`) ;\n\ninds:=[indices(EM)];\n\nprint(`(got indices...)`);\n\nfor i from \+ 1 to numops do\n# fprintf(fEMi,`%20.20g\\n`,ilist[i]);\n# fprintf( fEMj,`%20.20g\\n`,jlist[i]);\n# fprintf(fEMx,`%20.20g\\n`,elist[i]); \n\n# fprintf(fIPi,`%20.20g\\n`,ilist[i]);\n# fprintf(fIPj,`%20.20 g\\n`,jlist[i]);\n# fprintf(fIPx,`%20.20g\\n`,iplist[i]);\n\n fpri ntf(fEMi,`%20.20g\\n`,inds[i][1]);\n fprintf(fEMj,`%20.20g\\n`,inds[ i][2]);\n fprintf(fEMx,`%20.20g\\n`,EM[inds[i][1],inds[i][2]]);\n\n \+ fprintf(fIPi,`%20.20g\\n`,inds[i][1]);\n fprintf(fIPj,`%20.20g\\n` ,inds[i][2]);\n fprintf(fIPx,`%20.20g\\n`,IP[inds[i][1],inds[i][2]]) ;\nod; \n\nclose(fEMi);\nclose(fEMj);\nclose(fEMx);\nclose(fIPi);\nclo se(fIPj);\nclose(fIPx);\n\nod;\n\nend;\n\n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5generate_EM_IP_filesGf*6#%%maxkG6>%\"SG%#vmG%#EMG%#I PG%\"kG%(outfileG%\"iG%'numopsG%&ilistG%&jlistG%&elistG%'iplistG%,file name_EMG%,filename_IPG%-filename_EMiG%-filename_EMjG%-filename_EMxG%-f ilename_IPiG%-filename_IPjG%-filename_IPxG%%fEMxG%%fEMiG%%fEMjG%%fIPiG %%fIPjG%%fIPxG%*directoryG%%indsG6\"FE?(8(\"\"\"FH9$%%trueGCJ-%&printG 6#FG>8$-%'gasketGFN-%%evalG6#FP>8%-%)verticesGFN-%/make_neighborsG6%FP FGFWFS-FM6#%8calculating~matrices...G>8&-FT6#-%2gen_energy_matrixGFfn- FM6#%1calculated~EM...G>8'-FT6#-%1gen_iprod_matrixGFfn-FM6#%1calculate d~IP...G>8+-%%nopsG6#7#-%(indicesG6#F[o>80-%$catG6$Q#EMFE-%(convertG6$ FG%'stringG>81-Fhp6$Q#IPFEF[q>8>Q1/tahoe/home/raj/FE-FM6#%3saving~.m~f iles...GB6$F[o-Fhp6%FeqFfpQ#.mFEB6$Fdo-Fhp6%FeqF`qF^r>82-Fhp6&FeqFjpF[ qQ'_i.datFE>83-Fhp6&FeqFjpF[qQ'_j.datFE>84-Fhp6&FeqFjpF[qQ'_x.datFE>85 -Fhp6&FeqFcqF[qFgr>86-Fhp6&FeqFcqF[qF\\s>87-Fhp6&FeqFcqF[qFas>89-%%ope nG6$Fdr%&WRITEG>8:-Fat6$FirFct>88-Fat6$F^sFct>8;-Fat6$FcsFct>8<-Fat6$F gsFct>8=-Fat6$F[tFct-FM6#%9copying~to~text~files...G>8?Fap-FM6#%1(got~ indices...)G?(8*FHFHF]pFJC(-%(fprintfG6%F_t%)%20.20g|+G&&F\\v6#Fav6#FH -Fdv6%FetFfv&Fhv6#\"\"#-Fdv6%FitFfv&F[o6$FgvF]w-Fdv6%F]uFfvFgv-Fdv6%Fa uFfvF]w-Fdv6%FeuFfv&FdoFcw-%&closeG6#F_t-F\\x6#Fet-F\\x6#Fit-F\\x6#F]u -F\\x6#Fau-F\\x6#FeuFEFEFE" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "\n#generate_EM_IP_fi les(6);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 264 "#This \+ function multiplies func1 and func2 pointwise and stores \n#the result in outfunc...\n\nmult:=proc(T,m,func1,func2,outfunc) \n\nlocal inds,i ;\n\ninds:=vertices(m);\n\nfor i from 1 to nops(inds) do\n T[inds[i ],outfunc]:=T[inds[i],func1]*T[inds[i],func2];\nod;\n\nend; " }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%%multGf*6'%\"TG%\"mG%&func1G%&func2G %(outfuncG6$%%indsG%\"iG6\"F/C$>8$-%)verticesG6#9%?(8%\"\"\"F9-%%nopsG 6#F2%%trueG>&9$6$&F26#F89(*&&F@6$FB9&F9&F@6$FB9'F9F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 244 "#This function takes the sine of f unc1 (pointwise) and stores\n#it in outfunc.\nsinof:=proc(T,m,func1,ou tfunc)\n local inds,i;\n\n inds:=vertices(m); \n\n for i from 1 \+ to nops(inds) do\n T[inds[i],outfunc]:=sin(T[inds[i],func1]);\n \+ od;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%&sinofGf*6&%\"TG%\"mG% &func1G%(outfuncG6$%%indsG%\"iG6\"F.C$>8$-%)verticesG6#9%?(8%\"\"\"F8- %%nopsG6#F1%%trueG>&9$6$&F16#F79'-%$sinG6#&F?6$FA9&F.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 325 "#This procedure takes func1,func2, mult1, and mult2, and stores\n# mult1*func1+mult2*func2 in outfunc (ag ain, pointwise).\n\nadd_func_lin:=proc(T,m,func1,func2,mult1,mult2,out func) \nlocal inds,i;\n\ninds:=vertices(m);\n\nfor i from 1 to nops(in ds) do\n T[inds[i],outfunc]:=mult1*T[inds[i],func1]+mult2*T[inds[i] ,func2];\nod;\n\nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-add_func _linGf*6)%\"TG%\"mG%&func1G%&func2G%&mult1G%&mult2G%(outfuncG6$%%indsG %\"iG6\"F1C$>8$-%)verticesG6#9%?(8%\"\"\"F;-%%nopsG6#F4%%trueG>&9$6$&F 46#F:9*,&*&9(F;&FB6$FD9&F;F;*&9)F;&FB6$FD9'F;F;F1F1F1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2047 "#So this function stores the lapl acian of harm_funcno+(odd number up to 11)\n#and stores it in harm_fun cno+12+(odd number up to 11). It stores this by\n#using the laplacian s of the \"old\" basis functions, which are stored in \n#harm_funcno+( odd number up to 11) (values are overridden later).\n#See calc_biharms in splines.mws for full details on this process. \n\nf_g_lap:=proc(T, m,harm_funcno)\n\nlocal inds,i;\n\ninds:=vertices(m);\n\nfor i from 1 \+ to nops(inds) do\n\n#Uses the 30-15-15 rule...\n\nT[inds[i],harm_funcn o + 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];\n\nT[inds [i],harm_funcno + 12+5]:=T[inds[i],harm_funcno+3]-30*T[inds[i],harm_fu ncno+9] +15 * T[inds[i],harm_funcno + 7]+15 * T[inds[i],harm_funcno + \+ 11];\n\nT[inds[i],harm_funcno + 12+9]:=T[inds[i],harm_funcno+5]-30*T[i nds[i],harm_funcno+11] +15 * T[inds[i],harm_funcno + 9]+15 * T[inds[i] ,harm_funcno + 7];\n\nT[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];\n\nT[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];\n\nT[ 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];\n\nod;\nend;\n\n# This function creates the f and the g type biharmonic splines\n#from t he old basis splines by using the 30-15-15 rule and the 4-11-11 rule. \n#It puts the f-type vector in out_funcno, and then puts\n#the g-type vector in out_funcno+2.\n\nf_g_bih:=proc(T,m,num,harm_funcno,out_func no)\n\nlocal inds,i;\n\ninds:= vertices(m);\n\nfor i from 1 to nops(in ds) do\n\n T[inds[i],out_funcno]:=T[inds[i],harm_funcno+num]-30*T[i nds[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];\n 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];\nod;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6# >%(f_g_lapGf*6%%\"TG%\"mG%,harm_funcnoG6$%%indsG%\"iG6\"F-C$>8$-%)vert icesG6#9%?(8%\"\"\"F7-%%nopsG6#F0%%trueGC(>&9$6$&F06#F6,&9&F7\"#8F7,*& F?6$FA,&FDF7F7F7F7*&\"#IF7&F?6$FA,&FDF7\"\"(F7F7!\"\"*&\"#:F7&F?6$FA,& FDF7\"\"*F7F7F7*&FRF7&F?6$FA,&FDF7\"#6F7F7F7>&F?6$FA,&FDF7\"#&F?6$FA,&FDF7\" #@F7,*&F?6$FA,&FDF7\"\"&F7F7*&FKF7FXF7FP*&FRF7FSF7F7*&FRF7FLF7F7>&F?6$ FA,&FDF7FRF7,(FLFen*&\"\"%F7FSF7FP*&FfpF7FXF7FP>&F?6$FA,&FDF7\"#>F7,(F SFen*&FfpF7FLF7FP*&FfpF7FXF7FP>&F?6$FA,&FDF7\"#BF7,(FXFen*&FfpF7FSF7FP *&FfpF7FLF7FPF-F-F-" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%(f_g_bihGf*6' %\"TG%\"mG%$numG%,harm_funcnoG%+out_funcnoG6$%%indsG%\"iG6\"F/C$>8$-%) verticesG6#9%?(8%\"\"\"F9-%%nopsG6#F2%%trueGC$>&9$6$&F26#F89(,*&FA6$FC ,&9'F99&F9F9*&\"#IF9&FA6$FC,(FJF9\"\"'F9FKF9F9!\"\"*&\"#:F9&FA6$FC,(FJ F9-%$modG6$,&\"\")F9FKF9FQF9FQF9F9F9*&FTF9&FA6$FC,(FJF9-FY6$,&\"#5F9FK F9FQF9FQF9F9F9>&FA6$FC,&FEF9\"\"#F9,(&FA6$FC,(FJF9-FY6$,&FQF9FKF9FQF9F QF9\"#6*&\"\"%F9FUF9FR*&F]pF9FhnF9FRF/F/F/" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 517 "#This function implements the trapzoid rule for in tegration (out to level m).\n#It should be noted that simpson's rule i s significantly better\n#than the trapezoid rule, and simp_rule should be used whenever\n#possible.\n\ntrap_rule:=proc(T,m,funcno)\n\nlocal \+ i, func, V0, trap,Vother,Vother_e;\n \n\nV0:=vertices_m(0);\nVother_ e:=vertices(m); Vother:=Vother_e[4..nops(Vother_e)];\n\n trap := eva l(3^(-m-1) * ( sum(2 * T[Vother[i],funcno], i = 1..nops(Vother) ) + su m(T[V0[i],funcno], i = 1..nops(V0)) )); \n\n trap;\n \nend;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 837 "#This function implements simpson' s rule. This function gives\n#much better results than trap_rule. It should be noted that this algorithm is\n#used for all the numerical i ntegrations performed when implementing the FEM.\n#However, most of th e time, instead of using simp_rule, the algorithm is\n#implemented wit hin the procedure itself as a matter of convenience and efficiency.\n \nsimp_rule:=proc(T,m,funcno)\n\nlocal i, func, V0,Vm, simp,Vother_e,V other;\n \n#get V0, Vm, Veverythinginbetween...\nV0:=vertices_m(0); \nif m =1 then Vother:=[] else Vother_e:=vertices(m-1); Vother:=Vother _e[4..nops(Vother_e)]; fi;\nVm:=vertices_m(m);\n\n#and then perform th e sums...\nsimp := 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)) ); \n\nsimp;\n \nend;" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%*trap_ruleGf*6%%\"TG%\"mG%'funcnoG6(%\"iG%%fun cG%#V0G%%trapG%'VotherG%)Vother_eG6\"F1C'>8&-%+vertices_mG6#\"\"!>8)-% )verticesG6#9%>8(&F:6#;\"\"%-%%nopsG6#F:>8'-%%evalG6#*&)\"\"$,&F>!\"\" \"\"\"FQFR,&-%$sumG6$,$&9$6$&F@6#8$9&\"\"#/Fgn;FR-FF6#F@FR-FU6$&FY6$&F 4FfnFhn/Fgn;FR-FF6#F4FRFRFIF1F1F1" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#> %*simp_ruleGf*6%%\"TG%\"mG%'funcnoG6)%\"iG%%funcG%#V0G%#VmG%%simpG%)Vo ther_eG%'VotherG6\"F2C'>8&-%+vertices_mG6#\"\"!@%/9%\"\"\">8*7\"C$>8)- %)verticesG6#,&FF?&FC6#;\"\"%-%%nopsG6#FC>8'-F76#F<>8(*&-%% evalG6#,$*&F=F=)\"\"$F " 0 "" {MPLTEXT 1 0 890 "#This function takes the laplacian of funcno and sto res it in funcno_out.\n#Make_neighbors must be run on T before this fu nction is run.\n#Also, it should be noted that the laplacian doesn't e xist on the \n#boundaries of the gasket when using the pointwise formu la used here;\n#instead, I just set the values on the edges to 0 for e ase of programming.\n#It can be changed quite easily...\n\n\nfind_lap: =proc(T,m,funcno,funcno_out)\n\n local i, j, Vm;\n Vm:=vertices( m);\n\n#assign the values on the boundaries...\n for i from 1 to 3 \+ do\n T[Vm[i],funcno_out]:=0; od;\n\n # T[Vm[i], funcno _out] := (3/2) * (5^(m)) * (add( T[T[Vm[i],j],funcno] - (2 * T[Vm[i],f uncno]), j = 1..2)); od;\n\n\n#and now do the rest...\n for i fro m 4 to nops(Vm) do\n\n T[Vm[i], funcno_out] := (3/2) * (5^(m)) * add( T[T[Vm[i],j],funcno] - T[Vm[i],funcno], j = 1..4) ; od;\n \n end;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)find_lap Gf*6&%\"TG%\"mG%'funcnoG%+funcno_outG6%%\"iG%\"jG%#VmG6\"F/C%>8&-%)ver ticesG6#9%?(8$\"\"\"F9\"\"$%%trueG>&9$6$&F26#F89'\"\"!?(F8\"\"%F9-%%no psG6#F2F;>F=,$*&)\"\"&F6F9-%$addG6$,&&F>6$&F>6$F@8%9&F9&F>6$F@FW!\"\"/ FV;F9FEF9#F:\"\"#F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 163 "#This function returns the last 2 digits of address.\n#It is used to \+ determine orientation.\nlast_digits:=proc(address)\naddress[nops(addre ss)-1..nops(address)]; end;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,last _digitsGf*6#%(addressG6\"F(F(&9$6#;,&-%%nopsG6#F*\"\"\"F1!\"\"F.F(F(F( " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1129 "#This function return s the addresses of the boundaries of the 2\n#triangles of which basis[ 1] is a common vertex. It requires that the \n#spline basis table S h as been make_neighbored. If which_triangle=1 then\n#it returns the fi rst triangle (in clockwise orientation), otherwise it returns\n#the se cond triangle (in clockwise orientation). It returns the addresses in \n#clockwise fashion, with the first vertex being the one furthest up (ie, \n#largest y value in the euclidean sense). These addresses are often used in\n#coord_transfer.\n\nget_addresses:=proc(S,basis,m_basi s,which_triangle)\nlocal digits;\ndigits:=last_digits(basis[1]);\n\nif digits=[0,1] then\n if which_triangle=1 then RETURN([basis[1],S[bas is[1],1],S[basis[1],2]]);\n else RETURN([S[basis[1],4],basis[1],S[ba sis[1],3]]); fi;\nfi;\nif digits = [0,2] then\n if which_triangle=1 \+ then RETURN([S[basis[1],1],S[basis[1],2],basis[1]]);\n else RETURN([ basis[1],S[basis[1],3],S[basis[1],4]]); fi;\nfi;\nif digits = [1,2] th en\n if which_triangle=1 then RETURN([S[basis[1],2],basis[1],S[basis [1],1]]);\n else RETURN([S[basis[1],3],S[basis[1],4],basis[1]]); fi; \nfi;\nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.get_addressesGf* 6&%\"SG%&basisG%(m_basisG%/which_triangleG6#%'digitsG6\"F-C&>8$-%,last _digitsG6#&9%6#\"\"\"@$/F07$\"\"!F7@%/9'F7-%'RETURNG6#7%F4&9$6$F4F7&FD 6$F4\"\"#-F@6#7%&FD6$F4\"\"%F4&FD6$F4\"\"$@$/F07$F;FH@%F=-F@6#7%FCFFF4 -F@6#7%F4FOFL@$/F07$F7FH@%F=-F@6#7%FFF4FC-F@6#7%FOFLF4F-F-F-" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 496 "#This generates a function \+ with values 0 on (V1-V0), and sets the boundaries to value,\n#and then fills in the rest of the function harmonically, thus generating\n#a p iecewise harmonic function with full D3 symmetry.\n\nsym_func:=proc(T, m,value,funcno)\n\nT[[0],funcno]:=value;\nT[[1],funcno]:=value;\nT[[2] ,funcno]:=value;\nT[[0,1],funcno]:=0;\nT[[0,2],funcno]:=0;\nT[[1,2],fu ncno]:=0;\n\n\nharm(T,m-1,[2],[0,2], [1,2], funcno);\nharm(T,m-1,[0],[ 0,1], [0,2], funcno);\nharm(T,m-1,[1],[1,2], [0,1], funcno);\nend;" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%)sym_funcGf*6&%\"TG%\"mG%&valueG%'fu ncnoG6\"F+F+C+>&9$6$7#\"\"!9'9&>&F/6$7#\"\"\"F3F4>&F/6$7#\"\"#F3F4>&F/ 6$7$F2F9F3F2>&F/6$7$F2F>F3F2>&F/6$7$F9F>F3F2-%%harmG6(F/,&9%F9F9!\"\"F =FFFJF3-FL6(F/FNF1FBFFF3-FL6(F/FNF8FJFBF3F+F+F+" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 365 "#This procedure takes funcno, and (given IPin v = inverse(IP) where\n#IP is the inner product matrix created by gen_ iprod_matrix) expresses\n#the function in the spline basis.\n\n\nexpre ss_in_basis:=proc(T,IPinv,S,m_basis,m_func,funcno,harm_funcno)\n\nloca l vec,vm;\nvm:=vertices(m_basis);\nvec:=gen_f_vector(T,S,m_basis,m_fun c,vm,funcno,harm_funcno);\nevalm(IPinv&*vec);\nend;" }}{PARA 12 "" 1 " " {XPPMATH 20 "6#>%1express_in_basisGf*6)%\"TG%&IPinvG%\"SG%(m_basisG% 'm_funcG%'funcnoG%,harm_funcnoG6$%$vecG%#vmG6\"F1C%>8%-%)verticesG6#9' >8$-%-gen_f_vectorG6)9$9&F89(F49)9*-%&evalmG6#-%#&*G6$9%F:F1F1F1" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 295 "#This function takes a vect or of length n_elts\n#and writes it to a file. The subsequent\n#colum n vector can be read into matlab quite easily.\n\nwrite_vec:=proc(file name,vec,n_elts)\nlocal fd,i;\nfd:=fopen(filename,WRITE);\nfor i from \+ 1 to n_elts do\n fprintf(fd,\"%g\\n\",vec[i]);\nod;\nfclose(fd);\nen d;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%*write_vecGf*6%%)filenameG%$ve cG%'n_eltsG6$%#fdG%\"iG6\"F-C%>8$-%&fopenG6$9$%&WRITEG?(8%\"\"\"F89&%% trueG-%(fprintfG6%F0Q$%g|+F-&9%6#F7-%'fcloseG6#F0F-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 213 "#This function simply takes funcno 1 and copies it to funcno2\n\ncopy_function:=proc(T1,T2,m,funcno1,func no2)\nlocal i,vm;\n\nvm:=vertices(m);\nfor i from 1 to nops(vm) do\n \+ T2[vm[i],funcno2]:=T1[vm[i],funcno1];\nod;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.copy_functionGf*6'%#T1G%#T2G%\"mG%(funcno1G%(funcno2 G6$%\"iG%#vmG6\"F/C$>8%-%)verticesG6#9&?(8$\"\"\"F9-%%nopsG6#F2%%trueG >&9%6$&F26#F89(&9$6$FB9'F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 518 "#This function will read (out of filename), and put the \n#re sults in a list of lists. Each sublist will contain n_elements\n#and \+ corresponds to a row of the file.\n#ie, if your file looks like:\n#1 2 \n#3 4\n#5 6\n#then running read_basis_vals returns:\n#[[1,2],[3,4],[5 ,6]]\n#This is used for reading in time-series solutions to the wave\n #and heat equations.\n\nread_basis_vals:=proc(filename,n_elements)\n\n local i,fd,vals,line,outvals;\n\n\nfd:=fopen(filename,READ);\n\nvals:= readdata(fd,n_elements,float);\n\nfclose(fd);\n\nvals;\n\nend;" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%0read_basis_valsGf*6$%)filenameG%+n_ elementsG6'%\"iG%#fdG%%valsG%%lineG%(outvalsG6\"F/C&>8%-%&fopenG6$9$%% READG>8&-%)readdataG6%F29%%&floatG-%'fcloseG6#F2F9F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "with(plots):with(linalg):" }}{PARA 7 "" 1 "" {TEXT -1 50 "Warning, the name changecoords has been redefined\n" }}{PARA 7 "" 1 "" {TEXT -1 80 "Warning, the protected names norm and trace have \+ been redefined and unprotected\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 561 "#This procedure returns the matrix E(fi,fj) where fi , fj are spline basis functions.\n#it stores values in sparse matrix f orm, and also takes advantage of the fact that the\n#matrix is symmetr ic...\n\ngen_energy_matrix:=proc(S,m_basis,vm)\n\nlocal i,j,en_matrix, en;\n\n\nen_matrix:=array(sparse,1..3^(m_basis+1),1..3^(m_basis+1));\n \nfor i from 1 to 3^(m_basis+1) do\n for j from 1 to i do\n en: =energy(S,m_basis,get_basis(vm,i),get_basis(vm,j));\n if en<>0 th en\n en_matrix[i,j]:=en;\n en_matrix[j,i]:=en;\n fi ;\n od;\nod;\n\nen_matrix;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%2gen_energy_matrixGf*6%%\"SG%(m_basisG%#vmG6&%\"iG%\" jG%*en_matrixG%#enG6\"F/C%>8&-%&arrayG6%%'sparseG;\"\"\")\"\"$,&9%F8F8 F8F7?(8$F8F8F9%%trueG?(8%F8F8F>F?C$>8'-%'energyG6&9$F<-%*get_basisG6$9 &F>-FJ6$FLFA@$0FD\"\"!C$>&F26$F>FAFD>&F26$FAF>FDF2F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 466 "#This function returns the biharmo nic basis element corresponding to number. The bases are numbered by \+ the\n#ordering provided by vm: the first 3 elements of vm are the boun daries, and there is no f type basis vectors\n#at those points. The f irst bunch of basis vectors are basically f-type basis elements, and t he rest are g-type.\n\nget_basis:=proc(vm,number)\noption remember;\n \nif number <= nops(vm)-3 then [vm[number+3],f]\nelse [vm[number-nops( vm)+3],g] fi; end; \n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%*get_basi sGf*6$%#vmG%'numberG6\"6#%)rememberGF)@%19%,&-%%nopsG6#9$\"\"\"\"\"$! \"\"7$&F36#,&F.F4F5F4%\"fG7$&F36#,(F.F4F0F6F5F4%\"gGF)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 451 "#This operates just the same way a s gen_energy_matrix, except using inner products...\n\ngen_iprod_matri x:=proc(S,m_basis,vm)\n\nlocal i,j,i_matrix,ip;\n\n\ni_matrix:=array(s parse,1..3^(m_basis+1),1..3^(m_basis+1));\n\nfor i from 1 to 3^(m_basi s+1) do\n for j from 1 to i do\n ip:=inner_prod(S,m_basis,get_b asis(vm,i),get_basis(vm,j));\n if ip <> 0 then\n i_matrix[ i,j]:=ip;\n i_matrix[j,i]:=ip;\n fi;\n od;\nod;\n\ni_mat rix;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%1gen_iprod_mat rixGf*6%%\"SG%(m_basisG%#vmG6&%\"iG%\"jG%)i_matrixG%#ipG6\"F/C%>8&-%&a rrayG6%%'sparseG;\"\"\")\"\"$,&9%F8F8F8F7?(8$F8F8F9%%trueG?(8%F8F8F>F? C$>8'-%+inner_prodG6&9$F<-%*get_basisG6$9&F>-FJ6$FLFA@$0FD\"\"!C$>&F26 $F>FAFD>&F26$FAF>FDF2F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 499 "#This procedure gen erates the f vector, which is basically f integrates versus all the di fferent basis functions.\n#It calls int_vs_basis, which contains the m eat of the functionality.\n\ngen_f_vector:=proc(T,S,m_basis,m_func,vm, funcno,harm_funcno)\n\nlocal i,f_vec;\n\nf_vec:=array(1..3^(m_basis+1) );\n\nfor i from 1 to 3^(m_basis+1) do\n\n printf(`%d `,i);\n if ( (i mod 15)=1) then\n printf(`\\n`): fi:\n\n f_vec[i]:=int_vs_ba sis(T,S,m_basis,m_func,get_basis(vm,i),funcno,harm_funcno);\nod;\n\nf_ vec;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-gen_f_vectorGf*6)%\"T G%\"SG%(m_basisG%'m_funcG%#vmG%'funcnoG%,harm_funcnoG6$%\"iG%&f_vecG6 \"F1C%>8%-%&arrayG6#;\"\"\")\"\"$,&9&F9F9F9?(8$F9F9F:%%trueGC%-%'print fG6$%$%d~GF?@$/-%$modG6$F?\"#:F9-FC6#%\"|+G>&F46#F?-%-int_vs_basisG6)9 $9%F=9'-%*get_basisG6$9(F?9)9*F4F1F1F1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 "\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7581 "#Th is function creates a actual function (which can then be plotted and t he such) from the values\n#of the spline projection stored in C. It c reates the function simply by going through each basis\n#element and a dding in the appropriate amount of it. Sounds, simple, huh?\n\nmake_p retty_picture:=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno,C)\n\n\n \nlocal i,j,vm2,k,vmh,triangle1,triangle2,basis,add0,add1,add2,which_h arm1,which_harm2,coord1,coord2;\n\n\n#first, let's set the initial val ues of our function to 0.\nvm2:=vertices(m_func);\nfor i from 1 to nop s(vm2) do\n T[vm2[i],funcno]:=0;\nod;\n\nvmh:=vertices(m_func-m_basi s); #These are the vertices which we will need from our spline basis \+ elements\n\nfor i from 1 to 3^(m_basis+1) do\n\nbasis:=get_basis(vm,i) ; #This is the basis element we will be working with.\nprintf(`%d `,i );\nif ((i mod 15)=1) then\n printf(`\\n`): fi:\n\n\n#So if the basis element is in V0, we know it is a g-type basis vector.\n#the scheme h ere is to find the appropriate triangle (which will be some\n#sequence of 0's or 1's or 2's) and also get the coordinates of the boundaries \n#of this triangle (using the neighbors defined in S; like S[basis[1] ,2]).\n#also, depending on the orientation of the basis, we need to pi ck the appropriate\n#basis element (so for basis = [[0],g], we need a \+ basis element which has values\n#of zero everywhere, and has normal de rivative 1 at [0], and normal derivative 0\n#at the other boundaries. \+ This is stored in harm_funcno+3). This explains the assignments\n#of triangle1, add0, add1, add2, and which_harm1. The function coord_tra nsfer is then used\n#to take a coordinate from the basis function and \+ transform it to the appropriate part\n#of the function we are trying t o generate (basically, by adding some more digits to\n#the front). Th en add in the value of the basis function * C[i] and voila!\n\nif nops (basis[1]) = 1 then #this is if the basis is an element of V0...\n \+ if basis[1] = [0] then\n triangle1:=[seq(0,k=1..m_basis)];\n \+ add0:=basis[1]; add1:=S[basis[1],1]; add2:=S[basis[1],2];which_ha rm1:=3; \n fi;\n if basis[1] = [1] then\n triangle1:=[seq (1,k=1..m_basis)];\n add0:=S[basis[1],2]; add1:=basis[1]; add2: =S[basis[1],1];which_harm1:=7;\n fi;\n if basis[1] = [2] then\n \+ triangle1:=[seq(2,k=1..m_basis)];\n add0:=S[basis[1],1]; add1:=S[basis[1],2]; add2:=basis[1];which_harm1:=11;\n fi;\n fo r j from 1 to nops(vmh) do\n coord1:=coord_transfer(vmh[j],triangle 1,add0,add1,add2);\n T[coord1,funcno]:=T[coord1,funcno]+C[i]*T[vmh[ j],harm_funcno + which_harm1];\n od; \n \nelse\n\n#Here, the strat egy is the same, except that we have 2 triangles to deal with: one on \+ each side of the\n#basis vertex. Notice that the values we feed into \+ coord_transfer as the triangle coordinates\n#are different depending o n the orientation of the basis element (ie, what the last 2 digits are ).\n#Also, the basis functions we want to use change as well, again de pending on the last 2 digits.\n#All of this hinges around the fact tha t the neighbors are ordered in a clockwise fashion, so that the\n#orie ntation of the basis element tells you which neighbors are the vertice s of which triangle.\n\n#first let's deal with the basis element being an f...\nif basis[2]=f then\n if basis[1][nops(basis[1])-1..nops(b asis[1])] = [0,1] then \n triangle1:=get_triangle(S,basis,m_basis ,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n for j from 1 to nops(vmh) do\n\n T[coord_transfer(vmh[j],triangle1,bas is[1],S[basis[1],1],S[basis[1],2]),funcno]:= T[coord_transfer(vmh[j],t riangle1,basis[1],S[basis[1],1],S[basis[1],2]),funcno] + C[i] * (T[vmh [j],harm_funcno+1]);\n\n T[coord_transfer(vmh[j],triangle2,S[basi s[1],4],basis[1],S[basis[1],3]),funcno]:= T[coord_transfer(vmh[j],tria ngle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno] + C[i] * (T[vmh[j] ,harm_funcno+5]);\n \n od;\n T[basis[1],funcno]:=T[basis [1],funcno]-C[i]; #we added in the value of the basis here twice, so s ubtract one away\n\n fi;\n\n if basis[1][nops(basis[1])-1..nops( basis[1])] = [0,2] then\n triangle1:=get_triangle(S,basis,m_basis ,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n\n for j from 1 to nops(vmh) do\n\n 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[v mh[j],harm_funcno+9]);\n\n T[coord_transfer(vmh[j],triangle2,basi s[1],S[basis[1],3],S[basis[1],4]),funcno]:= T[coord_transfer(vmh[j],tr iangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno] + C[i] * (T[vmh[ j],harm_funcno+1]);\n\n od;\n T[basis[1],funcno]:=T[basis[1] ,funcno]-C[i];\n fi; \n\n if basis[1][nops(basis[1])-1..nops(bas is[1])] = [1,2] then\n triangle1:=get_triangle(S,basis,m_basis,1) ;\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n for j fr om 1 to nops(vmh) do\n\n T[coord_transfer(vmh[j],triangle1,S[basi s[1],2],basis[1],S[basis[1],1]),funcno]:= T[coord_transfer(vmh[j],tria ngle1,S[basis[1],2],basis[1],S[basis[1],1]),funcno] + C[i] * (T[vmh[j] ,harm_funcno+5]);\n\n T[coord_transfer(vmh[j],triangle2,S[basis[1 ],3],S[basis[1],4],basis[1]),funcno]:= T[coord_transfer(vmh[j],triangl e2,S[basis[1],3],S[basis[1],4],basis[1]),funcno] + C[i] * (T[vmh[j],ha rm_funcno+9]);\n\n od;\n T[basis[1],funcno]:=T[basis[1],func no]-C[i];\n fi; \nfi;\n\n#now let us deal with the basis element be ing a g. The situation here is identical to the f situation, except \+ \n#that the basis functions change (increase by 2), and in triangle 2, we must subtract instead of add\n#because the basis values are negati ve there...\nif basis[2]=g then\n if basis[1][nops(basis[1])-1..nop s(basis[1])] = [0,1] then \n triangle1:=get_triangle(S,basis,m_ba sis,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n fo r j from 1 to nops(vmh) do\n\n 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]);\n\n T[coord_transfer(vmh[j],triangle2,S[b asis[1],4],basis[1],S[basis[1],3]),funcno]:= T[coord_transfer(vmh[j],t riangle2,S[basis[1],4],basis[1],S[basis[1],3]),funcno] - C[i] * (T[vmh [j],harm_funcno+7]);\n od;\n\n fi;\n\n if basis[1][nops(bas is[1])-1..nops(basis[1])] = [0,2] then\n triangle1:=get_triangle( S,basis,m_basis,1);\n triangle2:=get_triangle(S,basis,m_basis,2); \n\n\n for j from 1 to nops(vmh) do\n\n 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]);\n\n T[coord_transfer(vmh[j ],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno]:= T[coord_tr ansfer(vmh[j],triangle2,basis[1],S[basis[1],3],S[basis[1],4]),funcno] \+ - C[i] * (T[vmh[j],harm_funcno+3]);\n od;\n\n fi; \n\n if b asis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then\n triangle 1:=get_triangle(S,basis,m_basis,1);\n triangle2:=get_triangle(S,b asis,m_basis,2);\n\n\n for j from 1 to nops(vmh) do\n\n T[co ord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),fu ncno]:= T[coord_transfer(vmh[j],triangle1,S[basis[1],2],basis[1],S[bas is[1],1]),funcno] + C[i] * (T[vmh[j],harm_funcno+7]);\n\n T[coord _transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],basis[1]),funcn o]:= T[coord_transfer(vmh[j],triangle2,S[basis[1],3],S[basis[1],4],bas is[1]),funcno] - C[i] * (T[vmh[j],harm_funcno+11]);\n od;\n\n \+ fi; \n\nfi;\n\nfi;\n\nod;\n\nend;\n \n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%4make_pretty_pictureGf*6*%\"TG%\"SG%(m_basisG%'m_f uncG%#vmG%'funcnoG%,harm_funcnoG%\"CG61%\"iG%\"jG%$vm2G%\"kG%$vmhG%*tr iangle1G%*triangle2G%&basisG%%add0G%%add1G%%add2G%,which_harm1G%,which _harm2G%'coord1G%'coord2G6\"F?C&>8&-%)verticesG6#9'?(8$\"\"\"FI-%%nops G6#FB%%trueG>&9$6$&FB6#FH9)\"\"!>8(-FD6#,&FFFI9&!\"\"?(FHFIFI)\"\"$,&F enFIFIFIFMC&>8+-%*get_basisG6$9(FH-%'printfG6$%$%d~GFH@$/-%$modG6$FH\" #:FI-Fco6#%\"|+G@%/-FK6#&F]o6#FIFIC&@$/Fcp7#FUC'>8)7#-%$seqG6$FU/8';FI Fen>8,Fcp>8-&9%6$FcpFI>8.&Fhq6$Fcp\"\"#>8/Fin@$/Fcp7#FIC'>F[q7#-F^q6$F IF`q>FdqF\\r>FfqFcp>F[rFgq>F`r\"\"(@$/Fcp7#F^rC'>F[q7#-F^q6$F^rF`q>Fdq Fgq>FfqF\\r>F[rFcp>F`r\"#6?(8%FIFI-FK6#FWFMC$>81-%/coord_transferG6'&F W6#F\\tF[qFdqFfqF[r>&FP6$FatFT,&FhtFI*&&9+FSFI&FP6$Fet,&9*FIF`rFIFIFIC $@$/&F]o6#F^r%\"fGC%@$/&Fcp6#;,&FapFIFIFfnFap7$FUFIC&>F[q-%-get_triang leG6&FhqF]oFenFI>8*-Fcv6&FhqF]oFenF^r?(F\\tFIFIF]tFMC$>&FP6$-Fct6'FetF [qFcpFgqF\\rFT,&F\\wFI*&F\\uFI&FP6$Fet,&FauFIFIFIFIFI>&FP6$-Fct6'FetFf v&Fhq6$Fcp\"\"%Fcp&Fhq6$FcpFinFT,&FfwFI*&F\\uFI&FP6$Fet,&FauFI\"\"&FIF IFI>&FP6$FcpFT,&FfxFIF\\uFfn@$/F[v7$FUF^rC&>F[qFbv>FfvFgv?(F\\tFIFIF]t FMC$>&FP6$-Fct6'FetF[qFgqF\\rFcpFT,&FbyFI*&F\\uFI&FP6$Fet,&FauFI\"\"*F IFIFI>&FP6$-Fct6'FetFfvFcpF]xFjwFT,&F]zFIFawFI>FfxFhx@$/F[v7$FIF^rC&>F [qFbv>FfvFgv?(F\\tFIFIF]tFMC$>&FP6$-Fct6'FetF[qF\\rFcpFgqFT,&F\\[lFIF` xFI>&FP6$-Fct6'FetFfvF]xFjwFcpFT,&Fb[lFIFgyFI>FfxFhx@$/Feu%\"gGC%@$Fju C%>F[qFbv>FfvFgv?(F\\tFIFIF]tFMC$>F\\w,&F\\wFI*&F\\uFI&FP6$Fet,&FauFIF inFIFIFI>Ffw,&FfwFI*&F\\uFI&FP6$Fet,&FauFIF]sFIFIFfn@$FjxC%>F[qFbv>Ffv Fgv?(F\\tFIFIF]tFMC$>Fby,&FbyFI*&F\\uFI&FP6$Fet,&FauFIFjsFIFIFI>F]z,&F ]zFIFd\\lFfn@$FdzC%>F[qFbv>FfvFgv?(F\\tFIFIF]tFMC$>F\\[l,&F\\[lFIFj\\l FI>Fb[l,&Fb[lFIFf]lFfnF?F?F?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 662 "#This function generates the matrix for non-constant q (ie, int egral(q*basisi*basisj) ).\n#Again, the meat of the functionality is in int_vs_basis12, which performs the actual integration\n\ngen_q_matrix :=proc(T,S,m_basis,m_func,vm,funcno,harm_funcno)\n\nlocal i,j,q_matrix ,q_val;\n\n\nq_matrix:=array(sparse,1..3^(m_basis+1),1..3^(m_basis+1)) ;\n\n\nfor i from 1 to 3^(m_basis+1) do\nprintf(`%d `,i);\nif ((i mod \+ 15)=1) then\n printf(`\\n`): fi:\n\n for j from 1 to i do\n q_ val:=int_vs_basis12(T,S,m_basis,m_func,get_basis(vm,i),get_basis(vm,j) ,funcno,harm_funcno);\n if q_val <> 0 then q_matrix[i,j]:=q_val; \+ q_matrix[j,i]:=q_val;fi;\n od;\nod;\n\nq_matrix;\nend; \n" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%-gen_q_matrixGf*6)%\"TG%\"SG%(m_basi sG%'m_funcG%#vmG%'funcnoG%,harm_funcnoG6&%\"iG%\"jG%)q_matrixG%&q_valG 6\"F3C%>8&-%&arrayG6%%'sparseG;\"\"\")\"\"$,&9&F 8'-%/int_vs_basis12G6*9$9%F@9'-%*get_basisG6$9(FB-Fhn6$FjnFS9)9*@$0FV \"\"!C$>&F66$FBFSFV>&F66$FSFBFVF6F3F3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 859 "#This function creates the laplacian of the function defined by the spline\n#projection C. Notice that the values at the \+ Vm_basis are completely wrong\n#because the laplacian is discontinuous at those points. See the posteriori\n#error function in error_funcs. mws for a method of storing the discontinuous\n#values of the laplacia n at the basis vertices. Notice that the way this\n#function works is just by calling make_pretty picture, but instead of sending\n#the bih armonic basis elements, we send it the laplacians of the biharmonic\n# basis elements (hence harm_funcno+12).\n\nmake_pretty_picture_lap:=pro c(T,S,m_basis,m_func,vm,funcno,harm_funcno,C)\n\nlocal i,vmfunc;\n\nvm func:=vertices(m_func);\n\nmake_pretty_picture(T,S,m_basis,m_func,vm,f uncno,harm_funcno+12,C);\n\nfor i from 1 to nops(vmfunc) do\n T[vmfu nc[i],funcno]:=T[vmfunc[i],funcno]*5^m_basis;\nod;\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%8make_pretty_picture_lapGf*6*%\"TG%\"SG%(m _basisG%'m_funcG%#vmG%'funcnoG%,harm_funcnoG%\"CG6$%\"iG%'vmfuncG6\"F2 C%>8%-%)verticesG6#9'-%4make_pretty_pictureG6*9$9%9&F99(9),&9*\"\"\"\" #7FD9+?(8$FDFD-%%nopsG6#F5%%trueG>&F=6$&F56#FHFA*&FNFD)\"\"&F?FDF2F2F2 " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 1249 "#This function returns a table R which conta ins the values of p and q as defined\n#by Mike Usher for determining t he values of n-harmonic extensions. The values\n#a and b are values o f inner products. The function basically implements the \n#formula fo r determining extensions based on an iterative method involving inner \n#products. See \"Splines on Fractals\" from more information (ie, d on't ask me,\n#I just work here...).\n\nps_n_qs:=proc(n)\n local i,l ,j,R,m;\n R:=table();\n\n R[p,0]:=2/5;\n R[q,0]:=1/5;\n R[a,0] :=2;\n R[b,0]:=-1;\n\n R[p,1]:=-3/25;\n R[q,1]:=-7/75;\n R[a,1 ]:=7/45;\n R[b,1]:=4/45;\n\n for i from 1 to n-1 do\n\n 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);\n 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);\n \n 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));\n 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));\n\n\n od;\n R;\nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%(ps_n_qsGf*6#%\"nG6'%\"iG%\"lG%\"jG%\"RG%\"mG6\" F.C->8'-%&tableGF.>&F16$%\"pG\"\"!#\"\"#\"\"&>&F16$%\"qGF8#\"\"\"F;>&F 16$%\"aGF8F:>&F16$%\"bGF8!\"\">&F16$F7FA#!\"$\"#D>&F16$F?FA#!\"(\"#v>& F16$FEFA#\"\"(\"#X>&F16$FIFA#\"\"%Ffn?(8$FAFA,&9$FAFAFJ%%trueGC&>&F16$ FE,&F]oFAFAFA,$*&,()F;,$F]oF:FA*&\"\"'FA)F;,&F]oFAFAFJFAFJF@FAFJ-%$sum G6$*&,&*&,&)F;F]oF[o#F[o\"#:FJFA&F16$F7,(F]oFAFAFA8&FJFAFA*&,&FfpFA#\" #8FhpFAFA&F16$F?F[qFAFAFA,&&F16$FEF\\qFA*&F:FA&F16$FIF\\qFAFAFA/F\\q;F AF]oFA#F:Fhp>&F16$FIFeo,$*&FhoFJ-F`p6$*&,&*&,&Ffp\"\"$#F`qFhpFJFAFipFA FA*&,&FfpF:#\"#9FhpFJFAFaqFAFAFAFcqFAFiqFAF[r>&F16$F7Feo,(Fco#!\"*FP*& #\"#=FPFAF]rFAFJ*&#FAF;FA-F`p6$,&*&,&&F16$F7,(F]oFAFAFA8(FJF[o&F16$F?F _tFAFA&F16$FEF`tFAFA*&,&F]tFgr*&F:FAFatFAFAFA&F16$FIF`tFAFA/F`tFjqFAFJ >&F16$F?Feo,(Fco#FUFP*&#F\\sFPFAF]rFAFJ*&#FAF;FA-F`p6$,&*&,&F]tF:*&Fgr FAFatFAFAFAFctFAFA*&F\\tFAFhtFAFAFjtFAFJF1F.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 140 "#As long as I'm here, why don't we just calc ulate this thing out to 5\n#so that we can do sex-harmonic functions ( Yeah, baby!)\nR:=ps_n_qs(5);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"RG %\"RG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "eval(R):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4276 "#This function creates a n -harmonic function, using the values for harmonic extensions calculate d\n#in R. It uses the recursive method of calculating the n-harmonic \+ extension based on Mike\n#Usher's work. In order to understand how th e recursion addressing works, I recommend you look\n#at the documentat ion of make_points and gen_points in SG.mws, where I fully explain how the \n#recursion works (in terms of the addressing). Basically, afte r that, this function\n#is quite straightforward; it simply calculates the next level based on the values of the previous\n#level, as well a s the necessary laplacians. This procedure is called by the wrapper f unction\n#n_harm.\n\n\nnharm:=proc(T,R,n,m,add0,add1,add2,funcno)\n \+ local i,j;\n if m <> 0 then\n #print(funcno);\n #print(add 0,add1,add2);\n\n if add1[nops(add1)-1..nops(add1)] = [0,1] then \n #This loop is because we need to calculate up to n-1 lapla cians of the function\n #at this level as well...\n \+ for i from 0 to n-1 do\n #the sums here determine the valu e at the next level based on the values of \n #various lap lacians and the such before. Again, see Mike Usher's work for a\n \+ #derivation of this formula. \n T[[op(add1[1..no ps(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[a dd2,funcno+n-1-i+j] ),j=0..i); \n 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,fun cno+n-1-i+j] ),j=0..i); \n 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); \n\n od;\n nharm(T,R,n,m-1,add0,[o p(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],func no);\n 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);\n 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);\n fi;\n if add1[nops(add1)-1..nops(add1)] = [1, 2] then\n for i from 0 to n-1 do\n 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); \n T[[op(add1[1..nops(a dd1)-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);\n 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);\n od;\n nharm(T,R,n,m-1,add0,[op(a dd1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno) ;\n nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[o p(add1[1..nops(add1)-2]),1,1,2],funcno);\n nharm(T,R,n,m-1,ad d2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2] ,funcno);\n fi;\n if add1[nops(add1)-1..nops(add1)] = [0,2] \+ then\n for i from 0 to n-1 do\n T[[op(add1[1..nop s(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[ad d2,funcno+n-1-i+j] ),j=0..i); \n 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,func no+n-1-i+j] ),j=0..i);\n 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,fun cno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+ j] ),j=0..i);\n od;\n 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);\n \+ nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(ad d1[1..nops(add1)-2]),2,0,2],funcno);\n 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],fun cno);\n fi;\n fi; \nend;\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%&nharmGf*6*%\"TG%\"RG%\"nG%\"mG%%add0G%%add1G%%add2G%'funcnoG6$% \"iG%\"jG6\"F2@$09'\"\"!C%@$/&9)6#;,&-%%nopsG6#F;\"\"\"FB!\"\"F?7$F6FB C&?(8$F6FB,&9&FBFBFC%%trueGC%>&9$6$7&-%#opG6#&F;6#;FB,&F?FB\"\"#FCF6F6 FB,*9+FBFIFBFBFCFGFC-%$sumG6$*&)\"\"&*&8%FBF?FBFC,&*&&9%6$%\"pGF\\oFB, &&FN6$9(,,FZFBFIFBFBFCFGFCF\\oFBFB&FN6$F;FgoFBFBFB*&&F`o6$%\"qGF\\oFB& FN6$9*FgoFBFBFB/F\\o;F6FG>&FN6$7&FQF6FBFXFY-Ffn6$*&FinFC,&*&F_oFB,&Fho FBF^pFBFBFB*&F[pFBFdoFBFBFBFap>&FN6$7&FQF6F6FXFY-Ffn6$*&FinFC,&*&F_oFB ,&FdoFBF^pFBFBFB*&F[pFBFhoFBFBFBFap-F$6*FNF`oFI,&F5FBFBFCFfoFPFaqFZ-F$ 6*FNF`oFIF[rF;FfpFPFZ-F$6*FNF`oFIF[rF`pFaqFfpFZ@$/F:7$FBFXC&?(FGF6FBFH FJC%>&FN6$7&FQFBFBFXFYFen>&FN6$7&FQFBF6FXFYFgp>&FN6$7&FQFBF6FBFYFbq-F$ 6*FNF`oFIF[rFfoFirFasFZ-F$6*FNF`oFIF[rF;F]sFirFZ-F$6*FNF`oFIF[rF`pFasF ]sFZ@$/F:7$F6FXC&?(FGF6FBFHFJC%>&FN6$7&FQFXF6FXFYFen>&FN6$7&FQFXF6FBFY Fgp>&FN6$7&FQFXFBFXFYFbq-F$6*FNF`oFIF[rFfoFatFitFZ-F$6*FNF`oFIF[rF;Fet FatFZ-F$6*FNF`oFIF[rF`pFitFetFZF2F2F2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1707 "#So this is basically a wrapper function around the function nharm. This function is \n#necessary because nharm won't wo rk until V1, due to idiosyncracies in the labeling system.\n#This func tion basically computes the first n-harmonic extension, and then sends the rest out to\n#nharm to compute. n stands for n-harmonic, and R i s the table of harmonic extension values\n#calculated using ps_and_qs. In order to make this function work, you have to put in the values\n #of the function at the boundaries in funcno, as well as the values of the laplacians of the\n#function in funcno+1, L^2 of the function in \+ funcno+2... and so forth. This uniquely determines\n#an n-harmonic fu nction (ie, you need to specify boundary values of L^(0->n-1)). Notic e that \n#the algorithm developed required that the function as well a s all of it's laplacians be \n#defined out to all levels. This has th e added bonus that while calculating the function, you have\n#calculat ed the laplacian of it in funcno+1, L^2(f) in funcno+2... Fun, isn't \+ it?\n\nn_harm:=proc(T,R,n,m,add0,add1,add2,funcno)\n\n local i,j;\n \n\n for i from 0 to n-1 do\n 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);\n 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);\n T[[0,2],funcn o+n-1-i]:=sum((1/(5^(j*1)))*(R[p,j]*(T[add0,funcno+n-1-i+j]+T[add2,fun cno+n-1-i+j])+R[q,j]*T[add1,funcno+n-1-i+j]),j=0..i); \n od;\n\n\n \+ nharm(T,R,n,m-1,[0],[0,1],[0,2],funcno);\n nharm(T,R,n,m-1,[1],[1,2 ],[0,1],funcno);\n nharm(T,R,n,m-1,[2],[0,2],[1,2],funcno); \n\n \n\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%'n_harmGf*6*%\"TG%\"RG% \"nG%\"mG%%add0G%%add1G%%add2G%'funcnoG6$%\"iG%\"jG6\"F2C&?(8$\"\"!\" \"\",&9&F7F7!\"\"%%trueGC%>&9$6$7$F6F7,*9+F7F9F7F7F:F5F:-%$sumG6$*&)\" \"&*&8%F7-%%nopsG6#9)F7F:,&*&&9%6$%\"pGFKF7,&&F?6$9(,,FCF7F9F7F7F:F5F: FKF7F7&F?6$FOFZF7F7F7*&&FS6$%\"qGFKF7&F?6$9*FZF7F7F7/FK;F6F5>&F?6$7$F7 \"\"#FB-FE6$*&)FIFKF:,&*&FRF7,&FenF7F[oF7F7F7*&FhnF7FWF7F7F7F^o>&F?6$7 $F6FdoFB-FE6$*&FhoF:,&*&FRF7,&FWF7F[oF7F7F7*&FhnF7FenF7F7F7F^o-%&nharm G6*F?FSF9,&9'F7F7F:7#F6FAF`pFC-Fip6*F?FSF9F[q7#F7FcoFAFC-Fip6*F?FSF9F[ q7#FdoF`pFcoFCF2F2F2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3580 "# This function recursively generates a harmonic function. It is a much simpler algorithm than\n#the n-harmonic algorithm, so the addressing \+ business is much more transparent here. Basically,\n#it implements th e 2/5, 2/5, 1/5 rule for generating the harmonic extension, and recurs ively \n#generates values all the way down into the gasket. The recur sion is a little bit complicated,\n#because each step of the way, you \+ need to know the orientation of the triangle in order to know\n#what t ype of points are assigned what values. Basically, we used a rotation al method of \n#determining what the orientation is: we made it so tha t the 0th address always was a vertex of the\n#previous triangle (ie, \+ the 0th address of the first triangle is [1]...). This way, we always \n#which address to assign to, based on the last digits of address 1 ( if address 1 is [1,2], then \n#we know the triangle is a 1 triangle. \+ Thus, you know what addresses to assign what values.\n#Anyway, it work s. If you want a somewhat more transparent example, see make_points a nd gen_points\n#in SG.mws. The recursion method used there is exactly the same as used here, and should be easily\n#generalizable to other \+ fractals (ie, pentagasket...).\n\nharm:=proc(T,m,add0,add1,add2,funcno )\n if m <> 0 then\n #print(funcno);\n #print(add0,add1,ad d2);\n\n if add1[nops(add1)-1..nops(add1)] = [0,1] then\n \+ 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]; \n T[[op(add1 [1..nops(add1)-2]),0,1,2],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[ad d2,funcno] + 1/5 * T[add0,funcno];\n 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];\n harm(T,m-1,add0,[op(add1[1..nops(add1) -2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno);\n harm (T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2 ]),0,0,1],funcno);\n harm(T,m-1,add2,[op(add1[1..nops(add1)-2 ]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2],funcno);\n fi;\n \+ if add1[nops(add1)-1..nops(add1)] = [1,2] then\n T[[op(add1[ 1..nops(add1)-2]),1,1,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add 1,funcno] + 1/5 * T[add2,funcno]; \n 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];\n T[[op(add1[1..nops(add1)-2]),1,0,1],f uncno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,f uncno];\n harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[ op(add1[1..nops(add1)-2]),1,0,1],funcno);\n harm(T,m-1,add1,[ op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],fun cno);\n harm(T,m-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op (add1[1..nops(add1)-2]),1,0,2],funcno);\n fi;\n if add1[nops (add1)-1..nops(add1)] = [0,2] then\n 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]; \n 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];\n 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];\n \+ harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..no ps(add1)-2]),2,1,2],funcno);\n harm(T,m-1,add1,[op(add1[1..no ps(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno);\n \+ harm(T,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops (add1)-2]),2,0,1],funcno);\n fi;\n fi; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%%harmGf*6(%\"TG%\"mG%%add0G%%add1G%%add2G%'funcn oG6\"F-F-@$09%\"\"!C%@$/&9'6#;,&-%%nopsG6#F6\"\"\"F=!\"\"F:7$F1F=C(>&9 $6$7&-%#opG6#&F66#;F=,&F:F=\"\"#F>F1F1F=9),(&FC6$9&FN#FM\"\"&*&FSF=&FC 6$F6FNF=F=*&#F=FTF=&FC6$9(FNF=F=>&FC6$7&FFF1F=FMFN,(FVFS*&FSF=FZF=F=*& FYF=FPF=F=>&FC6$7&FFF1F1FMFN,(FPFS*&FSF=FZF=F=*&FYF=FVF=F=-F$6(FC,&F0F =F=F>FRFEFaoFN-F$6(FCFgoF6FjnFEFN-F$6(FCFgoFfnFaoFjnFN@$/F57$F=FMC(>&F C6$7&FFF=F=FMFNFO>&FC6$7&FFF=F1FMFNF[o>&FC6$7&FFF=F1F=FNFbo-F$6(FCFgoF RFcpF[qFN-F$6(FCFgoF6FgpFcpFN-F$6(FCFgoFfnF[qFgpFN@$/F57$F1FMC(>&FC6$7 &FFFMF1FMFNFO>&FC6$7&FFFMF1F=FNF[o>&FC6$7&FFFMF=FMFNFbo-F$6(FCFgoFRFiq FarFN-F$6(FCFgoF6F]rFiqFN-F$6(FCFgoFfnFarF]rFNF-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 740 "#This function is basically a wrapper fo r harm\n#It generates a harmonic function with val0,val1,val2 on the b oundaries.\n#This function calculates the first harmonic extension, an d then calls harm\n#on the 3 parts to generate the rest of the harmoni c function. It can't do it directly due\n#to idiosyncracies in the no menclature of the V0 points.\n\nharmonic:=proc(T,m,val0,val1,val2,func no)\nT[[0],funcno]:=val0;\nT[[1],funcno]:=val1;\nT[[2],funcno]:=val2; \n\nT[[0,1],funcno]:=2/5 * val0 + 2/5 * val1 + 1/5 * val2;\nT[[0,2],fu ncno]:=2/5 * val0 + 2/5 * val2 + 1/5 * val1;\nT[[1,2],funcno]:=2/5 * v al1 + 2/5 * val2 + 1/5 * val0;\n\nharm(T,m-1,[2],[0,2], [1,2], funcno) ;\nharm(T,m-1,[0],[0,1], [0,2], funcno);\nharm(T,m-1,[1],[1,2], [0,1], funcno);\nend;\n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)harmonic Gf*6(%\"TG%\"mG%%val0G%%val1G%%val2G%'funcnoG6\"F-F-C+>&9$6$7#\"\"!9)9 &>&F16$7#\"\"\"F59'>&F16$7#\"\"#F59(>&F16$7$F4F;F5,(F6#FA\"\"&*&FHF;F< F;F;*&#F;FIF;FBF;F;>&F16$7$F4FAF5,(F6FH*&FHF;FBF;F;*&FLF;F&F16$7 $F;FAF5,(F " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 403 "# This procedure creates a function in out_func which is u-Pu\n#It also \+ creates a function in out_func + 1 which is (u-Pu)^2 (used for the L2_ norm)\nerror_fns:=proc(T,m,func_T,func_A,out_func)\n local i,vm;\n \+ vm:=vertices(m);\n\n for i from 1 to nops(vm) do\n T[vm[i],out _func] := T[vm[i],func_T] - T[vm[i],func_A];\n T[vm[i],out_func+1 ]:= (T[vm[i],func_T] - T[vm[i],func_A])^2;\n od;\n \nend; " }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%*error_fnsGf*6'%\"TG%\"mG%'func_TG%' func_AG%)out_funcG6$%\"iG%#vmG6\"F/C$>8%-%)verticesG6#9%?(8$\"\"\"F9-% %nopsG6#F2%%trueGC$>&9$6$&F26#F89(,&&FA6$FC9&F9&FA6$FC9'!\"\">&FA6$FC, &FEF9F9F9*$)FF\"\"#F9F/F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 322 "#This function takes the L2_norn of error_func, assuming that the square of error_func\n#is in error_func+1, as it would be if generate d by error_fns.\nL2_norm:=proc(T,m,error_func) #!!!!NOTE, ASSUMES THE \+ SQUARE OF THE FUNCTION IS IN ERROR_FUNC+1!!\n\nlocal err,i;\n\n\nerr:= simp_rule(T,m,error_func+1);\n \nsqrt(err);\nend;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(L2_normGf*6%%\"TG%\"mG%+error_funcG6$%$errG% \"iG6\"F-C$>8$-%*simp_ruleG6%9$9%,&9&\"\"\"F8F8-%%sqrtG6#F0F-F-F-" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 538 "#This procedure finds the s quare root of the energy of func_in.\n#Note that this function require s that make_neighbors has been run\n#on T out to level m...\n\nget_roo t_energy:=proc(T,m,func_in)\n local err, Vm,V0,j, i;\n\n err:=0;\n\n V m:=vertices(m);\n V0:=vertices(0);\n \n for i from 4 to nops(Vm) do \n err := err + (5/3)^m * (add( (T[T[Vm[i],j],func_in] - T[Vm[i] ,func_in])^2, j = 1..4))/2;\n od:\n for i from 1 to 3 do\n err := \+ err + (5/3)^m * (add( (T[T[Vm[i],j],func_in] - T[Vm[i],func_in])^2, j \+ = 1..2))/2;\n od; \n\nsqrt(err);\n \nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%0get_root_energyGf*6%%\"TG%\"mG%(func_inG6'%$errG%#Vm G%#V0G%\"jG%\"iG6\"F0C(>8$\"\"!>8%-%)verticesG6#9%>8&-F86#F4?(8(\"\"% \"\"\"-%%nopsG6#F6%%trueG>F3,&F3FB*(#FB\"\"#FB)#\"\"&\"\"$F:FB-%$addG6 $*$),&&9$6$&FW6$&F66#F@8'9&FB&FW6$FenFhn!\"\"FKFB/Fgn;FBFAFBFB?(F@FBFB FOFF>F3,&F3FB*(FJFBFLFB-FQ6$FS/Fgn;FBFKFBFB-%%sqrtG6#F3F0F0F0" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 295 "#This procedure finds the m aximum of a function (used to find max error...)\n\nfind_max:=proc(T,m ,func_in)\n local i,abs_diff,temp,Vm;\n Vm:=vertices(m);\n temp: =0;\n\nfor i from 1 to nops(Vm) do\n abs_diff := abs(T[Vm[i],func_in] ):\n if abs_diff > temp then temp := abs_diff; fi; od:\n\ntemp;\nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)find_maxGf*6%%\"TG%\"mG%(func_i nG6&%\"iG%)abs_diffG%%tempG%#VmG6\"F/C&>8'-%)verticesG6#9%>8&\"\"!?(8$ \"\"\"F<-%%nopsG6#F2%%trueGC$>8%-%$absG6#&9$6$&F26#F;9&@$2F8FC>F8FCF8F /F/F/" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7085 "#So this functio n estimates the posteriori error with q nonconstant (if q is constant, just stick in a \n#function with constant values). It calculates L2_ norm(-Lp+qp-f). This function requires you send in C,\n#C being the c oefficients of the spline representation of p. This allows us to calc ulate a theoretical\n#laplacian, since we know the exact values of the laplacian of the basis functions (they're biharmonic...).\n#Also, l_f uncno is just some random function number where you want the laplacian of the function stored.\n#Note that the laplacian stored in l_funcno \+ has incorrect values at the vertices of the spline basis, because\n#th e laplacian is discontinuous at those points. The solution to this pr oblem is documented within this procedure.\n\nposteriori_error_q_nonc: =proc(T,S,m_basis,m_func,C,p_funcno,q_funcno,f_funcno,harm_funcno,l_fu ncno)\n\nlocal vm,Vm,fact,V0,Vother,Vother_e,i,basis,triangle1,triangl e2,addfuncs,addresses1,addresses2,tris,total_int,Vmsum1,Vothersum1,V0s um1,j;\n\nvm:=vertices(m_basis);\n\nfact:=1/(6*3^m_func); #This is ju st the factor from simpson's rule.\n\n\n#Here, we are just getting the vertices of V0, Vm, and (in Vother) everything in between.\nV0:=verti ces_m(0);\nif m_basis = 1 then Vother:=[] else Vother_e:=vertices(m_fu nc-m_basis-1);Vother:=Vother_e[4..nops(Vother_e)]; fi;\nVm:=vertices_m (m_func-m_basis);\n\n#Here, we are getting all the triangles in our sp line basis and storing them in tris\nif m_basis=1 then tris:=[[0],[1], [2]] else tris:=p([0,1,2],m_basis) fi;\n\n#Okay, so here's the general idea for dealing with the discontinuous laplacian:\n#We will store th e addresses of the vertices of the triangles in tris in S[tris[i],10-1 2],\n#which is performed by get_triangle_adds. Then, in S[tris[i],7-9 ], we store the values of the\n#laplacians at the boundaries of tris[i ]. Thus, we have stored both the values of the laplacian\n#at each sp line basis vertex: once for each triangle. Then we will piecewise int egrate our function\n#by summing the integrals over all the tris[i], u sing the boundary values which correspond to that\n#particular triangl e. Pretty tricky, huh?\n\nfor 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...\n \nget_triangle_adds(S,m_basis); #get's the addresses of the boundarie s of tris and stores them as indicated\n\n#Here, we generate the lapla cian for all values EXCEPT for the spline vertices\nprint(`making lapl acian...`);\nmake_pretty_picture_lap(T,S,m_basis,m_func,vm,l_funcno,ha rm_funcno,C);\nprint(`finished laplacian...`);\n\nfor i from 1 to 3^(m _basis+1) do #First, let's put the values of the laplacians at the spl ine basis vertices (DISCONTINUOUS) in the table S... \n\n basis:=get _basis(vm,i);\n if basis[1]<>[0] and basis[1]<>[1] and basis[1]<>[2] then #ie, not a boundary point\n triangle1:=get_triangle(S,basis ,m_basis,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n \n #Here, we determine which basis functions we need, which again \+ \n #depend entirely on the last two digits of the basis elt (ie, \+ the orientation)\n if last_digits(basis[1])=[0,1] then addfuncs:= [1,5,3,7] fi;\n if last_digits(basis[1])=[1,2] then addfuncs:=[5, 9,7,11] fi;\n if last_digits(basis[1])=[0,2] then addfuncs:=[9,1, 11,3] fi;\n\n #And then just start adding the laplacian values in .\n #Notice that the +12 after harm_funcno means that we are\n \+ #accessing the laplacians of the basis functions, rather than the\n #basis functions themselves. \n if basis[2]=f then \n \+ S[triangle1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],harm_funcno+12+addf uncs[1]]*5^m_basis;\n S[triangle1,7+1]:=S[triangle1,7+1]+C[i]* T[[1],harm_funcno+12+addfuncs[1]]*5^m_basis;\n S[triangle1,7+2 ]:=S[triangle1,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs[1]]*5^m_basis; \n \n S[triangle2,7+0]:=S[triangle2,7+0]+C[i]*T[[0],harm_fun cno+12+addfuncs[2]]*5^m_basis;\n S[triangle2,7+1]:=S[triangle2 ,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs[2]]*5^m_basis;\n S[tr iangle2,7+2]:=S[triangle2,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs[2]]* 5^m_basis;\n else #DON'T FORGET THE NEGATIVE SIGN FOR G-TYPE BAS IS ELTS!!!\n S[triangle1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],har m_funcno+12+addfuncs[3]]*5^m_basis;\n S[triangle1,7+1]:=S[tria ngle1,7+1]+C[i]*T[[1],harm_funcno+12+addfuncs[3]]*5^m_basis;\n \+ S[triangle1,7+2]:=S[triangle1,7+2]+C[i]*T[[2],harm_funcno+12+addfuncs [3]]*5^m_basis;\n \n S[triangle2,7+0]:=S[triangle2,7+0]-C[i] *T[[0],harm_funcno+12+addfuncs[4]]*5^m_basis;\n S[triangle2,7+ 1]:=S[triangle2,7+1]-C[i]*T[[1],harm_funcno+12+addfuncs[4]]*5^m_basis; \n S[triangle2,7+2]:=S[triangle2,7+2]-C[i]*T[[2],harm_funcno+1 2+addfuncs[4]]*5^m_basis;\n fi;\n else #If we have a basis elt \+ which is on the boundary\n triangle1:=[seq(op(basis[1]),j=1..m_ba sis)]; #So the triangle we're dealing with is just a sequence of \n \+ #0's, 1's, or 2's, d epending on the value of the boundary\n if basis[1]=[0] then addf uncs:=3 fi;\n if basis[1]=[1] then addfuncs:=7 fi;\n if basi s[1]=[2] then addfuncs:=11 fi;\n\n #We know it's a g-type basis e lt, and that the values will always be positive...\n\n S[triangle 1,7+0]:=S[triangle1,7+0]+C[i]*T[[0],harm_funcno+12+addfuncs]*5^m_basis ;\n S[triangle1,7+1]:=S[triangle1,7+1]+C[i]*T[[1],harm_funcno+12+ addfuncs]*5^m_basis;\n S[triangle1,7+2]:=S[triangle1,7+2]+C[i]*T[ [2],harm_funcno+12+addfuncs]*5^m_basis;\n fi;\n \n\nod;\n\nprint (`integrating difference fn... `);\n\n#So now we're ready to integrate each triangle...\n\n\ntotal_int:=0;\n\n\nfor j from 1 to nops(tris) d o\n\nprintf(`%d `,j);\nif ((j mod 15)=1) then\n printf(`\\n`): fi:\n \n#So basically, we call func_val to get the value of the at the verte x defined by V0 (Vother,Vm) [i], when considered as\n#a vertex of the \+ subtriangle tris[j]. This basically encapsulates the functionality of coord_transfer. So we just\n#add up -Lp+qp-f. Notice that in the V0 sum1, instead of using -func_val(laplacian), we use -S[tris[j],7+i-1], which is\n#the value of the laplacian at the boundary when considerin g the triangle tris[j]. Thus, we use the value of the \n#laplacian ap propriate to each triangle, thereby dealing with the multiple values o f the laplacian.\n\n\n 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) );\n#print(`sum1`); \+ \n if m_basis <> 1 then Vothersum1:=add( ( -func_val(T,S,tris[j],Vot her[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;\n\n#print(`sum2`); \n \+ Vmsum1:=add( ( -func_val(T,S,tris[j],Vm[i],l_funcno)+func_val(T,S,tr is[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));\n#print(`sum3`); \n\n \+ total_int:=total_int + fact*(5*Vmsum1 + 2*Vothersum1 + V0sum1);\nod; \n\nsqrt(total_int); #Don't forget the square root!\nend;\n\n" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%8posteriori_error_q_noncGf*6,%\"TG% \"SG%(m_basisG%'m_funcG%\"CG%)p_funcnoG%)q_funcnoG%)f_funcnoG%,harm_fu ncnoG%)l_funcnoG65%#vmG%#VmG%%factG%#V0G%'VotherG%)Vother_eG%\"iG%&bas isG%*triangle1G%*triangle2G%)addfuncsG%+addresses1G%+addresses2G%%tris G%*total_intG%'Vmsum1G%+Vothersum1G%'V0sum1G%\"jG6\"FEC2>8$-%)vertices G6#9&>8&,$*&\"\"\"FQ)\"\"$9'!\"\"#FQ\"\"'>8'-%+vertices_mG6#\"\"!@%/FL FQ>8(7\"C$>8)-FJ6#,(FTFQFLFUFQFU>F[o&F_o6#;\"\"%-%%nopsG6#F_o>8%-Fen6# ,&FTFQFLFU@%Fin>817%7#Fgn7#FQ7#\"\"#>Fbp-%\"pG6$7%FgnFQFgpFL?(8*FQFQ-F io6#Fbp%%trueGC%>&9%6$&Fbp6#F^q\"\"(Fgn>&Feq6$Fgq\"\")Fgn>&Feq6$Fgq\" \"*Fgn-%2get_triangle_addsG6$FeqFL-%&printG6#%4making~laplacian...G-%8 make_pretty_picture_lapG6*9$FeqFLFTFH9-9,9(-Ffr6#%6finished~laplacian. ..G?(F^qFQFQ)FS,&FLFQFQFQFaqC$>8+-%*get_basisG6$FHF^q@%330&Fhs6#FQFdp0 F`tFep0F`tFfpC(>8,-%-get_triangleG6&FeqFhsFLFQ>8--Fht6&FeqFhsFLFgp@$/- %,last_digitsG6#F`t7$FgnFQ>8.7&FQ\"\"&FSFiq@$/F`u7$FQFgp>Feu7&FguFarFi q\"#6@$/F`u7$FgnFgp>Feu7&FarFQF]vFS@%/&Fhs6#Fgp%\"fGC(>&Feq6$FftFiq,&F jvFQ*(&F_sFhqFQ&F\\s6$Fdp,(F^sFQ\"#7FQ&FeuFatFQFQ)FguFLFQFQ>&Feq6$FftF ]r,&FfwFQ*(F^wFQ&F\\s6$FepFawFQFdwFQFQ>&Feq6$FftFar,&F]xFQ*(F^wFQ&F\\s 6$FfpFawFQFdwFQFQ>&Feq6$F[uFiq,&FdxFQ*(F^wFQ&F\\s6$Fdp,(F^sFQFbwFQ&Feu FfvFQFQFdwFQFQ>&Feq6$F[uF]r,&F]yFQ*(F^wFQ&F\\s6$FepFjxFQFdwFQFQ>&Feq6$ F[uFar,&FdyFQ*(F^wFQ&F\\s6$FfpFjxFQFdwFQFQC(>Fjv,&FjvFQ*(F^wFQ&F\\s6$F dp,(F^sFQFbwFQ&Feu6#FSFQFQFdwFQFQ>Ffw,&FfwFQ*(F^wFQ&F\\s6$FepF`zFQFdwF QFQ>F]x,&F]xFQ*(F^wFQ&F\\s6$FfpF`zFQFdwFQFQ>Fdx,&FdxFQ*(F^wFQ&F\\s6$Fd p,(F^sFQFbwFQ&Feu6#FgoFQFQFdwFQFU>F]y,&F]yFQ*(F^wFQ&F\\s6$FepFb[lFQFdw FQFU>Fdy,&FdyFQ*(F^wFQ&F\\s6$FfpFb[lFQFdwFQFUC)>Fft7#-%$seqG6$-%#opGFb u/86;FQFL@$/F`tFdp>FeuFS@$/F`tFep>FeuFiq@$/F`tFfp>FeuF]v>Fjv,&FjvFQ*(F ^wFQ&F\\s6$Fdp,(F^sFQFbwFQFeuFQFQFdwFQFQ>Ffw,&FfwFQ*(F^wFQ&F\\s6$FepFh ]lFQFdwFQFQ>F]x,&F]xFQ*(F^wFQ&F\\s6$FfpFh]lFQFdwFQFQ-Ffr6#%>integratin g~difference~fn...~G>82Fgn?(Fh\\lFQFQF_qFaqC(-%'printfG6$%$%d~GFh\\l@$ /-%$modG6$Fh\\l\"#:FQ-F[_l6#%\"|+G>85-%$addG6$*$),(&Feq6$&Fbp6#Fh\\l,& FWFQF^qFQFU*&-%)func_valG6'F\\sFeqFa`l&FYFhq9*FQ-Ff`l6'F\\sFeqFa`lFh`l 9)FQFQ-Ff`l6'F\\sFeqFa`lFh`l9+FUFgpFQ/F^q;FQ-Fio6#FY@%0FLFQ>84-Fj_l6$* $),(-Ff`l6'F\\sFeqFa`l&F[oFhqF]sFU*&-Ff`l6'F\\sFeqFa`lF_blFi`lFQ-Ff`l6 'F\\sFeqFa`lF_blF\\alFQFQ-Ff`l6'F\\sFeqFa`lF_blF_alFUFgpFQ/F^q;FQ-Fio6 #F[o>FgalFgn>83-Fj_l6$*$),(-Ff`l6'F\\sFeqFa`l&F\\pFhqF]sFU*&-Ff`l6'F\\ sFeqFa`lFeclFi`lFQ-Ff`l6'F\\sFeqFa`lFeclF\\alFQFQ-Ff`l6'F\\sFeqFa`lFec lF_alFUFgpFQ/F^q;FQ-Fio6#F\\p>Fg^l,&Fg^lFQ*&FNFQ,(F]clFgu*&FgpFQFgalFQ FQFh_lFQFQFQ-%%sqrtG6#Fg^lFEFEFE" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 607 "#This function takes the value of funcno at address \+ add when add is considered an address on the subtriangle\n#triangle. \+ For instance, if triangle is 0 (ie, the top 1/3 of the SG), then func_ val(T,S,triangle,[0,1],funcno)\n#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\n #that this requires that we have run get_triangle_adds, since coord_tr ansfer requires the addresses of the vertices\n#of the triangles to op erate...\n\nfunc_val:=proc(T,S,triangle,add,funcno)\nT[coord_transfer( add,triangle,S[triangle,10],S[triangle,11],S[triangle,12]),funcno] end ;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)func_valGf*6'%\"TG%\"SG%)tri angleG%$addG%'funcnoG6\"F,F,&9$6$-%/coord_transferG6'9'9&&9%6$F4\"#5&F 66$F4\"#6&F66$F4\"#79(F,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 702 "#So this function stores the addresses of the boundaries of the triangles\n#out to m_ba sis in S[triaddress,10-12], which 10,11,12 corresponding to add0,add1, add2 respectively.\n#It does V1 manually, and then calls get_triangle_ addresses, which recursively does the rest.\n\nget_triangle_adds:=proc (S,m_basis)\nif m_basis=1 then \n S[[0],10]:=[0];\n S[[0],11]:=[0, 1];\n S[[0],12]:=[0,2];\n S[[1],10]:=[0,1];\n S[[1],11]:=[1];\n \+ S[[1],12]:=[1,2];\n S[[2],10]:=[0,2];\n S[[2],11]:=[1,2];\n S[ [2],12]:=[2];\nelse\n get_triangle_addresses(S,[0],[0],[0,1],[0,2],m _basis-1);\n get_triangle_addresses(S,[1],[0,1],[1],[1,2],m_basis-1) ;\n get_triangle_addresses(S,[2],[0,2],[1,2],[2],m_basis-1);\nfi;\ne nd;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%2get_triangle_addsGf*6$% \"SG%(m_basisG6\"F)F)@%/9%\"\"\"C+>&9$6$7#\"\"!\"#5F3>&F16$F3\"#67$F4F ->&F16$F3\"#77$F4\"\"#>&F16$7#F-F5F:>&F16$FDF9FD>&F16$FDF>7$F-F@>&F16$ 7#F@F5F?>&F16$FOF9FK>&F16$FOF>FOC%-%7get_triangle_addressesG6(F1F3F3F: F?,&F,F-F-!\"\"-FX6(F1FDF:FDFKFZ-FX6(F1FOF?FKFOFZF)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1289 "#This is the recursive call which get_triangle_adds calls to generate the triangle_addresses for m_basi s>1...\n\nget_triangle_addresses:=proc(S,triangle,add0,add1,add2,m_bas is)\n\n#This is our base case. We simply assign the values according \+ to triangle\nif m_basis=1 then\n S[[op(triangle),0],10]:=add0;\n S [[op(triangle),0],11]:=[op(triangle),0,1];\n S[[op(triangle),0],12]: =[op(triangle),0,2];\n\n S[[op(triangle),1],10]:=[op(triangle),0,1]; \n S[[op(triangle),1],11]:=add1;\n S[[op(triangle),1],12]:=[op(tri angle),1,2];\n\n S[[op(triangle),2],10]:=[op(triangle),0,2];\n S[[ op(triangle),2],11]:=[op(triangle),1,2];\n S[[op(triangle),2],12]:=a dd2;\n\nelse\n #if we aren't in the base case, we simply call the pr ocedure again 3 times,\n #except that this time, we add on the appro priate extra digit onto triangle.\n #Notice that nothing gets assign ed until we reach the base case; the rest is \n #just setting it up \+ (primarily just figuring out which triangles to assign).\n get_trian gle_addresses(S,[op(triangle),0],add0,[op(triangle),0,1],[op(triangle) ,0,2],m_basis-1);\n get_triangle_addresses(S,[op(triangle),1],[op(tr iangle),0,1],add1,[op(triangle),1,2],m_basis-1);\n get_triangle_addr esses(S,[op(triangle),2],[op(triangle),0,2],[op(triangle),1,2],add2,m_ basis-1);\nfi;\nend;\n " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%7get_tri angle_addressesGf*6(%\"SG%)triangleG%%add0G%%add1G%%add2G%(m_basisG6\" F-F-@%/9)\"\"\"C+>&9$6$7$-%#opG6#9%\"\"!\"#59&>&F56$F7\"#67%F8F&F5 6$F7\"#77%F8F<\"\"#>&F56$7$F8F1F=FC>&F56$FMFB9'>&F56$FMFG7%F8F1FI>&F56 $7$F8FIF=FH>&F56$FYFBFU>&F56$FYFG9(C%-F$6(F5F7F>FCFH,&F0F1F1!\"\"-F$6( F5FMFCFQFUF^o-F$6(F5FYFHFUFjnF^oF-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 222 "#I don't know if this function is ever used, but it \+ returns the triangle\n#address of triangle of which add is a midpoint. \n\nget_tri_digits:=proc(add)\nif nops(add)=1 or nops(add) = 2 then [] else add[1..nops(add)-2] fi; end;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6# >%/get_tri_digitsGf*6#%$addG6\"F(F(@%5/-%%nopsG6#9$\"\"\"/F,\"\"#7\"&F /6#;F0,&F,F0F2!\"\"F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1290 "#This function creat es the basis functions for the harmonic\n#finite element method. The \+ basis functions are quite simple:\n#They are just 1 at one boundary, 0 at the other boundaries, \n#and filled in harmonically the rest of th e way. Here, I store\n#(1,0,0) in funcno, (0,1,0) in funcno+1, and (0 ,0,1) in funcno+2, \n#where the first number corresponds to the value \+ at [0], the second num\n#to the value at [1], and I think you can figu re out the last case.\n#It should be noted that funcno, funcno+1, and \+ funcno+2 are all used up\n#by this function.\n\ncalc_harms:=proc(T,m,f uncno)\n\nT[[0],funcno]:=1;\nT[[1],funcno]:=0;\nT[[2],funcno]:=0;\nT[[ 0,1],funcno]:=2/5;\nT[[1,2],funcno]:=1/5;\nT[[0,2],funcno]:=2/5;\nharm (T,m,[0],[0,1],[0,2],funcno);\nharm(T,m,[1],[1,2],[0,1],funcno);\nharm (T,m,[2],[0,2],[1,2],funcno);\n\n\nT[[0],funcno+1]:=0;\nT[[1],funcno+1 ]:=1;\nT[[2],funcno+1]:=0;\nT[[0,1],funcno+1]:=2/5;\nT[[1,2],funcno+1] :=2/5;\nT[[0,2],funcno+1]:=1/5;\nharm(T,m,[0],[0,1],[0,2],funcno+1);\n harm(T,m,[1],[1,2],[0,1],funcno+1);\nharm(T,m,[2],[0,2],[1,2],funcno+1 );\n\n\nT[[0],funcno+2]:=0;\nT[[1],funcno+2]:=0;\nT[[2],funcno+2]:=1; \nT[[0,1],funcno+2]:=1/5;\nT[[1,2],funcno+2]:=2/5;\nT[[0,2],funcno+2]: =2/5;\nharm(T,m,[0],[0,1],[0,2],funcno+2);\nharm(T,m,[1],[1,2],[0,1],f uncno+2);\nharm(T,m,[2],[0,2],[1,2],funcno+2);\n\nend;\n" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%+calc_harmsGf*6%%\"TG%\"mG%'funcnoG6\"F*F*C=>& 9$6$7#\"\"!9&\"\"\">&F.6$7#F3F2F1>&F.6$7#\"\"#F2F1>&F.6$7$F1F3F2#F<\" \"&>&F.6$7$F3F&F.6$7$F1F&F.6$F0,&F2F3F3F3F1>&F.6$F7FWF3>&F.6$F;FWF1>& F.6$F@FWFA>&F.6$FFFWFA>&F.6$FKFWFG-FM6(F.FOF0F@FKFW-FM6(F.FOF7FFF@FW-F M6(F.FOF;FKFFFW>&F.6$F0,&F2F3F&F.6$F7FjoF1>&F.6$F;FjoF3>&F.6$F@F joFG>&F.6$FFFjoFA>&F.6$FKFjoFA-FM6(F.FOF0F@FKFjo-FM6(F.FOF7FFF@Fjo-FM6 (F.FOF;FKFFFjoF*F*F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 "\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1053 "#This function computes the inner product of t wo harmonic basis elements.\n#It uses theoretical values for inner pro ducts of harmonic functions. It\n#then scales down the answer by a fa ctor of 3^m, because the triangles are\n#smaller by a factor of 3^m.\n \ninner_prod_h:=proc(S,m,basis_1,basis_2)\n\n local out,temp,basis1, basis2;\n\n basis1:=basis_1;\n basis2:=basis_2;\n\n#if they aren't neighbors, we get nothing...\n if basis1[1] <> basis2[1] and S[basi s1[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;\n\n \n#If the 2 basis elts are the same, then we get 7/45*2 (the two comes from the 2 neighboring triangles)\n#the value 7/45 is theoretical\n \+ if basis1[1]=basis2[1] then\n out:=7/45 * 2 / (3^m);\n fi;\n\n #If the 2 basis elts are neighbors, then we get 4/45 (theoretical valu e)\n if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or \+ S[basis1[1],3] = basis2[1] or\nS[basis1[1],4] = basis2[1]\nthen\n \+ out:=4/45 / (3^m); \n fi;\n \n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-inner_prod_hGf*6&%\"SG%\"mG%(basis_1G%(ba sis_2G6&%$outG%%tempG%'basis1G%'basis2G6\"F0C(>8&9&>8'9'@$33330&F36#\" \"\"&F6F?0&9$6$F>F@FA0&FD6$F>\"\"#FA0&FD6$F>\"\"$FA0&FD6$F>\"\"%FA-%'R ETURNG6#\"\"!@$/F>FA>8$,$*&F@F@)FM9%!\"\"#\"#9\"#X@$555/FCFA/FGFA/FKFA /FOFA>FY,$Fen#FQF[oFYF0F0F0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 778 "#This function returns the energy of two basis functions.\n#It wo rks just like inner_prod_h, except that it uses a scaling factor of\n# (5/3)^m, which is the energy scaling factor from kigami.\n\n\nenergy_h :=proc(S,m,basis_1,basis_2)\n\n local out,temp,basis1,basis2;\n\n \+ basis1:=basis_1;\n basis2:=basis_2;\n\n#if they aren't neighbors, we get nothing...\n if basis1[1] <> basis2[1] and S[basis1[1],1] <> ba sis2[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;\n\n if basis1[1 ]=basis2[1] then\n out:=2 * 2 * (5/3)^m;\n fi;\n if S[basis1[ 1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[basis1[1],3] = ba sis2[1] or\nS[basis1[1],4] = basis2[1]\nthen\n out:=-1 * (5/3)^m; \n fi;\n \n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)energy_hGf*6&%\"SG%\"mG%(basis_1G%(basis_2G6&%$outG%%tempG%'b asis1G%'basis2G6\"F0C(>8&9&>8'9'@$33330&F36#\"\"\"&F6F?0&9$6$F>F@FA0&F D6$F>\"\"#FA0&FD6$F>\"\"$FA0&FD6$F>\"\"%FA-%'RETURNG6#\"\"!@$/F>FA>8$, $)#\"\"&FM9%FQ@$555/FCFA/FGFA/FKFA/FOFA>FY,$Fen!\"\"FYF0F0F0" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5041 "#This function takes a bas is element and a function, and returns integral(basis*function).\n#The nice thing is that the basis element is zero everywhere except for 2 \+ triangles. Thus,\n#we can speed computations by just integrating over those particular triangles.\n#Note that the integration method used h ere is simpson's rule out to maximum depth in the SG.\n\nint_vs_basis_ h:=proc(T,S,m_basis,m_func,basis,funct,harm_funcno)\n\nlocal V0,Vother ,Vm,fact,triangle1,triangle2,V0sum1,V0sum2,Vothersum1,Vothersum2,Vmsum 1,Vmsum2,i,j,Vother_e;\n\nfact:= 1/(6*3^m_func); #This is the factor \+ for simpson's rule...\n\n\n#get the vertices for V0, Vm, and everythin g in between...\n#Notice that here we are taking m to be m_func-m_basi s, because what we want to do is\n#consider vertices scaled down by m_ basis in funct and multiply by vertices in the basis\n\nV0:=vertices_m (0);\nif m_func-m_basis = 1 then Vother:=[] else Vother_e:=vertices(m_ func-m_basis-1);Vother:=Vother_e[4..nops(Vother_e)]; fi;\nVm:=vertices _m(m_func-m_basis);\n\n#all basis elements are of type f...\nif basis[ 2]=f then\n\n #So basically, we will first check the orientation of the basis element (i.e, the last two digits).\n #Then, we get the \+ triangles associated with the basis element. Then, given a point in V m_func-m_basis,\n #we can use coord_transfer to determine (using th e triangle) which is the corresponding coordinate\n #in Vm_func, an d then multiply that by the corresponding value of the basis function, and sum.\n \n if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1 ] then \n triangle1:=get_triangle(S,basis,m_basis,1);\n tria ngle2:=get_triangle(S,basis,m_basis,2);\n\n#let's calculate sum's for \+ triangle 1, shall we?\n V0sum1:=add(T[coord_transfer(V0[i],triang le1,basis[1],S[basis[1],1],S[basis[1],2]),funct] *(T[V0[i],harm_funcno ]),i=1..nops(V0));\n\n 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)); \n \n 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));\n\n#and now on to triangle num ero 2!\n 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));\n\n 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)); \n\n 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));\n\n fi;\n\n if basis[1][nops(bas is[1])-1..nops(basis[1])] = [0,2] then\n triangle1:=get_triangle( S,basis,m_basis,1);\n triangle2:=get_triangle(S,basis,m_basis,2); \n\n#let's calculate sum's for triangle 1, shall we?\n V0sum1:=ad d(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));\n\n Vothersum 1:=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)); \+ \n \n 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));\n\n#and now on to triangle numero 2!\n 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));\n\n Vothersum2:=add(T[co ord_transfer(Vother[i],triangle2,basis[1],S[basis[1],3],S[basis[1],4]) ,funct] *(T[Vother[i],harm_funcno]),i=1..nops(Vother)); \n\n \+ 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));\n\n \+ fi; \n\n if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then \n triangle1:=get_triangle(S,basis,m_basis,1);\n triangle2:= get_triangle(S,basis,m_basis,2);\n\n#let's calculate sum's for triangl e 1, shall we?\n V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[b asis[1],2],basis[1],S[basis[1],1]),funct]*(T[V0[i],harm_funcno+1]),i=1 ..nops(V0));\n\n Vothersum1:=add(T[coord_transfer(Vother[i],trian gle1,S[basis[1],2],basis[1],S[basis[1],1]),funct] *(T[Vother[i],harm_f uncno+1]),i=1..nops(Vother)); \n \n Vmsum1:=add(T[coord_tran sfer(Vm[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*(T[V m[i],harm_funcno+1]),i=1..nops(Vm));\n\n#and now on to triangle numero 2!\n 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)) ;\n\n Vothersum2:=add(T[coord_transfer(Vother[i],triangle2,S[basi s[1],3],S[basis[1],4],basis[1]),funct] *(T[Vother[i],harm_funcno+2]),i =1..nops(Vother)); \n\n 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));\n\n fi; \nfi;\n\n\n#And then just return the formula for simpson's rule...\nRETURN(fact * (5*(Vmsum1+Vmsum2) + 2*(Vothersum1+Vothersum2) + (V0sum1+V0sum2)));\n\n\nend;\n \n \n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/int_vs_basis_hGf*6)%\"TG%\" SG%(m_basisG%'m_funcG%&basisG%&functG%,harm_funcnoG61%#V0G%'VotherG%#V mG%%factG%*triangle1G%*triangle2G%'V0sum1G%'V0sum2G%+Vothersum1G%+Voth ersum2G%'Vmsum1G%'Vmsum2G%\"iG%\"jG%)Vother_eG6\"F>C(>8',$*&\"\"\"FD) \"\"$9'!\"\"#FD\"\"'>8$-%+vertices_mG6#\"\"!@%/,&FGFD9&FHFD>8%7\"C$>82 -%)verticesG6#,(FGFDFTFHFDFH>FV&FZ6#;\"\"%-%%nopsG6#FZ>8&-FN6#FS@$/&9( 6#\"\"#%\"fGC%@$/&&Fho6#FD6#;,&-F_o6#F`pFDFDFHFep7$FPFDC*>8(-%-get_tri angleG6&9%FhoFTFD>8)-F\\q6&F^qFhoFTFjo>8*-%$addG6$*&&9$6$-%/coord_tran sferG6'&FL6#80FjpF`p&F^q6$F`pFD&F^q6$F`pFjo9)FD&Fjq6$F_r9*FD/Far;FD-F_ o6#FL>8,-Ffq6$*&&Fjq6$-F]r6'&FVF`rFjpF`pFbrFdrFfrFD&Fjq6$FgsFirFD/Far; FD-F_o6#FV>8.-Ffq6$*&&Fjq6$-F]r6'&FboF`rFjpF`pFbrFdrFfrFD&Fjq6$FgtFirF D/Far;FD-F_o6#Fbo>8+-Ffq6$*&&Fjq6$-F]r6'F_rF`q&F^q6$F`pF]oF`p&F^q6$F`p FFFfrFD&Fjq6$F_r,&FirFDFDFDFDFjr>8--Ffq6$*&&Fjq6$-F]r6'FgsF`qFguF`pFiu FfrFD&Fjq6$FgsF]vFDFjs>8/-Ffq6$*&&Fjq6$-F]r6'FgtF`qFguF`pFiuFfrFD&Fjq6 $FgtF]vFDFjt@$/F_p7$FPFjoC*>FjpF[q>F`qFaq>Fdq-Ffq6$*&&Fjq6$-F]r6'F_rFj pFbrFdrF`pFfrFD&Fjq6$F_r,&FirFDFjoFDFDFjr>F_s-Ffq6$*&&Fjq6$-F]r6'FgsFj pFbrFdrF`pFfrFD&Fjq6$FgsFdxFDFjs>F_t-Ffq6$*&&Fjq6$-F]r6'FgtFjpFbrFdrF` pFfrFD&Fjq6$FgtFdxFDFjt>F_u-Ffq6$*&&Fjq6$-F]r6'F_rF`qF`pFiuFguFfrFDFgr FDFjr>F_v-Ffq6$*&&Fjq6$-F]r6'FgsF`qF`pFiuFguFfrFDFhsFDFjs>Fjv-Ffq6$*&& Fjq6$-F]r6'FgtF`qF`pFiuFguFfrFDFhtFDFjt@$/F_p7$FDFjoC*>FjpF[q>F`qFaq>F dq-Ffq6$*&&Fjq6$-F]r6'F_rFjpFdrF`pFbrFfrFDF[vFDFjr>F_s-Ffq6$*&&Fjq6$-F ]r6'FgsFjpFdrF`pFbrFfrFDFgvFDFjs>F_t-Ffq6$*&&Fjq6$-F]r6'FgtFjpFdrF`pFb rFfrFDFbwFDFjt>F_u-Ffq6$*&&Fjq6$-F]r6'F_rF`qFiuFguF`pFfrFDFbxFDFjr>F_v -Ffq6$*&&Fjq6$-F]r6'FgsF`qFiuFguF`pFfrFDF]yFDFjs>Fjv-Ffq6$*&&Fjq6$-F]r 6'FgtF`qFiuFguF`pFfrFDFgyFDFjt-%'RETURNG6#*&FAFD,.F_t\"\"&*&F\\_lFDFjv FDFD*&FjoFDF_sFDFD*&FjoFDF_vFDFDFdqFDF_uFDFDF>F>F>" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 546 "#Here, we generate the matrix of E(basisi, basisj) out to level m_basis\n#The final matrix is a sparse, symmetric matrix, and we take that into\n#account here.\n\ngen_energy_matrix_h: =proc(S,m_basis,vm)\n\nlocal i,j,en_matrix,n_verts,en;\n\nn_verts:=nop s(vm)-3; sum(3^i,i=1..m_basis);\n\nen_matrix:=array(sparse,1..n_verts, 1..n_verts);\n\nfor i from 1 to n_verts do\n for j from 1 to n_verts do\n en:=energy_h(S,m_basis,get_basis(vm,i),get_basis(vm,j));\n \+ if en <> 0 then en_matrix[i,j]:=en; en_matrix[j,i]:=en; fi;\n o d;\nod;\n\nen_matrix;\nend; \n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 483 "#This function returns a basis element given any particular basis element number.\n#Notice that this function works for both harmonic a nd biharmonic splines; if the \n#number is greater than the number of \+ f-type vectors, we return a g-type vector,\n#but if we are using harmo nic splines, we will never get a number that high, and\n#so all we get are the f-type vectors.\n\nget_basis:=proc(vm,number)\n\nif number <= nops(vm)-3 then [vm[number+3],f]\nelse [vm[number-nops(vm)+3],g] fi; \+ end; \n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%4gen_energy_matrix_hGf* 6%%\"SG%(m_basisG%#vmG6'%\"iG%\"jG%*en_matrixG%(n_vertsG%#enG6\"F0C'>8 ',&-%%nopsG6#9&\"\"\"\"\"$!\"\"-%$sumG6$)F:8$/F@;F99%>8&-%&arrayG6%%'s parseG;F9F3FJ?(F@F9F9F3%%trueG?(8%F9F9F3FLC$>8(-%)energy_hG6&9$FC-%*ge t_basisG6$F8F@-FW6$F8FN@$0FQ\"\"!C$>&FE6$F@FNFQ>&FE6$FNF@FQFEF0F0F0" } }{PARA 12 "" 1 "" {XPPMATH 20 "6#>%*get_basisGf*6$%#vmG%'numberG6\"F)F )@%19%,&-%%nopsG6#9$\"\"\"\"\"$!\"\"7$&F16#,&F,F2F3F2%\"fG7$&F16#,(F,F 2F.F4F3F2%\"gGF)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 473 "#T his is just like gen_energy_matrix_h, except that creates the inner pr oduct matrix.\n\ngen_iprod_matrix_h:=proc(S,m_basis,vm)\n\nlocal i,j,i _matrix,n_verts,ip;\n\nn_verts:=nops(vm)-3; add(3^i,i=1..m_basis);\n\n i_matrix:=array(sparse,1..n_verts,1..n_verts);\n\nfor i from 1 to n_ve rts do\n for j from 1 to n_verts do\n ip:=inner_prod_h(S,m_basi s,get_basis(vm,i),get_basis(vm,j));\n if ip <> 0 then i_matrix[i, j]:=ip; i_matrix[j,i]:=ip; fi;\n od;\nod;\n\ni_matrix;\nend; \n " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3gen_iprod_matrix_hGf*6%%\"SG%(m _basisG%#vmG6'%\"iG%\"jG%)i_matrixG%(n_vertsG%#ipG6\"F0C'>8',&-%%nopsG 6#9&\"\"\"\"\"$!\"\"-%$addG6$)F:8$/F@;F99%>8&-%&arrayG6%%'sparseG;F9F3 FJ?(F@F9F9F3%%trueG?(8%F9F9F3FLC$>8(-%-inner_prod_hG6&9$FC-%*get_basis G6$F8F@-FW6$F8FN@$0FQ\"\"!C$>&FE6$F@FNFQ>&FE6$FNF@FQFEF0F0F0" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 414 "#This procedure generates a vector of integral(funcno*basisi).\n#It calls int_vs_basis_h, where m ost of the hard work is actually\n#done.\n\n\ngen_f_vector_h:=proc(T,S ,m_basis,m_func,vm,funcno,harm_funcno)\n\nlocal i,f_vec,n_verts;\n\n\n n_verts:=add(3^i,i=1..m_basis);\n\nf_vec:=array(1..n_verts);\n\nfor i \+ from 1 to n_verts do\n f_vec[i]:=int_vs_basis_h(T,S,m_basis,m_func,g et_basis(vm,i),funcno,harm_funcno);\nod;\n\nf_vec;\nend;" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%/gen_f_vector_hGf*6)%\"TG%\"SG%(m_basisG%'m_fu ncG%#vmG%'funcnoG%,harm_funcnoG6%%\"iG%&f_vecG%(n_vertsG6\"F2C&>8&-%$a ddG6$)\"\"$8$/F;;\"\"\"9&>8%-%&arrayG6#;F>F5?(F;F>F>F5%%trueG>&FA6#F;- %/int_vs_basis_hG6)9$9%F?9'-%*get_basisG6$9(F;9)9*FAF2F2F2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1813 "#This function is called by make_ pretty_picture_h. It basically goes through not assigning values to \+ \n#any vertices until it reaches the end of m_basis, and then it just \+ fills everything in harmonically,\n#basically making a piecewise harmo nic function with boundary values at the harmonic spline basis\n#verti ces.\n\nmake_h_answer:=proc(T,m_basis,m_func,add0,add1,add2,funcno);\n \nif m_basis = 0 then\n# print(add0,add1,add2);\n harm(T,m_func, add0,add1,add2,funcno);\nelse\n if add1[nops(add1)-1..nops(add1)] = [0,1] then\n make_h_answer(T,m_basis-1,m_func-1,add0,[op(a dd1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno) ;\n 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);\n \+ 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);\n fi;\n if a dd1[nops(add1)-1..nops(add1)] = [1,2] then\n 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);\n 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);\n 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],fu ncno);\n fi;\n if add1[nops(add1)-1..nops(add1)] = [0,2] the n\n 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);\n \+ 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);\n make_h_answ er(T,m_basis-1,m_func-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add 1[1..nops(add1)-2]),2,0,1],funcno);\n fi;\n\nfi;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.make_h_answerGf*6)%\"TG%(m_basisG%'m_func G%%add0G%%add1G%%add2G%'funcnoG6\"F.F.@%/9%\"\"!-%%harmG6(9$9&9'9(9)9* C%@$/&F96#;,&-%%nopsG6#F9\"\"\"FF!\"\"FC7$F2FFC%-F$6)F6,&F1FFFFFG,&F7F FFFFGF87&-%#opG6#&F96#;FF,&FCFF\"\"#FGF2F2FF7&FOF2F2FVF;-F$6)F6FLFMF97 &FOF2FFFVFNF;-F$6)F6FLFMF:FWFZF;@$/F?7$FFFVC%-F$6)F6FLFMF87&FOFFFFFV7& FOFFF2FFF;-F$6)F6FLFMF97&FOFFF2FVF]oF;-F$6)F6FLFMF:F^oFaoF;@$/F?7$F2FV C%-F$6)F6FLFMF87&FOFVF2FV7&FOFVFFFVF;-F$6)F6FLFMF97&FOFVF2FFFjoF;-F$6) F6FLFMF:F[pF^pF;F.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1043 "#This function basically copies the values of C (a projective represe ntation of a function\n#in the harmonic spline basis) into a function, and then calls make_h_answer to fill in the\n#rest of the function. \+ This procedure is far more efficient than the method used in\n#make_pr etty_picture because there, we have to add in each basis function manu ally.\n#Here, we can use the fact that the harmonic functions are comp letely defined\n#by knowing the values of the function on the boundary and then just using the harmonic\n#algorithm to determine the rest.\n \nmake_pretty_picture_h:=proc(T,m_basis,m_func,C,funcno)\n\nlocal vm,j ;\nvm:=vertices(m_basis); \n\n#first, set the initial values...\nT[[0] ,funcno]:=0 ;T[[1],funcno]:=0; T[[2],funcno]:=0; \n\nfor j from 1 to n ops([indices(C)]) do\nT[vm[j+3],funcno]:=C[j];\nod;\n\n#then call make _h_answer to do all the real work...\nmake_h_answer(T,m_basis-1,m_func -1,[0],[0,1],[0,2],funcno);\nmake_h_answer(T,m_basis-1,m_func-1,[1],[1 ,2],[0,1],funcno);\nmake_h_answer(T,m_basis-1,m_func-1,[2],[0,2],[1,2] ,funcno);\n\nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%6make_prett y_picture_hGf*6'%\"TG%(m_basisG%'m_funcG%\"CG%'funcnoG6$%#vmG%\"jG6\"F /C*>8$-%)verticesG6#9%>&9$6$7#\"\"!9(F<>&F96$7#\"\"\"F=F<>&F96$7#\"\"# F=F&F96$&F26#,&FIFB\"\"$FB F=&FQ6#FI-%.make_h_answerG6)F9,&F6FBFB!\"\",&9&FBFBFjnF;7$F " 0 "" {MPLTEXT 1 0 1158 "#This function returns a list of \+ which harmonic basis functions\n#should be used for integrating versus two basis elements.\n#It assumes that the basis elements are neighbor s, because otherwise\n#it's a fairly trivial question. It is called b y int_vs_basis12_h.\n#Notice that the values depend on orientation as \+ well as neighbors...\n\nwhich_funcs_h:=proc(S,basis1,basis2)\nlocal ou t;\n\n\nif last_digits(basis1[1]) = [0,1] then\n if S[basis1[1],1] = basis2[1] then out:=[0,1] fi;\n if S[basis1[1],2] = basis2[1] then \+ out:=[0,2] fi;\n if S[basis1[1],3] = basis2[1] then out:=[1,2] fi;\n if S[basis1[1],4] = basis2[1] then out:=[1,0] fi;\nfi;\nif last_dig its(basis1[1]) = [0,2] then\n if S[basis1[1],1] = basis2[1] then out :=[2,0] fi;\n if S[basis1[1],2] = basis2[1] then out:=[2,1] fi;\n \+ if S[basis1[1],3] = basis2[1] then out:=[0,1] fi;\n if S[basis1[1],4 ] = basis2[1] then out:=[0,2] fi;\nfi;\nif last_digits(basis1[1]) = [1 ,2] then\n if S[basis1[1],1] = basis2[1] then out:=[1,2] fi;\n if \+ S[basis1[1],2] = basis2[1] then out:=[1,0] fi;\n if S[basis1[1],3] = basis2[1] then out:=[2,0] fi;\n if S[basis1[1],4] = basis2[1] then \+ out:=[2,1] fi;\nfi;\n\nout;\nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.which_funcs_hGf*6%%\"SG%'basis1G%'basis2G6#%$outG6\"F,C&@$/-%,la st_digitsG6#&9%6#\"\"\"7$\"\"!F6C&@$/&9$6$F3F6&9&F5>8$F7@$/&F=6$F3\"\" #F?>FB7$F8FG@$/&F=6$F3\"\"$F?>FB7$F6FG@$/&F=6$F3\"\"%F?>FB7$F6F8@$/F0F IC&@$F;>FB7$FGF8@$FD>FB7$FGF6@$FK>FBF7@$FR>FBFI@$/F0FPC&@$F;>FBFP@$FD> FBFW@$FK>FBFgn@$FR>FBFjnFBF,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4108 "#This function is used to generate the matrix integ ral(q,basisi,basisj).\n#It basically has to go through the logic of in ner_prod_h, except that instead\n#of just returning a theoretical valu e, it returns the appropriate integral.\n\nint_vs_basis12_h:=proc(T,S, m_basis,m_func,basis_1,basis_2,funct,harm_funcno)\n\n local out,temp ,basis1,basis2,fact,triangle1,triangle2,addresses1,addresses2,j,which_ harms;\n\n fact:=1/(3^m_basis); #scaling factor due to m_basis\n \+ basis1:=basis_1;\n basis2:=basis_2;\n# print(basis1,basis2);\n\n#i f they aren't neighbors, we get nothing...\n 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 RETU RN(0) fi;\n\n if basis1[1]=basis2[1] then\n triangle1:=get_tria ngle(S,basis1,m_basis,1); #first, let's get the triangles, as well a s the addresses of the vertices of those triangles\n triangle2:=g et_triangle(S,basis1,m_basis,2);\n addresses1:=get_addresses(S,ba sis1,m_basis,1); #this returns the addresses of the coordinates of th e triangles\n addresses2:=get_addresses(S,basis1,m_basis,2);\n \+ \n #Notice that if basis1 and basis2 are in the same place, then we have to integrate over two triangles\n #Another impor tant thing is that the harmonic basis elements that you integrate agai nst dependent solely\n #upon the orientation of the basis elem ents, since they are in the same place.\n #Notice that a lot o f the work here is done in integrate_basis12, which is in splines.mws. \n if last_digits(basis1[1]) = [0,1] then\n out:=integ rate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],ad dresses1[3],funct,harm_funcno+0,harm_funcno+0)+\n integra te_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addr esses2[3],funct,harm_funcno+1,harm_funcno+1); fi;\n\n if last_ digits(basis1[1]) = [0,2] then\n out:=integrate_basis12(T,m_fu nc-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,h arm_funcno+2,harm_funcno+2)+\n integrate_basis12(T,m_func -m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,har m_funcno+0,harm_funcno+0); fi;\n\n if last_digits(basis1[1]) = [1,2] then\n out:=integrate_basis12(T,m_func-m_basis,triangle 1,addresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+1,harm_f uncno+1)+\n integrate_basis12(T,m_func-m_basis,triangle2, addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+2,harm_fun cno+2); fi;\n \n fi; #END OF CASE WHERE BASIS1 and BASIS2 are i n same place\n\n #And so now what if the basis elements are neighbor s \n 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\n \+ if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then\n #This case is if the overlapping triangle is the first t riangle of basis1 (in the clockwise sense).\n triangle1:= get_triangle(S,basis1,m_basis,1); #first, let's get the triangle... \n addresses1:=get_addresses(S,basis1,m_basis,1); #and th e corresponding address...\n\n which_harms:=which_funcs_h (S,basis1,basis2); #This function does most of the work for determin ing which basis functions to use.\n\n out:=integrate_basi s12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[ 3],funct, harm_funcno+which_harms[1],harm_funcno+which_harms[2]);\n \+ else\n\n #This case is if the overlapping triangl e is the second triangle of basis1 (in the clockwise sense).\n \+ #Otherwise, it's just like before.\n triangle1:=get _triangle(S,basis1,m_basis,2);\n addresses1:=get_addresse s(S,basis1,m_basis,2);\n which_harms:=which_funcs_h(S,bas is1,basis2);\n out:=integrate_basis12(T,m_func-m_basis,tr iangle1,addresses1[1],addresses1[2],addresses1[3],funct, harm_funcno+w hich_harms[1],harm_funcno+which_harms[2]);\n fi;\n\n fi;\n \n out*fact; #Don't forget the scaling factor! \nend;\n\n\n" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%1int_vs_basis12_hGf*6*%\"TG%\"SG%(m_ basisG%'m_funcG%(basis_1G%(basis_2G%&functG%,harm_funcnoG6-%$outG%%tem pG%'basis1G%'basis2G%%factG%*triangle1G%*triangle2G%+addresses1G%+addr esses2G%\"jG%,which_harmsG6\"F;C)>8(*&\"\"\"F@)\"\"$9&!\"\">8&9(>8'9)@ $33330&FF6#F@&FIFR0&9%6$FQF@FS0&FV6$FQ\"\"#FS0&FV6$FQFBFS0&FV6$FQ\"\"% FS-%'RETURNG6#\"\"!@$/FQFSC)>8)-%-get_triangleG6&FVFFFCF@>8*-Fgo6&FVFF FCFen>8+-%.get_addressesGFho>8,-F`pF\\p@$/-%,last_digitsG6#FQ7$F`oF@>8 $,&-%2integrate_basis12G6+9$,&9'F@FCFDFeo&F^pFR&F^p6#Fen&F^p6#FB9*9+Fi qF@-F^q6+F`qFaqFjo&FbpFR&FbpFeq&FbpFgqFhq,&FiqF@F@F@F_rF@@$/Ffp7$F`oFe n>F[q,&-F^q6+F`qFaqFeoFcqFdqFfqFhq,&FiqF@FenF@FgrF@-F^q6+F`qFaqFjoF\\r F]rF^rFhqFiqFiqF@@$/Ffp7$F@Fen>F[q,&-F^q6+F`qFaqFeoFcqFdqFfqFhqF_rF_rF @-F^q6+F`qFaqFjoF\\rF]rF^rFhqFgrFgrF@@$555/FUFS/FYFS/FgnFS/FjnFS@%FfsC &>FeoFfo>F^pF_p>8.-%.which_funcs_hG6%FVFFFI>F[q-F^q6+F`qFaqFeoFcqFdqFf qFhq,&FiqF@&F`tFRF@,&FiqF@&F`tFeqF@C&>FeoF[p>F^pFcp>F`tFat>F[qFet*&F[q F@F>F@F;F;F;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 761 "#This func tion generates the matrix integral(q*basisi,basisj) which is used for \+ q-nonconstant\n#in the finite element method. The matrix is sparse an d symmetrical, and both of these properties are\n#used here. Most of \+ the work is done in int_vs_basis12_h.\n\n\ngen_q_matrix_h:=proc(T,S,m_ basis,m_func,vm,funcno,harm_funcno)\n\nlocal i,j,q_matrix,q_val,n_vert s;\n\nn_verts:=nops(vm)-3;\n\nq_matrix:=array(sparse,1..n_verts,1..n_v erts);\n\n\nfor i from 1 to n_verts do\n \n printf(`%d `,i);\n i f ((i mod 15)=1) then\n printf(`\\n`): fi:\n \n for j from 1 to \+ i do\n q_val:=int_vs_basis12_h(T,S,m_basis,m_func,get_basis(vm,i) ,get_basis(vm,j),funcno,harm_funcno);\n if q_val <> 0 then q_matr ix[i,j]:=q_val; q_matrix[j,i]:=q_val;fi;\n od;\nod;\n\nq_matrix;\nen d; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/gen_q_matrix_hGf*6)% \"TG%\"SG%(m_basisG%'m_funcG%#vmG%'funcnoG%,harm_funcnoG6'%\"iG%\"jG%) q_matrixG%&q_valG%(n_vertsG6\"F4C&>8(,&-%%nopsG6#9(\"\"\"\"\"$!\"\">8& -%&arrayG6%%'sparseG;F=F7FF?(8$F=F=F7%%trueGC%-%'printfG6$%$%d~GFH@$/- %$modG6$FH\"#:F=-FL6#%\"|+G?(8%F=F=FHFIC$>8'-%1int_vs_basis12_hG6*9$9% 9&9'-%*get_basisG6$F&FA6$FHFYFfn>&FA6$F YFHFfnFAF4F4F4" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5025 "#This function returns the inner product of two basis vectors in the biharmonic basis.\n#It retu rns theoretical values from the splines on fractals paper (scaled down by the appropriate\n#factor of 3 because we are dealing with smaller \+ triangles...). Notice there are several\n#cases to deal with, and one must also be careful to use the fact that we have defined the \n#norm al derivative to be negative in the 2nd triangle...\n\ninner_prod:=pro c(S,m,basis_1,basis_2)\n\n local out,temp,basis1,basis2;\n\n basis 1:=basis_1;#We do this so that we can change basis1 and basis2\n bas is2:=basis_2;\n\n#if they aren't neighbors, we get nothing...\n if b asis1[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;\n\n#first deal with special case of eithe r basis1 or basis2 elt of v0:\n\n 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\n #first make basis2 to be definitely in V0 so that we can keep track of it.\n if nops(basis1[1])=1 then\n temp :=basis2;\n basis2:=basis1;\n basis1:=temp;\n fi; \n \n if basis1[1]=basis2[1] then #I have included other cases \+ here for reasons I'm not even sure of.\n \+ #but the only case which matters is g and g, because only g type \+ basis\n #elements can exist on the b oundaries...\n if (basis1[2] = f and basis2[2] = g) or (basis1 [2] = g and basis2[2] = f) then out:=-47/1395 / (3^m) fi;\n if basis1[2] = f and basis2[2] = f then out:=190/837 / (3^m) fi;\n \+ if basis1[2] = g and basis2[2] = g then out:=206/37665 / (3^m) fi; \n else\n 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\n #we will always make it inside of this if statement \+ because the basis elts are always neighbors...\n if basis1[ 2] = f and basis2[2] = f then out:=89/1674 / (3^m) fi; #This case nev er happens...\n if basis1[2] = g and basis2[2] = g then #If we get within this case, then we have to figure out the signs\n \+ if S[basis1[1],1] = basis2[1] then out:=83/37665 / (3^m) fi; #what helps is the fact that the sign of \n if S[basis 1[1],4] = basis2[1] then out:=-83/37665 / (3^m) fi; #the boundary is a lways positive...\n fi;\n if (basis1[2] = f and \+ basis2[2] = g) then out:=-61/5580 / (3^m) fi;\n if (basis1[ 2] = g and basis2[2] = f) then #This case never happens, because basi s2 can never be f\n if S[basis1[1],1] = basis2[1] then \+ out:=-61/5580 / (3^m) fi; \n if S[basis1[1],4] = basis2 [1] then out:=61/5580 / (3^m) fi;\n fi;\n\n fi;\n \+ fi; \n\n else\n#OTHERWISE, if neither are on the boundari es...\n\n if basis1[1]=basis2[1] then #This is the case of the basi s elts being in the same place...\n if (basis1[2] = f and basis2[ 2] = g) or (basis1[2] = g and basis2[2] = f) then out:=0 fi; #by odd s ymmetry of g, this integral is 0\n if basis1[2] = f and basis2[2] = f then out:=190/837 * 2 / (3^m) fi;\n if basis1[2] = g and bas is2[2] = g then out:=206/37665 * 2 / (3^m) fi;\n fi;\n if S[basis1 [1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[basis1[1],3] = b asis2[1] or S[basis1[1],4] = basis2[1] then\n #again, we will always enter this if statement because they are always neighbors...\n i f basis1[2] = f and basis2[2] = f then out:=89/1674 / (3^m) fi;\n \+ if basis1[2] = g and basis2[2] = g then #Here, we have to be careful \+ of signs, using the fact that\n \+ #the 1 and 2 neighbors are in the 1 triangle, which is positive ... \n if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basi s2[1] then\n if S[basis2[1],1] = basis1[1] or S[basis2[1], 2] = basis1[1] then out:=83/37665 / (3^m) fi; \n if S[basi s2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] then out:=-83/37665 / (3^m) fi;\n fi; \n if S[basis1[1],3] = basis2[1] \+ or S[basis1[1],4] = basis2[1] then # out:=-83/37665 / (3^m) fi;\n \+ if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] t hen out:=-83/37665 / (3^m) fi; \n if S[basis2[1],3] = basi s1[1] or S[basis2[1],4] = basis1[1] then out:=83/37665 / (3^m) fi;\n \+ fi; \n\n fi;\n if (basis1[2] = f and basis2[2] = g) \+ then \n if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = bas is1[1] then out:=-61/5580 / (3^m) fi; \n if S[basis2[1],3] = \+ basis1[1] or S[basis2[1],4] = basis1[1] then out:=61/5580 / (3^m) fi; \n fi;\n if (basis1[2] = g and basis2[2] = f) then\n \+ if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then out :=-61/5580 / (3^m) fi;\n if S[basis1[1],3] = basis2[1] or S[b asis1[1],4] = basis2[1] then out:=61/5580 / (3^m) fi;\n fi; \n\n \+ \n fi;\n fi; \n out; \nend;\n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+inner_prodGf*6& %\"SG%\"mG%(basis_1G%(basis_2G6&%$outG%%tempG%'basis1G%'basis2G6\"F0C' >8&9&>8'9'@$33330&F36#\"\"\"&F6F?0&9$6$F>F@FA0&FD6$F>\"\"#FA0&FD6$F>\" \"$FA0&FD6$F>\"\"%FA-%'RETURNG6#\"\"!@%55555/F>7#FU/F>7#F@/F>7#FI/FAFg n/FAFin/FAF[oC$@$/-%%nopsG6#F>F@C%>8%F6>F6F3>F3Fgo@%/F>FAC%@$53/&F36#F I%\"fG/&F6Fbp%\"gG3/FapFfp/FepFcp>8$,$*&F@F@)FM9%!\"\"#!#Z\"%&R\"@$3F` pFip>F[q,$F]q#\"$!>\"$P)@$3FhpFdp>F[q,$F]q#\"$1#\"&lw$@$555/FCFA/FGFA/ FKFA/FOFAC&@$Feq>F[q,$F]q#\"#*)\"%u;@$F\\rC$@$Ffr>F[q,$F]q#\"#$)Far@$F ir>F[q,$F]q#!#$)Far@$F_p>F[q,$F]q#!#h\"%!e&@$FgpC$@$Ffr>F[qF_t@$Fir>F[ q,$F]q#\"#hFbtC$@$F[pC%@$F^p>F[qFU@$Feq>F[q,$F]q#\"$!QFjq@$F\\r>F[q,$F ]q#\"$7%Far@$FcrC&@$Feq>F[qF]s@$F\\rC$@$FerC$@$5/&FD6$FAF@F>/&FD6$FAFI F>>F[qFes@$5/&FD6$FAFMF>/&FD6$FAFQF>>F[qFjs@$5FhrFirC$@$Fdv>F[qFjs@$F] w>F[qFes@$F_pC$@$Fdv>F[qF_t@$F]w>F[qFit@$FgpC$@$Fer>F[qF_t@$Ffw>F[qFit F[qF0F0F0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4205 "#This functi on calculates the energy of basis1, basis2. It works in exactly the s ame way as the inner_prod function above,\n#since the energy form is b ilinear, just like the inner product. In fact, the only things which \+ change are the theoretical\n#values (again, found in the splines paper ), and the scaling factor becomes 5/3...\n\nenergy:=proc(S,m,basis_1,b asis_2)\n\n local out,temp,basis1,basis2;\n\n basis1:=basis_1;\n \+ basis2:=basis_2;\n\n#if they aren't neighbors, we get nothing...\n \+ 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;\n\n#first deal with special case of e ither basis1 or basis2 elt of v0:\n\n if basis1[1]=[0] or basis1[1] \+ = [1] or basis1[1] = [2] or basis2[1]=[0] or basis2[1] = [1] or basis2 [1] = [2]\nthen\n # print(basis1,basis2);\n #first make basis 1 the one not in V0...\n if nops(basis1[1])=1 then\n temp :=basis2;\n basis2:=basis1;\n basis1:=temp;\n fi; \n\n#print (basis1,basis2); \n if basis1[1]=basis2[1] then\n \+ if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis 2[2] = f) then out:=-7/18 * ((5/3)^m) fi;\n if basis1[2] = f a nd basis2[2] = f then out:=19/6 * ((5/3)^m) fi;\n if basis1[2 ] = g and basis2[2] = g then out:=5/27 * ((5/3)^m) fi;\n else\n \+ if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or \+ S[basis1[1],3] = basis2[1] or\nS[basis1[1],4] = basis2[1]\nthen\n \+ if basis1[2] = f and basis2[2] = f then out:=-19/12 * ((5/3)^m) fi;\n if basis1[2] = g and basis2[2] = g then \n \+ if S[basis1[1],1] = basis2[1] then out:=-1/108 * ((5/3)^m) fi; \+ \n if S[basis1[1],4] = basis2[1] then out:=1/108 * ((5/ 3)^m) fi;\n fi;\n if (basis1[2] = f and basis2[2 ] = g) then out:=7/36 * ((5/3)^m) fi;\n if (basis1[2] = g a nd basis2[2] = f) then\n if S[basis1[1],1] = basis2[1] \+ then out:=7/36 * ((5/3)^m) fi; \n if S[basis1[1],4] = b asis2[1] then out:=-7/36 * ((5/3)^m) fi;\n fi;\n\n \+ fi;\n fi; \n # fi; \n\n else\n#OTHERWISE...\n\n if bas is1[1]=basis2[1] then\n if (basis1[2] = f and basis2[2] = g) or ( basis1[2] = g and basis2[2] = f) then out:=0 fi;\n if basis1[2] = f and basis2[2] = f then out:=19/6 * 2 * ((5/3)^m) fi;\n if basi s1[2] = g and basis2[2] = g then out:=5/27 * 2 * ((5/3)^m) fi;\n fi; \n if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] or S[ basis1[1],3] = basis2[1] or\nS[basis1[1],4] = basis2[1]\nthen\n i f basis1[2] = f and basis2[2] = f then out:=-19/12 * ((5/3)^m) fi;\n# \+ if basis1[2] = g and basis2[2] = g then \n# if S[basis1[ 1],3] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=-1/108 * ((5 /3)^m) fi; \n# if S[basis1[1],1] = basis2[1] or S[basis1[1],4 ] = basis2[1] then out:=1/108 * ((5/3)^m) fi;\n# fi;\n if ba sis1[2] = g and basis2[2] = g then \n if S[basis1[1],1] = bas is2[1] or S[basis1[1],2] = basis2[1] then\n if S[basis2[1] ,1] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=-1/108 * ((5/3 )^m) fi; \n if S[basis2[1],3] = basis1[1] or S[basis2[1],4 ] = basis1[1] then out:=1/108 * ((5/3)^m) fi;\n fi; \n \+ if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then # out:=-83/37665 / (3^m) fi;\n if S[basis2[1],1] = basis1[1 ] or S[basis2[1],2] = basis1[1] then out:=1/108 * ((5/3)^m) fi; \n \+ if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = basis1[1] t hen out:=-1/108 * ((5/3)^m) fi;\n fi; \n\n fi;\n\n \+ if (basis1[2] = f and basis2[2] = g) then \n if S[basis2[1],1 ] = basis1[1] or S[basis2[1],2] = basis1[1] then out:=7/36 * ((5/3)^m) fi; \n if S[basis2[1],3] = basis1[1] or S[basis2[1],4] = bas is1[1] then out:=-7/36 * ((5/3)^m) fi;\n fi;\n if (basis1[2 ] = g and basis2[2] = f) then\n if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then out:=7/36 * ((5/3)^m) fi;\n \+ if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] then ou t:=-7/36 * ((5/3)^m) fi;\n fi; \n\n \n fi;\n fi; \n \+ out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%'energyGf*6&%\"S G%\"mG%(basis_1G%(basis_2G6&%$outG%%tempG%'basis1G%'basis2G6\"F0C'>8&9 &>8'9'@$33330&F36#\"\"\"&F6F?0&9$6$F>F@FA0&FD6$F>\"\"#FA0&FD6$F>\"\"$F A0&FD6$F>\"\"%FA-%'RETURNG6#\"\"!@%55555/F>7#FU/F>7#F@/F>7#FI/FAFgn/FA Fin/FAF[oC$@$/-%%nopsG6#F>F@C%>8%F6>F6F3>F3Fgo@%/F>FAC%@$53/&F36#FI%\" fG/&F6Fbp%\"gG3/FapFfp/FepFcp>8$,$)#\"\"&FM9%#!\"(\"#=@$3F`pFip>F[q,$F ]q#\"#>\"\"'@$3FhpFdp>F[q,$F]q#F_q\"#F@$555/FCFA/FGFA/FKFA/FOFAC&@$Feq >F[q,$F]q#!#>\"#7@$F\\rC$@$Fer>F[q,$F]q#!\"\"\"$3\"@$Fhr>F[q,$F]q#F@Fg s@$F_p>F[q,$F]q#\"\"(\"#O@$FgpC$@$Fer>F[qF^t@$Fhr>F[q,$F]q#FbqFatC$@$F [pC%@$F^p>F[qFU@$Feq>F[q,$F]q#FiqFM@$F\\r>F[q,$F]q#\"#5F`r@$FbrC&@$Feq >F[qF\\s@$F\\rC$@$FdrC$@$5/&FD6$FAF@F>/&FD6$FAFIF>>F[qFds@$5/&FD6$FAFM F>/&FD6$FAFQF>>F[qFjs@$5FgrFhrC$@$Fav>F[qFjs@$Fjv>F[qFds@$F_pC$@$Fav>F [qF^t@$Fjv>F[qFht@$FgpC$@$Fdr>F[qF^t@$Fcw>F[qFhtF[qF0F0F0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2106 "# After running this function, yo u will end up having:\n# funcno, funcno+2,fucno+4 are (100),(010),(001 ) basis f's (old basis)\n# funcno+6,funcno+8,funcno+10 are Lap(f)=(100 ),etc. for the old basis\n# funcno+1,funcno+5,funcno+9 are f-type base s for the new basis\n# funcno+3,funcno+7,funcno+11 are g-type bases fo r the new basis\n# funcno+12+1,12+5,12+9,12+3...->12+11 are laplacian s of 1,5,9,3,7,11\n# (used for calculating laplacians of splines)\n\nc alc_biharms:=proc(T,R,m_harm,funcno)\n\n#first, let's set up to calcul ate the old basis elements...\n\nT[[0],funcno+0]:=1: T[[1],funcno+0]:= 0: T[[2],funcno+0]:=0;\nT[[0],funcno+1]:=0: T[[1],funcno+1]:=0: T[[2], funcno+1]:=0;\n\nT[[0],funcno+2]:=0: T[[1],funcno+2]:=1: T[[2],funcno+ 2]:=0;\nT[[0],funcno+3]:=0: T[[1],funcno+3]:=0: T[[2],funcno+3]:=0;\n \nT[[0],funcno+4]:=0: T[[1],funcno+4]:=0: T[[2],funcno+4]:=1;\nT[[0],f uncno+5]:=0: T[[1],funcno+5]:=0: T[[2],funcno+5]:=0;\n\n\nT[[0],funcno +6]:=0: T[[1],funcno+6]:=0: T[[2],funcno+6]:=0;\nT[[0],funcno+7]:=1: T [[1],funcno+7]:=0: T[[2],funcno+7]:=0;\n\nT[[0],funcno+8]:=0: T[[1],fu ncno+8]:=0: T[[2],funcno+8]:=0;\nT[[0],funcno+9]:=0: T[[1],funcno+9]:= 1: T[[2],funcno+9]:=0;\n\nT[[0],funcno+10]:=0: T[[1],funcno+10]:=0: T[ [2],funcno+10]:=0;\nT[[0],funcno+11]:=0: T[[1],funcno+11]:=0: T[[2],fu ncno+11]:=1;\n\n#Now, we'll actually calculate the bi-harmonic functio ns, leaving the laplacian of the function\n#in funcno+1... (hence the \+ need to space them out by 2).\nn_harm(T,R,2,m_harm,[0],[1],[2],funcno) ;\nn_harm(T,R,2,m_harm,[0],[1],[2],funcno+2);\nn_harm(T,R,2,m_harm,[0] ,[1],[2],funcno+4);\nn_harm(T,R,2,m_harm,[0],[1],[2],funcno+6);\nn_har m(T,R,2,m_harm,[0],[1],[2],funcno+8);\nn_harm(T,R,2,m_harm,[0],[1],[2] ,funcno+10);\n\n#Here, we will apply the transforms to the new basis t o the laplacians of the old basis, thereby\n#creating the laplacians o f the new basis functions...\nf_g_lap(T,m_harm,funcno);\n\n#and lastly , we will create the new basis elements, sticking their values in betw een the old basis elts (to save memory).\n\nf_g_bih(T,m_harm,0,funcno, funcno+1);\nf_g_bih(T,m_harm,2,funcno,funcno+5);\nf_g_bih(T,m_harm,4,f uncno,funcno+9);\n\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-calc_bi harmsGf*6&%\"TG%\"RG%'m_harmG%'funcnoG6\"F+F+CP>&9$6$7#\"\"!9'\"\"\">& F/6$7#F4F3F2>&F/6$7#\"\"#F3F2>&F/6$F1,&F3F4F4F4F2>&F/6$F8FAF2>&F/6$F&F/6$F1,&F3F4F=F4F2>&F/6$F8FKF4>&F/6$F&F/6$F1,&F3F4\"\"$F4F2 >&F/6$F8FUF2>&F/6$F&F/6$F1,&F3F4\"\"%F4F2>&F/6$F8FjnF2>&F/6$F&F/6$F1,&F3F4\"\"&F4F2>&F/6$F8FeoF2>&F/6$F&F/6$F1,&F3F4\"\" 'F4F2>&F/6$F8F`pF2>&F/6$F&F/6$F1,&F3F4\"\"(F4F4>&F/6$F8F[qF2>&F /6$F&F/6$F1,&F3F4\"\")F4F2>&F/6$F8FfqF2>&F/6$F&F/6$F1,&F 3F4\"\"*F4F2>&F/6$F8FarF4>&F/6$F&F/6$F1,&F3F4\"#5F4F2>&F/6$F8F \\sF2>&F/6$F&F/6$F1,&F3F4\"#6F4F2>&F/6$F8FgsF2>&F/6$F " 0 "" {MPLTEXT 1 0 701 "#This function is used almost everywhere. Basically, let's say y ou have an address\n#which is given relative to a subtriangle, and the n want to find out what the address\n#of that point is relative to the total triangle. Then you would call this function,\n#with input the \+ relative address in triangle, and with add0,add1,add2 being the addres ses\n#of the 3 vertices of triangle (relative to the main triangle).\n \ncoord_transfer:=proc(input, triangle,add0,add1,add2)\n\nif nops(inpu t) = 1 then\n if input = [0] then RETURN(add0); fi;\n if input = [1] then RETURN(add1); fi;\n if input = [2] then RETURN(add2); fi; \nfi;\n\n \nif nops(input) > 1 then \n RETURN([op(triangle),op(i nput)]); fi;\n\nend; \n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/coor d_transferGf*6'%&inputG%)triangleG%%add0G%%add1G%%add2G6\"F,F,C$@$/-%% nopsG6#9$\"\"\"C%@$/F37#\"\"!-%'RETURNG6#9&@$/F37#F4-F;6#9'@$/F37#\"\" #-F;6#9(@$2F4F0-F;6#7$-%#opG6#9%-FQF2F,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2123 "#Basically, given a basis element, it is often u seful to determine the\n#two neighboring triangles (often used in coor d_transfer). Thus, this function\n#takes a basis element as its input , and if which_tri=1, it returns the triangle\n#address of the first t riangle, and if which_tri=2, it returns the triangle address\n#of the \+ second triangle. Basically, if the basis element is an element of Vm_ basis,\n#then we can just cut off the last 2 digits and tack on a 0 or a 1 or a 2, depending\n#on the orientation of the basis element (ie, \+ the last 2 digits). Otherwise,\n#we can use the fact that the neighbo rs must be in Vm_basis, so we can choose the \n#appropriate neighbor a nd chop off some digits, and so forth...\n\nget_triangle:=proc(S,basis ,m_basis,which_tri)\n\nif which_tri = 1 then\n\nif basis[1][nops(basis [1])-1..nops(basis[1])] = [0,1] then \n if nops(basis[1]) = m_basis+ 1 then RETURN([op(basis[1][1..nops(basis[1])-2]),1]);\n else RETURN( [op(S[basis[1],1][1..nops(S[basis[1],1])-2]),0]); fi;\nfi;\n\nif basis [1][nops(basis[1])-1..nops(basis[1])] = [1,2] then \n if nops(basis[ 1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),2]);\n \+ else RETURN([op(S[basis[1],1][1..nops(S[basis[1],1])-2]),1]); fi;\nf i;\n\nif basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then \n \+ if nops(basis[1]) = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1 ])-2]),0]);\n else RETURN([op(S[basis[1],1][1..nops(S[basis[1],1])-2 ]),2]); fi;\nfi;\n\n\nelse\n\nif basis[1][nops(basis[1])-1..nops(basis [1])] = [0,1] then \n if nops(basis[1]) = m_basis+1 then RETURN([op( basis[1][1..nops(basis[1])-2]),0]);\n else RETURN([op(S[basis[1],4][ 1..nops(S[basis[1],4])-2]),1]); fi;\nfi;\n\nif basis[1][nops(basis[1]) -1..nops(basis[1])] = [1,2] then \n if nops(basis[1]) = m_basis+1 th en RETURN([op(basis[1][1..nops(basis[1])-2]),1]);\n else RETURN([op( S[basis[1],4][1..nops(S[basis[1],4])-2]),2]); fi;\nfi;\n\nif basis[1][ nops(basis[1])-1..nops(basis[1])] = [0,2] then \n if nops(basis[1]) \+ = m_basis+1 then RETURN([op(basis[1][1..nops(basis[1])-2]),2]);\n el se RETURN([op(S[basis[1],4][1..nops(S[basis[1],4])-2]),0]); fi;\nfi;\n \nfi;\nend; \n \n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-get _triangleGf*6&%\"SG%&basisG%(m_basisG%*which_triG6\"F+F+@%/9'\"\"\"C%@ $/&&9%6#F/6#;,&-%%nopsG6#F4F/F/!\"\"F:7$\"\"!F/@%/F:,&9&F/F/F/-%'RETUR NG6#7$-%#opG6#&F46#;F/,&F:F/\"\"#F=F/-FE6#7$-FI6#&&9$6$F4F/6#;F/,&-F;6 #FVF/FOF=F?@$/F37$F/FO@%FA-FE6#7$FHFO-FE6#7$FSF/@$/F37$F?FO@%FA-FE6#7$ FHF?-FE6#7$FSFOC%@$F2@%FAFfo-FE6#7$-FI6#&&FW6$F4\"\"%6#;F/,&-F;6#FepF/ FOF=F/@$Fin@%FAFD-FE6#7$FbpFO@$Fco@%FAF\\o-FE6#7$FbpF?F+F+F+" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10819 "#This procedure takes a f unction and returns integral(funct*basis).\n#the basic strategy is to \+ determine which triangles basis is non zero on\n#and then use coord_tr ansfer to determine which values of funct to multiply\n#the value of t he appropriate basis function by and then integrate.\n\n\nint_vs_basis :=proc(T,S,m_basis,m_func,basis,funct,harm_funcno)\n\nlocal V0,Vother, Vm,fact,triangle1,triangle2,V0sum1,V0sum2,Vothersum1,Vothersum2,Vmsum1 ,Vmsum2,i,j,Vother_e;\n\nfact:= 1/(6*3^m_func); #factor for simpson's \+ rule...\n\n#get the vertices (ie, V0, Vm, Veverything in between)...\n V0:=vertices_m(0);\nif m_func-m_basis=1 then Vother:=[] else Vother_e: =vertices(m_func-m_basis-1); Vother:=Vother_e[4..nops(Vother_e)];fi;\n Vm:=vertices_m(m_func-m_basis);\n\n\nif nops(basis[1]) = 1 then #this \+ is if the basis is an element of V0...\n#Here, we basically know the b asis element is a g type vector. Then, we use the fact that\n#the tri angle in question is just a sequence of 0's or 1's. Then, we multiply by the appropriate\n#biharmonic basis function (ie, harm_funcno+3 or \+ some such number), and sum...\n if basis[1] = [0] then\n tri angle1:=[seq(0,j=1..m_basis)];\n\n 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));\n\n Vothersum1:=add(T[coord_ transfer(Vother[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),fun ct] *( T[Vother[i],harm_funcno+3]),i=1..nops(Vother)); \n \n \+ 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));\n\n \+ fi;\n if basis[1] = [1] then\n triangle1:=[seq(1,j=1..m_bas is)];\n\n 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..no ps(V0));\n\n Vothersum1:=add(T[coord_transfer(Vother[i],triangl e1,S[basis[1],2],basis[1],S[basis[1],1]),funct] *( T[Vother[i],harm_fu ncno+7]),i=1..nops(Vother)); \n \n Vmsum1:=add(T[coord_tra nsfer(Vm[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*( T [Vm[i],harm_funcno+7]),i=1..nops(Vm));\n\n fi;\n if basis[1] = [ 2] then\n triangle1:=[seq(2,j=1..m_basis)];\n\n 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));\n\n Vothe rsum1:=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(Vothe r)); \n \n 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));\n fi; \n \nRETURN(fact * (5*(Vmsum1) + 2*(Vothersum 1) + (V0sum1))); #And then return the simpson's rule formula...\n\n\n else #So we will come here if the basis element is not on the boundary (could be either f or g...)\n\n#first let's deal with the basis eleme nt being an f...\nif basis[2]=f then\n#So now we have to get 2 differe nt triangles (the two triangles associated with the vertex)\n#Hence we have 2 calls to get_triangle. So we have 2 sums, one for triangle 1 \+ and one for triangle2.\n#We also need to keep track of the orientation (ie, the last 2 digits of the basis) in order to \n#determine which b asis functions we want to use.\n if basis[1][nops(basis[1])-1..nops (basis[1])] = [0,1] then \n triangle1:=get_triangle(S,basis,m_bas is,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n #le t's calculate sum's for triangle 1, shall we?\n V0sum1:=add(T[coo rd_transfer(V0[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),func t]*(T[V0[i],harm_funcno+1] ),i=1..nops(V0));\n\n 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)); \n \+ \n 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)); \n\n #and now on to triangle numero 2!\n 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));\n\n Vothersum2:=add(T[c oord_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)); \n\n 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));\n\n fi;\n\n if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] t hen\n triangle1:=get_triangle(S,basis,m_basis,1);\n triangle 2:=get_triangle(S,basis,m_basis,2);\n\n #let's calculate sum's fo r triangle 1, shall we?\n V0sum1:=add(T[coord_transfer(V0[i],tria ngle1,S[basis[1],1],S[basis[1],2],basis[1]),funct]*(T[V0[i],harm_funcn o+9]),i=1..nops(V0));\n\n 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)); \n \n Vmsum1:=add(T[c oord_transfer(Vm[i],triangle1,S[basis[1],1],S[basis[1],2],basis[1]),fu nct]*(T[Vm[i],harm_funcno+9]),i=1..nops(Vm));\n\n #and now on to \+ triangle numero 2!\n 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));\n\n Vothersum2:=add(T[coord_transfer(Vother[i],t riangle2,basis[1],S[basis[1],3],S[basis[1],4]),funct] *(T[Vother[i],ha rm_funcno+1]),i=1..nops(Vother)); \n\n 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));\n\n fi; \n\n if basis[ 1][nops(basis[1])-1..nops(basis[1])] = [1,2] then\n triangle1:=ge t_triangle(S,basis,m_basis,1);\n triangle2:=get_triangle(S,basis, m_basis,2);\n\n #let's calculate sum's for triangle 1, shall we? \n V0sum1:=add(T[coord_transfer(V0[i],triangle1,S[basis[1],2],bas is[1],S[basis[1],1]),funct]*(T[V0[i],harm_funcno+5]),i=1..nops(V0));\n \n 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)); \n \n Vmsum1:=add(T[coord_transfer(Vm[i],tri angle1,S[basis[1],2],basis[1],S[basis[1],1]),funct]*(T[Vm[i],harm_func no+5]),i=1..nops(Vm));\n\n #and now on to triangle numero 2!\n \+ 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));\n\n \+ 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)); \n\n Vmsum2:=add(T[coord_transfer(Vm[i],triangl e2,S[basis[1],3],S[basis[1],4],basis[1]),funct]*(T[Vm[i],harm_funcno+9 ]),i=1..nops(Vm));\n\n fi; \nfi;\n\n#now let us deal with the basis element being a g...\nif basis[2]=g then\n#The code is exactly the sa me as for the f-type vector, except that \n#we use g-type basis vector s (which basically means adding 2 to the\n#harm_funcno+k part).\n i f basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1] then \n tria ngle1:=get_triangle(S,basis,m_basis,1);\n triangle2:=get_triangle (S,basis,m_basis,2);\n\n #let's calculate sum's for triangle 1, s hall we?\n 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));\n\n Vothersum1:=add(T[coord_transfer(Vother[i],triangle1,b asis[1],S[basis[1],1],S[basis[1],2]),funct] *(T[Vother[i],harm_funcno+ 3]),i=1..nops(Vother)); \n \n Vmsum1:=add(T[coord_transfer(V m[i],triangle1,basis[1],S[basis[1],1],S[basis[1],2]),funct]*(T[Vm[i],h arm_funcno+3]),i=1..nops(Vm));\n\n #and now on to triangle numero 2!\n 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)) ;\n\n Vothersum2:=add(T[coord_transfer(Vother[i],triangle2,S[basi s[1],4],basis[1],S[basis[1],3]),funct] *(T[Vother[i],harm_funcno+7]),i =1..nops(Vother)); \n\n 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));\n\n fi;\n\n if basis[1][nops(basis[1] )-1..nops(basis[1])] = [0,2] then\n triangle1:=get_triangle(S,bas is,m_basis,1);\n triangle2:=get_triangle(S,basis,m_basis,2);\n\n \+ #let's calculate sum's for triangle 1, shall we?\n V0sum1:=a dd(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));\n\n Vothersu m1:=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)); \n \n 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..no ps(Vm));\n\n #and now on to triangle numero 2!\n 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));\n\n Vothersum2 :=add(T[coord_transfer(Vother[i],triangle2,basis[1],S[basis[1],3],S[ba sis[1],4]),funct] *(T[Vother[i],harm_funcno+3]),i=1..nops(Vother)); \+ \n\n 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));\n\n fi; \n\n if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then\n triangle1:=get_triangle(S,basis,m_basis,1);\n \+ triangle2:=get_triangle(S,basis,m_basis,2);\n\n #let's calculat e sum's for triangle 1, shall we?\n 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));\n\n Vothersum1:=add(T[coord_trans fer(Vother[i],triangle1,S[basis[1],2],basis[1],S[basis[1],1]),funct] * (T[Vother[i],harm_funcno+7]),i=1..nops(Vother)); \n \n Vmsum 1:=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));\n\n #and \+ now on to triangle numero 2!\n 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));\n\n 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));\n \n Vmsum2 :=add(T[coord_transfer(Vm[i],triangle2,S[basis[1],3],S[basis[1],4],bas is[1]),funct]*(T[Vm[i],harm_funcno+11]),i=1..nops(Vm));\n\n fi; \n \nfi;\n\n#Return simpson's rule. Notice that the second sum must be s ubtracted if it's a g-type vector\n#because the normal derivative is d efined as negative in the 2nd triangle...\nif basis[2] = f then RETURN (fact * (5*(Vmsum1+Vmsum2) + 2*(Vothersum1+Vothersum2) + (V0sum1+V0sum 2))); fi;\nif basis[2] = g then RETURN(fact * (5*(Vmsum1-Vmsum2) + 2*( Vothersum1-Vothersum2) + (V0sum1-V0sum2))); fi;\n\nfi;\n\n\nend;\n \+ \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-int_vs_basisGf*6)%\"TG% \"SG%(m_basisG%'m_funcG%&basisG%&functG%,harm_funcnoG61%#V0G%'VotherG% #VmG%%factG%*triangle1G%*triangle2G%'V0sum1G%'V0sum2G%+Vothersum1G%+Vo thersum2G%'Vmsum1G%'Vmsum2G%\"iG%\"jG%)Vother_eG6\"F>C'>8',$*&\"\"\"FD )\"\"$9'!\"\"#FD\"\"'>8$-%+vertices_mG6#\"\"!@%/,&FGFD9&FHFD>8%7\"C$>8 2-%)verticesG6#,(FGFDFTFHFDFH>FV&FZ6#;\"\"%-%%nopsG6#FZ>8&-FN6#FS@%/-F _o6#&9(6#FDFDC&@$/Fio7#FPC&>8(7#-%$seqG6$FP/81;FDFT>8*-%$addG6$*&&9$6$ -%/coord_transferG6'&FL6#80FbpFio&9%6$FioFD&Fjq6$Fio\"\"#9)FD&Faq6$Ffq ,&9*FDFFFDFD/Fhq;FD-F_o6#FL>8,-F]q6$*&&Faq6$-Fdq6'&FVFgqFbpFioFiqF\\rF _rFD&Faq6$FasFbrFD/Fhq;FD-F_o6#FV>8.-F]q6$*&&Faq6$-Fdq6'&FboFgqFbpFioF iqF\\rF_rFD&Faq6$FatFbrFD/Fhq;FD-F_o6#Fbo@$/Fio7#FDC&>Fbp7#-Fep6$FDFgp >F[q-F]q6$*&&Faq6$-Fdq6'FfqFbpF\\rFioFiqF_rFD&Faq6$Ffq,&FcrFD\"\"(FDFD Fdr>Fir-F]q6$*&&Faq6$-Fdq6'FasFbpF\\rFioFiqF_rFD&Faq6$FasFjuFDFds>Fis- F]q6$*&&Faq6$-Fdq6'FatFbpF\\rFioFiqF_rFD&Faq6$FatFjuFDFdt@$/Fio7#F^rC& >Fbp7#-Fep6$F^rFgp>F[q-F]q6$*&&Faq6$-Fdq6'FfqFbpFiqF\\rFioF_rFD&Faq6$F fq,&FcrFD\"#6FDFDFdr>Fir-F]q6$*&&Faq6$-Fdq6'FasFbpFiqF\\rFioF_rFD&Faq6 $FasFbxFDFds>Fis-F]q6$*&&Faq6$-Fdq6'FatFbpFiqF\\rFioF_rFD&Faq6$FatFbxF DFdt-%'RETURNG6#*&FAFD,(Fis\"\"&*&F^rFDFirFDFDF[qFDFDC&@$/&Fjo6#F^r%\" fGC%@$/&Fio6#;,&FgoFDFDFHFgo7$FPFDC*>Fbp-%-get_triangleG6&FjqFjoFTFD>8 )-F`[l6&FjqFjoFTF^r>F[q-F]q6$*&F`qFD&Faq6$Ffq,&FcrFDFDFDFDFdr>Fir-F]q6 $*&F]sFD&Faq6$FasF\\\\lFDFds>Fis-F]q6$*&F]tFD&Faq6$FatF\\\\lFDFdt>8+-F ]q6$*&&Faq6$-Fdq6'FfqFc[l&Fjq6$FioF]oFio&Fjq6$FioFFF_rFD&Faq6$Ffq,&Fcr FDF]zFDFDFdr>8--F]q6$*&&Faq6$-Fdq6'FasFc[lFb]lFioFd]lF_rFD&Faq6$FasFh] lFDFds>8/-F]q6$*&&Faq6$-Fdq6'FatFc[lFb]lFioFd]lF_rFD&Faq6$FatFh]lFDFdt @$/Fhz7$FPF^rC*>FbpF_[l>Fc[lFd[l>F[q-F]q6$*&F\\xFD&Faq6$Ffq,&FcrFD\"\" *FDFDFdr>Fir-F]q6$*&FhxFD&Faq6$FasF[`lFDFds>Fis-F]q6$*&FbyFD&Faq6$FatF [`lFDFdt>Fj\\l-F]q6$*&&Faq6$-Fdq6'FfqFc[lFioFd]lFb]lF_rFDFj[lFDFdr>Fj] l-F]q6$*&&Faq6$-Fdq6'FasFc[lFioFd]lFb]lF_rFDFa\\lFDFds>Fe^l-F]q6$*&&Fa q6$-Fdq6'FatFc[lFioFd]lFb]lF_rFDFg\\lFDFdt@$/Fhz7$FDF^rC*>FbpF_[l>Fc[l Fd[l>F[q-F]q6$*&FduFDFf]lFDFdr>Fir-F]q6$*&F`vFDFb^lFDFds>Fis-F]q6$*&Fj vFDF]_lFDFdt>Fj\\l-F]q6$*&&Faq6$-Fdq6'FfqFc[lFd]lFb]lFioF_rFDFi_lFDFdr >Fj]l-F]q6$*&&Faq6$-Fdq6'FasFc[lFd]lFb]lFioF_rFDFa`lFDFds>Fe^l-F]q6$*& &Faq6$-Fdq6'FatFc[lFd]lFb]lFioF_rFDFg`lFDFdt@$/Fbz%\"gGC%@$FgzC*>FbpF_ [l>Fc[lFd[l>F[qF\\q>FirFjr>FisFjs>Fj\\l-F]q6$*&F^]lFDFhuFDFdr>Fj]l-F]q 6$*&F^^lFDFdvFDFds>Fe^l-F]q6$*&Fi^lFDF^wFDFdt@$F`_lC*>FbpF_[l>Fc[lFd[l >F[qFiw>FirFex>FisF_y>Fj\\l-F]q6$*&F]alFDF`rFDFdr>Fj]l-F]q6$*&FealFDFb sFDFds>Fe^l-F]q6$*&F]blFDFbtFDFdt@$FbblC*>FbpF_[l>Fc[lFd[l>F[qFau>FirF ]v>FisFgv>Fj\\l-F]q6$*&FgclFDF`xFDFdr>Fj]l-F]q6$*&F_dlFDF\\yFDFds>Fe^l -F]q6$*&FgdlFDFfyFDFdt@$Faz-Fiy6#*&FAFD,.FisF]z*&F]zFDFe^lFDFD*&F^rFDF irFDFD*&F^rFDFj]lFDFDF[qFDFj\\lFDFD@$F\\el-Fiy6#*&FAFD,.FisF]z*&F]zFDF e^lFDFH*&F^rFDFirFDFD*&F^rFDFj]lFDFHF[qFDFj\\lFHFDF>F>F>" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11448 "#so this function returns integr al(funct*basis1*basis2).\n#it is used to calculate the q-matrix. This function basically is like\n#a combination of the logic of inner_prod with the integration of int_vs_basis.\n#several help functions are de fined to help do this (it's a pretty complicated function).\n\nint_vs_ basis12:=proc(T,S,m_basis,m_func,basis_1,basis_2,funct,harm_funcno)\n \n local out,temp,basis1,basis2,fact,triangle1,triangle2,addresses1, addresses2,j,which_harms;\n\n fact:=1/(3^m_basis); #integration scal ing factor...\n basis1:=basis_1; #again, so that we can switch arou nd basis1 and basis2, if need be.\n basis2:=basis_2;\n\n#if they are n't neighbors, we get nothing...\n if basis1[1] <> basis2[1] and S[b asis1[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; \n\n#first deal with special case of either basis1 or basis2 elt of v0 :\n\n if basis1[1]=[0] or basis1[1] = [1] or basis1[1] = [2] or basi s2[1]=[0] or basis2[1] = [1] or basis2[1] = [2] then\n #first mak e basis1 the one not in V0...\n #ie, make basis2 the one which is in V0 for sure\n if nops(basis1[1])=1 then\n temp:=basis 2;\n basis2:=basis1;\n basis1:=temp;\n fi;\n\n \+ if basis1[1]=basis2[1] then #The only way we'll ever get here is \+ if both elements are g-type basis vectors\n if basis1[2] = g a nd basis2[2] = g then \n if basis1[1] = [0] then \n \+ triangle1:=[seq(0,j=1..m_basis)];\n #The function integrate_basis12 basically does all the work; alls we have to do is \+ tell it which functions...\n out:=integrate_basis12(T,m_ func-m_basis,triangle1,basis1[1],S[basis1[1],1],S[basis1[1],2],funct,h arm_funcno+3,harm_funcno+3);\n fi;\n if basis1[1 ] = [1] then\n triangle1:=[seq(1,j=1..m_basis)];\n \+ 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); \n \+ fi;\n if basis1[1] = [2] then\n tria ngle1:=[seq(2,j=1..m_basis)];\n out:=integrate_basis12(T ,m_func-m_basis,triangle1,S[basis1[1],1],S[basis1[1],2],basis1[1],func t,harm_funcno+11,harm_funcno+11);\n fi; \n \+ fi;\n else\n 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] = basi s2[1] then\n # Notice that the only way we'll get in here is i f basis2 is a g-type vector; however, basis1 can still be either g or \+ f... \n if basis1[2] = g and basis2[2] = g then\n \+ if S[basis1[1],1] = basis2[1] then\n if basis2[1] = [0] then \n triangle1:=[seq(0,j=1..m_basis)];\n \+ out:=integrate_basis12(T,m_func-m_basis,triangle1,bas is2[1],S[basis2[1],1],S[basis2[1],2],funct,harm_funcno+3,harm_funcno+1 1);\n fi;\n if basis2[1] = [1] then\n \+ triangle1:=[seq(1,j=1..m_basis)];\n \+ out:=integrate_basis12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis 2[1],S[basis2[1],1],funct,harm_funcno+7,harm_funcno+3); \n \+ fi;\n if basis2[1] = [2] then\n t riangle1:=[seq(2,j=1..m_basis)];\n out:=integrate_ba sis12(T,m_func-m_basis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[ 1],funct,harm_funcno+11,harm_funcno+7);\n fi;\n \+ fi;\n \n if S[basis1[1],4] = basis2[1] the n\n if basis2[1] = [0] then \n triang le1:=[seq(0,j=1..m_basis)];\n out:=-integrate_basis1 2(T,m_func-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],f unct,harm_funcno+3,harm_funcno+7);\n fi;\n \+ if basis2[1] = [1] then\n triangle1:=[seq(1,j=1.. m_basis)];\n out:=-integrate_basis12(T,m_func-m_basi s,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+ 7,harm_funcno+11); \n fi;\n if basis2[1] = [2] then\n triangle1:=[seq(2,j=1..m_basis)];\n \+ out:=-integrate_basis12(T,m_func-m_basis,triangle1,S[b asis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcno+11,harm_funcno+ 3);\n fi;\n fi;\n fi;\n \+ \n if (basis1[2] = f and basis2[2] = g) then #Thi s is the only remaining case, because f-type vectors can never be at t he boundaries\n if S[basis1[1],1] = basis2[1] then\n \+ if basis2[1] = [0] then \n triangle1:=[ seq(0,j=1..m_basis)];\n out:=integrate_basis12(T,m_f unc-m_basis,triangle1,basis2[1],S[basis2[1],1],S[basis2[1],2],funct,ha rm_funcno+3,harm_funcno+9);\n fi;\n if b asis2[1] = [1] then\n triangle1:=[seq(1,j=1..m_basis )];\n out:=integrate_basis12(T,m_func-m_basis,triang le1,S[basis2[1],2],basis2[1],S[basis2[1],1],funct,harm_funcno+7,harm_f uncno+1); \n fi;\n if basis2[1] = [2] th en\n triangle1:=[seq(2,j=1..m_basis)];\n \+ 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);\n \+ fi;\n fi;\n \n if S[basis1[1] ,4] = basis2[1] then\n\n if basis2[1] = [0] then \n \+ triangle1:=[seq(0,j=1..m_basis)];\n o ut:=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);\n \+ fi;\n if basis2[1] = [1] then\n tri angle1:=[seq(1,j=1..m_basis)];\n out:=integrate_basi s12(T,m_func-m_basis,triangle1,S[basis2[1],2],basis2[1],S[basis2[1],1] ,funct,harm_funcno+7,harm_funcno+9); \n fi;\n \+ if basis2[1] = [2] then\n triangle1:=[seq(2,j= 1..m_basis)];\n out:=integrate_basis12(T,m_func-m_ba sis,triangle1,S[basis2[1],1],S[basis2[1],2],basis2[1],funct,harm_funcn o+11,harm_funcno+1);\n fi;\n fi;\n\n\n\n fi;\n\n fi;\n fi; \n # fi; \n\n else \+ #i.e, neither are in V0....\n\n if basis1[1]=basis2[1] then #they' re at the same vertex.\n #here, the only thing that matters is th e orientation to determine which basis functions to use...\n tria ngle1:=get_triangle(S,basis1,m_basis,1); #first, let's get the trian gles, as well as the addresses of the vertices of those triangles\n \+ triangle2:=get_triangle(S,basis1,m_basis,2);\n addresses1:=get _addresses(S,basis1,m_basis,1);\n addresses2:=get_addresses(S,bas is1,m_basis,2);\n if (basis1[2] = f and basis2[2] = g) or (basis1 [2] = g and basis2[2] = f) then #out:=0 fi;\n if last_digits(b asis1[1]) = [0,1] then\n out:=integrate_basis12(T,m_func-m_bas is,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_func no+1,harm_funcno+3)-\n integrate_basis12(T,m_func-m_basis ,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno +5,harm_funcno+7); fi;\n\n if last_digits(basis1[1]) = [0,2] t hen\n out:=integrate_basis12(T,m_func-m_basis,triangle1,addres ses1[1],addresses1[2],addresses1[3],funct,harm_funcno+9,harm_funcno+11 )-\n integrate_basis12(T,m_func-m_basis,triangle2,address es2[1],addresses2[2],addresses2[3],funct,harm_funcno+1,harm_funcno+3); fi;\n\n if last_digits(basis1[1]) = [1,2] then\n out: =integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addresses1 [2],addresses1[3],funct,harm_funcno+5,harm_funcno+7)-\n i ntegrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2 ],addresses2[3],funct,harm_funcno+9,harm_funcno+11); fi;\n fi;\n \n\n\n if basis1[2] = f and basis2[2] = f then #out:=190/837 * 2 \+ / (3^m) fi;\n if last_digits(basis1[1]) = [0,1] then\n \+ out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addre sses1[2],addresses1[3],funct,harm_funcno+1,harm_funcno+1)+\n \+ integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],address es2[2],addresses2[3],funct,harm_funcno+5,harm_funcno+5); fi;\n\n \+ if last_digits(basis1[1]) = [0,2] then\n out:=integrate_bas is12(T,m_func-m_basis,triangle1,addresses1[1],addresses1[2],addresses1 [3],funct,harm_funcno+9,harm_funcno+9)+\n integrate_basis 12(T,m_func-m_basis,triangle2,addresses2[1],addresses2[2],addresses2[3 ],funct,harm_funcno+1,harm_funcno+1); fi;\n\n if last_digits(b asis1[1]) = [1,2] then\n out:=integrate_basis12(T,m_func-m_bas is,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm_func no+5,harm_funcno+5)+\n integrate_basis12(T,m_func-m_basis ,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_funcno +9,harm_funcno+9); fi;\n fi;\n\n\n\n if basis1[2] = g and ba sis2[2] = g then #out:=206/37665 * 2 / (3^m) fi;\n if last_dig its(basis1[1]) = [0,1] then\n out:=integrate_basis12(T,m_func- m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct,harm _funcno+3,harm_funcno+3)+\n integrate_basis12(T,m_func-m_ basis,triangle2,addresses2[1],addresses2[2],addresses2[3],funct,harm_f uncno+7,harm_funcno+7); fi;\n\n if last_digits(basis1[1]) = [0 ,2] then\n out:=integrate_basis12(T,m_func-m_basis,triangle1,a ddresses1[1],addresses1[2],addresses1[3],funct,harm_funcno+11,harm_fun cno+11)+\n integrate_basis12(T,m_func-m_basis,triangle2,a ddresses2[1],addresses2[2],addresses2[3],funct,harm_funcno+3,harm_func no+3); fi;\n\n if last_digits(basis1[1]) = [1,2] then\n \+ out:=integrate_basis12(T,m_func-m_basis,triangle1,addresses1[1],addr esses1[2],addresses1[3],funct,harm_funcno+7,harm_funcno+7)+\n \+ integrate_basis12(T,m_func-m_basis,triangle2,addresses2[1],addres ses2[2],addresses2[3],funct,harm_funcno+11,harm_funcno+11); fi;\n \+ fi;\n fi; #END OF CASE WHERE BASIS1 and BASIS2 are in same place \+ \n 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\n #if they are neighbors, then we're in here. Notice that we just call\n \+ #which_funcs to determine which basis functions to use, as well as getting the negative sign to make the\n #sum negative if we are \+ in a negative normal derivative region...\n if S[basis1[1],1 ] = basis2[1] or S[basis1[1],2] = basis2[1] then\n triang le1:=get_triangle(S,basis1,m_basis,1); #first, let's get the triangl es, as well as the addresses of the vertices of those triangles\n \+ addresses1:=get_addresses(S,basis1,m_basis,1);\n \+ which_harms:=which_funcs(S,basis1,basis2);\n 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_fun cno+which_harms[2]);\n else\n triangle1:=get_tr iangle(S,basis1,m_basis,2);\n addresses1:=get_addresses(S ,basis1,m_basis,2);\n which_harms:=which_funcs(S,basis1,b asis2);\n out:=which_harms[3] * integrate_basis12(T,m_fun c-m_basis,triangle1,addresses1[1],addresses1[2],addresses1[3],funct, h arm_funcno+which_harms[1],harm_funcno+which_harms[2]);\n fi; \n\n fi;\n\n \n fi;\n \n out*fact; \nend;\n\n\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/int_vs_basis12Gf*6*%\"TG%\"SG%(m_ basisG%'m_funcG%(basis_1G%(basis_2G%&functG%,harm_funcnoG6-%$outG%%tem pG%'basis1G%'basis2G%%factG%*triangle1G%*triangle2G%+addresses1G%+addr esses2G%\"jG%,which_harmsG6\"F;C(>8(*&\"\"\"F@)\"\"$9&!\"\">8&9(>8'9)@ $33330&FF6#F@&FIFR0&9%6$FQF@FS0&FV6$FQ\"\"#FS0&FV6$FQFBFS0&FV6$FQ\"\"% FS-%'RETURNG6#\"\"!@%55555/FQ7#F`o/FQ7#F@/FQ7#Fen/FSFho/FSFjo/FSF\\pC$ @$/-%%nopsG6#FQF@C%>8%FI>FIFF>FFFhp@%/FQFS@$3/&FF6#Fen%\"gG/&FIFaqFbqC %@$FgoC$>8)7#-%$seqG6$F`o/8-;F@FC>8$-%2integrate_basis12G6+9$,&9'F@FCF DFiqFQFUFY9*,&9+F@FBF@Fjr@$FioC$>Fiq7#-F\\r6$F@F^r>Fbr-Fdr6+FfrFgrFiqF YFQFUFir,&F[sF@\"\"(F@Fes@$F[pC$>Fiq7#-F\\r6$FenF^r>Fbr-Fdr6+FfrFgrFiq FUFYFQFir,&F[sF@\"#6F@F`t@$555/FUFS/FYFS/FgnFS/FjnFSC$@$F^qC$@$FftC%@$ F]pC$>FiqFjq>Fbr-Fdr6+FfrFgrFiqFS&FV6$FSF@&FV6$FSFenFirFjrF`t@$F^pC$>F iqF_s>Fbr-Fdr6+FfrFgrFiqFguFSFeuFirFesFjr@$F_pC$>FiqFjs>Fbr-Fdr6+FfrFg rFiqFeuFguFSFirF`tFes@$FitC%@$F]pC$>FiqFjq>Fbr,$-Fdr6+FfrFgrFiqFSFeuFg uFirFjrFesFD@$F^pC$>FiqF_s>Fbr,$-Fdr6+FfrFgrFiqFguFSFeuFirFesF`tFD@$F_ pC$>FiqFjs>Fbr,$-Fdr6+FfrFgrFiqFeuFguFSFirF`tFjrFD@$3/F`q%\"fGFcqC$@$F ftC%@$F]pC$>FiqFjq>Fbr-Fdr6+FfrFgrFiqFSFeuFguFirFjr,&F[sF@\"\"*F@@$F^p C$>FiqF_s>Fbr-Fdr6+FfrFgrFiqFguFSFeuFirFes,&F[sF@F@F@@$F_pC$>FiqFjs>Fb r-Fdr6+FfrFgrFiqFeuFguFSFirF`t,&F[sF@\"\"&F@@$FitC%@$F]pC$>FiqFjq>Fbr- Fdr6+FfrFgrFiqFSFeuFguFirFjrFhy@$F^pC$>FiqF_s>Fbr-Fdr6+FfrFgrFiqFguFSF euFirFesFix@$F_pC$>FiqFjs>Fbr-Fdr6+FfrFgrFiqFeuFguFSFirF`tFayC$@$F\\qC )>Fiq-%-get_triangleG6&FVFFFCF@>8*-Fc[l6&FVFFFCFen>8+-%.get_addressesG Fd[l>8,-F\\\\lFh[l@$5F]x3F_q/FdqF_xC%@$/-%,last_digitsGFep7$F`oF@>Fbr, &-Fdr6+FfrFgrFiq&Fj[lFR&Fj[lFaq&Fj[l6#FBFirFayFjrF@-Fdr6+FfrFgrFf[l&F^ \\lFR&F^\\lFaq&F^\\lFa]lFirFhyFesFD@$/Fg\\l7$F`oFen>Fbr,&-Fdr6+FfrFgrF iqF^]lF_]lF`]lFirFixF`tF@-Fdr6+FfrFgrFf[lFd]lFe]lFf]lFirFayFjrFD@$/Fg \\l7$F@Fen>Fbr,&-Fdr6+FfrFgrFiqF^]lF_]lF`]lFirFhyFesF@-Fdr6+FfrFgrFf[l Fd]lFe]lFf]lFirFixF`tFD@$3F^xFc\\lC%@$Ff\\l>Fbr,&-Fdr6+FfrFgrFiqF^]lF_ ]lF`]lFirFayFayF@-Fdr6+FfrFgrFf[lFd]lFe]lFf]lFirFhyFhyF@@$Fh]l>Fbr,&-F dr6+FfrFgrFiqF^]lF_]lF`]lFirFixFixF@-Fdr6+FfrFgrFf[lFd]lFe]lFf]lFirFay FayF@@$Fa^l>Fbr,&-Fdr6+FfrFgrFiqF^]lF_]lF`]lFirFhyFhyF@-Fdr6+FfrFgrFf[ lFd]lFe]lFf]lFirFixFixF@@$F^qC%@$Ff\\l>Fbr,&-Fdr6+FfrFgrFiqF^]lF_]lF`] lFirFjrFjrF@-Fdr6+FfrFgrFf[lFd]lFe]lFf]lFirFesFesF@@$Fh]l>Fbr,&-Fdr6+F frFgrFiqF^]lF_]lF`]lFirF`tF`tF@-Fdr6+FfrFgrFf[lFd]lFe]lFf]lFirFjrFjrF@ @$Fa^l>Fbr,&-Fdr6+FfrFgrFiqF^]lF_]lF`]lFirFesFesF@-Fdr6+FfrFgrFf[lFd]l Fe]lFf]lFirF`tF`tF@@$Fct@%FetC&>FiqFb[l>Fj[lF[\\l>8.-%,which_funcsG6%F VFFFI>Fbr*&&F^blFa]lF@-Fdr6+FfrFgrFiqF^]lF_]lF`]lFir,&F[sF@&F^blFRF@,& F[sF@&F^blFaqF@F@C&>FiqFg[l>Fj[lF_\\l>F^blF_bl>FbrFcbl*&FbrF@F>F@F;F;F ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1031 "#This function basic ally integrates funcno*h_func1*h_func2 over triangle.\n#it works by us ing coordinates relative to the subtriangle triangle, and \n#then usin g coord_transfer to determine which address of funcno in the \n#main t riangle to multiply by. Used in both int_vs_basis12 and int_vs_basis1 2_h...\n\nintegrate_basis12:=proc(T,m,triangle,add0,add1,add2,funcno,h _func1,h_func2)\n\nlocal fact,V0,Vm, Vother_e,Vother,V0sum1,Vothersum1 ,Vmsum1,i; \n\nfact:=1/(6*3^m);\n\nV0:=vertices_m(0);\nif m = 1 then V other:=[] else Vother_e:=vertices(m-1);Vother:=Vother_e[4..nops(Vother _e)]; fi;\nVm:=vertices_m(m);\n\nV0sum1:=add(T[coord_transfer(V0[i],tr iangle,add0,add1,add2),funcno]*T[V0[i],h_func1]*T[V0[i],h_func2],i=1.. nops(V0));\n\nVothersum1:=add(T[coord_transfer(Vother[i],triangle,add0 ,add1,add2),funcno] *T[Vother[i],h_func1]*T[Vother[i],h_func2],i=1..no ps(Vother));\n\nVmsum1:=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));\n\nRET URN(fact*(5*Vmsum1 + 2 * Vothersum1 + V0sum1));\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%2integrate_basis12Gf*6+%\"TG%\"mG%)triangleG%%ad d0G%%add1G%%add2G%'funcnoG%(h_func1G%(h_func2G6+%%factG%#V0G%#VmG%)Vot her_eG%'VotherG%'V0sum1G%+Vothersum1G%'Vmsum1G%\"iG6\"F:C*>8$,$*&\"\" \"F@)\"\"$9%!\"\"#F@\"\"'>8%-%+vertices_mG6#\"\"!@%/FCF@>8(7\"C$>8'-%) verticesG6#,&FCF@F@FD>FP&FT6#;\"\"%-%%nopsG6#FT>8&-FJ6#FC>8)-%$addG6$* (&9$6$-%/coord_transferG6'&FH6#8,9&9'9(9)9*F@&Ffo6$F[p9+F@&Ffo6$F[p9,F @/F]p;F@-Fin6#FH>8*-Fbo6$*(&Ffo6$-Fio6'&FPF\\pF^pF_pF`pFapFbpF@&Ffo6$F fqFepF@&Ffo6$FfqFhpF@/F]p;F@-Fin6#FP>8+-Fbo6$*(&Ffo6$-Fio6'&F\\oF\\pF^ pF_pF`pFapFbpF@&Ffo6$FhrFepF@&Ffo6$FhrFhpF@/F]p;F@-Fin6#F\\o-%'RETURNG 6#*&F=F@,(F`r\"\"&*&\"\"#F@F^qF@F@F`oF@F@F:F:F:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 2537 "#This function is called by int_vs_basis12. \n#it basically assumes basis1 and basis2 are neighbors,\n#and using t hat, determines which harmonic basis functions\n#you would want to mul tiply against when you integrate.\n#it uses the orientation of the ver tex, as well as the \n#f or g status. Also, it includes, as the last \+ element of the\n#output list, a multiplicative factor of 1 or -1, depe nding on whether\n#you want the positive or the negative of the result ant integral. This\n#sign depends on the g-vectors and which triangle you are in, relatively\n#speaking...\n\nwhich_funcs:=proc(S,basis1,ba sis2)\nlocal out;\n\n\nif last_digits(basis1[1]) = [0,1] then\n if S [basis1[1],1] = basis2[1] then out:=[1,5] fi;\n if S[basis1[1],2] = \+ basis2[1] then out:=[1,9] fi;\n if S[basis1[1],3] = basis2[1] then o ut:=[5,9] fi;\n if S[basis1[1],4] = basis2[1] then out:=[5,1] fi;\nf i;\nif last_digits(basis1[1]) = [0,2] then\n if S[basis1[1],1] = bas is2[1] then out:=[9,1] fi;\n if S[basis1[1],2] = basis2[1] then out: =[9,5] fi;\n if S[basis1[1],3] = basis2[1] then out:=[1,5] fi;\n i f S[basis1[1],4] = basis2[1] then out:=[1,9] fi;\nfi;\nif last_digits( basis1[1]) = [1,2] then\n if S[basis1[1],1] = basis2[1] then out:=[5 ,9] fi;\n if S[basis1[1],2] = basis2[1] then out:=[5,1] fi;\n if S [basis1[1],3] = basis2[1] then out:=[9,1] fi;\n if S[basis1[1],4] = \+ basis2[1] then out:=[9,5] fi;\nfi;\n\n#if they are both f's...\nif bas is1[2] = f and basis2[2] = f then out:=[op(out),1]; fi;\n\n#if 1 is an f and the other is a g...\nif basis1[2] = f and basis2[2] = g then\n \+ if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then ou t:=[op(out),1];\n else out:=[op(out),-1]; fi;\nfi; \n\nif basis1[2 ] = g and basis2[2] = f then\n if S[basis1[1],1] = basis2[1] or S[b asis1[1],2] = basis2[1] then out:=[op(out),1];\n else out:=[op(out) ,-1]; fi;\nfi; \n\n#if they're both g's (word!)...\nif basis1[2] = g a nd basis2[2] = g then\n if S[basis1[1],1] = basis2[1] or S[basis1[1] ,2] = basis2[1] then\n if S[basis2[1],1] = basis1[1] or S[basis2[1 ],2] = basis1[1] then out:=[op(out),1] fi;\n if S[basis2[1],3] = b asis1[1] or S[basis2[1],4] = basis1[1] then out:=[op(out),-1] fi;\n \+ fi; \n if S[basis1[1],3] = basis2[1] or S[basis1[1],4] = basis2[1] t hen\n if S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] \+ then out:=[op(out),-1] fi; \n if S[basis2[1],3] = basis1[1] or S[b asis2[1],4] = basis1[1] then out:=[op(out),1] fi;\n fi; \nfi;\n\n\n \+ \n\nif basis1[2] = g then out[1]:=out[1] + 2; fi;\nif basis2[2] = g th en out[2]:=out[2] + 2; fi;\nout;\nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%,which_funcsGf*6%%\"SG%'basis1G%'basis2G6#%$outG6\"F, C,@$/-%,last_digitsG6#&9%6#\"\"\"7$\"\"!F6C&@$/&9$6$F3F6&9&F5>8$7$F6\" \"&@$/&F=6$F3\"\"#F?>FB7$F6\"\"*@$/&F=6$F3\"\"$F?>FB7$FDFL@$/&F=6$F3\" \"%F?>FB7$FDF6@$/F07$F8FIC&@$F;>FB7$FLF6@$FF>FB7$FLFD@$FN>FBFC@$FU>FBF K@$/F07$F6FIC&@$F;>FBFS@$FF>FBFZ@$FN>FBF[o@$FU>FBF^o@$3/&F46#FI%\"fG/& F@FcpFdp>FB7$-%#opG6#FBF6@$3Fap/Ffp%\"gG@%5/&F=6$F?F6F3/&F=6$F?FIF3>FB Fhp>FB7$Fip!\"\"@$3/FbpF_qFep@%5F;FF>FBFhp>FBFjq@$3F^rF^qC$@$F`rC$@$Fa q>FBFhp@$5/&F=6$F?FQF3/&F=6$F?FXF3>FBFjq@$5FNFUC$@$Faq>FBFjq@$F[s>FBFh p@$F^r>&FBF5,&F\\tF6FIF6@$F^q>&FBFcp,&F`tF6FIF6FBF,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 877 "#we are going to make a bunch of f unctions to be used as f's and q's...\n#5 -> x\n#6 -> y\n#7 -> constan t = 1\n#8 -> harm(5,0,0)\n#9 -> harm(.1,0,0)\n#10 -> harm(100,0,0)\n#1 1 -> (harm(5,0,0))^2\n#12 -> quad-harm(-100,40,-130)\n#16 -> tri-harm( -100,40,-130)\n# function number 7 should be plugged in for funcno... \n \nmake_f_and_q:=proc(T,m,funcno)\n \n make_points(T,m);\n\n ha rmonic(T,m,1,1,1,funcno);\n\n harmonic(T,m,5,0,0,funcno + 1);\n\n \+ harmonic(T,m,.1,0,0,funcno + 2);\n\n harmonic(T,m,100,0,0,funcno + 3 );\n\n mult(T,m,funcno,funcno,funcno + 4);\n\n T[[0],funcno+4]:=0: \n T[[0],funcno+5]:=0:\n T[[0],funcno+6]:=0:\n T[[0],funcno+7]:= -100:\n T[[1],funcno+4]:=0:\n T[[1],funcno+5]:=0:\n T[[1],funcno +6]:=0:\n T[[1],funcno+7]:=40:\n T[[2],funcno+4]:=0:\n T[[2],fun cno+5]:=0:\n T[[2],funcno+6]:=0:\n T[[2],funcno+7]:=-130:\n n_ha rm(T,R,4,m,[0],[1],[2],funcno+5); \nend; \n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-make_f_and_qGf*6%%\"TG%\"mG%'funcnoG6\"F*F*C5-%,make _pointsG6$9$9%-%)harmonicG6(F/F0\"\"\"F4F49&-F26(F/F0\"\"&\"\"!F9,&F5F 4F4F4-F26(F/F0$F4!\"\"F9F9,&F5F4\"\"#F4-F26(F/F0\"$+\"F9F9,&F5F4\"\"$F 4-%%multG6'F/F0F5F5,&F5F4\"\"%F4>&F/6$7#F9FIF9>&F/6$FN,&F5F4F8F4F9>&F/ 6$FN,&\"\"'F4F5F4F9>&F/6$FN,&F5F4\"\"(F4!$+\">&F/6$7#F4FIF9>&F/6$F[oFR F9>&F/6$F[oFVF9>&F/6$F[oFen\"#S>&F/6$7#F@FIF9>&F/6$FioFRF9>&F/6$FioFVF 9>&F/6$FioFen!$I\"-%'n_harmG6*F/%\"RGFJF0FNF[oFioFRF*F*F*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}}{MARK "101 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }