# How to generate 3D FEM mesh with inner points at the boundary

How to generate a mesh which contains given points at the boundary. The simplest representation would be a cube with a circle presented by this code:

<< NDSolveFEM
cubeCoordinates = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0,
1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
circle = {#[[1]], #[[2]], 0} & /@ CirclePoints[{0.5, 0.5}, 0.25, 10];
ListPointPlot3D[{cubeCoordinates, circle}]
QE = QuadElement[{{1, 2, 3, 4}, {5, 6, 7, 8}, {1, 2, 6, 5}, {2, 3, 7,
6}, {3, 4, 8, 7}, {4, 1, 5, 8}}];
bm = ToBoundaryMesh["Coordinates" -> cubeCoordinates,
"BoundaryElements" -> {QE}];
bm["Wireframe"]
em = ToElementMesh[bm];
em["Wireframe"]


Generating the mesh without including the circle is very simple, however when trying to include the circle points I always get an error. Also nothing about such case is in documentation (only 2D examples are present).

• You could use dm = DelaunayMesh[ Cases[Join[cubeCoordinates, circle], {_, _, 0}][[All, {1, 2}]]] as a basis for the lower part. Dec 1, 2016 at 17:53

Here is a way to do it:

<< NDSolveFEM
cubeCoordinates = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0,
1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
circle = {#[[1]], #[[2]], 0} & /@ CirclePoints[{0.5, 0.5}, 0.25, 50];
(*ListPointPlot3D[{cubeCoordinates,circle}]*)


Generate a Delaunay triangulation:

dm = DelaunayMesh[
Cases[Join[cubeCoordinates, circle], {_, _, 0}][[All, {1, 2}]]];


And generate the boundary mesh manually.

baseCoords = ArrayPad[MeshCoordinates[dm], {0, 1}, 0.];
restCoords = N[Cases[cubeCoordinates, Except[{_, _, 0}]]];
n = Length[baseCoords];
QE = QuadElement[{n + {1, 2, 3, 4}, {1, 2, 2 + n, 1 + n}, {2, 3,
3 + n, 2 + n}, {3, 4, 4 + n, 3 + n}, {4, 1, 1 + n, 4 + n}}];
bm = ToBoundaryMesh[
"Coordinates" -> Join[baseCoords, restCoords],
"BoundaryElements" ->
Join[MeshCells[dm, 2, Multicells -> True] /.
Polygon -> TriangleElement, {QE}]];
bm["Wireframe"]


em = ToElementMesh[bm];
em["Wireframe"]