Skip to content

Commit

Permalink
Use AdditionalData and fix headers inclusion
Browse files Browse the repository at this point in the history
  • Loading branch information
fdrmrc committed Jun 6, 2022
1 parent bc9dbf7 commit b5418bc
Show file tree
Hide file tree
Showing 7 changed files with 179 additions and 141 deletions.
105 changes: 105 additions & 0 deletions include/deal.II/cgal/point_conversion.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_cgal_point_conversion_h
#define dealii_cgal_point_conversion_h

#include <deal.II/base/config.h>

#include <deal.II/base/point.h>

#ifdef DEAL_II_WITH_CGAL

# include <CGAL/Cartesian.h>
# include <CGAL/Exact_predicates_exact_constructions_kernel.h>
# include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
# include <CGAL/Simple_cartesian.h>


DEAL_II_NAMESPACE_OPEN
namespace CGALWrappers
{
/**
* Convert from a deal.II Point to any compatible CGAL point.
*
* @tparam CGALPointType Any of the CGAL point types
* @tparam dim Dimension of the point
* @param [in] p An input deal.II Point<dim>
* @return CGALPointType A CGAL point
*/
template <typename CGALPointType, int dim>
inline CGALPointType
dealii_point_to_cgal_point(const dealii::Point<dim> &p);

/**
* Convert from various CGAL point types to deal.II Point.
*
* @tparam dim Dimension of the point
* @tparam CGALPointType Any of the CGAL point types
* @param p An input CGAL point type
* @return dealii::Point<dim> The corresponding deal.II point.
*/
template <int dim, typename CGALPointType>
inline dealii::Point<dim>
cgal_point_to_dealii_point(const CGALPointType &p);


# ifndef DOXYGEN
// Template implementations

template <typename CGALPointType, int dim>
inline CGALPointType
dealii_point_to_cgal_point(const dealii::Point<dim> &p)
{
constexpr int cdim = CGALPointType::Ambient_dimension::value;
static_assert(dim <= cdim, "Only dim <= cdim supported");
if constexpr (cdim == 1)
return CGALPointType(p[0]);
else if constexpr (cdim == 2)
return CGALPointType(p[0], dim > 1 ? p[1] : 0);
else if constexpr (cdim == 3)
return CGALPointType(p[0], dim > 1 ? p[1] : 0, dim > 2 ? p[2] : 0);
else
Assert(false, dealii::ExcNotImplemented());
return CGALPointType();
}



template <int dim, typename CGALPointType>
inline dealii::Point<dim>
cgal_point_to_dealii_point(const CGALPointType &p)
{
constexpr int cdim = CGALPointType::Ambient_dimension::value;
if constexpr (dim == 1)
return dealii::Point<dim>(CGAL::to_double(p.x()));
else if constexpr (dim == 2)
return dealii::Point<dim>(CGAL::to_double(p.x()),
cdim > 1 ? CGAL::to_double(p.y()) : 0);
else if constexpr (dim == 3)
return dealii::Point<dim>(CGAL::to_double(p.x()),
cdim > 1 ? CGAL::to_double(p.y()) : 0,
cdim > 2 ? CGAL::to_double(p.z()) : 0);
else
Assert(false, dealii::ExcNotImplemented());
}
} // namespace CGALWrappers

# endif

DEAL_II_NAMESPACE_CLOSE

#endif
#endif
4 changes: 1 addition & 3 deletions include/deal.II/cgal/surface_mesh.h
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,10 @@

#include <deal.II/grid/tria.h>

#include <deal.II/cgal/utilities.h>

#ifdef DEAL_II_WITH_CGAL
# include <CGAL/Polygon_mesh_processing/stitch_borders.h>
# include <CGAL/Surface_mesh.h>

# include <deal.II/cgal/point_conversion.h>

DEAL_II_NAMESPACE_OPEN

Expand Down
4 changes: 2 additions & 2 deletions include/deal.II/cgal/triangulation.h
Original file line number Diff line number Diff line change
Expand Up @@ -257,8 +257,8 @@ namespace CGALWrappers
// A face is identified by a cell and by the index within the cell
// of the opposite vertex. Loop over vertices, and retain only those
// that belong to this face.
unsigned int j = 0;
for (unsigned int i = 0; i < 4; ++i)
int j = 0;
for (int i = 0; i < 4; ++i)
if (i != cgal_vertex_face_index)
dealii_face.vertices[j++] =
cgal_to_dealii_vertex_map[cgal_cell->vertex(i)];
Expand Down
63 changes: 0 additions & 63 deletions include/deal.II/cgal/utilities.h
Original file line number Diff line number Diff line change
Expand Up @@ -209,30 +209,6 @@ namespace CGALWrappers



/**
* Convert from a deal.II Point to any compatible CGAL point.
*
* @tparam CGALPointType Any of the CGAL point types
* @tparam dim Dimension of the point
* @param [in] p An input deal.II Point<dim>
* @return CGALPointType A CGAL point
*/
template <typename CGALPointType, int dim>
inline CGALPointType
dealii_point_to_cgal_point(const dealii::Point<dim> &p);

/**
* Convert from various CGAL point types to deal.II Point.
*
* @tparam dim Dimension of the point
* @tparam CGALPointType Any of the CGAL point types
* @param p An input CGAL point type
* @return dealii::Point<dim> The corresponding deal.II point.
*/
template <int dim, typename CGALPointType>
inline dealii::Point<dim>
cgal_point_to_dealii_point(const CGALPointType &p);

/**
* Given a closed CGAL::Surface_mesh, this function fills the
* region bounded by the surface with tets.
Expand Down Expand Up @@ -375,45 +351,6 @@ namespace CGALWrappers
// Template implementations
namespace CGALWrappers
{
template <typename CGALPointType, int dim>
inline CGALPointType
dealii_point_to_cgal_point(const dealii::Point<dim> &p)
{
constexpr int cdim = CGALPointType::Ambient_dimension::value;
static_assert(dim <= cdim, "Only dim <= cdim supported");
if constexpr (cdim == 1)
return CGALPointType(p[0]);
else if constexpr (cdim == 2)
return CGALPointType(p[0], dim > 1 ? p[1] : 0);
else if constexpr (cdim == 3)
return CGALPointType(p[0], dim > 1 ? p[1] : 0, dim > 2 ? p[2] : 0);
else
Assert(false, dealii::ExcNotImplemented());
return CGALPointType();
}



template <int dim, typename CGALPointType>
inline dealii::Point<dim>
cgal_point_to_dealii_point(const CGALPointType &p)
{
constexpr int cdim = CGALPointType::Ambient_dimension::value;
if constexpr (dim == 1)
return dealii::Point<dim>(CGAL::to_double(p.x()));
else if constexpr (dim == 2)
return dealii::Point<dim>(CGAL::to_double(p.x()),
cdim > 1 ? CGAL::to_double(p.y()) : 0);
else if constexpr (dim == 3)
return dealii::Point<dim>(CGAL::to_double(p.x()),
cdim > 1 ? CGAL::to_double(p.y()) : 0,
cdim > 2 ? CGAL::to_double(p.z()) : 0);
else
Assert(false, dealii::ExcNotImplemented());
}



template <typename C3t3>
void
cgal_surface_mesh_to_cgal_triangulation(
Expand Down
72 changes: 7 additions & 65 deletions include/deal.II/grid/grid_generator.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,25 +30,18 @@

#ifdef DEAL_II_WITH_CGAL
// All functions needed by the CGAL mesh generation utilities
<<<<<<< HEAD
# include <CGAL/Complex_2_in_triangulation_3.h>
# include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>
# include <CGAL/Implicit_surface_3.h>
=======
>>>>>>> 4b1e1f45a5 (Surface deal.II Tria to volumetric deal.II Tria)
# include <CGAL/Labeled_mesh_domain_3.h>
# include <CGAL/Mesh_complex_3_in_triangulation_3.h>
# include <CGAL/Mesh_criteria_3.h>
# include <CGAL/Mesh_triangulation_3.h>
<<<<<<< HEAD
# include <CGAL/Surface_mesh.h>
# include <CGAL/Surface_mesh_default_triangulation_3.h>
# include <CGAL/make_mesh_3.h>
# include <CGAL/make_surface_mesh.h>
=======
# include <CGAL/make_mesh_3.h>
# include <deal.II/cgal/surface_mesh.h>
>>>>>>> 4b1e1f45a5 (Surface deal.II Tria to volumetric deal.II Tria)
# include <deal.II/cgal/triangulation.h>
# include <deal.II/cgal/utilities.h>
#endif
Expand Down Expand Up @@ -1831,8 +1824,8 @@ namespace GridGenerator
CGALWrappers::AdditionalData<dim>{},
const Point<3> &interior_point = Point<3>(),
const double & outer_ball_radius = 1.0);
/**

/**
* Create a deal.II Triangulation<3> out of a deal.II Triangulation<2,3>
* by filling it with tetrahedrons.
*
Expand All @@ -1855,14 +1848,14 @@ namespace GridGenerator
*
* @param [in] surface_tria The input deal.II Triangulation<2,3>.
* @param [out] vol_tria The output deal.II Triangulation<3>.
* @param [in] cgal_args Additional parameters to pass to the CGAL::make_mesh_3
* function.
* @param[in] data Additional parameters to pass to the CGAL::make_mesh_3
* function and to the CGAL::make_surface_mesh functions
*/
template <typename... Args>
void
surface_mesh_to_volumetric_mesh(const Triangulation<2, 3> &surface_tria,
Triangulation<3> &vol_tria,
Args... cgal_args);
Triangulation<3> & vol_tria,
const CGALWrappers::AdditionalData<3> &data =
CGALWrappers::AdditionalData<3>{});
#endif
///@}

