6
$\begingroup$

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

Region Plot

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

2D Mesh Example

Any help would be much appreciated.

$\endgroup$
1
  • 1
    $\begingroup$ One of the nicer features in the 12.+ versions, is the OpenCascadeLink. OpenCascade is an open source mesher and does very good job snapping to sharp features (an area that is problematic for standard Region functions). Here is an example 223631 of where I used it (at the end is a quarter symmetry case). You should be able to use OpenCascade and a MeshRefinementFunction to achieve most of what you want in the newer versions. $\endgroup$ – Tim Laska Aug 26 '20 at 8:57
8
$\begingroup$

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

Boundary Mesh Region

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

Volume Mesh

$\endgroup$
1
  • $\begingroup$ Thanks Tim for the comments and work. You nailed exactly what I was seeking - only now I need to consider whether to invest in a MM upgrade. $\endgroup$ – Greenasnz Aug 27 '20 at 15:35

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.