From: Martin Kronbichler Date: Wed, 1 Mar 2017 17:35:14 +0000 (+0100) Subject: Fix manifold code by introducing a new flag X-Git-Tag: v8.5.0-rc1~17^2~7 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6fdbe51d91162146d6ea6c4bf9ef3f21b7768dd5;p=dealii.git Fix manifold code by introducing a new flag --- diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 7787f93fa9..2694e8521d 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -69,6 +69,49 @@ DEAL_II_NAMESPACE_OPEN * cells; again this class is then equivalent to using MappingQGeneric right * away. * + *

Behavior along curved boundaries and with different manifolds

+ * + * As described above, one often only knows a manifold description of a + * surface but not the interior of the computational domain. Thus, a + * StraightBoundary object will be assigned to the interior entities that + * describes a usual planar coordinate system where the additional points of + * this class are placed exactly according to a bi-/trilinear mapping. When + * combined with a non-flat manifold on the boundary, for example a circle, + * the two manifold descriptions are in general incompatible. For example, a + * StraightBoundary would put an interior point located at some small distance + * epsilon away from the boundary along a flat line and thus in general + * outside the concave part of a circle. If the polynomial degree of + * MappingQGeneric is sufficiently high, the transformation from the reference + * cell to such a cell would in general contain inverted regions close to the + * boundary. + * + * In order to avoid this situation, this class applies a smoothing on cells + * adjacent to the boundary by using so-called Laplace smoothing. In the + * algorithm computing new points, all the entities of the cells are passed + * through hierarchically, starting from the lines to the quads and finally + * hexes. The elements higher up in the hierarchy that sit on a + * StraightBoundary will then get their points interpolated from all the + * surrounding points and not just the corner points. If only a line is + * associated to a curved boundary but the adjacent quad is on a flat + * manifold, the points inside the quad will be computed according to the + * deformed line and thus always result in a well-defined transformation. If + * this smoothing is undesired, the optional argument @p + * smooth_support_points_on_flat_manifold can be used to disable the + * smoothing, placing the additional points strictly according to the + * manifold. + * + * While this smoothing approach works well for low and medium convergence + * orders up to approximately three to four, this mechanism has an inherent + * shortcoming because it switches from a curved manifold to a flat manifold + * within one layer of elements. This will cause the Jacobian transformation + * to have jumps between the first and second element layer that can reduce + * the order of convergence. For example, the convergence rates for solving + * the Laplacian on a circle where only the boundary is deformed and the above + * mesh smoothing algorithm is applied will typically not exceed 3.5, even for + * fourth or fifth degree polynomials. In such a case, the curved manifold + * needs to be switched to a flat manifold in a smooth way that does not + * depend on the mesh size and eventuell covers a whole layer of cells. + * * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001, * Wolfgang Bangerth, 2015 */ @@ -92,9 +135,14 @@ public: * The value of @p use_mapping_q_on_all_cells is ignored if @p dim is not * equal to @p spacedim, i.e., if we are considering meshes on surfaces * embedded into higher dimensional spaces. + * + * The optional parameter @p smooth_support_points_on_flat_manifold controls + * whether smoothing adjancent to StraightBoundary objects according to the + * general class description should be enabled (default) or not. */ MappingQ (const unsigned int polynomial_degree, - const bool use_mapping_q_on_all_cells = false); + const bool use_mapping_q_on_all_cells = false, + const bool smooth_support_points_on_flat_manifold = true); /** * Copy constructor. diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index e9ec6483fa..1832e37786 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -81,6 +81,48 @@ template class MappingQ; * a manifold object to all cells that can then also be used by the mapping * classes for higher order mappings. * + *

Behavior along curved boundaries and with different manifolds

+ * + * As described above, one often only knows a manifold description of a + * surface but not the interior of the computational domain. Thus, a + * StraightBoundary object will be assigned to the interior entities that + * describes a usual planar coordinate system where the additional points of + * this class are placed exactly according to a bi-/trilinear mapping. When + * combined with a non-flat manifold on the boundary, for example a circle, + * the two manifold descriptions are in general incompatible. For example, a + * StraightBoundary would put an interior point located at some small distance + * epsilon away from the boundary along a flat line and thus in general + * outside the concave part of a circle. If the polynomial degree of + * MappingQGeneric is sufficiently high, the transformation from the reference + * cell to such a cell would in general contain inverted regions close to the + * boundary. + * + * In order to avoid this situation, this class applies a smoothing on cells + * adjacent to the boundary by using so-called Laplace smoothing. In the + * algorithm computing new points, all the entities of the cells are passed + * through hierarchically, starting from the lines to the quads and finally + * hexes. The elements higher up in the hierarchy that sit on a + * StraightBoundary will then get their points interpolated from all the + * surrounding points and not just the corner points. If only a line is + * associated to a curved boundary but the adjacent quad is on a flat + * manifold, the points inside the quad will be computed according to the + * deformed line and thus always result in a well-defined transformation. If + * this smoothing is undesired, the optional argument @p + * smooth_support_points_on_flat_manifold can be used to disable the + * smoothing, placing the additional points strictly according to the + * manifold. + * + * While this smoothing approach works well for low and medium convergence + * orders up to approximately three to four, this mechanism has an inherent + * shortcoming because it switches from a curved manifold to a flat manifold + * within one layer of elements. This will cause the Jacobian transformation + * to have jumps between the first and second element layer that can reduce + * the order of convergence. For example, the convergence rates for solving + * the Laplacian on a circle where only the boundary is deformed and the above + * mesh smoothing algorithm is applied will typically not exceed 3.5, even for + * fourth or fifth degree polynomials. In such a case, the curved manifold + * needs to be switched to a flat manifold in a smooth way that does not + * depend on the mesh size and eventuell covers a whole layer of cells. * * @author Wolfgang Bangerth, 2015 */ @@ -92,8 +134,13 @@ public: * Constructor. @p polynomial_degree denotes the polynomial degree of the * polynomials that are used to map cells from the reference to the real * cell. + * + * The optional parameter @p smooth_support_points_on_flat_manifold controls + * whether smoothing adjancent to StraightBoundary objects according to the + * general class description should be enabled (default) or not. */ - MappingQGeneric (const unsigned int polynomial_degree); + MappingQGeneric (const unsigned int polynomial_degree, + const bool smooth_support_points_on_flat_manifold = true); /** * Copy constructor. @@ -505,6 +552,12 @@ protected: */ const unsigned int polynomial_degree; + /** + * Stores whether we want to smooth the placement of interior points on flat + * manifolds. + */ + const bool smooth_support_points_on_flat_manifold; + /* * The default line support points. These are used when computing * the location in real space of the support points on lines and diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 12d4ea6e03..87ed327bb8 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -54,7 +54,8 @@ MappingQ::InternalData::memory_consumption () const template MappingQ::MappingQ (const unsigned int degree, - const bool use_mapping_q_on_all_cells) + const bool use_mapping_q_on_all_cells, + const bool smooth_support_points_on_flat_manifold) : polynomial_degree (degree), @@ -76,7 +77,7 @@ MappingQ::MappingQ (const unsigned int degree, // created via the shared_ptr objects qp_mapping (this->polynomial_degree>1 ? - std_cxx11::shared_ptr >(new MappingQGeneric(degree)) + std_cxx11::shared_ptr >(new MappingQGeneric(degree, smooth_support_points_on_flat_manifold)) : q1_mapping) {} diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index b9c6c65dc5..27a1fbb99f 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -1206,9 +1206,11 @@ namespace template -MappingQGeneric::MappingQGeneric (const unsigned int p) +MappingQGeneric::MappingQGeneric (const unsigned int p, + const bool smooth_support_points_on_flat_manifold) : polynomial_degree(p), + smooth_support_points_on_flat_manifold(smooth_support_points_on_flat_manifold), line_support_points(this->polynomial_degree+1), fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : 0), support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), @@ -1224,6 +1226,7 @@ template MappingQGeneric::MappingQGeneric (const MappingQGeneric &mapping) : polynomial_degree(mapping.polynomial_degree), + smooth_support_points_on_flat_manifold(mapping.smooth_support_points_on_flat_manifold), line_support_points(mapping.line_support_points), fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : 0), support_point_weights_on_quad (mapping.support_point_weights_on_quad), @@ -3812,9 +3815,13 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, // On a quad, we have to check whether the manifold should determine the // point distribution or rather a weighted sum should be created. This - // is the same logic as in the compute_mapping_support_points - // function below. - if (dynamic_cast *>(&face->get_manifold()) == NULL) + // is the same logic as in the compute_mapping_support_points function + // below: if we have a StraightBoundary but no "real" boundary + // description doing something else, we should go for the smoothing of + // the support points through add_weighted_interior_points + if (smooth_support_points_on_flat_manifold == false || + std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") == + std::string::npos) { // ask the boundary/manifold object to return intermediate points on it get_intermediate_points_on_object(face->get_manifold(), line_support_points, @@ -3901,13 +3908,17 @@ compute_mapping_support_points(const typename Triangulation::cell_ // (outer) points add_line_support_points(cell, a); - // then get the support points on the quad if we are on a - // manifold, otherwise compute them from the points around it + // then get the support points on the quad. if we have a + // StraightBoundary but no "real" boundary description doing something + // else, we should go for the smoothing of the support points through + // add_weighted_interior_points. we need to go through typeid because + // curved boundaries are derived from StraightBoundary and thus a + // plain dynamic_cast will not suffice. if (dim != spacedim) add_quad_support_points(cell, a); - // TODO: use get_intermediate_points_on_object for the else case - // unconditionally as soon as the interface of Boundary is fixed - else if (dynamic_cast *>(&cell->get_manifold()) == NULL) + else if (smooth_support_points_on_flat_manifold == false || + std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") == + std::string::npos) { std::vector > quad_points (Utilities::fixed_power(polynomial_degree-1)); get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points); @@ -3924,9 +3935,15 @@ compute_mapping_support_points(const typename Triangulation::cell_ add_line_support_points (cell, a); add_quad_support_points (cell, a); - // then compute the interior points - // TODO: remove else case as soon as boundary is fixed - if (dynamic_cast *>(&cell->get_manifold()) == NULL) + // then compute the interior points. if we have a StraightBoundary but + // no "real" boundary description doing something else, we should go + // for the smoothing of the support points through + // add_weighted_interior_points. we need to go through typeid because + // curved boundaries are derived from StraightBoundary and thus a + // plain dynamic_cast will not suffice. + if (smooth_support_points_on_flat_manifold == false || + std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") == + std::string::npos) { std::vector > hex_points (Utilities::fixed_power(polynomial_degree-1)); get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points); diff --git a/tests/mappings/mapping_q_mixed_manifolds_01.cc~ b/tests/mappings/mapping_q_mixed_manifolds_01.cc~ deleted file mode 100644 index fa0840c718..0000000000 --- a/tests/mappings/mapping_q_mixed_manifolds_01.cc~ +++ /dev/null @@ -1,201 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 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 at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -// Shows the shape functions implemented and computes the area of cells. - -#include "../tests.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -const double D = 0.1; -const double R = D/2.0; -const double R_1 = 1.2*R; -const double R_2 = 1.7*R; -const double H = 0.41; -const double X_0 = 0.0; -const double X_1 = 0.3; -const double X_C = 0.5; // center -const double X_2 = 0.7; - -const double Y_0 = 0.0; -const double Y_C = 0.2; // center - -const unsigned int MANIFOLD_ID = 1; - - -void create_triangulation(Triangulation<2> &tria) -{ - AssertThrow(std::abs((X_2-X_1) - 2.0*(X_C-X_1))<1.0e-12, ExcMessage("Geometry parameters X_1,X_2,X_C invalid!")); - SphericalManifold<2> spherical_manifold(Point<2>(X_C,Y_C)); - - Triangulation<2> circle_1, circle_2, circle_tmp, middle, middle_tmp, middle_tmp2, tmp_3D; - std::vector ref_1(2, 2); - ref_1[1] = 2; - - // create middle part first as a hyper shell - const double outer_radius = (X_2-X_1)/2.0; - const unsigned int n_cells = 4; - GridGenerator::hyper_shell(middle, Point<2>(X_C, Y_C), R_2, outer_radius, n_cells, true); - middle.set_all_manifold_ids(MANIFOLD_ID); - middle.set_manifold(MANIFOLD_ID, spherical_manifold); - middle.refine_global(1); - - // two inner circles in order to refine towards the cylinder surface - const unsigned int n_cells_circle = 8; - GridGenerator::hyper_shell(circle_1, Point<2>(X_C, Y_C), R, R_1, n_cells_circle, true); - circle_1.set_all_manifold_ids(MANIFOLD_ID); - circle_1.set_manifold(MANIFOLD_ID,spherical_manifold); - - GridGenerator::hyper_shell(circle_2, Point<2>(X_C, Y_C), R_1, R_2, n_cells_circle, true); - circle_2.set_all_manifold_ids(MANIFOLD_ID); - circle_2.set_manifold(MANIFOLD_ID,spherical_manifold); - - // then move the vertices to the points where we want them to be to create a slightly asymmetric cube with a hole - for (Triangulation<2>::cell_iterator cell = middle.begin(); cell != middle.end(); ++cell) - { - for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v) - { - Point<2> &vertex = cell->vertex(v); - if (std::abs(vertex[0] - X_2) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10) - { - vertex = Point<2>(X_2, H/2.0); - } - else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) - { - vertex = Point<2>(X_2, H); - } - else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) - { - vertex = Point<2>(X_2, Y_0); - } - else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C +(X_2-X_1)/2.0)) < 1e-10) - { - vertex = Point<2>(X_C, H); - } - else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C-(X_2-X_1)/2.0)) < 1e-10) - { - vertex = Point<2>(X_C, Y_0); - } - else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) - { - vertex = Point<2>(X_1, H); - } - else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) - { - vertex = Point<2>(X_1, Y_0); - } - else if (std::abs(vertex[0] - X_1) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10) - { - vertex = Point<2>(X_1, H/2.0); - } - } - } - - // must copy the triangulation because we cannot merge triangulations with refinement... - GridGenerator::flatten_triangulation(middle, middle_tmp); - GridGenerator::merge_triangulations(circle_1,circle_2,circle_tmp); - GridGenerator::merge_triangulations(middle_tmp,circle_tmp,tria); - - // Set the cylinder boundary to 2, outflow to 1, the rest to 0. - //tria.set_all_manifold_ids(0); - for (Triangulation<2>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell) - { - if (Point<2>(X_C,Y_C).distance(cell->center())<= R_2) - cell->set_all_manifold_ids(MANIFOLD_ID); - } -} - -void create_triangulation(Triangulation<3> &tria) -{ - Triangulation<2> tria_2d; - create_triangulation(tria_2d); - GridGenerator::extrude_triangulation(tria_2d, 3, H, tria); - - // Set the cylinder boundary to 2, outflow to 1, the rest to 0. - tria.set_all_manifold_ids(0); - for (Triangulation<3>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell) - { - if (Point<3>(X_C,Y_C,cell->center()[2]).distance(cell->center())<= R_2) - cell->set_all_manifold_ids(MANIFOLD_ID); - } -} - - - -template -void test() -{ - Point center; - center[0] = X_C; - center[1] = Y_C; - Point direction; - direction[dim-1] = 1.; - - std_cxx11::shared_ptr > cylinder_manifold - (dim == 2 ? static_cast*>(new SphericalManifold(center)) : - static_cast*>(new CylindricalManifold(direction, center))); - Triangulation tria; - create_triangulation(tria); - tria.set_manifold(MANIFOLD_ID, *cylinder_manifold); - - FE_Nothing fe; - for (unsigned int degree = 1; degree < 7; ++degree) - { - MappingQGeneric mapping(degree); - QGauss quad(degree+1); - FEValues fe_values(mapping, fe, quad, update_JxW_values); - double sum = 0.; - for (typename Triangulation::active_cell_iterator cell=tria.begin_active(); - cell != tria.end(); ++cell) - { - fe_values.reinit(cell); - double local_sum = 0; - for (unsigned int q=0; q(); - test<3>(); -}