Is there a way to generate 2D stretched mesh using CGAL?

Question

I currently use CGAL to generate 2D Delaunay triangulation.One of the mesh control parameter is the maximum length of the triangle edge. The examples suggests that this parameter is a constant. I would like to know how this parameter be made function of some thing else, for example spatial location.

Answer 1

I think Delaunay meshing with variable density is not directly supported by CGAL although you could mesh your regions independently. Alternatively you may have a look at: http://www.geom.at/advanced-mesh-generation/ where I have implemented that as a callback function.

Answer 2

It doesn't look like CGAL provides an example of this but they machinery is all there. The details get a little complicated since the objects that control if triangles need to be refined also have to understand the priority under which triangles get refined.

To do this, I copied Delaunay_mesh_size_criteria_2 to create a new class (Delaunay_mesh_user_criteria_2) that has a spatially varying sizing field. Buried in the class is a function (user_sizing_field) that can be implemented with a varying size field based on location. The code below compares the size of the longest edge of the triangle to the minimum of the sizing field at the three vertices, but you could use a size at the barycenter or circumcenter or even send the entire triangle to the sizing function if you have a good way to compute the smallest allowable size on the triangle altogether.

This is a starting point, although a better solution would,

• refactor some things to avoid so much duplication with with existing Delaunay_mesh_size_criteria,
• allow the user to pass in the sizing function as an argument to the criteria object, and
• be shipped with CGAL.

template <class CDT>
class Delaunay_mesh_user_criteria_2 : 
    public virtual Delaunay_mesh_criteria_2<CDT>
{
protected:
  typedef typename CDT::Geom_traits Geom_traits;
  double sizebound;

public:
  typedef Delaunay_mesh_criteria_2<CDT> Base;

  Delaunay_mesh_user_criteria_2(const double aspect_bound = 0.125,
                                const Geom_traits& traits = Geom_traits())
    : Base(aspect_bound, traits){}

  // first: squared_minimum_sine
  // second: size
  struct Quality : public std::pair<double, double>
  {
    typedef std::pair<double, double> Base;

    Quality() : Base() {};
    Quality(double _sine, double _size) : Base(_sine, _size) {}

    const double& size() const { return second; }
    const double& sine() const { return first; }

    // q1<q2 means q1 is prioritised over q2
    // ( q1 == *this, q2 == q )
    bool operator<(const Quality& q) const
    {
      if( size() > 1 )
        if( q.size() > 1 )
          return ( size() > q.size() );
        else
          return true; // *this is big but not q
      else
        if( q.size() >  1 )
          return false; // q is big but not *this
      return( sine() < q.sine() );
    }

    std::ostream& operator<<(std::ostream& out) const
    {
      return out << "(size=" << size()
         << ", sine=" << sine() << ")";
    }
  };

  class Is_bad: public Base::Is_bad
  {
  public:
    typedef typename Base::Is_bad::Point_2 Point_2;

    Is_bad(const double aspect_bound,
           const Geom_traits& traits)
      : Base::Is_bad(aspect_bound, traits) {}

    Mesh_2::Face_badness operator()(const Quality q) const
    {
      if( q.size() > 1 )
        return Mesh_2::IMPERATIVELY_BAD;
      if( q.sine() < this->B )
        return Mesh_2::BAD;
      else
    return Mesh_2::NOT_BAD;
    }

    double user_sizing_function(const Point_2 p) const
    {
      // IMPLEMENT YOUR CUSTOM SIZING FUNCTION HERE.
      // BUT MAKE SURE THIS RETURNS SOMETHING LARGER
      // THAN ZERO TO ALLOW THE ALGORITHM TO TERMINATE
      return std::abs(p.x()) + .025;
    }

    Mesh_2::Face_badness operator()(const typename CDT::Face_handle& fh,
                    Quality& q) const
    {
      typedef typename CDT::Geom_traits Geom_traits;
      typedef typename Geom_traits::Compute_area_2 Compute_area_2;
      typedef typename Geom_traits::Compute_squared_distance_2 Compute_squared_distance_2;

      Geom_traits traits; /** @warning traits with data!! */

      Compute_squared_distance_2 squared_distance = 
        traits.compute_squared_distance_2_object();

      const Point_2& pa = fh->vertex(0)->point();
      const Point_2& pb = fh->vertex(1)->point();
      const Point_2& pc = fh->vertex(2)->point();

      double size_bound = std::min(std::min(user_sizing_function(pa),
                                            user_sizing_function(pb)),
                                   user_sizing_function(pc));
      double
        a = CGAL::to_double(squared_distance(pb, pc)),
        b = CGAL::to_double(squared_distance(pc, pa)),
        c = CGAL::to_double(squared_distance(pa, pb));

      double max_sq_length; // squared max edge length
      double second_max_sq_length;

      if(a<b)
      {
        if(b<c) {
          max_sq_length = c;
          second_max_sq_length = b;
        }
        else { // c<=b
          max_sq_length = b;
          second_max_sq_length = ( a < c ? c : a );
        }
      }
      else // b<=a
      {
        if(a<c) {
          max_sq_length = c;
          second_max_sq_length = a;
        }
        else { // c<=a
          max_sq_length = a;
          second_max_sq_length = ( b < c ? c : b );
        }
      }

      q.second = 0;

      q.second = max_sq_length / (size_bound*size_bound);
      // normalized by size bound to deal
      // with size field
      if( q.size() > 1 )
      {
        q.first = 1; // (do not compute sine)
        return Mesh_2::IMPERATIVELY_BAD;
      }

      Compute_area_2 area_2 = traits.compute_area_2_object();

      double area = 2*CGAL::to_double(area_2(pa, pb, pc));

      q.first = (area * area) / (max_sq_length * second_max_sq_length); // (sine)

      if( q.sine() < this->B )
        return Mesh_2::BAD;
      else
        return Mesh_2::NOT_BAD;
    }
  };

  Is_bad is_bad_object() const
  { return Is_bad(this->bound(), this->traits /* from the bad class */); }
};

Answer 3

I am also interested for variable mesh criteria on the domaine with CGAL. I have found an alternative many years ago : https://www.cs.cmu.edu/~quake/triangle.html

But i am still interested to do the same things with CGAL ... I don't know if it is possible ...