Mesh refinement along grain boundaries

Question

I have again a mesh related question. I want to do simulations on a microstructure and therefore I have to mesh it. That works fine, but the refinement of the grain boundaries does not work and let "ToElementMesh" crash and/or takes very long time.

Needs["FEMAddOns`"]
<< NDSolve`FEM`

(*Points*)
pts = {{0, 0}, {3333, 0}, {3333, 2166}, {0, 2166}, {1150, 
    0}, {1100, 900}, {1050, 950}, {0, 1333}, {3333, 1000}, {3066, 
    1500}, {2900, 1700}, {2800, 1750}, {1120, 900}, {3000, 2166}};

(*Create the BMeshes from pts*)
densBound = 10;(*Is intended to refine the boundaries, but is extremely time consuming and/or let ToElementMesh abort*)
a1 = 
 ToBoundaryMesh[
  Line[{pts[[1]], pts[[5]], pts[[6]], pts[[7]], pts[[8]], pts[[1]]}]
  (*  ,"MaxBoundaryCellMeasure"->densBound*)]
a2 = ToBoundaryMesh[
  Line[{pts[[5]], pts[[2]], pts[[9]], pts[[10]], pts[[11]], pts[[12]],
     pts[[13]], pts[[6]], pts[[5]]}]
  (*  ,"MaxBoundaryCellMeasure"->densBound*)]
a3 = ToBoundaryMesh[
  Line[{pts[[9]], pts[[3]], pts[[14]], pts[[12]], pts[[11]], 
    pts[[10]], pts[[9]]}]
  (*  ,"MaxBoundaryCellMeasure"->densBound*)]
a4 = ToBoundaryMesh[
  Line[{pts[[8]], pts[[7]], pts[[6]], pts[[13]], pts[[12]], pts[[14]],
     pts[[4]], pts[[8]]}]
  (*  ,"MaxBoundaryCellMeasure"->densBound*)]

(*Join all BMeshes to only one BMesh*)
mesh1 = BoundaryElementMeshJoin[a1, a2];
mesh2 = BoundaryElementMeshJoin[mesh1, a3];
mesh3 = BoundaryElementMeshJoin[mesh2, a4];
mesh3["Wireframe"]

Thats what the boundary mesh looks like:

(*Create the Element mesh*)
DenseGrain = 5000;(*max.size of the elements*)
mesh = 
 ToElementMesh[mesh3, 
  "RegionMarker" -> {{{100, 100}, 1, DenseGrain}, {{2000, 500}, 2, 
     DenseGrain}, {{1000, 1500}, 3, DenseGrain}, {{3200, 1800}, 4, 
     DenseGrain}}, "MeshOrder" -> 2]

plotmesh = 
 mesh["Wireframe"[
   "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red], 
     FaceForm[Blue], FaceForm[Orange], FaceForm[LightBlue]}]]

Thats what the element mesh without refinement looks like:

The refinement should follow the grain boundaries. I am sure that there are more elegant ways to get the result, but I am stuck at this point and glad for all help or suggestions you can provide.

Thanks in advances

Max

Answer 1

The ratio of your largest feature to your smallest feature is about 200. Such a large feature ratio can often cause problems with automatic meshing. You can use DistanceMatrix to get an idea of the feature ratio.

DistanceMatrix[pts] // N // MatrixForm

The following approach is not likely to be robust or flexible, but it works on my machine. I set densBound=5 (i.e., the 1/4 the minimum feature distance) and also added a refinement zone about pts[[6]] (i.e., the point where the smallest features are).

Needs["FEMAddOns`"]
<< NDSolve`FEM`

(*Points*)
pts = {{0, 0}, {3333, 0}, {3333, 2166}, {0, 2166}, {1150, 0}, {1100, 
    900}, {1050, 950}, {0, 1333}, {3333, 1000}, {3066, 1500}, {2900, 
    1700}, {2800, 1750}, {1120, 900}, {3000, 2166}};

(*Create the BMeshes from pts*)
densBound = 5;(*Is intended to refine the boundaries,but is extremely \
time consuming and/or let ToElementMesh abort*)a1 = 
 ToBoundaryMesh[
  Line[{pts[[1]], pts[[5]], pts[[6]], pts[[7]], pts[[8]], pts[[1]]}], 
  "MaxBoundaryCellMeasure" -> densBound];
a2 = ToBoundaryMesh[
   Line[{pts[[5]], pts[[2]], pts[[9]], pts[[10]], pts[[11]], 
     pts[[12]], pts[[13]], pts[[6]], pts[[5]]}], 
   "MaxBoundaryCellMeasure" -> densBound];
a3 = ToBoundaryMesh[
   Line[{pts[[9]], pts[[3]], pts[[14]], pts[[12]], pts[[11]], 
     pts[[10]], pts[[9]]}], "MaxBoundaryCellMeasure" -> densBound];
a4 = ToBoundaryMesh[
   Line[{pts[[8]], pts[[7]], pts[[6]], pts[[13]], pts[[12]], 
     pts[[14]], pts[[4]], pts[[8]]}], 
   "MaxBoundaryCellMeasure" -> densBound];

(*Join all BMeshes to only one BMesh*)
mesh1 = BoundaryElementMeshJoin[a1, a2];
mesh2 = BoundaryElementMeshJoin[mesh1, a3];
mesh3 = BoundaryElementMeshJoin[mesh2, a4];
Show[mesh3["Wireframe"], 
 mesh3["Wireframe"["MeshElement" -> "PointElements"]]]
DenseGrain = 5000;(*max.size of the elements*)
mrf = 
 With[{rmf = 
    RegionMember[Disk[pts[[6]], 80]]}, 
  Function[{vertices, area}, 
   Block[{x, y}, {x, y} = Mean[vertices]; 
    If[rmf[{x, y}], area > DenseGrain/128, 
     area > DenseGrain]]]]; mesh = 
 ToElementMesh[mesh3, 
  "RegionMarker" -> {{{100, 100}, 1, DenseGrain}, {{2000, 500}, 2, 
     DenseGrain}, {{1000, 1500}, 3, DenseGrain}, {{3200, 1800}, 4, 
     DenseGrain}}, "MeshOrder" -> 2, MeshRefinementFunction -> mrf];

plotmesh = 
 mesh["Wireframe"[
   "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red], 
     FaceForm[Blue], FaceForm[Orange], FaceForm[LightBlue]}]]