7
$\begingroup$

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]

enter image description here

$\endgroup$
8
$\begingroup$

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["NDSolve`FEM`"]
(* 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[0] /@ 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]]

Strucured Mesh

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[
    Length /@ mr["ConnectivityMatrix"[1, 2]]["AdjacencyLists"], 1];
HighlightMesh[mr, Style[{1, boundaryedgeindices}, Thick, Red]]

Free Edges

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[0] /@ 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.

2D Mesh

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[[2]]];
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]}]]]

Top Surface Mesh

You can see that the free surface nodes align well with the base mesh.

| improve this answer | |
$\endgroup$
  • $\begingroup$ That looks great, Tim! Do you know how to close the ends of this irregular tube with a mesh? $\endgroup$ – mh2000 Oct 17 at 17:45
  • 1
    $\begingroup$ @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. $\endgroup$ – Tim Laska Oct 18 at 19:05
  • $\begingroup$ 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. $\endgroup$ – mh2000 Oct 19 at 17:03
  • 1
    $\begingroup$ I updated the answer to show how to create an equivalenced cap mesh in 2D. $\endgroup$ – Tim Laska 2 days ago
  • 1
    $\begingroup$ @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. $\endgroup$ – Tim Laska 2 days ago
7
$\begingroup$

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]

mesh hole edges

Note there is a seam in the mesh because the mesh is not connected there.

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – mh2000 Oct 16 at 21:49
7
$\begingroup$

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[#[[1]] == 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 @ 
  Position[Length /@ mr["ConnectivityMatrix"[1, 2]]["AdjacencyLists"],  1];

be2 == boundaryedges
True

Original answer:

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

enter image description here

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

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Nice. Where did you find this property? I can't find it in the documentation. $\endgroup$ – flinty Oct 16 at 22:55
  • 4
    $\begingroup$ @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.) $\endgroup$ – kglr Oct 16 at 23:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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