Generating a boundary mesh from a list of points

Question

I have a list of points from which I am trying to generate a mesh:

 {{{261.324, 937.805}, {1139.98, 937.843}}, {{261.324, 
   937.805}, {261.014, 57.1923}}, {{1140.95, 505.675}, {705.906, 
   57.3194}}, {{1578.03, 507.712}, {1139.98, 937.843}}, {{261.014, 
   57.1923}, {705.906, 57.3194}}, {{1140.95, 505.675}, {1578.03, 
   507.712}}}

fldd is a list of coordinates (without duplicates) generated by flattening the given array and using the DeleteDuplicates command. That definitely works. However, the BoundaryElements part of the command is giving me issues.

I use the following code to generate the connectivity array which tells Mathematica in what order to connect the points:

cn = {}
For[c1 = 1, c1 <= len, c1++,
 For[c2 = 1, c2 <= 2, c2++,
  holder = Position[lines, lines[[c1, c2]]];
  sholder = Sort[holder];
  Print[sholder];
  If[FreeQ[cn, {sholder[[1, 1]], sholder[[2, 1]]}],
   Print["In If"];
   AppendTo[cn, {sholder[[1, 1]], sholder[[2, 1]]}]
   ]
  ]
 ]

When I run the code below to actually generate the boundary mesh, I get an error:

bmesh = ToBoundaryMesh["Coordinates" -> fldd[[1 ;; Length[fldd]]], 
  "BoundaryElements" -> {LineElement[cn[[1 ;; Length[cn]]]]}]

ToBoundaryMesh::fememib: The input has or generated an intersecting boundary and cannot be processed. >>

I can't really see what's going wrong. I know that the loop above generates a list of numbers of the form:

{{startpoint1, endpoint1}, ..., {startpointN, endpointN}} 

This looks identical to the examples I've seen. I also know that the points I want connected are connected, and each point is listed as a vertex for 2 lines, which is what I want geometrically. Any help would be greatly appreciated.

Full code is:

lines = {
   {{261.324, 937.805}, {1139.98, 937.843}},
   {{261.324, 937.805}, {261.014, 57.1923}}, 
   {{1140.95, 505.675}, {705.906, 57.3194}}, 
   {{1578.03, 507.712}, {1139.98, 937.843}}, 
   {{261.014, 57.1923}, {705.906, 57.3194}},
   {{1140.95, 505.675}, {1578.03, 507.712}}
  };


len = Length[lines];
fl = Flatten[lines, 1];
fldd = DeleteDuplicates[fl];

(*To find connectivity (Uggh)*)
cn = {};
For[c1 = 1, c1 <= len, c1++,
 For[c2 = 1, c2 <= 2, c2++,
  holder = Position[lines, lines[[c1, c2]]];
  sholder = Sort[holder];
  Print[sholder];
  If[FreeQ[cn, {sholder[[1, 1]], sholder[[2, 1]]}],
   Print["In If"];
   AppendTo[cn, {sholder[[1, 1]], sholder[[2, 1]]}]
  ]
 ]
]

bmesh = ToBoundaryMesh[
    "Coordinates" -> fldd[[1 ;; Length[fldd]]], 
    "BoundaryElements" -> {LineElement[cn[[1 ;; Length[cn]]]]}
  ]
bmesh["Wireframe"]
emesh = ToElementMesh[bmesh]
emesh["Wireframe"]

Answer 1

Your problem is the fact that the lines in your list of lines are not contiguous, or in other words the points in your fldd list generated from the lines variable are not each other's nearest neighbors.

Another issue may be the nested For loops. If you find yourself using procedural constructs such as looping constructs, you might want to rethink your approach, since there will typically be a more efficient and often simpler one using a functional approach instead. In other words, you want to get into the mindset of telling Mathematica what you want done, and not how to do it! It is brain-twisting at first, but it pays off in the long run.

The problem: First of all, let's see why ToBoundaryMesh was complaining in your original implementation. I will use the definition of lines that you provided in the following. This is along the lines of what you were trying to do:

Needs["NDSolve`FEM`"]

points = DeleteDuplicates[Flatten[lines, 1]];

ToBoundaryMesh[
  "Coordinates" -> points,
  "BoundaryElements" -> {LineElement@Partition[Range[Length[points]], 2, 1]}
]["Wireframe"]

As you can see, the mesh is generated from the points in the order in which you provide them. When you attempted to close that mesh region by connecting the top leftmost point with the rightmost point, the function found that this last LineElement would intersect one of the existing lines, and threw an error.

The solution: We need to order your points in such a way that they can be traversed in a continuous path that goes through each of them once. Sounds familiar? It may, because this is more or less the classical traveling salesman problem. Fortunately, Mathematica already has a built in function to accomplish that: FindShortestTour. This function takes a list of points, and returns the length of the tour (we don't care about that), and a list of indices that indicates in what order one should take the original points to generate this closed path: this last piece of information is exactly what we need.

boundarymesh =
  ToBoundaryMesh[
   "Coordinates" -> points,
   "BoundaryElements" -> {LineElement@ Partition[Last@FindShortestTour[points], 2, 1]}
   ];

boundarymesh["Wireframe"]

Now that we have the boundary mesh, we can create an element mesh from it:

elementmesh = ToElementMesh[boundarymesh];

elementmesh["Wireframe"]

The functional approach: Note how much easier it is to produce your desired connectivity array using a more idiomatic approach than For loops. Let's dissect what is going on in the generation of the BoundaryElements list:

ordering = Last@FindShortestTour[points]
Partition[ordering, 2, 1]

In this case, Partition picks elements out of the ordering list two at a time (second argument), with an overlap between the generated lists of one element (the third argument). The output of the above is:

output of FindShortestTour: {1, 3, 5, 4, 6, 2, 1}
output of Partition:        {{1, 3}, {3, 5}, {5, 4}, {4, 6}, {6, 2}, {2,1}}