Cgal corefinement #111

fdrmrc · 2022-05-06T20:27:44Z

Just a draft to describe things a little.

The "Surface with edge"d part in the following image the is result of the corefinement test. It is just a CGAL::Surface_mesh.
Here, that precise surface is filled with tets (use medit visualize it). This has still to be added. Notice that this could be used to generate tetrahedral meshes, not to compute Quadratures over the result of BooleanOperations.

fdrmrc · 2022-05-07T15:50:15Z

Fixed output files
Added boolean operations between Surface_meshes
Added Quadrature formulas over volumetric region bounded by a Surface_mesh and simple test
Use coarse volumetric mesh to compute intersections(i.e. just compute the ConvexHull)

So far this works in 3D only. I'll try to make this work for <2,2> objects too.

luca-heltai

Can you make the PR on top of cgal-triangulation, instead of on top of master?

luca-heltai · 2022-05-07T18:05:35Z

include/deal.II/cgal/surface_mesh.h

+  enum class BooleanOperation
+  {
+    only_corefinement = 1 << 0,
+    union_op          = 1 << 2,
+    intersection_op   = 1 << 3,
+  };


Suggested change

enum class BooleanOperation

{

only_corefinement = 1 << 0,

union_op = 1 << 2,

intersection_op = 1 << 3,

};

enum class BooleanOperation

{

none = 1 << 0,

union = 1 << 1,

intersection = 1 << 2,

};

Shalle we support also difference?

Sure, thanks.

luca-heltai · 2022-05-07T18:06:30Z

include/deal.II/cgal/surface_mesh.h

 #  include <CGAL/Surface_mesh.h>
