AceFEM: Conversion of triangular to quadrilateral mesh

Question

I am using AceFEM command SMTTriangularToQuad to convert 2D triangular mesh (created by Mathematica) to unstructured quadrilateral mesh. This mesh has internal "domains" (e.g. regions inside the mesh that have different material properties) and conversion failed. Is there a way to solve the problem illustrated in the example bellow?

First I define some filled (primitive) geometric regions and plot their edges.

r1 = Rectangle[{0, 0}, {4, 4}];
r2 = Rectangle[{0.25, 0.25}, {3.75, 3.75}];
d1 = Disk[{2, 2}, 1.25];
d2 = Disk[{2, 2}, 0.5];

Graphics[{
  {FaceForm[], EdgeForm[Black], r1},
  {FaceForm[], EdgeForm[Blue], r2},
  {FaceForm[], EdgeForm[Red], d1},
  {FaceForm[], EdgeForm[Orange], d2}
  }, ImageSize -> 150]

Then I create a derived geometric region with RegionUnion and mesh it with ToElementMesh. Innermost disk is not meshed, it is defined to be a hole. Rectangular and circular internal border is still visible as a path over element edges.

<< NDSolve`FEM`;
(* For conversion to quadrilaterals it is essential that "MeshOrder"->1. *)
mesh = ToElementMesh[
   RegionUnion[
    RegionDifference[r1, r2],
    RegionDifference[d1, d2]
    ],
   "RegionHoles" -> {{2, 2}}, "MeshOrder" -> 1];
Show[mesh["Wireframe"], ImageSize -> 250]

Conversion to quadrilaterals does not respect internal borders between regions. This is visible as a lack of circle shaped path over element edges inside the domain.

<< AceFEM`;
(* SMTTriangularToQuad doesn't load with the package. You have to \
search for it in documentation and evaluate the cell with its \
definition manually. *)
quadMesh = SMTTriangularToQuad[mesh];
Show[quadMesh["Wireframe"], ImageSize -> 250]

Answer 1

"RegionMarker" option has to be specified when you are converting boundary mesh to element mesh. The first element in a list is the coordinates of the point contained in the region to which "RegionMarker" refers to. The second argument is the marker (some kind of region ID) and the third is the maximal cell size in that region. This way you can control mesh density in particular regions of your domain. Using the OP example:

pts = {{0.1, 0.1}, {0.5, 0.5}, {1.5, 1.5}};
mesh = ToElementMesh[
   RegionUnion[
    RegionDifference[r1, r2],
    RegionDifference[d1, d2]
    ],
   "RegionHoles" -> {{2, 2}},
   "MeshOrder" -> 1,
   "RegionMarker" -> {{pts[[1]],1,1},{pts[[2]],2,1},{pts[[3]],3,0.01}}, 
   MeshQualityGoal -> "Maximal"
   ];

(* Showing the location of points specifying distinct sub-regions. *)  
Show[{
  mesh["Wireframe"],
  Graphics[{Red, PointSize[0.025], Point[pts]}]
  }, ImageSize -> 250]

Conversion to quadrilaterals is now successful and respects borders of internal sub-regions. It would be interesting to see how much the quality of the mesh can be further improved by slight repositioning of nodes?

quadMesh = SMTTriangularToQuad[mesh];
Show[
 quadMesh[
  "Wireframe"[
   "MeshElementStyle" -> {FaceForm[Yellow], FaceForm[Orange], 
     FaceForm[Red]}]],
 ImageSize -> 250]

---

EDIT

In the comments, User21 suggested adding some more information about meshes in question. Conversion to quadrilaterals approximately doubles the number of elements in this case.

(Length@*Flatten) /@ {mesh["Quality"], quadMesh["Quality"]}
(* {1276, 2708} *)

Distribution of quality of element shape is worse for quadrilaterals, but I don't know how significant this effect is.

Histogram[
 {Flatten@mesh["Quality"], Flatten@quadMesh["Quality"]},
 {0.02},
 "Probability",
 ChartLegends -> {"Triangles", "Quadrilaterals"}
 ]

---

EDIT 2

Quality of the quadrilateral mesh can be improved by using Laplacian smoothing (function LaplacianElementMeshSmoothing) from this great answer.

smooth = LaplacianElementMeshSmoothing[quadMesh];

Histogram[
 {Flatten@quadMesh["Quality"], Flatten@smooth["Quality"]},
 {0.02},
 "Probability",
 ChartLegends -> {"Original quad mesh", "Smoothed quad mesh"}
 ]

And this is how both meshes look like.