#C24.txt, April 21, 2014, starting the Finite Element Method Help:=proc(): print(` ET1(T,A,x,y), Le(A,B,x,y)`): print(` InsT(T,x,y), ET(T,A,x,y), Tg(N) , TestF(T,c,x,y)`): end: #ET1(T,A,x,y): inputs a triangle T #(given as a set of three vertices), and #A, a member of T, and variables x,y #outputs the (uniqe) linear expression a+b*x+c*y #that is 1 at A and 0 at the other two vertices ET1:=proc(T,A,x,y) local L,a,b,c,var,eq,T1: L:=a+b*x+c*y: if not member(A,T) then RETURN(FAIL): fi: var:={a,b,c}: T1:= T minus {A}: eq:={subs({x=A[1],y=A[2]},L)-1, seq(subs({x=T1[i][1],y=T1[i][2]},L),i=1..nops(T1))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): else subs(var,L): fi: end: #Le(A,B,x,y): inputs two points A and B and #outputs the line through them Le:=proc(A,B,x,y) local m: if B[1]=A[1] then RETURN(x-A[1]): fi: m:=(B[2]-A[2])/(B[1]-A[1]): expand((y-B[2])-m*(x-B[1])): end: #InsT(T,x,y): the inside of the triangle T InsT:=proc(T,x,y) local C,L1,L2,L3: C:=convert(T,`+`)/3: L1:=Le(T[1],T[2],x,y): L2:=Le(T[2],T[3],x,y): L3:=Le(T[3],T[1],x,y): expand(sign(subs({x=C[1],y=C[2]},L1))*L1)>0 and expand(sign(subs({x=C[1],y=C[2]},L2))*L2)>0 and expand(sign(subs({x=C[1],y=C[2]},L3))*L3)>0: end: #ET(T,A,x,y): inputs a triangle T #(given as a set of three vertices), and #A, a member of T, and variables x,y #outputs the (uniqe) piecewise expression #that equals something like a+b*x+c*y #INSIDE the triangle, such #that is 1 at A and 0 at the other two vertices #and 0 outside ET:=proc(T,A,x,y) local f,c: f:=ET1(T,A,x,y): c:=InsT(T,x,y): piecewise(c,f, not c, 0): end: #Tg(N): tirangulates into right-angled triangles #of edge-side 1 the square [0,N]x[0,N] Tg:=proc(N) local i,j: #{seq(seq({[i,j],[i+1,j],[i,j+1],[i+1,j+1]},i=0..N-1),j=0..N-1)}: {seq(seq({[i,j],[i+1,j],[i,j+1]},i=0..N-1),j=0..N-1)} union {seq(seq({[i+1,j],[i,j+1],[i+1,j+1]},i=0..N-1),j=0..N-1)}: end: #TestF(T,c,x,y): the generic linear combination #of elements that defines a piece-wise continuous #function on the domain that is the union of the #insides of the triangles of the set T TestF:=proc(T,c,x,y) local V,t,F,v,var: V:={seq(op(t), t in T)}: var:={seq(c[v],v in V)}: F:=0: for t in T do for v in t do F:=F+c[v]*ET(t,v,x,y): od: od: F,var: end: