6
$\begingroup$

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["NDSolve`FEM`"];
(*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:

M

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

enter image description here

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

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

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

Needs["NDSolve`FEM`"];
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]
Head[reg12]
(*Polygon*)

If we compare this to the disk case we get

Needs["NDSolve`FEM`"];
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2]
Head[reg12]
(* 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}]]]

enter image description here

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"]

enter image description here

Sorry for the trouble.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – james Apr 6 at 8:39
  • $\begingroup$ @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. $\endgroup$ – user21 Apr 6 at 8:42
  • 1
    $\begingroup$ 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. $\endgroup$ – james Apr 6 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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