FEMClass`

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

Example 4 :
Solve
D[ f[x,y], {x,2} ] + D[ f[x,y], {y,2} ] - 1 = 0,
f = 0 on x = 0, y = 0, y = 1
D[f,NormalVector] = 1 over x = 1

Needs["FEMClass`"]

Pre-Process

myMalla = 
      NewMalla[T3, {{1, 2, 5}, {5, 2, 3}, {5, 3, 4}, {1, 5, 4}}, {{0, 0}, {1, 
            0}, {1, 1}, {0, 1}, {(1/2), (1/2)}}]

Refinamela[myMalla]; Refinamela[myMalla];

GrafMalla[myMalla]

this is the boundary of the mesh

fm = Frontera[myMalla]

myRes = 
      NewRestricciones[{g[D0, 1] == 0, g[D0, 2] == 0, g[D0, 3] == 0, 
          g[D0, 4] == 0, g[D0, 5] == 0, g[D0, 10] == 0, g[D0, 14] == 0, 
          g[D0, 19] == 0, g[D0, 23] == 0, g[D0, 28] == 0, g[D0, 32] == 0}, {
          g}]

myCampo = NewCampoEF[T3C0, myMalla]

Principal Process

ao = CreaArgOmega[myCampo, myRes];

ag1 = CreaArgGama[myCampo, myRes, ao, fm];

ag2 = 
      CreaArgGama[myCampo, myRes, ao, 
        fm[[{5, 6, 7, 8}]]];

ieo = 
      IntegralesElementalesOmega[myCampo, myRes, 
        ao, {(∂_x) #2 (∂_x) #1 + 
              (∂_y) #2 (∂_y) #1&}, {x, y}];

this perform a integral over all boundary,

ieg1 = 
      IntegralesElementalesGama[myCampo, myRes, ag1, 
        fm, {#1 {(∂_x) #2, (∂_y) #2}.NormalVector&}, {x, 
          y}];

this perform integrals over the elements 5,6,7 and 8 of boundary, the functional include a normal vector " n", to define the normal derivative

\(\(ieg2 = 
      IntegralesElementalesGama[myCampo, myRes, ag2, 
        fm[[{5, 6, 7, 8}]], 

        				{#1\ 
              \(({∂\_x #2, ∂\_y #2}.NormalVector - 1)
                \)&}, {x, y}]; \)\)

\(ito = IntegralesTotales[ieo]; 
\(itg1 = IntegralesTotales[ieg1]; \)\)itg2 = IntegralesTotales[ieg2];

it = itg1 - ito + itg2;

si = SistemaIncompleto[it];

sc = SistemaCompleto[si, myRes];

this expresion form a expresion for the initial values for the variables to find

\(varIni = 
      Transpose[

        		{Incognitas[myCampo, myRes], 
			
          Table[Random[], {Length[Incognitas[myCampo, myRes]]}]
		}]\)

s = FindRoot[Evaluate[sc],  Evaluate[Sequence@@varIni]];

ColumnForm[s]

Post-Process

GrafCampo[myCampo, g, s]

Show[%, ViewPoint -> {-2.738, 0.268, 1.970}]