From b50b86d6551f27dbbd072bd8fcd7b8c8f52dfe95 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 17 Jun 2016 10:31:39 -0400 Subject: [PATCH] Allow Triangulation to set spacedim different than dim. --- include/deal.II/python/point_wrapper.h | 2 +- .../deal.II/python/triangulation_wrapper.h | 31 ++ source/python/export_triangulation.cc | 3 +- source/python/point_wrapper.cc | 6 +- source/python/triangulation_wrapper.cc | 270 ++++++++++++------ tests/python/triangulation_wrapper.py | 34 +-- 6 files changed, 239 insertions(+), 107 deletions(-) diff --git a/include/deal.II/python/point_wrapper.h b/include/deal.II/python/point_wrapper.h index 2206515a20..7fbeafc112 100644 --- a/include/deal.II/python/point_wrapper.h +++ b/include/deal.II/python/point_wrapper.h @@ -45,7 +45,7 @@ namespace PyDealII /** * Assignment operator. The dimension of the point is changed if it is - * different that the one of @p other. + * different than the one of @p other. */ PointWrapper &operator= (const PointWrapper &other); diff --git a/include/deal.II/python/triangulation_wrapper.h b/include/deal.II/python/triangulation_wrapper.h index 2ac14b3bfc..df4cb28617 100644 --- a/include/deal.II/python/triangulation_wrapper.h +++ b/include/deal.II/python/triangulation_wrapper.h @@ -31,6 +31,14 @@ namespace PyDealII */ TriangulationWrapper(const std::string &dim); + /** + * Constructor. Takes a string @p dim with one of the following values + * "2D", "2d", "3D", or "3d" and a string @p spacedim with one of the + * following values "2D", "2d", "3D", or "3d". The dimension of @p spacedim + * must be larger than the dimension of @p dim + */ + TriangulationWrapper(const std::string &dim, const std::string &spacedim); + /** * Destructor. */ @@ -119,6 +127,11 @@ namespace PyDealII */ int get_dim() const; + /** + * Return the space dimension of the underlying Triangulation object. + */ + int get_spacedim() const; + /** * Return a pointer that can be casted to the underlying Triangulation * object. @@ -126,11 +139,21 @@ namespace PyDealII void *get_triangulation(); private: + /** + * Helper function for the contructors. + */ + void setup(const std::string &dimension, const std::string &spacedimension); + /** * Dimension of the underlying Triangulation object. */ int dim; + /** + * Space dimension of the underlying Triangulation object. + */ + int spacedim; + /** * Pointer that can be casted to the underlying Triangulation object. */ @@ -150,6 +173,14 @@ namespace PyDealII + inline + int TriangulationWrapper::get_spacedim() const + { + return spacedim; + } + + + inline void *TriangulationWrapper::get_triangulation() { diff --git a/source/python/export_triangulation.cc b/source/python/export_triangulation.cc index f6a316fa3b..27628033d8 100644 --- a/source/python/export_triangulation.cc +++ b/source/python/export_triangulation.cc @@ -20,7 +20,7 @@ namespace PyDealII { // Macro to enable default arguments - BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_hyper_cube_overloads, generate_hyper_cube, 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, @@ -33,6 +33,7 @@ namespace PyDealII void export_triangulation() { boost::python::class_("Triangulation", boost::python::init()) + .def(boost::python::init()) .def("n_active_cells", &TriangulationWrapper::n_active_cells, "Return the number of active cells", diff --git a/source/python/point_wrapper.cc b/source/python/point_wrapper.cc index 87cc48090e..cf17461706 100644 --- a/source/python/point_wrapper.cc +++ b/source/python/point_wrapper.cc @@ -156,6 +156,9 @@ namespace PyDealII { dim = other.dim; + AssertThrow(other.point != nullptr, + dealii::ExcMessage("Underlying point does not exist.")); + if (dim == 2) { dealii::Point<2> *other_point = static_cast*>(other.point); @@ -167,6 +170,7 @@ namespace PyDealII point = new dealii::Point<3> ((*other_point)[0], (*other_point)[1], (*other_point)[2]); } else - point = nullptr; + AssertThrow(false, + dealii::ExcMessage("The dimension of the point should be 2 or 3.")); } } diff --git a/source/python/triangulation_wrapper.cc b/source/python/triangulation_wrapper.cc index 915a830553..6f16f30067 100644 --- a/source/python/triangulation_wrapper.cc +++ b/source/python/triangulation_wrapper.cc @@ -25,14 +25,14 @@ namespace PyDealII { namespace internal { - template + template void generate_hyper_cube(const double left, const double right, const bool colorize, void *triangulation) { - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); tria->clear(); dealii::GridGenerator::hyper_cube(*tria, left, right, colorize); } @@ -56,21 +56,21 @@ namespace PyDealII - template + template void generate_subdivided_hyper_cube(const unsigned int repetitions, const double left, const double right, void *triangulation) { - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); tria->clear(); dealii::GridGenerator::subdivided_hyper_cube(*tria, repetitions, left, right); } - template + template void generate_hyper_rectangle(PointWrapper &p1, PointWrapper &p2, const bool colorize, @@ -84,15 +84,15 @@ namespace PyDealII dealii::Point point_1 = *(static_cast*>(p1.get_point())); dealii::Point point_2 = *(static_cast*>(p2.get_point())); - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); tria->clear(); dealii::GridGenerator::hyper_rectangle(*tria, point_1, point_2, colorize); } - template + template void generate_subdivided_hyper_rectangle(const std::vector &repetitions, PointWrapper &p1, PointWrapper &p2, @@ -107,8 +107,8 @@ namespace PyDealII dealii::Point point_1 = *(static_cast*>(p1.get_point())); dealii::Point point_2 = *(static_cast*>(p2.get_point())); - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); tria->clear(); dealii::GridGenerator::subdivided_hyper_rectangle(*tria, repetitions, point_1, point_2, colorize); @@ -133,34 +133,34 @@ namespace PyDealII - template + template void shift(boost::python::list &shift_list, void *triangulation) { // Extract the shift vector from the python list - dealii::Tensor<1,dim> shift_vector; - for (int i=0; i shift_vector; + for (int i=0; i(shift_list[i]); - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); dealii::GridTools::shift(shift_vector, *tria); } - template + template void merge_triangulations(TriangulationWrapper &triangulation_1, TriangulationWrapper &triangulation_2, void *triangulation) { - dealii::Triangulation *tria = - static_cast*>(triangulation); + dealii::Triangulation *tria = + static_cast*>(triangulation); tria->clear(); - dealii::Triangulation *tria_1 = - static_cast*>(triangulation_1.get_triangulation()); - dealii::Triangulation *tria_2 = - static_cast*>(triangulation_2.get_triangulation()); + dealii::Triangulation *tria_1 = + static_cast*>(triangulation_1.get_triangulation()); + dealii::Triangulation *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. @@ -171,43 +171,53 @@ namespace PyDealII TriangulationWrapper::TriangulationWrapper(const std::string &dimension) - : - triangulation(nullptr) { if ((dimension.compare("2D")==0) || (dimension.compare("2d")==0)) - { - dim = 2; - triangulation = new dealii::Triangulation<2>(); - } + setup("2D", "2D"); else if ((dimension.compare("3D")==0) || (dimension.compare("3d")==0)) - { - dim = 3; - triangulation = new dealii::Triangulation<3>(); - } + setup("3D", "3D"); else AssertThrow(false, dealii::ExcMessage("Dimension needs to be 2D or 3D")); } + TriangulationWrapper::TriangulationWrapper(const std::string &dimension, + const std::string &spacedimension) + { + setup(dimension, spacedimension); + } + + + TriangulationWrapper::~TriangulationWrapper() { if (triangulation != nullptr) + { + if (dim == 2) { - 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; - } + if (spacedim == 2) + { + // We cannot call delete on a void pointer so cast the void pointer back + // first. + dealii::Triangulation<2,2> *tmp = + static_cast*>(triangulation); + delete tmp; + } else - { - dealii::Triangulation<3> *tmp = static_cast*>(triangulation); - delete tmp; - } - triangulation = nullptr; + { + dealii::Triangulation<2,3> *tmp = + static_cast*>(triangulation); + delete tmp; + } } + else + { + dealii::Triangulation<3,3> *tmp = static_cast*>(triangulation); + delete tmp; + } + triangulation = nullptr; + } dim = -1; } @@ -215,10 +225,12 @@ namespace PyDealII unsigned int TriangulationWrapper::n_active_cells() const { - if (dim == 2) - return (*static_cast*>(triangulation)).n_active_cells(); + if ((dim == 2) && (spacedim == 2)) + return (*static_cast*>(triangulation)).n_active_cells(); + else if ((dim == 2) && (spacedim == 3)) + return (*static_cast*>(triangulation)).n_active_cells(); else - return (*static_cast*>(triangulation)).n_active_cells(); + return (*static_cast*>(triangulation)).n_active_cells(); } @@ -227,10 +239,12 @@ namespace PyDealII const double right, const bool colorize) { - if (dim == 2) - internal::generate_hyper_cube<2>(left, right, colorize, triangulation); + if ((dim == 2) && (spacedim == 2)) + internal::generate_hyper_cube<2,2>(left, right, colorize, triangulation); + else if ((dim == 2) && (spacedim == 3)) + internal::generate_hyper_cube<2,3>(left, right, colorize, triangulation); else - internal::generate_hyper_cube<3>(left, right, colorize, triangulation); + internal::generate_hyper_cube<3,3>(left, right, colorize, triangulation); } @@ -239,6 +253,8 @@ namespace PyDealII { AssertThrow(boost::python::len(vertices) == dim+1, dealii::ExcMessage("The number of vertices should be equal to dim+1.")); + AssertThrow(dim == spacedim, + dealii::ExcMessage("This function is only implemented for dim equal to spacedim.")); // Extract the PointWrapper object from the python list std::vector wrapped_points(dim+1); for (int i=0; i(repetitions, left, right, triangulation); + if ((dim == 2) && (spacedim == 2)) + internal::generate_subdivided_hyper_cube<2,2>(repetitions, left, right, triangulation); + else if ((dim == 2) && (spacedim == 3)) + internal::generate_subdivided_hyper_cube<2,3>(repetitions, left, right, triangulation); else - internal::generate_subdivided_hyper_cube<3>(repetitions, left, right, triangulation); + internal::generate_subdivided_hyper_cube<3,3>(repetitions, left, right, triangulation); } @@ -272,10 +290,12 @@ namespace PyDealII PointWrapper &p2, const bool colorize) { - if (dim == 2) - internal::generate_hyper_rectangle<2>(p1, p2, colorize, triangulation); + if ((dim == 2) && (spacedim == 2)) + internal::generate_hyper_rectangle<2,2>(p1, p2, colorize, triangulation); + else if ((dim == 2) && (spacedim == 3)) + internal::generate_hyper_rectangle<2,3>(p1, p2, colorize, triangulation); else - internal::generate_hyper_rectangle<3>(p1, p2, colorize, triangulation); + internal::generate_hyper_rectangle<3,3>(p1, p2, colorize, triangulation); } @@ -293,12 +313,15 @@ namespace PyDealII for (int i=0; i(repetition_list[i]); - if (dim == 2) - internal::generate_subdivided_hyper_rectangle<2>(repetitions, p1, p2, colorize, - triangulation); + if ((dim == 2) && (spacedim == 2)) + internal::generate_subdivided_hyper_rectangle<2,2>(repetitions, p1, p2, colorize, + triangulation); + else if ((dim == 2) && (spacedim == 3)) + internal::generate_subdivided_hyper_rectangle<2,3>(repetitions, p1, p2, colorize, + triangulation); else - internal::generate_subdivided_hyper_rectangle<3>(repetitions, p1, p2, colorize, - triangulation); + internal::generate_subdivided_hyper_rectangle<3,3>(repetitions, p1, p2, colorize, + triangulation); } @@ -306,6 +329,8 @@ namespace PyDealII void TriangulationWrapper::generate_hyper_ball(PointWrapper ¢er, const double radius) { + AssertThrow(dim == spacedim, + dealii::ExcMessage("This function is only implemented for dim equal to spacedim.")); if (dim == 2) internal::generate_hyper_ball<2>(center, radius, triangulation); else @@ -316,10 +341,14 @@ namespace PyDealII void TriangulationWrapper::shift(boost::python::list &shift_list) { - if (dim == 2) - internal::shift<2>(shift_list, triangulation); + AssertThrow(boost::python::len(shift_list) == spacedim, + dealii::ExcMessage("Size of the shift vector is not equal to spacedim.")); + if ((dim == 2) && (spacedim == 2)) + internal::shift<2,2>(shift_list, triangulation); + else if ((dim ==2) && (spacedim == 3)) + internal::shift<2,3>(shift_list, triangulation); else - internal::shift<3>(shift_list, triangulation); + internal::shift<3,3>(shift_list, triangulation); } @@ -329,11 +358,43 @@ namespace PyDealII { AssertThrow(triangulation_1.get_dim() == triangulation_2.get_dim(), dealii::ExcMessage("Triangulation_1 and Triangulation_2 should have the same dimension.")); + AssertThrow(dim == triangulation_2.get_dim(), + dealii::ExcMessage("Triangulation and Triangulation_2 should have the same dimension.")); + AssertThrow(triangulation_1.get_spacedim() == triangulation_2.get_spacedim(), + dealii::ExcMessage("Triangulation_1 and Triangulation_2 should have the same space dimension.")); + AssertThrow(spacedim == triangulation_2.get_spacedim(), + dealii::ExcMessage("Triangulation and Triangulation_2 should have the same space dimension.")); + + if ((triangulation_1.get_dim() == 2) && (triangulation_1.get_spacedim() == 2)) + internal::merge_triangulations<2,2>(triangulation_1, triangulation_2, triangulation); + else if ((triangulation_1.get_dim() == 2) && (triangulation_1.get_spacedim() == 3)) + internal::merge_triangulations<2,3>(triangulation_1, triangulation_2, triangulation); + else + internal::merge_triangulations<3,3>(triangulation_1, triangulation_2, triangulation); + } + - if (triangulation_1.get_dim() == 2) - internal::merge_triangulations<2>(triangulation_1, triangulation_2, triangulation); + + void TriangulationWrapper::refine_global(const unsigned int n) + { + if ((dim == 2) && (spacedim == 2)) + { + dealii::Triangulation<2,2> *tria = + static_cast*>(triangulation); + tria->refine_global(n); + } + else if ((dim == 2) && (spacedim == 3)) + { + dealii::Triangulation<2,3> *tria = + static_cast*>(triangulation); + tria->refine_global(n); + } else - internal::merge_triangulations<3>(triangulation_1, triangulation_2, triangulation); + { + dealii::Triangulation<3,3> *tria = + static_cast*>(triangulation); + tria->refine_global(n); + } } @@ -343,17 +404,26 @@ namespace PyDealII std::ofstream ofs(filename); boost::archive::text_oarchive oa(ofs); - if (dim == 2) + if ((dim == 2) && (spacedim == 2)) + { + dealii::Triangulation<2,2> *tria = + static_cast*>(triangulation); + + oa << *tria; + } + else if ((dim == 2) && (spacedim == 3)) + { { - dealii::Triangulation<2> *tria = - static_cast*>(triangulation); + dealii::Triangulation<2,3> *tria = + static_cast*>(triangulation); oa << *tria; } + } else { - dealii::Triangulation<3> *tria = - static_cast*>(triangulation); + dealii::Triangulation<3,3> *tria = + static_cast*>(triangulation); oa << *tria; } @@ -366,17 +436,24 @@ namespace PyDealII std::ifstream ifs(filename); boost::archive::text_iarchive ia(ifs); - if (dim == 2) + if ((dim == 2) && (spacedim == 2)) { - dealii::Triangulation<2> *tria = - static_cast*>(triangulation); + dealii::Triangulation<2,2> *tria = + static_cast*>(triangulation); + + ia >> *tria; + } + else if ((dim == 2) && (spacedim == 3)) + { + dealii::Triangulation<2,3> *tria = + static_cast*>(triangulation); ia >> *tria; } else { dealii::Triangulation<3> *tria = - static_cast*>(triangulation); + static_cast*>(triangulation); ia >> *tria; } @@ -384,19 +461,38 @@ namespace PyDealII - void TriangulationWrapper::refine_global(const unsigned int n) + void TriangulationWrapper::setup(const std::string &dimension, + const std::string &spacedimension) { - if (dim == 2) + if ((dimension.compare("2D") == 0) || (dimension.compare("2d") == 0)) { - dealii::Triangulation<2> *tria = - static_cast*>(triangulation); - tria->refine_global(n); + dim = 2; + + if ((spacedimension.compare("2D") == 0) || + (spacedimension.compare("2d") == 0)) + { + spacedim = 2; + triangulation = new dealii::Triangulation<2,2>(); + } + else if ((spacedimension.compare("3D") == 0) || + (spacedimension.compare("3d") == 0)) + { + spacedim = 3; + triangulation = new dealii::Triangulation<2,3>(); + } + else + AssertThrow(false, dealii::ExcMessage("Spacedimension needs to be 2D or 3D.")); } - else + else if ((dimension.compare("3D") == 0) || (dimension.compare("3d") == 0)) { - dealii::Triangulation<3> *tria = - static_cast*>(triangulation); - tria->refine_global(n); + if ((spacedimension.compare("3D") != 0) && + (spacedimension.compare("3d") != 0)) + AssertThrow(false, dealii::ExcMessage("Spacedimension needs to be 3D.")); + dim = 3; + spacedim = 3; + triangulation = new dealii::Triangulation<3,3>(); } + else + AssertThrow(false, dealii::ExcMessage("Dimension needs to be 2D or 3D.")); } } diff --git a/tests/python/triangulation_wrapper.py b/tests/python/triangulation_wrapper.py index 42cc7d4216..04eabbe905 100644 --- a/tests/python/triangulation_wrapper.py +++ b/tests/python/triangulation_wrapper.py @@ -20,16 +20,16 @@ from pydealii.PyDealII import * class TestTriangulationWrapper(unittest.TestCase): def setUp(self): - self.dim = ['2D', '3D'] + self.dim = [['2D', '2D'], ['2D', '3D'], ['3D', '3D']] def build_hyper_cube_triangulation(self, dim): - triangulation_1 = Triangulation(dim) + triangulation_1 = Triangulation(dim[0], dim[1]) triangulation_1.generate_hyper_cube() return triangulation_1 def build_hyper_rectangle_triangulation(self, dim): - triangulation_2 = Triangulation(dim) - if (dim == '2D'): + triangulation_2 = Triangulation(dim[0], dim[1]) + if (dim[0] == '2D'): point_1 = Point([0., 0.]) point_2 = Point([1., 1.]) else: @@ -46,8 +46,8 @@ class TestTriangulationWrapper(unittest.TestCase): def test_simplex(self): for dim in self.dim: - triangulation = Triangulation(dim) - if (dim == '2D'): + triangulation = Triangulation(dim[0]) + if (dim[0] == '2D'): point_1 = Point([0., 0.]) point_2 = Point([1., 0.]) point_3 = Point([1., 1.]) @@ -64,11 +64,11 @@ class TestTriangulationWrapper(unittest.TestCase): def test_subdivided_hyper_cube(self): for dim in self.dim: - triangulation = Triangulation(dim) + triangulation = Triangulation(dim[0], dim[1]) repetitions = 2 triangulation.generate_subdivided_hyper_cube(repetitions) n_cells = triangulation.n_active_cells() - if (dim == '2D'): + if (dim[0] == '2D'): self.assertEqual(n_cells, 4) else: self.assertEqual(n_cells, 8) @@ -81,8 +81,8 @@ class TestTriangulationWrapper(unittest.TestCase): def test_subdivided_hyper_rectangle(self): for dim in self.dim: - triangulation = Triangulation(dim) - if (dim == '2D'): + triangulation = Triangulation(dim[0], dim[1]) + if (dim[0] == '2D'): repetitions = [6, 4] point_1 = Point([0., 0.]) point_2 = Point([1., 1.]) @@ -98,8 +98,8 @@ class TestTriangulationWrapper(unittest.TestCase): def test_hyper_ball(self): for dim in self.dim: - triangulation = Triangulation(dim) - if (dim == '2D'): + triangulation = Triangulation(dim[0]) + if (dim[0] == '2D'): center = Point([0., 0.]) n_cells_ref = 5 else: @@ -113,11 +113,11 @@ class TestTriangulationWrapper(unittest.TestCase): for dim in self.dim: triangulation_1 = self.build_hyper_cube_triangulation(dim) triangulation_2 = self.build_hyper_rectangle_triangulation(dim) - if (dim == '2D'): + if (dim[1] == '2D'): triangulation_2.shift([1., 0.]) else: triangulation_2.shift([1., 0., 0.]) - triangulation = Triangulation(dim) + triangulation = Triangulation(dim[0], dim[1]) triangulation.merge_triangulations(triangulation_1, triangulation_2) n_cells = triangulation.n_active_cells() @@ -128,7 +128,7 @@ class TestTriangulationWrapper(unittest.TestCase): triangulation = self.build_hyper_cube_triangulation(dim) triangulation.refine_global(2) n_cells = triangulation.n_active_cells() - if (dim == '2D'): + if (dim[0] == '2D'): self.assertEqual(n_cells, 16) else: self.assertEqual(n_cells, 64) @@ -138,10 +138,10 @@ class TestTriangulationWrapper(unittest.TestCase): triangulation_1 = self.build_hyper_cube_triangulation(dim) triangulation_1.refine_global(1) triangulation_1.save('mesh.output') - triangulation_2 = Triangulation(dim) + triangulation_2 = Triangulation(dim[0], dim[1]) triangulation_2.load('mesh.output') n_cells = triangulation_2.n_active_cells() - if (dim == '2D'): + if (dim[0] == '2D'): self.assertEqual(n_cells, 4) else: self.assertEqual(n_cells, 8) -- 2.39.5