Skip to content

Commit

Permalink
Move implementation of implicit_function() into .cc
Browse files Browse the repository at this point in the history
  • Loading branch information
@fdrmrc
fdrmrc committed Jun 8, 2022
1 parent 190dda6 commit ada3eb1
Show file tree
Hide file tree
Showing 4 changed files with 157 additions and 111 deletions.
37 changes: 32 additions & 5 deletions include/deal.II/cgal/triangulation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,16 +18,30 @@

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

#ifdef DEAL_II_WITH_CGAL
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>

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

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

# include <boost/hana.hpp>

# 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>
# 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>
# include <CGAL/Polyhedron_3.h>
# include <CGAL/Surface_mesh.h>
# include <CGAL/Surface_mesh_default_triangulation_3.h>
# include <CGAL/Triangulation_2.h>
# include <CGAL/Triangulation_3.h>
# include <CGAL/make_mesh_3.h>
# include <CGAL/make_surface_mesh.h>
# include <deal.II/cgal/surface_mesh.h>
# include <deal.II/cgal/utilities.h>

DEAL_II_NAMESPACE_OPEN
Expand Down Expand Up @@ -256,8 +270,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 Expand Up @@ -545,7 +559,20 @@ namespace CGALWrappers
}
triangulation.create_triangulation(vertices, cells, subcell_data);
}
# endif
} // namespace CGALWrappers
# endif // doxygen

DEAL_II_NAMESPACE_CLOSE
#else

DEAL_II_NAMESPACE_OPEN
namespace CGALWrappers
{
template <int dim>
struct AdditionalData
{};


} // namespace CGALWrappers

DEAL_II_NAMESPACE_CLOSE
Expand Down
107 changes: 1 addition & 106 deletions include/deal.II/grid/grid_generator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1765,7 +1765,6 @@ namespace GridGenerator
const std::string & grid_generator_function_name,
const std::string & grid_generator_function_arguments);

