Solving Laplace Equation using Finite Elements

Questions
, ,
1

Hi there,

I was excited to read that the newer versions of Mathematica have the capacity of carrying out Finite Element calculations for solving Partial Differential Equations. I bought Mathematica v.10 and was all set to go … or so I thought …

For starters, I tried a problem that’s so simple I don’t even need a PC to solve it: the Laplace equation in spherical coordinates in a domain consisting of the shell between two concentric spheres with simple (constant) Dirichlet conditions on the surface of the two spheres.

In the attachment below is the notebook with my 4 lines of code. As far as I can figure out, I did everything just as in the examples given in the relevant Chapter of the Wolfram Documentation section, but this is apparently not so, since NDSolve doesn’t do what it’s supposed to do.The error message in the output leaves me nonplussed.

Can some kind spirit tell me what I’m doing wrong?

Thanks already,
René Samson

2

Non-success messages are a kind of stupid, because everybody is inapt to do this or that, if you only search a bit; but nevertheless I cannot get this to work. A few remarks

  • in spherical co-ordinates 1. <= r <= 100. && 0 <= r1 <= Pi && 0 <= r2 <= 2 Pi is not an ImplicitRegion, but an explicit one
  • Csc[r1] has singularities at each multiple of Pi
  • even this NDSolveValue::femcnmd message is not found in the Message folder S:\Program Files\WolframResearch\Mathematica\10.0.1\Documentation\English\System\ReferencePages\Messages
  • it might be rewarding to use cartesian co-ordinates and as ImplicitRegion in that case the Ball
  • where exactly did you find this example in the Wolfram
    documentation?
3

Dear Dr. Krause,

Thanks a great deal for your reaction.

To explain better what I did and why I did it:
I am a relative newcomer to Mathematica and this is the first time I tried to use the 3D-PDE module of Mathematica. I tried my hand on this problem because:

  1. the domain is finite so that it lends itself to a Finite Element treatment;
  2. the solution is a very simple analytical expression. It’s nice to have a check on the Mathematica outcome.

To stay on the safe side I tried - as much as possible - to copy code I read in the Mathematica documentation on this subject (Solving Partial Differential Equations with Finite Elements—Wolfram Documentation). I agree with you that there’s nothing implicit in my first line of code. I also tried the region definition without the ImplicitRegion statement; e.g. as a RegionDifference between a Ball of radius 100 and a Ball of radius 1. That didn’t help, unfortunately: I got the same error message.

I will try out your suggestion to do everything consistently in Cartesian coordinates although it “feels” very weird to do that in a geometry that so obviously has spherical symmetry. I will let you know whether that works.

Thanks, René Samson

4

FEM is used to have as much as possible flexibility about the regions, the region is triangulated and the basic tiles (triangles, tetrahedra, …) do usually not preserve any symmetry of the original problem: they are by refinement procedures fitted to the region; let’s do it for the shell in 3D (see notebook), it ends in an

During evaluation of In[49]:= NDSolveValue::fememib: The input has or generated an intersecting boundary and cannot be processed. >>
During evaluation of In[49]:= NDSolveValue::femtemnm: A mesh could not be generated. >>

that did not work. The intersecting boundary is again some artefact. But it works in 2D

comm353406.PNG
comm353406.PNG1281×833 44.4 KB

5

It has to do with the ToElementMesh[] function

In[1]:= Needs["NDSolve`FEM`"]

In[2]:= Clear[\[CapitalOmega]]
\[CapitalOmega] = ImplicitRegion[2 <= (x^2 + y^2 + z^2) <= 4, {x, y, z}]

Out[3]= ImplicitRegion[2 <= x^2 + y^2 + z^2 <= 4, {x, y, z}]

In[4]:= mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 0.01];
During evaluation of In[4]:= ToElementMesh::fememib: The input has or generated an intersecting boundary and cannot be processed. >>
During evaluation of In[4]:= ToElementMesh::femtemnm: A mesh could not be generated. >>

this function has interesting options,

In[5]:= Options[ToElementMesh]

