Examples for the Package BVP2

\(\( (*\ it\ might\ be\ necessary\ to\ set\ the\ current\ directory\ *) 
      
 (*\ SetDirectory["~stahel/tex/wahlfach/Variationsrechnung"]\ *) 
      \)\)

Loading the package and displaying the minimal documentation

<<  ./BVP2.m

?ReadMesh

?MeshGraphics

?FEMSolve

?FEMGradient

?FEMSurface

?FEMSurfaceNoMesh

?LevelGraphics2D

?LevelGraphics3D

Example 1

\({nodes, elements, segments} = ReadMesh["EasyMesh/test1"]; 

    Show[MeshGraphics[nodes, elements]]; \)

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, \(-1\)]; 
gDirichlet = Function[{x, y}, y/10]; 

    gNeumann = Function[{x, y}, \(-1\)]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = 
      Table[Graphics3D[
          Polygon[{points[\([elements[\([k, 2]\)]]\)], 
              points[\([elements[\([k, 3]\)]]\)], 
              points[\([elements[\([k, 4]\)]]\)]}]], {k, Length[elements]}]; 
    
gg1 = \(solution = 
        Show[g, 
	PlotRange\  -> All, 
	Axes -> True, 
	
          AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 

          \ ViewPoint -> {4, \(-4\), 2}]\); \)

\({nodes, elements, segments} = ReadMesh["EasyMesh/test3"]; 

    gMesh = Show[MeshGraphics[nodes, elements]]; \)

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, \(-1\)]; 
gDirichlet = Function[{x, y}, y/10]; 

    gNeumann = Function[{x, y}, \(-1\)]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = \ FEMSurface[points, elements]; 

    gg3 = \(solution = 
        Show[g, 
	PlotRange\  -> All, 
	Axes -> True, 
	
          AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 

          \ ViewPoint -> {4, \(-4\), 2}]\); \)

\(g = \ FEMSurfaceNoMesh[points, elements]; 

    glevel = LevelGraphics3D[points, elements, Table[le, {le, 0, 2, 0.1}]]; 

    gg4 = Show[{g, glevel}, 
	PlotRange\  -> All, 
	Axes -> True, 
	
        AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 
\ 
        ViewPoint -> {4, \(-4\), 4}]; \)

Example 2

{nodes, elements, segments} = ReadMesh["EasyMesh/test2"];

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, 0]; 
gDirichlet = Function[{x, y}, 4 - x]; 

    gNeumann = Function[{x, y}, 0]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = 
      Table[Graphics3D[
          Polygon[{points[\([elements[\([k, 2]\)]]\)], 
              points[\([elements[\([k, 3]\)]]\)], 
              points[\([elements[\([k, 4]\)]]\)]}]], {k, Length[elements]}]; 
    
solution = 
      Show[g, 
	Axes -> True, AxesLabel -> {"x", "y", "u"}, 

        AspectRatio\  -> 1, 
ViewPoint -> {4, \(-4\), 1}]; \)

Example 3, computing the gradient and a flux field

\({nodes, elements, segments} = ReadMesh["EasyMesh/test3"]; 

    Show[MeshGraphics[nodes, elements]]; \)

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, 0]; 
gDirichlet = Function[{x, y}, 1]; 

    gNeumann = Function[{x, y}, 1]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = 
      Table[Graphics3D[
          Polygon[{points[\([elements[\([k, 2]\)]]\)], 
              points[\([elements[\([k, 3]\)]]\)], 
              points[\([elements[\([k, 4]\)]]\)]}]], {k, Length[elements]}]; 
    
solution = 
      Show[g, 
	Axes -> True, AxesLabel -> {"x", "y", "u"}, 

        AspectRatio\  -> 1, 
ViewPoint -> {4, \(-4\), 1}]; \)

gradient = FEMGradient[points, elements];

Needs["Graphics`PlotField`"];
field = Table[{points[[k,{1,2}]],-0.5*gradient[[k]]}
   ,{k,Length[points]}];
gFlux=ListPlotVectorField[field];

\(gMesh = MeshGraphics[nodes, elements]; 

    Show[gMesh, gFlux, AspectRatio -> Automatic]; \)

\(\(
glevel = 
      LevelGraphics2D[points, elements, Table[le, {le, 1, 2.5, 0.1}]]; 

    Show[glevel, gFlux, AspectRatio -> Automatic]; \)\)

Example 4, including comparison with exact solution

{nodes, elements, segments} = ReadMesh["EasyMesh/test3"];

Choose your solution of the equation:

Clear[a, b, f, gNeumann, gDirichlet]fsol[x_, y_]  = x*(y - 3) + Sin[(y/5)^4*6];

Now the terms in the PDE have to be computed

Clear[dxfsol]

\(lapfsol = D[fsol[x, y], {x, 2}] + D[fsol[x, y], {y, 2}]; 

    dxfsolFunc = \ D[fsol[x, y], x]; 

    dxfsol[x_, y_] = \ Evaluate[dxfsolFunc]; 
dyfsol = \ D[fsol[x, y], y]; 

    a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, Evaluate[lapfsol]]; 

    gDirichlet = Function[{x, y}, Evaluate[\ fsol[x, y]]]; 

    gNeumann = Function[{x, y}, Evaluate[dxfsol[x, y]*Sign[x - 2]]]; \)

points = 
      FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, segments];

\(\(Show[FEMSurface[points, elements], 
Axes -> True, 
      AspectRatio\  -> 1]; \)\)

\(dpoints = 
      Table[{points[\([k, 1]\)], points[\([k, 2]\)], 
          points[\([k, 3]\)] - fsol[points[\([k, 1]\)], points[\([k, 2]\)]]}, 
        
	{k, Length[points]}]; 

    maxError = Max[Abs[\(Transpose[dpoints]\)[\([3]\)]]]\)maxValue = Max[Abs[Transpose[points][[3]]]]maxError/maxValue

\(g1 = FEMSurface[dpoints, elements]; 

    g2 = FEMSurface[points, elements]; 

    Show[{g1, g2}, 
	Axes -> True, 
	AspectRatio\  -> \ 1, 
	
      ViewPoint -> {0.921, \ \(-3.248\), \ 0.01}]; \)

\(g3 = FEMSurfaceNoMesh[points, elements]; 

    glevel = LevelGraphics3D[points, elements, 
        Table[le, {le, \(-15\), 5, 1}]]; 

    Show[g3, glevel, 
Axes -> True, AspectRatio\  -> 1]; \)

Compare the x component of the gradient with the correct exact solution

pointsGrad = FEMGradient[points, elements];

\(dxpoints = 
      Table[{points[\([k, 1]\)], points[\([k, 2]\)], 
          pointsGrad[\([k, 1]\)] - 
            dxfsol[points[\([k, 1]\)], points[\([k, 2]\)]]}, 

        	{k, Length[points]}]; 

    maxError = Max[Abs[\(Transpose[dxpoints]\)[\([3]\)]]]\)maxValue = Max[Abs[Transpose[pointsGrad][[1]]]]maxError/maxValue

Example 5

\({nodes, elements, segments} = ReadMesh["EasyMesh/test5"]; 

    Show[MeshGraphics[nodes, elements]]; \)

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, \(-1\)]; 
gDirichlet = Function[{x, y}, 0]; 

    gNeumann = Function[{x, y}, 0]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = FEMSurface[points, elements]; 

    solution = 
      Show[g, 
	Axes -> True, AxesLabel -> {"x", "y", "u"}, 

        AspectRatio\  -> 1, 
ViewPoint -> {\(-4\), 4, 1}]; \)

\(g = FEMSurfaceNoMesh[points, elements]; 

    levels = Table[level, {level, 0, 0.15, 0.01}]; 

    glevel = LevelGraphics3D[points, elements, levels]; 

    solution = 
      Show[{g, glevel}, 
	Axes -> True, 
        AxesLabel -> {"x", "y", "u"}, 
AspectRatio\  -> 1, 

        ViewPoint -> {\(-4\), 4, 1}]; \)

Trivial Mesh Example

\(\(nodes = {{1, 0, 0, 1}, {2, 0.25, 0, 1}, {3, 0.5, 0, 1}, {4, 0.75, 0, 
          1}, {5, 1, 0, 1}, 

        		{6, 0, 0.25, 2}, {7, 0.25, 0.25, 0}, {8, 0.5, 0.25, 0}, {9, 0.75, 
          0.25, 0}, {10, 1, 0.25, 0}, 

        		{11, 0, 0.5, 2}, {12, 0.25, 0.5, 2}, {13, 0.5, 0.5, 0}, {14, 
          0.75, 0.5, 0}, {15, 1, 0.5, 2}, 

        		{16, 0, 0.75, 2}, {17, 0.25, 0.75, 0}, {18, 0.5, 0.75, 0}, {19, 
          0.75, 0.75, 0}, {20, 1, 0.75, 2}, 

        		{21, 0, 1, 1}, {22, 0.25, 1, 1}, {23, 0.5, 1, 1}, {24, 0.75, 1, 
          1}, {25, 1, 1, 1}
	}; \)\)

\(\(elements = 

      	{{1, 1, 7, 6}, {2, 1, 7, 2}, {3, 2, 8, 7}, {4, 2, 3, 8}, 

        	{5, 3, 9, 8}, {6, 3, 4, 9}, {7, 4, 10, 9}, {8, 4, 5, 10}, 

        	{9, 6, 12, 11}, {10, 6, 12, 7}, {11, 7, 13, 12}, {12, 7, 8, 13}, 

        	{13, 8, 14, 13}, {14, 8, 9, 14}, {15, 9, 15, 14}, {16, 9, 10, 15}, 
        
	{17, 11, 17, 16}, {18, 11, 17, 12}, {19, 12, 18, 17}, {20, 12, 
          13, 18}, 

        	{21, 13, 19, 18}, {22, 13, 14, 19}, {23, 14, 20, 19}, {24, 14, 15, 
          20}, 
	{25, 16, 22, 21}, {26, 16, 22, 17}, {27, 17, 23, 22}, {28, 
          17, 18, 23}, 

        	{29, 18, 24, 23}, {30, 18, 19, 24}, {31, 19, 25, 24}, {32, 19, 20, 
          25}}; \)\)

\(\(segments = {
{1, 2, 1}, {2, 3, 1}, {3, 4, 1}, {4, 5, 1}, 
        
{21, 22, 1}, {22, 23, 1}, {23, 24, 1}, {24, 25, 1}, 
        
{1, 6, 2}, {6, 11, 2}, {11, 16, 2}, {16, 21, 2}, 
        
{5, 10, 2}, {10, 15, 2}, {15, 20, 2}, {20, 25, 2}		
	}; \)\)

Show[MeshGraphics[nodes, elements]];

\(a = Function[{x, y}, 1]; 
b = Function[{x, y}, 0]; 

    f = Function[{x, y}, \(-1\)]; 
gDirichlet = Function[{x, y}, 0*y/10]; 

    gNeumann = Function[{x, y}, \(-1\)]; 


    points = FEMSolve[a, b, f, gDirichlet, gNeumann, nodes, elements, 
        segments]; \)

\(g = \ FEMSurface[points, elements]; 

    gg3 = \(solution = 
        Show[g, 
	PlotRange\  -> All, 
	Axes -> True, 
	
          AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 

          \ ViewPoint -> {4, \(-4\), 2}]\); \)

\(g = \ FEMSurfaceNoMesh[points, elements]; 

    glevel = LevelGraphics3D[points, elements, Table[le, {le, 0, 2, 0.1}]]; 

    gg4 = Show[{g, glevel}, 
	PlotRange\  -> All, 
	Axes -> True, 
	
        AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 
\ 
        ViewPoint -> {4, \(-4\), 2}]; \)