]> https://gitweb.dealii.org/ - dealii.git/commitdiff
extend MappingQEulerian to work with multi-grid solution vectors
authorDenis Davydov <davydden@gmail.com>
Wed, 27 Sep 2017 14:25:27 +0000 (16:25 +0200)
committerDenis Davydov <davydden@gmail.com>
Wed, 27 Sep 2017 14:25:27 +0000 (16:25 +0200)
include/deal.II/fe/mapping_q_eulerian.h
source/fe/mapping_q_eulerian.cc
tests/mappings/mapping_q_eulerian_08.cc [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output [new file with mode: 0644]

index 5c1a7e130a98ac29adb159f5c764ea0d619f4257..1b3117f153371af20f7c2791547fa8d0e3f9c2e9 100644 (file)
@@ -43,11 +43,11 @@ template <typename> class Vector;
  * <h3>Usage</h3>
  *
  * The constructor of this class takes three arguments: the polynomial degree
- * of the desire Qp mapping, a reference to the vector that defines the
+ * of the desired Qp mapping, a reference to the vector that defines the
  * mapping from the initial configuration to the current configuration, and a
  * reference to the DoFHandler. The most common case is to use the solution
  * vector for the problem under consideration as the shift vector. The key
- * requirement is that the number of components of the given vector field be
+ * requirement is that the number of components of the given vector field must be
  * equal to (or possibly greater than) the number of space dimensions. If
  * there are more components than space dimensions (for example, if one is
  * working with a coupled problem where there are additional solution
@@ -57,7 +57,7 @@ template <typename> class Vector;
  * the triangulation and associate the desired shift vector to it.
  *
  * Typically, the DoFHandler operates on a finite element that is constructed
- * as a system element (FESystem) from continuous FE_Q() objects. An example
+ * as a system element (FESystem) from continuous FE_Q objects. An example
  * is shown below:
  * @code
  *    FESystem<dim> fe(FE_Q<dim>(2), dim, FE_Q<dim>(1), 1);
@@ -80,7 +80,7 @@ template <typename> class Vector;
  * that whenever you use this object, the given objects still represent valid
  * data.
  *
- * To enable the use of the MappingQ1Eulerian class also in the context of
+ * To enable the use of the MappingQEulerian class also in the context of
  * parallel codes using the PETSc or Trilinos wrapper classes, the type
  * of the vector can be specified as template parameter <tt>VectorType</tt>.
  *
@@ -102,10 +102,14 @@ public:
    * the second argument. The first dim components of this function will be
    * interpreted as the displacement we use in defining the mapping, relative
    * to the location of cells of the underlying triangulation.
+   * @param[in] level. Is the multi-grid level at which the mapping will
+   * be used. It is mainly used to check if the size of the @p euler_vector
+   * is consistent with the @p euler_dof_handler .
    */
   MappingQEulerian (const unsigned int             degree,
                     const DoFHandler<dim,spacedim> &euler_dof_handler,
-                    const VectorType               &euler_vector);
+                    const VectorType               &euler_vector,
+                    const unsigned int             level = numbers::invalid_unsigned_int);
 
   /**
    * Return the mapped vertices of the cell. For the current class, this
@@ -125,15 +129,16 @@ public:
   Mapping<dim,spacedim> *clone () const;
 
   /**
-   * Always returns @p false because MappingQ1Eulerian does not in general
+   * Always returns @p false because MappingQEulerian does not in general
    * preserve vertex locations (unless the translation vector happens to
-   * provide for zero displacements at vertex locations).
+   * provide zero displacements at vertex locations).
    */
   virtual
   bool preserves_vertex_locations () const;
 
   /**
-   * Exception
+   * Exception which is thrown when the mapping is being evaluated at
+   * non-active cell.
    */
   DeclException0 (ExcInactiveCell);
 
@@ -167,6 +172,11 @@ protected:
 
 private:
 
+  /**
+   * Multigrid level at which the mapping is to be used.
+   */
+  const unsigned int level;
+
   /**
    * A class derived from MappingQGeneric that provides the generic mapping
    * with support points on boundary objects so that the corresponding Q3
@@ -202,7 +212,7 @@ private:
     compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
 
     /**
-     * Always returns @p false because MappingQ1Eulerian does not in general
+     * Always returns @p false because MappingQEulerianGeneric does not in general
      * preserve vertex locations (unless the translation vector happens to
      * provide for zero displacements at vertex locations).
      */
@@ -220,14 +230,12 @@ private:
      * Special quadrature rule used to define the support points in the
      * reference configuration.
      */
-
     class SupportQuadrature : public Quadrature<dim>
     {
     public:
       /**
        * Constructor, with an argument defining the desired polynomial degree.
        */
-
       SupportQuadrature (const unsigned int map_degree);
 
     };
index b5a3296689f0125d0275c1d99a1cea2c6569653d..d06ea777af89c8c02117f47317e46995c61b77f0 100644 (file)
@@ -58,11 +58,13 @@ template <int dim, class VectorType, int spacedim>
 MappingQEulerian<dim, VectorType, spacedim>::
 MappingQEulerian (const unsigned int              degree,
                   const DoFHandler<dim,spacedim> &euler_dof_handler,
-                  const VectorType               &euler_vector)
+                  const VectorType               &euler_vector,
+                  const unsigned int              level)
   :
   MappingQ<dim,spacedim>(degree, true),
   euler_vector(&euler_vector),
