Mathematica Asked by Andrew Green on June 27, 2021
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.
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"]
Answered by Tim Laska on June 27, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP