]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add FieldTransfer for the transfer of field with FE_Nothing
authorBruno Turcksin <bruno.turcksin@gmail.com>
Tue, 14 Sep 2021 18:17:50 +0000 (18:17 +0000)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Fri, 15 Jul 2022 12:56:16 +0000 (12:56 +0000)
include/deal.II/distributed/field_transfer.h [new file with mode: 0644]
source/distributed/CMakeLists.txt
source/distributed/field_transfer.cc [new file with mode: 0644]
source/distributed/field_transfer.inst.in [new file with mode: 0644]

diff --git a/include/deal.II/distributed/field_transfer.h b/include/deal.II/distributed/field_transfer.h
new file mode 100644 (file)
index 0000000..6b5ef36
--- /dev/null
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_distributed_field_transfer_h
+#define dealii_distributed_field_transfer_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/distributed/cell_data_transfer.templates.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <memory>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace parallel
+{
+  namespace distributed
+  {
+    /**
+     * Classes and functions in the experimental namespace can be modified or
+     * removed without being deprecated first.
+     */
+    namespace experimental
+    {
+      /**
+       * This class is similar to SolutionTransfer but it supports the case
+       * where elements have been activated during refinement, i.e., FE_Nothing
+       * elements have been associated with a finite elements during refinement.
+       */
+      template <int dim, typename VectorType, int spacedim = dim>
+      class FieldTransfer
+      {
+      private:
+        using Number = typename VectorType::value_type;
+
+      public:
+        /**
+         * 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.
+         */
+        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
+         * @param[in] fe_nothing_index The finite element index associated with
+         * FE_Nothing
+         */
+        void
+        prepare_for_coarsening_and_refinement(
+          const VectorType & in,
+          const unsigned int fe_nothing_index);
+
+        /**
+         * Interpolate the data previously stored in this object before the mesh
+         * was refined or coarsenend 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.
+         */
+        void
+        interpolate(const Number &                   new_value,
+                    const AffineConstraints<Number> &affine_constraints,
+                    VectorType &                     out);
+
+      private:
+        /**
+         * DoFHandler associated with the object.
+         */
+        const DoFHandler<dim, spacedim> &dof_handler;
+
+        /**
+         * Data transfered by cell_data_transfer.
+         */
+        std::vector<Vector<Number>> data_to_transfer;
+
+        /**
+         * CellDataTransfer used to perform the field transfer.
+         */
+        std::unique_ptr<
+          CellDataTransfer<dim, spacedim, std::vector<Vector<Number>>>>
+          cell_data_transfer;
+      };
+    } // namespace experimental
+  }   // namespace distributed
+} // namespace parallel
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 836602ef8dffcfad2a8b2417354ca0a6a04342bf..7754496a0e166e9f3e11c1d5b08a2121f17e749c 100644 (file)
@@ -26,6 +26,7 @@ SET(_unity_include_src
   tria_base.cc
   shared_tria.cc
   p4est_wrappers.cc
+  field_transfer.cc
   )
 
 SET(_separate_src
@@ -44,6 +45,7 @@ SET(_inst
   grid_refinement.inst.in
   cell_weights.inst.in
   cell_data_transfer.inst.in
+  field_transfer.inst.in
   fully_distributed_tria.inst.in
   repartitioning_policy_tools.inst.in
   solution_transfer.inst.in
diff --git a/source/distributed/field_transfer.cc b/source/distributed/field_transfer.cc
new file mode 100644 (file)
index 0000000..75b2d32
--- /dev/null
@@ -0,0 +1,242 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/distributed/field_transfer.h>
+
+#ifdef DEAL_II_WITH_P4EST
+
+#  include <deal.II/lac/lapack_full_matrix.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace parallel
+{
+  namespace distributed
+  {
+    namespace experimental
+    {
+      template <int dim, typename VectorType, int spacedim>
+      FieldTransfer<dim, VectorType, spacedim>::FieldTransfer(
+        const DoFHandler<dim, spacedim> &dof)
+        : dof_handler(dof)
+      {
+        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())));
+      }
+
+
+
+      template <int dim, typename VectorType, int spacedim>
+      void
+      FieldTransfer<dim, VectorType, spacedim>::
+        prepare_for_coarsening_and_refinement(
+          const VectorType & in,
+          const unsigned int fe_nothing_index)
+      {
+        const unsigned int dofs_per_cell =
+          dof_handler.get_fe_collection().max_dofs_per_cell();
+
+        Vector<Number> cell_solution(dofs_per_cell);
+        Vector<Number> dummy_cell_solution(dofs_per_cell);
+
+        for (auto &val : dummy_cell_solution)
+          {
+            val = std::numeric_limits<Number>::infinity();
+          }
+
+        in.update_ghost_values();
+
+        std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          {
+            if ((cell->is_locally_owned()) &&
+                (cell->active_fe_index() != fe_nothing_index))
+              {
+                cell->get_dof_values(in, cell_solution);
+                data_to_transfer.push_back(cell_solution);
+              }
+            else
+              {
+                data_to_transfer.push_back(dummy_cell_solution);
+              }
+          }
+
+        cell_data_transfer->prepare_for_coarsening_and_refinement(
+          data_to_transfer);
+      }
+
+
+
+      template <int dim, typename VectorType, int spacedim>
+      void
+      FieldTransfer<dim, VectorType, spacedim>::interpolate(
+        const Number &                   new_value,
+        const AffineConstraints<Number> &affine_constraints,
+        VectorType &                     out)
+      {
+        const unsigned int dofs_per_cell =
+          dof_handler.get_fe_collection().max_dofs_per_cell();
+        std::vector<Vector<Number>> transferred_data(
+          dof_handler.get_triangulation().n_active_cells(),
+          Vector<Number>(dofs_per_cell));
+
+        cell_data_transfer->unpack(transferred_data);
+
+        // Free the memory allocated by data_to_transfer
+        data_to_transfer.clear();
+
+        for (unsigned int i = 0; i < out.locally_owned_size(); ++i)
+          out.local_element(i) = std::numeric_limits<Number>::infinity();
+
+        unsigned int cell_i = 0;
+        for (auto const &cell : dof_handler.active_cell_iterators())
+          {
+            if ((cell->is_locally_owned()) &&
+                (transferred_data[cell_i][0] !=
+                 std::numeric_limits<Number>::infinity()))
+              {
+                cell->set_dof_values(transferred_data[cell_i], out);
+              }
+            ++cell_i;
+          }
+
+
+        // Communicate the results.
+        out.compress(dealii::VectorOperation::min);
+
+        // Treat hanging nodes
+        std::vector<types::global_dof_index> dof_indices;
+        std::vector<types::global_dof_index> dofs_map;
+        std::vector<std::vector<std::pair<types::global_dof_index, Number>>>
+                            constraint_lines;
+        std::vector<Number> constraint_values;
+        IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+        for (auto constrained_dof : locally_owned_dofs)
+          if (affine_constraints.is_constrained(constrained_dof))
+            {
+              auto *constraint =
+                affine_constraints.get_constraint_entries(constrained_dof);
+              const unsigned int line_size = constraint->size();
+              bool               add_line  = false;
+              for (unsigned int i = 0; i < line_size; ++i)
+                {
+                  types::global_dof_index constraining_dof =
+                    (*constraint)[i].first;
+                  // If one of the constraining value is infinity, we need to
+                  // reverse the relationship
+                  if (out[constraining_dof] ==
+                      std::numeric_limits<Number>::infinity())
+                    {
+                      add_line = true;
+                      break;
+                    }
+                }
+
+              if (add_line)
+                {
+                  std::vector<std::pair<types::global_dof_index, Number>> line;
+                  Number val = out[constrained_dof];
+
+                  for (unsigned int i = 0; i < line_size; ++i)
+                    {
+                      types::global_dof_index constraining_dof =
+                        (*constraint)[i].first;
+
+                      if (out[constraining_dof] ==
+                          std::numeric_limits<Number>::infinity())
+                        {
+                          auto constraining_dof_map_it =
+                            std::find(dofs_map.begin(),
+                                      dofs_map.end(),
+                                      constraining_dof);
+                          if (constraining_dof_map_it == dofs_map.end())
+                            {
+                              dofs_map.push_back(constraining_dof);
+                            }
+                          line.push_back((*constraint)[i]);
+                        }
+                      else
+                        {
+                          val -=
+                            out[constraining_dof] * (*constraint)[i].second;
+                        }
+                    }
+                  constraint_lines.push_back(line);
+                  constraint_values.push_back(val);
+                }
+            }
+
+        // Build a constraint matrix that we invert
+        const unsigned int n_rows = constraint_lines.size();
+        if (n_rows > 0)
+          {
+            const unsigned int               n_cols = dofs_map.size();
+            dealii::LAPACKFullMatrix<Number> constraints_matrix(n_rows, n_cols);
+            dealii::Vector<Number>           constraint_values_vec(n_rows);
+
+            for (unsigned int i = 0; i < n_rows; ++i)
+              {
+                for (unsigned int j = 0; j < n_cols; ++j)
+                  {
+                    if (j < constraint_lines[i].size())
+                      {
+                        auto constraint_it =
+                          std::find(dofs_map.begin(),
+                                    dofs_map.end(),
+                                    constraint_lines[i][j].first);
+                        constraints_matrix(i,
+                                           constraint_it - dofs_map.begin()) =
+                          constraint_lines[i][j].second;
+                      }
+                  }
+
+                constraint_values_vec[i] = constraint_values[i];
+              }
+
+            constraints_matrix.compute_inverse_svd();
+
+            dealii::Vector<Number> new_constrained_values(n_cols);
+            constraints_matrix.vmult(new_constrained_values,
+                                     constraint_values_vec);
+
+            for (unsigned int i = 0; i < n_cols; ++i)
+              {
+                out[dofs_map[i]] = new_constrained_values[i];
+              }
+          }
+
+        // Set the value to the newly create DoFs.
+        std::for_each(out.begin(), out.end(), [&](Number &val) {
+          if (val == std::numeric_limits<Number>::infinity())
+            {
+              val = new_value;
+            }
+        });
+      }
+    } // namespace experimental
+  }   // namespace distributed
+} // namespace parallel
+
+// explicit instantiations
+#  include "field_transfer.inst"
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/source/distributed/field_transfer.inst.in b/source/distributed/field_transfer.inst.in
new file mode 100644 (file)
index 0000000..3dcbc62
--- /dev/null
@@ -0,0 +1,39 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+    namespace parallel
+    \{
+      namespace distributed
+      \{
+        namespace experimental
+        \{
+#if deal_II_dimension <= deal_II_space_dimension
+          template class FieldTransfer<
+            deal_II_dimension,
+            ::dealii::LinearAlgebra::distributed::Vector<double>,
+            deal_II_space_dimension>;
+          template class FieldTransfer<
+            deal_II_dimension,
+            ::dealii::LinearAlgebra::distributed::Vector<float>,
+            deal_II_space_dimension>;
+#endif
+        \}
+      \}
+    \}
+  }

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.