7
$\begingroup$

Can anyone explain the following behavior? I create a mesh on a Disk in two different ways. First I do it in one step with ToElementMesh. Then, in two steps using ToBoundaryMesh followed by ToElementMesh. Why are the total areas so different?

Seems like the indirect method introduces errors of a first order mesh, even though the order in both cases is 2.

Thanks!

In[201]:= mesh1 = ToElementMesh[Disk[]];

In[194]:= Pi - Total[mesh1["MeshElementMeasure"], 2]

Out[194]= 2.00118*10^-6

In[202]:= bmesh = ToBoundaryMesh[Disk[]];

In[203]:= mesh2 = ToElementMesh[bmesh];

In[197]:= Pi - Total[mesh2["MeshElementMeasure"], 2]

Out[197]= 0.00908652

In[200]:= mesh1["BoundaryElements"] == mesh2["BoundaryElements"]

Out[200]= True
Improve this question
$\endgroup$
3
  • $\begingroup$ Not to mention mesh1["PointElements"] == mesh2["PointElements"] $\endgroup$
    – aardvark2012
    Commented Jun 27, 2017 at 11:31
  • $\begingroup$ @aardvark2012, those are just the IDs of the node coordinates and it does not say anything about how similar a mesh is. Look at the mesh coordinates instead. $\endgroup$
    – user21
    Commented Jun 27, 2017 at 14:46
  • $\begingroup$ @user21 Ah, yes. Of course. My mistake. $\endgroup$
    – aardvark2012
    Commented Jun 28, 2017 at 0:36

1 Answer 1

Reset to default
6
$\begingroup$

Most importantly I'd like to recommend to take a look at the Element Mesh Generation tutorial. That tutorial explains mesh generation for numerical applications like the Finite Element Method and covers you question. If anything is unclear there let me know and it can be improved.

I'll try do give a different explanation than given in the tutorial next. Let's consider for a minute that we have a boundary element mesh like the following:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh["Coordinates" -> {
{1., 0.}, {0.9125378206934781, 0.4089923297618155},{0.6654505497212123, 0.7464419373774067}, {0.32914518683708227,  0.9442793262493796}, {2.8415758474179748*^-8, 0.9999999999999996}, {-0.40899232976181543, 0.9125378206934781}, {-0.7464419373774067, 0.6654505497212122}, {-0.9442793262493796, 0.3291451868370823}, {-0.9999999999999996, 2.8415758313367482*^-8}, {-0.9125378206934783, -0.40899232976181493}, {-0.6654505497212126, -0.7464419373774064}, {-0.3291451868370832, -0.9442793262493793}, {-2.841576059504587*^-8, -0.9999999999999996}, {0.40899232976181327, -0.9125378206934791}, {0.7464419373774028, -0.6654505497212166}, {0.9442793262493757, -0.3291451868370933}}, 
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 10}, {10, 11}, {11, 12}, {12, 13}, {13, 14}, {14, 15}, {15, 16}, {16, 1}}]}];

When we pass that to ToElementMesh it will generate a full mesh that approximates the boundary given. ToElementMesh has no way of knowing what the original input to ToBoundaryMesh was when it is given a boundary mesh representation.

mesh = ToElementMesh[bmesh];

So how well does ToElementMesh approximate bmesh? We can not really tell because we do not know what bmesh is.

Now, I am telling you that bmesh is supposed to represent a Disk[]. Then and only then we can check:

Pi - Total[mesh["MeshElementMeasure"], 2]
0.0819661

And it is a poor presentation. If you do have a symbolic representation of your region it's a good idea to pass that along. That is what ToNumericalRegion is for. Let's look at an example. This generates a numerical region of a Disk[]:

nr = ToNumericalRegion[Disk[]];

We can now 'fill in' a boundary mesh like so:

bem2 = ToBoundaryMesh[nr, "MaxBoundaryCellMeasure" -> .5, 
  AccuracyGoal -> 1]

These options are in fact the ones I used to generate the above example boundary mesh. Note that now the NumericalRegion has a boundary mesh - the same as bem2

bem2 === nr["BoundaryMesh"]
True

When you pass the numerical region to ToElementMesh things are very different as now ToElementMesh has access to the boundary representation and the symbolic representation of the region and can thus generate a better mesh.

mesh2 = ToElementMesh[nr];
Pi - Total[mesh2["MeshElementMeasure"], 2]
3.893310501545955`*^-6

When you call

ToElementMesh[Disk[]]

then ToElementMesh does have access to the symbolic region. In fact internally, it generates a NumericalRegion just as in this post and proceeds like shown here.

You can also set a (boundary) mesh to a NumericalRegion:

SetNumericalRegionElementMesh[nr, bmesh]
mesh2 = ToElementMesh[nr]
Pi - Total[mesh2["MeshElementMeasure"], 2]

Hope that helps.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks for your thoughtful reply. The very last comment is exactly what I need. I understand that this function SetNumericalRegionElementMesh can pass the information that allows ToElementMesh to construct a good representation of something like a disk. The information contained in a boundary mesh is insufficient. The tutorial you mention does not seem to cover SetNumericalRegionElementMesh. This seems to provide the functionality needed to create an element mesh with a custom boundary mesh -- which was the problem that motivated my question. So that was indeed helpful! $\endgroup$
    – Will.Mo
    Commented Jun 28, 2017 at 11:47
  • 1
    $\begingroup$ @Will.Mo, SetNumericalRegionElementMesh is actually in the mentioned tutorial. Look for the section called 'Numerical Regions'. Since SetNumericalRegionElementMesh is there only as program code, it might not be searchable in a browser. It is, however, searchable in the tutorial notebook under FEMDocumentation/tutorial/ElementMeshCreation. If you have other suggestions for improvement, let me know. $\endgroup$
    – user21
    Commented Jun 28, 2017 at 13:39
  • $\begingroup$ Sorry, I missed that. Thanks for the tips, and thanks for teaching me how to search the Mathematica tutorials correctly :). $\endgroup$
    – Will.Mo
    Commented Jun 28, 2017 at 14:27
  • 1
    $\begingroup$ @Will.Mo, I must admit that when I searched for it on the web page I did not find it either only when I wanted to add it to the tutorial I did see that I had already done it.... haha. Out of curiosity and if you can say, what are you working on? $\endgroup$
    – user21
    Commented Jun 28, 2017 at 14:29
  • $\begingroup$ my project involves solving an elliptic PDE with an "oblique boundary condition". I've managed to get something that works for the "regular" case (where the direction of the derivative is nowhere tangential to the boundary) and now I'm working on the irregular case, and I believe I will need a bit more control with the boundary mesh. $\endgroup$
    – Will.Mo
    Commented Jun 28, 2017 at 14:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.