-  euler_dof_handler(&euler_dof_handler)
+  euler_dof_handler(&euler_dof_handler),
+  level(level)
 {
   // reset the q1 mapping we use for interior cells (and previously
   // set by the MappingQ constructor) to a MappingQ1Eulerian with the
@@ -168,10 +170,13 @@ std::vector<Point<spacedim> >
 MappingQEulerian<dim, VectorType, spacedim>::MappingQEulerianGeneric::
 compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
 {
-  // first, basic assertion with respect to vector size,
+  const bool mg_vector = mapping_q_eulerian.level != numbers::invalid_unsigned_int;
 
-  const types::global_dof_index n_dofs  = mapping_q_eulerian.euler_dof_handler->n_dofs();
+  const types::global_dof_index n_dofs  = mg_vector ?
+                                          mapping_q_eulerian.euler_dof_handler->n_dofs(mapping_q_eulerian.level) :
+                                          mapping_q_eulerian.euler_dof_handler->n_dofs();
   const types::global_dof_index vector_size = mapping_q_eulerian.euler_vector->size();
+
   (void)n_dofs;
   (void)vector_size;
 
@@ -182,20 +187,19 @@ compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell
   typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell,
                                                             mapping_q_eulerian.euler_dof_handler);
 
-  Assert (dof_cell->active() == true, ExcInactiveCell());
+  Assert (mg_vector || dof_cell->active() == true, ExcInactiveCell());
 
   // our quadrature rule is chosen so that each quadrature point corresponds
-  // to a support point in the undeformed configuration.  we can then query
+  // to a support point in the undeformed configuration. We can then query
   // the given displacement field at these points to determine the shift
   // vector that maps the support points to the deformed configuration.
 
-  // we assume that the given field contains dim displacement components, but
+  // We assume that the given field contains dim displacement components, but
   // that there may be other solution components as well (e.g. pressures).
   // this class therefore assumes that the first dim components represent the
   // actual shift vector we need, and simply ignores any components after
-  // that.  this implies that the user should order components appropriately,
+  // that. This implies that the user should order components appropriately,
   // or create a separate dof handler for the displacements.
-
   const unsigned int n_support_pts = support_quadrature.size();
   const unsigned int n_components  = mapping_q_eulerian.euler_dof_handler->get_fe(0).n_components();
 
@@ -205,12 +209,21 @@ compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell
   shift_vector(n_support_pts,
                Vector<typename VectorType::value_type>(n_components));
 
+  std::vector<types::global_dof_index> dof_indices(mapping_q_eulerian.euler_dof_handler->get_fe(0).dofs_per_cell);
   // fill shift vector for each support point using an fe_values object. make
   // sure that the fe_values variable isn't used simultaneously from different
   // threads
   Threads::Mutex::ScopedLock lock(fe_values_mutex);
   fe_values.reinit(dof_cell);
-  fe_values.get_function_values(*mapping_q_eulerian.euler_vector, shift_vector);
+  if (mg_vector)
+    {
+      dof_cell->get_mg_dof_indices(dof_indices);
+      fe_values.get_function_values(*mapping_q_eulerian.euler_vector,
+                                    dof_indices,
+                                    shift_vector);
+    }
+  else
+    fe_values.get_function_values(*mapping_q_eulerian.euler_vector, shift_vector);
 
   // and finally compute the positions of the support points in the deformed
   // configuration.
diff --git a/tests/mappings/mapping_q_eulerian_08.cc b/tests/mappings/mapping_q_eulerian_08.cc
new file mode 100644 (file)
index 0000000..2f1bc4f
--- /dev/null
@@ -0,0 +1,320 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test that MappingQEulerian works in parallel with geometric multigrids.
+// We apply a simple linear deformation which can be represented exactly
+// at the coarse level.
+//
+// inspired by _07
+
+#include "../tests.h"
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/identity_matrix.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/filtered_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_reordering.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <iostream>
+#include <vector>
+
+
+
+
+
+
+using namespace dealii;
+
+
+template <int dim>
+class Displacement : public Function<dim>
+{
+public:
+  Displacement() :
+    Function<dim>(dim)
+  {}
+
+  double value (const Point<dim> &p,
+                const unsigned int component) const
+  {
+    return (0.1 + 2.*component) + (0.2 + 3.*component) * p[component];
+  }
+
+  template <typename NumberType>
+  Tensor<1,dim, VectorizedArray<NumberType> >
+  shift_value(const Point<dim, VectorizedArray<NumberType>> &p_vec) const
+  {
+    Tensor<1,dim, VectorizedArray<NumberType>> shift_vec;
+    Point<dim> p;
+    for (unsigned int v = 0; v < VectorizedArray<NumberType>::n_array_elements; ++v)
+      {
+        for (unsigned int d = 0; d < dim; ++d)
+          p[d] = p_vec[d][v];
+
+        for (unsigned int d = 0; d < dim; ++d)
+          shift_vec[d][v] = this->value(p,d);
+      }
+
+    return shift_vec;
+
+  }
+};
+
+
+template <int dim, int fe_degree=2, int n_q_points=fe_degree+1, typename NumberType=double, typename LevelNumberType=NumberType>
+void test (const unsigned int n_ref = 0)
+{
+  Displacement<dim> displacement_function;
+  const unsigned int euler_fe_degree=2;
+
+  deallog << "dim=" << dim << std::endl;
+  MPI_Comm mpi_communicator(MPI_COMM_WORLD);
+  unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator);
+  unsigned int numproc = Utilities::MPI::n_mpi_processes (mpi_communicator);
+
+  deallog << "numproc=" << numproc << std::endl;
+
+  parallel::distributed::Triangulation<dim>
+  triangulation (mpi_communicator,
+                 Triangulation<dim>::limit_level_difference_at_vertices,
+                 parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global(3);
+
+  // do some adaptive refinement
+  for (unsigned int ref=0; ref<n_ref; ++ref)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator cell=triangulation.begin_active();
+           cell != triangulation.end(); ++cell)
+        if (cell->is_locally_owned() &&
+            (cell->center().norm() < 0.5 && (cell->level() < 5 ||
+                                             cell->center().norm() > 0.45)
+             ||
+             (dim == 2 && cell->center().norm() > 1.2)))
+          cell->set_refine_flag();
+      triangulation.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof_handler (triangulation);
+  dof_handler.distribute_dofs(fe);
+  dof_handler.distribute_mg_dofs();
+
+  // quadrature for MatrixFree, not related to the degree in the
+  // FE displacement field.
+  QGauss<1> quadrature_formula(n_q_points);
+
+
+  FESystem<dim> fe_euler(FE_Q<dim>(euler_fe_degree), dim);
+  DoFHandler<dim> dof_handler_euler(triangulation);
+  dof_handler_euler.distribute_dofs(fe_euler);
+  dof_handler_euler.distribute_mg_dofs ();
+
+  const IndexSet &locally_owned_dofs_euler = dof_handler_euler.locally_owned_dofs ();
+  IndexSet locally_relevant_dofs_euler;
+  DoFTools::extract_locally_relevant_dofs (dof_handler_euler,locally_relevant_dofs_euler);
+
+  const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs ();
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs);
+
+  // constraints:
+  ConstraintMatrix constraints_euler;
+  constraints_euler.reinit (locally_relevant_dofs_euler);
+  DoFTools::make_hanging_node_constraints  (dof_handler_euler, constraints_euler);
+  constraints_euler.close ();
+
+  ConstraintMatrix constraints;
+  constraints.reinit (locally_relevant_dofs);
+  DoFTools::make_hanging_node_constraints  (dof_handler, constraints);
+  constraints_euler.close ();
+
+  // MG constraints:
+  MGConstrainedDoFs mg_constrained_dofs_euler;
+  mg_constrained_dofs_euler.initialize(dof_handler_euler);
+
+  // Displacement vector
+  LinearAlgebra::distributed::Vector<NumberType> displacement;
+  displacement.reinit(locally_owned_dofs_euler,
+                      locally_relevant_dofs_euler,
+                      mpi_communicator);
+
+  VectorTools::project<dim,LinearAlgebra::distributed::Vector<NumberType>,dim>
+  (dof_handler_euler,
+   constraints_euler,
+   QGauss<dim>(n_q_points),
+   displacement_function,
+   displacement);
+  displacement.update_ghost_values();
+
+  MGTransferMatrixFree<dim, LevelNumberType> mg_transfer_euler(mg_constrained_dofs_euler);
+  mg_transfer_euler.build(dof_handler);
+
+  // now the core of the test:
+  const unsigned int max_level = dof_handler.get_triangulation().n_global_levels()-1;
+  const unsigned int min_level = 0;
+  MGLevelObject<LinearAlgebra::distributed::Vector<LevelNumberType>> displacement_level(min_level, max_level);
+  mg_transfer_euler.interpolate_to_mg(dof_handler_euler,displacement_level, displacement);
+
+  // First, check fine-level only:
+  {
+    MappingQEulerian<dim,LinearAlgebra::distributed::Vector<NumberType>>
+        euler_fine(euler_fe_degree, dof_handler_euler, displacement);
+
+
+
+    MatrixFree<dim,NumberType> matrix_free_euler;
+    typename MatrixFree<dim,NumberType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim,NumberType>::AdditionalData::partition_color;
+    data.mapping_update_flags = update_values | update_gradients | update_JxW_values | update_quadrature_points;
+    matrix_free_euler.reinit (euler_fine, dof_handler, constraints, quadrature_formula, data);
+
+    MatrixFree<dim,NumberType> matrix_free;
+    matrix_free.reinit (dof_handler, constraints, quadrature_formula, data);
+
+
+    // test fine-level mapping:
+    {
+      FEEvaluation<dim,fe_degree,n_q_points,1,NumberType> fe_eval_euler(matrix_free_euler);
+      FEEvaluation<dim,fe_degree,n_q_points,1,NumberType> fe_eval(matrix_free);
+      const unsigned int n_cells = matrix_free_euler.n_macro_cells();
+      Assert (matrix_free_euler.n_macro_cells() == matrix_free.n_macro_cells(),
+              ExcInternalError());
+      const unsigned int nqp = fe_eval.n_q_points;
+      for (unsigned int cell=0; cell<n_cells; ++cell)
+        {
+          fe_eval_euler.reinit(cell);
+          fe_eval.reinit(cell);
+          for (unsigned int q=0; q<nqp; ++q)
+            {
+              const auto &v1 = fe_eval_euler.quadrature_point(q);
+              const auto &qp = fe_eval.quadrature_point(q);
+              const auto v2 = qp + displacement_function.shift_value(qp);
+              VectorizedArray<NumberType> dist = v1.distance(v2);
+              for (unsigned int v = 0; v < VectorizedArray<NumberType>::n_array_elements; ++v)
+                AssertThrow (dist[v] < 1e-8,
+                             ExcMessage("distance: " + std::to_string(dist[v])));
+            }
+        }
+    }
+
+  }
+
+  // now go through all GMG levels:
+  std::set<types::boundary_id> dirichlet_boundary_ids;
+  dirichlet_boundary_ids.insert(0);
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof_handler);
+  mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,dirichlet_boundary_ids);
+
+  for (unsigned int level = min_level; level <= max_level; ++level)
+    {
+      typename MatrixFree<dim,LevelNumberType>::AdditionalData mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme = MatrixFree<dim,LevelNumberType>::AdditionalData::partition_color;
+      mg_additional_data.level_mg_handler = level;
+      mg_additional_data.mapping_update_flags = update_values | update_gradients | update_JxW_values | update_quadrature_points;
+
+      ConstraintMatrix level_constraints;
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof_handler, level,
+                                                    relevant_dofs);
+      level_constraints.reinit(relevant_dofs);
+      level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
+      level_constraints.close();
+
+      MatrixFree<dim,LevelNumberType> mg_level_euler;
+
+      MappingQEulerian<dim,LinearAlgebra::distributed::Vector<LevelNumberType>>
+          euler_level(euler_fe_degree, dof_handler_euler, displacement_level[level], level);
+
+      mg_level_euler.reinit (euler_level, dof_handler, level_constraints, quadrature_formula,
+                             mg_additional_data);
+
+      MatrixFree<dim,LevelNumberType> mg_level;
+      mg_level.reinit (dof_handler, level_constraints, quadrature_formula,
+                       mg_additional_data);
+
+      // go through all cells and quadrature points:
+      {
+        FEEvaluation<dim,fe_degree,n_q_points,1,NumberType> fe_eval_euler(mg_level_euler);
+        FEEvaluation<dim,fe_degree,n_q_points,1,NumberType> fe_eval(mg_level);
+        const unsigned int n_cells = mg_level_euler.n_macro_cells();
+        Assert (mg_level_euler.n_macro_cells() == mg_level.n_macro_cells(),
+                ExcInternalError());
+        const unsigned int nqp = fe_eval.n_q_points;
+        for (unsigned int cell=0; cell<n_cells; ++cell)
+          {
+            fe_eval_euler.reinit(cell);
+            fe_eval.reinit(cell);
+            for (unsigned int q=0; q<nqp; ++q)
+              {
+                const auto &v1 = fe_eval_euler.quadrature_point(q);
+                const auto &qp = fe_eval.quadrature_point(q);
+                const auto v2 = qp + displacement_function.shift_value(qp);
+                VectorizedArray<NumberType> dist = v1.distance(v2);
+                for (unsigned int v = 0; v < VectorizedArray<NumberType>::n_array_elements; ++v)
+                  AssertThrow (dist[v] < 1e-8,
+                               ExcMessage("Level " + std::to_string(level) +
+                                          " distance: " + std::to_string(dist[v])));
+              }
+          }
+      }
+    }
+
+  deallog << "Ok" << std::endl;
+}
+
+
+
+int main (int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  MPILogInitAll log;
+
+  test<2> ();
+}
+
+
diff --git a/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..047c058
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:0::dim=2
+DEAL:0::numproc=1
+DEAL:0::Ok
diff --git a/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output b/tests/mappings/mapping_q_eulerian_08.with_p4est=true.mpirun=3.output
new file mode 100644 (file)
index 0000000..d171703
--- /dev/null
@@ -0,0 +1,14 @@
+
+DEAL:0::dim=2
+DEAL:0::numproc=3
+DEAL:0::Ok
+
+DEAL:1::dim=2
+DEAL:1::numproc=3
+DEAL:1::Ok
+
+
+DEAL:2::dim=2
+DEAL:2::numproc=3
+DEAL:2::Ok
+

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.