Calculating a potential function using the finite element method

Question

This is my first attempt to use the Finite Element method available in version 10. There are questions and I am very open to suggestions. My example is flow around a cylinder which is a well known example of a potential problem. I wish to extend this to more general shapes. I first set up the region which is straightforward.

Needs["NDSolve`FEM`"];
x2 = 2; y2 = 1; r = 1/8;
reg = ImplicitRegion[0 <= x <= x2 &&0 <= y <= y2 && 
(x - x2/2)^2 + (y - y2/2)^2 >= r^2,{x, y}];

I make the mesh.

Question 1: can I just put the refined boundary on the inner circle and not on the outer boundary?

mesh = ToElementMesh[reg, MaxCellMeasure -> 0.0005, "MaxBoundaryCellMeasure" ->0.01];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

There is a problem around the circle at about 1 o'clock—expanded view below.

Question 2: is this a problem with the grid generation or the plotting of the grid?

Show[mesh["Wireframe"], Frame -> True,PlotRange -> {{0.6, 1.4}, {0.3, 0.7}}]

The differential equation is just the Laplacian and a solution is found easily.

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
  NeumannValue[1., x == 0] + NeumannValue[-1, x == 2]}, 
 u, Element[{x, y}, reg]];

The solution looks good except for the artifact close to {0, 0}. This artifact disappears if PlotPoints is increased to {50,25}.

Question 3: Is this just a plotting bug that can be ignored?

Column[{
  Plot3D[sol[x, y], Element[{x, y}, reg], BoxRatios -> {2, 1, 1},
   ImageSize -> 300],
  Plot3D[sol[x, y], Element[{x, y}, reg], BoxRatios -> {2, 1, 1},
   ImageSize -> 300, PlotPoints -> {50, 25}]
}]

A plot of the potential function around the circle, below, looks good.

Question 4: if the grid is made up as a polygon how does the solution find the values that are on the circle but not on the grid?

Plot[Evaluate[sol[1 + (r) Cos[θ], 0.5 + (r) Sin[θ]]],
 {θ, 0,2 π}]

I am interested in the velocity which is the gradient of the potential function. I calculate the velocity using Grad.

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]]

The following should plot the speed of the fluid around the circle.

Question 5: what is going wrong?

ϵ = 0.0;
Plot[Evaluate[Norm[f[1 + (r + ϵ) Cos[θ], 
0.5 + (r + ϵ) Sin[θ]]]], {θ, 0, 2 π}]

Moving off the circle does not help much.

ϵ = 0.1;
Plot[Evaluate[  Norm[f[1 + (r + ϵ) Cos[θ], 
0.5 + (r + ϵ) Sin[θ]]]], {θ, 0, 2 π}]

As usual StreamPlot has a mind of its own and does not produce the streamlines I ask for.

Question 6: is there a way of using the region defined at the start rather than a new RegionFunction?

spts = Table[{0.001, y}, {y, 0, 1, 0.02}];
StreamPlot[f[x, y], {x, 0, 2}, {y, 0, 1}, AspectRatio -> Automatic, 
RegionFunction -> 
 Function[{x, y, vx, vy, n}, (x - x2/2)^2 + (y - y2/2)^2 >= r^2], 
StreamPoints -> spts, 
Epilog -> {Point[spts], Circle[{x2/2, y2/2}, r]}]

This is an edit to my problem which may help to establish the issue as being one of a bad mesh. Starting again but making a finer mesh we see that it is sill bad but good in parts. Question: Are there other methods for making a mesh?

mesh = ToElementMesh[reg, MaxCellMeasure -> 0.00001,
 "MaxBoundaryCellMeasure" -> 0.01];
 Show[mesh["Wireframe"], Frame -> True,PlotRange -> {{0.95, 1.2}, {0.5, 0.7}}]

If we now solve using this faulty mesh rather than in the original question by using a default mesh we have:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
 NeumannValue[1., x == 0] + NeumannValue[-1, x == 2]}, u, 
  Element[{x, y}, mesh]];

Calculating the gradient as before

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]]

The plot of the speed is now

 ϵ = 0.0;
Plot[Evaluate[Norm[f[1 + (r + ϵ) Cos[θ], 
0.5 + (r + ϵ) Sin[θ]]]], {θ, 0, 2 π}]

