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

Question

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:

<< NDSolve`FEM`
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).

Answer 1

Here is a way to do it:

<< NDSolve`FEM`
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"]