3
$\begingroup$

In a finite element code I whish to use "IncludePoints" in ToElementMesh and I get a failure. Here is a minimum working example.

Needs["NDSolve`FEM`"];
Needs["OpenCascadeLink`"];
ClearAll[u, v, w];
Lbeam = 0.15; (*Length of beam *)
wth = 0.04;  (*Width of beams*)
tk = 0.008; (*Thickness of beams*)
rhole = 11/(1000. 2); (* Radius of hole*)

beam = Cuboid[{-Lbeam/2, -wth/2, 0}, {Lbeam/2, wth/2, tk}];
hole = Cylinder[{{0, 0, 0}, {0, 0, tk}}, rhole];
geom = RegionDifference[beam, hole];
rcp = 0.019;
extraPoints = 
  Table[{rcp Cos[\[Theta]], rcp Sin[\[Theta]], 0}, {\[Theta], 0, 
    2 \[Pi] - (2 \[Pi])/32, (2 \[Pi])/32}];
mesh = ToElementMesh[geom, "IncludePoints" -> extraPoints];

The error message is

TetGenTetrahedralize::reterr: Tetrahedralize returned an error, 3.

ToElementMesh::femtemnm: A mesh could not be generated.

For a mesh with quads this works. I tried running the code in this post and got the same failure.

Version 14.3 for Windows.

I would like to do this in OpenCascade but this does not seem to be an option.

Is there a workaround? Thanks for any help.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
6
$\begingroup$

Currently "IncludePoints" can not be on the boundary, they need to be strictly inside the region. And we can see that your include points are, in fact, on the boundary:

bmesh = ToBoundaryMesh[geom];
Show[
 bmesh["Edgeframe"],
 Graphics3D[{Red, PointSize[0.02], Point[extraPoints]}],
 bmesh["Wireframe"[
   "MeshElementStyle" -> 
    Directive[Opacity[0.2], FaceForm[Blue], EdgeForm[]]]]
 ]

enter image description here

There may be cases where it works to put them on the boundary but currently, that's not generally possible. This is explained in this section of ToBoundaryMesh. ToElementMesh calls ToBoundaryMesh for that part of the job and that is explained here.

That being said, what can you do about it? Probably your best bet is to use OpenCascade for that. The idea is the following: OpenCascade can generate boolean operations like the union of the shapes given to it are of the same dimensional. We extract the faces of the geometry geom and add additional faces of the based on the include points.

Get the faces of the geom:

of = OpenCascadeShapeFaces[OpenCascadeShape[geom]];

Next, we construct a crude 'circle' with the extra points in a z==0 plane and then extrude that to until the z==tk plane.

st = FindShortestTour[extraPoints][[2]];
ols = OpenCascadeShape[Line[extraPoints[[st]]]];
sweep = OpenCascadeShapeLinearSweep[ols, {{0, 0, 0}, {0, 0, tk}}];

Now, we find the union of all of those faces:

ns = OpenCascadeShapeUnion[Flatten[{of, sweep}]];

And we can generate the boundary mesh:

(bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[ns])["Wireframe"]

enter image description here

And that can then be meshed:

mesh = ToElementMesh[bmesh]
(* ElementMesh[{{-0.075, 0.075}, {-0.02, 0.02}, {0., 
   0.008}}, {TetrahedronElement["<" 16374 ">"]}] *)

Hope that helps.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Works perfectly. Always a good and helpful answer. $\endgroup$ Commented Feb 13 at 13:03

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.