# Identify the Vertices and Lines around the edges of a free surface - Mesh/DiscreteGraphics

To write a mesher that closes a mesh surface using an implied surface from a boundary spline curve and mesh generated using DiscreteGraphics, I need to identify the Vertices and Lines around the edges of a free surface. Is there a simple way to extract just this data using MeshCells?

Below is a sample mesh with two free edges/openings that I would like to extract this information from.

curv1 = {{3, 0, 0}, {1, 1, 0}, {0, 2, 0}, {-2, 0, 0}, {0, -2, 0}, {3, 0, 0}};
curv2 = {{2, 0, 5}, {1, 1, 2}, {0, 2, 2}, {-1, 0, 5}, {0, -2, 5}, {3, 0, 5}};
sur1 = BSplineSurface[{curv1, curv2}, SplineClosed -> {False, True},     SplineDegree -> 3];
DiscretizeGraphics[sur1] The presence of a seam in both @flinty's and @kglr's answer suggest that an invalid mesh (for FEM purposes) is being created by the OP's DiscretizeGraphics approach. Rather than troubleshooting the meshing approach, I will present a structured meshing approach that eliminates the spurious edge artifact.

First, we will use a BSplineFunction to map a structured UV map to the curved surface.

curv1 = {{3, 0, 0}, {1, 1, 0}, {0, 2, 0}, {-2, 0, 0}, {0, -2, 0}, {3,
0, 0}};
curv2 = {{2, 0, 5}, {1, 1, 2}, {0, 2, 2}, {-1, 0, 5}, {0, -2, 5}, {3,
0, 5}};
bsf1 = BSplineFunction[{curv1, curv2}, SplineClosed -> {False, True},
SplineDegree -> 3];


The following workflow will create structure MeshRegion:

(* Import Required Package *)
Needs["NDSolveFEM"]
(* Create a UV Tensor Product Grid *)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
seg = pointsToMesh@Subdivide[0, 1, 36]
rp = RegionProduct[seg, seg]
(* Extract Coords from RegionProduct *)
crd = MeshCoordinates[rp];
(* Map coordinates to BSPlineFunction *)
crd2 = crd /. {{x_, y_} -> Chop[bsf1[x, y], 1*^-7]};
(*grab incidents RegionProduct mesh*)
inc = Delete /@ MeshCells[rp, 2];
(* Convert Quads to Triangles *)
inc2 = Partition[
Flatten[ inc /. {{i_, j_, k_, l_} -> {{i, j, k}, {i, k, l}}}], 3];
mrkrs = ConstantArray[1, First@Dimensions@inc2];
(* FEM Create BoundaryMesh *)
bm = ToBoundaryMesh["Coordinates" -> crd2,
"BoundaryElements" -> {TriangleElement[inc2, mrkrs]}];
(* Convert BoundaryMesh to MeshRegion *)
mr = MeshRegion[bm];
HighlightMesh[mr, Style[1, Orange]] The resulting mesh looks pretty good.

Now, we can apply @kglr's approach to see that we removed the spurious edge:

