CSGCGALMeshGenerator2D: Keep location of vertices given in polygon in generated mesh

[Tom Andreas Nærland]

In the old mesh generator in Dolfin it was possible to retain the location of vertices given in a CSG polygon in the generated mesh. This does not seem to be possible now due to how the PSLG class in CSGCGALDomain2D uses the CGAL::snap_rounding_2 which sweeps "pixels" to snap together close vertices. I propose to make the snap_rounding_2 routine optional by skipping the routine if a negative parameter value is given (as is the case with some of the other features).

[Benjamin Dam Kehlet]

Yes, it may be a good idea to optionally skip the snap rounding. I'll take a look.

The use of CGAL::snap_rounding_2 is suboptimal for mshr for at least a couple of reasons: * It is quite slow (I had to patch CGAL to work around a bug that they seem unwilling to fix, because the fix slows down the algorithm) * It snaps vertices to the pixel centers even if there is only one vertex within the pixel. * It is not designed for this use, so the output format (a list of polylines) is not optimal.

Unfortunately it is a little more complicated than just skipping the entire routine since this is what computes the pslg. This may introduce completely new vertices if there are subdomains defined. Skipping the snap rounding routine will then crash the mesh generator. I think this can happen only if there are subdomains defined (otherwise the boolean operations package will have introduced all the vertices).

A better solution in the long term would be a proper implementation of the PSLG class which only moves vertices if they need to be merged (and maybe issues a warning).

[Tom Andreas Nærland]

Hi! Many thanks for mshr, by the way!

And yeah, snap_rounding_2 seems very slow. In the old mesh generator the PSLG was just naively given by the order of the input vertices. This was handy since it was cheap and sufficient for simply connected domains without subdomains. I suggest that such a naive PSLG is given if the pixel size is negative (since collapsing some vertices kind of indicates that one is willing to lose the information about the input vertices) or if one supplies some kind of parameter indicating that the PSLG is to be derived from the order of the given vertices.

Below is how I've edited the PSLG to accomodate my needs. It is just skipping the whole snap_rounding_2. This gives the same behaviour as the old mesh generator. (Why I need this behaviour can be seen in the video: http://vimeo.com/100802906 where an initial generated mesh is deformed by a changing-in-time boundary mesh (via Mesh::move(BoundaryMesh)) matching at the boundary vertices in the generated mesh).

PSLG::PSLG(const std::vector<std::pair<std::size_t, CSGCGALDomain2D> >& domains, 
           double pixel_size, 
           double edge_truncate_tolerance)
{
  // Collect all segments from all domains to send to snap rounding
  std::vector<Segment_2> segments;
  const Quotient truncate_tolerance(edge_truncate_tolerance);
  const Quotient truncate_tolerance_squared(truncate_tolerance*truncate_tolerance);

  for (auto it = domains.begin(); it != domains.end(); it++)
  {
    const Polygon_set_2& p = it->second.impl->polygon_set;

    std::list<Polygon_with_holes_2> polygon_list;
    p.polygons_with_holes(std::back_inserter(polygon_list));

    for (std::list<Polygon_with_holes_2>::const_iterator pit = polygon_list.begin();
         pit != polygon_list.end(); ++pit)
    {
      add_simple_polygon(segments, pit->outer_boundary(), truncate_tolerance_squared);

      // Add holes
      Hole_const_iterator hit;
      for (hit = pit->holes_begin(); hit != pit->holes_end(); ++hit)
      {
        add_simple_polygon(segments, *hit, truncate_tolerance_squared);
      }
    }
  }

  std::map<Point_2, std::size_t> point_to_index;

  for (std::vector<Segment_2>::const_iterator iter1 = segments.begin();
       iter1 != segments.end(); ++iter1)
  {
    std::size_t prev = get_vertex_index(point_to_index,
                                               vertices,
                                               (*iter1).source());
    std::size_t current = get_vertex_index(point_to_index,
                                             vertices,
                                             (*iter1).target());

      edges.push_back(std::make_pair(prev, current));
      prev = current;
  }
}

[Benjamin Dam Kehlet]

This should be fixed now. I disabled the snap rounding when there are no subdomains. I also decreased the pixel size of the snap rounding should it be within the tolerance of dolfin::near(). This means that no vertices should be moved more than 1e-16. It is still not optimal (one would expect that vertices are kept exact at their original location), but acceptable for now.

I'm closing this one now, since there is already a registered issue about reimplementing the PSLG code (https://bitbucket.org/benjamik/mshr/issue/6/implement-pslg)

Nice work with the video! Cool to see the ALE functionality in Dolfin at work!