#ifdef DEAL_II_WITH_CGAL
/**
* Generate a Triangulation from the zero level set of an implicit function,
* using the CGAL library.
Expand Down Expand Up @@ -1823,7 +1822,6 @@ namespace GridGenerator
CGALWrappers::AdditionalData<dim>{},
const Point<3> &interior_point = Point<3>(),
const double & outer_ball_radius = 1.0);
#endif
///@}

/**
Expand Down Expand Up @@ -2733,110 +2731,7 @@ namespace GridGenerator
const double,
const bool);

# ifdef DEAL_II_WITH_CGAL
template <int dim>
void
implicit_function(Triangulation<dim, 3> &tria,
const Function<3> & dealii_implicit_function,
const CGALWrappers::AdditionalData<dim> &data,
const Point<3> & interior_point,
const double & outer_ball_radius)
{
Assert(dealii_implicit_function.n_components == 1,
ExcMessage(
"The implicit function must have exactly one component."));
Assert(dealii_implicit_function.value(interior_point) < 0,
ExcMessage(
"The implicit function must be negative at the interior point."));
Assert(outer_ball_radius > 0,
ExcMessage("The outer ball radius must be positive."));
Assert(tria.n_active_cells() == 0,
ExcMessage("The triangulation must be empty."));

if constexpr (dim == 3)
{
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using NumberType = K::FT;
using Point_3 = K::Point_3;
using Sphere_3 = K::Sphere_3;

using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
using Tr =
CGAL::Mesh_triangulation_3<Mesh_domain,
CGAL::Default,
CGALWrappers::ConcurrencyTag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;


auto cgal_implicit_function = [&](const Point_3 &p) {
return NumberType(
dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
};

Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
cgal_implicit_function,
K::Sphere_3(
Point_3(interior_point[0], interior_point[1], interior_point[2]),
outer_ball_radius * outer_ball_radius));

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);

auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
CGALWrappers::cgal_triangulation_to_dealii_triangulation(
cgal_triangulation, tria);
}
else if constexpr (dim == 2)
{
// default triangulation for Surface_mesher
using Tr = CGAL::Surface_mesh_default_triangulation_3;
using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
using GT = Tr::Geom_traits;
using Sphere_3 = GT::Sphere_3;
using Point_3 = GT::Point_3;
using FT = GT::FT;
typedef FT (*Function)(Point_3);
using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;


auto cgal_implicit_function = [&](const Point_3 &p) {
return FT(
dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
};

Surface_3 surface(cgal_implicit_function,
Sphere_3(Point_3(interior_point[0],
interior_point[1],
interior_point[2]),
outer_ball_radius * outer_ball_radius));

Tr tr;
C2t3 c2t3(tr);
Surface_mesh mesh;

CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
data.radius_bound,
data.distance_bound);
CGAL::make_surface_mesh(c2t3,
surface,
criteria,
CGAL::Non_manifold_tag());
CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
}
else
{
Assert(false, ExcImpossibleInDim(dim));
}
}
# endif


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



template <int dim>
void
implicit_function(Triangulation<dim, 3> &tria,
const Function<3> & dealii_implicit_function,
const CGALWrappers::AdditionalData<dim> &data,
const Point<3> & interior_point,
const double & outer_ball_radius)
{
# ifdef DEAL_II_WITH_CGAL
Assert(dealii_implicit_function.n_components == 1,
ExcMessage(
"The implicit function must have exactly one component."));
Assert(dealii_implicit_function.value(interior_point) < 0,
ExcMessage(
"The implicit function must be negative at the interior point."));
Assert(outer_ball_radius > 0,
ExcMessage("The outer ball radius must be positive."));
Assert(tria.n_active_cells() == 0,
ExcMessage("The triangulation must be empty."));

if constexpr (dim == 3)
{
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using NumberType = K::FT;
using Point_3 = K::Point_3;
using Sphere_3 = K::Sphere_3;

using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
using Tr =
CGAL::Mesh_triangulation_3<Mesh_domain,
CGAL::Default,
CGALWrappers::ConcurrencyTag>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr, int, int>;
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;


auto cgal_implicit_function = [&](const Point_3 &p) {
return NumberType(
dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
};

Mesh_domain domain = Mesh_domain::create_implicit_mesh_domain(
cgal_implicit_function,
K::Sphere_3(
Point_3(interior_point[0], interior_point[1], interior_point[2]),
outer_ball_radius * outer_ball_radius));

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);

auto cgal_triangulation = CGAL::make_mesh_3<C3t3>(domain, criteria);
CGALWrappers::cgal_triangulation_to_dealii_triangulation(
cgal_triangulation, tria);
}
else if constexpr (dim == 2)
{
// default triangulation for Surface_mesher
using Tr = CGAL::Surface_mesh_default_triangulation_3;
using C2t3 = CGAL::Complex_2_in_triangulation_3<Tr>;
using GT = Tr::Geom_traits;
using Sphere_3 = GT::Sphere_3;
using Point_3 = GT::Point_3;
using FT = GT::FT;
typedef FT (*Function)(Point_3);
using Surface_3 = CGAL::Implicit_surface_3<GT, Function>;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;


auto cgal_implicit_function = [&](const Point_3 &p) {
return FT(
dealii_implicit_function.value(Point<3>(p.x(), p.y(), p.z())));
};

Surface_3 surface(cgal_implicit_function,
Sphere_3(Point_3(interior_point[0],
interior_point[1],
interior_point[2]),
outer_ball_radius * outer_ball_radius));

Tr tr;
C2t3 c2t3(tr);
Surface_mesh mesh;

CGAL::Surface_mesh_default_criteria_3<Tr> criteria(data.angular_bound,
data.radius_bound,
data.distance_bound);
CGAL::make_surface_mesh(c2t3,
surface,
criteria,
CGAL::Non_manifold_tag());
CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, mesh);
CGALWrappers::cgal_surface_mesh_to_dealii_triangulation(mesh, tria);
}
else
{
Assert(false, ExcImpossibleInDim(dim));
}

# else

(void)tria;
(void)implicit_function;
(void)data;
(void)interior_point;
(void)outer_ball_radius;
AssertThrow(false, ExcMessage("This function needs CGAL to be installed."));

# endif
}
} // namespace GridGenerator

// explicit instantiations
Expand Down
8 changes: 8 additions & 0 deletions source/grid/grid_generator.inst.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,6 +228,14 @@ for (deal_II_dimension : DIMENSIONS)
hyper_ball_balanced<deal_II_dimension>(Triangulation<deal_II_dimension> &,
const Point<deal_II_dimension> &,
const double);

template void
implicit_function<deal_II_dimension>(
Triangulation<deal_II_dimension, 3> &,
const Function<3> &,
const CGALWrappers::AdditionalData<deal_II_dimension> &,
const Point<3> &,
const double &);
#endif
\}
}
Expand Down

0 comments on commit ada3eb1

Please sign in to comment.