Suggest improvement for the API of Mesh_criteria_3

[nschloe]

I would like to create a mesh with variable cell size, and only the cell size, not edge size, angles etc., should influence the mesh. (See https://github.com/nschloe/pygalmesh/issues/114.) This is using pygalmesh, but the underlying issue in in CGAL.

Setting the edge_size to 0.0 does not work:

import pygalmesh

class Field(pygalmesh.SizingFieldBase):
    def eval(self, x):
        return (1 - x[0]) * 0.025 + 0.025

mesh = pygalmesh.generate_mesh(
    pygalmesh.Cuboid([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]),
    cell_size=Field(),
    edge_size=0.0,
    # edge_size=100000000.0
)

mesh.write("cgal_cuboid.vtk")

File: /usr/include/CGAL/Mesh_3/Protect_edges_sizing_field.h
Line: 448
Explanation: Error: the sizing field is null at corner (0 0 0)

Fair enough. Let's set it to some really large value (which is the default, too). Result:

This isn't good either, as it seems that the cell size is not respected at the edges.

In pygalmesh, the mesh criteria are set as

    auto criteria = Mesh_criteria(
        CGAL::parameters::edge_size=edge_size,
        CGAL::parameters::facet_angle=facet_angle,
        CGAL::parameters::facet_size=facet_size,
        CGAL::parameters::facet_distance=facet_distance,
        CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio,
        CGAL::parameters::cell_size=size
        );

[lrineau]

edge_size cannot be null: In the documentation of the constructor CGAL::Mesh_criteria_3< Tr >::Mesh_criteria_3, it is said:

edge_size: a scalar field (resp. a constant) providing a space varying (resp. a uniform) upper bound for the lengths of curve edges. This parameter has to be set to a positive value when 1-dimensional features protection is used.

In your Python example, you can use edge_size=Field(), and it will have the effect you want, I think.

[nschloe]

I'd have to update pygalmesh for it I think. So edge_size can also eat a std::function à la

const auto size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
    return 0.5;  // some evaluation
};

?

[lrineau]

That is right: it can be a sizing field.

Note that the second argument is the dimension of the element:

const auto size = [&](K::Point_3 p, const int dimension, const Mesh_domain::Index&) {
    if(dimension <= 1) { // for edges and corners
        return 0.8; 
    } else if(dimension == 2) { // for facets
        return 1.0;
    } else { // for cells
        return 0.5;  // some evaluation
    }
};

The third argument is actually the index of the input feature the element is on. Its meaning depends on your mesh domain definition, and I will not talk more about that for the moment.

[nschloe]

Alright, thank you. One minor gripe I have is that cell_size use circumradii as a bound, but edge size uses the edge length.

cell_size: a scalar field (resp. a constant) describing a space varying (resp. a uniform) upper-bound for the circumradii of the mesh tetrahedra.

Using the same Field() for both will result in some distortion. Well, let's assume it won't be too awful.

[nschloe]

Question: How does CGAL implement a "float or function" interface? I'm getting cursed by combinatorics. (Link would already be helpful.)

[lrineau]

Alright, thank you. One minor gripe I have is that cell_size use circumradii as a bound, but edge size uses the edge length.

You are right. And in fact the edge size being a bound on the edge length is actually a lie, or at least imprecise. The precise implementation is complicated to explain without the reading the scientific articles.

Using the same Field() for both will result in some distortion. Well, let's assume it won't be too awful.

There will always be a kind of distortion. The way the sharp edges are protected pushes away new vertices from the sharp features. The triangles incident to the sharp features will always have an aspect ratio slightly degraded from the other ones. For example, on a cube, if I use the same constant for edge_size and facet_size, then: Note that the triangles incident to the sharp edges have one edge on the sharp edge that is slightly smaller than the rest of the mesh, but their two other edges have edge lengths similar to the triangles that are not-incident to sharp edges.

On the contrary, I can use a 50% bigger edge_size: Then the edges on the sharp edges have roughly the same size as the other ones, but the two other edges of triangles incident to sharp edges are bigger.

In my opinion, the mesh is better-looking when the edge_size and facet_size have the same sizing field, or sizing constant.

[nschloe]

In my opinion, the mesh is better-looking when the edge_size and facet_size have the same sizing field, or sizing constant.

Very good, thanks for looking into it.

[lrineau]

Question: How does CGAL implement a "float or function" interface? I'm getting cursed by combinatorics. (Link would already be helpful.)

There is a virtual base class template: https://github.com/CGAL/cgal/blob/afc0abb82f341e5dadaae545dc3ac51d30a8f5dc/Mesh_3/include/CGAL/Mesh_edge_criteria_3.h#L33-L48

And then two constructors for Mesh_edge_criteria: https://github.com/CGAL/cgal/blob/afc0abb82f341e5dadaae545dc3ac51d30a8f5dc/Mesh_3/include/CGAL/Mesh_edge_criteria_3.h#L99-L121

The combinatorics is simple, because Mesh_criteria_3 actually contains three members: https://github.com/CGAL/cgal/blob/afc0abb82f341e5dadaae545dc3ac51d30a8f5dc/Mesh_3/include/CGAL/Mesh_criteria_3.h#L142-L147 and each of the three classes implements the two constructors like Mesh_edge_critera_3.

[nschloe]

Cool. I'll figure out the rest.

[nschloe]

