]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FieldTransfer: enable coarsening 14802/head
authorSebastian Proell <sebastian.proell@tum.de>
Tue, 14 Feb 2023 18:10:22 +0000 (19:10 +0100)
committerSebastian Proell <sebastian.proell@tum.de>
Wed, 15 Feb 2023 17:25:17 +0000 (18:25 +0100)
include/deal.II/distributed/field_transfer.h
source/distributed/field_transfer.cc
tests/hp/field_transfer_05.cc [new file with mode: 0644]
tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output [new file with mode: 0644]

index 0793e4fcce7ce6d99c8211872716b40b89583638..95c0ec8f8fb6eb61ad54d2cf27afcdc38e1b0a9a 100644 (file)
@@ -52,17 +52,11 @@ namespace parallel
          * Constructor.
          *
          * @param[in] dof_handler The DoFHandler on which all the operations
-         * will happen. At the time when this constructor is call, the
-         * DoFHandler still points to the Triangulation before the refinement in
-         * question happens.
+         * will happen. This constructor must be called before the underlying
+         * Triangulation is coarsened/refined.
          */
         FieldTransfer(const DoFHandler<dim, spacedim> &dof_handler);
 
-        /**
-         * Destructor.
-         */
-        ~FieldTransfer() = default;
-
         /**
          * Prepare the current object for coarsening and refinement.
          * @param[in] in The vector that will be interpolated
@@ -76,10 +70,10 @@ namespace parallel
 
         /**
          * Interpolate the data previously stored in this object before the mesh
-         * was refined or coarsenend onto the current set of cells. @p new_value
+         * was refined or coarsened onto the current set of cells. @p new_value
          * is the value associated to the new degrees of freedom that where
          * created during the element activation. @p affine_constraints is the
-         * AffinieConstraints after refinement.
+         * AffineConstraints after refinement.
          */
         void
         interpolate(const Number &                   new_value,
index 2c8366d38ae47940c67d0318c7b5fa503f3c9f6f..2f633debc22473a9bfd45055ecbddec3e840bda5 100644 (file)
@@ -34,12 +34,96 @@ namespace parallel
         const DoFHandler<dim, spacedim> &dof)
         : dof_handler(dof)
       {
+        // When coarsening, we want to mimic the behavior of SolutionTransfer
+        // and interpolate from child cells to parent. Define this strategy here
+        // since it is not readily available
+        const auto coarsening_strategy =
+          [this](
+            const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+              &                                parent,
+            const std::vector<Vector<Number>> &children_values) {
+            // get the equivalent DoFCellAccessor
+            typename DoFHandler<dim, spacedim>::cell_iterator dof_cell_iterator(
+              &dof_handler.get_triangulation(),
+              parent->level(),
+              parent->index(),
+              &dof_handler);
+
+            int fe_index = 0;
+            if (dof_handler.has_hp_capabilities())
+              fe_index = dealii::internal::hp::DoFHandlerImplementation::
+                dominated_future_fe_on_children<dim, spacedim>(
+                  dof_cell_iterator);
+
+            const auto &fe = dof_handler.get_fe(fe_index);
+            Assert(fe.n_dofs_per_cell() > 0,
+                   ExcMessage(
+                     "Cannot coarsen onto a FiniteElement with no DoFs."));
+            AssertDimension(dof_cell_iterator->n_children(),
+                            children_values.size());
+
+            const auto child_iterators = dof_cell_iterator->child_iterators();
+            const unsigned int n_children_with_fe_nothing =
+              std::count_if(child_iterators.begin(),
+                            child_iterators.end(),
+                            [](const auto &child_cell) {
+                              return child_cell->get_fe().n_dofs_per_cell() ==
+                                     0;
+                            });
+
+            Assert(
+              n_children_with_fe_nothing == 0 ||
+                n_children_with_fe_nothing == dof_cell_iterator->n_children(),
+              ExcMessage(
+                "Coarsening is only supported for parent cells where either all"
+                " or none of the child cells are FE_Nothing."));
+
+            // in case all children are FE_Nothing there is nothing to
+            // interpolate and we just return the first entry from the children
+            // values (containing invalid entries)
+            if (n_children_with_fe_nothing == dof_cell_iterator->n_children())
+              {
+                return children_values[0];
+              }
+
+            const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+            Vector<Number>     tmp(dofs_per_cell);
+            Vector<Number>     interpolated_values(dofs_per_cell);
+
+            // Otherwise, perform the actual interpolation here. Due to the
+            // assert above, we know that all child cells have data to
+            // interpolate.
+            for (unsigned int child = 0;
+                 child < dof_cell_iterator->n_children();
+                 ++child)
+              {
+                // interpolate the previously stored values on a child to the
+                // mother cell
+                fe.get_restriction_matrix(child,
+                                          dof_cell_iterator->refinement_case())
+                  .vmult(tmp, children_values[child]);
+
+                // and add up or set them in the output vector
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  if (fe.restriction_is_additive(i))
+                    interpolated_values(i) += tmp(i);
+                  else if (tmp(i) != Number())
+                    interpolated_values(i) = tmp(i);
+              }
+
+            return interpolated_values;
+          };
+
         cell_data_transfer = std::make_unique<
           CellDataTransfer<dim, spacedim, std::vector<Vector<Number>>>>(
           dynamic_cast<
             dealii::parallel::distributed::Triangulation<dim, spacedim> &>(
             const_cast<dealii::Triangulation<dim, spacedim> &>(
-              dof_handler.get_triangulation())));
+              dof_handler.get_triangulation())),
+          false,
+          &dealii::AdaptationStrategies::Refinement::
+            preserve<dim, spacedim, Vector<Number>>,
+          coarsening_strategy);
       }
 
 
