From: Alexander Grayver Date: Mon, 2 Dec 2019 14:47:22 +0000 (+0100) Subject: Add more methods to work with manifolds. Add python tests X-Git-Tag: v9.2.0-rc1~811^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9109%2Fhead;p=dealii.git Add more methods to work with manifolds. Add python tests --- diff --git a/contrib/python-bindings/include/cell_accessor_wrapper.h b/contrib/python-bindings/include/cell_accessor_wrapper.h index 10e368a9d3..fd2d2bb284 100644 --- a/contrib/python-bindings/include/cell_accessor_wrapper.h +++ b/contrib/python-bindings/include/cell_accessor_wrapper.h @@ -125,6 +125,13 @@ namespace python int get_manifold_id() const; + /** + * Do as set_manifold_id() but also set the manifold indicators of + * the objects that bound the current object. + */ + void + set_all_manifold_ids(const int manifold_id); + /** * Return the ith neighbor of a cell. If the neighbor does not exist, * i.e., if the ith face of the current object is at the boundary, @@ -154,6 +161,18 @@ namespace python boost::python::list faces() const; + /** + * Compute the dim-dimensional measure of the object. + * For a dim-dimensional cell in dim-dimensional space, + * this equals its volume. On the other hand, for a 2d + * cell in 3d space, or if the current object pointed to + * is a 2d face of a 3d cell in 3d space, then the function + * computes the area the object occupies. For a + * one-dimensional object, return its length. + */ + double + measure() const; + /** * Exception. */ diff --git a/contrib/python-bindings/include/manifold_wrapper.h b/contrib/python-bindings/include/manifold_wrapper.h index 54b848d34c..c1e6cf2113 100644 --- a/contrib/python-bindings/include/manifold_wrapper.h +++ b/contrib/python-bindings/include/manifold_wrapper.h @@ -68,6 +68,17 @@ namespace python void * get_manifold() const; + /** + * Return the dimension of the underlying object + */ + int + get_dim() const; + + /** + * Return the space dimension of the underlying object + */ + int + get_spacedim() const; private: /** diff --git a/contrib/python-bindings/include/tria_accessor_wrapper.h b/contrib/python-bindings/include/tria_accessor_wrapper.h index b34d711087..ebdd8aa152 100644 --- a/contrib/python-bindings/include/tria_accessor_wrapper.h +++ b/contrib/python-bindings/include/tria_accessor_wrapper.h @@ -102,6 +102,18 @@ namespace python bool at_boundary() const; + /** + * Compute the dim-dimensional measure of the object. + * For a dim-dimensional cell in dim-dimensional space, + * this equals its volume. On the other hand, for a 2d + * cell in 3d space, or if the current object pointed to + * is a 2d face of a 3d cell in 3d space, then the function + * computes the area the object occupies. For a + * one-dimensional object, return its length. + */ + double + measure() const; + /** * Exception. */ diff --git a/contrib/python-bindings/include/triangulation_wrapper.h b/contrib/python-bindings/include/triangulation_wrapper.h index a599e62ed5..6910143fc7 100644 --- a/contrib/python-bindings/include/triangulation_wrapper.h +++ b/contrib/python-bindings/include/triangulation_wrapper.h @@ -278,7 +278,8 @@ namespace python generate_hyper_shell(PointWrapper & center, const double inner_radius, const double outer_radius, - const unsigned n_cells = 0); + const unsigned n_cells = 0, + bool colorize = false); /** * Shift each vertex of the Triangulation by the given @p shift_list. @@ -316,11 +317,26 @@ namespace python flatten_triangulation(TriangulationWrapper &tria_out); /** - * + * Assign a manifold object to a certain part of the triangulation. + * If an object with manifold number is refined, this object + * is used to find the location of new vertices (see the results + * section of step-49 for a more in-depth discussion of this, with + * examples). It is also used for non-linear (i.e.: non-Q1) + * transformations of cells to the unit cell in shape function + * calculations. */ void set_manifold(const int number, ManifoldWrapper &manifold); + /** + * Reset those parts of the triangulation with the given manifold_number to + * use a FlatManifold object. This is the default state of a non-curved + * triangulation, and undoes assignment of a different Manifold object by + * the function Triangulation::set_manifold(). + */ + void + reset_manifold(const int number); + /** * Refine all the cells @p n times. */ diff --git a/contrib/python-bindings/source/cell_accessor_wrapper.cc b/contrib/python-bindings/source/cell_accessor_wrapper.cc index bff328a579..73991421da 100644 --- a/contrib/python-bindings/source/cell_accessor_wrapper.cc +++ b/contrib/python-bindings/source/cell_accessor_wrapper.cc @@ -251,6 +251,17 @@ namespace python + template + void + set_all_manifold_ids(const int manifold_id, void *cell_accessor) + { + const CellAccessor *cell = + static_cast *>(cell_accessor); + return cell->set_all_manifold_ids(manifold_id); + } + + + template bool at_boundary(const void *cell_accessor) @@ -304,6 +315,18 @@ namespace python return faces_list; } + + + + template + double + measure(const void *cell_accessor) + { + const CellAccessor *cell = + static_cast *>(cell_accessor); + + return cell->measure(); + } } // namespace internal @@ -558,6 +581,19 @@ namespace python + void + CellAccessorWrapper::set_all_manifold_ids(const int manifold_id) + { + if ((dim == 2) && (spacedim == 2)) + internal::set_all_manifold_ids<2, 2>(manifold_id, cell_accessor); + else if ((dim == 2) && (spacedim == 3)) + internal::set_all_manifold_ids<2, 3>(manifold_id, cell_accessor); + else + internal::set_all_manifold_ids<3, 3>(manifold_id, cell_accessor); + } + + + bool CellAccessorWrapper::at_boundary() const { @@ -610,6 +646,19 @@ namespace python return internal::faces<3, 3>(cell_accessor); } + + + double + CellAccessorWrapper::measure() const + { + if ((dim == 2) && (spacedim == 2)) + return internal::measure<2, 2>(cell_accessor); + else if ((dim == 2) && (spacedim == 3)) + return internal::measure<2, 3>(cell_accessor); + else + return internal::measure<3, 3>(cell_accessor); + } + } // namespace python DEAL_II_NAMESPACE_CLOSE diff --git a/contrib/python-bindings/source/export_cell_accessor.cc b/contrib/python-bindings/source/export_cell_accessor.cc index 9235bdcbdd..877da21b98 100644 --- a/contrib/python-bindings/source/export_cell_accessor.cc +++ b/contrib/python-bindings/source/export_cell_accessor.cc @@ -57,6 +57,12 @@ namespace python + const char set_all_manifold_ids_docstring[] = + "Do as set_manifold_id() but also set the manifold indicators of \n" + "the objects that bound the current object. \n"; + + + const char barycenter_docstring[] = "Return the barycenter of the current cell \n"; @@ -97,6 +103,11 @@ namespace python + const char measure_docstring[] = + " Compute the dim-dimensional measure of the object. \n"; + + + void export_cell_accessor() { @@ -119,6 +130,10 @@ namespace python &CellAccessorWrapper::get_manifold_id, &CellAccessorWrapper::set_manifold_id, manifold_id_docstring) + .def("set_all_manifold_ids", + &CellAccessorWrapper::set_all_manifold_ids, + set_all_manifold_ids_docstring, + boost::python::args("self", "number")) .def("barycenter", &CellAccessorWrapper::get_barycenter, barycenter_docstring, @@ -146,6 +161,10 @@ namespace python .def("faces", &CellAccessorWrapper::faces, faces_docstring, + boost::python::args("self")) + .def("measure", + &CellAccessorWrapper::measure, + measure_docstring, boost::python::args("self")); } } // namespace python diff --git a/contrib/python-bindings/source/export_tria_accessor.cc b/contrib/python-bindings/source/export_tria_accessor.cc index 237aa481c2..9ac2c0f66f 100644 --- a/contrib/python-bindings/source/export_tria_accessor.cc +++ b/contrib/python-bindings/source/export_tria_accessor.cc @@ -44,6 +44,9 @@ namespace python const char at_boundary_docstring[] = " Return whether the face is at the boundary \n"; + const char measure_docstring[] = + " Compute the dim-dimensional measure of the object. \n"; + void export_tria_accessor() { @@ -77,7 +80,11 @@ namespace python .def("set_all_boundary_ids", &TriaAccessorWrapper::set_all_boundary_ids, set_all_boundary_ids_docstring, - boost::python::args("self", "boundary_id")); + boost::python::args("self", "boundary_id")) + .def("measure", + &TriaAccessorWrapper::measure, + measure_docstring, + boost::python::args("self")); } } // namespace python diff --git a/contrib/python-bindings/source/export_triangulation.cc b/contrib/python-bindings/source/export_triangulation.cc index 8992688d91..5919562899 100644 --- a/contrib/python-bindings/source/export_triangulation.cc +++ b/contrib/python-bindings/source/export_triangulation.cc @@ -93,11 +93,10 @@ namespace python generate_half_hyper_ball, 1, 2) - BOOST_PYTHON_MEMBER_FUNCTION_OVERLOADS(generate_hyper_shell_overloads, generate_hyper_shell, 3, - 4) + 5) const char n_active_cells_docstring[] = "Return the number of active cells \n"; @@ -360,6 +359,12 @@ namespace python + const char reset_manifold_docstring[] = + "Reset those parts of the triangulation with the given manifold_number \n" + "to use a FlatManifold object. \n"; + + + void export_triangulation() { @@ -468,7 +473,8 @@ namespace python "center", "inner_radius", "outer_radius", - "n_cells"), + "n_cells", + "colorize"), generate_hyper_shell_docstring)) .def("shift", &TriangulationWrapper::shift, @@ -513,7 +519,11 @@ namespace python .def("set_manifold", &TriangulationWrapper::set_manifold, set_manifold_docstring, - boost::python::args("self", "number", "manifold")); + boost::python::args("self", "number", "manifold")) + .def("reset_manifold", + &TriangulationWrapper::reset_manifold, + reset_manifold_docstring, + boost::python::args("self", "number")); } } // namespace python diff --git a/contrib/python-bindings/source/manifold_wrapper.cc b/contrib/python-bindings/source/manifold_wrapper.cc index 239ba9dfec..14d4888967 100644 --- a/contrib/python-bindings/source/manifold_wrapper.cc +++ b/contrib/python-bindings/source/manifold_wrapper.cc @@ -215,12 +215,30 @@ namespace python } + void * ManifoldWrapper::get_manifold() const { return manifold_ptr; } + + + int + ManifoldWrapper::get_dim() const + { + return dim; + } + + + + int + ManifoldWrapper::get_spacedim() const + { + return spacedim; + } + + } // 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 4f5daf5c14..196daad426 100644 --- a/contrib/python-bindings/source/tria_accessor_wrapper.cc +++ b/contrib/python-bindings/source/tria_accessor_wrapper.cc @@ -140,6 +140,17 @@ namespace python tria_accessor); return accessor->at_boundary(); } + + + template + double + measure(const void *tria_accessor) + { + const TriaAccessor *accessor = + static_cast *>( + tria_accessor); + return accessor->measure(); + } } // namespace internal @@ -336,6 +347,19 @@ namespace python return internal::at_boundary<2, 3, 3>(tria_accessor); } + + + double + TriaAccessorWrapper::measure() const + { + if ((dim == 2) && (spacedim == 2) && (structdim == 1)) + return internal::measure<1, 2, 2>(tria_accessor); + else if ((dim == 2) && (spacedim == 3) && (structdim == 1)) + return internal::measure<1, 2, 3>(tria_accessor); + else + return internal::measure<2, 3, 3>(tria_accessor); + } + } // namespace python DEAL_II_NAMESPACE_CLOSE diff --git a/contrib/python-bindings/source/triangulation_wrapper.cc b/contrib/python-bindings/source/triangulation_wrapper.cc index 5457d2f00e..d029380b1e 100644 --- a/contrib/python-bindings/source/triangulation_wrapper.cc +++ b/contrib/python-bindings/source/triangulation_wrapper.cc @@ -380,6 +380,7 @@ namespace python const double inner_radius, const double outer_radius, const unsigned n_cells, + bool colorize, void * triangulation) { // Cast the PointWrapper object to Point @@ -390,7 +391,7 @@ namespace python static_cast *>(triangulation); tria->clear(); GridGenerator::hyper_shell( - *tria, center_point, inner_radius, outer_radius, n_cells); + *tria, center_point, inner_radius, outer_radius, n_cells, colorize); } @@ -1057,7 +1058,8 @@ namespace python TriangulationWrapper::generate_hyper_shell(PointWrapper & center, const double inner_radius, const double outer_radius, - const unsigned n_cells) + const unsigned n_cells, + bool colorize) { AssertThrow( dim == spacedim, @@ -1065,10 +1067,10 @@ namespace python "This function is only implemented for dim equal to spacedim.")); if (dim == 2) internal::generate_hyper_shell<2>( - center, inner_radius, outer_radius, n_cells, triangulation); + center, inner_radius, outer_radius, n_cells, colorize, triangulation); else internal::generate_hyper_shell<3>( - center, inner_radius, outer_radius, n_cells, triangulation); + center, inner_radius, outer_radius, n_cells, colorize, triangulation); } @@ -1354,6 +1356,15 @@ namespace python TriangulationWrapper::set_manifold(const int number, ManifoldWrapper &manifold) { + AssertThrow( + dim == manifold.get_dim(), + ExcMessage( + "The Triangulation and Manifold should have the same dimension.")); + AssertThrow( + spacedim == manifold.get_spacedim(), + ExcMessage( + "The Triangulation and Manifold should have the same space dimension.")); + if ((dim == 2) && (spacedim == 2)) { Triangulation<2, 2> *tria = @@ -1382,6 +1393,31 @@ namespace python + void + TriangulationWrapper::reset_manifold(const int number) + { + if ((dim == 2) && (spacedim == 2)) + { + Triangulation<2, 2> *tria = + static_cast *>(triangulation); + tria->reset_manifold(number); + } + else if ((dim == 2) && (spacedim == 3)) + { + Triangulation<2, 3> *tria = + static_cast *>(triangulation); + tria->reset_manifold(number); + } + else + { + Triangulation<3, 3> *tria = + static_cast *>(triangulation); + tria->reset_manifold(number); + } + } + + + void TriangulationWrapper::setup(const std::string &dimension, const std::string &spacedimension) diff --git a/contrib/python-bindings/tests/manifold_wrapper.py b/contrib/python-bindings/tests/manifold_wrapper.py new file mode 100644 index 0000000000..d3270ccca2 --- /dev/null +++ b/contrib/python-bindings/tests/manifold_wrapper.py @@ -0,0 +1,71 @@ +# --------------------------------------------------------------------- +# +# 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 math +import unittest +from PyDealII.Debug import * + +class TestManifoldWrapperShell(unittest.TestCase): + + def setUp(self): + self.triangulation = Triangulation('2D') + p_center = Point([0, 0]) + self.triangulation.generate_hyper_shell(center = p_center, inner_radius = 0.5, outer_radius = 1., n_cells = 0, colorize = True) + + self.manifold = Manifold(dim = 2, spacedim = 2) + self.manifold.create_polar(p_center) + + def test_manifold(self): + self.triangulation.set_manifold(0, self.manifold) + for cell in self.triangulation.active_cells(): + cell.manifold_id = 0 + + self.triangulation.refine_global(3) + + circumference = 0 + for cell in self.triangulation.active_cells(): + for face in cell.faces(): + if face.at_boundary() and face.boundary_id == 1: + circumference += face.measure() + + self.assertTrue(abs(circumference - 2*math.pi)/(2*math.pi) < 1e-2) + +class TestManifoldWrapperBall(unittest.TestCase): + + def setUp(self): + self.triangulation = Triangulation('3D') + p_center = Point([0., 0., 0.]) + self.triangulation.generate_hyper_ball(center = p_center, radius = 1.) + + self.manifold = Manifold(dim = 3, spacedim = 3) + self.manifold.create_spherical(p_center) + + def test_manifold(self): + self.triangulation.reset_manifold(number = 0) + self.triangulation.set_manifold(number = 0, manifold = self.manifold) + for cell in self.triangulation.active_cells(): + if cell.at_boundary(): + cell.manifold_id = 0 + + self.triangulation.refine_global(3) + + volume = 0 + for cell in self.triangulation.active_cells(): + volume += cell.measure() + + self.assertTrue(abs(volume - 4./3. * math.pi) / (4./3.*math.pi) < 2e-2) + +if __name__ == '__main__': + unittest.main() diff --git a/contrib/python-bindings/tests/mapping_wrapper.py b/contrib/python-bindings/tests/mapping_wrapper.py index ef18ced7e5..b264e7bac5 100644 --- a/contrib/python-bindings/tests/mapping_wrapper.py +++ b/contrib/python-bindings/tests/mapping_wrapper.py @@ -13,12 +13,6 @@ # # --------------------------------------------------------------------- -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 *