# FEM Simulation: Meshing two Arbitrary objects in an “air” mesh

Based on this excellent answer to my previous question: FEM: Electric Field between two arbitrary defined shapes

I wanted to try out other shapes:

Needs["NDSolveFEM"];
(*Define Boundaries*)
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Rectangle[{-2.5, 2.5}, {2.5, 2}];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12]

mesh = ToElementMesh[reg, {{-5, 5}, {-5, 5}},
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]


But all I am getting is this:

When I use the Diskas object 1, everything works correctly:

So, I was wondering: Why does it not work with other shapes ?

• Weird. Most surprisingly, it works when leaving away the bounding box information {{-5, 5}, {-5, 5}}... – Henrik Schumacher Apr 6 at 7:16
• @HenrikSchumacher Good point. Yes, very weird. – james Apr 6 at 7:34

At some point in time the geometry team decided to return a Polygon object in cases like this:

Needs["NDSolveFEM"];
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Rectangle[{-2.5, 2.5}, {2.5, 2}];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2]
(*Polygon*)


If we compare this to the disk case we get

Needs["NDSolveFEM"];
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2]
(* BooleanRegion *)


My opinion is that this is not a good change, since for a Polygon you do not know if it is an exact representation of the original geometry or only an approximation. For example is this:

Graphics[Polygon[
Table[{Cos[2 \[Pi] k/6], Sin[2 \[Pi] k/6]}, {k, 0, 5}]]]


a crude approximation to a disk, or is this the intended shape? You can not tell. For FEM not being able to tell makes a difference; for example for a second order mesh for an inexact region (like Polygon) one does not know where to move the mid side nodes to. This is different for the BooleanRegion object returned in the disk case. This is an exact symbolic representation of the region and thus preferable for FEM.

That being said, it seems that the boundary intersection algorithm (= giving region bounds) does not work for (this?) Polygon. Whether this is a bug or the bounds algorithm needs improvements or if this is as designed I'd need to investigate.

Luckily, as has been pointed out in the comments, the workaround is simple. Omit the bounding box:

mesh = ToElementMesh[reg,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]];
mesh["Wireframe"]


Sorry for the trouble.

• Thank you very much for your detailed investigation and explanation of the source of error! Here the mesh refinement function makes no sens, right ? It would need to be finer around the edges of the polygons not in the center. – james Apr 6 at 8:39
• @james, depends on what you want to do. If you have some sort of a source at the center of the refinement function, then perhaps. Use a mesh refinement in a specific area if you see a the change in a solution after refining the mesh with MaxCellMeasure. – user21 Apr 6 at 8:42
• Okay, thanks ! In this particular case, I just have two charged plates (Parallel Plate Capacitor) hence I could omit the mesh refinement in the center. – james Apr 6 at 8:45