Expand Down Expand Up @@ -2876,57 +2869,6 @@ namespace GridGenerator
Assert(false, ExcImpossibleInDim(dim));
}
}



template <typename... Args>
void
surface_mesh_to_volumetric_mesh(const Triangulation<2, 3> &surface_tria,
Triangulation<3> &vol_tria,
Args... cgal_args)
{
Assert(
surface_tria.n_cells() > 0,
ExcMessage(
"The input triangulation cannot be empty when calling this function."));
Assert(
vol_tria.n_cells() == 0,
ExcMessage(
"The output triangulation must be empty when calling this function."));
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;

using Mesh_domain =
CGAL::Polyhedral_mesh_domain_with_features_3<K,
CGAL::Surface_mesh<Point_3>>;
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
CGAL::Default,
CGALWrappers::ConcurrencyTag>::type;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
using C3t3 =
CGAL::Mesh_complex_3_in_triangulation_3<Tr,
Mesh_domain::Corner_index,
Mesh_domain::Curve_index>;

CGAL::Surface_mesh<Point_3> mesh;
// This function "fills" the missing arrow of the following diagram.
// Tria<2,3> Tria<3>
// | ^
// | |
// | |
// | |
// V |
// CGAL::Surface_mesh -----------> CGAL::C3t3
CGALWrappers::dealii_tria_to_cgal_surface_mesh(surface_tria, mesh);
CGAL::Polygon_mesh_processing::triangulate_faces(mesh);
CGAL::Polygon_mesh_processing::stitch_borders(mesh);
Mesh_domain domain(mesh);
domain.detect_features();
Mesh_criteria criteria(cgal_args...);
const auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
CGALWrappers::cgal_triangulation_to_dealii_triangulation(cgal_triangulation,
vol_tria);
}
# endif

#endif
Expand Down
57 changes: 57 additions & 0 deletions source/grid/grid_generator.cc
Original file line number Diff line number Diff line change
Expand Up @@ -8639,6 +8639,63 @@ namespace GridGenerator
}
}


# ifdef DEAL_II_WITH_CGAL
void
surface_mesh_to_volumetric_mesh(const Triangulation<2, 3> &surface_tria,
Triangulation<3> & vol_tria,
const CGALWrappers::AdditionalData<3> &data)
{
Assert(
surface_tria.n_cells() > 0,
ExcMessage(
"The input triangulation cannot be empty when calling this function."));
Assert(
vol_tria.n_cells() == 0,
ExcMessage(
"The output triangulation must be empty when calling this function."));
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = K::Point_3;

using Mesh_domain =
CGAL::Polyhedral_mesh_domain_with_features_3<K,
CGAL::Surface_mesh<Point_3>>;
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain,
CGAL::Default,
CGALWrappers::ConcurrencyTag>::type;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
using C3t3 =
CGAL::Mesh_complex_3_in_triangulation_3<Tr,
Mesh_domain::Corner_index,
Mesh_domain::Curve_index>;

CGAL::Surface_mesh<Point_3> mesh;
// This function "fills" the missing arrow of the following diagram.
// Tria<2,3> Tria<3>
// | ^
// | |
// | |
// | |
// V |
// CGAL::Surface_mesh -----------> CGAL::C3t3
CGALWrappers::dealii_tria_to_cgal_surface_mesh(surface_tria, mesh);
CGAL::Polygon_mesh_processing::triangulate_faces(mesh);
CGAL::Polygon_mesh_processing::stitch_borders(mesh);
Mesh_domain domain(mesh);
domain.detect_features();
Mesh_criteria criteria(CGAL::parameters::facet_size = data.facet_size,
CGAL::parameters::facet_angle = data.facet_angle,
CGAL::parameters::facet_distance =
data.facet_distance,
CGAL::parameters::cell_radius_edge_ratio =
data.cell_radius_edge_ratio,
CGAL::parameters::cell_size = data.cell_size);
const auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
CGALWrappers::cgal_triangulation_to_dealii_triangulation(cgal_triangulation,
vol_tria);
}
# endif

} // namespace GridGenerator

// explicit instantiations
Expand Down

0 comments on commit b5418bc

Please sign in to comment.