8
$\begingroup$

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}]

enter image description here

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"]

enter image description here

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.

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
12
$\begingroup$

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]

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
5
$\begingroup$

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}]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for your answer. I have problem of reproducing your results. I got error RegionUnion::argm: RegionUnion called with 0 arguments; 1 or more arguments are expected. $\endgroup$
    – Moonwalk
    Jul 22, 2022 at 8:24
1
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.