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.
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 aMeshRefinementFunction
to achieve most of what you want in the newer versions. $\endgroup$ – Tim Laska Aug 26 '20 at 8:57