From: Alexander Grayver Date: Mon, 18 Nov 2019 14:09:07 +0000 (+0100) Subject: Add mapping wrappers X-Git-Tag: v9.2.0-rc1~844^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9055%2Fhead;p=dealii.git Add mapping wrappers --- diff --git a/contrib/python-bindings/include/cell_accessor_wrapper.h b/contrib/python-bindings/include/cell_accessor_wrapper.h index 159349691a..2802b9458c 100644 --- a/contrib/python-bindings/include/cell_accessor_wrapper.h +++ b/contrib/python-bindings/include/cell_accessor_wrapper.h @@ -190,6 +190,8 @@ namespace python * Pointer that can be casted to the underlying CellAccessor object. */ void *cell_accessor; + + friend class MappingQGenericWrapper; }; diff --git a/contrib/python-bindings/include/mapping_wrapper.h b/contrib/python-bindings/include/mapping_wrapper.h new file mode 100644 index 0000000000..b80752ce14 --- /dev/null +++ b/contrib/python-bindings/include/mapping_wrapper.h @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 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_mapping_wrapper_h +#define dealii_mapping_wrapper_h + +#include + +#include + +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace python +{ + class MappingQGenericWrapper + { + public: + /** + * Copy constructor. + */ + MappingQGenericWrapper(const MappingQGenericWrapper &other); + + /** + * Constructor. + */ + MappingQGenericWrapper(const int dim, const int spacedim, const int degree); + + /** + * Destructor. + */ + ~MappingQGenericWrapper(); + + /** + * Map the point p on the unit cell to the corresponding point on the real + * cell. + */ + PointWrapper + transform_unit_to_real_cell(CellAccessorWrapper &cell, PointWrapper &point); + + /** + * Map the point p on the real cell to the corresponding point on the unit + * cell. + */ + PointWrapper + transform_real_to_unit_cell(CellAccessorWrapper &cell, PointWrapper &point); + + private: + /** + * Dimension of the underlying Mapping object. + */ + int dim; + + /** + * Space dimension of the underlying Mapping object. + */ + int spacedim; + + /** + * Degree of the underlying mapping. + */ + int degree; + + /** + * Pointer that can be casted to the underlying Mapping object. + */ + void *mapping_ptr; + }; + +} // namespace python + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/contrib/python-bindings/include/point_wrapper.h b/contrib/python-bindings/include/point_wrapper.h index 57d24db69a..f2b8ae14ac 100644 --- a/contrib/python-bindings/include/point_wrapper.h +++ b/contrib/python-bindings/include/point_wrapper.h @@ -230,6 +230,8 @@ namespace python * Pointer to the underlying Point object. */ void *point; + + friend class MappingQGenericWrapper; }; diff --git a/contrib/python-bindings/source/CMakeLists.txt b/contrib/python-bindings/source/CMakeLists.txt index 9a6ac24d88..376dbe35b0 100644 --- a/contrib/python-bindings/source/CMakeLists.txt +++ b/contrib/python-bindings/source/CMakeLists.txt @@ -94,10 +94,12 @@ SET(_src export_cell_accessor.cc export_point.cc export_triangulation.cc + export_mapping.cc cell_accessor_wrapper.cc tria_accessor_wrapper.cc point_wrapper.cc triangulation_wrapper.cc + mapping_wrapper.cc ) FOREACH(_build ${DEAL_II_BUILD_TYPES}) diff --git a/contrib/python-bindings/source/cell_accessor_wrapper.cc b/contrib/python-bindings/source/cell_accessor_wrapper.cc index 3a5ab80b50..958bb52ee9 100644 --- a/contrib/python-bindings/source/cell_accessor_wrapper.cc +++ b/contrib/python-bindings/source/cell_accessor_wrapper.cc @@ -167,7 +167,7 @@ namespace python static_cast *>(cell_accessor); Point barycenter = cell->barycenter(); boost::python::list barycenter_list; - for (int i = 0; i < dim; ++i) + for (int i = 0; i < spacedim; ++i) barycenter_list.append(barycenter[i]); return PointWrapper(barycenter_list); diff --git a/contrib/python-bindings/source/export_mapping.cc b/contrib/python-bindings/source/export_mapping.cc new file mode 100644 index 0000000000..e93a007063 --- /dev/null +++ b/contrib/python-bindings/source/export_mapping.cc @@ -0,0 +1,52 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace python +{ + const char transform_unit_to_real_cell_docstring[] = + " Map the point p on the unit cell to the corresponding point \n" + " on the real cell. \n"; + + const char transform_real_to_unit_cell_docstring[] = + " Map the point p on the real cell to the corresponding point \n" + " on the unit cell. \n"; + + + + void + export_mapping() + { + boost::python::class_( + "MappingQGeneric", + boost::python::init( + boost::python::args("dim", "spacedim", "degree"))) + .def("transform_real_to_unit_cell", + &MappingQGenericWrapper::transform_real_to_unit_cell, + transform_real_to_unit_cell_docstring, + boost::python::args("self", "cell_wrapper", "point_wrapper")) + .def("transform_unit_to_real_cell", + &MappingQGenericWrapper::transform_unit_to_real_cell, + transform_unit_to_real_cell_docstring, + boost::python::args("self", "cell_wrapper", "point_wrapper")); + } +} // namespace python + +DEAL_II_NAMESPACE_CLOSE diff --git a/contrib/python-bindings/source/mapping_wrapper.cc b/contrib/python-bindings/source/mapping_wrapper.cc new file mode 100644 index 0000000000..0c7ad75c67 --- /dev/null +++ b/contrib/python-bindings/source/mapping_wrapper.cc @@ -0,0 +1,227 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 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. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace python +{ + namespace internal + { + template + PointWrapper + transform_unit_to_real_cell(void *mapping_ptr, + void *cell_accessor_ptr, + void *point_ptr) + { + const MappingQGeneric *mapping = + static_cast *>(mapping_ptr); + + const CellAccessor *cell_accessor = + static_cast *>(cell_accessor_ptr); + + const Point *point = static_cast *>(point_ptr); + + typename Triangulation::active_cell_iterator cell( + &cell_accessor->get_triangulation(), + cell_accessor->level(), + cell_accessor->index()); + + Point p_real = + mapping->transform_unit_to_real_cell(cell, *point); + + boost::python::list coord_list; + for (int i = 0; i < spacedim; ++i) + coord_list.append(p_real[i]); + + return PointWrapper(coord_list); + } + + + + template + PointWrapper + transform_real_to_unit_cell(void *mapping_ptr, + void *cell_accessor_ptr, + void *point_ptr) + { + const MappingQGeneric *mapping = + static_cast *>(mapping_ptr); + + const CellAccessor *cell_accessor = + static_cast *>(cell_accessor_ptr); + + const Point *point = + static_cast *>(point_ptr); + + typename Triangulation::active_cell_iterator cell( + &cell_accessor->get_triangulation(), + cell_accessor->level(), + cell_accessor->index()); + + Point p_real = mapping->transform_real_to_unit_cell(cell, *point); + + boost::python::list coord_list; + for (int i = 0; i < dim; ++i) + coord_list.append(p_real[i]); + + return PointWrapper(coord_list); + } + } // namespace internal + + + + MappingQGenericWrapper::MappingQGenericWrapper(const int dim, + const int spacedim, + const int degree) + : dim(dim) + , spacedim(spacedim) + , degree(degree) + { + if ((dim == 2) && (spacedim == 2)) + { + mapping_ptr = new MappingQGeneric<2, 2>(degree); + } + else if ((dim == 2) && (spacedim == 3)) + { + mapping_ptr = new MappingQGeneric<2, 3>(degree); + } + else if ((dim == 3) && (spacedim == 3)) + { + mapping_ptr = new MappingQGeneric<3, 3>(degree); + } + else + AssertThrow(false, ExcMessage("Wrong dim-spacedim combination.")); + } + + + + MappingQGenericWrapper::MappingQGenericWrapper( + const MappingQGenericWrapper &other) + { + dim = other.dim; + spacedim = other.spacedim; + degree = other.degree; + + AssertThrow(other.mapping_ptr != nullptr, + ExcMessage("Underlying mapping does not exist.")); + + if ((dim == 2) && (spacedim == 2)) + { + mapping_ptr = new MappingQGeneric<2, 2>(other.degree); + } + else if ((dim == 2) && (spacedim == 3)) + { + mapping_ptr = new MappingQGeneric<2, 3>(other.degree); + } + else if ((dim == 3) && (spacedim == 3)) + { + mapping_ptr = new MappingQGeneric<3, 3>(other.degree); + } + else + AssertThrow(false, ExcMessage("Wrong dim-spacedim combination.")); + } + + + + MappingQGenericWrapper::~MappingQGenericWrapper() + { + if (dim != -1) + { + if ((dim == 2) && (spacedim == 2)) + { + // We cannot call delete on a void pointer so cast the void pointer + // back first. + MappingQGeneric<2, 2> *tmp = + static_cast *>(mapping_ptr); + delete tmp; + } + else if ((dim == 2) && (spacedim == 3)) + { + MappingQGeneric<2, 3> *tmp = + static_cast *>(mapping_ptr); + delete tmp; + } + else + { + MappingQGeneric<3, 3> *tmp = + static_cast *>(mapping_ptr); + delete tmp; + } + + dim = -1; + spacedim = -1; + degree = -1; + mapping_ptr = nullptr; + } + } + + + + PointWrapper + MappingQGenericWrapper::transform_unit_to_real_cell(CellAccessorWrapper &cell, + PointWrapper & p) + { + AssertThrow( + dim == p.get_dim(), + ExcMessage( + "Dimension of the point is not equal to the dimension of the mapping.")); + + if ((dim == 2) && (spacedim == 2)) + return internal::transform_unit_to_real_cell<2, 2>(mapping_ptr, + cell.cell_accessor, + p.point); + else if ((dim == 2) && (spacedim == 3)) + return internal::transform_unit_to_real_cell<2, 3>(mapping_ptr, + cell.cell_accessor, + p.point); + else + return internal::transform_unit_to_real_cell<3, 3>(mapping_ptr, + cell.cell_accessor, + p.point); + } + + + + PointWrapper + MappingQGenericWrapper::transform_real_to_unit_cell(CellAccessorWrapper &cell, + PointWrapper & p) + { + AssertThrow( + spacedim == p.get_dim(), + ExcMessage( + "Dimension of the point is not equal to the space dimension of the mapping.")); + + if ((dim == 2) && (spacedim == 2)) + return internal::transform_real_to_unit_cell<2, 2>(mapping_ptr, + cell.cell_accessor, + p.point); + else if ((dim == 2) && (spacedim == 3)) + return internal::transform_real_to_unit_cell<2, 3>(mapping_ptr, + cell.cell_accessor, + p.point); + else + return internal::transform_real_to_unit_cell<3, 3>(mapping_ptr, + cell.cell_accessor, + p.point); + } + +} // namespace python + +DEAL_II_NAMESPACE_CLOSE diff --git a/contrib/python-bindings/source/tria_accessor_wrapper.cc b/contrib/python-bindings/source/tria_accessor_wrapper.cc index 49188e9512..7dfc296937 100644 --- a/contrib/python-bindings/source/tria_accessor_wrapper.cc +++ b/contrib/python-bindings/source/tria_accessor_wrapper.cc @@ -34,7 +34,7 @@ namespace python tria_accessor); Point barycenter = accessor->barycenter(); boost::python::list barycenter_list; - for (int i = 0; i < dim; ++i) + for (int i = 0; i < spacedim; ++i) barycenter_list.append(barycenter[i]); return PointWrapper(barycenter_list); diff --git a/contrib/python-bindings/source/wrappers.cc b/contrib/python-bindings/source/wrappers.cc index d19257597c..5d26719cd0 100644 --- a/contrib/python-bindings/source/wrappers.cc +++ b/contrib/python-bindings/source/wrappers.cc @@ -31,6 +31,8 @@ namespace python export_point(); void export_triangulation(); + void + export_mapping(); } // namespace python DEAL_II_NAMESPACE_CLOSE @@ -63,6 +65,7 @@ BOOST_PYTHON_MODULE(Debug) dealii::python::export_cell_accessor(); dealii::python::export_point(); dealii::python::export_triangulation(); + dealii::python::export_mapping(); } #else @@ -80,6 +83,7 @@ BOOST_PYTHON_MODULE(Release) dealii::python::export_cell_accessor(); dealii::python::export_point(); dealii::python::export_triangulation(); + dealii::python::export_mapping(); } #endif diff --git a/contrib/python-bindings/tests/mapping_wrapper.py b/contrib/python-bindings/tests/mapping_wrapper.py new file mode 100644 index 0000000000..ef18ced7e5 --- /dev/null +++ b/contrib/python-bindings/tests/mapping_wrapper.py @@ -0,0 +1,60 @@ +# --------------------------------------------------------------------- +# +# 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.md at +# the top level directory of deal.II. +# +# --------------------------------------------------------------------- + +import os +import sys +module_path = os.path.abspath('/home/agrayver/lib/dealii/build_gofem/lib/python3.7/site-packages') +if module_path not in sys.path: + sys.path.append(module_path) + +import unittest +from PyDealII.Debug import * + +class TestMappingWrapperCube(unittest.TestCase): + + def setUp(self): + self.triangulation = Triangulation('2D') + self.triangulation.generate_hyper_cube() + self.triangulation.refine_global(1) + + def test_mapping(self): + mapping = MappingQGeneric(dim = 2, spacedim = 2, degree = 1) + p_unit = Point([0.5, 0.5]) + + for cell in self.triangulation.active_cells(): + p_real = mapping.transform_unit_to_real_cell(cell, p_unit) + p_unit_back = mapping.transform_real_to_unit_cell(cell, p_real) + + self.assertTrue(p_unit.distance(p_unit_back) < 1e-10) + +class TestMappingWrapperSphere(unittest.TestCase): + + def setUp(self): + self.triangulation = Triangulation('2D', '3D') + p_center = Point([0., 0., 0.]) + self.triangulation.generate_hyper_sphere(p_center) + + def test_mapping(self): + mapping = MappingQGeneric(dim = 2, spacedim = 3, degree = 4) + p_unit = Point([0.5, 0.5]) + + for cell in self.triangulation.active_cells(): + p_real = mapping.transform_unit_to_real_cell(cell, p_unit) + p_unit_back = mapping.transform_real_to_unit_cell(cell, p_real) + + self.assertTrue(p_unit.distance(p_unit_back) < 1e-10) + +if __name__ == '__main__': + unittest.main() diff --git a/doc/news/changes/minor/20191118Grayver b/doc/news/changes/minor/20191118Grayver new file mode 100644 index 0000000000..5b89534f77 --- /dev/null +++ b/doc/news/changes/minor/20191118Grayver @@ -0,0 +1,3 @@ +Added: New python wrappers for the MappingQGeneric class +
+(Alexander Grayver, 2019/11/18)