Extract 2D quad mesh from 3D hexahedral mesh

Question

I would like to extract 2D mesh of outer surface of 3D meshed object.

Let's say I have 3D mesh data from here and I import the data into Mathematica.

    Needs["NDSolve`FEM`"];
    SetDirectory[NotebookDirectory[]];
    nodes3Dmesh = Import["nodes_3D_mesh.txt", "Table"];
    conn3Dmesh = Import["connection_3D_mesh.txt", "Table"] + 1;
    mesh3D = ToElementMesh["Coordinates" -> nodes3Dmesh, "MeshElements" -> 
    {HexahedronElement[conn3Dmesh]}];(*3D hex mesh*)

    Graphics3D[{ElementMeshToGraphicsComplex[mesh3D]}, Axes -> True, 
    AxesLabel -> {x, y, z}]

Now, I want to extract 2D mesh of surface lying in the x-z plane at y = 0 (this reddish surface in the above picture). My approach to this problem was:

    PosNodesY0 = Flatten[Position[nodes3Dmesh, _?(#[[2]] == 0. &)]];(*positions of nodes which 
    have y coordinate equal to zero*)
    conn2Dmesh = Select[Map[Select[#, MemberQ[PosNodesY0, #] &] &, conn3Dmesh], 
    UnsameQ[#, {}] &];(*computed connections for the 2D mesh*)
    mesh2D = ToElementMesh["Coordinates" -> Drop[nodes3Dmesh, None, {2}], 
    "MeshElements" -> {QuadElement[conn2Dmesh]}];(*resulting 2D quad mesh*)
    mesh2D["Wireframe"]

Questions:

1. My approach of extracting 2D mesh takes a long time for large 3D meshed object. Specifically the conn2Dmesh implementation. It is possible to make it faster?
2. My approach does not work for surface with more complicated geometry (e.g. white surface in the first picture). Would it be possible to generalize this for any outer surface of 3D object?

I would appreciate any help.

Answer 1

Here's the different approach (not sure how fast it will be though).

First, construct the mesh region and convert it to a boundary mesh:

mesh = MeshRegion[nodes3Dmesh, Hexahedron[conn3Dmesh]];
bmesh = BoundaryMesh[mesh];

Compute normal vectors of polygons:

enormal = Chop[Region`Mesh`MeshCellNormals[bmesh, 2]];

compare normals of adjacent polygons of lines and find corner edges:

mg = MeshConnectivityGraph[bmesh, {1, 2}];    
edges = VertexList[mg, {1, _}];    
adj = (AdjacencyList[mg, #] & /@ edges)[[All, All, 2]];    
corneredges = Pick[edges, Dot @@ enormal[[#]] & /@ adj, x_ /; x < 1];

HighlightMesh[bmesh, corneredges]

Partition the graph with corner edges and construct submeshes:

partition = WeaklyConnectedComponents[VertexDelete[mg, corneredges]];
meshes = 
  MeshRegion[MeshCoordinates[bmesh], MeshCells[bmesh, #]] & /@ 
   partition;

then project resulting meshes onto 2d:

projMesh[mesh_] :=
 Block[{coords, p, m},
  coords = MeshPrimitives[mesh, {2, 1}][[1]];
  p = First[coords];
  m = Select[
    Orthogonalize[Transpose[Transpose[Rest[coords]] - p], Dot], 
    AnyTrue[#, # != 0 &] &];
  m = AffineTransform[{m, -Dot[m, p]}];
  MeshRegion[m[MeshCoordinates[mesh]], MeshCells[mesh, 2]]
  ]

projMesh /@ meshes

Answer 2

Another possible approach. Use ConvexHullMesh to get the 7 surfaces and collect the polygons.

Clear[bmesh, chmesh, bmeshNormal, chmeshNormal, indexs, regs, meshs];
bmesh = BoundaryMeshRegion[mesh3D];
chmesh = ConvexHullMesh[bmesh];
bmeshNormal = Region`Mesh`MeshCellNormals[bmesh, 2];
chmeshNormal = Region`Mesh`MeshCellNormals[chmesh, 2];
indexs = 
  Table[Position[bmeshNormal, 
     x_?VectorQ /; 
      EuclideanDistance[x, chmeshNormal[[i]]] < .001], {i, 
     Length@chmeshNormal}][[;; , ;; , 1]];
regs = Table[
  RegionUnion@MeshPrimitives[bmesh, 2][[indexs[[i]]]], {i, 
   Length@chmeshNormal}];
meshs = 
 Table[reg = 
   TransformedRegion[regs[[i]], 
    If[Norm@Cross[chmeshNormal[[i]], {0, 0, 1}] > 0, 
     RotationTransform[{chmeshNormal[[i]], {0, 0, 1}}], Identity]]; 
  MeshRegion[
   MeshCoordinates[reg] /. {x_Real, y_Real, z_Real} :> {x, y}, 
   MeshCells[reg, 2]], {i, Length@chmeshNormal}]

Answer 3

Yet another way to use Region`Mesh`MeshCellNormals:

nodes3D = Import[$HomeDirectory <> "/Downloads/nodes_3D_mesh.txt", "Data"];
cons3D = Import[$HomeDirectory <> "/Downloads/connection_3D_mesh.txt", "Data"];

bmesh = BoundaryMesh[MeshRegion[nodes3D, Hexahedron[1 + cons3D]]]

faces = Values @ GroupBy[MeshCellIndex[bmesh, 2], 
   Round[Region`Mesh`MeshCellNormals[bmesh, #], 10^-4] &, 
   MeshRegion[MeshCoordinates @ bmesh, MeshCells[bmesh, #]] &]

Multicolumn[faces, 4]