Unfortunately, I'll have to reopen this. To create a Mesh_criteria object,

    const auto criteria = Mesh_criteria(
        CGAL::parameters::edge_size=edge_size,
        CGAL::parameters::facet_angle=facet_angle,
        CGAL::parameters::facet_size=facet_size,
        CGAL::parameters::facet_distance=facet_distance,
        CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio,
        CGAL::parameters::cell_size=cell_size
        );

from either input functions or input floats, I create functions no matter what. For cell_size:

    std::function<double(K::Point_3, const int, const Mesh_domain::Index&)> cell_size;
    if (cell_size_field) {
      cell_size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
        return cell_size_field->eval({p.x(), p.y(), p.z()});
      };
    } else {
      cell_size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
        return cell_size_value;
      };
    }

(This could be shortened a bit with C++17 and the ternary operator.) Passing this as mesh criterion works well for cell_size and edge_size. When doing it for facet_size and facet_distance, the mesh generation process takes forever. I have no idea what the difference between

const auto criteria = Mesh_criteria(
    CGAL::parameters::facet_size=0.0
    );

and

auto facet_size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
  return 0.0;
};
const auto criteria = Mesh_criteria(
    CGAL::parameters::facet_size=facet_size
    );

could be. Any hints? I would really love to see something like

auto criteria = Mesh_criteria()
criteria.add_facet_size_criterion(facet_size);

which would alleviate me of such problems.

[janetournois]

IIRC, the implementation of the facet_size criterion is such that :

• facet_size = 0 is equivalent to no size criterion at all
• facet_size = a_function_that_always_returns_0 is actually taken into account as a "real" sizing field, and tries to create 0-sized faces... forever. Can we test that the input function always returns 0?

[nschloe]

Hm yes, that'd make sense. It's different for edge and cell size, right? I think I'll just scrap field specification for facets in pygalmesh then. To reiterate, it'd be great if one could set the criteria individually.

[lrineau]

I agree we should improve the API of Mesh_criteria_3. @nschloe's suggestions are interesting.

[nschloe]

I just noticed that cell_size behaves the same way: If set the 0, it's treated as nonexistent, if set to a function that always returns 0, it's taken for face-value, resulting in weird behavior (in my case, a segfault). Before the new API lands, do you have an idea of how to implement float/function selection? The following approach does not work for the above reason:

  // cell_size_field is a ptr, take it if not nullptr
  // if nullptr, fall back to cell_size_value

  std::function<double(K::Point_3, const int, const Mesh_domain::Index&)> cell_size;
  if (cell_size_field) {
    cell_size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
      return cell_size_field->eval({p.x(), p.y(), p.z()});
    };
  } else {
    cell_size = [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
      return cell_size_value;
    };
  }

  // do the same for edge_size etc.

  const auto criteria = Mesh_criteria(
      CGAL::parameters::edge_size=edge_size,
      CGAL::parameters::cell_size=cell_size
      );

[nschloe]

Unfortunately, this isn't legal C++:

    const auto criteria = Mesh_criteria(
        CGAL::parameters::cell_size=cell_size_field ?
          [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
              return cell_size_field->eval({p.x(), p.y(), p.z()});
          } : cell_size_value
        );

error: incompatible operand types ('(lambda at src/generate.cpp:168:9)' and 'const double')
      CGAL::parameters::cell_size=cell_size_field ?
                                                  ^

[nschloe]

I had some hopes for the following, but doesn't work either.

typedef Mesh_criteria::Cell_criteria Cell_criteria;

        CGAL::parameters::cell_size=cell_size_field ?
          Cell_criteria(0.0, [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
              return cell_size_field->eval({p.x(), p.y(), p.z()});
          }) : Cell_criteria(0.0, cell_size_value)

  Mesh_cell_criteria_3(const FT& radius_edge_bound,
  ^
/usr/include/CGAL/Mesh_cell_criteria_3.h:34:7: note: candidate constructor (the implicit copy constructor) not viable: requires 1 argument, but 2 were provided
class Mesh_cell_criteria_3
      ^

Some template-magic going on here; I'm out of my depth.

[lrineau]

@nschloe

Use the ternary operator to construct a Cell_criteria, that way:

  typedef Mesh_criteria::Cell_criteria Cell_criteria;
  Cell_criteria cell_criteria = cell_size_field ?
     Mesh_criteria(cell_radius_edge_ratio,
                  [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
                    return cell_size_field->eval({p.x(), p.y(), p.z()});
                  }) :
     Mesh_criteria(cell_radius_edge_ratio, cell_size_value);

Do the same for the edge_criteria and facet_criteria, and then:

  Mesh_criteria criteria = Mesh_criteria(edge_criteria, facet_criteria, cell_criteria);

[nschloe]

This should be

  Cell_criteria cell_criteria = cell_size_field ?
     Cell_criteria(cell_radius_edge_ratio,
                  [&](K::Point_3 p, const int, const Mesh_domain::Index&) {
                    return cell_size_field->eval({p.x(), p.y(), p.z()});
                  }) :
     Cell_criteria(cell_radius_edge_ratio, cell_size_value);

but yeah, it works. For facets it's getting messy because there are two arguments which could be fields (facet_size and facet_distance), but a nested ternary solves it.

[lrineau]

Good. I am glad that you managed to setup the code as you wished.

But I agree there should be another API that ease that sort of thing. That is why we will keep that issue open.

[nschloe]

Has there been any work on cell size criteria yet?