]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Do not reuse CellSimilarity for higher order MappingQGeneric maps. 3738/head
authorDavid Wells <wellsd2@rpi.edu>
Sun, 8 Jan 2017 01:26:59 +0000 (20:26 -0500)
committerDavid Wells <wellsd2@rpi.edu>
Sun, 8 Jan 2017 05:24:50 +0000 (00:24 -0500)
This commit fixes bug where two cells, both alike in vertices (but varying in
curvature) would sometimes be treated as translations of each-other by
MappingQGeneric. This lead to an artificial limit of second order convergence
since Jacobian weights (and other values) were not correctly updated.

doc/news/changes/minor/20170105DavidWells [new file with mode: 0644]
source/fe/mapping_q_generic.cc
tests/fe/fe_values_function_manifold.cc [new file with mode: 0644]
tests/fe/fe_values_function_manifold.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170105DavidWells b/doc/news/changes/minor/20170105DavidWells
new file mode 100644 (file)
index 0000000..297cc39
--- /dev/null
@@ -0,0 +1,3 @@
+Fixed: MappingQGeneric, when computing values on a given cell, sometimes reused information computed on a different cell whose faces were curved in a different way. MappingQGeneric now recomputes all values when its order is greater than one (that is, when it can represent curved cells).
+<br>
+(David Wells, 2017/01/05)
index ae48043d156c8df10f76b58a6c8931d6e7190e91..333dae05e70ea037e976832e63277edb69358bec 100644 (file)
@@ -2589,10 +2589,17 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   data.mapping_support_points = this->compute_mapping_support_points(cell);
   data.cell_of_current_support_points = cell;
 
+  // if the order of the mapping is greater than 1, then do not reuse any cell
+  // similarity information. This is necessary because the cell similarity
+  // value is computed with just cell vertices and does not take into account
+  // cell curvature.
+  const CellSimilarity::Similarity computed_cell_similarity =
+    (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none);
+
   internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
                                                   data,
                                                   output_data.quadrature_points);
-  internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
+  internal::maybe_update_Jacobians<dim,spacedim> (computed_cell_similarity,
                                                   QProjector<dim>::DataSetDescriptor::cell (),
                                                   data);
 
@@ -2612,7 +2619,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
               ExcDimensionMismatch(output_data.normal_vectors.size(), n_q_points));
 
 
