From 071be319945e6742a6adf133d3f9949e07e9f3cf Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 10 Jun 2016 08:02:10 -0400 Subject: [PATCH] Add 3D for Triangulation and Point interface. --- include/deal.II/python/point_wrapper.h | 63 ++- .../deal.II/python/triangulation_wrapper.h | 165 +++++++ source/python/CMakeLists.txt | 1 + source/python/export_point.cc | 2 +- source/python/export_triangulation.cc | 221 +++------- source/python/point_wrapper.cc | 60 ++- source/python/triangulation_wrapper.cc | 412 ++++++++++++++++++ 7 files changed, 750 insertions(+), 174 deletions(-) create mode 100644 include/deal.II/python/triangulation_wrapper.h create mode 100644 source/python/triangulation_wrapper.cc diff --git a/include/deal.II/python/point_wrapper.h b/include/deal.II/python/point_wrapper.h index b3f081a786..73c40ac938 100644 --- a/include/deal.II/python/point_wrapper.h +++ b/include/deal.II/python/point_wrapper.h @@ -22,18 +22,79 @@ namespace PyDealII class PointWrapper { public: + /** + * This constructor exists only so that a std::vector of PointWrapper can be + * created. This constructor should not be used. + */ + PointWrapper(); + + /** + * Constructor. Use a list of double to initialize the Point. + */ PointWrapper(boost::python::list list); + + /** + * Destructpr. + */ + ~PointWrapper(); + + /** + * Return the first component of the Point. + */ double get_x(); + + /** + * Set the first component of the Point. + */ void set_x(double x); + + /** + * Return the second component of the Point. + */ double get_y(); + + /** + * Set the second component of the Point. + */ void set_y(double y); + + /** + * Return the third component of the Point. + */ double get_z(); + + /** + * Set the third component of the Point. + */ void set_z(double z); + /** + * Return a pointer that can be casted to the underlying Point. + */ + void *get_point(); + private: + /** + * Dimension of the Point. + */ int dim; - std::shared_ptr point; + + /** + * Pointer to the underlying Point object. + */ + void *point; }; + + + //--------------------- Inline functions ----------------------// + + + + inline + void *PointWrapper::get_point() + { + return point; + } } #endif diff --git a/include/deal.II/python/triangulation_wrapper.h b/include/deal.II/python/triangulation_wrapper.h new file mode 100644 index 0000000000..cfcd16daf5 --- /dev/null +++ b/include/deal.II/python/triangulation_wrapper.h @@ -0,0 +1,165 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__triangulation_wrapper_h +#define dealii__triangulation_wrapper_h + +#include +#include +#include + +namespace PyDealII +{ + class TriangulationWrapper + { + public: + /** + * Constructor. Takes a string @p dim with one of the following values + * "2D", "2d", "3D", or "3d". + */ + TriangulationWrapper(const std::string &dim); + + /** + * Destructor. + */ + ~TriangulationWrapper(); + + /** + * Return the number of active cells. + */ + unsigned int n_active_cells(); + + /** + * Generate a hyper cube (square in 2D and cube in 3D) with exactly one + * cell. + */ + void generate_hyper_cube(const double left = 0., + const double right = 1., + const bool colorize = false); + + /** + * Generate a simplex with (dim+1) vertices and mesh cells. + */ + void generate_simplex(boost::python::list &vertices); + + /** + * Same as hyper_cube but not only one cell is created but each coordinate + * direction is subdivided in @p repetitions cells. + */ + void generate_subdivided_hyper_cube(const unsigned int repetitions, + const double left = 0., + const double right = 1.); + + /** + * Generate a coordinate-parallel brick from the two diagonally opposite + * corners points @p1 and @p2. + */ + void generate_hyper_rectangle(PointWrapper &p1, + PointWrapper &p2, + const bool colorize = false); + + /** + * Generate a coordinate-parallel brick from the two diagonally opposite + * corners points @p1 and @p2. In direction i, repetitions[i] cells are + * created. + */ + void generate_subdivided_hyper_rectangle(boost::python::list &repetitions, + PointWrapper &p1, + PointWrapper &p2, + const bool colorize = false); + + /** + * Generate a hyperball, i.e. a circle or a ball around @p center with + * given @p radius. + */ + void generate_hyper_ball(PointWrapper ¢er, + const double radius = 1.); + + /** + * Shift each vertex of the Triangulation by the given @p shift_list. + */ + void shift(boost::python::list &shift_list); + + /** + * Given two triangulations, create the triangulation that contains the + * cells of both triangulations. + */ + void merge_triangulations(TriangulationWrapper &triangulation_1, + TriangulationWrapper &triangulation_2); + + /** + * Refine all the cells @p times times. + */ + void refine_global(const unsigned int times); + + /** + * Write the Triangulation in file. + */ + void save(const std::string &filename); + + /** + * Load the Triangulation from a file. + */ + void load(const std::string &filename); + + /** + * Return the dimension of the underlying Triangulation object. + */ + int get_dim(); + + /** + * Return a pointer that can be casted to the underlying Triangulation + * object. + */ + void *get_triangulation(); + + private: + /** + * Reset the underlying Triangulation object. + */ + void reset_triangulation(); + + /** + * Dimension of the underlying Triangulation object. + */ + int dim; + + /** + * Pointer that can be casted to the underlying Triangulation object. + */ + void *triangulation; + }; + + +//-------------------- inline functions -----------------------// + + + + inline + int TriangulationWrapper::get_dim() + { + return dim; + } + + + + inline + void *TriangulationWrapper::get_triangulation() + { + return triangulation; + } +} + +#endif diff --git a/source/python/CMakeLists.txt b/source/python/CMakeLists.txt index 7607ad8bec..6ee5c7261b 100644 --- a/source/python/CMakeLists.txt +++ b/source/python/CMakeLists.txt @@ -24,6 +24,7 @@ SET(_src export_point.cc export_triangulation.cc point_wrapper.cc + triangulation_wrapper.cc ) FILE(GLOB _header diff --git a/source/python/export_point.cc b/source/python/export_point.cc index b2f9d42cde..f908f1673f 100644 --- a/source/python/export_point.cc +++ b/source/python/export_point.cc @@ -22,7 +22,7 @@ namespace PyDealII { void export_point() { - boost::python::class_ ("Point",boost::python::init()) + boost::python::class_ ("Point", boost::python::init()) .add_property("x", &PointWrapper::get_x, &PointWrapper::set_x, "Get the x component of the point.") .add_property("y", &PointWrapper::get_y, &PointWrapper::set_y, "Get the y component of the point.") .add_property("z", &PointWrapper::get_z, &PointWrapper::set_z, "Get the z component of the point."); diff --git a/source/python/export_triangulation.cc b/source/python/export_triangulation.cc index 9fe9c4accc..9ee2fd177c 100644 --- a/source/python/export_triangulation.cc +++ b/source/python/export_triangulation.cc @@ -13,172 +13,79 @@ // // --------------------------------------------------------------------- -#include -#include #include -#include -#include -#include +#include namespace PyDealII { - template - unsigned int n_active_cells(const dealii::Triangulation &triangulation) - { - return triangulation.n_active_cells(); - } - - - - template - void generate_hyper_cube(dealii::Triangulation &triangulation, - const double left = 0., - const double right = 0., - const bool colorize = false) - { - dealii::GridGenerator::hyper_cube(triangulation, left, right, colorize); - } - - - - void generate_simplex(dealii::Triangulation<2> &triangulation, - dealii::Point<2> vertex_1, - dealii::Point<2> vertex_2, - dealii::Point<2> vertex_3) - { - std::vector> vertices(3); - vertices[0] = vertex_1; - vertices[1] = vertex_2; - vertices[2] = vertex_3; - - dealii::GridGenerator::simplex(triangulation, vertices); - } - - - - template - void generate_subdivided_hyper_cube(dealii::Triangulation &triangulation, - const unsigned int repetitions, - const double left = 0., - const double right = 1.) - { - dealii::GridGenerator::subdivided_hyper_cube(triangulation, repetitions, - left, right); - } - - - - template - void generate_hyper_rectangle(dealii::Triangulation &triangulation, - const dealii::Point &p1, - const dealii::Point &p2, - const bool colorize = false) - { - dealii::GridGenerator::hyper_rectangle(triangulation, p1, p2, colorize); - } - - - - void generate_subdivided_hyper_rectangle(dealii::Triangulation<2> &triangulation, - const unsigned int repetition_x, - const unsigned int repetition_y, - const dealii::Point<2> &p1, - const dealii::Point<2> &p2, - const bool colorize = false) - { - std::vector repetitions(2); - repetitions[0] = repetition_x; - repetitions[1] = repetition_y; - dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, repetitions, - p1, p2, colorize); - } - - - - template - void generate_hyper_ball(dealii::Triangulation &triangulation, - const dealii::Point ¢er = dealii::Point(), - const double radius = 1.) - { - dealii::GridGenerator::hyper_ball(triangulation, center, radius); - } - - - - template - void save(const dealii::Triangulation &triangulation, - const std::string filename) - { - std::ofstream ofs(filename); - boost::archive::text_oarchive oa(ofs); - oa << triangulation; - } - - - - template - void load(dealii::Triangulation &triangulation, - const std::string filename) - { - std::ifstream ifs(filename); - boost::archive::text_iarchive ia(ifs); - ia >> triangulation; - } - - - // Macro to enable default arguments - BOOST_PYTHON_FUNCTION_OVERLOADS(generate_hyper_cube_overloads, generate_hyper_cube, 1, 4) - BOOST_PYTHON_FUNCTION_OVERLOADS(generate_subdivided_hyper_cube_overloads, - generate_subdivided_hyper_cube, 2, 4) - BOOST_PYTHON_FUNCTION_OVERLOADS(generate_hyper_rectangle_overloads, - generate_hyper_rectangle, 3, 4) - BOOST_PYTHON_FUNCTION_OVERLOADS(generate_subdivided_hyper_rectangle_overloads, - generate_subdivided_hyper_rectangle, 5, 6) - BOOST_PYTHON_FUNCTION_OVERLOADS(generate_hyper_ball_overloads, generate_hyper_ball, 1, 3) + BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_hyper_cube_overloads, generate_hyper_cube, 0, 3) + BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_subdivided_hyper_cube_overloads, + generate_subdivided_hyper_cube, 1, 3) + BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_hyper_rectangle_overloads, + generate_hyper_rectangle, 2, 3) + BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_subdivided_hyper_rectangle_overloads, + generate_subdivided_hyper_rectangle, 3, 4) + BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_hyper_ball_overloads, generate_hyper_ball, 1, 2) void export_triangulation() { - boost::python::class_> ("Triangulation") - .def("n_active_cells", n_active_cells<2>, - "Return the number of active cells", - boost::python::args("self")) - .def("generate_hyper_cube", generate_hyper_cube<2>, - generate_hyper_cube_overloads( - boost::python::args("self", "left", "right", "colorize"), - "Generate a hyper_cube.")) - .def("generate_simplex", generate_simplex, - "Generate a simplex.", - boost::python::args("self", "vertex_1", "vertex_2", "vertex_3")) - .def("generate_subdivided_hyper_cube", generate_subdivided_hyper_cube<2>, - generate_subdivided_hyper_cube_overloads( - boost::python::args("self", "repetitions", "left", "right"), - "Generate a subdivided hyper_cube.")) - .def("generate_hyper_rectangle", generate_hyper_rectangle<2>, - generate_hyper_rectangle_overloads( - boost::python::args("self", "p1", "p2", "colorize"), - "Generate a hyper_rectangle.")) - .def("generate_subdivided_hyper_rectangle", - generate_subdivided_hyper_rectangle, - generate_subdivided_hyper_rectangle_overloads( - boost::python::args("self", "repetition_x", - "repetition_y","p1", "p2", - "colorize"), - "Generate a subdivided hyper_rectangle.")) - .def("generate_hyper_ball", - generate_hyper_ball<2>, - generate_hyper_ball_overloads( - boost::python::args("self", "center", "radius"), - "Generate a hyper_ball.")) - .def("refine_global", &dealii::Triangulation<2>::refine_global, - "Refine the mesh uniformly.", - boost::python::args("self", "times")) - .def("save", save<2>, "Serialize and save the triangulation.", - boost::python::args("self", "filename")) - .def("load", load<2>, "Load and deserialize a triangulation.", - boost::python::args("self", "filename")); + boost::python::class_("Triangulation", boost::python::init()) + .def("n_active_cells", + &TriangulationWrapper::n_active_cells, + "Return the number of active cells", + boost::python::args("self")) + .def("generate_hyper_cube", + &TriangulationWrapper::generate_hyper_cube, + generate_hyper_cube_overloads( + boost::python::args("self", "left", "right", "colorize"), + "Generate a hyper_cube (square in 2D and cube in 3D) with exactly one cell.")) + .def("generate_simplex", + &TriangulationWrapper::generate_simplex, + "Generate a simplex with (dim+1) vertices and mesh cells.", + boost::python::args("self", "vertices")) + .def("generate_subdivided_hyper_cube", + &TriangulationWrapper::generate_subdivided_hyper_cube, + generate_subdivided_hyper_cube_overloads( + boost::python::args("self", "repetitions", "left", "right"), + "Same as hyper_cube but not only one cells is created but each coordinate direction is subdivdided in repetitions cells.")) + .def("generate_hyper_rectangle", + &TriangulationWrapper::generate_hyper_rectangle, + generate_hyper_rectangle_overloads( + boost::python::args("self", "p1", "p2", "colorize"), + "Generate a coordinate-parallel brick from the two diagonally opposite corners points p1 and p2.")) + .def("generate_subdivided_hyper_rectangle", + &TriangulationWrapper::generate_subdivided_hyper_rectangle, + generate_subdivided_hyper_rectangle_overloads( + boost::python::args("self", "repetitions", + "p1", "p2", "colorize"), + "Generate a coordinate-parallel brick from the two diagonally opposite corners point @p1 and @p2. In direction i, repetitions[i] cells are created.")) + .def("generate_hyper_ball", + &TriangulationWrapper::generate_hyper_ball, + generate_hyper_ball_overloads( + boost::python::args("self", "center", "radius"), + "Generate a hyperball, i.e. a circle or a ball around @p center with given @p radius.")) + .def("shift", + &TriangulationWrapper::shift, + "Shift every vertex of the Triangulation by the given shift vector.", + boost::python::args("self", "shift")) + .def("merge_triangulations", + &TriangulationWrapper::merge_triangulations, + "Given two triangulations, create the triangulation that contains the cells of both triangulations.", + boost::python::args("self", "triangulation_1", "triangulation_2")) + .def("refine_global", + &TriangulationWrapper::refine_global, + "Refine all the cells times times.", + boost::python::args("self", "times")) + .def("save", + &TriangulationWrapper::save, + "Write the Triangulation to a file.", + boost::python::args("self", "filename")) + .def("load", + &TriangulationWrapper::load, + "Load the Triangulation from a file.", + boost::python::args("self", "filename")); } - } diff --git a/source/python/point_wrapper.cc b/source/python/point_wrapper.cc index 07529ac726..8459bd5ddf 100644 --- a/source/python/point_wrapper.cc +++ b/source/python/point_wrapper.cc @@ -18,16 +18,23 @@ namespace PyDealII { + PointWrapper::PointWrapper() + : + dim(-1), + point(nullptr) + {} + + PointWrapper::PointWrapper(boost::python::list coord) { dim = boost::python::len(coord); if (dim == 2) - point.reset(new dealii::Point<2> (boost::python::extract(coord[0]), - boost::python::extract(coord[1]))); + point = new dealii::Point<2> (boost::python::extract(coord[0]), + boost::python::extract(coord[1])); else if (dim == 3) - point.reset(new dealii::Point<3> (boost::python::extract(coord[0]), - boost::python::extract(coord[1]), - boost::python::extract(coord[2]))); + point = new dealii::Point<3> (boost::python::extract(coord[0]), + boost::python::extract(coord[1]), + boost::python::extract(coord[2])); else AssertThrow(false, dealii::ExcMessage("The list of coordinates must contain two or three elements.")); @@ -35,12 +42,35 @@ namespace PyDealII + PointWrapper::~PointWrapper() + { + if (point != nullptr) + { + if (dim == 2) + { + // We cannot call delete on a void pointer so cast the void pointer back + // first. + dealii::Point<2> *tmp = static_cast*>(point); + delete tmp; + } + else + { + dealii::Point<3> *tmp = static_cast*>(point); + delete tmp; + } + point = nullptr; + } + dim = -1; + } + + + double PointWrapper::get_x() { if (dim == 2) - return (*std::static_pointer_cast>(point))(0); + return (*static_cast*>(point))(0); else - return (*std::static_pointer_cast>(point))(0); + return (*static_cast*>(point))(0); } @@ -48,9 +78,9 @@ namespace PyDealII void PointWrapper::set_x(double x) { if (dim == 2) - (*std::static_pointer_cast>(point))(0) = x; + (*static_cast*>(point))(0) = x; else - (*std::static_pointer_cast>(point))(0) = x; + (*static_cast*>(point))(0) = x; } @@ -58,9 +88,9 @@ namespace PyDealII double PointWrapper::get_y() { if (dim == 2) - return (*std::static_pointer_cast>(point))(1); + return (*static_cast*>(point))(1); else - return (*std::static_pointer_cast>(point))(1); + return (*static_cast*>(point))(1); } @@ -68,9 +98,9 @@ namespace PyDealII void PointWrapper::set_y(double y) { if (dim == 2) - (*std::static_pointer_cast>(point))(1) = y; + (*static_cast*>(point))(1) = y; else - (*std::static_pointer_cast>(point))(1) = y; + (*static_cast*>(point))(1) = y; } @@ -78,7 +108,7 @@ namespace PyDealII double PointWrapper::get_z() { if (dim == 3) - return (*std::static_pointer_cast>(point))(2); + return (*static_cast*>(point))(2); else AssertThrow(false, dealii::ExcMessage("The z coordinate is only available for three-dimensional points")); @@ -89,7 +119,7 @@ namespace PyDealII void PointWrapper::set_z(double z) { if (dim == 3) - (*std::static_pointer_cast>(point))(2) = z; + (*static_cast*>(point))(2) = z; else AssertThrow(false, dealii::ExcMessage("The z coordinate is only available for three-dimensional points")); diff --git a/source/python/triangulation_wrapper.cc b/source/python/triangulation_wrapper.cc new file mode 100644 index 0000000000..98717593b4 --- /dev/null +++ b/source/python/triangulation_wrapper.cc @@ -0,0 +1,412 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include + +namespace PyDealII +{ + TriangulationWrapper::TriangulationWrapper(const std::string &dimension) + : + triangulation(nullptr) + { + if ((dimension.compare("2D")==0) || (dimension.compare("2d")==0)) + { + dim = 2; + } + else if ((dimension.compare("3D")==0) || (dimension.compare("3d")==0)) + { + dim = 3; + } + else + AssertThrow(false, dealii::ExcMessage("Dimension needs to be 2D or 3D")); + } + + + + TriangulationWrapper::~TriangulationWrapper() + { + if (triangulation != nullptr) + { + if (dim == 2) + { + // We cannot call delete on a void pointer so cast the void pointer back + // first. + dealii::Triangulation<2> *tmp = static_cast*>(triangulation); + delete tmp; + } + else + { + dealii::Triangulation<3> *tmp = static_cast*>(triangulation); + delete tmp; + } + triangulation = nullptr; + } + dim = -1; + } + + + + unsigned int TriangulationWrapper::n_active_cells() + { + if (dim == 2) + { + return (*static_cast*>(triangulation)).n_active_cells(); + } + else + return (*static_cast*>(triangulation)).n_active_cells(); + } + + + + void TriangulationWrapper::generate_hyper_cube(const double left, + const double right, + const bool colorize) + { + // Create a new Triangulation + reset_triangulation(); + + if (dim == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_cube(*tria, left, right, colorize); + } + else + { + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_cube(*tria, left, right, colorize); + } + } + + + + void TriangulationWrapper::generate_simplex(boost::python::list &vertices) + { + // Create a new Triangulation + reset_triangulation(); + + // Extract the PointWrapper object from the python list + std::vector wrapped_points(dim); + for (int i=0; i(vertices[i]); + + if (dim == 2) + { + // Cast the PointWrapper objects to Point + std::vector> points(dim); + for (int i=0; i*>((wrapped_points[i]).get_point())); + + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::simplex(*tria, points); + } + else + { + std::vector> points(dim); + for (int i=0; i*>((wrapped_points[i]).get_point())); + + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::simplex(*tria, points); + } + } + + + + void TriangulationWrapper::generate_subdivided_hyper_cube(const unsigned int repetitions, + const double left, + const double right) + { + // Create a new Triangulation + reset_triangulation(); + + if (dim == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::subdivided_hyper_cube(*tria, repetitions, left, right); + } + else + { + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::subdivided_hyper_cube(*tria, repetitions, left, right); + } + } + + + + void TriangulationWrapper::generate_hyper_rectangle(PointWrapper &p1, + PointWrapper &p2, + const bool colorize) + { + // Create a new Triangulation + reset_triangulation(); + + if (dim == 2) + { + // Cast the PointWrapper object to Point + dealii::Point<2> point_1 = *(static_cast*>(p1.get_point())); + dealii::Point<2> point_2 = *(static_cast*>(p2.get_point())); + + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_rectangle(*tria, point_1, point_2, colorize); + } + else + { + dealii::Point<3> point_1 = *(static_cast*>(p1.get_point())); + dealii::Point<3> point_2 = *(static_cast*>(p2.get_point())); + + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_rectangle(*tria, point_1, point_2, colorize); + } + + } + + + + void TriangulationWrapper::generate_subdivided_hyper_rectangle(boost::python::list &repetition_list, + PointWrapper &p1, + PointWrapper &p2, + const bool colorize) + { + // Create a new Triangulation + reset_triangulation(); + + // Extract the repetitions from the python list + std::vector repetitions(dim); + for (int i=0; i(repetition_list[i]); + + if (dim == 2) + { + // Cast the PointWrapper object to Point + dealii::Point<2> point_1 = *(static_cast*>(p1.get_point())); + dealii::Point<2> point_2 = *(static_cast*>(p2.get_point())); + + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::subdivided_hyper_rectangle(*tria, repetitions, point_1, + point_2, colorize); + } + else + { + dealii::Point<3> point_1 = *(static_cast*>(p1.get_point())); + dealii::Point<3> point_2 = *(static_cast*>(p2.get_point())); + + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::subdivided_hyper_rectangle(*tria, repetitions, point_1, + point_2, colorize); + } + } + + + + void TriangulationWrapper::generate_hyper_ball(PointWrapper ¢er, + const double radius) + { + // Create a new Triangulation + reset_triangulation(); + + if (dim == 2) + { + // Cast the PointWrapper object to Point + dealii::Point<2> center_point = *(static_cast*>( + center.get_point())); + + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_ball(*tria, center_point, radius); + } + else + { + dealii::Point<3> center_point = *(static_cast*>( + center.get_point())); + + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridGenerator::hyper_ball(*tria, center_point, radius); + } + } + + + + void TriangulationWrapper::shift(boost::python::list &shift_list) + { + if (dim == 2) + { + // Extract the shift vector from the python list + dealii::Tensor<1,2> shift_vector; + for (int i=0; i(shift_list[i]); + + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + dealii::GridTools::shift(shift_vector, *tria); + } + else + { + dealii::Tensor<1,3> shift_vector; + for (int i=0; i(shift_list[i]); + + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + dealii::GridTools::shift(shift_vector, *tria); + } + } + + + + void TriangulationWrapper::merge_triangulations(TriangulationWrapper &triangulation_1, + TriangulationWrapper &triangulation_2) + { + AssertThrow(triangulation_1.get_dim() == triangulation_2.get_dim(), + dealii::ExcMessage("Triangulation_1 and Triangulation_2 should have the same dimension.")); + + // Create a new Triangulation + reset_triangulation(); + + if (triangulation_1.get_dim() == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + tria = new dealii::Triangulation<2>(); + dealii::Triangulation<2> *tria_1 = + static_cast*>(triangulation_1.get_triangulation()); + dealii::Triangulation<2> *tria_2 = + static_cast*>(triangulation_2.get_triangulation()); + dealii::GridGenerator::merge_triangulations(*tria_1, *tria_2, *tria); + // We need to reassign tria to triangulation because tria was cleared + // inside merge_triangulations. + triangulation = tria; + } + else + { + TriangulationWrapper merged_triangulation("3D"); + dealii::Triangulation<3> *tria = + static_cast*>(merged_triangulation.get_triangulation()); + tria = new dealii::Triangulation<3>(); + dealii::Triangulation<3> *tria_1 = + static_cast*>(triangulation_1.get_triangulation()); + dealii::Triangulation<3> *tria_2 = + static_cast*>(triangulation_2.get_triangulation()); + dealii::GridGenerator::merge_triangulations(*tria_1, *tria_2, *tria); + triangulation = tria; + } + } + + + + void TriangulationWrapper::save(const std::string &filename) + { + std::ofstream ofs(filename); + boost::archive::text_oarchive oa(ofs); + + if (dim == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + + oa << *tria; + } + else + { + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + + oa << *tria; + } + } + + + + void TriangulationWrapper::load(const std::string &filename) + { + std::ifstream ifs(filename); + boost::archive::text_iarchive ia(ifs); + + if (dim == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + + ia >> *tria; + } + else + { + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + + ia >> *tria; + } + } + + + + void TriangulationWrapper::refine_global(const unsigned int times) + { + if (dim == 2) + { + dealii::Triangulation<2> *tria = + static_cast*>(triangulation); + tria->refine_global(times); + } + else + { + dealii::Triangulation<3> *tria = + static_cast*>(triangulation); + tria->refine_global(times); + } + } + + + + void TriangulationWrapper::reset_triangulation() + { + if (dim == 2) + { + if (triangulation != nullptr) + { + // We cannot call delete on a void pointer so cast the void pointer back + // first. + dealii::Triangulation<2> *tmp = static_cast*>(triangulation); + delete tmp; + } + triangulation = new dealii::Triangulation<2>(); + } + else + { + if (triangulation != nullptr) + { + dealii::Triangulation<3> *tmp = static_cast*>(triangulation); + delete tmp; + } + triangulation = new dealii::Triangulation<3> (); + } + } +} -- 2.39.5