(* 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, then
*       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
*       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.
*    5) change definition statements to read :
*       A(M,3),B(M,3),C(M,3),XX(LM),YY(LM)
*       definition LINE(NL,2),LL(LM),II(LM,MAX LL(LM)),
*            NV(LM,MAX LL(LM))
*       definition NTR(M,3),GAMMA(LM),ALPHA(LN,LN+1), use
*       the appropriate numbers for the variables M, LM,
*       NL.
*)
P[x_,y_] := Evaluate[1];
Q[x_,y_] := Evaluate[1];
R[x_,y_] := Evaluate[0];
F[x_,y_] := Evaluate[0];
G[x_,y_] := Evaluate[4];
G1[x_,y_] := Evaluate[0];
G2[X_,Y_] := Module[{T,Z1,g2},
   Z1 = 0;
   If[0.0 <= X && X <= 0.2,
      Z1 = X;
   ];
   If[0.2 < X && X <= 0.4,
      Z1 = (X+Y)/Sqrt[2];
   ];
   g2 = Z1;
   g2
];
RR[X_,Y_,A_,B_,c_,i_,I1_,I2_] := Module[{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);
   rr
];
FFF[X_,Y_,A_,B_,c_,i_,I1_] := Module[{fff},
   fff = F[X,Y]*(A[i-1,I1-1]+B[i-1,I1-1]*X+c[i-1,I1-1]*Y);
   fff
];
GG1[X_,Y_,A_,B_,c_,i_,I1_,I2_] := Module[{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);
   gg1
];
GG2[X_,Y_,A_,B_,c_,i_,I1_] := Module[{gg2},
   gg2 = G2[X,Y]*(A[i-1,I1-1]+B[i-1,I1-1]*X+c[i-1,I1-1]*Y);
   gg2
];
GG3[X_,Y_,A_,B_,c_,i_,I2_] := Module[{gg3},
   gg3 = G2[X,Y]*(A[i-1,I2-1]+B[i-1,I2-1]*X+c[i-1,I2-1]*Y);
   gg3
];
GG4[X_,Y_,A_,B_,c_,i_,I1_] := Module[{gg4},
   gg4 = G1[X,Y]*(A[i-1,I1-1]+B[i-1,I1-1]*X+c[i-1,I1-1]*Y)^2;
   gg4
];
GG5[X_,Y_,A_,B_,c_,i_,I2_] := Module[{gg5},
   gg5 = G1[X,Y]*(A[i-1,I2-1]+B[i-1,I2-1]*X+c[i-1,I2-1]*Y)^2;
   gg5
];
QQ[L_,XX_,YY_,DELTA_,J1_,J2_,J3_,I1_,I2_,A_,B_,c_] := Module[{X,Y,i,T1,T2,T3,QQQ,qq},
   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,
      For[i = 1, i <= 7, i++,
	 S[i-1] = P[X[i-1],Y[i-1]];
      ];
   ];
   If[L == 2,
      For[i = 1, i <= 7, i++,
	 S[i-1] = Q[X[i-1],Y[i-1]];
      ];
   ];
   If[L == 3,
      For[i = 1, i <= 7, i++,
     	 S[i-1] = RR[X[i-1],Y[i-1],A,B,c,i,I1,I2];
      ];
   ];
   If[L == 4,
      For[i = 1, i <= 7, i++,
	 S[i-1] = FFF[X[i-1],Y[i-1],A,B,c,i,I1];
      ];
   ];
   T1 = 0;
   For[i = 1, i <= 3, i++,  
      T1 = T1+S[i-1];
   ];
   T2 = 0;
   For[i = 4, i <= 6, i++,
      T2 = T2+S[i-1];
   ];
   T3 = S[6];
   QQQ = 0.5*(T1/20+2*T2/15+9*T3/20)*Abs[DELTA];
   qq = QQQ;
   qq
];
SQ[L_,XX_,YY_,J1_,J2_,i_,I1_,I2_,A_,B_,c_] := Module[{X1,Y1,X2,Y2,H,T1,T2,T3,K,X,Q3,Q1,Q2,SSQ,sq},
   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 = 1, K <= 101, K++,
      X = (K-1)*H;
      If[L == 1,
	 S[K-1] = T3*GG1[T1*X+X1,T2*X+Y1,A,B,c,i,I1,I2];
      ];
      If[L == 2,
	 S[K-1] = T3*GG2[T1*X+X1,T2*X+Y1,A,B,c,i,I1];
      ];
      If[L == 3,
	 S[K-1] = T3*GG3[T1*X+X1,T2*X+Y1,A,B,c,i,I2];
      ];
      If[L == 4,
	 S[K-1] = T3*GG4[T1*X+X1,T2*X+Y1,A,B,c,i,I1];
      ];
      If[L == 5,
	 S[K-1] = T3*GG5[T1*X+X1,T2*X+Y1,A,B,c,i,I2];
      ];
   ];
   Q3 = S[0]+S[100];
   Q1 = 0;
   Q2 = S[99];
   For[K = 1, K <= 49, K++,
      Q1 = Q1+S[2*K];
      Q2 = Q2+S[2*K-1];
   ];
   SSQ = (Q3+2*Q1+4*Q2)*H/3;
   sq = SSQ;
   sq
];
OK = 0;
Print["The input will be from a text file in the following form:\n"];
Print["K N M n m n1\n"];
Print["\n"];
Print["Each of the above is an integer - \n"];
Print["  separate with at least one blank\n"];
Print["\n"];
Print["Follow with the input for each node in the form:\n"];
Print["x-coord., y-coord., number of triangles in which the\n"];
Print["node is a vertex.\n"];
Print["\n"];
Print["Continue with the triangle number and the vertex number for\n"];
Print["each triangle in which the node is a vertex.\n"];
Print["Separate each entry with at least one blank\n"];
Print["\n"];
Print["After all of the nodes have been entered follow with information\n"];
Print["on the lines over which line integrals must be computed\n"];
Print["The format of this data will be the node number of the\n"];
Print["starting node, followed by the node number of the ending\n"];
Print["node for each line, taken in the positive direction.\n"];
Print["There should be 2 * nl such entries,each an integer\n"];
Print["separated by a blank.\n"];
Print["\n"];
Print["make sure the functions P,Q,R,F,G,G1,G2 have been created \n"];
Print["along with the input file \n"];
AA = InputString["This is the Finite Element Method\n
\n
The input will be from a text file in the\n
following form (see screen):\n
Has the input file been created? 'yes' or 'no'\n"];
If[AA == "yes" || AA == "Y" || AA == "y",
   NAME = InputString["Input the file name\nIn the form - drive:\\name.ext\nFor example:   a:\output.dta\n"]; 
   INP = OpenRead[NAME];
   OK = 1;
   K = Read[INP,Number];
   n = Read[INP,Number];
   M = Read[INP,Number];
   LN1 = Read[INP,Number];
   LM = Read[INP,Number];
   NL = Read[INP,Number];
   For[KK = 1, KK <= LM, KK++,
      XX[KK-1] = Read[INP,Number];
      YY[KK-1] = Read[INP,Number];
      LLL = Read[INP,Number];
      For[J = 1, J <= LLL, J++,
	 II[KK-1,J-1] = Read[INP,Number];
	 NV[KK-1,J-1] = Read[INP,Number];
      ];
      LL[KK-1] = LLL;
      For[J = 1, J <= LLL, J++,
	 N1 = II[KK-1,J-1];
	 N2 = NV[KK-1,J-1];
	 NTR[N1-1,N2-1] = KK;
      ];
   ];
   If[NL > 0,
      For[i = 1, i <= NL, i++,
	 For[J = 1, J <= 2, J++,
	    LINE[i-1,J-1] = Read[INP,Number];
         ];
      ];
   ];
   Close[INP],
   Input["The program will end so that the input file\n
   can be created\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"FINITE ELEMENT METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"\n"];
   K1 = K+1;
   N1 = LN1+1;
   Write[OUP,"Verticies and nodes of Triangles\n"];
   Write[OUP,"Triangle-node number for vertex 1 to 3\n"];
   For[i = 1, i <= M, i++,
      Write[OUP,i];
      For[J = 1, J <= 3, J++,
	 Write[OUP,NTR[i-1,J-1]];
      ];
      Write[OUP,"\n"];
   ];
   Write[OUP,"x and y coordinates of nodes\n"];
   For[KK = 1, KK <= LM, KK++,
      Write[OUP,KK,"      ",SetPrecision[XX[KK-1],9],"      ",
            SetPrecision[YY[KK-1],9]];
   ];
   Write[OUP,"Lines of the Domain\n"];
   For[KK = 1, KK <= NL, KK++,
      Write[OUP,KK,"      ",LINE[KK-1,0],"      ",LINE[KK-1,1]];
   ];
   (*STEP 1*)
   If[LM != LN1,
      For[L = N1, L <= LM, L++,
	 GAMMA[L-1] = G[XX[L-1],YY[L-1]];
      ];
   ];
   (* step 2 - initialization of ALPHA - note that ALPHA[i,LN1+1] = BETA *)
   For[i = 1,i <= LN1,i++,
      For[J = 1,J <= N1,J++,
         ALPHA[i-1,J-1] = 0;
      ];
   ];
   (* 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 br vertex 3 *)
   (* step 3 *)
   For[i = 1,i <= M,i++,
      J1 = NTR[i-1,0];
      J2 = NTR[i-1,1];
      J3 = NTR[i-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[i-1,0] = (XX[J2-1]*YY[J3-1]-YY[J2-1]*XX[J3-1])/DELTA;
      B[i-1,0] = (YY[J2-1]-YY[J3-1])/DELTA;
      c[i-1,0] = (XX[J3-1]-XX[J2-1])/DELTA;
      A[i-1,1] = (XX[J3-1]*YY[J1-1]-YY[J3-1]*XX[J1-1])/DELTA;
      B[i-1,1] = (YY[J3-1]-YY[J1-1])/DELTA;
      c[i-1,1] = (XX[J1-1]-XX[J3-1])/DELTA;
      A[i-1,2] = (XX[J1-1]*YY[J2-1]-YY[J1-1]*XX[J2-1])/DELTA;
      B[i-1,2] = (YY[J1-1]-YY[J2-1])/DELTA;
      c[i-1,2] = (XX[J2-1]-XX[J1-1])/DELTA;
      (* step 4 *)
      (* I1 = J for step 4 and I1 = K for step 7 *)
      For[I1 = 1,I1 <= 3,I1++,
         (* step 8 *)
         JJ1 = NTR[i-1,I1-1];
	 (* I2 = K for step 4 and I2 = J for step 9 *)
         For[I2 = 1,I2 <= I1,I2++,
	    (* step 4 and step 10 *)
            JJ2 = NTR[i-1,I2-1];
            ZZ = B[i-1,I1-1]*B[i-1,I2-1]*QQ[1,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,c]+c[i-1,I1-1]*c[i-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 12 *)
            If[JJ1 <= LN1,
               If[JJ2 <= LN1,
                  ALPHA[JJ1-1,JJ2-1] = ALPHA[JJ1-1,JJ2-1]+ZZ;
                  If[JJ1 != JJ2,
                     ALPHA[JJ2-1,JJ1-1] = ALPHA[JJ2-1,JJ1-1]+ZZ;
                  ],
                  ALPHA[JJ1-1,N1-1] = ALPHA[JJ1-1,N1-1]-ZZ*GAMMA[JJ2-1];
               ],
               If[JJ2 <= LN1,
                  ALPHA[JJ2-1,N1-1] = ALPHA[JJ2-1,N1-1]-ZZ*GAMMA[JJ1-1];
               ];
            ];
         ];
         HH = -QQ[4,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,c];
         If[JJ1 <= LN1,
            ALPHA[JJ1-1,N1-1] = ALPHA[JJ1-1,N1-1]+HH;
         ];
      ];
   ];
   (* output the basis fuctions *)
   Write[OUP,"Basis Functions\n"];
   Write[OUP,"Triangle - Vertix - Node - Function\n"];
   For[i = 1,i <= M,i++,
      For[J = 1,J <= 3,J++,
         Write[OUP,i," ",J," ",NTR[i-1,J-1]," ",
           SetPrecision[A[i-1,J-1],9]," ",SetPrecision[B[i-1,J-1],9]," ",
            SetPrecision[c[i-1,J-1],9]];
      ];
   ];
   (* step 5 *)
   (* for each line segment J1 = 1,...,NL and for each triangle I,I = K1,...,N
      which may contain line JI search for all 3 vertices foe possible correspondences *)
   (* step 5 and steps 13-19 *)
   If[NL != 0 && n != K,
      For[JI = 1,JI <= NL,JI++,
         For[i = K1,i <= n,i++,
            For[I1 = 1,I1 <= 3,I1++,
		(* I1 = J in step 5 and I1 = K in step 4 *)
                (* step 15 *)
                J1 = NTR[i-1,I1-1];
                If[LINE[JI-1,0] == J1,
                   For[I2 = 1,I2 <= 3,I2++,
		      (* I2 = K in step 5 and I2 = J in step 16 *)
	 	      (* step 17 *)
                      J2 = NTR[i-1,I2-1];
                      If[LINE[JI-1,1] == J2,
(* there is a correspondence of vertex I1 in triangle I with node J1 as the 
   start of line JI an vertex I2 with node J2 as the end of line JI *)
			 (* step 5 *)
                         XJ = SQ[1,XX,YY,J1,J2,i,I1,I2,A,B,c];
                         XJ1 = SQ[4,XX,YY,J1,J2,i,I1,I2,A,B,c];
                         XJ2 = SQ[5,XX,YY,J1,J2,i,I1,I2,A,B,c];
                         XI1 = SQ[2,XX,YY,J1,J2,i,I1,I2,A,B,c];
                         XI2 = SQ[3,XX,YY,J1,J2,i,I1,I2,A,B,c];
	  		 (* steps 8 and 19 *)
                         If[J1 <= LN1,
                            If[J2 <= LN1,
                               ALPHA[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,
                               ALPHA[J1-1,N1-1] = ALPHA[J1-1,N1-1]-XJ*GAMMA[J2-1]+XI1;
                               ALPHA[J1-1,J1-1] = ALPHA[J1-1,J1-1]+XJ1;
                            ],
                            If[J2 <= LN1,
                               ALPHA[J2-1,N1-1] = ALPHA[J2-1,N1-1]-XJ*GAMMA[J1-1]+XI2;
                               ALPHA[J2-1,J2-1] = ALPHA[J2-1,J2-1]+XJ2;
                            ];
                         ];
                      ];
                   ];
                ];
             ];
          ];
       ];
    ];
   (* output ALPHA *)
   Write[OUP,"Matrix ALPHA follows:\n"];
   For[i = 1,i <= LN1,i++,
      Write[OUP,"Row ",i];
      For[J = 1,J <= N1,J++,
         Write[OUP,SetPrecision[ALPHA[i-1,J-1],9]];
      ];
   ];
   Write[OUP,"\n"];
   (* step 20 *)
   If[LN1 > 1,
      INN = LN1-1;
      For[i = 1,i <= INN,i++,
         I1 = i+1;    
            For[J = I1,J <= LN1,J++,
            CCC := ALPHA[J-1,i-1]/ALPHA[i-1,i-1];
            For[J1 = I1,J1 <= N1,J1++,
               ALPHA[J-1,J1-1] = ALPHA[J-1,J1-1]-CCC*ALPHA[i-1,J1-1];
            ];
            ALPHA[J-1,i-1] = 0;
         ];
      ];
   ];
   GAMMA[LN1-1] = ALPHA[LN1-1,N1-1]/ALPHA[LN1-1,LN1-1];
   If[LN1 > 1,
      For[i = 1,i <= INN,i++,
         J = LN1-i;
         CCC = ALPHA[J-1,N1-1];
         J1 = J+1;
         For[KK = J1,KK <= LN1,KK++,
            CCC = CCC-ALPHA[J-1,KK-1]*GAMMA[KK-1];
         ];
         GAMMA[J-1] = CCC/ALPHA[J-1,J-1];
      ];
   ];
   (* step 21 *)
   (* output GAMMA *)
   Write[OUP,"Coefficients of Basis Functions:\n"];
   For[i = 1,i <= LM,i++,
      LLL = LL[i-1];
      Write[OUP,i," ",SetPrecision[GAMMA[i-1],9]," ",i];
      For[J = 1,J <= LLL,J++,
         Write[OUP,II[i-1,J-1]];
      ];
   ];
   Write[OUP,"\n"];
];
Quit[];