]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix manifold code by introducing a new flag
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 1 Mar 2017 17:35:14 +0000 (18:35 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 18 Mar 2017 17:27:34 +0000 (18:27 +0100)
include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_generic.h
source/fe/mapping_q.cc
source/fe/mapping_q_generic.cc
tests/mappings/mapping_q_mixed_manifolds_01.cc~ [deleted file]

index 7787f93fa912a9cb01b872806ddfa13e010e8c60..2694e8521d040c4fa60ec7749bdc4f5feb5247ea 100644 (file)
@@ -69,6 +69,49 @@ DEAL_II_NAMESPACE_OPEN
  * cells; again this class is then equivalent to using MappingQGeneric right
  * away.
  *
+ * <h4>Behavior along curved boundaries and with different manifolds</h4>
+ *
+ * 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.
index e9ec6483fa9b53d16b0812e2f700334d6af4d25a..1832e377865186f47378bbe442565bb79f061460 100644 (file)
@@ -81,6 +81,48 @@ template <int,int> class MappingQ;
  * a manifold object to all cells that can then also be used by the mapping
  * classes for higher order mappings.
  *
+ * <h4>Behavior along curved boundaries and with different manifolds</h4>
+ *
+ * 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
index 12d4ea6e03d272481823fea99afd256685aa1780..87ed327bb8d5251a7fac24b7e0fd9798924f80bb 100644 (file)
@@ -54,7 +54,8 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 
 template<int dim, int spacedim>
 MappingQ<dim,spacedim>::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<dim,spacedim>::MappingQ (const unsigned int degree,
   // created via the shared_ptr objects
   qp_mapping (this->polynomial_degree>1
               ?
-              std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree))
+              std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree, smooth_support_points_on_flat_manifold))
               :
               q1_mapping)
 {}
index b9c6c65dc5ef86b9b671a52dab1145755bde643e..27a1fbb99f0bf42d043f028e89be7301519015b4 100644 (file)
@@ -1206,9 +1206,11 @@ namespace
 
 
 template<int dim, int spacedim>
-MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
+MappingQGeneric<dim,spacedim>::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<dim>(this->polynomial_degree) : 0),
   support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
@@ -1224,6 +1226,7 @@ template<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &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<dim>(*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<const Boundary<3,3> *>(&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<dim,spacedim>::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<const Boundary<dim,spacedim> *>(&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<Point<spacedim> > quad_points (Utilities::fixed_power<dim>(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<dim,spacedim>::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<const Boundary<dim,spacedim> *>(&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<Point<spacedim> > hex_points (Utilities::fixed_power<dim>(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 (file)
index fa0840c..0000000
+++ /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 <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/fe/mapping_q_generic.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe_nothing.h>
-#include <deal.II/grid/manifold_lib.h>
-#include <vector>
-#include <fstream>
-#include <iomanip>
-#include <string>
-#include <sstream>
-
-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<unsigned int> 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 <int dim>
-void test()
-{
-  Point<dim> center;
-  center[0] = X_C;
-  center[1] = Y_C;
-  Point<dim> direction;
-  direction[dim-1] = 1.;
-
-  std_cxx11::shared_ptr<Manifold<dim> > cylinder_manifold
-  (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(center)) :
-   static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, center)));
-  Triangulation<dim> tria;
-  create_triangulation(tria);
-  tria.set_manifold(MANIFOLD_ID, *cylinder_manifold);
-
-  FE_Nothing<dim> fe;
-  for (unsigned int degree = 1; degree < 7; ++degree)
-    {
-      MappingQGeneric<dim> mapping(degree);
-      QGauss<dim> quad(degree+1);
-      FEValues<dim> fe_values(mapping, fe, quad, update_JxW_values);
-      double sum = 0.;
-      for (typename Triangulation<dim>::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<quad.size(); ++q)
-            local_sum += fe_values.JxW(q);
-          sum += local_sum;
-        }
-      const double reference = (dim==2 ? 1. : H) * (H*(X_2-X_1)-D*D*numbers::PI*0.25);
-      deallog << "Volume " << dim << "D mapping degree " << degree << ": "
-              << sum << " error: " << (sum-reference)/reference << std::endl;
-    }
-}
-
-
-int main()
-{
-  std::ofstream logfile ("output");
-  deallog.attach(logfile);
-  deallog.threshold_double(1.e-14);
-
-  test<2>();
-  test<3>();
-}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.