# How can I make a region for boundary conditions for a 3D geometry?

I sometimes work with pipes which can take up complicated positions in 3D geometries. Here is a minimum working example of what I am trying to do.

Suppose I want to mesh and calculate vibration frequencies for this solid curved tube.

R0 = 1;
r = 0.05;
ϕ = 0.654 π;
g = Graphics3D[{CapForm["Butt"],
Tube[Table[{R0 Cos[θ], R0 Sin[θ], 0}, {θ,
0, ϕ, ϕ/25}], r]}, Axes -> True,
AxesLabel -> {"x", "y", "z"}]


I can make a mesh for my tube easily:

Needs["NDSolveFEM"];
reg = BoundaryDiscretizeGraphics[g,
MaxCellMeasure -> {"Area" -> 0.001} ]


Now I need a stress operator and boundary conditions. I give the stress operator below and add boundary conditions on one end of the tube which is easy. However, on the other end of the tube I also want to add boundary conditions which is difficult and the reason for my question. To get going I don't add boundary conditions on the "difficult" end.

ClearAll[stressOperator, u, v, w, x, y, z, Y, ν];
stressOperator[
Y_, ν_] := {Inactive[
Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0,
0}, {-Y/(2*(1 + ν)), 0, 0}} .
Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))),
0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}} .
Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0,
0}, {0, -Y/(2*(1 + ν)), 0}, {0, 0, -Y/(2*(1 + ν))}} .
Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -((Y*ν)/((1 -
2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)), 0}} .
Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -Y/(2*(1 + ν)),
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0,
0}} . Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0,
0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0,
0, -Y/(2*(1 + ν))}} .
Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}],
Inactive[Div][{{0, 0, 0}, {0,
0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 -
2*ν)*(1 + ν))), 0}} .
Inactive[Grad][v[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0,
0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}} .
Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,
0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}} .
Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}]};


Here is the code for the boundary conditions on the easy end;

{vals, funs} = NDEigensystem[{stressOperator[10^3, 33/100],
DirichletCondition[u[x, y, z] == 0,
y == 0 && (x - R0)^2 + z^2 <= r^2],
DirichletCondition[v[x, y, z] == 0,
y == 0 && (x - R0)^2 + z^2 <= r^2],
DirichletCondition[w[x, y, z] == 0,
y == 0 && (x - R0)^2 + z^2 <= r^2]},
{u, v, w}, {x, y, z} \[Element] mesh, 12];


Here are the vibration frequencies:

TableForm[Sqrt[vals], TableHeadings -> {Automatic, None}]


This is an example of a deflected vibration eigenmode:

n = 6;
fac = 0.02;
uif = funs[[n, 1]];
vif = funs[[n, 2]];
wif = funs[[n, 3]];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> fac];
Show[mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]],
Axes -> True, AxesLabel -> {"x", "y", "z"}]


So everything seems to be working and the issue is: How do I add a DirichletCondition boundary condition on the difficult end? I have to be able to specify the region. My first thought was to identify nodes, making the end region I wish to specify, and then make a MeshRegion. This is not straightforward, particularly in the general case, because you have to define a region to select the coordinates. So the most general question would be: How do you define a region of a mesh where you wish to add a boundary condition?

If you increase the resolution of your BoundaryDiscretizeGraphics, BoundaryElementMarkerUnion will separate your shape into three nice features as shown below.

Needs["NDSolveFEM"];
reg = BoundaryDiscretizeGraphics[g,
MaxCellMeasure -> {"Area" -> 0.001/2}];
mesh = ToElementMesh[reg];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = {ColorData["BrightBands"][#]} & /@ temp
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]


As you can see, the blue and yellow features represent the ends of the shape.

Now, you can refer to your DirichletCondition by ElementMarker like so:

{vals, funs} =
With[{em2 = (ElementMarker == 2), em3 = (ElementMarker == 3)},
NDEigensystem[{stressOperator[10^3, 33/100],
DirichletCondition[u[x, y, z] == 0, em2],
DirichletCondition[v[x, y, z] == 0, em2],
DirichletCondition[w[x, y, z] == 0, em2],
DirichletCondition[u[x, y, z] == 0, em3],
DirichletCondition[v[x, y, z] == 0, em3],
DirichletCondition[w[x, y, z] == 0, em3]}, {u, v,
w}, {x, y, z} \[Element] mesh, 12]
];


When you plot the solution, you can see that both ends are fixed.

n = 6;
fac = 0.02;
uif = funs[[n, 1]];
vif = funs[[n, 2]];
wif = funs[[n, 3]];
dmesh = ElementMeshDeformation[mesh, {uif, vif, wif},
"ScalingFactor" -> fac];
Show[mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]],
Axes -> True, AxesLabel -> {"x", "y", "z"}]


• Excellent solution. I will have to explore the markers more. Thanks.
– Hugh
Commented Jul 7, 2021 at 13:44
• Is BoundaryElementMarkerUnion Written up anywhere? I can find a little under ElemetMesh. I notice that as you change the area of the mesh elements the number of features you get changes. How does the code look for features? Is it to do with change in normal direction? Should I ask a new question on this?
– Hugh
Commented Jul 7, 2021 at 14:30
• @Hugh Thank you! They use it quite a bit in the OpenCascadeLink Tutorial. That is where I learned about it. Commented Jul 7, 2021 at 14:56
• I can't resist writing that I find both the question and this answer quite impressive. I imagine they must have involved a lot of time and persistence. Commented Jul 14, 2021 at 4:53

As Tim points out, using markers is the way to go. Here is another way to generate that torus:

Needs["OpenCascadeLink"]
OpenCascadeTorus[{{0, 0, 0}, {0, 0, 1}}, 1, 0.05, 0.654 \[Pi]]];
"ShapeSurfaceMeshOptions" -> {(*"AngularDeflection"->0.5,*)
"LinearDeflection" -> 0.001}];
bmesh["Wireframe"]


You can then visualize the mesh with its markers

groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
Show[
bmesh["Edgeframe"],
bmesh["Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> (Directive[EdgeForm[], FaceForm[#]] & /@
colors)]]
, Epilog -> Inset[LineLegend[colors, groups], Scaled[{0.85, 0.8}]]
]
`

• Thanks useful as usual. There is a lot to learn here. Is the OpendCascade Tutorial the only place to learn from?
– Hugh
Commented Jul 7, 2021 at 21:31
• @Hugh, yes, currently (V12.3) the tutorial is the only place. Commented Jul 8, 2021 at 4:54
• (its = possessive, it's = "it is" or "it has". See for example youtube.com/watch?v=8Gv0H-vPoDc&t=1m20s and wikihow.com/Use-Its-and-It%27s) Commented Jul 8, 2021 at 8:51
• @PeterMortensen, fixed. thanks. Commented Jul 8, 2021 at 14:44