#C25.txt, April 24, 2014, continuing the Finite Element #method Help:=proc(): print(`TestF(T,c,x,y), IntTr(f,x,y,t)`): print(` Iw1(w,p,q,r,f,x,y,t) `): print(` Iw(F,p,q,r,f,x,y), FEM(p,q,r,f,g,x,y,N), EvalF1(F,x,y,pt)`): end: ####old stuff #C24.txt, April 21, 2014, starting the Finite Element Method Help24:=proc(): print(` ET1(T,A,x,y), Le(A,B,x,y)`): print(` InsT(T,x,y), ET(T,A,x,y), Tg(N) , TestFstupid(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,i: 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: #TestFstupid(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 TestFstupid:=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: ###end old stuff #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,i: V:={seq(op(t), t in T)}: var:={seq(c[v],v in V)}: F:=[]: for t in T do F:=[op(F),[t, add(c[t[i]]* ET1(t,t[i],x,y),i=1..3)]]: od: F,var: end: #IntTr(f,x,y,t): double-integral of f #over the interior of the triangle t IntTr:=proc(f,x,y,t) local i: VectorCalculus[int](f,[x,y]=Triangle(seq(`<,>`(op(t[i])),i=1..3))); end: #Iw1(w,p,q,r,f,x,y,t): inputs an expression #w, and expressions p,q,r,f in the variables #x,y, and a triangle t, finds I[w] as given #in (12.44) of the book Iw1:=proc(w,p,q,r,f,x,y,t) local w1: w1:=1/2*(p*diff(w,x)^2+ q*diff(w,y)^2-r*w^2)+ f*w: expand(IntTr(w1,x,y,t)): end: #Iw(F,p,q,r,f,x,y): inputs a piece-wise function #F (given a list of pairs [t,ExpIn(x,y)] # and expressions p,q,r,f in the variables #x,y, and a triangle t, finds I[w] as given #in (12.44) of the book Iw:=proc(F,p,q,r,f,x,y) local i: #region:=F[i][1]: expr:=F[i][2] expand(add(Iw1(F[i][2],p,q,r,f,x,y,F[i][1]),i=1..nops(F))): end: #FEM(p,q,r,f,g,x,y,N): the finite-element apprx. solution #to the elliptic PDE Dirichlet BVP #(12.41), u=g on the boundary using the #triang. Tg(N) FEM:=proc(p,q,r,f,g,x,y,N) local w,c,var,B,T,b,N1244,v1,eq,i: T:=Tg(N): B:={seq([0,i],i=0..N),seq([i,0],i=0..N), seq([N,i],i=0..N),seq([i,N],i=0..N)}: w:=TestF(T,c,x,y): var:=w[2]: w:=w[1]: for b in B do w:=subs(c[b]=subs({x=b[1],y=b[2]},g),w): var:=var minus {c[b]}: od: if var={} then RETURN(w): fi: N1244:=Iw(w,p,q,r,f,x,y): eq:={seq(diff(N1244,v1), v1 in var)}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,w): end: #EvalF1(F,x,y,pt): evaluate the point pt #in the piecewise function F EvalF1:=proc(F,x,y,pt) local i: for i from 1 to nops(F) do if evalb(subs({x=pt[1],y=pt[2]},InsT(F[i][1],x,y))) then RETURN(subs({x=pt[1],y=pt[2]},F[i][2])): fi: od: FAIL: end: