10
$\begingroup$

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]

enter image description here

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]

enter image description here

$\endgroup$
11
$\begingroup$

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

triangles

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]

quadrilaterals


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

histogram


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

histogram2

And this is how both meshes look like.

mesh_comparison

$\endgroup$
  • 2
    $\begingroup$ Nice. You could add two histograms Histogram[mesh["Quality"]] and Histogram[quadMesh["Quality"]]. Also interesting is the number of elements. Do you happen to know what algorithm is used for the triangle to quad conversion? Thanks. $\endgroup$ – user21 Apr 4 '17 at 16:36
  • 5
    $\begingroup$ Copied from ´AceFEM´ Help for ´SMTTriangularToQuad´ function: "ToElementMesh directly creates only triangular meshes. Unstructured mesh composed of linear triangles can be converted to unstructured mesh composed of quadrilateral with the help of an additional function defined here. Algorithm is base on procedures described in HOUMAN BOROUCHAKI AND PASCAL J. FRE, ADAPTIVE TRIANGULAR QUADRILATERAL MESH GENERATION, INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 41, 915-934 (1998)." Module for ´SMTTriangularToQuad´ is shown in Help, you can instal AceGen and AceFEM packages. $\endgroup$ – NZupan Apr 5 '17 at 6:53
  • $\begingroup$ @NZupan, thanks for the reference. $\endgroup$ – user21 Apr 5 '17 at 8:24
  • $\begingroup$ Could you do me a favor and add the picture of the smoothed quad mesh? I am curious to see it. Thanks a lot. $\endgroup$ – user21 Oct 2 '17 at 11:50
  • $\begingroup$ Thanks! Now we 'just' need someone to implement the tri to quad algorithm for Mathematica :-P $\endgroup$ – user21 Oct 2 '17 at 12:24

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.