This is good at locations away from the bad mesh. Thus the problem seems to be a mesh problem with improved mesh resolution giving better results.

Answer 1

Here are a couple of suggestions. First let's look at the mesh generation with is documented with ToElementMesh and in the mesh generation tutorial.

Needs["NDSolve`FEM`"];
x2 = 2; y2 = 1; r = 1/8;
reg = ImplicitRegion[
   0 <= x <= x2 && 
    0 <= y <= y2 && (x - x2/2)^2 + (y - y2/2)^2 >= r^2, {x, y}];

For the mesh generation you have (in V10) two choices, one is a RegionPlot based boundary mesh generator and the other is a continuation based method. The RegionPlot is 'quick and dirty' and the continuation based method is exact but can be slow.

mesh = ToElementMesh[reg,
   "BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 33}
   , "MaxBoundaryCellMeasure" -> 0.01
   ];

By increasing the sample points (which is the same as increasing the plot points) a better mesh is generated. And here is the same with the continuation method.

mesh = ToElementMesh[reg
   , "BoundaryMeshGenerator" -> {"Continuation"}
   , "RegionHoles" -> {{1, 1/2}}
   , "MaxBoundaryCellMeasure" -> 0.01
   ];

A necessary condition that the mesh is approximating the region is that the are is about the same as that from the original region.

a = Integrate[1, {x, y} ∈ reg]
(128 - π)/64

We can then compute:

Total[mesh["MeshElementMeasure"], {2}] - a
{1.67415*10^-9}

Currently, the "MaxBoundaryCellMeasure" will address all boundaries. If you can suggest a convenient syntax for selecting a sub-boundary let me know. DiscretizeRegion only takes the "MaxCellMeasure" option and specifying something like "MaxCellMeasure"->{"Area"->..., "Length"->...} will apply the smallest of both values too all elements, not just the boundary elements.

You could play with a 'MeshRefinementFunction' but careful, as this will subdivide a full mesh and not a boundary mesh. Another option would be to generate a boundary mesh with a Circle, extract the coordinates and boundary elements and add the four corner points and four boundary elements.

After a visual inspection of the mesh, we can solve the PDE:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 
     NeumannValue[1., x == 0] + NeumannValue[-1, x == 2]}, u, 
   Element[{x, y}, mesh], 
   "ExtrapolationHandler" -> {(Indeterminate &), 
     "WarningMessage" -> True}];

A few comments here: First, this is not going to work. You need a DirichletCondition or at the very least a NeumannValue with a Robin part. Please see the documentation to NeumannValue possible issues section.

If you are like me, and do not like that InterpolatingFunction extrapolates, you can switch that off by specifying an "ExtrapolationHandler", please see the Solving PDE with FEM tutorial.

Something like this will then give a message and return Indeterminate

sol[1, 1/2]
InterpolatingFunction::dmval: "Input value {1.,0.5} lies outside the range of data in the interpolating function. Extrapolation will be used."
Indeterminate

To plot over the mesh which we have generated anyways one can use:

Plot3D[sol[x, y], Element[{x, y}, mesh], BoxRatios -> {2, 1, 1}, 
 ImageSize -> 300, Boxed -> False]

This advantageous since the mesh is already generated, using reg the mesh would need to generated again.

The following plot will 'float' as long as no essential conditions like a DirichletCondition are specified. Did you notice the result changes depending on refinement? Maybe NDSolve should give a warning here.

Plot[Evaluate[
  sol[1 + (r) Cos[θ], 0.5 + (r) Sin[θ]]], {θ, 0, 
  2 π}]

If it's wiggly then it is most likely not well approximated with the mesh, this looks fine though:

ClearAll[f];
f[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}, "Cartesian"]]
ϵ = 0.0;
Plot[Evaluate[
  Norm[f[1 + (r + ϵ) Cos[θ], 
    0.5 + (r + ϵ) Sin[θ]]]], {θ, 0, 2 π}]

Summary:

• Question 1: no.
• Question 2: grid generation.
• Question 3: plotting bug.
• Question 4: interpolation.
• Question 5: bad mesh and computing derivatives is always a challenge.
• Question 6: not yet, an upcoming release will hopefully do that.

Have a look at the FEM documentation, many of your questions are addressed there in much greater detail than can be done here. I hope that helps to get you started into the FEM world.