-      if (cell_similarity != CellSimilarity::translation)
+      if (computed_cell_similarity != CellSimilarity::translation)
         for (unsigned int point=0; point<n_q_points; ++point)
           {
 
@@ -2649,7 +2656,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                 output_data.JxW_values[point]
                   = sqrt(determinant(G)) * weights[point];
 
-                if (cell_similarity == CellSimilarity::inverted_translation)
+                if (computed_cell_similarity == CellSimilarity::inverted_translation)
                   {
                     // we only need to flip the normal
                     if (update_flags & update_normal_vectors)
@@ -2693,7 +2700,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   if (update_flags & update_jacobians)
     {
       AssertDimension (output_data.jacobians.size(), n_q_points);
-      if (cell_similarity != CellSimilarity::translation)
+      if (computed_cell_similarity != CellSimilarity::translation)
         for (unsigned int point=0; point<n_q_points; ++point)
           output_data.jacobians[point] = data.contravariant[point];
     }
@@ -2702,42 +2709,42 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   if (update_flags & update_inverse_jacobians)
     {
       AssertDimension (output_data.inverse_jacobians.size(), n_q_points);
-      if (cell_similarity != CellSimilarity::translation)
+      if (computed_cell_similarity != CellSimilarity::translation)
         for (unsigned int point=0; point<n_q_points; ++point)
           output_data.inverse_jacobians[point] = data.covariant[point].transpose();
     }
 
-  internal::maybe_update_jacobian_grads<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_grads<dim,spacedim> (computed_cell_similarity,
                                                        QProjector<dim>::DataSetDescriptor::cell (),
                                                        data,
                                                        output_data.jacobian_grads);
 
-  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_pushed_forward_grads<dim,spacedim> (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_pushed_forward_grads);
 
-  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_2nd_derivatives<dim,spacedim> (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_2nd_derivatives);
 
-  internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,spacedim> (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_pushed_forward_2nd_derivatives);
 
-  internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_3rd_derivatives);
 
-  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (cell_similarity,
+  internal::maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,spacedim> (computed_cell_similarity,
       QProjector<dim>::DataSetDescriptor::cell (),
       data,
       output_data.jacobian_pushed_forward_3rd_derivatives);
 
-  return cell_similarity;
+  return computed_cell_similarity;
 }
 
 
diff --git a/tests/fe/fe_values_function_manifold.cc b/tests/fe/fe_values_function_manifold.cc
new file mode 100644 (file)
index 0000000..f277963
--- /dev/null
@@ -0,0 +1,376 @@
+#include <deal.II/base/function.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "../tests.h"
+
+// Test that, for a single thread, CellSimilarity calculations return 'none'
+// when we are using curved cells. This is a regression test for a bug where,
+// for this particular manifold, some cells would be marked as translations
+// even though they were not (they varied in curvature of faces). This (small)
+// error causes the Laplace solver to converge at order 2 instead of order 5,
+// resulting in much higher L2 errors when the grid is refined.
+
+using namespace dealii;
+
+static const unsigned int fe_order = 4;
+static const dealii::types::boundary_id boundary_id = 0;
+static const dealii::types::manifold_id cubic_manifold_id = 1;
+static const double pi = M_PI;
+
+// ----------------------------------------------------------------------------
+// Manufactured solution and manufactured forcing
+// ----------------------------------------------------------------------------
+
+template <int dim>
+class HardManufacturedSolution : public Function<dim>
+{
+public:
+  virtual double value(const Point<dim> &point,
+                       const unsigned int) const
+  {
+    const double &x = point[0];
+    const double &y = point[1];
+
+    return (Utilities::fixed_power<3>(y) + std::exp(-Utilities::fixed_power<2>(y))
+            + std::sin(4.5*Utilities::fixed_power<2>(y))
+            + std::sin(20*y))*(20*std::cos(4*pi*x) + 0.1*std::sin(20*pi*x)
+                               - 80*std::sin(6*pi*x));
+  }
+};
+
+
+
+template <int dim>
+class HardManufacturedForcing : public Function<dim>
+{
+public:
+  virtual double value(const Point<dim> &point,
+                       const unsigned int) const
+  {
+    const double &x = point[0];
+    const double &y = point[1];
+    return -40.0*(Utilities::fixed_power<3>(y) + std::exp(-Utilities::fixed_power<2>(y)) + std::sin(4.5*Utilities::fixed_power<2>(y)) + std::sin(20*y))*(-8.0*pi*pi*std::cos(4*pi*x) - 1.0*pi*pi*std::sin(20*pi*x) + 72.0*pi*pi*std::sin(6*pi*x)) - (4*Utilities::fixed_power<2>(y)*std::exp(-Utilities::fixed_power<2>(y)) - 81.0*Utilities::fixed_power<2>(y)*std::sin(4.5*Utilities::fixed_power<2>(y)) + 6*y + 9.0*std::cos(4.5*Utilities::fixed_power<2>(y)) - 2*std::exp(-Utilities::fixed_power<2>(y)) - 400*std::sin(20*y))*(20*std::cos(4*pi*x) + 0.1*std::sin(20*pi*x) - 80*std::sin(6*pi*x));
+  }
+};
+
+
+// ----------------------------------------------------------------------------
+// Description of the curved geometry
+// ----------------------------------------------------------------------------
+
+template <int dim>
+class PushForward : public Function<dim>
+{
+public:
+  PushForward() : Function<dim>(dim)
+  {}
+
+  double value(const Point<dim>   &point,
+               const unsigned int  component = 0) const
+  {
+    switch (component)
+      {
+      case 0:
+        return point[0];
+      case 1:
+      {
+        const double &x = point[0];
+        return point[1] + 0.25*(2*x - 1.0)*(x - 1.0)*x;
+      }
+      default:
+        Assert(false, ExcNotImplemented());
+      }
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+};
+
+template <int dim>
+class PullBack : public Function<dim>
+{
+public:
+  PullBack() : Function<dim>(dim)
+  {}
+
+  double value(const Point<dim>   &point,
+               const unsigned int  component = 0) const
+  {
+    switch (component)
+      {
+      case 0:
+        return point[0];
+      case 1:
+      {
+        const double &x = point[0];
+        return point[1] - 0.25*(2*x - 1.0)*(x - 1.0)*x;
+      }
+      default:
+        Assert(false, ExcNotImplemented());
+      }
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+};
+
+/**
+ * A struct containing the push forward and pull back functions. This is
+ * only used with CubicRoofManifold: the only reason this struct exists is
+ * so that CubicRoofManifold holds instances of these functions (which are
+ * fed into FunctionManifold, so they must exist at that point in the
+ * constructor call) and destroys them after destroying its FunctionManifold
+ * part.
+ */
+template <int dim>
+struct CubicRoofFunctions
+{
+  PushForward<dim> forward;
+  PullBack<dim> backward;
+};
+
+/**
+ * A manifold describing the cubic roof.
+ */
+template <int dim>
+class CubicRoofManifold : private CubicRoofFunctions<dim>,
+  public FunctionManifold<dim>
+{
+public:
+  CubicRoofManifold() : FunctionManifold<dim>(this->forward, this->backward)
+  {}
+};
+
+template <int dim>
+std_cxx11::shared_ptr<Manifold<dim> >
+cubic_roof(Triangulation<dim> &triangulation)
+{
+  std_cxx11::shared_ptr<Manifold<dim> > boundary(new CubicRoofManifold<dim>());
+  GridGenerator::hyper_cube(triangulation);
+
+  triangulation.set_all_manifold_ids(cubic_manifold_id);
+  triangulation.set_manifold(cubic_manifold_id, *boundary);
+  return boundary;
+}
+
+// ----------------------------------------------------------------------------
+// The Laplace solver that has trouble with 1 thread, but works well with 2
+// ----------------------------------------------------------------------------
+template <int dim>
+class JxWError
+{
+public:
+  JxWError(const unsigned int n_global_refines);
+
+  double run();
+
+protected:
+  std_cxx11::shared_ptr<Function<dim> > manufactured_solution;
+  std_cxx11::shared_ptr<Function<dim> > manufactured_forcing;
+
+  std_cxx11::shared_ptr<Manifold<dim> > boundary_manifold;
+  Triangulation<dim> triangulation;
+  FE_Q<dim> finite_element;
+  DoFHandler<dim> dof_handler;
+  QGauss<dim> cell_quadrature;
+  MappingQGeneric<dim> cell_mapping;
+
+  ConstraintMatrix all_constraints;
+  SparsityPattern sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+  Vector<double> system_rhs;
+  Vector<double> solution;
+
+  void setup_dofs();
+  void setup_matrices();
+  double solve();
+};
+
+
+
+template <int dim>
+JxWError<dim>::JxWError(const unsigned int n_global_refines) :
+  manufactured_solution(new HardManufacturedSolution<dim>()),
+  manufactured_forcing(new HardManufacturedForcing<dim>()),
+  finite_element(fe_order),
+  dof_handler(triangulation),
+  cell_quadrature(fe_order + 1),
+  cell_mapping(fe_order)
+
+{
+  boundary_manifold = cubic_roof(triangulation);
+  triangulation.refine_global(n_global_refines);
+}
+
+
+
+template <int dim>
+void JxWError<dim>::setup_dofs()
+{
+  dof_handler.distribute_dofs(finite_element);
+  VectorTools::interpolate_boundary_values(cell_mapping,
+                                           dof_handler,
+                                           boundary_id,
+                                           *manufactured_solution,
+                                           all_constraints);
+  all_constraints.close();
+
+  DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
+                                                  dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler,
+                                  dynamic_sparsity_pattern,
+                                  all_constraints,
+                                  /*keep_constrained_dofs=*/false);
+  sparsity_pattern.copy_from(dynamic_sparsity_pattern);
+}
+
+
+template <int dim>
+void JxWError<dim>::setup_matrices()
+{
+  system_matrix.reinit(sparsity_pattern);
+  system_rhs.reinit(dof_handler.n_dofs());
+
+  const UpdateFlags flags = update_values | update_gradients | update_JxW_values
+                            | update_quadrature_points;
+  FEValues<dim> fe_values(cell_mapping, finite_element, cell_quadrature, flags);
+
+  const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
+  FullMatrix<double> cell_system(dofs_per_cell, dofs_per_cell);
+  Vector<double> cell_rhs(dofs_per_cell);
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  for (; cell != endc; ++cell)
+    {
+      cell->get_dof_indices(local_dof_indices);
+      cell_system = 0.0;
+      cell_rhs = 0.0;
+      fe_values.reinit(cell);
+
+      for (unsigned int q_point_n = 0;
+           q_point_n < fe_values.n_quadrature_points; ++q_point_n)
+        {
+          const double point_forcing =
+            manufactured_forcing->value(fe_values.quadrature_point(q_point_n));
+
+          for (unsigned int test_n = 0; test_n < dofs_per_cell; ++test_n)
+            {
+              for (unsigned int trial_n = 0; trial_n < dofs_per_cell; ++trial_n)
+                {
+                  cell_system(test_n, trial_n) +=
+                    fe_values.JxW(q_point_n)*
+                    (fe_values.shape_grad(test_n, q_point_n)
+                     *fe_values.shape_grad(trial_n, q_point_n));
+                }
+
+              cell_rhs[test_n] += fe_values.JxW(q_point_n)
+                                  *fe_values.shape_value(test_n, q_point_n)
+                                  *point_forcing;
+            }
+        }
+
+      all_constraints.distribute_local_to_global(cell_system, cell_rhs,
+                                                 local_dof_indices,
+                                                 system_matrix, system_rhs);
+    }
+}
+
+
+
+template <int dim>
+double JxWError<dim>::solve()
+{
+  {
+    SolverControl solver_control(std::max(types::global_dof_index(100),
+                                          static_cast<types::global_dof_index>(system_rhs.size())),
+                                 1e-14*system_rhs.l2_norm(),
+                                 false,
+                                 false);
+    SolverCG<> solver(solver_control);
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
+
+    solution.reinit(system_rhs);
+    solver.solve(system_matrix, solution, system_rhs, preconditioner);
+    all_constraints.distribute(solution);
+  }
+
+  Vector<double> cell_l2_error(triangulation.n_cells());
+  VectorTools::integrate_difference(cell_mapping,
+                                    dof_handler,
+                                    solution,
+                                    *manufactured_solution,
+                                    cell_l2_error,
+                                    // use QIterated to avoid spurious superconvergence
+                                    QIterated<dim>(QGauss<1>(finite_element.degree), 2),
+                                    VectorTools::L2_norm);
+
+  const double l2_error = VectorTools::compute_global_error
+                          (triangulation, cell_l2_error, VectorTools::L2_norm);
+  return l2_error;
+}
+
+
+
+template <int dim>
+double JxWError<dim>::run()
+{
+  setup_dofs();
+  setup_matrices();
+  const double error = solve();
+  return error;
+}
+
+
+
+int main(int argc, char **argv)
+{
+  // Use exactly one thread so that the CellSimilarity checks are not disabled.
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  static const int dim = 2;
+
+  std::ofstream logfile ("output");
+  deallog << std::setprecision(10);
+  deallog.attach(logfile);
+  for (unsigned int n_global_refines = 3; n_global_refines < 6; ++n_global_refines)
+    {
+      JxWError<dim> solver(n_global_refines);
+      deallog << "L2 error: "
+              << solver.run()
+              << std::endl;
+    }
+}
diff --git a/tests/fe/fe_values_function_manifold.output b/tests/fe/fe_values_function_manifold.output
new file mode 100644 (file)
index 0000000..7e92c1b
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::L2 error: 0.1233048027
+DEAL::L2 error: 0.004171790017
+DEAL::L2 error: 0.0001334762670

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.