4
$\begingroup$

I am trying to create a finer mesh in a subregion of a solid cylinder in order to improve the resolution there. I don't want to have too many elements for the entire cylinder because it slows down the computation. From other questions here and the tutorial on Element Mesh Creation https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html, I see that this can be done for 2D shapes and composite shapes (different materials).

  1. I'd like a finer mesh along the axis of the cylinder, with a smaller radius (e.g. 1/3 of the radius).

I suppose I can break it in two separate cylinders--the inner solid cylinder and the outer "pipe" and then combine them with RegionUnion. But, then I am not sure how to define the boundary surfaces for the boundary conditions (outer cylindrical surface and the whole end faces). I hope there is a more elegant solution.

  1. How to create finer mesh just for the end faces? The boundary surfaces are defined with ElementMarker. I tried defining separate boundary mesh for each face and then turning it ToElementMesh, but it didn't work. Below is the code for the whole solid cylinder.
OpenCascadeShapeType[cyl]; 
cyl = OpenCascadeShape[c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003]];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[cyl]
mesh = ToElementMesh[bmesh]
Show[mesh["Wireframe"], Graphics3D[c1]]

When I visualize the solution with SliceContourPlot3D, the resolution is insufficient to capture small features near the axis.

enter image description here

$\endgroup$
1
  • $\begingroup$ You'll probably need to find an appropriate MeshRefinementFunction. Perhaps start by searching the site with that keyword. $\endgroup$
    – MarcoB
    Jul 17, 2023 at 0:17

1 Answer 1

8
$\begingroup$

Version 1

One way to do it could be to generate a 2D disk mesh with an appropriate refinement and the extrude that to 3D.

(m2 = ToElementMesh[Disk[{0, 0}, 0.003], 
    MeshRefinementFunction -> 
     Function[{vertices, area}, 
      area > 0.00000005 (0.001 + 
          500 Norm[Mean[vertices]])]])["Wireframe"]

enter image description here

mesh = ElementMeshRegionProduct[m2, 
  ToElementMesh[Line[{{0}, {0.01}}]]]


mesh["Wireframe"]

enter image description here

Here is a view inside:

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{-0.003, 0}, All, {0, 0.005}}]]

enter image description here

Version 2

A variant of this approach is to use a graded mesh for the element mesh region product, like so:

m1 = ToGradedMesh[{Line[{{0}, {0.01}}], <|
     "Alignment" -> "BothEnds"|>}];
MeshRegion[m1]

enter image description here

mesh = ElementMeshRegionProduct[m2, m1];

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{-0.003, 0}, All, {0, 0.005}}]]

enter image description here

Note, how the mesh is becoming denser towards the bottom (and top not shown) Here is a view of the boundaries:

mesh["Wireframe"[
  "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red], 
    FaceForm[Blue]}, PlotRange -> {{-0.003, 0}, All, {0, 0.01}}]]

enter image description here

Version 3

Here is a version that has a refined boundary but a coarse interior:

Needs["NDSolve`FEM`"]
c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
mesh = ToElementMesh[c1, "MaxBoundaryCellMeasure" -> 0.1, 
   "MaxCellMeasure" -> Infinity];
Show[mesh["Wireframe"], Graphics3D[c1]]

enter image description here

Cross section:

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "ElementMeshDirective" -> 
   Directive[EdgeForm[Black], FaceForm[LightGray]], 
  PlotRange -> {{0, 0.005}, {0, 0.003}, {0, 0.003}}]]

enter image description here

Version 4

Yet another way to do it would be by using OpenCascade and create an inner cylinder that is refined to a different level.

Needs["NDSolve`FEM`"]
Needs["OpenCascadeLink`"]
c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
c2 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003/2];

s1 = OpenCascadeShape[c1];
s2 = OpenCascadeShape[c2];

faces = Flatten[OpenCascadeShapeFaces /@ {s1, s2}];
union = OpenCascadeShapeUnion[faces];

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];

bmesh["Wireframe"]

enter image description here

Create to material markers p1 and p2:

p1 = {0.01/2, 0, 0};
p2 = {0.01/2, 0, 0.005/2};
Show[
 Graphics3D[{{Red, PointSize[0.02], Point[p1]}, {Green, 
    PointSize[0.02], Point[p2]}}],
 bmesh["Wireframe"],
 Graphics3D[{Opacity[0.2], c1}]]

enter image description here

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{p1, 1, 10^-12}, {p2, 2, 10^-10}}];
Show[mesh["Wireframe"], Graphics3D[c1]]

enter image description here

Version 5

Here is a version with a MeshRefinementFunction. By testing I found that the mesh changes when a "MaxCellMeasure" around 10^-11 was given. So I create a function that re-scales around this value from 0 to the boundary at 0.003.

Needs["NDSolve`FEM`"]
rescale = Rescale[Sqrt[y^2 + z^2], {0, 0.003}, {10^-14, 10^-11}]

The Sqrt[y^2 + z^2] is just the radial distance from the x-axis. This I then injected into the mesh refinement function:

c1 = Cylinder[{{0, 0, 0}, {0.01, 0, 0}}, 0.003];
cf = With[{r = rescale},
   Compile[{{coordinates, _Real, 2}, {vol, _Real, 0}},
    Block[{com, x, y, z, rvol},
     com = Mean[coordinates];
     {x, y, z} = com;
     rvol = r;
     If[vol > rvol, True, False]
     ]
    ]
   ];
mesh = ToElementMesh[c1, MeshRefinementFunction -> cf]
Show[mesh["Wireframe"], Graphics3D[c1]]

enter image description here

$\endgroup$
2
  • $\begingroup$ how would you define the boundaries--cylindrical surface, top-end face, and bottom-end face? I tried mesh["Wireframe"[ "MeshElementStyle" -> {FaceForm[Green], FaceForm[Red], FaceForm[Blue]}, PlotRange -> {{-0.003, 0}, All, {0, 0.005}}]] but it gives a series of cross-sections along the length (instead of top and bottom). In addition, I cannot assign names to the boundaries in order to use them in the boundary conditions for solving the PDE. $\endgroup$
    – user93372
    Jul 17, 2023 at 22:08
  • $\begingroup$ @user93372, I hope i have addressed the first part of your question, I do not understand what you mean with 'assign names' though. $\endgroup$
    – user21
    Jul 18, 2023 at 18:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.