restart:# Finite Element Algorithm 12.5## To approximate the solution to an elliptic partial-differential# equation subject to Dirichlet, mixed, or Neumann boundary# conditions:## Input: see STEP 0## Output: description of triangles, nodes, line integrals, basis# functions, linear system to be solved, and the# coefficients of the basis functions### Step 0# General outline## 1. Triangles numbered: 1 to K for triangles with no edges on# Script-S-1 or Script-S-2, K+1 to N for triangles with# edges on script-S-2, N+1 to M for remaining triangles.# Note: K=0 implies that no triangle is interior to D.# Note: M=N implies that all triangles have edges on# script-S-2.## 2. Nodes numbered: 1 to LN for interior nodes and nodes on# script-S-2, LN+1 to LM for nodes on script-S-1.# Note: LM and LN represent lower case m and n resp.# Note: LN=LM implies that script-S-1 contains no nodes.# Note: If a node is on both script-S-1 and script-S-2, # it should be treated as if it were only on script-S-1.## 3. NL=number of line segments on script-S-2# line(I,J) is an NL by 2 array where# line(I,1)= first node on line I and# line(I,2)= second node on line I taken# in positive direction along line I## 4. For the node labelled KK,KK=1,...,LM we have:# A) coordinates XX(KK),YY(KK)# B) number of triangles in which KK is a vertex= LL(KK)# C) II(KK,J) labels the triangles KK is in and# NV(KK,J) labels which vertex node KK is for# each J=1,...,LL(KK)## 5. NTR is an M by 3 array where# NTR(I,1)=node number of vertex 1 in triangle I# NTR(I,2)=node number of vertex 2 in triangle I# NTR(I,3)=node number of vertex 3 in triangle I## 6. Function subprograms:# A) P,Q,R,F,G,G1,G2 are the functions specified by# the particular differential equation# B) RR is the integrand# R*N(J)*N(K) on triangle I in step 4# C) FFF is the integrand F*N(J) on triangle I in step 4# D) GG1=G1*N(J)*N(K)# GG2=G2*N(J)# GG3=G2*N(K)# GG4=G1*N(J)*N(J)# GG5=G1*N(K)*N(K)# integrands in step 5# E) QQ(FF) computes the double integral of any# integrand FF over a triangle with vertices given by# nodes J1,J2,J3 - the method is an O(H**2) approximation# for triangles# F) SQ(PP) computes the line integral of any integrand PP# along the line from (XX(J1),YY(J1)) to (XX(J2),YY(J2))# by using a parameterization given by:# X=XX(J1)+(XX(J2)-XX(J1))*T# Y=YY(J1)+(YY(J2)-YY(J1))*T# for 0 <= t <= 1# and applying Simpson's composite method with H=.01## 7. Arrays:# A) A,B,C are M by 3 arrays where the basis function N# for the Ith triangle, Jth vertex is# N(X,Y)=A(I,J)+B(I,J)*X+C(I,J)*Y# for J=1,2,3 and I=1,2,...,M# B) XX,YY are LM by 1 arrays to hold coordinates of nodes# C) line,LL,II,NV,NTR have been explained above# D) GAMMA1 and Alpha are clear## 8. Note that A,B,C,XX,YY,I,I1,I2,J1,J2,J3,Delta are reserved# storage so that in any subprogram we know that# triangle I has vertices (XX(J1),YY(J1)),(XX(J2),YY(J2)),# (XX(J3),YY(J3)). That is, vertex 1 is node J1, vertex 2 # is node J2, vertex 3 is node J3 unless a line integral is# involved in which case the line integral goes from node J1# to node J2 in triangle I or unless vertex I1 is node J1# and vertex I2 is node J2 - the basis functions involve# A(I,I1)+B(I,I1)*X+C(I,I1)*Y for vertex I1 in triangle I# and A(I,I2)+B(I,I2)*X+C(I,I2)*Y for vertex I2 in triangle# I delta is 1/2 the area of triangle I## To change problems:# 1) change function subprograms P,Q,R,F,G,G1,G2# 2) change data input for K,N,M,LN,LM,NL.# 3) change data input for XX,YY,LLL,II,NV.# 4) change data input for line.global P,Q,R,F,G,G1,G2,G3,G4,G5,RR,FFF,GG1,GG2,GG3,GG4,GG5,QQ,SQ,S:# The following code is commented out. It could be used to include# the functions P, Q, R, F, G, G1 as part of the code.#P := proc(X,Y) local p:#p := 1:#RETURN(p):#Q := proc(X,Y) local q:#q := 1:#RETURN(q):#R := proc(X,Y) local r:#r := 0:#RETURN(r):#F := proc(X,Y) local f:#f := 0:#RETURN(f):#G := proc(X,Y) local g:#g := 4:#RETURN(g):#end:#G1 := proc(X,Y) local g1:#g1 := 0:#RETURN(g1):#end:# Function G2 is defined in code for this example instead of inputG2 := proc(X,Y) local T, Z1, g2:T := 1.0E-05:Z1 := 0:if 0.2-T <= X and X <= 0.4+T and abs(Y-0.2) <= T thenZ1 := X:fi:if 0.5-T <= X and X <= (0.6+T) and abs(Y-0.1) <= T thenZ1 := X:fi:if -T <= Y and Y <= 0.1+T and abs(X-0.6) <= T thenZ1 := Y:fi:if -T <= X and X <= 0.2+T and abs(Y+X-0.4) <= T thenZ1 := (X+Y)/sqrt(2):fi:if 0.4 -T <= X and X <= 0.5+T and abs(Y+X-0.6) <= T thenZ1 := (X+Y)/sqrt(2):fi:g2 := Z1:RETURN(g2):end:RR := proc(X,Y,A,B,C,I3,I1,I2) local rr:rr := R(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):RETURN(rr):end:FFF := proc(X,Y,A,B,C,I3,I1) local fff:fff := F(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):RETURN(fff):end:GG1 := proc(X,Y,A,B,C,I3,I1,I2) local gg1:gg1 := G1(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):RETURN(gg1):end:GG2 := proc(X,Y,A,B,C,I3,I1) local gg2:gg2 := G2(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):RETURN(gg2):end:GG3 := proc(X,Y,A,B,C,I3,I2) local gg3:gg3 := G2(X,Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):RETURN(gg3):end:GG4 := proc(X,Y,A,B,C,I3,I1) local gg4:gg4 := G1(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):RETURN(gg4):end:GG5 := proc(X,Y,A,B,C,I3,I2) local gg5:gg5 := G1(X,Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):RETURN(gg5):end:QQ := proc(L,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C) local X, Y, I3, T1, T2, T3, QQQ, qq: global S:X[0] := XX[J1-1]:Y[0] := YY[J1-1]:X[1] := XX[J2-1]:Y[1] := YY[J2-1]:X[2] := XX[J3-1]:Y[2] := YY[J3-1]:X[3] := 0.5*(X[0]+X[1]):Y[3] := 0.5*(Y[0]+Y[1]):X[4] := 0.5*(X[0]+X[2]):Y[4] := 0.5*(Y[0]+Y[2]):X[5] := 0.5*(X[1]+X[2]):Y[5] := 0.5*(Y[1]+Y[2]):X[6] := (X[0]+X[1]+X[2])/3.0:Y[6] := (Y[0]+Y[1]+Y[2])/3.0:if L = 1 thenfor I3 from 1 to 7 doS[I3-1] := P(X[I3-1],Y[I3-1]):od: fi:if L = 2 thenfor I3 from 1 to 7 doS[I3-1] := Q(X[I3-1],Y[I3-1]):od: fi:if L = 3 thenfor I3 from 1 to 7 doS[I3-1] := RR(X[I3-1],Y[I3-1],A,B,C,I3,I1,I2):od: fi:if L = 4 thenfor I3 from 1 to 7 doS[I3-1] := FFF(X[I3-1],Y[I3-1],A,B,C,I3,I1):od: fi:T1 := 0:for I3 from 1 to 3 doT1 := T1+S[I3-1]:od:T2 := 0:for I3 from 4 to 6 doT2 := T2+S[I3-1]:od:T3 := S[6]:QQQ := 0.5*(T1/20+2*T2/15+9*T3/20)*abs(DELTA):qq := QQQ:RETURN(qq):end:SQ := proc(L,XX,YY,J1,J2,I3,I1,I2,A,B,C) local X1, Y1, X2, Y2, H, T1, T2, T3, K, X, Q3, Q1, Q2, SSQ, sq: global S:X1 := XX[J1-1]:Y1 := YY[J1-1]:X2 := XX[J2-1]:Y2 := YY[J2-1]:H := 0.01:T1 := X2-X1:T2 := Y2-Y1:T3 := sqrt(T1*T1+T2*T2):for K from 1 to 101 doX := (K-1)*H:if L = 1 thenS[K-1] := T3*GG1(T1*X+X1,T2*X+Y1,A,B,C,I3,I1,I2):fi:if L = 2 thenS[K-1] := T3*GG2(T1*X+X1,T2*X+Y1,A,B,C,I3,I1):fi:if L = 3 thenS[K-1] := T3*GG3(T1*X+X1,T2*X+Y1,A,B,C,I3,I2):fi:if L = 4 thenS[K-1] := T3*GG4(T1*X+X1,T2*X+Y1,A,B,C,I3,I1):fi:if L = 5 thenS[K-1] := T3*GG5(T1*X+X1,T2*X+Y1,A,B,C,I3,I2):fi:od:Q3 := S[0]+S[100]:Q1 := 0:Q2 := S[99]:for K from 1 to 49 doQ1 := Q1+S[2*K]:Q2 := Q2+S[2*K-1]:od:SSQ := (Q3+2*Q1+4*Q2)*H/3:sq := SSQ:RETURN(sq):end:print(`This is the Finite Element Method.`):OK := FALSE:print(`The input will be from a text file in the following form:`):print(`K N M n m nl`):print(`Each of the above is an integer -`):print(`separate with at least one blank.`):print(`Follow with the input for each node in the form:`):print(`x-coord., y-coord., number of triangles in which the`):print(` node is a vertex.`):print(`Continue with the triangle number and vertex number for`):print(`each triangle in which the node is a vertex.`):print(`Separate each entry with at least one blank.`):print(`After all nodes have been entered follow with information`):print(`on the lines over which line integrals must be computed.`):print(`The format of this data will be the node number of the`):print(`starting node, followed by the node number of the ending`):print(`node for each line, taken in the positive direction.`):print(`There should be 2 * nl such entries, each an integer`):print(`separated by a blank.`):print(`Functions can be input or coded as procedures.`):print(`The example has G2 as a procedure and the rest`):print(`as input.`):print(`Are the functions P,Q,R,F,G,G1,G2 ready for input or`):print(`have they been included in code and has the input `):print(`file been created? Answer 1 for yes or 2 for no.`):AA := scanf(`%d`)[1]: print(`Input is `): print(AA):if AA = 1 thenprint(`Input function P(X,Y) in terms of x and y`):P := scanf(`%a`)[1]: print(`P(x,y) = `): print(P):P := unapply(P,x,y):print(`Input function Q(X,Y) in terms of x and y`):Q := scanf(`%a`)[1]: print(`Q(x,y) = `): print(Q):Q := unapply(Q,x,y):print(`Input function R(X,Y) in terms of x and y`):R := scanf(`%a`)[1]: print(`R(x,y) = `): print(R):R := unapply(R,x,y):print(`Input function F(X,Y) in terms of x and y`):F := scanf(`%a`)[1]: print(`F(x,y) = `): print(F):F := unapply(F,x,y):print(`Input function G(X,Y) in terms of x and y`):G := scanf(`%a`)[1]: print(`G(x,y) = `): print(G):G := unapply(G,x):print(`Input function G1(X,Y) in terms of x and y`):G1 := scanf(`%a`)[1]: print(`G1(x,y) = `): print(G1):G1 := unapply(G1,x):#print(`Input function G2(X) in terms of x`):#G2 := scanf(`%a`)[1]: print(`G2(x,y) = `): print(G2):#G2 := unapply(G2,x):print(`Input the file name in the form - drive:\134\134name.ext`):print(`for example: A:\134\134DATA.DTA`):NAME := scanf(`%s`)[1]: print(`Input file is `): print(NAME):INP := fopen(NAME,READ,TEXT):OK := TRUE:K := fscanf(INP, `%d`)[1]:N := fscanf(INP, `%d`)[1]:M := fscanf(INP, `%d`)[1]:LN1 := fscanf(INP, `%d`)[1]:LM := fscanf(INP, `%d`)[1]:NL := fscanf(INP, `%d`)[1]:for KK from 1 to LM doXX[KK-1] := fscanf(INP, `%f`)[1]:YY[KK-1] := fscanf(INP, `%f`)[1]:LLL := fscanf(INP, `%d`)[1]:for J from 1 to LLL doII[KK-1,J-1] := fscanf(INP, `%d`)[1]:NV[KK-1,J-1] := fscanf(INP, `%d`)[1]:od:LL[KK-1] := LLL:for J from 1 to LLL doN1 := II[KK-1,J-1]:N2 := NV[KK-1,J-1]:NTR[N1-1,N2-1] := KK:od:od:if NL > 0 thenfor I3 from 1 to NL dofor J from 1 to 2 doLINE1[I3-1,J-1] := fscanf(INP, `%d`)[1]:od:od:fi:fclose(INP):elseprint(`The program will end so that the input file `):print(`can be created and the functions defined.`):fi:if OK = TRUE thenprint(`Choice of output method:`):print(`1. Output to screen`):print(`2. Output to text file`):print(`Please enter 1 or 2.`):FLAG := scanf(`%d`)[1]: print(`Input is `): print(FLAG):if FLAG = 2 thenprint(`Input the file name in the form - drive:\134\134name.ext`):print(`for example: A:\134\134OUTPUT.DTA`):NAME := scanf(`%s`)[1]: print(`Output file is `): print(NAME):OUP := fopen(NAME,WRITE,TEXT):elseOUP := default:fi:fprintf(OUP, `FINITE ELEMENT METHOD\134n\134n\134n`):K1 := K+1:N1 := LN1+1:fprintf(OUP, `Vertices and Nodes of Triangles\134n`):fprintf(OUP, `Triangle-node number for vertex 1 to 3\134n`):for I3 from 1 to M dofprintf(OUP, ` %5d`, I3):for J from 1 to 3 dofprintf(OUP, ` %4d`, NTR[I3-1,J-1]):od:fprintf(OUP, `\134n`):od:fprintf(OUP, `x and y coordinates of nodes\134n`):for KK from 1 to LM dofprintf(OUP, `%5d %11.8f %11.8f\134n`, KK, XX[KK-1], YY[KK-1]):od:fprintf(OUP, `Lines of the Domain\134n`):for KK from 1 to NL dofprintf(OUP, `%5d %4d %4d\134n`, KK, LINE1[KK-1,0], LINE1[KK-1,1]):od:# Step 1if LM <> LN1 thenfor L from N1 to LM doGAMMA1[L-1] := G(XX[L-1],YY[L-1]):od:fi:# Step 2 - initialization of ALPHA# Note that ALHPA[I,LN1+1] = BETA[I]for I3 from 1 to LN1 dofor J from 1 to N1 doALPHA[I3-1,J-1] := 0:od:# Steps 3, 4 and 6 - 12 are within the next loop# For each triangle I let node J1 be vertex 1, node J2 be vertex 2,# and node J3 be vertex 3# Step 3od:for I3 from 1 to M doJ1 := NTR[I3-1,0]:J2 := NTR[I3-1,1]:J3 := NTR[I3-1,2]:DELTA := XX[J2-1]*YY[J3-1]-XX[J3-1]*YY[J2-1]-XX[J1-1]*(YY[J3-1]-YY[J2-1])+YY[J1-1]*(XX[J3-1]-XX[J2-1]):A[I3-1,0] := (XX[J2-1]*YY[J3-1]-YY[J2-1]*XX[J3-1])/DELTA:B[I3-1,0] := (YY[J2-1]-YY[J3-1])/DELTA:C[I3-1,0] := (XX[J3-1]-XX[J2-1])/DELTA:A[I3-1,1] := (XX[J3-1]*YY[J1-1]-YY[J3-1]*XX[J1-1])/DELTA:B[I3-1,1] := (YY[J3-1]-YY[J1-1])/DELTA:C[I3-1,1] := (XX[J1-1]-XX[J3-1])/DELTA:A[I3-1,2] := (XX[J1-1]*YY[J2-1]-YY[J1-1]*XX[J2-1])/DELTA:B[I3-1,2] := (YY[J1-1]-YY[J2-1])/DELTA:C[I3-1,2] := (XX[J2-1]-XX[J1-1])/DELTA:# Step 4# I1 = J for Step 4 and I1 = K for Step 7for I1 from 1 to 3 do# Step 8JJ1 := NTR[I3-1,I1-1]:# I2 = K for Step 4 and I2 = J for Step 9for I2 from 1 to I1 do# Steps 4 and 10JJ2 := NTR[I3-1,I2-1]:ZZ := B[I3-1,I1-1]*B[I3-1,I2-1]*QQ(1,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)+C[I3-1,I1-1]*C[I3-1,I2-1]*QQ(2,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)-QQ(3,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C):# Steps 11 and 12if JJ1 <= LN1 thenif JJ2 <= LN1 thenALPHA[JJ1-1,JJ2-1] := ALPHA[JJ1-1,JJ2-1]+ZZ:if JJ1 <> JJ2 thenALPHA[JJ2-1,JJ1-1] := ALPHA[JJ2-1,JJ1-1]+ZZ:fi:elseALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]-ZZ*GAMMA1[JJ2-1]:fi:elseif JJ2 <= LN1 thenALPHA[JJ2-1,N1-1] := ALPHA[JJ2-1,N1-1]-ZZ*GAMMA1[JJ1-1]:fi:fi:od:HH := -QQ(4,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C):if JJ1 <= LN1 thenALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]+HH:fi:od:od:# Output the basis functionsfprintf(OUP, `Basis Functions\134n`):fprintf(OUP, `Triangle - Vertex - Node - Function\134n`):for I3 from 1 to M dofor J from 1 to 3 dofprintf(OUP, `%4d %3d %3d %13.8f %13.8f %13.8f\134n`, I3, J, NTR[I3-1,J-1], A[I3-1,J-1],B[I3-1,J-1], C[I3-1,J-1]):od:od:# Step 5# For each line segment JI = 1, ..., NL and for each triangle I,# I = K1, ..., N which may contain line JI search all 3# vertices for possible correspondences.# Step 5 and Steps 13 - 19if NL <> 0 and N <> K thenfor JI from 1 to NL dofor I3 from K1 to N dofor I1 from 1 to 3 do# I1 = J in Step 5 and I1 = K in Step 14# Step 15J1 := NTR[I3-1,I1-1]:if LINE1[JI-1,0] = J1 thenfor I2 from 1 to 3 do# I2 = K in Step 5 and I2 = J in Step 16# Step 17J2 := NTR[I3-1,I2-1]:if LINE1[JI-1,1] = J2 then# There is a correspondence of vertix I1 in triangle I with node J1# as the start of line JI and vertex I2 with node J2 as the# end of line JI# Step 5XJ := SQ(1,XX,YY,J1,J2,I3,I1,I2,A,B,C):XJ1 := SQ(4,XX,YY,J1,J2,I3,I1,I2,A,B,C):XJ2 := SQ(5,XX,YY,J1,J2,I3,I1,I2,A,B,C):XI1 := SQ(2,XX,YY,J1,J2,I3,I1,I2,A,B,C):XI2 := SQ(3,XX,YY,J1,J2,I3,I1,I2,A,B,C):# Steps 8 and 19if J1 <= LN1 thenif J2 <= LN1 thenALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1:ALPHA[J1-1,J2-1] := ALPHA[J1-1,J2-1]+XJ:ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2:ALPHA[J2-1,J1-1] := ALPHA[J2-1,J1-1]+XJ:ALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]+XI1:ALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]+XI2:elseALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]-XJ*GAMMA1[J2-1]+XI1:ALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1:fi:elseif J2 <= LN1 thenALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]-XJ*GAMMA1[J1-1]+XI2:ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2:fi:fi:fi:od:fi:od:od:od:fi:# Output ALPHAfprintf(OUP, `Matrix ALPHA follows:\134n`):for I3 from 1 to LN1 dofprintf(OUP, `Row %4d\134n`, I3):for J from 1 to N1 dofprintf(OUP, ` %13.10e\134n `, evalf(ALPHA[I3-1,J-1])):od:od:fprintf(OUP, `\134n`):# Step 20if LN1 > 1 thenINN := LN1-1:for I3 from 1 to INN doI1 := I3+1:for J from I1 to LN1 doCCC := ALPHA[J-1,I3-1]/ALPHA[I3-1,I3-1]:for J1 from I1 to N1 doALPHA[J-1,J1-1] := ALPHA[J-1,J1-1]-CCC*ALPHA[I3-1,J1-1]:od:ALPHA[J-1,I3-1] := 0:od:od:fi:GAMMA1[LN1-1] := ALPHA[LN1-1,N1-1]/ALPHA[LN1-1,LN1-1]:if LN1 > 1 thenfor I3 from 1 to INN doJ := LN1-I3:CCC := ALPHA[J-1,N1-1]:J1 := J+1:for KK from J1 to LN1 doCCC := CCC-ALPHA[J-1,KK-1]*GAMMA1[KK-1]:od:GAMMA1[J-1] := CCC/ALPHA[J-1,J-1]:od:fi:# Step 21# Output GAMMA1fprintf(OUP, `Coefficients of Basis Functions:\134n`):for I3 from 1 to LM doLLL := LL[I3-1]:fprintf(OUP, `%3d %13.8f %d`, I3, evalf(GAMMA1[I3-1]), I3):for J from 1 to LLL dofprintf(OUP, ` %4d`, II[I3-1,J-1]):od:fprintf(OUP, `\134n`):od:if OUP <> default thenfclose(OUP):print(`Output file `,NAME,` created successfully`):fi:# Step 23fi:JSFH