alg125.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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) Gamma 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.
> alg125 := proc()  local OK, AA, NAME, INP, K, N, M, LN1, LM, NL, KK, XX, YY, LLL, J, II, NV, LL, N1, N2, NTR, I, LINE, FLAG, OUP, K1, L, GAMMA, ALPHA, J1, J2, J3, DELTA, A, B, C, I1, JJ1, I2, JJ2, ZZ, HH, JI, XJ, XJ1, XJ2, XI1, XI2, INN, CCC; 
> 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 input
> G2 := 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 then
> Z1 := X;
> fi;
> if 0.5-T <= X and X <= (0.6+T) and abs(Y-0.1) <= T then
> Z1 := X;
> fi;
> if -T <= Y and Y <= 0.1+T and abs(X-0.6) <= T then
> Z1 := Y;
> fi;
> if -T <= X and X <= 0.2+T and abs(Y+X-0.4) <= T then
> Z1 := (X+Y)/sqrt(2);
> fi;
> if 0.4 -T <= X and X <= 0.5+T and abs(Y+X-0.6) <= T then
> Z1 := (X+Y)/sqrt(2);
> fi;
> g2 := Z1;
> RETURN(g2);
> end;
> RR := proc(X,Y,A,B,C,I,I1,I2) local rr;
> rr := R(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(rr);
> end;
> FFF := proc(X,Y,A,B,C,I,I1) local fff;
> fff := F(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(fff);
> end;
> GG1 := proc(X,Y,A,B,C,I,I1,I2) local gg1;
> gg1 := G1(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg1);
> end;
> GG2 := proc(X,Y,A,B,C,I,I1) local gg2;
> gg2 := G2(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(gg2);
> end;
> GG3 := proc(X,Y,A,B,C,I,I2) local gg3;
> gg3 := G2(X,Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg3);
> end;
> GG4 := proc(X,Y,A,B,C,I,I1) local gg4;
> gg4 := G1(X,Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y)*(A[I-1,I1-1]+B[I-1,I1-1]*X+C[I-1,I1-1]*Y);
> RETURN(gg4);
> end;
> GG5 := proc(X,Y,A,B,C,I,I2) local gg5;
> gg5 := G1(X,Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y)*(A[I-1,I2-1]+B[I-1,I2-1]*X+C[I-1,I2-1]*Y);
> RETURN(gg5);
> end;
> QQ := proc(L,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)  local X, Y, I, 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 then
> for I from 1 to 7 do
> S[I-1] := P(X[I-1],Y[I-1]);
> od;    
> fi;
> if L = 2 then
> for I from 1 to 7 do
> S[I-1] := Q(X[I-1],Y[I-1]);
> od;    
> fi;
> if L = 3 then
> for I from 1 to 7 do
> S[I-1] := RR(X[I-1],Y[I-1],A,B,C,I,I1,I2);
> od;    
> fi;
> if L = 4 then
> for I from 1 to 7 do
> S[I-1] := FFF(X[I-1],Y[I-1],A,B,C,I,I1);
> od;    
> fi;
> T1 := 0;
> for I from 1 to 3 do
> T1 := T1+S[I-1];
> od;
> T2 := 0;
> for I from 4 to 6 do
> T2 := T2+S[I-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,I,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 do
> X := (K-1)*H;
> if L = 1 then
> S[K-1] := T3*GG1(T1*X+X1,T2*X+Y1,A,B,C,I,I1,I2);
> fi;
> if L = 2 then
> S[K-1] := T3*GG2(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
> fi;
> if L = 3 then
> S[K-1] := T3*GG3(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
> fi;
> if L = 4 then
> S[K-1] := T3*GG4(T1*X+X1,T2*X+Y1,A,B,C,I,I1);
> fi;
> if L = 5 then
> S[K-1] := T3*GG5(T1*X+X1,T2*X+Y1,A,B,C,I,I2);
> fi;
> od;
> Q3 := S[0]+S[100];
> Q1 := 0;
> Q2 := S[99];
> for K from 1 to 49 do
> Q1 := Q1+S[2*K];
> Q2 := Q2+S[2*K-1];
> od;
> SSQ := (Q3+2*Q1+4*Q2)*H/3;
> sq := SSQ;
> RETURN(sq);
> end;
> printf(`This is the Finite Element Method.\n`);
> OK := FALSE;
> printf(`The input will be from a text file in the following form:\n`);
> printf(`K  N  M  n  m  nl\n\n`);
> printf(`Each of the above is an integer -`);
> printf(`separate with at least one blank.\n\n`);
> printf(`Follow with the input for each node in the form:\n`);
> printf(`x-coord., y-coord., number of triangles in which the`);
> printf(` node is a vertex.\n\n`);
> printf(`Continue with the triangle number and vertex number for\n`);
> printf(`each triangle in which the node is a vertex.\n`);
> printf(`Separate each entry with at least one blank.\n\n`);
> printf(`After all nodes have been entered follow with information\n`);
> printf(`on the lines over which line integrals must be computed.\n`);
> printf(`The format of this data will be the node number of the\n`);
> printf(`starting node, followed by the node number of the ending\n`);
> printf(`node for each line, taken in the positive direction.\n`);
> printf(`There should be 2 * nl such entries, each an integer\n`);
> printf(`separated by a blank.\n\n`);
> printf(`Functions can be input or coded as procedures.\n`);
> printf(`The example has G2 as a procedure and the rest\n`);
> printf(`as input.\n`);