# FEM in 3D: assigning proper mesh for NDSolveValue

I am trying to compute electric potential in a simple 3D geometry containing (in a simplied version) 3 conductive rings (2 outer rings at 0 V, and the inner at 1 V) in different planes and I don't get a proper mesh. I define the rings ok and can visualize a mesh for all of them and also I define an Air region around them to compute the potential in all this region. But, if I try to define a region which is the Difference between the Air-region and the Cylinders-region, it does not generate a proper mesh. Here is my trial. Defining the geometry:

ClearAll["Global*"]
Needs["NDSolveFEM"];
(* Define regions *)
cylin4KLength = 1.0;
cylin4KOffset = 3.0;
cylin40KLength = 1.0;
cylin40KOffset = -3.0;
cylinAnelLength = 1.0;
cylinAnelOffset = 0.0;

air = Cylinder[{{0, 0, -2 cylin4KOffset}, {0, 0, + 2 cylin4KOffset}},

cylin4Kout =
Cylinder[{{0, 0, -cylin4KLength/2 + cylin4KOffset}, {0,
cylin4Kin =
Cylinder[{{0, 0, -cylin4KLength/2 + cylin4KOffset}, {0,
cylin4K = RegionDifference[cylin4Kout, cylin4Kin];

cylin40Kout =
Cylinder[{{0, 0, -cylin40KLength/2 + cylin40KOffset}, {0,
cylin40Kin =
Cylinder[{{0, 0, -cylin40KLength/2 + cylin40KOffset}, {0,
cylin40K = RegionDifference[cylin40Kout, cylin40Kin];

cylinAnelout =
Cylinder[{{0, 0, -cylinAnelLength/2 + cylinAnelOffset}, {0,
cylinAnelin =
Cylinder[{{0, 0, -cylinAnelLength/2 + cylinAnelOffset}, {0,
0, +cylinAnelLength/2 + cylinAnelOffset}}, (cylinAnelRadius - 2)];
cylinAnel = RegionDifference[cylinAnelout, cylinAnelin];

cylins = RegionUnion[cylin4K, cylin40K, cylinAnel];


Visualizing this geometry:

bmeshair = ToBoundaryMesh[air]["Wireframe"];
bmeshCylin4K = ToBoundaryMesh[cylin4K]["Wireframe"];
bmeshCylin40K = ToBoundaryMesh[cylin40K]["Wireframe"];
bmeshCylinAnel = ToBoundaryMesh[cylinAnel]["Wireframe"];
Show[bmeshair, bmeshCylin4K, bmeshCylin40K, bmeshCylinAnel]


Alternatively (this is good for visualizing the boundary conditions for the Potential in NDSolveValue), the rings could be defined as:

regionGroundPlates =
ImplicitRegion[((x^2 + y^2 > cylinRadiusIn^2) && (x^2 + y^2 <
cylin40KRadius^2)) && ((z - cylin4KOffset)^2 < (cylin4KLength/
2)^2  || (z - cylin40KOffset)^2 < (cylin40KLength/2)^2), {x,
y, z}];
regionAnel =
ImplicitRegion[((x^2 + y^2 > (cylinAnelRadius - 2)^2) && (x^2 +
y^2 < cylinAnelRadius^2)) && ((z -
cylinAnelOffset)^2 < (cylinAnelLength/2)^2), {x, y, z}];
RegPlates = RegionUnion[regionAnel, regionGroundPlates];
Region[RegPlates]


Now to the mesh: if I try to define a mesh from a RegionDifference as in:

regMesh = RegionDifference[Air, cylins];
bmeshtotal = ToBoundaryMesh[regMesh]["Wireframe"];
Show[bmeshtotal]


I do get and error: "...is not a type of graphics". Once I get a proper mesh (that's my problem), I'd like to use NDSolveValue to solve for the potential using the boundary conditions (OV in the outer rings and 1V in the inner one) using:

eq = Laplacian[u[x, y, z], {x, y, z}]; Vanel = 1;
bc = {DirichletCondition[
u[x, y, z] ==
0, ((x^2 + y^2 > cylinRadiusIn^2) && (x^2 + y^2 <
cylin4KOffset)^2 < (cylin40KLength/2)^2  || (z -
cylin40KOffset)^2 < (cylin40KLength/2)^2) ],
DirichletCondition[
u[x, y, z] ==
Vanel, ((x^2 + y^2 > (cylinAnelRadius - 2)^2) && (x^2 + y^2 <
cylinAnelOffset)^2 < (cylinAnelLength/2)^2)  ]};
solMV = NDSolveValue[{eq == 0, bc}, u, {x, y, z} \[Element] regMesh];


Can someone help me with this mesh? Claudio

Use something like this to generate the mesh:

object = Fold[RegionDifference, {air, cylin4K, cylin40K, cylinAnel}];
(mesh = ToElementMesh[object,
"RegionHoles" -> {{30, 0, 3}, {30, 0, -3}, {11, 0, 0}},
"MaxCellMeasure" -> Infinity])["Wireframe"]


I have added region holes, as I think you will want to set boundary condition inside the air enclosed geometry on those rings.

You can either use predicates to address the (presumably) Dirihclet conditions you wish to set. You can also use ElementMarker to do so.

Here is a way to scroll through the surfaces to get the ElementMarker for NeumannValue

bmesh = ToBoundaryMesh[mesh];
bIDs = bmesh["BoundaryElementMarkerUnion"];
edgeframe = bmesh["Edgeframe"];
outline = bmesh["Wireframe"["MeshElementStyle" ->
Directive[Opacity[0.2], FaceForm[LightBlue], EdgeForm[]]]];
Manipulate[Show[
outline,
edgeframe,
bmesh["Wireframe"[ElementMarker == bIDs[[id]],
"MeshElementStyle" -> Directive[FaceForm[Green], EdgeForm[]]]]
], {{id, 1, "ElementMarker ID"}, 1, Length[bIDs], 1,
Appearance -> "Open"}, SaveDefinitions -> True]


A similar code can also be used for DirichletCondition, for that I'd refer you to the documentation: Paste this into the documentation system FEMDocumentation/tutorial/ElementMeshVisualization#2089232889` to get to the relevant section.

• Excellent, thanks user21. When I first tried your solution above, it didn't work: gave the same problem of an infinite boundary. Then I decided to upgrade my Mathematica from 12.0 to the 12.3 version and then it worked beautifully. Nov 23, 2021 at 17:46