FEMClass`

By Edgar Camacho Cota
Universidad Nacional Autonoma de Mexico (UNAM)
ecc@mira.labvis.unam.mx
June 2000.

Example 1 :
Solve
y''[x] + y'[x] + y[x] = 0, y'[0] = 1, y[0] = 0

Pre-Process

Load the package

Needs["FEMClass`"]

This create an instance of the class L2, this is a mesh of linear elements, with two nodes, the second argument of new, is a list of the nodes in each element, for example, the third element have the nodes 3 and 4; the third argument is a list of cordenates.

myMalla = 
      NewMalla[L2, {{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}}, {{0}, {(2/10)}, {
            (4/10)}, {(6/10)}, {(8/10)}, {1}}]

This graph the mesh

GrafMalla[myMalla]

This stands for the restricctions, the function value of f (cero derivative, D0) in the node 1, is cero; and the firs derivative (Dx)  in node 1 is 1.

myRes = NewRestricciones[{f[D0, 1] == 0, f[Dx, 1] == 1}, {f}]

This create an instance of the class L2C2, that is one dimensional element, with linear geometry, with C2 continuity.

myCampo = NewCampoEF[L2C2, myMalla]

Principal Process

ao is a matrix of arguments for the calculation of the elements Integrals over the domain

ao = CreaArgOmega[myCampo, myRes];

ieo is the matrix of the elements integrals; #2 reprecent the function f, and the #1 reprecent the basis function, need for implement the Galerkin Method.

\(\(ieo = 
      IntegralesElementalesOmega[myCampo, myRes, 
        ao, {#1\ \((∂\_{t, 2}#2 + ∂\_t #2 + #2)\)&}, {
          t}]; \)\)

this sum for each basis

ito = IntegralesTotales[ieo];

it = ito;

si = SistemaIncompleto[it];

sc is the sistem of algebraic equations that solve the problem for each f[D0, n], where D0 is the stands for de cero derivatives, and n is the number of the node, then f[D0,n] is the value of function in the "n" node. Dx, is the stand for one derivative and Dxx the stand for the second derivative.

sc = SistemaCompleto[si, myRes];

s = 
      Solve[sc, Incognitas[myCampo, myRes]][[1
        ]];

ColumnForm[s // N]

this command are not defined, for this element, is more easy with the built in mathematica functions.

GrafCampo[myCampo, g, s]

Post Process

\(\(ps = 
      With[{mm = Coords[myMalla]}, 
			
          Table[{mm[[n, 1]], {f[D0, n], 
                f[Dx, n], f[Dxx, n]}}, {n, Length[mm]}]] /. s; \)\)

nf is the quest function

nf = Interpolation[ps, InterpolationOrder -> 5]

this is the graph of the solution

femSolution = Plot[nf[x], {x, 0, 1}]

this is the residual

Plot[nf''[x]+nf'[x]+nf[x],{x, 0, 1},PlotRange -> All]

this is the solution with built in mathematica fuctions

DSolve[{y''[x]+y'[x]+y[x]==0,y[0] == 0,y'[0]==1},y[x],x];dSolveSolution = Plot[Evaluate[y[x] /. %], {x, 0, 1}]

this is the comparations of FEMClass solution and DSolve Solution

Show[femSolution, dSolveSolution]