{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 "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" -1 12 1 {CSTYLE "" -1 -1 "Tim es" 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 705 "# This is a < compa rison operator for words (where a word is a list of 0's, 1's,\n# and 2 's, i.e. [0,1,0,2]). This function is designed to be used to sort a li st\n# of words. Shorter words < longer words, and words of the same le ngth are \n# essentially compared like their numeric equivalents (i.e. [0,0,1] < [0,0,2])\nword_compare := proc(w1, w2)\n local i;\n\n \+ if nops(w1) < nops(w2) then \n RETURN(true);\n elif nops(w2 ) < nops(w1) then\n RETURN(false);\n else\n for i fro m 1 to nops(w1) do\n if w1[i] < w2[i] then\n \+ RETURN(true);\n elif w2[i] < w1[i] then\n R ETURN(false);\n fi;\n od;\n fi;\n RETURN(true) ;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-word_compareGf*6$%#w1G%# w2G6#%\"iG6\"F+C$@'2-%%nopsG6#9$-F06#9%-%'RETURNG6#%%trueG2F3F/-F76#%& falseG?(8$\"\"\"F@F/F9@&2&F26#F?&F5FDF62FEFCF;F6F+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1061 "# This takes a word list (which r epresents a subdivision of SG) and divides\n# everything one level fur ther.\nsubdivide_irreg := proc(word_list)\n local i, out;\n\n ou t := word_list;\n\n # \"base case\" - add \"first\" words if they a re not present\n if not member([0], word_list) then\n out := [op(out), [0]];\n fi;\n if not member([1], word_list) then\n \+ out := [op(out), [1]];\n fi;\n if not member([2], word_list ) then\n out := [op(out), [2]];\n fi;\n \n # take each word and add 3 new ones based on it\n for i from 1 to nops(word_li st) do\n if not member([op(word_list[i]),0], word_list) then\n \+ out := [op(out), [op(word_list[i]),0]];\n fi;\n \+ if not member([op(word_list[i]),1], word_list) then\n ou t := [op(out), [op(word_list[i]),1]];\n fi;\n if not mem ber([op(word_list[i]),2], word_list) then\n out := [op(out) , [op(word_list[i]),2]];\n fi;\n od;\n\n # make sure word list is sorted (see word_compare above)\n sort(out, word_compare); \nend;\n " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%0subdivide_irregGf*6#%*word_listG6$%\"iG%$outG6 \"F+C(>8%9$@$4-%'memberG6$7#\"\"!F/>F.7$-%#opG6#F.F5@$4-F36$7#\"\"\"F/ >F.7$F9F@@$4-F36$7#\"\"#F/>F.7$F9FH?(8$FAFA-%%nopsG6#F/%%trueGC%@$4-F3 6$7$-F:6#&F/6#FMF6F/>F.7$F9FW@$4-F36$7$FXFAF/>F.7$F9F\\o@$4-F36$7$FXFI F/>F.7$F9Fco-%%sortG6$F.%-word_compareGF+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 475 "# This function generates all of the vertices o f the subdivision represented\n# by word_list in one long list. \n# NO TE: these come out in nicely sorted order, but only if word_list is al ready\n# sorted.\nvertices_irreg := proc(word_list)\n local i, outl ist;\n\n outlist := [[0],[1],[2],[0,1],[0,2],[1,2]];\n\n for i f rom 1 to nops(word_list) do\n outlist := [op(outlist),[op(word_ list[i]),0,1],[op(word_list[i]),0,2],[op(word_list[i]),1,2]];\n od; \n\n outlist;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/vertices_ irregGf*6#%*word_listG6$%\"iG%(outlistG6\"F+C%>8%7(7#\"\"!7#\"\"\"7#\" \"#7$F1F37$F1F57$F3F5?(8$F3F3-%%nopsG6#9$%%trueG>F.7&-%#opG6#F.7%-FC6# &F>6#F:F1F37%FFF1F57%FFF3F5F.F+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 595 "# This is a boolean function that determines if word is of the format:\n# prefix + one occurrence of num1 + all num2\n# Th e first two cases are for words that are of the same length as prefix \n# or one more (they do not need to have the rest of the \"ideal form at\" to match).\ndesired_format := proc(word, prefix, num1, num2)\n \+ if nops(word) = nops(prefix) then\n evalb(word = [op(prefix)]) ;\n elif nops(word) = nops(prefix) + 1 then\n evalb(word = [ op(prefix), num1]);\n else\n evalb(word = [op(prefix), num1, seq(num2, i = 1..(nops(word) - nops(prefix) - 1))]);\n fi;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/desired_formatGf*6&%%wordG%'prefi xG%%num1G%%num2G6\"F+F+@'/-%%nopsG6#9$-F/6#9%-%&evalbG6#/F17#-%#opGF3/ F.,&F2\"\"\"F>F>-F66#/F17$F:9&-F66#/F17%F:FC-%$seqG6$9'/%\"iG;F>,(F.F> F2!\"\"F>FPF+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 478 "# Thi s function simply creates a gasket table based on a word list, and \n# initializes all neighbors to 0, so we can tell when a point hasn't ye t\n# been assigned its proper neighbor.\ngasket_irreg := proc(word_lis t)\n local i, vm, T;\n\n vm := vertices_irreg(word_list);\n T := table();\n\n for i from 1 to nops(vm) do\n T[vm[i], 0] : = vm[i];\n T[vm[i], 1] := 0;\n T[vm[i], 2] := 0;\n \+ T[vm[i], 3] := 0;\n T[vm[i], 4] := 0;\n od;\n\n T;\nend ;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-gasket_irregGf*6#%*word_listG6 %%\"iG%#vmG%\"TG6\"F,C&>8%-%/vertices_irregG6#9$>8&-%&tableGF,?(8$\"\" \"F:-%%nopsG6#F/%%trueGC'>&F56$&F/6#F9\"\"!FC>&F56$FCF:FE>&F56$FC\"\"# FE>&F56$FC\"\"$FE>&F56$FC\"\"%FEF5F,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8892 "# This function creates the entries 1, 2, 3, and 4 \+ in the table for every \n# vertex in the basis. Entries 1, 2, 3, and 4 represent a vertex's first, second, \n# third, and fourth neighbors, \+ respectively (numbered counterclockwise).\n# Notice that we start \"sh allow\" and go \"deep,\" coding each neighbor relationship\n# in both \+ directions as soon as we establish it.\nmake_neighbors_irreg := proc(T , word_list, vm)\n\nlocal i;\n\n# if word_list is empty, we have V1 an d we hardcode everything\nif nops(word_list) = 0 then\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], 1] := [1];\n T[[0,1], 2] := [1,2];\n T[[0,1], \+ 3] := [0,2];\n T[[0,1], 4] := [0];\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 T[[1,2], 1] := [2];\n T[[1,2], 2] := [0,2];\n T[[1,2] , 3] := [0,1];\n T[[1,2], 4] := [1];\n\nelse\n # hardcode the bo undary\n # Notice how we use longest_word, because we don't how dee p in SG\n # the neighbors will be.\n T[[0], 1] := [op(longest_wo rd(word_list, [], 0, 0)),0,1];\n T[[0], 2] := [op(T[[0], 1][1..(nop s(T[[0], 1]) - 2)]), 0, 2];\n T[[1], 1] := [op(longest_word(word_li st, [], 1, 1)),1,2];\n T[[1], 2] := [op(T[[1], 1][1..(nops(T[[1], 1 ]) - 2)]), 0, 1];\n T[[2], 1] := [op(longest_word(word_list, [], 2, 2)),0,2];\n T[[2], 2] := [op(T[[2], 1][1..(nops(T[[2], 1]) - 2)]), 1, 2];\n\n T[T[[0], 1], 4] := [0];\n T[T[[0], 2], 1] := [0];\n \+ T[T[[1], 1], 4] := [1];\n T[T[[1], 2], 1] := [1];\n T[T[[2], \+ 1], 4] := [2];\n T[T[[2], 2], 1] := [2]:\n\n #hardcode V1 \\ V0 \n # Structure of this is:\n # for each vertex in V1 \\ V0\n \+ # for each of its neighbors\n # if it hasn't been coded yet\n # find its neighbor\n # set relat ionship both ways\n\n if T[[0,1], 1] = 0 then\n T[[0,1], 1] \+ := [op(longest_word(word_list, [], 1, 0)),0,1];\n T[T[[0,1], 1] , 4] := [0,1];\n fi;\n if T[[0,1], 2] = 0 then\n if membe r([1], word_list) then \n T[[0,1], 2] := [op(longest_word(w ord_list, [], 1, 0)),0,2];\n T[T[[0,1], 2], 1] := [0,1];\n \+ else\n T[[0,1], 2] := [1,2];\n T[[1,2], 3 ] := [0,1];\n fi; \n fi;\n if T[[0,1], 3] = 0 then\n \+ if member([0], word_list) then \n T[[0,1], 3] := [op(lon gest_word(word_list, [], 0, 1)),1,2];\n T[T[[0,1], 3], 4] : = [0,1];\n else\n T[[0,1], 3] := [0,2];\n \+ T[[0,2], 2] := [0,1];\n fi; \n fi;\n if T[[0,1], 4] = 0 \+ then \n T[[0,1], 4] := [op(longest_word(word_list, [], 0, 1)),0 ,1];\n T[T[[0,1], 4], 1] := [0,1]; \n fi;\n \n if T[[0 ,2], 1] = 0 then \n T[[0,2], 1] := [op(longest_word(word_list, \+ [], 0, 2)),0,2];\n T[T[[0,2], 1], 4] := [0,2]; \n fi;\n i f T[[0,2], 2] = 0 then\n if member([0], word_list) then \n \+ T[[0,2], 2] := [op(longest_word(word_list, [], 0, 2)),1,2];\n \+ T[T[[0,2], 2], 1] := [0,2];\n else\n T[[0, 2], 2] := [0,1];\n T[[0,1], 3] := [0,2];\n fi; \n \+ fi;\n if T[[0,2], 3] = 0 then\n if member([2], word_list) t hen \n T[[0,2], 3] := [op(longest_word(word_list, [], 2, 0) ),0,1];\n T[T[[0,2], 3], 4] := [0,2];\n else\n \+ T[[0,2], 3] := [1,2];\n T[[1,2], 2] := [0,2];\n \+ fi; \n fi;\n if T[[0,2], 4] = 0 then \n T[[0,2], 4] := \+ [op(longest_word(word_list, [], 2, 0)),0,2];\n T[T[[0,2], 4], 1 ] := [0,2]; \n fi;\n \n if T[[1,2], 1] = 0 then \n T[[ 1,2], 1] := [op(longest_word(word_list, [], 2, 1)),1,2];\n T[T[ [1,2], 1], 4] := [1,2]; \n fi;\n if T[[1,2], 2] = 0 then\n \+ if member([2], word_list) then \n T[[1,2], 2] := [op(long est_word(word_list, [], 2, 1)),0,1];\n T[T[[1,2], 2], 1] := [1,2];\n else\n T[[1,2], 2] := [0,2];\n \+ T[[0,2], 3] := [1,2];\n fi; \n fi;\n if T[[1,2], 3] = 0 t hen\n if member([1], word_list) then \n T[[1,2], 3] \+ := [op(longest_word(word_list, [], 1, 2)),0,2];\n T[T[[1,2] , 3], 4] := [1,2];\n else\n T[[1,2], 3] := [0,1];\n \+ T[[0,1], 2] := [1,2];\n fi; \n fi;\n if T[[1,2 ], 4] = 0 then \n T[[1,2], 4] := [op(longest_word(word_list, [] , 1, 2)),1,2];\n T[T[[1,2], 4], 1] := [1,2]; \n fi;\n\n # end of hardcoding\n \n # Starting with V2 vertices\n for i fr om 7 to nops(vm) do\n \n # Here we have to handl e each ending separately \n # To really understand this, the on ly thing to do is work out\n # some examples to see the pattern s in the numbers.\n # It all falls out of the addressing system .\n if vm[i][(nops(vm[i])-1)..(nops(vm[i]))] = [0,1] then\n \+ \n if T[vm[i], 1] = 0 then\n T[vm[i] , 1] := [op(longest_word(word_list, vm[i][1..(nops(vm[i])-2)], 1, 0)), 0,1];\n T[T[vm[i], 1], 4] := vm[i];\n fi;\n \+ \n if T[vm[i], 2] = 0 then\n if m ember([op(vm[i][1..(nops(vm[i])-2)]), 1], word_list) then\n \+ T[vm[i], 2] := [op(longest_word(word_list, vm[i][1..(nops(vm[ i])-2)], 1, 0)),0,2];\n T[T[vm[i], 2], 1] := vm[i]; \n else\n T[vm[i], 2] := [op(vm[i][1 ..(nops(vm[i])-2)]),1,2];\n T[T[vm[i], 2], 3] := vm [i];\n fi;\n fi;\n if T[vm[i], 3] = 0 then\n if member([op(vm[i][1..(nops(vm[i])-2)]), 0 ], word_list) then\n T[vm[i], 3] := [op(longest_wor d(word_list, vm[i][1..(nops(vm[i])-2)], 0, 1)),1,2];\n \+ T[T[vm[i], 3], 4] := vm[i];\n else\n \+ T[vm[i], 3] := [op(vm[i][1..(nops(vm[i])-2)]),0,2];\n \+ T[T[vm[i], 3], 2] := vm[i];\n fi;\n \+ fi;\n if T[vm[i], 4] = 0 then\n T[vm[i], 4] \+ := [op(longest_word(word_list, vm[i][1..(nops(vm[i])-2)], 0, 1)),0,1]; \n T[T[vm[i], 4], 1] := vm[i];\n fi;\n \+ fi;\n if vm[i][(nops(vm[i])-1)..(nops(vm[i]))] = [0,2] then\n if T[vm[i], 1] = 0 then\n T[vm[i], 1] := [o p(longest_word(word_list, vm[i][1..(nops(vm[i])-2)], 0, 2)),0,2];\n \+ T[T[vm[i], 1], 4] := vm[i];\n fi;\n \+ if T[vm[i], 2] = 0 then\n if member([op(vm[i][1..(nops (vm[i])-2)]), 0], word_list) then\n T[vm[i], 2] := \+ [op(longest_word(word_list, vm[i][1..(nops(vm[i])-2)], 0, 2)),1,2];\n \+ T[T[vm[i], 2], 1] := vm[i];\n else\n T[vm[i], 2] := [op(vm[i][1..(nops(vm[i])-2)]),0,1] ;\n T[T[vm[i], 2], 3] := vm[i];\n fi ;\n fi;\n if T[vm[i], 3] = 0 then\n \+ if member([op(vm[i][1..(nops(vm[i])-2)]), 2], word_list) then\n \+ T[vm[i], 3] := [op(longest_word(word_list, vm[i][1..(n ops(vm[i])-2)], 2, 0)),0,1];\n T[T[vm[i], 3], 4] := vm[i];\n else\n T[vm[i], 3] := [op( vm[i][1..(nops(vm[i])-2)]),1,2];\n T[T[vm[i], 3], 2 ] := vm[i];\n fi;\n fi;\n if T[vm [i], 4] = 0 then\n T[vm[i], 4] := [op(longest_word(word _list, vm[i][1..(nops(vm[i])-2)], 2, 0)),0,2];\n T[T[vm [i], 4], 1] := vm[i];\n fi;\n fi;\n if vm[i][ (nops(vm[i])-1)..(nops(vm[i]))] = [1,2] then\n if T[vm[i], \+ 1] = 0 then\n T[vm[i], 1] := [op(longest_word(word_list , vm[i][1..(nops(vm[i])-2)], 2, 1)),1,2];\n T[T[vm[i], \+ 1], 4] := vm[i];\n fi;\n if T[vm[i], 2] = 0 then \n if member([op(vm[i][1..(nops(vm[i])-2)]), 2], word_l ist) then\n T[vm[i], 2] := [op(longest_word(word_li st, vm[i][1..(nops(vm[i])-2)], 2, 1)),0,1];\n T[T[v m[i], 2], 1] := vm[i];\n else\n T[vm [i], 2] := [op(vm[i][1..(nops(vm[i])-2)]),0,2];\n T [T[vm[i], 2], 3] := vm[i];\n fi;\n fi;\n \+ if T[vm[i], 3] = 0 then\n if member([op(vm[i][1 ..(nops(vm[i])-2)]), 1], word_list) then\n T[vm[i], 3] := [op(longest_word(word_list, vm[i][1..(nops(vm[i])-2)], 1, 2)),0 ,2];\n T[T[vm[i], 3], 4] := vm[i];\n \+ else\n T[vm[i], 3] := [op(vm[i][1..(nops(vm[i])-2) ]),0,1];\n T[T[vm[i], 3], 2] := vm[i];\n \+ fi;\n fi;\n if T[vm[i], 4] = 0 then\n \+ T[vm[i], 4] := [op(longest_word(word_list, vm[i][1..(nops(vm [i])-2)], 1, 2)),1,2];\n T[T[vm[i], 4], 1] := vm[i];\n \+ fi;\n fi;\n od;\n\nfi;\n\nT;\n \nend;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5make_neighbors_irregGf*6%%\"TG%*word_listG%#vmG 6#%\"iG6\"F,C$@%/-%%nopsG6#9%\"\"!C4>&9$6$7#F4\"\"\"7$F4F;>&F86$F:\"\" #7$F4F@>&F86$7#F;F;7$F;F@>&F86$FEF@F<>&F86$7#F@F;FA>&F86$FMF@FF>&F86$F &F86$F&F86$F<\"\"$FA>&F86$F<\"\"%F:>&F86$FAF;F:>&F86$FAF@F <>&F86$FAFZFF>&F86$FAFhnFM>&F86$FFF;FM>&F86$FFF@FA>&F86$FFFZF<>&F86$FF FhnFEC;>F77%-%#opG6#-%-longest_wordG6&F37\"F4F4F4F;>F>7%-Fep6#&F76#;F; ,&-F16#F7F;F@!\"\"F4F@>FC7%-Fep6#-Fhp6&F3FjpF;F;F;F@>FH7%-Fep6#&FC6#;F ;,&-F16#FCF;F@FeqF4F;>FK7%-Fep6#-Fhp6&F3FjpF@F@F4F@>FO7%-Fep6#&FK6#;F; ,&-F16#FKF;F@FeqF;F@>&F86$F7FhnF:>&F86$F>F;F:>&F86$FCFhnFE>&F86$FHF;FE >&F86$FKFhnFM>&F86$FOF;FM@$/FRF4C$>FR7%-Fep6#-Fhp6&F3FjpF;F4F4F;>&F86$ FRFhnF<@$/FUF4@%-%'memberG6$FEF3C$>FU7%F]uF4F@>&F86$FUF;FFUFF>F\\p F<@$/FXF4@%-Fhu6$F:F3C$>FX7%-Fep6#-Fhp6&F3FjpF4F;F;F@>&F86$FXFhnFF XFA>F]oF<@$/FfnF4C$>Ffn7%F[wF4F;>&F86$FfnF;F<@$/FjnF4C$>Fjn7%-Fep6#-Fh p6&F3FjpF4F@F4F@>&F86$FjnFhnFA@$/F]oF4@%FfvC$>F]o7%FbxF;F@>&F86$F]oF;F AC$>F]oF<>FXFA@$/F`oF4@%-Fhu6$FMF3C$>F`o7%-Fep6#-Fhp6&F3FjpF@F4F4F;>&F 86$F`oFhnFAC$>F`oFF>FioFA@$/FcoF4C$>Fco7%F]zF4F@>&F86$FcoF;FA@$/FfoF4C $>Ffo7%-Fep6#-Fhp6&F3FjpF@F;F;F@>&F86$FfoFhnFF@$/FioF4@%FhyC$>Fio7%Fd[ lF4F;>&F86$FioF;FFC$>FioFA>F`oFF@$/F\\pF4@%FguC$>F\\p7%-Fep6#-Fhp6&F3F jpF;F@F4F@>&F86$F\\pFhnFFC$>F\\pF<>FUFF@$/F_pF4C$>F_p7%F]]lF;F@>&F86$F _pF;FF?(8$\"\"(F;-F16#9&%%trueGC%@$/&&Fd^l6#F`^l6#;,&-F16#Fj^lF;F;FeqF __lFFd_l7%-Fep6#-Fhp6&F3&Fj^l6#;F;,&F__lF;F@FeqF ;F4F4F;>&F86$Fd_lFhnFj^l@$/&F86$Fj^lF@F4@%-Fhu6$7$-Fep6#F]`lF;F3C$>Ff` l7%Fi_lF4F@>&F86$Ff`lF;Fj^lC$>Ff`l7%F\\alF;F@>&F86$Ff`lFZFj^l@$/&F86$F j^lFZF4@%-Fhu6$7$F\\alF4F3C$>F\\bl7%-Fep6#-Fhp6&F3F]`lF4F;F;F@>&F86$F \\blFhnFj^lC$>F\\bl7%F\\alF4F@>&F86$F\\blF@Fj^l@$/&F86$Fj^lFhnF4C$>Fdc l7%FeblF4F;>&F86$FdclF;Fj^l@$/Fi^lFAC&@$Fc_lC$>Fd_l7%-Fep6#-Fhp6&F3F]` lF4F@F4F@>Fb`lFj^l@$Fe`l@%F_blC$>Ff`l7%FcdlF;F@>FbalFj^lC$>Ff`l7%F\\al F4F;>FhalFj^l@$F[bl@%-Fhu6$7$F\\alF@F3C$>F\\bl7%-Fep6#-Fhp6&F3F]`lF@F4 F4F;>FjblFj^lC$>F\\blFfal>F`clFj^l@$FcclC$>Fdcl7%FjelF4F@>FjclFj^l@$/F i^lFFC&@$Fc_lC$>Fd_l7%-Fep6#-Fhp6&F3F]`lF@F;F;F@>Fb`lFj^l@$Fe`l@%FdelC $>Ff`l7%F^glF4F;>FbalFj^lC$>Ff`lF^cl>FhalFj^l@$F[bl@%Fi`lC$>F\\bl7%-Fe p6#-Fhp6&F3F]`lF;F@F4F@>FjblFj^lC$>F\\blF`el>F`clFj^l@$FcclC$>Fdcl7%Fa hlF;F@>FjclFj^lF8F,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 458 "# This function simply finds the longest word in the list that has th e \"desired \n# format\" (see documentation for desired_format procedu re above)\n# CAUTION: This function assumes word_list is sorted!\nlong est_word := proc(word_list, prefix, num1, num2)\n local valid_words ;\n\n valid_words := select(desired_format, word_list, prefix, num1 , num2);\n\n if nops(valid_words) = 0 then\n [];\n else\n valid_words[nops(valid_words)];\n fi;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-longest_wordGf*6&%*word_listG%'prefixG%%num1G%% num2G6#%,valid_wordsG6\"F-C$>8$-%'selectG6'%/desired_formatG9$9%9&9'@% /-%%nopsG6#F0\"\"!7\"&F06#F;F-F-F-" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 908 "# This function finds the energy between two harmoni c basis elements. It is\n# based entirely on scaling theoretical value s\nenergy_irreg_h:=proc(S,basis1,basis2)\n\n local out,m1,m2;\n\n \+ # if they aren't neighbors or identical, we get nothing...\n if basi s1[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] <> bas is2[1] then RETURN(0) fi;\n\n # If they are the same, there are 2 tr iangles contributing, so scale by them\n # appropriately and add\n \+ if basis1[1] = basis2[1] then\n m1 := max(nops(basis1[1]), nops (S[basis1[1], 1])) - 1;\n m2 := max(nops(basis1[1]), nops(S[basi s1[1], 4])) - 1;\n out := 2 * ((5/3)^m1 + (5/3)^m2);\n # other wise there is only one triangle to deal with\n else\n m1 := ma x(nops(basis1[1]), nops(basis2[1])) - 1;\n out := -1 * (5/3)^m1; \n fi;\n \n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%/energy_irreg_hGf*6%%\"SG%'basis1G%'basis2G6%%$outG%#m1G%#m2G6 \"F.C&@$33330&9%6#\"\"\"&9&F80&9$6$F6F9F:0&F>6$F6\"\"#F:0&F>6$F6\"\"$F :0&F>6$F6\"\"%F:-%'RETURNG6#\"\"!@$/F6F:C%>8%,&-%$maxG6$-%%nopsG6#F6-F Z6#F=F9F9!\"\">8&,&-FW6$FY-FZ6#FIF9F9Fhn>8$,&)#\"\"&FGFTFC*&FCF9)FdoFj nF9F9@$555/F=F:/FAF:/FEF:/FIF:C$>FT,&-FW6$FY-FZ6#F:F9F9Fhn>Fao,$FcoFhn FaoF.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 598 "# This functi on creates the matrix of energies between all pairs of vertices\n# in \+ the basis. It is a sparse, symmetric matrix, and we take advantage of \+ this\n# in the coding.\ngen_energy_matrix_irreg_h := proc(S,vm)\n\nloc al i, j, en_matrix, n_verts, en;\n\nn_verts := nops(vm) - 3;\n\nen_mat rix := array(sparse, 1..n_verts, 1..n_verts);\n\nfor i from 1 to n_ver ts do\n for j from 1 to n_verts do\n en := energy_irreg_h(S, get_basis(vm,i), get_basis(vm,j));\n if en <> 0 then \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#>%:gen_energy_matrix_irreg_hGf*6$%\"SG%#vmG6'%\"iG%\"jG %*en_matrixG%(n_vertsG%#enG6\"F/C&>8',&-%%nopsG6#9%\"\"\"\"\"$!\"\">8& -%&arrayG6%%'sparseG;F8F2FA?(8$F8F8F2%%trueG?(8%F8F8F2FDC$>8(-%/energy _irreg_hG6%9$-%*get_basisG6$F7FC-FO6$F7FF@$0FI\"\"!C$>&F<6$FCFFFI>&F<6 $FFFCFIF " 0 "" {MPLTEXT 1 0 933 "# This \+ function finds the inner product (int(basis1*basis2)) between two harm onic\n# basis elements. It works the same way as energy_irreg_h, but w ith different\n# theoretical values.\ninner_prod_irreg_h := proc(S,bas is1,basis2)\n\n local out, m1, m2;\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] <> bas is2[1] and S[basis1[1],4] <> basis2[1] then RETURN(0) fi;\n\n\n # i f they are the same, we have 2 triangles, so scale by them and add\n \+ if basis1[1] = basis2[1] then\n m1 := max(nops(basis1[1]), no ps(S[basis1[1], 1])) - 1;\n m2 := max(nops(basis1[1]), nops(S[b asis1[1], 4])) - 1;\n \n out := 7/45 * ((1/3)^m1 + (1/3)^m2) ;\n # if neighbors, only one triangle\n else\n m1 := max( nops(basis1[1]), nops(basis2[1])) - 1;\n\n out := 4/45 / (3^m1) ; \n fi;\n \n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3inner_prod_irreg_hGf*6%%\"SG%'basis1G%'basis2G6%%$ou tG%#m1G%#m2G6\"F.C&@$33330&9%6#\"\"\"&9&F80&9$6$F6F9F:0&F>6$F6\"\"#F:0 &F>6$F6\"\"$F:0&F>6$F6\"\"%F:-%'RETURNG6#\"\"!@$/F6F:C%>8%,&-%$maxG6$- %%nopsG6#F6-FZ6#F=F9F9!\"\">8&,&-FW6$FY-FZ6#FIF9F9Fhn>8$,&)#F9FGFT#\" \"(\"#X*&FeoF9)FdoFjnF9F9@$555/F=F:/FAF:/FEF:/FIF:C$>FT,&-FW6$FY-FZ6#F :F9F9Fhn>Fao,$*&F9F9)FGFTFhn#FKFgoFaoF.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 505 "# generates the inner product matrix (just like g en_energy_matrix_irreg_h above).\ngen_iprod_matrix_irreg_h:=proc(S,vm) \n\nlocal i, j, i_matrix, n_verts, ip;\n\nn_verts := nops(vm) - 3;\n\n i_matrix := array(sparse, 1..n_verts, 1..n_verts);\n\nfor i from 1 to \+ n_verts do\n for j from 1 to i do\n ip := inner_prod_irreg_h (S, get_basis(vm,i), get_basis(vm,j));\n \n if ip <> 0 t hen \n i_matrix[i,j] := ip; \n i_matrix[j,i] := \+ ip; \n fi;\n od;\nod;\n\ni_matrix;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%9gen_iprod_matrix_irreg_hGf*6$%\"SG%#vmG6' %\"iG%\"jG%)i_matrixG%(n_vertsG%#ipG6\"F/C&>8',&-%%nopsG6#9%\"\"\"\"\" $!\"\">8&-%&arrayG6%%'sparseG;F8F2FA?(8$F8F8F2%%trueG?(8%F8F8FCFDC$>8( -%3inner_prod_irreg_hG6%9$-%*get_basisG6$F7FC-FO6$F7FF@$0FI\"\"!C$>&F< 6$FCFFFI>&F<6$FFFCFIF " 0 "" {MPLTEXT 1 0 5508 "# This function uses Simpson's rule to approximate the integra l of q * basis, \n# where q is stored in slot funct of T. \nint_vs_bas is_irreg_h:=proc(T,S,m_func,basis,funct,harm_funcno)\n\nlocal V0, Voth er1, Vother2, Vm1, Vm2, fact, triangle1, triangle2, V0sum1, V0sum2, Vo ther1sum, Vother2sum, Vm1sum, Vm2sum, i, j, Vother_e, m1, m2;\n\nfact \+ := 1/(6*3^m_func); #This is the factor for simpson's rule...\n\n# \"l evel\" to either side of basis vertex\nm1 := max(nops(basis[1]), nops( S[basis[1], 1])) - 1;\nm2 := max(nops(basis[1]), nops(S[basis[1], 4])) - 1;\n\n# set up lists of vertices ...\n# we will imagine the triangl es of the basis which contribute as entire SG's.\nV0 := vertices_m(0); \nif m_func - m1 <= 1 then \n Vother1 := []; \nelse \n Vother_e \+ := vertices(m_func - m1 - 1);\n Vother1 := Vother_e[4..nops(Vother_ e)]; \nfi;\nif m_func - m2 <= 1 then\n Vother2 := [];\nelse\n Vo ther_e := vertices(m_func - m2 - 1);\n Vother2 := Vother_e[4..nops( Vother_e)];\nfi;\nif m_func = m1 then\n Vm1 := [];\nelse\n Vm1 : = vertices_m(m_func - m1);\nfi;\nif m_func = m2 then\n Vm2 := [];\n else\n Vm2 := vertices_m(m_func - m2);\nfi;\n\n# establish the two \+ triangles contributing\ntriangle1 := get_triangle(S, basis, m1, 1);\nt riangle2 := get_triangle(S, basis, m2, 2);\n\n\n # Here we make use of the coord_transfer function, which\n # takes an address relativ e to some subtriangle of SG and tells us what\n # its absolute addr ess is in the real SG. Since we imagined the \n # contributing tria ngles as whole SG's, we can take each point in them\n # and use coo rd_transfer to get the q function value, and use the relative\n # a ddress to get the basis function value. It's by doing this that we\n \+ # avoid actually constructing all of the basis functions.\n\n # W e need to do cases based on the ending so that we orient everything\n \+ # correctly.\n if basis[1][nops(basis[1])-1..nops(basis[1])] = [ 0,1] then \n #let's calculate sums for triangle 1, shall we?\n \+ V0sum1 := add(T[coord_transfer(V0[i], triangle1, basis[1], S[ba sis[1], 1], S[basis[1], 2]), funct] * (T[V0[i], harm_funcno]), i = 1.. nops(V0));\n\n Vother1sum := add(T[coord_transfer(Vother1[i], t riangle1, basis[1], S[basis[1], 1], S[basis[1], 2]), funct] * (T[Vothe r1[i], harm_funcno]), i = 1..nops(Vother1)); \n \n Vm1sum \+ := add(T[coord_transfer(Vm1[i], triangle1, basis[1], S[basis[1], 1], S [basis[1], 2]), funct] * (T[Vm1[i], harm_funcno]), i = 1..nops(Vm1)); \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 + 1] ), i = 1..nops(V0));\n\n \+ Vother2sum := add(T[coord_transfer(Vother2[i], triangle2, S[basis [1], 4], basis[1], S[basis[1], 3]), funct] * (T[Vother2[i], harm_funcn o + 1]), i = 1..nops(Vother2)); \n\n Vm2sum := add(T[coo rd_transfer(Vm2[i], triangle2, S[basis[1], 4], basis[1], S[basis[1], 3 ]), funct] * (T[Vm2[i], harm_funcno + 1]), i = 1..nops(Vm2));\n\n f i;\n\n if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then\n \n #let's calculate sums for triangle 1, shall we?\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 + 2] ), i = 1. .nops(V0));\n\n Vother1sum := add(T[coord_transfer(Vother1[i], \+ triangle1, S[basis[1], 1], S[basis[1], 2], basis[1]), funct] * (T[Voth er1[i], harm_funcno + 2]), i = 1..nops(Vother1)); \n \n Vm 1sum := add(T[coord_transfer(Vm1[i], triangle1, S[basis[1], 1], S[basi s[1], 2], basis[1]), funct] * (T[Vm1[i], harm_funcno + 2]), i = 1..nop s(Vm1));\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 Vother2sum := add(T[coord_transfer(Vother2[i], triangle2, basi s[1], S[basis[1], 3], S[basis[1], 4]), funct] * (T[Vother2[i], harm_fu ncno]), i = 1..nops(Vother2)); \n\n Vm2sum := add(T[coor d_transfer(Vm2[i], triangle2, basis[1], S[basis[1], 3], S[basis[1], 4] ), funct] * (T[Vm2[i], harm_funcno]), i = 1..nops(Vm2));\n\n fi; \n \n if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] then\n \+ \n #let's calculate sums for triangle 1, shall we?\n V 0sum1 := add(T[coord_transfer(V0[i], triangle1, S[basis[1], 2], basis[ 1], S[basis[1], 1]), funct] * (T[V0[i], harm_funcno + 1]), i = 1..nops (V0));\n\n Vother1sum := add(T[coord_transfer(Vother1[i], trian gle1, S[basis[1], 2], basis[1], S[basis[1], 1]), funct] * (T[Vother1[i ], harm_funcno + 1]), i = 1..nops(Vother1)); \n \n Vm1sum \+ := add(T[coord_transfer(Vm1[i], triangle1, S[basis[1], 2], basis[1], S [basis[1], 1]), funct] * (T[Vm1[i], harm_funcno + 1]), i = 1..nops(Vm1 ));\n\n #and now on to triangle numero 2!\n V0sum2 := ad d(T[coord_transfer(V0[i], triangle2, S[basis[1], 3], S[basis[1], 4], b asis[1]), funct] * (T[V0[i], harm_funcno + 2]), i = 1..nops(V0));\n\n \+ Vother2sum := add(T[coord_transfer(Vother2[i], triangle2, S[bas is[1], 3], S[basis[1], 4], basis[1]), funct] * (T[Vother2[i], harm_fun cno + 2]), i = 1..nops(Vother2)); \n\n Vm2sum := add(T[c oord_transfer(Vm2[i], triangle2, S[basis[1], 3], S[basis[1], 4], basis [1]), funct] * (T[Vm2[i], harm_funcno + 2]), i = 1..nops(Vm2));\n\n \+ fi; \n\n\n#And then just return the formula for simpson's rule...\nRE TURN(fact * (5 * (Vm1sum + Vm2sum) + 2 * (Vother1sum + Vother2sum) + ( V0sum1 + V0sum2)));\n\n\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5 int_vs_basis_irreg_hGf*6(%\"TG%\"SG%'m_funcG%&basisG%&functG%,harm_fun cnoG65%#V0G%(Vother1G%(Vother2G%$Vm1G%$Vm2G%%factG%*triangle1G%*triang le2G%'V0sum1G%'V0sum2G%+Vother1sumG%+Vother2sumG%'Vm1sumG%'Vm2sumG%\"i G%\"jG%)Vother_eG%#m1G%#m2G6\"FAC.>8),$*&\"\"\"FG)\"\"$9&!\"\"#FG\"\"' >85,&-%$maxG6$-%%nopsG6#&9'6#FG-FU6#&9%6$FWFGFGFGFK>86,&-FR6$FT-FU6#&F gn6$FW\"\"%FGFGFK>8$-%+vertices_mG6#\"\"!@%1,&FJFGFOFKFG>8%7\"C$>84-%) verticesG6#,(FJFGFOFKFGFK>F]p&Fap6#;Fbo-FU6#Fap@%1,&FJFGFjnFKFG>8&F^pC $>Fap-Fcp6#,(FJFGFjnFKFGFK>F`qFgp@%/FJFO>8'F^p>Fjq-Ffo6#F[p@%/FJFjn>8( F^p>Far-Ffo6#F^q>8*-%-get_triangleG6&FgnFXFOFG>8+-Fhr6&FgnFXFjn\"\"#@$ /&FX6#F^s%\"fGC%@$/&FW6#;,&FTFGFGFKFT7$FhoFGC(>8,-%$addG6$*&&9$6$-%/co ord_transferG6'&Fdo6#82FfrFWFfn&Fgn6$FWF^s9(FG&Fdt6$Fit9)FG/F[u;FG-FU6 #Fdo>8.-F`t6$*&&Fdt6$-Fgt6'&F]pFjtFfrFWFfnF\\uF^uFG&Fdt6$F_vFauFG/F[u; FG-FU6#F]p>80-F`t6$*&&Fdt6$-Fgt6'&FjqFjtFfrFWFfnF\\uF^uFG&Fdt6$F_wFauF G/F[u;FG-FU6#Fjq>8--F`t6$*&&Fdt6$-Fgt6'FitF[sF`oFW&Fgn6$FWFIF^uFG&Fdt6 $Fit,&FauFGFGFGFGFbu>8/-F`t6$*&&Fdt6$-Fgt6'&F`qFjtF[sF`oFWF_xF^uFG&Fdt 6$F]yFcxFG/F[u;FG-FU6#F`q>81-F`t6$*&&Fdt6$-Fgt6'&FarFjtF[sF`oFWF_xF^uF G&Fdt6$F]zFcxFG/F[u;FG-FU6#Far@$/Fgs7$FhoF^sC(>F^t-F`t6$*&&Fdt6$-Fgt6' FitFfrFfnF\\uFWF^uFG&Fdt6$Fit,&FauFGF^sFGFGFbu>Fgu-F`t6$*&&Fdt6$-Fgt6' F_vFfrFfnF\\uFWF^uFG&Fdt6$F_vFb[lFGFbv>Fgv-F`t6$*&&Fdt6$-Fgt6'F_wFfrFf nF\\uFWF^uFG&Fdt6$F_wFb[lFGFbw>Fgw-F`t6$*&&Fdt6$-Fgt6'FitF[sFWF_xF`oF^ uFGF_uFGFbu>Fex-F`t6$*&&Fdt6$-Fgt6'F]yF[sFWF_xF`oF^uFG&Fdt6$F]yFauFGF` y>Fey-F`t6$*&&Fdt6$-Fgt6'F]zF[sFWF_xF`oF^uFG&Fdt6$F]zFauFGF`z@$/Fgs7$F GF^sC(>F^t-F`t6$*&&Fdt6$-Fgt6'FitFfrF\\uFWFfnF^uFGFaxFGFbu>Fgu-F`t6$*& &Fdt6$-Fgt6'F_vFfrF\\uFWFfnF^uFG&Fdt6$F_vFcxFGFbv>Fgv-F`t6$*&&Fdt6$-Fg t6'F_wFfrF\\uFWFfnF^uFG&Fdt6$F_wFcxFGFbw>Fgw-F`t6$*&&Fdt6$-Fgt6'FitF[s F_xF`oFWF^uFGF`[lFGFbu>Fex-F`t6$*&&Fdt6$-Fgt6'F]yF[sF_xF`oFWF^uFG&Fdt6 $F]yFb[lFGF`y>Fey-F`t6$*&&Fdt6$-Fgt6'F]zF[sF_xF`oFWF^uFG&Fdt6$F]zFb[lF GF`z-%'RETURNG6#*&FDFG,.Fgv\"\"&*&FdblFGFeyFGFG*&F^sFGFguFGFG*&F^sFGFe xFGFGF^tFGFgwFGFGFAFAFA" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 376 "# generates a vector of int(f * basis[i]), using simpson's rule\ngen_ f_vector_irreg_h := proc(T, S, m_func, vm, funcno, harm_funcno)\n\n \+ local i, f_vec, n_verts;\n\n n_verts := nops(vm) - 3;\n\n f_vec := array(1..n_verts);\n\n for i from 1 to n_verts do\n f_ve c[i] := int_vs_basis_irreg_h(T, S, m_func, get_basis(vm, i), funcno, h arm_funcno);\n od;\n\n f_vec;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5gen_f_vector_irreg_hGf*6(%\"TG%\"SG%'m_funcG%#vmG%'f uncnoG%,harm_funcnoG6%%\"iG%&f_vecG%(n_vertsG6\"F1C&>8&,&-%%nopsG6#9' \"\"\"\"\"$!\"\">8%-%&arrayG6#;F:F4?(8$F:F:F4%%trueG>&F>6#FD-%5int_vs_ basis_irreg_hG6(9$9%9&-%*get_basisG6$F9FD9(9)F>F1F1F1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 710 "# This function constructs an actu al function (which can be plotted, etc.) from\n# the list of coefficie nts C which results from the FEM.\nmake_pretty_picture_irreg_h:=proc(T ,S,vm,m_func,C,funcno)\n\n local j;\n\n # set the boundary\n \+ T[[0], funcno] := 0;\n T[[1], funcno] := 0; \n T[[2], funcno] := 0; \n\n # since they are harmonic, we can just set the values dire ctly ...\n for j from 1 to (nops(vm) - 3) do\n T[vm[j+3], fu ncno] := C[j];\n od;\n\n # ... then harmonically extend\n mak e_h_answer_irreg(T, S, m_func-1, [0], [0,1], [0,2], funcno);\n make _h_answer_irreg(T, S, m_func-1, [1], [1,2], [0,1], funcno);\n make_ h_answer_irreg(T, S, m_func-1, [2], [0,2], [1,2], funcno);\n\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%&9$6$7#\"\"!9)F6>& F36$7#\"\"\"F7F6>&F36$7#\"\"#F7F6?(8$F&F36$&FH6#,&FCF " 0 "" {MPLTEXT 1 0 1729 "# recursive helper to harmonically extend in irregular basis\nmake_h_answer_irreg := proc(T ,S,m_func,add0,add1,add2,funcno);\n\n if S[add1, 3] = add2 then\n \+ harm(T,m_func,add0,add1,add2,funcno);\n else\n if add1 [nops(add1)-1..nops(add1)] = [0,1] then\n make_h_answer_irr eg(T, S, m_func - 1, add0, [op(add1[1..nops(add1) - 2]), 0, 0, 1], [op (add1[1..nops(add1) - 2]), 0, 0, 2], funcno);\n make_h_answ er_irreg(T, S, 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_irreg(T, S, 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 add1[nops(add1)-1..nops(add1)] = [1,2] then\n \+ make_h_answer_irreg(T, S, 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_irreg(T, S, m_func - 1, add1, [op(add1[1..nops(add 1) - 2]), 1, 0, 2], [op(add1[1..nops(add1) - 2]), 1, 1, 2], funcno);\n make_h_answer_irreg(T, S, m_func - 1, add2, [op(add1[1..no ps(add1) - 2]), 1, 0, 1], [op(add1[1..nops(add1) - 2]), 1, 0, 2], func no);\n fi;\n if add1[nops(add1)-1..nops(add1)] = [0,2] t hen\n make_h_answer_irreg(T, S, 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_irreg(T, S, 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_answer_irreg(T, S, m_func - 1, add 2, [op(add1[1..nops(add1) - 2]), 2, 1, 2], [op(add1[1..nops(add1) - 2] ), 2, 0, 1], funcno);\n fi;\n\n fi;\nend;" }}{PARA 12 "" 1 " " {XPPMATH 20 "6#>%4make_h_answer_irregGf*6)%\"TG%\"SG%'m_funcG%%add0G %%add1G%%add2G%'funcnoG6\"F.F.@%/&9%6$9(\"\"$9)-%%harmG6(9$9&9'F4F69*C %@$/&F46#;,&-%%nopsG6#F4\"\"\"FH!\"\"FE7$\"\"!FHC%-F$6)F:F2,&F;FHFHFIF <7&-%#opG6#&F46#;FH,&FEFH\"\"#FIFKFKFH7&FQFKFKFXF=-F$6)F:F2FOF47&FQFKF HFXFPF=-F$6)F:F2FOF6FYFfnF=@$/FA7$FHFXC%-F$6)F:F2FOF<7&FQFHFHFX7&FQFHF KFHF=-F$6)F:F2FOF47&FQFHFKFXF_oF=-F$6)F:F2FOF6F`oFcoF=@$/FA7$FKFXC%-F$ 6)F:F2FOF<7&FQFXFKFX7&FQFXFHFXF=-F$6)F:F2FOF47&FQFXFKFHF\\pF=-F$6)F:F2 FOF6F]pF`pF=F.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3799 "# T his function uses Simpson's rule to approximate the integral of\n# q * basis1 * basis2, where q is a function stored in slot funct of T.\nin t_vs_basis12_irreg_h:=proc(T,S,m_func,basis1,basis2,funct,harm_funcno) \n\n local out, temp, fact, triangle1, triangle2, addresses1, addre sses2, j, which_harms, m1, m2;\n\n #if they aren't neighbors, we ge t nothing...\n if basis1[1] <> basis2[1] and S[basis1[1],1] <> basi s2[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 # first, let's get the triangles, the addre sses of the vertices of \n # those triangles, and the scaling f actors for them\n m1 := max(nops(basis1[1]), nops(S[basis1[1], \+ 1])) - 1;\n m2 := max(nops(basis1[1]), nops(S[basis1[1], 4])) - 1;\n\n triangle1 := get_triangle(S, basis1, m1, 1); \n \+ triangle2 := get_triangle(S, basis1, m2, 2);\n\n addresses1 : = get_addresses(S, basis1, m1, 1); \n addresses2 := get_addres ses(S, basis1, m2, 2);\n \n #Notice that if basis1 and bas is2 are in the same place, then we have to integrate over two triangle s\n #Another important thing is that the harmonic basis element s that you integrate against depend solely\n #upon the orientat ion of the basis elements, since they are in the same place.\n \+ #Notice that a lot of the work here is done in integrate_basis12, whic h is in splines.mws.\n if last_digits(basis1[1]) = [0,1] then\n out := (1 / 3^m1) * integrate_basis12(T, m_func - m1, tria ngle1, addresses1[1], addresses1[2], addresses1[3], funct, harm_funcno + 0, harm_funcno + 0) + (1 / 3^m2) * integrate_basis12(T, m_func - m2 , triangle2, addresses2[1], addresses2[2], addresses2[3], funct, harm_ funcno + 1, harm_funcno + 1); \n fi;\n\n if last_digits( basis1[1]) = [0,2] then\n out := (1 / 3^m1) * integrate_bas is12(T, m_func - m1, triangle1, addresses1[1], addresses1[2], addresse s1[3], funct, harm_funcno + 2, harm_funcno + 2) + (1 / 3^m2) * integra te_basis12(T, m_func - m2, triangle2, addresses2[1], addresses2[2], ad dresses2[3], funct, harm_funcno + 0, harm_funcno + 0); \n fi;\n \n if last_digits(basis1[1]) = [1,2] then\n out := ( 1 / 3^m1) * integrate_basis12(T, m_func - m1, triangle1, addresses1[1] , addresses1[2], addresses1[3], funct, harm_funcno + 1, harm_funcno + \+ 1) + (1 / 3^m2) * integrate_basis12(T, m_func - m2, triangle2, address es2[1], addresses2[2], addresses2[3], funct, harm_funcno + 2, harm_fun cno + 2); \n fi;\n \n #END OF CASE WHERE BASIS1 and BASI S2 are in same place\n\n # basis1 and basis2 are therefore neighbor s\n else\n m1 := max(nops(basis1[1]), nops(basis2[1])) - 1; \n\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 f irst triangle of basis1 (in the clockwise sense).\n triangl e1 := get_triangle(S, basis1, m1, 1); #first, let's get the triangle ...\n addresses1 := get_addresses(S, basis1, m1, 1); #and t he corresponding address...\n else\n #This case is i f the overlapping triangle is the second triangle of basis1 (in the cl ockwise sense).\n #Otherwise, it's just like before.\n \+ triangle1 := get_triangle(S, basis1, m1, 2);\n addre sses1 := get_addresses(S, basis1, m1, 2);\n fi;\n\n #Thi s function does most of the work for determining which basis functions to use.\n which_harms := which_funcs_h(S,basis1,basis2);\n \+ \n out := (1 / 3^m1) * integrate_basis12(T, m_func - m1, triang le1, addresses1[1], addresses1[2], addresses1[3], funct, harm_funcno + which_harms[1], harm_funcno + which_harms[2]);\n fi;\n\n out; \+ \nend;\n\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%7int_vs_basis12_irreg _hGf*6)%\"TG%\"SG%'m_funcG%'basis1G%'basis2G%&functG%,harm_funcnoG6-%$ outG%%tempG%%factG%*triangle1G%*triangle2G%+addresses1G%+addresses2G% \"jG%,which_harmsG%#m1G%#m2G6\"F:C&@$33330&9'6#\"\"\"&9(FD0&9%6$FBFEFF 0&FJ6$FB\"\"#FF0&FJ6$FB\"\"$FF0&FJ6$FB\"\"%FF-%'RETURNG6#\"\"!@$/FBFFC +>8-,&-%$maxG6$-%%nopsG6#FB-F`o6#FIFEFE!\"\">8.,&-F]o6$F_o-F`o6#FUFEFE Fdo>8'-%-get_triangleG6&FJFCFjnFE>8(-F_p6&FJFCFfoFO>8)-%.get_addresses GF`p>8*-FhpFdp@$/-%,last_digitsGFao7$FenFE>8$,&*&)FSFjnFdo-%2integrate _basis12G6+9$,&9&FEFjnFdoF]p&FfpFD&Ffp6#FO&Ffp6#FS9)9*FbrFEFE*&)FSFfoF do-Fgq6+Fiq,&F[rFEFfoFdoFbp&FjpFD&FjpF^r&FjpF`rFar,&FbrFEFEFEF[sFEFE@$ /F^q7$FenFO>Fbq,&*&FeqFdo-Fgq6+FiqFjqF]pF\\rF]rF_rFar,&FbrFEFOFEFdsFEF E*&FdrFdo-Fgq6+FiqFgrFbpFhrFirFjrFarFbrFbrFEFE@$/F^q7$FEFO>Fbq,&*&FeqF do-Fgq6+FiqFjqF]pF\\rF]rF_rFarF[sF[sFEFE*&FdrFdo-Fgq6+FiqFgrFbpFhrFirF jrFarFdsFdsFEFE@$555/FIFF/FMFF/FQFF/FUFFC&>Fjn,&-F]o6$F_o-F`o6#FFFEFEF do@%FftC$>F]pF^p>FfpFgpC$>F]p-F_p6&FJFCFjnFO>Ffp-FhpFiu>8,-%.which_fun cs_hG6%FJFCFG>Fbq*&FeqFdo-Fgq6+FiqFjqF]pF\\rF]rF_rFar,&FbrFE&F]vFDFE,& FbrFE&F]vF^rFEFEFbqF:F:F:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 589 "# generates the matrix of int(q * basis[i] * basis[j])\ngen_q_mat rix_irreg_h := proc(T,S,m_func,vm,funcno,harm_funcno)\n\n local i, \+ j, q_matrix, q_val, n_verts;\n\n n_verts := nops(vm) - 3;\n\n q_ matrix := array(sparse, 1..n_verts, 1..n_verts);\n\n\n for i from 1 to n_verts do\n for j from 1 to i do\n q_val := int _vs_basis12_irreg_h(T, S, m_func, get_basis(vm, i), get_basis(vm, j), \+ funcno, harm_funcno);\n if q_val <> 0 then \n \+ q_matrix[i,j] := q_val; \n q_matrix[j,i] := q_val;\n \+ fi;\n od;\n od;\n\nq_matrix;\nend;" }}{PARA 12 " " 1 "" {XPPMATH 20 "6#>%5gen_q_matrix_irreg_hGf*6(%\"TG%\"SG%'m_funcG% #vmG%'funcnoG%,harm_funcnoG6'%\"iG%\"jG%)q_matrixG%&q_valG%(n_vertsG6 \"F3C&>8(,&-%%nopsG6#9'\"\"\"\"\"$!\"\">8&-%&arrayG6%%'sparseG;F8'-%7int_vs_basis12_irreg_hG6)9$9%9&-%* get_basisG6$F;FG-FU6$F;FJ9(9)@$0FM\"\"!C$>&F@6$FGFJFM>&F@6$FJFGFMF@F3F 3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1157 "# This function si mply encapsulates all of the FEM functions together\n# CAUTION: This f unction uses T, assuming it has been created properly.\nharm_FEM := pr oc(W, level, baseno)\n local S, vm, EM, QM, EM_QM, f_v, C;\n \n \+ print(\"making harms\");\n calc_harms(T, level, 20):\n print( \"making basis gasket\");\n S := gasket_irreg(W):\n print(\"maki ng basis vertices\");\n vm := vertices_irreg(W):\n print(\"# of \+ basis vertices:\", nops(vm));\n print(\"making basis neighbors\"); \n make_neighbors_irreg(S,W,vm):\n print(\"generating energy mat rix\");\n EM := gen_energy_matrix_irreg_h(S,vm):\n print(\"gener ating q matrix\");\n QM := gen_q_matrix_irreg_h(T,S,level,vm,14,20) :\n print(\"adding matrices\");\n EM_QM := evalm(EM + QM):\n \+ print(\"calculating f\");\n mult(T,level,14,8,16): \n add_func_ lin(T,level,9,16,-1,1,11):\n print(\"generating f vector\");\n f _v := gen_f_vector_irreg_h(T,S,level,vm,11,20):\n print(\"solving e quation\");\n C := evalf(linsolve(EM_QM,f_v),20):\n print(\"maki ng points\");\n make_points(T,level):\n print(\"making pretty pi cture\");\n make_pretty_picture_irreg_h(T,S,vm,level,C,baseno+15): \n \nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%)harm_FEMGf*6%%\"WG%&le velG%'basenoG6)%\"SG%#vmG%#EMG%#QMG%&EM_QMG%$f_vG%\"CG6\"F2C<-%&printG 6#Q-making~harmsF2-%+calc_harmsG6%%\"TG9%\"#?-F56#Q4making~basis~gaske tF2>8$-%-gasket_irregG6#9$-F56#Q6making~basis~verticesF2>8%-%/vertices _irregGFE-F56$Q5#~of~basis~vertices:F2-%%nopsG6#FK-F56#Q7making~basis~ neighborsF2-%5make_neighbors_irregG6%FBFFFK-F56#Q9generating~energy~ma trixF2>8&-%:gen_energy_matrix_irreg_hG6$FBFK-F56#Q4generating~q~matrix F2>8'-%5gen_q_matrix_irreg_hG6(F;FBF 8(-%&evalmG6#,&F`o\"\"\"FhnF^p-F56#Q.calculating~fF2-%%multG6'F;F8)-%5gen_f_vector_irreg_hG6(F;FBF8*-%&evalfG6$-%)linsolveG6$FioFaqF=-F56#Q.making~pointsF2-%,ma ke_pointsG6$F;F<-F56#Q6making~pretty~pictureF2-% " 0 "" {MPLTEXT 1 0 4549 "#This function returns the inner product of two bas is vectors in the biharmonic basis.\n#It returns theoretical values fr om the splines on fractals paper (scaled down by the appropriate\n#fac tor of 3 because we are dealing with smaller triangles...). Notice th ere are several\n#cases to deal with, and one must also be careful to \+ use the fact that we have defined the \n#normal derivative to be negat ive in the 2nd triangle...\n\n# We also must be careful to scale value s for g-type basis functions, something\n# which only makes a differen ce on uneven grids. \n\ninner_prod_irreg := proc(S, basis_1, basis_2) \n\n local fact_ff, fact_gg, fact_fg, out, temp, basis1, basis2, m1 , m2;\n \n #We do this so that we can change basis1 and basis2\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 \n RETURN( 0);\n 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 basi s1[1] = [2] or basis2[1] = [0] or basis2[1] = [1] or basis2[1] = [2] t hen\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 te mp := basis2;\n basis2 := basis1;\n basis1 := te mp;\n fi;\n \n if basis1[1] = basis2[1] then\n \+ out := 206/37665 * (3/25) ^ (nops(S[basis1[1], 1]) - 1); \n \+ else\n fact_gg := (3/25) ^ (nops(basis1[1]) - 1); \+ \n \n if basis1[2] = g then \n \+ if S[basis1[1],1] = basis2[1] then \n ou t := 83/37665 * fact_gg; \n else \n \+ out := -83/37665 * fact_gg; \n fi; \n else \+ \n out := -61/5580 * (1/5) ^ (nops(basis1[1]) - 1); \n fi;\n fi;\n #OTHERWISE, if neither are on the boundaries...\n else\n\n #This is the case of the basis elt s being in the same place...\n if basis1[1] = basis2[1] then\n \+ if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and \+ basis2[2] = f) then \n out := 47/1395 * ((1/5) ^ (max(n ops(basis1[1]), nops(S[basis1[1], 4])) - 1) -\n \+ (1/5) ^ (max(nops(basis1[1]), nops(S[basis1[1], 1])) - 1)); \n elif basis1[2] = f and basis2[2] = f then \n \+ out := 190/837 * ((1/3) ^ (max(nops(basis1[1]), nops(S[basis1[1], 1])) - 1) +\n (1/3) ^ (max(nops(basi s1[1]), nops(S[basis1[1], 4])) - 1)); \n else \n \+ out := 206/37665 * ((3/25) ^ (max(nops(basis1[1]), nops(S[basis1[ 1], 1])) - 1) +\n (3/25) ^ (max(nop s(basis1[1]), nops(S[basis1[1], 4])) - 1)); \n fi;\n \+ #This is the case of the basis elts being neighbors\n else\n \+ fact_ff := (1/3) ^ (max(nops(basis1[1]), nops(basis2[1])) - \+ 1);\n fact_gg := (3/25) ^ (max(nops(basis1[1]), nops(basis2 [1])) - 1);\n fact_fg := (1/5) ^ (max(nops(basis1[1]), nops (basis2[1])) - 1);\n\n if basis1[2] = f and basis2[2] = f t hen \n out := 89/1674 * fact_ff; \n elif bas is1[2] = g and basis2[2] = g then \n if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then\n i f S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then \n \+ out := 83/37665 * fact_gg; \n e lse \n out := -83/37665 * fact_gg; \n \+ fi;\n else \n if S[basis2 [1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then \n \+ out := -83/37665 * fact_gg; \n else \n \+ out := 83/37665 * fact_gg; \n \+ fi;\n fi;\n elif (basis1[2] = f and basis2[2 ] = g) then \n if S[basis2[1],1] = basis1[1] or S[basis 2[1],2] = basis1[1] then \n out := -61/5580 * fact_ fg; \n else \n out := 61/5580 * fac t_fg; \n fi;\n else\n if S[ba sis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then \n \+ out := -61/5580 * fact_fg; \n else \n \+ out := 61/5580 * fact_fg; \n fi;\n \+ fi; \n fi;\n fi;\n\n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%1inner_prod_irregGf*6%%\"SG%(basis_1G%(basis_2G6+% (fact_ffG%(fact_ggG%(fact_fgG%$outG%%tempG%'basis1G%'basis2G%#m1G%#m2G 6\"F4C'>8)9%>8*9&@$33330&F76#\"\"\"&F:FC0&9$6$FBFDFE0&FH6$FB\"\"#FE0&F H6$FB\"\"$FE0&FH6$FB\"\"%FE-%'RETURNG6#\"\"!@%55555/FB7#FY/FB7#FD/FB7# FM/FEF[o/FEF]o/FEF_oC$@$/-%%nopsG6#FBFDC%>8(F:>F:F7>F7F[p@%/FBFE>8',$) #FQ\"#D,&-Fgo6#FGFDFD!\"\"#\"$1#\"&lw$C$>8%)Fdp,&FfoFDFDFip@%3/&F76#FM %\"gG/&F:FfqFgq@%/FGFE>Fap,$F_q#\"#$)F\\q>Fap,$F_q#!#$)F\\q>Fap,$)#FD \"\"&Faq#!#h\"%!e&@%F_p@'53/Feq%\"fGFhq3Fdq/FiqFas>Fap,&)Fgr,&-%$maxG6 $Ffo-Fgo6#FSFDFDFip#\"#Z\"%&R\"*&#F^tF_tFD)Fgr,&-Fis6$FfoFgpFDFDFipFDF ip3F`sFcs>Fap,&)#FDFQFct#\"$!>\"$P)*&F[uFD)FjtFgsFDFD>Fap,&)FdpFctFjp* &FjpFD)FdpFgsFDFDC&>8$)Fjt,&-Fis6$Ffo-Fgo6#FEFDFDFip>F_q)FdpFiu>8&)Fgr Fiu@)Fft>Fap,$Fgu#\"#*)\"%u;Fcq@%5F[r/FKFE@%5/&FH6$FEFDFB/&FH6$FEFMFB> FapF]r>FapFar@%F]w>FapFar>FapF]rF_s@%F]w>Fap,$FavFir>Fap,$Fav#\"#hF[s@ %Fjv>FapF[x>FapF]xFapF4F4F4" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 528 "# Same as gen_iprod_matrix_irreg_h but for biharmonic splines\nge n_iprod_matrix_irreg := proc(S, vm)\n\n local i, j, n_verts, i_matr ix, ip;\n \n n_verts := 2 * nops(vm) - 3;\n i_matrix := array(sp arse, 1..n_verts, 1..n_verts);\n\n for i from 1 to n_verts do\n \+ for j from 1 to i do\n ip := inner_prod_irreg(S, get_ba sis(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;\n od;\n\n i_matrix;\nend; \n" }} {PARA 12 "" 1 "" {XPPMATH 20 "6#>%7gen_iprod_matrix_irregGf*6$%\"SG%#v mG6'%\"iG%\"jG%(n_vertsG%)i_matrixG%#ipG6\"F/C&>8&,&-%%nopsG6#9%\"\"# \"\"$!\"\">8'-%&arrayG6%%'sparseG;\"\"\"F2FA?(8$FBFBF2%%trueG?(8%FBFBF DFEC$>8(-%1inner_prod_irregG6%9$-%*get_basisG6$F7FD-FP6$F7FG@$0FJ\"\"! C$>&F<6$FDFGFJ>&F<6$FGFDFJF " 0 "" {MPLTEXT 1 0 4185 "#This function returns the energy of two basis vect ors in the biharmonic basis.\n#It returns 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#normal derivative to be negative in \+ the 2nd triangle...\n\n\n# We also must be careful to scale the functi on values of the g-type basis\n# functions - this only makes a differe nce on irregular grids\n\nenergy_irreg := proc(S, basis_1, basis_2)\n \n local fact_ff, fact_gg, out, temp, basis1, basis2, m1, m2;\n \+ \n #We do this so that we can change basis1 and basis2\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] <> basis 2[1] and S[basis1[1],4] <> basis2[1] then \n RETURN(0);\n 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 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 := basis 2;\n basis2 := basis1;\n basis1 := temp;\n \+ fi;\n \n if basis1[1] = basis2[1] then\n out := 5/27 * (3/5) ^ (nops(S[basis1[1], 1]) - 1); \n else\n \+ fact_gg := (3/5) ^ (nops(basis1[1]) - 1); \+ \n \n if basis1[2] = g and basis2[2] = g then \n \+ if S[basis1[1],1] = basis2[1] then \n o ut := -1/108 * fact_gg; \n else \n o ut := 1/108 * fact_gg; \n fi; \n else \n \+ out := 7/36; \n fi;\n fi;\n #OTHE RWISE, if neither are on the boundaries...\n else\n\n #This \+ is the case of the basis elts being in the same place...\n if b asis1[1] = basis2[1] then\n if (basis1[2] = f and basis2[2] = g) or (basis1[2] = g and basis2[2] = f) then \n out \+ := 0; \n elif basis1[2] = f and basis2[2] = f then \n \+ out := 19/6 * ((5/3) ^ (max(nops(basis1[1]), nops(S[basis1[ 1], 1])) - 1) +\n (5/3) ^ (max(nops(basi s1[1]), nops(S[basis1[1], 4])) - 1)); \n else \n \+ out := 5/27 * ((3/5) ^ (max(nops(basis1[1]), nops(S[basis1[1], 1] )) - 1) +\n (3/5) ^ (max(nops(basis1[1]) , nops(S[basis1[1], 4])) - 1)); \n fi;\n #This is th e case of the basis elts being neighbors\n else\n fa ct_ff := (5/3) ^ (max(nops(basis1[1]), nops(basis2[1])) - 1);\n \+ fact_gg := (3/5) ^ (max(nops(basis1[1]), nops(basis2[1])) - 1);\n \n if basis1[2] = f and basis2[2] = f then \n \+ out := -19/12 * fact_ff; \n elif basis1[2] = g and basis2 [2] = g then \n if S[basis1[1],1] = basis2[1] or S[bas is1[1],2] = basis2[1] then\n if S[basis2[1],1] = ba sis1[1] or S[basis2[1],2] = basis1[1] then \n o ut := -1/108 * fact_gg; \n else \n \+ out := 1/108 * fact_gg; \n fi;\n \+ else \n if S[basis2[1],1] = basis1[1] or S[bas is2[1],2] = basis1[1] then \n out := 1/108 * fa ct_gg; \n else \n out := -1/ 108 * fact_gg; \n fi;\n fi;\n \+ elif (basis1[2] = f and basis2[2] = g) then \n if \+ S[basis2[1],1] = basis1[1] or S[basis2[1],2] = basis1[1] then \n \+ out := 7/36; \n else \n \+ out := -7/36; \n fi;\n else\n \+ if S[basis1[1],1] = basis2[1] or S[basis1[1],2] = basis2[1] then \+ \n out := 7/36; \n else \n \+ out := -7/36; \n fi;\n fi; \n \+ fi;\n fi;\n\n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%-energy_irregGf*6%%\"SG%(basis_1G%(basis_2G6*%(fact_ffG%(fact_ ggG%$outG%%tempG%'basis1G%'basis2G%#m1G%#m2G6\"F3C'>8(9%>8)9&@$33330&F 66#\"\"\"&F9FB0&9$6$FAFCFD0&FG6$FA\"\"#FD0&FG6$FA\"\"$FD0&FG6$FA\"\"%F D-%'RETURNG6#\"\"!@%55555/FA7#FX/FA7#FC/FA7#FL/FDFjn/FDF\\o/FDF^oC$@$/ -%%nopsG6#FAFCC%>8'F9>F9F6>F6Fjo@%/FAFD>8&,$)#FP\"\"&,&-Ffo6#FFFCFC!\" \"#Fdp\"#FC$>8%)Fcp,&FeoFCFCFhp@%3/&F66#FL%\"gG/&F9FdqFeq@%/FFFD>F`p,$ F]q#Fhp\"$3\">F`p,$F]q#FCF]r>F`p#\"\"(\"#O@%F^p@'53/Fcq%\"fGFfq3Fbq/Fg qFjr>F`pFX3FirF\\s>F`p,&)#FdpFP,&-%$maxG6$FeoFfpFCFCFhp#\"#>\"\"'*&Fgs FC)Fbs,&-Fes6$Feo-Ffo6#FRFCFCFhpFCFC>F`p,&)FcpFcsFip*&FipFC)FcpF\\tFCF CC%>8$)Fbs,&-Fes6$Feo-Ffo6#FDFCFCFhp>F]q)FcpFjt@)F^s>F`p,$Fht#!#>\"#7F aq@%5Fiq/FJFD@%5/&FG6$FDFCFA/&FG6$FDFLFA>F`pF[r>F`pF_r@%F[v>F`pF_r>F`p F[rFhr@%F[v>F`pFbr>F`p#!\"(Fdr@%Fhu>F`pFbr>F`pFjvF`pF3F3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 658 "#This procedure returns the matrix E(fi,fj) where fi, fj are spline basis functions.\n#it stores values \+ in sparse matrix form, and also takes advantage of the fact that the\n #matrix is symmetric...\n\ngen_energy_matrix_irreg := proc(S, vm)\n\n \+ local i, j, n_verts, en_matrix, en;\n\n n_verts := 2 * nops(vm) \+ - 3;\n en_matrix := array(sparse, 1..n_verts, 1..n_verts);\n\n f or i from 1 to n_verts do\n for j from 1 to i do\n e n := energy_irreg(S, get_basis(vm,i), get_basis(vm,j));\n i f en<>0 then\n en_matrix[i,j] := en;\n e n_matrix[j,i] := en;\n fi;\n od;\n od;\n\n en_ matrix;\nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%8gen_energy_matr ix_irregGf*6$%\"SG%#vmG6'%\"iG%\"jG%(n_vertsG%*en_matrixG%#enG6\"F/C&> 8&,&-%%nopsG6#9%\"\"#\"\"$!\"\">8'-%&arrayG6%%'sparseG;\"\"\"F2FA?(8$F BFBF2%%trueG?(8%FBFBFDFEC$>8(-%-energy_irregG6%9$-%*get_basisG6$F7FD-F P6$F7FG@$0FJ\"\"!C$>&F<6$FDFGFJ>&F<6$FGFDFJF " 0 "" {MPLTEXT 1 0 8418 "# integrates (function * basis) using S impson's rule ... similar to harmonic\n# version, but of course more c omplicated\nint_vs_basis_irreg := proc(T, S, m_func, basis, funct, har m_funcno)\n\n local m1, m2, V0, Vother1, Vother2, Vm1, Vm2, fact, t riangle1, triangle2, V0sum1, V0sum2, Vothersum1, Vothersum2, Vmsum1, V msum2, i, j, Vother_e, blah;\n\n fact := 1 / (6 * 3^m_func); #facto r for simpson's rule...\n\n m1 := max(nops(basis[1]), nops(S[basis[ 1], 1])) - 1;\n m2 := max(nops(basis[1]), nops(S[basis[1], 4])) - 1 ;\n\n V0 := vertices_m(0);\n if m_func - m1 <= 1 then \n \+ Vother1 := []; \n else \n Vother_e := vertices(m_func - m1 - 1);\n Vother1 := Vother_e[4..nops(Vother_e)]; \n fi;\n i f m_func - m2 <= 1 then\n Vother2 := [];\n else\n Vot her_e := vertices(m_func - m2 - 1);\n Vother2 := Vother_e[4..no ps(Vother_e)];\n fi;\n if m_func = m1 then\n Vm1 := [];\n else\n Vm1 := vertices_m(m_func - m1);\n fi;\n if m_f unc = m2 then\n Vm2 := [];\n else\n Vm2 := vertices_m (m_func - m2);\n fi;\n \n \n if nops(basis[1]) = 1 then #t his is if the basis is an element of V0...\n #Here, we know the basis element is a g type vector. Then, we use the fact that\n \+ #the triangle in question is just a sequence of 0's or 1's. Then, w e multiply by the appropriate\n #biharmonic basis function (ie, harm_funcno+3 or some such number), and sum...\n \n\n t riangle1 := [seq(basis[1][1], j = 1..m1)];\n \n if basis [1] = [0] then\n V0sum1 := add(T[coord_transfer(V0[i], tria ngle1, basis[1], S[basis[1], 1], S[basis[1], 2]), funct] * (T[V0[i], h arm_funcno + 3]), i = 1..nops(V0));\n\n Vothersum1 := add(T [coord_transfer(Vother1[i], triangle1, basis[1], S[basis[1], 1], S[bas is[1], 2]), funct] * (T[Vother1[i], harm_funcno + 3]), i = 1..nops(Vot her1)); \n \n Vmsum1 := add(T[coord_transfer(Vm1[i], t riangle1, basis[1], S[basis[1], 1], S[basis[1], 2]), funct] * (T[Vm1[i ], harm_funcno + 3]), i = 1..nops(Vm1));\n fi;\n if basi s[1] = [1] then\n V0sum1 := add(T[coord_transfer(V0[i], tri angle1, 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_transfer(Vother1[i], triangle1, S[basis[1], 2], basis[1], S[ba sis[1], 1]), funct] * (T[Vother1[i], harm_funcno + 7]), i = 1..nops(Vo ther1)); \n \n Vmsum1 := add(T[coord_transfer(Vm1[i], \+ triangle1, S[basis[1], 2], basis[1], S[basis[1], 1]), funct] * (T[Vm1[ i], harm_funcno + 7]), i = 1..nops(Vm1));\n fi;\n if bas is[1] = [2] then\n V0sum1 := add(T[coord_transfer(V0[i], tr iangle1, S[basis[1], 1], S[basis[1], 2], basis[1]), funct] * (T[V0[i], harm_funcno + 11]), i = 1..nops(V0));\n\n Vothersum1 := ad d(T[coord_transfer(Vother1[i], triangle1, S[basis[1], 1], S[basis[1], \+ 2], basis[1]), funct] * (T[Vother1[i], harm_funcno + 11]), i = 1..nops (Vother1)); \n \n Vmsum1 := add(T[coord_transfer(Vm1[i ], triangle1, S[basis[1], 1], S[basis[1], 2], basis[1]), funct] * (T[V m1[i], harm_funcno + 11]), i = 1..nops(Vm1));\n fi;\n \n #normal derivative scaling\n V0sum1 := V0sum1 * (3/5) ^ m1;\n Vothersum1 := Vothersum1 * (3/5) ^ m1;\n Vmsum1 : = Vmsum1 * (3/5) ^ m1;\n \n \n RETURN(fact * (5 *(Vmsum1) + 2*(Vothersum1) + (V0sum1))); #And then return the simpson 's rule formula...\n\n\n else #So we will come here if the basis el ement is not on the boundary (could be either f or g...)\n\n tr iangle1 := get_triangle(S, basis, m1, 1);\n triangle2 := get_tr iangle(S, basis, m2, 2);\n\n\n # set flag (determines which har m_func values to use later)\n if basis[2]=f then\n b lah := 0;\n else\n blah := 2;\n fi;\n\n \+ if basis[1][nops(basis[1])-1..nops(basis[1])] = [0,1] then \n \+ \n #let's calculate sums for triangle 1, shall we?\n \+ V0sum1 := add(T[coord_transfer(V0[i], triangle1, basis[1], S[basi s[1], 1], S[basis[1], 2]), funct] * (T[V0[i], harm_funcno + 1 + blah]) , i = 1..nops(V0));\n\n Vothersum1 := add(T[coord_transfer( Vother1[i], triangle1, basis[1], S[basis[1], 1], S[basis[1], 2]), func t] * (T[Vother1[i], harm_funcno + 1 + blah]), i = 1..nops(Vother1)); \+ \n \n Vmsum1 := add(T[coord_transfer(Vm1[i], triangle1 , basis[1], S[basis[1], 1], S[basis[1], 2]), funct] * (T[Vm1[i], harm_ funcno + 1 + blah]), i = 1..nops(Vm1));\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 + blah]), i = 1..nops(V0));\n\n Vothersu m2 := add(T[coord_transfer(Vother2[i], triangle2, S[basis[1], 4], basi s[1], S[basis[1], 3]), funct] * (T[Vother2[i], harm_funcno + 5 + blah] ), i = 1..nops(Vother2)); \n\n Vmsum2 := add(T[coord _transfer(Vm2[i], triangle2, S[basis[1], 4], basis[1], S[basis[1], 3]) , funct] * (T[Vm2[i], harm_funcno + 5 + blah]), i = 1..nops(Vm2));\n\n elif basis[1][nops(basis[1])-1..nops(basis[1])] = [0,2] then\n \n #let's calculate sums for triangle 1, shall we?\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 + 9 + \+ blah]), i = 1..nops(V0));\n\n Vothersum1 := add(T[coord_tra nsfer(Vother1[i], triangle1, S[basis[1], 1], S[basis[1], 2], basis[1]) , funct] * (T[Vother1[i], harm_funcno + 9 + blah]), i = 1..nops(Vother 1)); \n \n Vmsum1 := add(T[coord_transfer(Vm1[i], tria ngle1, S[basis[1], 1], S[basis[1], 2], basis[1]), funct] * (T[Vm1[i], \+ harm_funcno + 9 + blah]), i = 1..nops(Vm1));\n\n #and now o n 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 + blah]), i = 1..nops(V0));\n\n Vot hersum2 := add(T[coord_transfer(Vother2[i], triangle2, basis[1], S[bas is[1], 3], S[basis[1], 4]), funct] * (T[Vother2[i], harm_funcno + 1 + \+ blah]), i = 1..nops(Vother2)); \n\n Vmsum2 := add(T[ coord_transfer(Vm2[i], triangle2, basis[1], S[basis[1], 3], S[basis[1] , 4]), funct] * (T[Vm2[i], harm_funcno + 1 + blah]), i = 1..nops(Vm2)) ;\n\n elif basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] t hen\n \n #let's calculate sum's for triangle 1, shall \+ we?\n V0sum1 := add(T[coord_transfer(V0[i], triangle1, S[ba sis[1], 2], basis[1], S[basis[1], 1]), funct] * (T[V0[i], harm_funcno \+ + 5 + blah]), i = 1..nops(V0));\n\n Vothersum1 := add(T[coo rd_transfer(Vother1[i], triangle1, S[basis[1], 2], basis[1], S[basis[1 ], 1]), funct] * (T[Vother1[i], harm_funcno + 5 + blah]), i = 1..nops( Vother1)); \n \n Vmsum1 := add(T[coord_transfer(Vm1[i] , triangle1, S[basis[1], 2], basis[1], S[basis[1], 1]), funct] * (T[Vm 1[i], harm_funcno + 5 + blah]), i = 1..nops(Vm1));\n\n #and now on to triangle numero 2!\n V0sum2 := add(T[coord_trans fer(V0[i], triangle2, S[basis[1], 3], S[basis[1], 4], basis[1]), funct ] * (T[V0[i], harm_funcno + 9 + blah]), i = 1..nops(V0));\n\n \+ Vothersum2 := add(T[coord_transfer(Vother2[i], triangle2, S[basis[1 ], 3], S[basis[1], 4], basis[1]), funct] * (T[Vother2[i], harm_funcno \+ + 9 + blah]), i = 1..nops(Vother2)); \n\n Vmsum2 := \+ add(T[coord_transfer(Vm2[i], triangle2, S[basis[1], 3], S[basis[1], 4] , basis[1]), funct] * (T[Vm2[i], harm_funcno + 9 + blah]), i = 1..nops (Vm2));\n\n fi; \n\n #normal derivative scaling\n \+ if basis[2] = g then\n V0sum1 := V0sum1 * (3/5) ^ m1;\n \+ Vothersum1 := Vothersum1 * (3/5) ^ m1;\n Vmsum1 := Vmsum1 * (3/5) ^ m1;\n V0sum2 := V0sum2 * (3/5) ^ m2;\n \+ Vothersum2 := Vothersum2 * (3/5) ^ m2;\n Vmsum2 := Vmsum2 * (3/5) ^ m2;\n fi;\n \n #Return simpson' s rule. Notice that the second sum must be subtracted if it's a g-typ e vector\n #because the normal derivative is defined as negativ e in the 2nd triangle...\n if basis[2] = f then RETURN(fact * ( 5*(Vmsum1+Vmsum2) + 2*(Vothersum1+Vothersum2) + (V0sum1+V0sum2))); fi; \n if basis[2] = g then RETURN(fact * (5*(Vmsum1-Vmsum2) + 2*(V othersum1-Vothersum2) + (V0sum1-V0sum2))); fi;\n\n fi;\n\n\nend;" } }{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3int_vs_basis_irregGf*6(%\"TG%\"SG% 'm_funcG%&basisG%&functG%,harm_funcnoG66%#m1G%#m2G%#V0G%(Vother1G%(Vot her2G%$Vm1G%$Vm2G%%factG%*triangle1G%*triangle2G%'V0sum1G%'V0sum2G%+Vo thersum1G%+Vothersum2G%'Vmsum1G%'Vmsum2G%\"iG%\"jG%)Vother_eG%%blahG6 \"FBC+>8+,$*&\"\"\"FH)\"\"$9&!\"\"#FH\"\"'>8$,&-%$maxG6$-%%nopsG6#&9'6 #FH-FV6#&9%6$FXFHFHFHFL>8%,&-FS6$FU-FV6#&Fhn6$FX\"\"%FHFHFL>8&-%+verti ces_mG6#\"\"!@%1,&FKFHFPFLFH>8'7\"C$>86-%)verticesG6#,(FKFHFPFLFHFL>F^ p&Fbp6#;Fco-FV6#Fbp@%1,&FKFHF[oFLFH>8(F_pC$>Fbp-Fdp6#,(FKFHF[oFLFHFL>F aqFhp@%/FKFP>8)F_p>F[r-Fgo6#F\\p@%/FKF[o>8*F_p>Fbr-Fgo6#F_q@%/FUFHC*>8 ,7#-%$seqG6$&FXFZ/85;FHFP@$/FX7#FioC%>8.-%$addG6$*&&9$6$-%/coord_trans ferG6'&Feo6#84FjrFXFgn&Fhn6$FX\"\"#9(FH&F^t6$Fct,&9)FHFJFHFH/Fet;FH-FV 6#Feo>80-Fjs6$*&&F^t6$-Fat6'&F^pFdtFjrFXFgnFftFitFH&F^t6$F[vF\\uFH/Fet ;FH-FV6#F^p>82-Fjs6$*&&F^t6$-Fat6'&F[rFdtFjrFXFgnFftFitFH&F^t6$F[wF\\u FH/Fet;FH-FV6#F[r@$/FX7#FHC%>Fhs-Fjs6$*&&F^t6$-Fat6'FctFjrFftFXFgnFitF H&F^t6$Fct,&F]uFH\"\"(FHFHF^u>Fcu-Fjs6$*&&F^t6$-Fat6'F[vFjrFftFXFgnFit FH&F^t6$F[vF`xFHF^v>Fcv-Fjs6$*&&F^t6$-Fat6'F[wFjrFftFXFgnFitFH&F^t6$F[ wF`xFHF^w@$/FX7#FhtC%>Fhs-Fjs6$*&&F^t6$-Fat6'FctFjrFgnFftFXFitFH&F^t6$ Fct,&F]uFH\"#6FHFHF^u>Fcu-Fjs6$*&&F^t6$-Fat6'F[vFjrFgnFftFXFitFH&F^t6$ F[vFdzFHF^v>Fcv-Fjs6$*&&F^t6$-Fat6'F[wFjrFgnFftFXFitFH&F^t6$F[wFdzFHF^ w>Fhs*&FhsFH)#FJ\"\"&FPFH>Fcu*&FcuFHF\\\\lFH>Fcv*&FcvFHF\\\\lFH-%'RETU RNG6#*&FEFH,(FcvF^\\l*&FhtFHFcuFHFHFhsFHFHC)>Fjr-%-get_triangleG6&FhnF YFPFH>8--F\\]l6&FhnFYF[oFht@%/&FY6#Fht%\"fG>87Fio>Fh]lFht@(/&FX6#;,&FU FHFHFLFU7$FioFHC(>Fhs-Fjs6$*&F]tFH&F^t6$Fct,(F]uFHFHFHFh]lFHFHF^u>Fcu- Fjs6$*&FguFH&F^t6$F[vFh^lFHF^v>Fcv-Fjs6$*&FgvFH&F^t6$F[wFh^lFHF^w>8/-F js6$*&&F^t6$-Fat6'FctF_]lFaoFX&Fhn6$FXFJFitFH&F^t6$Fct,(F]uFHF^\\lFHFh ]lFHFHF^u>81-Fjs6$*&&F^t6$-Fat6'&FaqFdtF_]lFaoFXF^`lFitFH&F^t6$F\\alFb `lFH/Fet;FH-FV6#Faq>83-Fjs6$*&&F^t6$-Fat6'&FbrFdtF_]lFaoFXF^`lFitFH&F^ t6$F\\blFb`lFH/Fet;FH-FV6#Fbr/F\\^l7$FioFhtC(>Fhs-Fjs6$*&F^zFH&F^t6$Fc t,(F]uFH\"\"*FHFh]lFHFHF^u>Fcu-Fjs6$*&FjzFH&F^t6$F[vF\\clFHF^v>Fcv-Fjs 6$*&Fd[lFH&F^t6$F[wF\\clFHF^w>Ff_l-Fjs6$*&&F^t6$-Fat6'FctF_]lFXF^`lFao FitFHFf^lFHF^u>Fd`l-Fjs6$*&&F^t6$-Fat6'F\\alF_]lFXF^`lFaoFitFH&F^t6$F \\alFh^lFHF_al>Fdal-Fjs6$*&&F^t6$-Fat6'F\\blF_]lFXF^`lFaoFitFH&F^t6$F \\blFh^lFHF_bl/F\\^l7$FHFhtC(>Fhs-Fjs6$*&FjwFHF``lFHF^u>Fcu-Fjs6$*&Ffx FH&F^t6$F[vFb`lFHF^v>Fcv-Fjs6$*&F`yFH&F^t6$F[wFb`lFHF^w>Ff_l-Fjs6$*&&F ^t6$-Fat6'FctF_]lF^`lFaoFXFitFHFjblFHF^u>Fd`l-Fjs6$*&&F^t6$-Fat6'F\\al F_]lF^`lFaoFXFitFH&F^t6$F\\alF\\clFHF_al>Fdal-Fjs6$*&&F^t6$-Fat6'F\\bl F_]lF^`lFaoFXFitFH&F^t6$F\\blF\\clFHF_bl@$/Fd]l%\"gGC(>FhsF[\\l>FcuF` \\l>FcvFb\\l>Ff_l*&Ff_lFH)F]\\lF[oFH>Fd`l*&Fd`lFHF^ilFH>Fdal*&FdalFHF^ ilFH@$Fc]l-Fd\\l6#*&FEFH,.FcvF^\\l*&F^\\lFHFdalFHFH*&FhtFHFcuFHFH*&Fht FHFd`lFHFHFhsFHFf_lFHFH@$Ffhl-Fd\\l6#*&FEFH,.FcvF^\\l*&F^\\lFHFdalFHFL *&FhtFHFcuFHFH*&FhtFHFd`lFHFLFhsFHFf_lFLFHFBFBFB" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 359 "# same as gen_f_vector_irreg_h but for bihar monic splines\ngen_f_vector_irreg := proc(T,S,m_func,vm,funcno,harm_fu ncno)\n\n local i, f_vec, n_verts;\n\n n_verts := 2 * nops(vm) - 3;\n\n f_vec := array(1..n_verts);\n\n for i from 1 to n_verts \+ do\n f_vec[i] := int_vs_basis_irreg(T,S,m_func,get_basis(vm,i), funcno,harm_funcno);\n od;\n\n f_vec;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3gen_f_vector_irregGf*6(%\"TG%\"SG%'m_funcG%#vmG%'fun cnoG%,harm_funcnoG6%%\"iG%&f_vecG%(n_vertsG6\"F1C&>8&,&-%%nopsG6#9'\" \"#\"\"$!\"\">8%-%&arrayG6#;\"\"\"F4?(8$FCFCF4%%trueG>&F>6#FE-%3int_vs _basis_irregG6(9$9%9&-%*get_basisG6$F9FE9(9)F>F1F1F1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7670 "#This function creates an actual functi on (which can then be plotted and such) from the values\n#of the splin e projection stored in C. It creates the function simply by going thr ough each basis\n#element and adding in the appropriate amount of it. \+ Sounds, simple, huh?\n\n# for the irreg version of this function, the convention i've adopted is to \n# plot to some uniform depth past the basis, which is why we use W and depth in all of our plotting functio ns ... \n\nmake_pretty_picture_irreg := proc(T, S, W, depth, vm, funcn o, harm_funcno, C)\n\n local i, j, vm2, k, vmh, triangle1, triangle 2, basis, add0, add1, add2, which_harm1, coord1, fg1, fg2, fg3_1, fg3_ 2, m1, m2, tempW;\n\n\n #first, let's set the initial values of our function to 0.\n tempW := W;\n for i from 1 to depth do\n \+ tempW := subdivide_irreg(tempW);\n od;\n vm2 := vertices_irreg (tempW);\n for i from 1 to nops(vm2) do\n T[vm2[i], funcno] \+ := 0;\n od;\n\n for i from 1 to 2 * nops(vm) - 3 do\n\n b asis:=get_basis(vm,i); #This is the basis element we will be working \+ with.\n \n m1 := max(nops(basis[1]), nops(S[basis[1], 1] )) - 1;\n m2 := max(nops(basis[1]), nops(S[basis[1], 4])) - 1; \n\n #So if the basis element is in V0, we know it is a g-type \+ basis vector.\n #the scheme here is to find the appropriate tri angle (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 #als o, depending on the orientation of the basis, we need to pick the appr opriate\n #basis element (so for basis = [[0],g], we need a bas is element which has values\n #of zero everywhere, and has norm al derivative 1 at [0], and normal derivative 0\n #at the other boundaries. This is stored in harm_funcno+3). This explains the ass ignments\n #of triangle1, add0, add1, add2, and which_harm1. T he function coord_transfer 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 to generate (basically, by adding \+ some more digits to\n #the front). Then add in the value of th e basis function * C[i] and voila!\n\n if nops(basis[1]) = 1 th en #this is if the basis is an element of V0...\n if basis[ 1] = [0] then\n triangle1 := [seq(0, k = 1..m1)];\n \+ add0 := basis[1]; \n add1 := S[basis[1], 1] ; \n add2 := S[basis[1], 2];\n which_har m1 := 3; \n fi;\n if basis[1] = [1] then\n \+ triangle1 := [seq(1, k = 1..m1)];\n add0 := S [basis[1], 2]; \n add1 := basis[1]; \n a dd2 := S[basis[1], 1];\n which_harm1 := 7;\n \+ fi;\n if basis[1] = [2] then\n triangle1 := [seq(2, k = 1..m1)];\n add0 := S[basis[1], 1]; \n \+ add1 := S[basis[1], 2]; \n add2 := basis[1]; \n which_harm1 := 11;\n fi;\n \n \+ vmh := vertices(depth);\n for j from 1 to nops(vmh) do\n coord1 := coord_transfer(vmh[j], triangle1, add0, add1 , add2);\n T[coord1,funcno] := T[coord1, funcno] + C[i] * ((3/5) ^ m1) * T[vmh[j], harm_funcno + which_harm1];\n o d; \n\n else\n\n #Here, the strategy is the same, ex cept 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_t ransfer as the triangle coordinates\n #are different depend ing on the orientation of the basis element (ie, what the last 2 digit s are).\n #Also, the basis functions we want to use change \+ as well, again depending on the last 2 digits.\n #All of th is hinges around the fact that the neighbors are ordered in a clockwis e fashion, so that the\n #orientation of the basis element \+ tells you which neighbors are the vertices of which triangle.\n\n \+ triangle1 := get_triangle(S, basis, m1, 1);\n triang le2 := get_triangle(S, basis, m2, 2);\n\n vmh := vertices(d epth);\n\n #set flags for f vs g ...\n # fg1 is \+ because the normal derivative is negative in the second\n # triangle, so g-types need a - instead of a + at times\n # \+ fg2 is to select the correct slot in the table to get basis\n \+ # function values\n # fg3's are scaling for g-type value s \n if basis[2] = f then\n fg1 := 1;\n \+ fg2 := 0;\n fg3_1 := 1;\n fg3 _2 := 1;\n else\n fg1 := -1;\n \+ fg2 := 2;\n fg3_1 := (3/5) ^ m1;\n fg3 _2 := (3/5) ^ m2;\n fi;\n\n if basis[1][nops(bas is[1])-1..nops(basis[1])] = [0,1] then \n \n \+ for j from 1 to nops(vmh) do\n T[coord_transfer(vm h[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] * fg3_1 * (T[vmh[j], harm_funcno + 1 + fg2]); \n od;\n \n for j from 1 \+ to nops(vmh) do\n T[coord_transfer(vmh[j], triangle 2, S[basis[1], 4], basis[1], S[basis[1], 3]), funcno] := T[coord_trans fer(vmh[j], triangle2, S[basis[1], 4], basis[1], S[basis[1], 3]), func no] + fg1 * C[i] * fg3_2 * (T[vmh[j], harm_funcno + 5 + fg2]);\n \+ od;\n \n\n #we added in \+ the value of the basis here twice, so subtract one away (f case only) \n if basis[2] = f then\n T[basis[1] , funcno] := T[basis[1], funcno] - C[i];\n fi;\n \+ \n fi;\n\n if basis[1][nops(basis[1])- 1..nops(basis[1])] = [0,2] then\n \n for j from 1 \+ to nops(vmh) do\n T[coord_transfer(vmh[j], triangle 1, S[basis[1], 1], S[basis[1], 2], basis[1]), funcno] := T[coord_trans fer(vmh[j], triangle1, S[basis[1], 1], S[basis[1], 2], basis[1]), func no] + C[i] * fg3_1 * (T[vmh[j], harm_funcno + 9 + fg2]);\n \+ od;\n for j from 1 to nops(vmh) do\n \+ T[coord_transfer(vmh[j], triangle2, basis[1], S[basis[1], 3], S[b asis[1], 4]), funcno] := T[coord_transfer(vmh[j], triangle2, basis[1], S[basis[1], 3], S[basis[1], 4]), funcno] + fg1 * C[i] * fg3_2 * (T[vm h[j], harm_funcno + 1 + fg2]);\n od;\n \n \+ if basis[2] = f then\n T[basis[1], funcno] : = T[basis[1], funcno] - C[i];\n fi;\n fi; \n \n if basis[1][nops(basis[1])-1..nops(basis[1])] = [1,2] th en\n \n for j from 1 to nops(vmh) do\n \+ T[coord_transfer(vmh[j], triangle1, S[basis[1], 2], basis [1], S[basis[1], 1]), funcno] := T[coord_transfer(vmh[j], triangle1, S [basis[1], 2], basis[1], S[basis[1], 1]), funcno] + C[i] * fg3_1 * (T[ vmh[j], harm_funcno + 5 + fg2]);\n od;\n \+ for j from 1 to nops(vmh) do\n T[coord_transfer(vm h[j], triangle2, S[basis[1], 3], S[basis[1], 4], basis[1]), funcno] := T[coord_transfer(vmh[j], triangle2, S[basis[1], 3], S[basis[1], 4], b asis[1]), funcno] + fg1 * C[i] * fg3_2 * (T[vmh[j], harm_funcno + 9 + \+ fg2]);\n od;\n \n if basis[2] = f t hen\n T[basis[1], funcno] := T[basis[1], funcno] - \+ C[i];\n fi;\n fi; \n fi;\n od;\nen d;\n \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%:make_pretty_pict ure_irregGf*6*%\"TG%\"SG%\"WG%&depthG%#vmG%'funcnoG%,harm_funcnoG%\"CG 66%\"iG%\"jG%$vm2G%\"kG%$vmhG%*triangle1G%*triangle2G%&basisG%%add0G%% add1G%%add2G%,which_harm1G%'coord1G%$fg1G%$fg2G%&fg3_1G%&fg3_2G%#m1G%# m2G%&tempWG6\"FDC'>879&?(8$\"\"\"FK9'%%trueG>FG-%0subdivide_irregG6#FG >8&-%/vertices_irregGFQ?(FJFKFK-%%nopsG6#FSFM>&9$6$&FS6#FJ9)\"\"!?(FJF KFK,&-FX6#9(\"\"#\"\"$!\"\"FMC&>8+-%*get_basisG6$F`oFJ>85,&-%$maxG6$-F X6#&Ffo6#FK-FX6#&9%6$FbpFKFKFKFco>86,&-F^p6$F`p-FX6#&Fgp6$Fbp\"\"%FKFK Fco@%/F`pFKC'@$/Fbp7#F[oC'>8)7#-%$seqG6$F[o/8';FKF[p>8,Fbp>8-Ffp>8.&Fg p6$FbpFao>8/Fbo@$/Fbp7#FKC'>F[r7#-F^r6$FKF`r>FdrFir>FfrFbp>FhrFfp>F\\s \"\"(@$/Fbp7#FaoC'>F[r7#-F^r6$FaoF`r>FdrFfp>FfrFir>FhrFbp>F\\s\"#6>8(- %)verticesG6#FL?(8%FKFK-FX6#FhtFMC$>80-%/coord_transferG6'&Fht6#F]uF[r FdrFfrFhr>&Ffn6$FbuFjn,&FiuFK*(&9+FinFK)#Fbo\"\"&F[pFK&Ffn6$Ffu,&9*FKF \\sFKFKFKC)>F[r-%-get_triangleG6&FgpFfoF[pFK>8*-Fiv6&FgpFfoFjpFao>FhtF it@%/&Ffo6#Fao%\"fGC&>81FK>82F[o>83FK>84FKC&>FgwFco>FiwFao>F[xF_v>F]x) F`vFjp@$/&Fbp6#;,&F`pFKFKFcoF`p7$F[oFKC%?(F]uFKFKF^uFM>&Ffn6$-Fdu6'Ffu F[rFbpFfpFirFjn,&F^yFK*(F]vFKF[xFK&Ffn6$Ffu,(FevFKFKFKFiwFKFKFK?(F]uFK FKF^uFM>&Ffn6$-Fdu6'FfuF\\wF`qFbp&Fgp6$FbpFboFjn,&FiyFK**FgwFKF]vFKF]x FK&Ffn6$Ffu,(FevFKFavFKFiwFKFKFK@$Faw>&Ffn6$FbpFjn,&FfzFKF]vFco@$/Ffx7 $F[oFaoC%?(F]uFKFKF^uFM>&Ffn6$-Fdu6'FfuF[rFfpFirFbpFjn,&F_[lFK*(F]vFKF [xFK&Ffn6$Ffu,(FevFK\"\"*FKFiwFKFKFK?(F]uFKFKF^uFM>&Ffn6$-Fdu6'FfuF\\w FbpF]zF`qFjn,&F[\\lFK**FgwFKF]vFKF]xFKFdyFKFK@$Faw>FfzFhz@$/Ffx7$FKFao C%?(F]uFKFKF^uFM>&Ffn6$-Fdu6'FfuF[rFirFbpFfpFjn,&Fi\\lFK*(F]vFKF[xFKFa zFKFK?(F]uFKFKF^uFM>&Ffn6$-Fdu6'FfuF\\wF]zF`qFbpFjn,&Fa]lFK**FgwFKF]vF KF]xFKFe[lFKFK@$Faw>FfzFhzFDFDFD" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12472 "# Rather huge function to integrate a function aga inst 2 biharmonic basis\n# functions. \nint_vs_basis12_irreg := proc(T , S, m_func, basis_1, basis_2, funct, harm_funcno)\n\n local out, t emp, basis1, basis2, triangle1, triangle2, addresses1, addresses2, j, \+ which_harms, m1, m2, fact1, fact2, g_fact1, g_fact2;\n\n basis1 := \+ basis_1; #again, so that we can switch around basis1 and basis2, if n eed be.\n basis2 := basis_2;\n\n # set various scaling factors a nd such\n m1 := max(nops(basis1[1]), nops(S[basis1[1], 1])) - 1;\n \+ m2 := max(nops(basis1[1]), nops(S[basis1[1], 4])) - 1;\n fact1 : = (1/3)^m1;\n fact2 := (1/3)^m2;\n g_fact1 := (3/5)^m1;\n g_f act2 := (3/5)^m2;\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[bas is1[1],4] <> basis2[1] then RETURN(0) fi;\n\n #first deal with spec ial case of either 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 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 := basis2;\n \+ basis2 := basis1;\n basis1 := temp;\n fi;\n\n if basis1[1] = basis2[1] then #The only way we'll ever get h ere is if both elements are g-type basis vectors\n \n \+ if basis1[1] = [0] then \n triangle1 := [seq(0, j = 1..m1)];\n out := g_fact1^2 * fact1 * integrate_basis1 2(T, m_func - m1, triangle1, basis1[1], S[basis1[1], 1], S[basis1[1], \+ 2], funct, harm_funcno + 3, harm_funcno + 3);\n fi;\n \+ if basis1[1] = [1] then\n triangle1 := [seq(1, j \+ = 1..m1)];\n out := g_fact1^2 * fact1 * integrate_basis 12(T, m_func - m1, 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 triangle1 := [seq(2, \+ j = 1..m1)];\n out := g_fact1^2 * fact1 * integrate_bas is12(T, m_func - m1, triangle1, S[basis1[1], 1], S[basis1[1], 2], basi s1[1], funct, harm_funcno + 11, harm_funcno + 11);\n fi; \+ \n elif basis1[2] = g then\n if S[basis1[1],1] = basis2[1] then\n if basis2[1] = [0] then \n \+ triangle1 := [seq(0, j = 1..m1)];\n out \+ := g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, ba sis2[1], S[basis2[1], 1], S[basis2[1], 2], funct, harm_funcno + 3, har m_funcno + 11);\n fi;\n if basis2[1] = [ 1] then\n triangle1 := [seq(1, j = 1..m1)];\n \+ out := g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, S[basis2[1], 2], basis2[1], S[basis2[1], 1], funct, ha rm_funcno + 7, harm_funcno + 3); \n fi;\n \+ if basis2[1] = [2] then\n triangle1 := [seq(2, j \+ = 1..m1)];\n out := g_fact1^2 * fact1 * integrate_b asis12(T, m_func - m1, triangle1, S[basis2[1], 1], S[basis2[1], 2], ba sis2[1], funct, harm_funcno + 11,harm_funcno + 7);\n fi ;\n fi;\n \n if S[basis1[1],4] = basis2[1] \+ then\n if basis2[1] = [0] then \n tr iangle1 := [seq(0, j = 1..m1)];\n out := -g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, basis2[1], S[b asis2[1], 1], S[basis2[1], 2], funct, harm_funcno + 3, harm_funcno + 7 );\n fi;\n if basis2[1] = [1] then\n \+ triangle1 := [seq(1, j = 1..m1)];\n \+ out := -g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangl e1, 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..m1)]; \n out := -g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, S[basis2[1], 1], S[basis2[1], 2], basis2[1], \+ funct, harm_funcno + 11, harm_funcno + 3);\n fi;\n \+ fi;\n #last remaining case: basis1 is f, basis2 is g\n \+ else\n \n if S[basis1[1],1] = basis2[1] th en\n if basis2[1] = [0] then \n tria ngle1 := [seq(0, j = 1..m1)];\n out := g_fact1 * fa ct1 * integrate_basis12(T, m_func - m1, triangle1, basis2[1], S[basis2 [1], 1], S[basis2[1], 2], funct, harm_funcno + 3, harm_funcno + 9);\n \+ fi;\n if basis2[1] = [1] then\n \+ triangle1 := [seq(1, j = 1..m1)];\n out \+ := g_fact1 * fact1 * integrate_basis12(T, m_func - m1, triangle1, S[ba sis2[1], 2], basis2[1], S[basis2[1], 1], funct, harm_funcno + 7, harm_ funcno + 1); \n fi;\n if basis2[1] = [2] then\n triangle1 := [seq(2,j = 1 ..m1)];\n \+ out := g_fact1 * fact1 * integrate_basis12(T, m_func - m1, triangle1, S[basis2[1], 1], S[basis2[1], 2], basis2[1], funct, harm_f uncno + 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 ..m1)];\n out := g_fact1 * fact1 * integrate_basis1 2(T, m_func - m1, 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 triangle1 : = [seq(1, j = 1..m1)];\n out := g_fact1 * fact1 * i ntegrate_basis12(T, m_func - m1, 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..m1)];\n out := g_f act1 * fact1 * integrate_basis12(T, m_func - m1, triangle1, S[basis2[1 ], 1], S[basis2[1], 2], basis2[1], funct, harm_funcno + 11, harm_funcn o + 1);\n fi;\n fi;\n\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 o nly thing that matters is the orientation to determine which basis fun ctions to use...\n #first, let's get the triangles, as well as the addresses of the vertices of those triangles\n tria ngle1 := get_triangle(S, basis1, m1, 1); \n triangle2 := \+ get_triangle(S, basis1, m2, 2);\n addresses1 := get_address es(S, basis1, m1, 1);\n addresses2 := get_addresses(S, basi s1, m2, 2);\n \n if (basis1[2] = f and basis2[2] = g) \+ or (basis1[2] = g and basis2[2] = f) then\n if last_dig its(basis1[1]) = [0,1] then\n out := g_fact1 * fact 1 * integrate_basis12(T, m_func - m1, triangle1, addresses1[1], addres ses1[2], addresses1[3], funct, harm_funcno + 1, harm_funcno + 3) - g_f act2 * fact2 * integrate_basis12(T, m_func - m2, triangle2, addresses2 [1], addresses2[2], addresses2[3], funct, harm_funcno + 5, harm_funcno + 7); \n fi;\n\n if last_digits(basis1[ 1]) = [0,2] then\n out := g_fact1 * fact1 * integra te_basis12(T, m_func - m1, triangle1, addresses1[1], addresses1[2], ad dresses1[3], funct, harm_funcno + 9, harm_funcno + 11) - g_fact2 * fac t2 * integrate_basis12(T, m_func - m2, triangle2, addresses2[1], addre sses2[2], addresses2[3], funct, harm_funcno + 1, harm_funcno + 3); \n \+ fi;\n\n if last_digits(basis1[1]) = [1,2 ] then\n out := g_fact1 * fact1 * integrate_basis12 (T, m_func - m1, triangle1, addresses1[1], addresses1[2], addresses1[3 ], funct, harm_funcno + 5, harm_funcno + 7) - g_fact2 * fact2 * integr ate_basis12(T, m_func - m2, triangle2, addresses2[1], addresses2[2], a ddresses2[3], funct, harm_funcno + 9, harm_funcno + 11); \n \+ fi;\n fi;\n\n if basis1[2] = f and basis2[2 ] = f then \n if last_digits(basis1[1]) = [0,1] then\n \+ out := fact1 * integrate_basis12(T, m_func - m1, tr iangle1, addresses1[1], addresses1[2], addresses1[3], funct, harm_func no + 1, harm_funcno + 1) + fact2 * integrate_basis12(T, m_func - m2, t riangle2, addresses2[1], addresses2[2], addresses2[3], funct, harm_fun cno + 5, harm_funcno + 5); \n fi;\n\n if last_digits(basis1[1]) = [0,2] then\n out := fact1 * integrate_basis12(T, m_func - m1, triangle1, addresses1[1], address es1[2], addresses1[3], funct, harm_funcno + 9, harm_funcno + 9) + fact 2 * integrate_basis12(T, m_func - m2, triangle2, addresses2[1], addres ses2[2], addresses2[3], funct, harm_funcno + 1, harm_funcno + 1); \n \+ fi;\n\n if last_digits(basis1[1]) = [1,2] then\n out := fact1 * integrate_basis12(T, m_func \+ - m1, triangle1, addresses1[1], addresses1[2], addresses1[3], funct, h arm_funcno + 5, harm_funcno + 5) + fact2 * integrate_basis12(T, m_func - m2, triangle2, addresses2[1], addresses2[2], addresses2[3], funct, \+ harm_funcno + 9, harm_funcno + 9); \n fi;\n \+ fi;\n\n if basis1[2] = g and basis2[2] = g then\n \+ if last_digits(basis1[1]) = [0,1] then\n out \+ := g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, ad dresses1[1], addresses1[2], addresses1[3], funct, harm_funcno + 3, har m_funcno + 3) + g_fact2^2 * fact2 * integrate_basis12(T, m_func - m2, \+ triangle2, addresses2[1], addresses2[2], addresses2[3], funct, harm_fu ncno + 7, harm_funcno + 7); \n fi;\n\n i f last_digits(basis1[1]) = [0,2] then\n out := g_fa ct1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, addresses 1[1], addresses1[2], addresses1[3], funct, harm_funcno + 11, harm_func no + 11) + g_fact2^2 * fact2 * integrate_basis12(T, m_func - m2, trian gle2, addresses2[1], addresses2[2], addresses2[3], funct, harm_funcno \+ + 3, harm_funcno + 3); \n fi;\n\n if las t_digits(basis1[1]) = [1,2] then\n out := g_fact1^2 * fact1 * integrate_basis12(T, m_func - m1, triangle1, addresses1[1], addresses1[2], addresses1[3], funct, harm_funcno + 7, harm_funcno + 7 ) + g_fact2^2 * fact2 * integrate_basis12(T, m_func - m2, triangle2, a ddresses2[1], addresses2[2], addresses2[3], funct, harm_funcno + 11, h arm_funcno + 11); \n fi;\n fi;\n #END OF CASE WHERE BASIS1 and BASIS2 are in same place \n \n\n el se\n #if they are neighbors, then we're in here. Notice th at we just call\n #which_funcs to determine which basis fun ctions 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] = basi s2[1] then\n #first, let's get the triangles, as well a s the addresses of the vertices of those triangles\n tr iangle1 := get_triangle(S, basis1, m1, 1); \n address es1 := get_addresses(S, basis1, m1, 1);\n which_harms : = which_funcs(S, basis1, basis2);\n out := which_harms[ 3] * fact1 * integrate_basis12(T, m_func - m1, triangle1, addresses1[1 ], addresses1[2], addresses1[3], funct, harm_funcno + which_harms[1], \+ harm_funcno + which_harms[2]);\n if basis1[2] = g then \n out := g_fact1 * out;\n fi;\n \+ if basis2[2] = g then\n out := g_fact1 \+ * out;\n fi;\n else\n triangl e1 := get_triangle(S, basis1, m2, 2);\n addresses1 := g et_addresses(S, basis1, m2, 2);\n which_harms := which_ funcs(S, basis1, basis2);\n out := which_harms[3] * fac t2 * integrate_basis12(T, m_func - m2, triangle1, addresses1[1], addre sses1[2], addresses1[3], funct, harm_funcno + which_harms[1], harm_fun cno + which_harms[2]);\n if basis1[2] = g then\n \+ out := g_fact2 * out;\n fi;\n \+ if basis2[2] = g then\n out := g_fact2 * out;\n \+ fi;\n fi;\n\n fi;\n \n fi;\n \+ \n out; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5int_vs_ basis12_irregGf*6)%\"TG%\"SG%'m_funcG%(basis_1G%(basis_2G%&functG%,har m_funcnoG62%$outG%%tempG%'basis1G%'basis2G%*triangle1G%*triangle2G%+ad dresses1G%+addresses2G%\"jG%,which_harmsG%#m1G%#m2G%&fact1G%&fact2G%(g _fact1G%(g_fact2G6\"F?C->8&9'>8'9(>8.,&-%$maxG6$-%%nopsG6#&FB6#\"\"\"- FN6#&9%6$FPFRFRFR!\"\">8/,&-FK6$FM-FN6#&FV6$FP\"\"%FRFRFX>80)#FR\"\"$F H>81)F`oFZ>82)#Fao\"\"&FH>83)FhoFZ@$33330FP&FEFQ0FUFcp0&FV6$FP\"\"#Fcp 0&FV6$FPFaoFcp0FjnFcp-%'RETURNG6#\"\"!@%55555/FP7#F`q/FP7#FR/FP7#Fhp/F cpFhq/FcpFjq/FcpF\\rC$@$/FMFRC%>8%FE>FEFB>FBFer@'/FPFcpC%@$FgqC$>8(7#- %$seqG6$F`q/8,;FRFH>8$*()FfoFhpFRF^oFR-%2integrate_basis12G6+9$,&9&FRF HFXF^sFPFUFfp9),&9*FRFaoFRFatFR@$FiqC$>F^s7#-Fas6$FRFcs>Fgs*(FisFRF^oF R-F[t6+F]tF^tF^sFfpFPFUF`t,&FbtFR\"\"(FRF]uFR@$F[rC$>F^s7#-Fas6$FhpFcs >Fgs*(FisFRF^oFR-F[t6+F]tF^tF^sFUFfpFPF`t,&FbtFR\"#6FRFiuFR/&FB6#Fhp% \"gGC$@$/FUFcpC%@$F]rC$>F^sF_s>Fgs*(FisFRF^oFR-F[t6+F]tF^tF^sFcp&FV6$F cpFR&FV6$FcpFhpF`tFatFiuFR@$F^rC$>F^sFft>Fgs*(FisFRF^oFR-F[t6+F]tF^tF^ sF\\wFcpFjvF`tF]uFatFR@$F_rC$>F^sFbu>Fgs*(FisFRF^oFR-F[t6+F]tF^tF^sFjv F\\wFcpF`tFiuF]uFR@$/FjnFcpC%@$F]rC$>F^sF_s>Fgs,$*(FisFRF^oFR-F[t6+F]t F^tF^sFcpFjvF\\wF`tFatF]uFRFX@$F^rC$>F^sFft>Fgs,$*(FisFRF^oFR-F[t6+F]t F^tF^sF\\wFcpFjvF`tF]uFiuFRFX@$F_rC$>F^sFbu>Fgs,$*(FisFRF^oFR-F[t6+F]t F^tF^sFjvF\\wFcpF`tFiuFatFRFXC$@$FavC%@$F]rC$>F^sF_s>Fgs*(FfoFRF^oFR-F [t6+F]tF^tF^sFcpFjvF\\wF`tFat,&FbtFR\"\"*FRFR@$F^rC$>F^sFft>Fgs*(FfoFR F^oFR-F[t6+F]tF^tF^sF\\wFcpFjvF`tF]u,&FbtFRFRFRFR@$F_rC$>F^sFbu>Fgs*(F foFRF^oFR-F[t6+F]tF^tF^sFjvF\\wFcpF`tFiu,&FbtFRFioFRFR@$F]xC%@$F]rC$>F ^sF_s>Fgs*(FfoFRF^oFR-F[t6+F]tF^tF^sFcpFjvF\\wF`tFatFb[lFR@$F^rC$>F^sF ft>Fgs*(FfoFRF^oFR-F[t6+F]tF^tF^sF\\wFcpFjvF`tF]uFazFR@$F_rC$>F^sFbu>F gs*(FfoFRF^oFR-F[t6+F]tF^tF^sFjvF\\wFcpF`tFiuFjzFRC$@$FirC)>F^s-%-get_ triangleG6&FVFBFHFR>8)-F_]l6&FVFBFZFhp>8*-%.get_addressesGF`]l>8+-Fh]l Fd]l@$53/F\\v%\"fG/&FEF]vF^v3F[v/Fb^lF`^lC%@$/-%,last_digitsGFO7$F`qFR >Fgs,&*(FfoFRF^oFR-F[t6+F]tF^tF^s&Ff]lFQ&Ff]lF]v&Ff]l6#FaoF`tFjzFatFRF R*(F[pFRFcoFR-F[t6+F]t,&F_tFRFZFXFb]l&Fj]lFQ&Fj]lF]v&Fj]lFc_lF`tFb[lF] uFRFX@$/Fh^l7$F`qFhp>Fgs,&*(FfoFRF^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tFa zFiuFRFR*(F[pFRFcoFR-F[t6+F]tFg_lFb]lFh_lFi_lFj_lF`tFjzFatFRFX@$/Fh^l7 $FRFhp>Fgs,&*(FfoFRF^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tFb[lF]uFRFR*(F[p FRFcoFR-F[t6+F]tFg_lFb]lFh_lFi_lFj_lF`tFazFiuFRFX@$3F_^lFd^lC%@$Fg^l>F gs,&*&F^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tFjzFjzFRFR*&FcoFR-F[t6+F]tFg_ lFb]lFh_lFi_lFj_lF`tFb[lFb[lFRFR@$F\\`l>Fgs,&*&F^oFR-F[t6+F]tF^tF^sF`_ lFa_lFb_lF`tFazFazFRFR*&FcoFR-F[t6+F]tFg_lFb]lFh_lFi_lFj_lF`tFjzFjzFRF R@$Fg`l>Fgs,&*&F^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tFb[lFb[lFRFR*&FcoFR- F[t6+F]tFg_lFb]lFh_lFi_lFj_lF`tFazFazFRFR@$3F[vFa^lC%@$Fg^l>Fgs,&*(Fis FRF^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tFatFatFRFR*()F[pFhpFRFcoFR-F[t6+F ]tFg_lFb]lFh_lFi_lFj_lF`tF]uF]uFRFR@$F\\`l>Fgs,&*(FisFRF^oFR-F[t6+F]tF ^tF^sF`_lFa_lFb_lF`tFiuFiuFRFR*(FiclFRFcoFR-F[t6+F]tFg_lFb]lFh_lFi_lFj _lF`tFatFatFRFR@$Fg`l>Fgs,&*(FisFRF^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`tF ]uF]uFRFR*(FiclFRFcoFR-F[t6+F]tFg_lFb]lFh_lFi_lFj_lF`tFiuFiuFRFR@$555F av/FfpFcp/FjpFcpF]x@%FaelC(>F^sF^]l>Ff]lFg]l>8--%,which_funcsG6%FVFBFE >Fgs*(&FielFc_lFRF^oFR-F[t6+F]tF^tF^sF`_lFa_lFb_lF`t,&FbtFR&FielFQFR,& FbtFR&FielF]vFRFR@$F[v>Fgs*&FfoFRFgsFR@$Fa^l>FgsFhflC(>F^sFc]l>Ff]lF[^ l>FielFjel>Fgs*(F_flFRFcoFR-F[t6+F]tFg_lF^sF`_lFa_lFb_lF`tFbflFdflFR@$ F[v>Fgs*&F[pFRFgsFR@$Fa^l>FgsFeglFgsF?F?F?" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 598 "# same as gen_q_matrix_irreg_h but for biharmonic \+ splines\ngen_q_matrix_irreg := proc(T, S, m_func, vm, funcno, harm_fun cno)\n\n local i, j, q_matrix, q_val, n_verts;\n\n n_verts := 2 \+ * nops(vm) - 3;\n q_matrix := array(sparse, 1..n_verts, 1..n_verts) ;\n\n for i from 1 to n_verts do\n for j from 1 to i do\n \+ q_val := int_vs_basis12_irreg(T, S, m_func, get_basis(vm,i), \+ get_basis(vm,j), funcno, harm_funcno);\n if q_val <> 0 then \n q_matrix[i,j] := q_val; \n q_matrix[ j,i] := q_val;\n fi;\n od;\n od;\n\n q_matrix; \nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%3gen_q_matrix_irregGf*6(% \"TG%\"SG%'m_funcG%#vmG%'funcnoG%,harm_funcnoG6'%\"iG%\"jG%)q_matrixG% &q_valG%(n_vertsG6\"F3C&>8(,&-%%nopsG6#9'\"\"#\"\"$!\"\">8&-%&arrayG6% %'sparseG;\"\"\"F6FE?(8$FFFFF6%%trueG?(8%FFFFFHFIC$>8'-%5int_vs_basis1 2_irregG6)9$9%9&-%*get_basisG6$F;FH-FV6$F;FK9(9)@$0FN\"\"!C$>&F@6$FHFK FN>&F@6$FKFHFNF@F3F3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 711 " # for a q which is the analog of a square well (a triangular well in o ur case, \n# i guess), we can compute the integrals theoretically, so \+ we have separate \n# functions.\ngen_square_well_q_matrix_irreg := pro c(S, vm, add0, add1, add2, a)\n\n local i, j, q_matrix, q_val, n_ve rts;\n\n n_verts := 2 * nops(vm) - 3;\n q_matrix := array(sparse , 1..n_verts, 1..n_verts);\n\n for i from 1 to n_verts do\n \+ for j from 1 to i do\n q_val := int_square_well_vs_basis12_ irreg(S, get_basis(vm,i), get_basis(vm,j), add0, add1, add2, a);\n \+ if q_val <> 0 then \n q_matrix[i,j] := q_val; \+ \n q_matrix[j,i] := q_val;\n fi;\n od ;\n od;\n\n q_matrix;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#> %?gen_square_well_q_matrix_irregGf*6(%\"SG%#vmG%%add0G%%add1G%%add2G% \"aG6'%\"iG%\"jG%)q_matrixG%&q_valG%(n_vertsG6\"F3C&>8(,&-%%nopsG6#9% \"\"#\"\"$!\"\">8&-%&arrayG6%%'sparseG;\"\"\"F6FE?(8$FFFFF6%%trueG?(8% FFFFFHFIC$>8'-%Aint_square_well_vs_basis12_irregG6)9$-%*get_basisG6$F; FH-FT6$F;FK9&9'9(9)@$0FN\"\"!C$>&F@6$FHFKFN>&F@6$FKFHFNF@F3F3F3" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1117 "# same as harm_FEM above: \+ encapsulates all the main stuff into one function\n# good for iteratin g.\nbiharm_FEM := proc(W, level, baseno)\n local S, vm, EM, QM, EM_ QM, f_v, C;\n \n print(\"making biharms\");\n calc_biharms(T, R, level, 20):\n print(\"making basis gasket\");\n S := gasket_ irreg(W):\n print(\"making basis vertices\");\n vm := vertices_i rreg(W):\n print(\"# of basis vertices:\", nops(vm));\n print(\" making basis neighbors\");\n make_neighbors_irreg(S,W,vm):\n pri nt(\"generating energy matrix\");\n EM := gen_energy_matrix_irreg(S ,vm):\n print(\"generating q matrix\");\n QM := gen_q_matrix_irr eg(T,S,level,vm,14,20):\n print(\"adding matrices\");\n EM_QM := evalm(EM + QM):\n print(\"calculating f\");\n mult(T,level,14,8 ,16): \n add_func_lin(T,level,9,16,-1,1,11):\n print(\"generati ng f vector\");\n f_v := gen_f_vector_irreg(T,S,level,vm,11,20):\n \+ print(\"solving equation\");\n C := evalf(linsolve(EM_QM,f_v),20 ):\n print(\"making points\");\n make_points(T,level):\n prin t(\"making pretty picture\");\n make_pretty_picture_irreg(T,S,level ,vm,baseno+15,20,C):\n \nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+bi harm_FEMGf*6%%\"WG%&levelG%'basenoG6)%\"SG%#vmG%#EMG%#QMG%&EM_QMG%$f_v G%\"CG6\"F2C<-%&printG6#Q/making~biharmsF2-%-calc_biharmsG6&%\"TG%\"RG 9%\"#?-F56#Q4making~basis~gasketF2>8$-%-gasket_irregG6#9$-F56#Q6making ~basis~verticesF2>8%-%/vertices_irregGFF-F56$Q5#~of~basis~vertices:F2- %%nopsG6#FL-F56#Q7making~basis~neighborsF2-%5make_neighbors_irregG6%FC FGFL-F56#Q9generating~energy~matrixF2>8&-%8gen_energy_matrix_irregG6$F CFL-F56#Q4generating~q~matrixF2>8'-%3gen_q_matrix_irregG6(F;FCF=FL\"#9 F>-F56#Q0adding~matricesF2>8(-%&evalmG6#,&Fao\"\"\"FinF_p-F56#Q.calcul ating~fF2-%%multG6'F;F=Feo\"\")\"#;-%-add_func_linG6)F;F=\"\"*Fgp!\"\" F_p\"#6-F56#Q4generating~f~vectorF2>8)-%3gen_f_vector_irregG6(F;FCF=FL F]qF>-F56#Q1solving~equationF2>8*-%&evalfG6$-%)linsolveG6$FjoFbqF>-F56 #Q.making~pointsF2-%,make_pointsG6$F;F=-F56#Q6making~pretty~pictureF2- %:make_pretty_picture_irregG6)F;FCF=FL,&9&F_p\"#:F_pF>FjqF2F2F2" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1566 "# generates the matrices f or harmonic basis, and stores them in files\n# that we can load into s parse matrices in matlab\ngen_matrix_files_irreg_h := proc(T, S, vm, m _func, funcno)\n local EM, IP, QM, inds, fEMi, fEMj, fEMx, fIPi, fI Pj, fIPx, fQMi, fQMj, fQMx, i;\n\n print(\"generating energy matrix \");\n EM := gen_energy_matrix_irreg_h(S, vm):\n print(\"generat ing inner product matrix\");\n IP := gen_iprod_matrix_irreg_h(S, vm ):\n print(\"generating q matrix\");\n QM := gen_q_matrix_irreg_ h(T,S,m_func,vm,funcno,120):\n\n fEMi := open(\"EMi.dat\", WRITE): \n fEMj := open(\"EMj.dat\", WRITE):\n fEMx := open(\"EMx.dat\", WRITE):\n fIPi := open(\"IPi.dat\", WRITE):\n fIPj := open(\"IP j.dat\", WRITE):\n fIPx := open(\"IPx.dat\", WRITE):\n fQMi := o pen(\"QMi.dat\", WRITE):\n fQMj := open(\"QMj.dat\", WRITE):\n f QMx := open(\"QMx.dat\", WRITE):\n\n inds := [indices(EM)]:\n\n \+ for i from 1 to nops(inds) do\n fprintf(fEMi, `%20.20g\\n`, ind s[i][1]):\n fprintf(fEMj, `%20.20g\\n`, inds[i][2]):\n f printf(fEMx, `%20.20g\\n`, EM[inds[i][1], inds[i][2]]):\n fprin tf(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], in ds[i][2]]):\n fprintf(fQMi, `%20.20g\\n`, inds[i][1]):\n \+ fprintf(fQMj, `%20.20g\\n`, inds[i][2]):\n fprintf(fQMx, `%20. 20g\\n`, QM[inds[i][1], inds[i][2]]): \n od;\n\n close(fEMi):\n close(fEMj):\n close(fEMx):\n close(fIPi):\n close(fIPj): \n close(fIPx):\n close(fQMi):\n close(fQMj):\n close(fQMx ):\n\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%9gen_matrix_files_irre g_hGf*6'%\"TG%\"SG%#vmG%'m_funcG%'funcnoG60%#EMG%#IPG%#QMG%%indsG%%fEM iG%%fEMjG%%fEMxG%%fIPiG%%fIPjG%%fIPxG%%fQMiG%%fQMjG%%fQMxG%\"iG6\"F;C< -%&printG6#Q9generating~energy~matrixF;>8$-%:gen_energy_matrix_irreg_h G6$9%9&-F>6#Q@generating~inner~product~matrixF;>8%-%9gen_iprod_matrix_ irreg_hGFE-F>6#Q4generating~q~matrixF;>8&-%5gen_q_matrix_irreg_hG6(9$F F9'FG9(\"$?\">8(-%%openG6$Q(EMi.datF;%&WRITEG>8)-Fhn6$Q(EMj.datF;F[o>8 *-Fhn6$Q(EMx.datF;F[o>8+-Fhn6$Q(IPi.datF;F[o>8,-Fhn6$Q(IPj.datF;F[o>8- -Fhn6$Q(IPx.datF;F[o>8.-Fhn6$Q(QMi.datF;F[o>8/-Fhn6$Q(QMj.datF;F[o>80- Fhn6$Q(QMx.datF;F[o>8'7#-%(indicesG6#FB?(81\"\"\"F\\r-%%nopsG6#Feq%%tr ueGC+-%(fprintfG6%Ffn%)%20.20g|+G&&Feq6#F[r6#F\\r-Fcr6%F]oFer&Fgr6#\" \"#-Fcr6%FboFer&FB6$FfrF\\s-Fcr6%FgoFerFfr-Fcr6%F\\pFerF\\s-Fcr6%FapFe r&FLFbs-Fcr6%FfpFerFfr-Fcr6%F[qFerF\\s-Fcr6%F`qFer&FSFbs-%&closeG6#Ffn -Fbt6#F]o-Fbt6#Fbo-Fbt6#Fgo-Fbt6#F\\p-Fbt6#Fap-Fbt6#Ffp-Fbt6#F[q-Fbt6# F`qF;F;F;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 390 "# reads in li sts of coefficients from files (which are generated by matlab)\nread_v ect := proc(filename,n_elements)\n\n local i, fd, vals, val;\n\n \+ fd := fopen(filename, READ);\n\n #vals := readdata(fd, n_elements , float);\n\n vals := [];\n\n for i from 1 to n_elements do\n \+ val := fscanf(fd, \"%g\\n\")[1];\n vals := [op(vals), val] ;\n od;\n\n fclose(fd);\n\n vals;\n\nend;" }}{PARA 12 "" 1 " " {XPPMATH 20 "6#>%*read_vectGf*6$%)filenameG%+n_elementsG6&%\"iG%#fdG %%valsG%$valG6\"F.C'>8%-%&fopenG6$9$%%READG>8&7\"?(8$\"\"\"F<9%%%trueG C$>8'&-%'fscanfG6$F1Q$%g|+F.6#F<>F87$-%#opG6#F8FA-%'fcloseG6#F1F8F.F.F ." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 266 "# takes a list of coe fficients and normalizes them using the first one\nnormalize_coefficie nts := proc(C)\n local i, fact, temp;\n\n temp = array(1..nops(C ));\n fact := C[1];\n\n for i from 1 to nops(C) do\n temp [i] := C[i] / fact;\n od;\n\n temp;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%7normalize_coefficientsGf*6#%\"CG6%%\"iG%%factG%%temp G6\"F,C&/8&-%&arrayG6#;\"\"\"-%%nopsG6#9$>8%&F86#F4?(8$F4F4F5%%trueG>& F/6#F>*&&F8FBF4F:!\"\"F/F,F,F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 900 "# This function comes in handy when you're only examining cha nges in q, since\n# the energy and iprod matrices only depend on the b asis\ngen_matrix_files_irreg_q_only := proc(T, S, vm, m_func, funcno, \+ inds, sqflag, add0, add1, add2, a)\n local EM, IP, QM, fEMi, fEMj, \+ fEMx, fIPi, fIPj, fIPx, fQMi, fQMj, fQMx, i;\n\n print(\"generating q matrix\");\n if sqflag = true then\n QM := gen_square_wel l_q_matrix_irreg(S, vm, add0, add1, add2, a);\n else\n QM := gen_q_matrix_irreg(T,S,m_func,vm,funcno,20):\n fi;\n \n fQMi := open(\"QMi.dat\", WRITE):\n fQMj := open(\"QMj.dat\", WRITE):\n fQMx := open(\"QMx.dat\", WRITE):\n \n for i from 1 to nops( inds) do\n fprintf(fQMi, `%20.20g\\n`, inds[i][1]):\n fp rintf(fQMj, `%20.20g\\n`, inds[i][2]):\n fprintf(fQMx, `%20.20g \\n`, QM[inds[i][1], inds[i][2]]):\n od;\n\n close(fQMi):\n c lose(fQMj):\n close(fQMx):\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6 #>%>gen_matrix_files_irreg_q_onlyGf*6-%\"TG%\"SG%#vmG%'m_funcG%'funcno G%%indsG%'sqflagG%%add0G%%add1G%%add2G%\"aG6/%#EMG%#IPG%#QMG%%fEMiG%%f EMjG%%fEMxG%%fIPiG%%fIPjG%%fIPxG%%fQMiG%%fQMjG%%fQMxG%\"iG6\"F@C+-%&pr intG6#Q4generating~q~matrixF@@%/9*%%trueG>8&-%?gen_square_well_q_matri x_irregG6(9%9&9+9,9-9.>FK-%3gen_q_matrix_irregG6(9$FO9'FP9(\"#?>8--%%o penG6$Q(QMi.datF@%&WRITEG>8.-Fjn6$Q(QMj.datF@F]o>8/-Fjn6$Q(QMx.datF@F] o?(80\"\"\"Fjo-%%nopsG6#9)FIC%-%(fprintfG6%Fhn%)%20.20g|+G&&F^p6#Fio6# Fjo-Fap6%F_oFcp&Fep6#\"\"#-Fap6%FdoFcp&FK6$FdpFjp-%&closeG6#Fhn-Fbq6#F _o-Fbq6#FdoF@F@F@" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1472 "# ta kes a word list and extends it by N zooms about vert\nzoom_wordlist := proc(W, vert, N)\n local temp, w1, w2, x1, x2, w_vert, i;\n\n i f vert <> [0] and vert <> [1] and vert <> [2] and vert <> [0,1] and ve rt <> [0,2] and vert <> [1,2] and not member(vert[1..(nops(vert) - 2)] , W) then\n print(\"vertex is not a member of the subdivision\" );\n RETURN(W);\n fi;\n\n temp := W;\n\n # special cas e of boundary\n if nops(vert) = 1 then\n w1 := longest_word( W, [], vert[1], vert[1]);\n for i from 1 to N do\n w 1 := [op(w1), vert[1]];\n temp := [op(temp), w1];\n \+ od;\n else\n\n w_vert := vert[1 .. (nops(vert) - 2)];\n \+ if vert[(nops(vert) - 1) .. nops(vert)] = [0,1] then\n w 1 := longest_word(W, w_vert, 0, 1);\n w2 := longest_word(W, w_vert, 1, 0);\n x1 := 1;\n x2 := 0;\n e lif vert[(nops(vert) - 1) .. nops(vert)] = [0,2] then\n w1 \+ := longest_word(W, w_vert, 0, 2);\n w2 := longest_word(W, w _vert, 2, 0);\n x1 := 2;\n x2 := 0;\n els e\n w1 := longest_word(W, w_vert, 2, 1);\n w2 := longest_word(W, w_vert, 1, 2);\n x1 := 1;\n x2 \+ := 2;\n fi;\n \n for i from 1 to N do\n w1 := [op(w1), x1];\n w2 := [op(w2), x2];\n temp : = [op(temp), w1];\n temp := [op(temp), w2];\n od;\n \+ fi;\n\n sort(temp, word_compare);\nend; " }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.zoom_wordlistGf*6%%\"WG%%vertG%\"NG6)%%tempG%#w1G%#w 2G%#x1G%#x2G%'w_vertG%\"iG6\"F2C&@$33333309%7#\"\"!0F<7#\"\"\"0F<7#\" \"#0F<7$F>FA0F<7$F>FD0F<7$FAFD4-%'memberG6$&F<6#;FA,&-%%nopsG6#F8$FW@%/FSFAC$>8%-%-longest_wordG6&FW7\"&F<6#FAFeo?(8*FAFA9&%% trueGC$>F`o7$-%#opG6#F`oFeo>F[o7$-F_p6#F[oF`oC%>8)FO@'/&F<6#;,&FSFAFAF VFSFFC&>F`o-Fbo6&FWFgpF>FA>8&-Fbo6&FWFgpFAF>>8'FA>8(F>/FjpFHC&>F`o-Fbo 6&FWFgpF>FD>Fcq-Fbo6&FWFgpFDF>>FgqFD>FiqF>C&>F`o-Fbo6&FWFgpFDFA>Fcq-Fb o6&FWFgpFAFD>FgqFA>FiqFD?(FhoFAFAFioFjoC&>F`o7$F^pFgq>Fcq7$-F_p6#FcqFi q>F[oFbp>F[o7$FcpFcq-%%sortG6$F[o%-word_compareGF2F2F2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 196 "# simply produces a wordlist to re present level m\nreg_wordlist := proc(m)\n local out, i;\n\n out := [];\n for i from 2 to m do\n out := subdivide_irreg(out) ;\n od;\n\n out;\nend; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>% -reg_wordlistGf*6#%\"mG6$%$outG%\"iG6\"F+C%>8$7\"?(8%\"\"#\"\"\"9$%%tr ueG>F.-%0subdivide_irregG6#F.F.F+F+F+" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 379 "# returns integral of a biharmonic basis function\n# note: F-TYPE BASIS FUNCTION ONLY\nint_biharm_basis_irreg := proc(S, b asis)\n local m1, m2;\n\n # implicitly adding one to traditional values of \"m's\"\n # because formula is (1/3) * (1/3^m1 + 1/3^m2) \n m1 := max(nops(basis), nops(S[basis, 1]));\n m2 := max(nops(b asis), nops(S[basis, 4]));\n\n (1/3)^m1 + (1/3)^m2;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%7int_biharm_basis_irregGf*6$%\"SG%&basisG6 $%#m1G%#m2G6\"F,C%>8$-%$maxG6$-%%nopsG6#9%-F46#&9$6$F6\"\"\">8%-F16$F3 -F46#&F:6$F6\"\"%,&)#F<\"\"$F/F<)FHF>F " 0 "" {MPLTEXT 1 0 1024 "# constructs an array of plots based on an arr ay of coefficients\n# Contains two flags, one is to restrict the z ran ge displayed, and \n# the other is to zoom in on a point\nmake_animati on_irreg := proc(T, S, W, depth, vm, vals, funcno, harm_funcno, range_ flag, zmin, zmax, zoom_flag, pt, m)\n\n local SG_plots, final_point s, inds, i;\n\n SG_plots := [];\n\n \n for i from 1 to nops(v als) do\n\n print(`computing plot number `,i);\n make_pr etty_picture_irreg(T, S, W, depth, vm, funcno + i, harm_funcno, vals[i ]);\n if zoom_flag = true then\n final_points := zoo m_plot(T, W, depth, funcno + i, pt, m);\n else\n fin al_points := SG_plot_irreg(T, W, depth, funcno + i);\n fi;\n \+ if range_flag = true then\n SG_plots := [op(SG_plots), pointplot3d(final_points, color = BLUE, symbol = POINT, view = zmin.. zmax)];\n else\n SG_plots := [op(SG_plots), pointplo t3d(final_points, color = BLUE, symbol = POINT)];\n fi;\n od ;\n\n SG_plots;\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%5make_ animation_irregGf*60%\"TG%\"SG%\"WG%&depthG%#vmG%%valsG%'funcnoG%,harm _funcnoG%+range_flagG%%zminG%%zmaxG%*zoom_flagG%#ptG%\"mG6&%)SG_plotsG %-final_pointsG%%indsG%\"iG6\"F:C%>8$7\"?(8'\"\"\"FA-%%nopsG6#9)%%true GC&-%&printG6$%7computing~plot~number~GF@-%:make_pretty_picture_irregG 6*9$9%9&9'9(,&9*FAF@FA9+&FE6#F@@%/9/FF>8%-%*zoom_plotG6(FOFQFRFT9091>F gn-%.SG_plot_irregG6&FOFQFRFT@%/9,FF>F=7$-%#opG6#F=-%,pointplot3dG6&Fg n/%&colorG%%BLUEG/%'symbolG%&POINTG/%%viewG;9-9.>F=7$Ffo-Fjo6%FgnF\\pF _pF=F:F:F:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1944 "# Generates all the matrices for the biharmonic basis (and q integrals if qflag\n # is true) and stores them in files so we can read them into matlab. \+ \ngen_matrix_files_irreg := proc(T, S, vm, m_func, funcno, qflag, sqfl ag, add0, add1, add2, a)\n local EM, IP, QM, inds, fEMi, fEMj, fEMx , fIPi, fIPj, fIPx, fQMi, fQMj, fQMx, i;\n\n print(\"generating ene rgy matrix\");\n EM := gen_energy_matrix_irreg(S, vm):\n print( \"generating inner product matrix\");\n IP := gen_iprod_matrix_irre g(S, vm):\n if qflag = true then\n print(\"generating q matr ix\");\n if sqflag = true then\n QM := gen_square_we ll_q_matrix_irreg(S, vm, add0, add1, add2, a);\n else\n \+ QM := gen_q_matrix_irreg(T,S,m_func,vm,funcno,20):\n fi;\n \+ fi;\n\n fEMi := open(\"EMi.dat\", WRITE):\n fEMj := open(\"EM j.dat\", WRITE):\n fEMx := open(\"EMx.dat\", WRITE):\n fIPi := o pen(\"IPi.dat\", WRITE):\n fIPj := open(\"IPj.dat\", WRITE):\n f IPx := open(\"IPx.dat\", WRITE):\n if qflag = true then\n fQ Mi := open(\"QMi.dat\", WRITE):\n fQMj := open(\"QMj.dat\", WRI TE):\n fQMx := open(\"QMx.dat\", WRITE):\n fi;\n\n inds : = [indices(EM)]:\n\n for i from 1 to nops(inds) do\n fprintf (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 fprintf(fIPi, `%20.20g\\n`, inds[i][1]):\n f printf(fIPj, `%20.20g\\n`, inds[i][2]):\n fprintf(fIPx, `%20.20 g\\n`, IP[inds[i][1], inds[i][2]]):\n if qflag = true then\n \+ fprintf(fQMi, `%20.20g\\n`, inds[i][1]):\n fprintf (fQMj, `%20.20g\\n`, inds[i][2]):\n fprintf(fQMx, `%20.20g \\n`, QM[inds[i][1], inds[i][2]]):\n fi; \n od;\n\n clos e(fEMi):\n close(fEMj):\n close(fEMx):\n close(fIPi):\n cl ose(fIPj):\n close(fIPx):\n if qflag = true then\n close( fQMi):\n close(fQMj):\n close(fQMx):\n fi;\n\nend;" } }{PARA 12 "" 1 "" {XPPMATH 20 "6#>%7gen_matrix_files_irregGf*6-%\"TG% \"SG%#vmG%'m_funcG%'funcnoG%&qflagG%'sqflagG%%add0G%%add1G%%add2G%\"aG 60%#EMG%#IPG%#QMG%%indsG%%fEMiG%%fEMjG%%fEMxG%%fIPiG%%fIPjG%%fIPxG%%fQ MiG%%fQMjG%%fQMxG%\"iG6\"FAC7-%&printG6#Q9generating~energy~matrixFA>8 $-%8gen_energy_matrix_irregG6$9%9&-FD6#Q@generating~inner~product~matr ixFA>8%-%7gen_iprod_matrix_irregGFK@$/9)%%trueGC$-FD6#Q4generating~q~m atrixFA@%/9*FX>8&-%?gen_square_well_q_matrix_irregG6(FLFM9+9,9-9.>F[o- %3gen_q_matrix_irregG6(9$FL9'FM9(\"#?>8(-%%openG6$Q(EMi.datFA%&WRITEG> 8)-F^p6$Q(EMj.datFAFap>8*-F^p6$Q(EMx.datFAFap>8+-F^p6$Q(IPi.datFAFap>8 ,-F^p6$Q(IPj.datFAFap>8--F^p6$Q(IPx.datFAFap@$FVC%>8.-F^p6$Q(QMi.datFA Fap>8/-F^p6$Q(QMj.datFAFap>80-F^p6$Q(QMx.datFAFap>8'7#-%(indicesG6#FH? (81\"\"\"Fds-%%nopsG6#F]sFXC)-%(fprintfG6%F\\p%)%20.20g|+G&&F]s6#Fcs6# Fds-Fjs6%FcpF\\t&F^t6#\"\"#-Fjs6%FhpF\\t&FH6$F]tFct-Fjs6%F]qF\\tF]t-Fj s6%FbqF\\tFct-Fjs6%FgqF\\t&FRFit@$FVC%-Fjs6%F^rF\\tF]t-Fjs6%FcrF\\tFct -Fjs6%FhrF\\t&F[oFit-%&closeG6#F\\p-F[v6#Fcp-F[v6#Fhp-F[v6#F]q-F[v6#Fb q-F[v6#Fgq@$FVC%-F[v6#F^r-F[v6#Fcr-F[v6#FhrFAFAFA" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 2127 "# If we are integrating against a square w ell (triangular well?) we can compute\n# the value theoretically, sinc e q is constant. The only nontrivial part of this\n# is determining if the basis elements are within the well.\nint_square_well_vs_basis12_i rreg := proc(S, basis1, basis2, add0, add1, add2, a)\n local m, g_f lag, triangle, out;\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[basis 1[1],4] <> basis2[1] then \n RETURN(0);\n fi;\n \n #fo r our purposes we don't deal with square wells that include the bounda ry\n if nops(basis1[1]) = 1 or nops(basis2[1]) = 1 then\n RE TURN(0);\n fi;\n\n triangle := add1[1..(nops(add1) - 2)];\n t riangle := [op(triangle), add0[nops(add0)]];\n\n # if either one is exterior to the square well, return 0\n if (basis1[1] <> add0 and \+ basis1[1] <> add1 and basis1[1] <> add2 and\n basis1[1][1 .. no ps(triangle)] <> triangle) or\n (basis2[1] <> add0 and basis2[1] <> add1 and basis2[1] <> add2 and\n basis2[1][1 .. nops(triang le)] <> triangle) then\n RETURN(0);\n fi;\n\n # if both a re at one boundary vertex of square well, code special case of IP\n \+ if (basis1[1] = add0 and basis2[1] = add0) or \n (basis1[1] = a dd1 and basis2[1] = add1) or \n (basis1[1] = add2 and basis2[1] \+ = add2) then\n if S[basis1[1], 1][1 .. nops(triangle)] = triang le then\n m := max(nops(basis1[1]), nops(S[basis1[1], 1])) \+ - 1;\n g_flag := 1;\n else\n m := max(nop s(basis1[1]), nops(S[basis1[1], 4])) - 1;\n g_flag := -1;\n fi;\n if (basis1[2] = f and basis2[2] = g) or (basis1[2 ] = g and basis2[2] = f) then \n out := g_flag * 47/1395 * \+ (1/5) ^ m; \n elif basis1[2] = f and basis2[2] = f then \n \+ out := 190/837 * (1/3) ^ m; \n else \n out := 206/37665 * (3/25) ^ m; \n fi;\n # else return a * IP(basis 1, basis2) \n else\n out := a * inner_prod_irreg(S, basis1, \+ basis2);\n fi;\n\n\n out;\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%Aint_square_well_vs_basis12_irregGf*6)%\"SG%'basis1G%'basis2G%%ad d0G%%add1G%%add2G%\"aG6&%\"mG%'g_flagG%)triangleG%$outG6\"F3C)@$33330& 9%6#\"\"\"&9&F=0&9$6$F;F>F?0&FC6$F;\"\"#F?0&FC6$F;\"\"$F?0&FC6$F;\"\"% F?-%'RETURNG6#\"\"!@$5/-%%nopsG6#F;F>/-FY6#F?F>FQ>8&&9(6#;F>,&-FY6#F[o F>FH!\"\">Fin7$-%#opG6#Fin&9'6#-FY6#Fho@$53330F;Fho0F;F[o0F;9)0&F;6#;F >-FYFfoFin3330F?Fho0F?F[o0F?Fdp0&F?FgpFinFQ@%553/F;Fho/F?Fho3/F;F[o/F? F[o3/F;Fdp/F?FdpC$@%/&FBFgpFinC$>8$,&-%$maxG6$FX-FY6#FBF>F>Fao>8%F>C$> Fdr,&-Fgr6$FX-FY6#FNF>F>Fao>F\\sFao@'53/&F<6#FH%\"fG/&F@Fjs%\"gG3/FisF ^t/F]tF[t>8',$*&F\\sF>)#F>\"\"&FdrF>#\"#Z\"%&R\"3FhsFat>Fct,$)#F>FLFdr #\"$!>\"$P)>Fct,$)#FL\"#DFdr#\"$1#\"&lw$>Fct*&9*F>-%1inner_prod_irregG 6%FCFFctF3F3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 688 "# w rapper to create x and y values of basis vertices for plotting\n# actu ally computes uniform depth past basis.\nmake_points_irreg := proc(T, \+ W, depth)\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_points recursively creates the rest of \+ the coordinates out to Vm... \n\n gen_points_irreg(T, W, depth, [0] , [0], [0,1], [0,2]);\n gen_points_irreg(T, W, depth, [1], [1], [1, 2], [0,1]);\n gen_points_irreg(T, W, depth, [2], [2], [0,2], [1,2]) ; \nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%2make_points_irregGf*6%% \"TG%\"WG%&depthG6\"F*F*C1>&9$6$7#\"\"!\"\"&F1>&F.6$7#\"\"\"F2$F2!\"\" >&F.6$7#\"\"#F2$!\"&F9>&F.6$F0\"\"',$-%%sqrtG6#\"\"$#F7F>>&F.6$F6FDF1> &F.6$F=FDF1>&F.6$7$F1F7F2$\"#D!\"#>&F.6$FTFD,$FF#F7\"\"%>&F.6$7$F7F>F2 F1>&F.6$F[oFDF1>&F.6$7$F1F>F2$!#DFW>&F.6$FboFDFen-%1gen_points_irregG6 )F.9%9&F0F0FTFbo-Fio6)F.F[pF\\pF6F6F[oFT-Fio6)F.F[pF\\pF=F=FboF[oF*F*F *" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3310 "# recursively genera tes x and y values of points up to depth past basis\ngen_points_irreg \+ := proc(T, W, depth, curtriangle, add0, add1, add2)\n if member(cur triangle, W) then\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[ad d0,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_poi nts_irreg(T, W, depth, [op(curtriangle), 0], add0, [op(add1[1..nops(ad d1)-2]),0,0,1], [op(add1[1..nops(add1)-2]),0,0,2]);\n gen_p oints_irreg(T, W, depth, [op(curtriangle), 1], add1, [op(add1[1..nops( add1)-2]),0,1,2], [op(add1[1..nops(add1)-2]),0,0,1]);\n gen _points_irreg(T, W, depth, [op(curtriangle), 2], add2, [op(add1[1..nop s(add1)-2]),0,0,2], [op(add1[1..nops(add1)-2]),0,1,2]);\n fi;\n \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(add1)-2]),1,1,2],6] := .5 * T[add0,6] + .5 * T[add1,6]; \n \+ T[[op(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_irreg(T, W, d epth, [op(curtriangle), 1], add0, [op(add1[1..nops(add1)-2]),1,1,2], [ op(add1[1..nops(add1)-2]),1,0,1]);\n gen_points_irreg(T, W, depth, [op(curtriangle), 2], add1, [op(add1[1..nops(add1)-2]),1,0,2], [op(add1[1..nops(add1)-2]),1,1,2]);\n gen_points_irreg(T, \+ W, depth, [op(curtriangle), 0], add2, [op(add1[1..nops(add1)-2]),1,0,1 ], [op(add1[1..nops(add1)-2]),1,0,2]);\n fi;\n \n i f add1[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(add1)-2]),2,0,1],6] := .5 * T[add1,6] + .5 * T[add2,6]; \n T[[op(add1[1..nops(add1)-2]),2,1,2],6] := .5 * T[add0,6] + .5 * T[add2,6];\n gen_points_irreg(T, W, depth, [op(curt riangle), 2], add0, [op(add1[1..nops(add1)-2]),2,0,2], [op(add1[1..nop s(add1)-2]),2,1,2]);\n gen_points_irreg(T, W, depth, [op(cu rtriangle), 0], add1, [op(add1[1..nops(add1)-2]),2,0,1], [op(add1[1..n ops(add1)-2]),2,0,2]);\n gen_points_irreg(T, W, depth, [op( curtriangle), 1], add2, [op(add1[1..nops(add1)-2]),2,1,2], [op(add1[1. .nops(add1)-2]),2,0,1]);\n fi;\n else\n gen_points(T, depth, add0, add1, add2);\n fi; \nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%1gen_points_irregGf*6)%\"TG%\"WG%&depthG%,curtriangle G%%add0G%%add1G%%add2G6\"F.F.@%-%'memberG6$9'9%C%@$/&9)6#;,&-%%nopsG6# F9\"\"\"F@!\"\"F=7$\"\"!F@C+>&9$6$7&-%#opG6#&F96#;F@,&F=F@\"\"#FAFCFCF @\"\"&,&&FG6$9(FR$FRFA*&FWF@&FG6$F9FRF@F@>&FG6$7&FJFCF@FQFR,&FYFW*&FWF @&FG6$9*FRF@F@>&FG6$7&FJFCFCFQFR,&FTFW*&FWF@F[oF@F@>&FG6$FI\"\"',&&FG6 $FVFgoFW*&FWF@&FG6$F9FgoF@F@>&FG6$FhnFgo,&F\\pFW*&FWF@&FG6$F]oFgoF@F@> &FG6$FaoFgo,&FioFW*&FWF@FcpF@F@-F$6)FGF49&7$-FK6#F3FCFVFIFao-F$6)FGF4F \\q7$F^qF@F9FhnFI-F$6)FGF4F\\q7$F^qFQF]oFaoFhn@$/F87$F@FQC+>&FG6$7&FJF @F@FQFRFS>&FG6$7&FJF@FCFQFRFin>&FG6$7&FJF@FCF@FRFbo>&FG6$F]rFgoFho>&FG 6$FarFgoFap>&FG6$FerFgoFhp-F$6)FGF4F\\qFbqFVF]rFer-F$6)FGF4F\\qFeqF9Fa rF]r-F$6)FGF4F\\qF]qF]oFerFar@$/F87$FCFQC+>&FG6$7&FJFQFCFQFRFS>&FG6$7& FJFQFCF@FRFin>&FG6$7&FJFQF@FQFRFbo>&FG6$F\\tFgoFho>&FG6$F`tFgoFap>&FG6 $FdtFgoFhp-F$6)FGF4F\\qFeqFVF\\tFdt-F$6)FGF4F\\qF]qF9F`tF\\t-F$6)FGF4F \\qFbqF]oFdtF`t-%+gen_pointsG6'FGF\\qFVF9F]oF.F.F." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 393 "# creates our irregular plot (at a fixed d epth past the basis)\nSG_plot_irreg := proc(T, W, depth, funcno)\n \+ local i, tempW, plot_pts, func;\n\n tempW := W;\n for i from 1 t o depth do\n tempW := subdivide_irreg(tempW);\n od;\n\n f unc := vertices_irreg(tempW);\n \n plot_pts := [seq([T[func[i], 5] , T[func[i], 6], T[func[i], funcno]], i = 1..nops(func))];\n \n pl ot_pts;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%.SG_plot_irregGf *6&%\"TG%\"WG%&depthG%'funcnoG6&%\"iG%&tempWG%)plot_ptsG%%funcG6\"F0C' >8%9%?(8$\"\"\"F79&%%trueG>F3-%0subdivide_irregG6#F3>8'-%/vertices_irr egGF=>8&7#-%$seqG6$7%&9$6$&F?6#F6\"\"&&FJ6$FL\"\"'&FJ6$FL9'/F6;F7-%%no psG6#F?FCF0F0F0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 902 "# write s a series of functions (which must be spaced regularly starting at \n # func1) to a file for later retrieval\nsave_funcs := proc(T, W, depth , func1, total, step, nfname)\n\n local fd, i, j, vm, tempW, fname; \n\n fname := nfname;\n fd := fopen(fname, WRITE);\n\n fprint f(fd, \"%d\\n\", nops(W));\n for i from 1 to nops(W) do\n fp rintf(fd, \"%d\\n\", nops(W[i]));\n for j from 1 to nops(W[i]) \+ do\n fprintf(fd, \"%d\\n\", W[i][j]);\n od;\n od; \n fprintf(fd, \"%d\\n\", depth);\n fprintf(fd, \"%d\\n\", total );\n\n tempW := W;\n for i from 1 to depth do\n tempW := \+ subdivide_irreg(tempW);\n od;\n\n vm := vertices_irreg(tempW);\n for j from func1 by step to step*(total-1)+func1 do\n for i from 1 to nops(vm) do\n #print(j, i, evalf(T[vm[i], j])); \n fprintf(fd, \"%5.5g\\n\", evalf(T[vm[i], j]));\n \+ od;\n od;\n\n fclose(fd);\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+save_funcsGf*6)%\"TG%\"WG%&depthG%&func1G%&totalG%%stepG%'nfname G6(%#fdG%\"iG%\"jG%#vmG%&tempWG%&fnameG6\"F5C->8)9*>8$-%&fopenG6$F8%&W RITEG-%(fprintfG6%F;Q$%d|+F5-%%nopsG6#9%?(8%\"\"\"FJFD%%trueGC$-FA6%F; FC-FE6#&FG6#FI?(8&FJFJFOFK-FA6%F;FC&FQ6#FT-FA6%F;FC9&-FA6%F;FC9(>8(FG? (FIFJFJFenFK>Fjn-%0subdivide_irregG6#Fjn>8'-%/vertices_irregGF_o?(FT9' 9),&*&FfoFJ,&FhnFJFJ!\"\"FJFJFeoFJFK?(FIFJFJ-FE6#FaoFK-FA6%F;Q'%5.5g|+ F5-%&evalfG6#&9$6$&FaoFRFT-%'fcloseG6#F;F5F5F5" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 919 "# retrieves the functions written by save_fun cs above\nread_funcs := proc(T, funcno, nfname)\n\n local fd, i, j, vm, W, depth, tempW, len_W, len_Wi, curword, temp, total, fname;\n \+ \n fname := nfname;\n fd := fopen(fname, READ);\n\n len_W := \+ fscanf(fd, \"%d\\n\")[1];\n W := [];\n for i from 1 to len_W do \n len_Wi := fscanf(fd, \"%d\\n\")[1];\n curword := []; \n for j from 1 to len_Wi do\n temp := fscanf(fd, \" %d\\n\")[1];\n curword := [op(curword), temp];\n od; \n W := [op(W), curword];\n od;\n\n depth := fscanf(fd, \+ \"%d\\n\")[1];\n total := fscanf(fd, \"%d\\n\")[1];\n\n tempW := W;\n for i from 1 to depth do\n tempW := subdivide_irreg(te mpW);\n od;\n vm := vertices_irreg(tempW);\n\n for j from fun cno to total+funcno-1 do\n for i from 1 to nops(vm) do\n \+ T[vm[i], j] := fscanf(fd, \"%g\\n\")[1];\n od;\n od;\n \+ fclose(fd);\nend;" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%+read_funcsG f*6%%\"TG%'funcnoG%'nfnameG6/%#fdG%\"iG%\"jG%#vmG%\"WG%&depthG%&tempWG %&len_WG%'len_WiG%(curwordG%%tempG%&totalG%&fnameG6\"F8C.>809&>8$-%&fo penG6$F;%%READG>8+&-%'fscanfG6$F>Q$%d|+F86#\"\"\">8(7\"?(8%FKFKFD%%tru eGC&>8,FE>8-FN?(8&FKFKFTFQC$>8.FE>FV7$-%#opG6#FVFen>FM7$-Fin6#FMFV>8)F E>8/FE>8*FM?(FPFKFKF`oFQ>Fdo-%0subdivide_irregG6#Fdo>8'-%/vertices_irr egGFio?(FX9%FK,(FboFKF_pFKFK!\"\"FQ?(FPFKFK-%%nopsG6#F[pFQ>&9$6$&F[p6# FPFX&-FG6$F>Q$%g|+F8FJ-%'fcloseG6#F>F8F8F8" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 1042 "# creates a plot along a line of the gasket.\n# t he line is defined by a starting point, a neighbor direction to go the first\n# time, and a neighbor direction to go all subsequent times. A lways plots along\n# the line until the line ends.\nline_plot_irreg := proc(T, S, startpt, first_neighbor, rest_neighbor, funcno)\n local line_pts, plot_pts, ending, last_pt, numpts;\n\n line_pts := [star tpt, S[startpt, first_neighbor]];\n\n ending := last_digits(line_pt s[2]);\n last_pt := line_pts[2];\n\n while last_digits(last_pt) \+ = ending do\n line_pts := [op(line_pts), S[last_pt,rest_neighbo r]];\n last_pt := S[last_pt, rest_neighbor];\n if nops(l ast_pt) = 1 then\n last_pt := [op(last_pt), op(last_pt)]\n \+ fi;\n od;\n\n numpts := nops(line_pts);\n\n #use dista nce from first line point as x coordinate, since irregular\n plot_p ts := [seq([sqrt(( T[line_pts[i], 5] - T[line_pts[1], 5] ) ^ 2 + ( T[l ine_pts[i], 6] - T[line_pts[1], 6] ) ^ 2), T[line_pts[i], funcno]], i \+ = 1..numpts)];\n\n\n plot_pts;\nend;\n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%0line_plot_irregGf*6(%\"TG%\"SG%(startptG%/first_neig hborG%.rest_neighborG%'funcnoG6'%)line_ptsG%)plot_ptsG%'endingG%(last_ ptG%'numptsG6\"F3C)>8$7$9&&9%6$F89'>8&-%,last_digitsG6#&F66#\"\"#>8'FB ?(F3\"\"\"FHF3/-F@6#FFF>C%>F67$-%#opG6#F6&F:6$FF9(>FFFR@$/-%%nopsGFKFH >FF7$-FPFKFfn>8(-FYFQ>8%7#-%$seqG6$7$-%%sqrtG6#,&*$),&&9$6$&F66#%\"iG \"\"&FH&Fio6$&F66#FHF^p!\"\"FDFHFH*$),&&Fio6$F[p\"\"'FH&Fio6$FapFipFcp FDFHFH&Fio6$F[p9)/F]p;FHFhnF[oF3F3F3" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1266 "# creates a list of points that represent the funct ion in funcno restricted\n# to the two level m triangles around pt\n# \+ CAUTION: not intended for boundary vertices - although that's a simple function\n# to write (just restrict to first m digits are 0 \+ / 1 / 2\nzoom_plot := proc(T, W, depth, funcno, pt, m)\n local i, t empW, S, vm, prefix, ending, beg1, beg2, plot_pts, point_list;\n\n \+ # too computationally intensive for high m\n # I'm currently too la zy to hardcode finding just the four neighbors\n #S := gasket(m);\n #vm := vertices(m);\n #make_neighbors(S, m, vm);\n\n tempW : = W;\n for i from 1 to depth do\n tempW := subdivide_irreg(t empW);\n od;\n vm := vertices_irreg(tempW);\n\n prefix := pt[ 1..nops(pt)-2];\n ending := pt[nops(pt)-1..nops(pt)];\n\n beg1 : = [op(prefix), ending[1], seq(ending[2], i=1..m-nops(prefix)-1)];\n \+ beg2 := [op(prefix), ending[2], seq(ending[1], i=1..m-nops(prefix)-1) ];\n\n point_list := select(begins_with, vm, beg1);\n point_list := [op(point_list), op(select(begins_with, vm, beg2))];\n #point_l ist := [op(point_list), seq(S[pt, i], i=1..4)];\n \n plot_pts := [ seq([T[point_list[i], 5], T[point_list[i], 6], T[point_list[i], funcno ]], i = 1..nops(point_list))];\n \n plot_pts;\nend; \n" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#>%*zoom_plotGf*6(%\"TG%\"WG%&depthG%'funcnoG %#ptG%\"mG6,%\"iG%&tempWG%\"SG%#vmG%'prefixG%'endingG%%beg1G%%beg2G%)p lot_ptsG%+point_listG6\"F8C->8%9%?(8$\"\"\"F?9&%%trueG>F;-%0subdivide_ irregG6#F;>8'-%/vertices_irregGFE>8(&9(6#;F?,&-%%nopsG6#FMF?\"\"#!\"\" >8)&FM6#;,&FQF?F?FUFQ>8*7%-%#opG6#FK&FW6#F?-%$seqG6$&FW6#FT/F>;F?,(9)F ?-FRF[oFUF?FU>8+7%FinFao-F_o6$F\\oFco>8--%'selectG6%%,begins_withGFGFg n>F^p7$-Fjn6#F^p-Fjn6#-F`p6%FbpFGFio>8,7#-F_o6$7%&9$6$&F^p6#F>\"\"&&Fb q6$Fdq\"\"'&Fbq6$Fdq9'/F>;F?-FRFfpF\\qF8F8F8" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 184 "# just checks if l begins with beg - used in zo om_plot\nbegins_with := proc(l, beg)\n if nops(l) < nops(beg) then \n RETURN(false);\n fi;\n evalb(l[1..nops(beg)] = beg);\n end;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,begins_withGf*6$%\"lG%$begG 6\"F)F)C$@$2-%%nopsG6#9$-F.6#9%-%'RETURNG6#%&falseG-%&evalbG6#/&F06#; \"\"\"F1F3F)F)F)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "45 0 0" 284 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }