Examples for the Package BVP2

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

Loading the package

<< BVP2.m

?ReadMesh

?MeshGraphics

?FEMSolve

?FEMSurface

?FEMSurfaceNoMesh

?LevelGraphics2D

?LevelGraphics3D

Example 0

{nodes, elements, segments} = ReadMesh["test1"];

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}]\);\)\)

\(\(Show[
        FEMSurface[points, elements], 
Axes -> True, 
        AxesLabel -> {"x", "y", "u"}, 
\ AspectRatio\  -> 1, 
\ 
        ViewPoint -> {4, \(-4\), 2}];\)\)

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

\(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}]\); \)

Example 1

{nodes, elements, segments} = ReadMesh["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 2

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

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}];\)\)

Example 3, including comparison with exact solution

{nodes, elements, segments} = ReadMesh["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

lapfsol = D[fsol[x, y], {x, 2}] + D[fsol[x, y], {y, 2}]
dxfsol = D[fsol[x, y], x]
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*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];\)\)

Example 4

{nodes, elements, segments} = ReadMesh["test5"];

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}];\)\)