FEM Simulation: Meshing two Arbitrary objects in an "air" mesh

Question

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:

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

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

Answer 1

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

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.