3D FEM Meshing Internal Regions

Question

I am trying to set up a 3D element mesh with intersecting regions having different mesh densities. I am having difficulty setting up the defining boundary meshes from which I will then apply ToElementMesh. I understand how to do it in 2D but I do not know the best way to do it for 3D. My code below has been cut down to try and show the basic problem I have. I need to set up the boundary mesh on the green problem volume so the intersections with the "e-core" region on the x=z=0 axis can be meshed consistent with the finer mesh to be used in the e-core region volume. Although I have shown the full core, due to symmetry in the problem I will only use 1/4 of it, i.e, that which intersects with the green volume.

Please note I only have MM 10.4, so I do not have access to FEMAddons. However, I would also be interested to see how it could be done if I upgraded in the future.

Clear["Global`*"];
Needs["NDSolve`FEM`"];

eCore[cw_, ch_, cd_, ww_, wh_] := 
  Module[(*cw = core width, ch = core height, cd = core depth, www = 
   window width, w = window height*){vertices, topFace, reg},
   vertices = {{-cw/2, 0}, {-cw/4 - ww/2, 0}, {-cw/4 - ww/2, 
      wh}, {-cw/4 + ww/2, wh}, {-cw/4 + ww/2, 0}, {cw/4 - ww/2, 
      0}, {cw/4 - ww/2, wh}, {cw/4 + ww/2, wh}, {cw/4 + ww/2, 
      0}, {cw/2, 0}, {cw/2, ch}, {-cw/2, ch}};
   topFace = 
    BoundaryMeshRegion[vertices, 
     Line[{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1}]];
   reg = RegionProduct[topFace, 
     MeshRegion[{{-ch/2}, {ch/2}}, Line[{1, 2}]]]; reg];

(*Create an e-core using above function and rotate/translate position \
as required*)
regCore1 = 
  TransformedRegion[
   TransformedRegion[eCore[0.065, 0.033, .027, .013, .022], 
    RotationTransform[0, {0, 0, 1}]], 
   TranslationTransform[{0, 0.002, 0}]] ;
bmeshCore1 = 
  BoundaryDiscretizeRegion[regCore1, 
   MaxCellMeasure -> {"Length" -> 0.005}, Axes -> True, 
   AxesLabel -> {x, y, z}];
(*get coordinates of 1/4 core1 mesh in problem volume*)
core1Coord = 
  Cases[DeleteDuplicates[MeshCoordinates[bmeshCore1]], {x_, y_, z_} /;
     x \[GreaterSlantEqual] 0 && z \[LessSlantEqual] 0];

(*Create air region that defines the problem boundaries allowing for \
symmetry in the problem*)
radiusAir = 0.15;
regAir1 = 
  RegionIntersection[
   Cuboid[{0, 0, -radiusAir}, {radiusAir, radiusAir, 0}], 
   Ball[{0, 0, 0}, radiusAir]];
bmeshAir1 = 
  BoundaryDiscretizeRegion[regAir1, 
   MaxCellMeasure -> {"Length" -> 0.01}, Axes -> True, 
   AxesLabel -> {x, y, z}];
RegionPlot3D[{regCore1, regAir1}, Axes -> True, 
 AxesLabel -> {x, y, z}, PlotStyle -> {Blue, Green}]

I guess I want the 3D equivalent of the Wolfram 2D example given under Element Mesh Generation. Here I have modified it to have a higher mesh density on the internal line boundary.

(*2D Example of open line boundary within a closed rectangular \
boundary - modified from Wolfram FEM Meshing example*)
n = 20; \
lineCoord = 
 DeleteDuplicates[
  Join[Table[{1/6. + (i - 1)*4/(6.*(n - 1)), 1/6.}, {i, 1, n}], 
   Table[{5/6., 1/6. + (i - 1)*4/(6.*(n - 1))}, {i, 1, n}]]];
bmesh = ToBoundaryMesh[
   "Coordinates" -> Join[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, lineCoord],
    "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        1}}], LineElement[
      Partition[Delete[Last[FindShortestTour[lineCoord]], 1], 2, 1] + 
       4]}];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> {"Length" -> 0.5}];
mesh["Wireframe"]

Any help would be much appreciated.

Answer 1

Here is an answer based version 12.1.1 for Microsoft Windows (64-bit) (June 19, 2020) as alluded to in the comments.

Here is the workflow to create Computational Solid Geometry (CSG) with OpenCascadeLink:

Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]
(* Geometry Parameters *)
{cw, ch, cd, ww, wh} = {0.065, 0.033, .027, .013, .022};
yoff = 0.002;
radiusAir = 0.15;
(* Use CSG to Create Core Shape *)
shape0 = OpenCascadeShape[
   Cuboid[{-cw/2, 0 + yoff, -cd/2}, {cw/2, ch + yoff, cd/2}]];
shape1 = OpenCascadeShape[
   Cuboid[{-cw/4 - ww/2, 0 + yoff, -cd/2}, {-cw/4 + ww/2, wh + yoff, 
     cd/2}]];
shape2 = OpenCascadeShape[
   Cuboid[{cw/4 - ww/2, 0 + yoff, -cd/2}, {cw/4 + ww/2, wh + yoff, 
     cd/2}]];
core = OpenCascadeShapeDifference[shape0, shape1];
core = OpenCascadeShapeDifference[core, shape2];
(* Create Air Sphere *)
shapea = OpenCascadeShape[Ball[{0, 0, 0}, radiusAir]];
(* Create Quarter Symmetry *)
(* Create Quarter Symmetry Cube *)
shapeq = OpenCascadeShape[
   Cuboid[{0, 0, -radiusAir}, {radiusAir, radiusAir, 0}]];
(* Create Quarter Symmetry Regions *)
shapeinta = OpenCascadeShapeIntersection[shapeq, shapea];
shapeintcore = OpenCascadeShapeIntersection[shapeq, core];
(* Create Shape with Internal Boundaries *)
(* https://wolfram.com/xid/0bxz9t5u18ulek5jqypwwj4nro1wg77bu-xj0w1m*)


union = OpenCascadeShapeUnion[shapeinta, shapeintcore];
intersection = OpenCascadeShapeIntersection[shapeinta, shapeintcore];
shape = OpenCascadeShapeSewing[{union, intersection}];
(* Create Boundary Mesh *)
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[shape];
(* Visualize Surfaces *)
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = {Opacity[0.75], ColorData["BrightBands"][#]} & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]

Now, we can set up a refinement region based on the core and create a volume mesh like so:

(* Define Core as Refinement Region *)
refinementRegion = 
  MeshRegion@
   ToElementMesh[
    OpenCascadeShapeSurfaceMeshToBoundaryMesh[shapeintcore], 
    MaxCellMeasure -> Infinity];
(* Create Mesh Refinement Function *)
mrf = With[{rmf = RegionMember[refinementRegion]}, 
   Function[{vertices, volume}, 
    Block[{x, y, z}, {x, y, z} = Mean[vertices]; 
     If[rmf[{x, y, z}], volume > 1.25`*^-7/8^2, 
      volume > 1.0`*^-6/8]]]];
(* Create and Display Volumetric Mesh *)
(mesh = ToElementMesh[bmesh, 
    MeshRefinementFunction -> mrf])["Wireframe"]