Cgal surface to deal.II triangulation

[fdrmrc]

This allows to go from a CGAL::Surface_mesh or a CGAL::Polyhedron_3 to a Triangulation<2,3>. Currently this is supporting only surfaces made by quads.

Depends on #13661. Relevant changes to this PR can be seen at: luca-heltai/dealii@c965b7f...0dce7af

[fdrmrc]

The relative .cc test is done by reading two CGAL surfaces from eight.off and creating the corresponding Triangulation<2,3>.

[luca-heltai]

I removed a check for the orientation in the simplex case, since that was hardcoded for quads. The check fails also on some quad meshes.

With this PR we can now finally compute BEM scattering problems around gummy bears. :)

[luca-heltai]

/rebuild

[peterrum]

Why is this not an issue in other cases?

[luca-heltai]

It is an issue. But this loop won't work in that case (and it doesn't work also in the case of quadrilaterals for some perfectly valid grids, so I'm not sure what the loops actually does).

[peterrum]

Who is "I"?

[luca-heltai]

I have no idea. In this whole comment, I only added "and made of quadrilaterals", to clarify that the check is skipped in the case of triangular meshes.

[peterrum]

Maybe you can skip this for triangulation generation but are the results obtained on such meshes correct?

[luca-heltai]

If we check somewhere else that the orientation is correct in the case of triangles, then we should not have problems here. CGAL generates consistently oriented meshes. We should definitely check that they are consistently oriented for tringle meshes. How is this done in two dimensions + simplices now? If what we do when we create the grid is enough in 2d+simplices, then it should be ok also in this case. I think the issues with quads is that we have a highly non-standard orientation in our triangulation class (i.e., we do not base our triangulation on the half-edge datastructure, but on lexycographical ordering, and consistent orientation of edges across neighbors). This gives issues when extracting boundary triangulation (for example), since some faces will have the right orientation, and some other won't. I believe this is the reason for this check to be here in the first place.

This is not the case for triangles. Triangles follow the same orientation of most other libraries (anti-clockwise), and so does our tetra-hedron classes, so each face is oriented consistently w.r.t. to each cell, and you are guaranteed that even if you extract the bundary of a tet-only grid, you get a consistently oriented grid. I believe the check for triangles (and tetra) should be that, on each edge cell->face(f)->vertex_index(0) == cell->neighbor(f)->face(cell->neighbor(f)->neighbor_of_neighbor(f))->vertex_index(1) is true, but I have to double check this.

[peterrum]

We should definitely check that they are consistently oriented for tringle meshes. How is this done in two dimensions + simplices now?

I assumed that the meshes coming from a CAD tool are correctly oriented. But unfortunately, this is not always given as the issue #13346 showed. It was fixed in #13579. For the normal case dim==spacedim, we have the checks:

dealii/source/grid/tria.cc

Lines 2309 to 2323 in 5cdf7bd

	if (dim == spacedim)
	for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
	{
	// If we should check for distorted cells, then we permit them
	// to exist. If a cell has negative measure, then it must be
	// distorted (the converse is not necessarily true); hence
	// throw an exception if no such cells should exist.
	if (tria.check_for_distorted_cells)
	{
	const double cell_measure = GridTools::cell_measure<spacedim>(
	vertices,
	ArrayView<const unsigned int>(cells[cell_no].vertices));
	AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
	}
	}

and

dealii/source/fe/mapping_fe.cc

Lines 1233 to 1240 in 5cdf7bd

	// TODO: this allows for anisotropies of up to 1e6 in 3D and
	// 1e12 in 2D. might want to find a finer
	// (dimension-independent) criterion
	Assert(det >
	1e-12 * Utilities::fixed_power<dim>(
	cell->diameter() / std::sqrt(double(dim))),
	(typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
	cell->center(), det, point)));

that will definitely be triggered if the mesh is not correctly orientated.

I believe the check for triangles (and tetra) should be that, on each edge cell->face(f)->vertex_index(0) == cell->neighbor(f)->face(cell->neighbor(f)->neighbor_of_neighbor(f))->vertex_index(1) is true, but I have to double check this.

Do you want to open an issue? I guess one can come up with a minimal test with a hand full of cells?

[luca-heltai]

It is unclear to me what we want to do. For the codimension zero case, it is clear that if a mesh has negative determinant, it should not be allowed. For co-dimension one grids, however, the J in JxW is defined as sqrt(det(DF^T.DF)), and it is positive regardless of the orientation of the normal (i.e, the previous expression is the same as computing the norm of the cross product of the tangent vectors, or the columns of the Jacobian). The orientation of the mesh in the codimension one case does not influence the sign of the determinant but just the orientation of the normal to the cell. Do we want to enforce a consistent orientation, or should we allow discontinuous normals? A moebious strip is not orientable, but it is a perfectly valid co-dimension one mesh... It won't make sense to ask for a continuous normal in many cases. I think we should remove the check, and let the user be responsible for the orientation.

[luca-heltai]

I know my previous comments partly contradicts what I said earlier, but I did not have the moebious strip example in mind when I wrota that

If we check somewhere else that the orientation is correct in the case of triangles, then we should not have problems here. CGAL generates consistently oriented meshes.

[peterrum]

I am not sure what the influence of discontinuous normals are so I cannot really tell when/where it is the best to deal with them. Maybe you could post a simple example?

[peterrum]

Fine from my side to merge! @luca-heltai You are blocking the PR (changes requested) ;)

[luca-heltai]

Thanks @fdrmrc ! I really like where we are going... :)