(* Apply kglr's Edge Extraction Method *)
boundaryedgeindices =
Flatten@Position[
HighlightMesh[mr, Style[{1, boundaryedgeindices}, Thick, Red]] # Update: Top Surface 2D Mesh

In the comments, the OP had a question about capping the ends of the mesh. Because the projected curve is not convex, a simple capping is not generally possible. One possibility, is to create a minimal surface.

The following workflow shows how to create a 2D mesh with nodes that are equivalenced with the 3D mesh by setting the Mesh Order to 1 and the SteinerPoints option to False.

(* Extract Coords from segment *)
crd2d = MeshCoordinates[seg];
(* Map coordinates to BSPlineFunction *)
crd2d2 = Flatten[
crd2d /. {{x_} :> Chop[{bsf1[1, x][[1 ;; 2]]}, 1*^-7]}, 1];
(* grab incidents segmentr mesh *)
inc2d = Delete /@ MeshCells[seg, 1];
(* Create Boundary Mesh *)
bm2d = ToBoundaryMesh["Coordinates" -> crd2d2,
"BoundaryElements" -> {LineElement[inc2d]}];
bm2d["Wireframe"]
(* Create 2D element mesh *)
m2d = ToElementMesh[bm2d, "MeshOrder" -> 1, "SteinerPoints" -> False];
m2d["Wireframe"]


To use the Mathematica example to create a minimal surface will require a little thought to specify the DirichletCondition because the curve is specified parametrically. Since the OP has routines for calculating minimal surfaces, I will not go into it here. # Update: Minimal Surface

The following workflow will solve for the minimal surface using NDSolveValue.

(* convert bsf1 to x,y,z components *)
ztop0[v_?NumericQ] := Module[{val}, val = bsf1[1, v]; Last@val];
zmin = First@NMinimize[ztop0[t], {t, 0, 1}];
zmax = First@NMaximize[ztop0[t], {t, 0, 1}];
zmid = Mean[{zmin, zmax}];
xtop[v_?NumericQ] := Module[{val}, val = bsf1[1, v]; First@val];
ytop[v_?NumericQ] := Module[{val}, val = bsf1[1, v]; val[]];
ztop[v_?NumericQ] := Module[{val}, val = bsf1[1, v]; Last@val - zmid];
(* Use Nearest to find v given x,y *)
nf = Nearest[
Table[{xtop[t], ytop[t]}, {t, 0, 1, .0001}] ->
Table[t, {t, 0, 1, .0001}]];
(* calculate z given x,y for DirichletCondition *)
fz[x_?NumericQ, y_?NumericQ] := ztop[First@nf[{x, y}]]
(* Minimal Surface https://wolfram.com/xid/0bdpx7hg6-hvook1 *)
ufun = NDSolveValue[{-Inactive[Div][(1/Sqrt[1 + \!$$\*SubscriptBox[\(∇$$, $${x, y}$$]$$u[x, y]$$\).\!$$\*SubscriptBox[\(∇$$, $${x, y}$$]$$u[x, y]$$\)]) Inactive[Grad][
u[x, y], {x, y}], {x, y}] == 0,
DirichletCondition[u[x, y] == fz[x, y], True]},
u, {x, y} ∈ m2d];


Now, we can convert the 2D mesh to a 3D boundary mesh using the minimal surface solution for the z coordinates:

(* create and display minimal surface boundary mesh *)
c3d = Join[m2d["Coordinates"], List /@ (ufun["ValuesOnGrid"] + zmid),
2];
bmtop = ToBoundaryMesh["Coordinates" -> c3d,
"BoundaryElements" -> m2d["MeshElements"]];
Show[bm["Wireframe"["MeshElementStyle" -> {FaceForm[Green]}]],
bmtop["Wireframe"["MeshElementStyle" -> {FaceForm[Red]}]]] You can see that the free surface nodes align well with the base mesh.

• That looks great, Tim! Do you know how to close the ends of this irregular tube with a mesh? – mh2000 Oct 17 at 17:45
• @mh2000 Sorry for the delay. I was at my niece's wedding yesterday. I don't think there is a unique way to cap the surfaces or a unique surface for that matter. If the enclosing curve is convex, we could attempt a simple cap with the centroid. We could attempt a minimal surface as shown here. – Tim Laska Oct 18 at 19:05
• Thanks, Tim! That is exactly what I was thinking about doing! If you have an idea for creating the base mesh, I'd love to see it. I have a number of routines for doing the minimal surface calculations. – mh2000 Oct 19 at 17:03
• I updated the answer to show how to create an equivalenced cap mesh in 2D. – Tim Laska 2 days ago
• @mh2000 Although I used RegionProduct to create a "structured mesh", it is fortuitous that it was the first 36 elements since I don't know the inner workings of the internal functions. You probably need an actual test to be sure. Also, I added a minimal surface example for the top surface. – Tim Laska 2 days ago

Find the lines on the mesh polygons that are not shared with other polygons:

curv1 = {{3, 0, 0}, {1, 1, 0}, {0, 2, 0}, {-2, 0, 0}, {0, -2, 0}, {3, 0, 0}};
curv2 = {{2, 0, 5}, {1, 1, 2}, {0, 2, 2}, {-1, 0, 5}, {0, -2, 5}, {3, 0, 5}};
sur1 = BSplineSurface[{curv1, curv2}, SplineClosed -> {False, True}, SplineDegree -> 3];
mesh = DiscretizeGraphics[sur1];
cells = MeshCells[mesh, 2][[All, 1]];
lines = Join @@ ((Sort /@ Subsets[#, {2}]) & /@ cells);
edgecells = Select[Tally[lines], Last[#] == 1 &][[All, 1]];
coords = MeshCoordinates[mesh];
edgelines = Line[{coords[[First[#]]], coords[[Last[#]]]}] & /@ edgecells;
Graphics3D[{{EdgeForm[None], Opacity[.5], mesh}, Red, Thick,
edgelines}, Boxed -> False] Note there is a seam in the mesh because the mesh is not connected there.

• Thank you! Much smarter than the road I was going down (comparing vertices to the splines with a tolerance). Is there an easy way to close that seam? Also, how to separate the upper and lower edges? – mh2000 Oct 16 at 21:49

Update: Two alternative, more direct, ways to get the boundary edges:

1. Use the property "EdgeFaceConnectivityRules" and select the edges connected to 0:

be1 = Keys @ Select[#[] == 0 &] @ Association[ mr["EdgeFaceConnectivityRules"]]

be1 == boundaryedges

True


2. Use the property "ConnectivityMatrix"[1, 2]" (which gives a SparseArray where entry $$ij$$ is 1 iff 1-dimensional element $$i$$ is connected to 2-dimensional element $$j$$) and select the rows that contain a single element:

be2 = Flatten @

be2 == boundaryedges

True


You can process mr["FaceEdgeConnectivityRules"] to identify edges connected to a single face:

mr = DiscretizeGraphics[sur1];

boundaryedges = Keys @ Select[EqualTo @ 1] @
Counts @ Flatten @ Values @ mr["FaceEdgeConnectivityRules"];

HighlightMesh[mr, {1, boundaryedges}, PlotTheme -> "Lines"] Then we can use mr["EdgeVertexConnectivityRules"] to identify the vertices incident to boundaryedges:

boundaryvertices = Union @@ (boundaryedges /. mr["EdgeVertexConnectivityRules"]);

HighlightMesh[mr,
Style[{0, boundaryvertices}, PointSize[Medium], Red],
PlotTheme -> "Lines"] • Nice. Where did you find this property? I can't find it in the documentation. – flinty Oct 16 at 22:55
• @flinty, mr["Properties"] list a number of intriguing properties. None are documented; so trial/error is the only way to explore. (Over the last few versions/updates more and more of these properties are working.) – kglr Oct 16 at 23:02