Out[5]= {"BoundaryMarkerFunction" -> None, 
 "BoundaryMeshGenerator" -> Automatic, 
 "CheckIncidentsCompletness" -> True, "CheckIntersections" -> True, 
 "CheckQuality" -> Automatic, "DeleteDuplicateCoordinates" -> True, 
 "ElementMeshGenerator" -> Automatic, 
 "ImproveBoundaryPosition" -> Automatic, "IncludePoints" -> {}, 
 "MaxBoundaryCellMeasure" -> Automatic, "MeshElementBlocks" -> 1, 
 "MeshElementConstraint" -> Automatic, "MeshElementType" -> Automatic,
  "MeshOrder" -> Automatic, "MessageHead" -> Automatic, 
 "NodeReordering" -> Automatic, "PointMarkerFunction" -> None, 
 "RegionHoles" -> Automatic, "RegionMarker" -> None, 
 "SteinerPoints" -> Automatic, AccuracyGoal -> Automatic, 
 MaxCellMeasure -> Automatic, MeshQualityGoal -> Automatic, 
 MeshRefinementFunction -> None, PrecisionGoal -> Automatic, 
 PerformanceGoal :> $PerformanceGoal}

some of them, like RegionHoles are still undocumented. Nevertheless, the solid ball works in 3D out of the box

comm352406_1.PNG
comm352406_1.PNG1130×953 130 KB

only the result has to be understood.

6

Dear Dr.Krause, Thanks so much for your comment above. The essence of the problem is apparently that the grid-constructing algorithm, while doing fine on certain geometries, has problems coping with certain other geometries. For example: the interior of a ball is no problem at all, but the shell between two concentric shells is problematic. I find it difficult to understand why. If the gap between the two balls is many orders of magnitude larger than the maimum size of the grid element, then this should not be a problematic geometry for a smart grid algorithm. In the attachment below you will find an example where this goes wrong (radius ball #1 = 1; radius ball #2 = 100; size of maximum grid element = 0.01).

What’s the problem here? There simply MUST be a simple trick to overcome this problem. I wonder whether the Wolfram developers have an opinion on this.

Once again, my most sincere thanks for your time and energy,
René Samson

7

In contrast to the case with two concentric Balls which generate the errors when executing ToElementMesh on the RegionDifference, using Cuboids does not generate the error:

Needs["NDSolve`FEM`"]; 

inner = 10;
outer = 20;

c1 = Cuboid[{-inner, -inner, -inner}, {inner, inner, inner}];

c2 = Cuboid[{-outer, -outer, -outer}, {outer, outer, outer}];

\[CapitalOmega] = RegionDifference[c2, c1];

RegionPlot3D[\[CapitalOmega], PlotStyle -> Opacity[.5]] 

mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> {"Volume" -> 1}];

The generation of the mesh is time consuming and especially so as the “Volume” is decreased. One can see that the number of mesh elements in this case is already quite significant:

In[54]:= mesh

Out[54]= ElementMesh[{{-20., 20.}, {-20., 20.}, {-20., 
   20.}}, {TetrahedronElement["<" 106250 ">"]}]

All of this suggests a bug in the Ball case…

I also encountered Udo’s bug ( http://community.wolfram.com/groups/-/m/t/357364 ) when exploring plots of the above RegionPlot3D.

8

It seems that this problem has been solved in the most recent release of Mathematica (10.0.2). The mesh generation which failed before in the shell region between two concentric balls now seems to work fine. I am happy that this problem has been solved and look forward to start working with this tool. (René Samson).

9

Dear Rene,
Would you be able to post your code of a working solution of the Laplace equation in the spherical shell domain? I would find that useful.
Kind regards, WCC

10

Dear Craig, In the attachment please find my “working solution” for the Laplace eqn in a spherical shell domain. I put “working solution” in quotation marks because obviously this isn’t yet a correct solution. I suspect that I have to refine the boundary mesh to get a decent solution and I didn’t have the time to do that yet, but at least the error messages I was getting in the earlier Mathematica release which told me - basically - that the programme failed to set up a mesh in the shell domain, have disappeared, and I do get a solution to the NDSolve-command. Progress …, but still some work left to be done. I hope you find this useful. (René Samson).