+#  include <CGAL/Triangulation_3.h>
+#  include <CGAL/convex_hull_3.h>
+

 DEAL_II_NAMESPACE_OPEN

 namespace CGALWrappers
 {


Add a bit of documentation, with some examples from the CGAL web site, or computed with the tests.

luca-heltai · 2022-05-07T18:10:26Z

include/deal.II/cgal/surface_mesh.h

+   * Extensive examples can be found in the CGAL documentation at
+   * https://doc.cgal.org/latest/Polygon_mesh_processing/index.html#coref_coref_subsec
+   *


I would put here some examples compute with this function, similar to the one you have in the comments.

luca-heltai · 2022-05-07T18:11:00Z

include/deal.II/cgal/surface_mesh.h

+   * @tparam CGALPointType
+   * @tparam dim
+   * @tparam spacedim


Remove the template paramaters comments.

luca-heltai · 2022-05-07T18:11:39Z

include/deal.II/cgal/surface_mesh.h

+   * @param[in] cell1
+   * @param[in] cell2
+   * @param[out] outsm


Suggested change

* @param[in] cell1

* @param[in] cell2

* @param[out] outsm

* @param[in] cell1 First input cell

* @param[in] cell2 Second input cell

* @param[out] outsm Output surface mesh

luca-heltai · 2022-05-07T18:17:04Z

include/deal.II/cgal/surface_mesh.h

+  /**
+   * Build a Quadrature<spacedim> formula to integrate over the volumetric
+   * region described by a CGAL::Surface_mesh, by filling the interior with
+   * simplices and collecting points and weights to initegrate over it.


What happens if the surface mesh is not closed?

I cannot think to a case where having a non-closed surface makes sense, after all we're intersecting deal.II cells in that function.

I know, but does the function check that surface mesh is closed or not? You should document what happens in that case, and if the function throws an exception, you should state this here.

luca-heltai · 2022-05-07T18:17:33Z

include/deal.II/cgal/surface_mesh.h

+   * @tparam spacedim
+   * @tparam CGALPointType


Remove parameters documentation

luca-heltai · 2022-05-07T18:17:41Z

include/deal.II/cgal/surface_mesh.h

+   * @param sm
+   * @param degree
+   * @return Quadrature<spacedim>


Comments input and output

luca-heltai · 2022-05-07T18:19:42Z

include/deal.II/cgal/surface_mesh.h

+   * @return Quadrature<spacedim>
+   */
+  template <int spacedim, typename CGALPointType, typename CGALTriangulation>
+  Quadrature<spacedim>


Suggested change

Quadrature<spacedim>

Quadrature<CGALTriangulation::Space_dimension()>

I'm not sure if this works, but this would allow us to remove the template parameter.

That is not constexpr. I used CGALPointType::Ambient_dimension::value and this is constexpr. It works. Thanks!

luca-heltai · 2022-05-07T18:20:42Z

include/deal.II/cgal/surface_mesh.h

+   */
+  template <int spacedim, typename CGALPointType, typename CGALTriangulation>
+  Quadrature<spacedim>
+  quadrature_inside_region(const CGAL::Surface_mesh<CGALPointType> &sm,


Do we like the name?

Honestly I don't, but I was not able to find a meaningful name. Maybe collect_quadratures_inside_surface ?

luca-heltai · 2022-05-07T18:31:42Z

Do you think we could make this more general by having it return a polygon soup? I would really love to have this as we did with parmoonolith, i.e., arbitrary cell VS arbitrary cell. Of course if it is not possible, then we simply need to specialize this for each possibility.

fdrmrc · 2022-05-08T11:36:20Z

Thanks for the comments! I addressed them and added a bit of documentation, except for the non-closed surface case.
For what concerns the Polygon soup, I'm looking into it right now.

luca-heltai · 2022-05-08T13:13:16Z

include/deal.II/cgal/surface_mesh.h

+    NONE         = 1 << 0,
+    UNION        = 1 << 2,
+    INTERSECTION = 1 << 3,
+    DIFFERENCE   = 1 << 4,


deal.II notation: lowercase. If none doesn't work, use corefine as you had at the beginning. Aso, each one should be commented (using the ///< doxygen tag):

union = 1<<2, ///< Compute the union of the input arguments intersection = 1<<2, ///< Compute the intersection...

luca-heltai · 2022-05-08T13:14:32Z

include/deal.II/cgal/surface_mesh.h

+   * Performs corefinement and boolean operations on two deal.II cells and put
+   * the result in a CGAL::Surface_mesh object. If the cells are disjoint and
+   * you try to perform a boolean operation which makes no sense, like their
+   * intersection, an exception is thrown.


Do we want to throw an exception, or simply return an empty mesh? In the end, an empty set is a valid result for an intersection from a mathematical standpoint.

luca-heltai · 2022-05-08T13:15:47Z

include/deal.II/cgal/surface_mesh.h

+   * @param[in] degree The degree of the quadrature formula
+   * @param[out] tr A CGAL triangulation storing the convex_hull of the volume
+   * bounded by the Surface_mesh
+   * @return Quadrature<spacedim> Collection of Quadratures. The sum of the weights equals to the volume of the polygonal region bounded by sm


Suggested change

* @return Quadrature<spacedim> Collection of Quadratures. The sum of the weights equals to the volume of the polygonal region bounded by sm

* @return Quadrature<spacedim> Collection of Quadratures. The sum of the weights

* equals to the volume of the polygonal region bounded by sm

luca-heltai · 2022-05-08T13:17:11Z

include/deal.II/cgal/surface_mesh.h

+    CGAL::Surface_mesh<CGALPointType> &                         outsm,
+    const BooleanOperation &boolean_operation = BooleanOperation::NONE);
+
+  /**


You no longer have spacedim at your disposal. You should either document what spacedim is, or use directly CGALPointType::Ambient_dimension::value in place of spacedim in the documentation.

luca-heltai

Very nice!

luca-heltai · 2022-05-08T13:18:53Z

include/deal.II/cgal/surface_mesh.h

+   * cube and a tetrahedron, converted to a deal.II triangulation
+   *
+   * @image html intersection_cube_with_tetrahedron.png
+   *


Perfect!

Can you also add the same picture computin the union and the difference?

luca-heltai · 2022-05-08T13:20:35Z

include/deal.II/cgal/surface_mesh.h

+   * the result in a CGAL::Surface_mesh object. If the cells are disjoint and
+   * you try to perform a boolean operation which makes no sense, like their
+   * intersection, an exception is thrown.
+   *


You should add a line in which you describe who you subtract from whom. I'm imagining you do cell1 \setminus cell2, but you should document this.

Also in this case, you could have an empty result (i.e., cell1 is entirely included in cell2).

luca-heltai · 2022-05-08T13:21:17Z

include/deal.II/cgal/surface_mesh.h

+   * Build a Quadrature<spacedim> formula to integrate over the volumetric
+   * region described by a CGAL::Surface_mesh, by filling the interior with
+   * simplices and collecting points and weights to integrate over it.
+   *


Add a comment that explains what happens if the operation results in a empty set: you return an empty quadrature.

luca-heltai · 2022-05-08T13:23:19Z

There is an issue with the image:

[100%] Built target documentation
/home/runner/work/dealii/dealii/include/deal.II/cgal/surface_mesh.h:51: warning: image file intersection_cube_with_tetrahedron.png is not found in IMAGE_PATH: assuming external image.
Error: Process completed with exit code 1.

luca-heltai · 2022-05-08T14:41:33Z

source/cgal/surface_mesh.cc

+  boolean_operation(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell1,
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell2,
+    CGAL::Surface_mesh<CGALPointType> &                         outsm,


Add mapping to the input argument list.

luca-heltai · 2022-05-08T14:43:54Z

source/cgal/surface_mesh.cc

+        case BooleanOperation::DIFFERENCE:
+          res = PMP::corefine_and_compute_difference(sm1, sm2, outsm);
+          break;
+        case BooleanOperation::NONE:


This is not going to produce anything visible on the outside. outsm is not going to be touched by this function, and sm1 and sm2 will be destroied at the end of the call. We should document what will happen in this case, i.e., we return sm1, i.e., the corefinement of the first cell, therefore here you need to do outsm = std::move(sm1);

luca-heltai · 2022-05-08T14:50:48Z

source/cgal/surface_mesh.cc

+
+  template <typename CGALPointType, typename CGALTriangulation>
+  Quadrature<CGALPointType::Ambient_dimension::value>
+  collect_quadratures_inside_surface(


maybe a name like build_quadrature_on_convex_hull?

Also, remove GCALTriangulation from the input argument. Since you are only using it inside the function, just build it inside.

We should also provide a different version of this function for non-convex polyhedral meshes. Maybe adding a boolean argument like integrate_on_convex_hull which defaults to false that uses the algorithm you showed me yesterday, i.e., create a domain, then a mesh, then last lines.

luca-heltai · 2022-05-08T14:54:31Z

source/cgal/surface_mesh.cc

+    const unsigned int                       degree,
+    CGALTriangulation &                      tr)
+  {
+    constexpr unsigned int spacedim = 3;


Suggested change

constexpr unsigned int spacedim = 3;

constexpr unsigned int spacedim = CGALPointType::Ambient_dimension::value;

luca-heltai · 2022-05-08T14:55:44Z

source/cgal/surface_mesh.cc

+    QGaussSimplex<spacedim>              quad(degree);
+    std::vector<dealii::Point<spacedim>> pts;
+    std::vector<double>                  wts;
+    for (const auto &f : tr.finite_cell_handles())


this is different if spacedim is 3 or 2

fdrmrc · 2022-06-06T14:09:59Z

See upstream

luca-heltai and others added 7 commits May 4, 2022 09:21

Insert points in CGAL triangulation.

1684628

Test Triangulation_2

5cfa394

Triangulation_3.

50ff618

Test sphere in 3d.

d01c884

Update changelog.

bda0de1

Minor fixes.

db13e84

Fix indent

b87f33b

fdrmrc force-pushed the cgal-corefinement branch 5 times, most recently from bc3b17d to 98970b3 Compare May 7, 2022 15:44

Add corefinement and integration over region bounded by a surface_mesh

2eaa7b6

fdrmrc force-pushed the cgal-corefinement branch from a1008c3 to 2eaa7b6 Compare May 7, 2022 18:08

luca-heltai requested changes May 7, 2022

View reviewed changes

luca-heltai mentioned this pull request May 7, 2022

Cgal triangulation dealii/dealii#13661

Merged

code review

50c74c2

fdrmrc force-pushed the cgal-corefinement branch from d093d3b to 50c74c2 Compare May 8, 2022 11:37

luca-heltai requested changes May 8, 2022

View reviewed changes

fdrmrc closed this Jun 6, 2022

	Quadrature<spacedim>
	Quadrature<CGALTriangulation::Space_dimension()>

	* @return Quadrature<spacedim> Collection of Quadratures. The sum of the weights equals to the volume of the polygonal region bounded by sm
	* @return Quadrature<spacedim> Collection of Quadratures. The sum of the weights
	* equals to the volume of the polygonal region bounded by sm

	constexpr unsigned int spacedim = 3;
	constexpr unsigned int spacedim = CGALPointType::Ambient_dimension::value;

Cgal corefinement #111

Cgal corefinement #111

Conversation

fdrmrc commented May 6, 2022 • edited

fdrmrc commented May 7, 2022 • edited

luca-heltai left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fdrmrc May 7, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luca-heltai commented May 7, 2022

fdrmrc commented May 8, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luca-heltai left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luca-heltai commented May 8, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fdrmrc commented Jun 6, 2022

fdrmrc commented May 6, 2022 •

edited

fdrmrc commented May 7, 2022 •

edited

fdrmrc May 7, 2022 •

edited