diff --git a/tests/hp/field_transfer_05.cc b/tests/hp/field_transfer_05.cc
new file mode 100644 (file)
index 0000000..c48f5f9
--- /dev/null
@@ -0,0 +1,170 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+// Check that we can use FieldTransfer when coarsening cells with non-uniform
+// data
+
+#include <deal.II/distributed/field_transfer.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+class TestFunction : public Function<2>
+{
+public:
+  double
+  value(const dealii::Point<2> &p, const unsigned int component) const override
+  {
+    return std::exp(10 - 10 * p[1]);
+  }
+};
+
+void
+output(const DoFHandler<2> &                             dof_handler,
+       const LinearAlgebra::distributed::Vector<double> &vector,
+       const std::string &                               file_name)
+{
+  DataOut<2> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(vector, "solution");
+  data_out.build_patches();
+  data_out.write_vtu_in_parallel(file_name, MPI_COMM_WORLD);
+}
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MPILogInitAll mpi_log;
+
+  parallel::distributed::Triangulation<2> triangulation(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(triangulation);
+  triangulation.refine_global(3);
+
+  hp::FECollection<2> fe_collection;
+  fe_collection.push_back(FE_Q<2>(1));
+  fe_collection.push_back(FE_Nothing<2>());
+
+  DoFHandler<2> dof_handler(triangulation);
+
+  // Assign FE_Nothing to half of the domain
+  for (auto cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->is_locally_owned())
+        {
+          if (cell->center()[0] < 0.5)
+            {
+              cell->set_active_fe_index(0);
+            }
+          else
+            {
+              cell->set_active_fe_index(1);
+            }
+        }
+    }
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  // Initialize solution
+  auto locally_relevant_dofs =
+    DoFTools::extract_locally_relevant_dofs(dof_handler);
+  LinearAlgebra::distributed::Vector<double> solution(
+    dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
+
+  // Initialize with a non-constant function
+  VectorTools::interpolate(dof_handler, TestFunction(), solution);
+
+  {
+    if (false)
+      output(dof_handler, solution, "result_0.vtu");
+
+    deallog << "solution before: " << std::endl;
+    std::stringstream ss;
+    solution.print(ss);
+    deallog << ss.str() << std::endl;
+  }
+
+  parallel::distributed::experimental::
+    FieldTransfer<2, LinearAlgebra::distributed::Vector<double>>
+      field_transfer(dof_handler);
+  // Assign FE_Q to all the cells
+  for (auto cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->is_locally_owned())
+        {
+          cell->set_future_fe_index(0);
+        }
+    }
+
+  // coarsen some cells
+  for (auto cell : dof_handler.active_cell_iterators())
+    {
+      if (cell->is_locally_owned())
+        {
+          if (cell->center()[0] < 0.25 && cell->center()[1] < 0.25)
+            {
+              cell->set_coarsen_flag();
+            }
+        }
+    }
+
+  triangulation.prepare_coarsening_and_refinement();
+
+  field_transfer.prepare_for_coarsening_and_refinement(solution, 1);
+
+  triangulation.execute_coarsening_and_refinement();
+
+  dof_handler.distribute_dofs(fe_collection);
+
+  locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
+  LinearAlgebra::distributed::Vector<double> new_solution(
+    dof_handler.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
+
+  AffineConstraints<double> affine_constraints;
+  DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints);
+  affine_constraints.close();
+
+  const double new_value = 11.;
+
+  field_transfer.interpolate(new_value, affine_constraints, new_solution);
+
+  affine_constraints.distribute(new_solution);
+
+  {
+    if (false)
+      output(dof_handler, new_solution, "result_1.vtu");
+
+    deallog << "solution after coarsening and activation: " << std::endl;
+    std::stringstream ss;
+    new_solution.print(ss);
+    deallog << ss.str() << std::endl;
+  }
+}
diff --git a/tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output b/tests/hp/field_transfer_05.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..5bde4cb
--- /dev/null
@@ -0,0 +1,27 @@
+
+DEAL:0::solution before: 
+DEAL:0::Process #0
+Local range: [0, 25), global size: 45
+Vector data:
+2.203e+04 2.203e+04 6.311e+03 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 2.203e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 
+
+DEAL:0::solution after coarsening and activation: 
+DEAL:0::Process #0
+Local range: [0, 42), global size: 78
+Vector data:
+2.203e+04 2.203e+04 1.808e+03 1.808e+03 2.203e+04 1.192e+04 6.311e+03 2.203e+04 6.311e+03 1.808e+03 1.808e+03 1.808e+03 5.180e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.484e+02 5.180e+02 5.180e+02 1.484e+02 1.484e+02 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 
+
+
+DEAL:1::solution before: 
+DEAL:1::Process #1
+Local range: [25, 45), global size: 45
+Vector data:
+4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 
+
+DEAL:1::solution after coarsening and activation: 
+DEAL:1::Process #1
+Local range: [42, 78), global size: 78
+Vector data:
+4.252e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 1.218e+01 4.252e+01 4.252e+01 1.218e+01 1.218e+01 3.490e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.000e+00 3.490e+00 3.490e+00 1.000e+00 1.000e+00 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 1.100e+01 
+
+

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.