]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Improve asserts in parallel::distributed::SolutionTransfer 12463/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 25 Feb 2022 18:18:04 +0000 (19:18 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 10 May 2022 09:33:15 +0000 (11:33 +0200)
doc/news/changes/minor/20220225Munch [new file with mode: 0644]
include/deal.II/base/partitioner.templates.h
include/deal.II/dofs/dof_accessor.h
source/distributed/solution_transfer.cc
source/dofs/dof_accessor_set.cc
source/dofs/dof_accessor_set.inst.in
source/numerics/solution_transfer.cc
tests/codim_one/solution_transfer_01.cc

diff --git a/doc/news/changes/minor/20220225Munch b/doc/news/changes/minor/20220225Munch
new file mode 100644 (file)
index 0000000..1e461f0
--- /dev/null
@@ -0,0 +1,9 @@
+Improved: The asserts in parallel::distributed::SolutionTransfer have
+been improved. Before this change, an assert has been only called
+in the parallel case for distributed vectors for cells neighboring
+interprocess boundaries. Now, we check for all DoFs that the 
+contributions from all cells have the same values. This condition
+might be not given if the solution has high gradients at hanging
+nodes.
+<br>
+(Peter Munch, Magdalena Schreter, 2022/02/25)
index 221b3e55420e38935d3161b86b0d61d92b5fdab6..34fe115dacfd9b7581ffe94f7b4654edfd3757bc 100644 (file)
@@ -620,7 +620,10 @@ namespace Utilities
                 // p::d::SolutionTransfer. The rationale is that during
                 // interpolation on two elements sharing the face, values on
                 // this face obtained from each side might be different due to
-                // additions being done in different order.
+                // additions being done in different order. If the local
+                // value is zero, it indicates that the local process has not
+                // set the value during the cell loop and its value can be
+                // safely overridden.
                 Assert(*read_position == Number() ||
                          internal::get_abs(locally_owned_array[j] -
                                            *read_position) <=
index 39ee53ef313e868b4a8e05dfdb4df22a4ef953e5..2f561c44ed455574c7987130a6769feb52e01721 100644 (file)
@@ -1743,6 +1743,15 @@ public:
    * of the respective finite element class for a description of what the
    * prolongation matrices represent in this case.
    *
+   * @note Cells set the values of DoFs independently and might overwrite
+   * previously set values in the global vector, for example when calling the
+   * same function earlier from a different cell. By setting @p perform_check,
+   * you can enable a check that the previous value and the one to be set here
+   * are at least roughly the same. In practice, they might be slightly
+   * different because they are computed in a way that theoretically ensures
+   * that they are the same, but in practice they are only equal up to
+   * round-off.
+   *
    * @note Unlike the get_dof_values() function, this function is only
    * available on cells, rather than on lines, quads, and hexes, since
    * interpolation is presently only provided for cells by the finite element
@@ -1754,7 +1763,8 @@ public:
     const Vector<number> &local_values,
     OutputVector &        values,
     const unsigned int    fe_index =
-      DoFHandler<dimension_, space_dimension_>::invalid_fe_index) const;
+      DoFHandler<dimension_, space_dimension_>::invalid_fe_index,
+    const bool perform_check = false) const;
 
   /**
    * Distribute a local (cell based) vector to a global one by mapping the
index 99dcf55b42e741f6e0b66e7eb876277d803ec769..889658aad0ce08e4ebf412203d108f802d09dcf7 100644 (file)
@@ -424,7 +424,8 @@ namespace parallel
       for (; it_input != dof_values.cend(); ++it_input, ++it_output)
         cell->set_dof_values_by_interpolation(*it_input,
                                               *(*it_output),
-                                              fe_index);
+                                              fe_index,
+                                              true);
     }
   } // namespace distributed
 } // namespace parallel
index ca5dc3cf9633ed1988fbbbbe3daff49795516761..bd2f50e5c68fe30648f8c424305a7dfbb0c48ec0 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+template <typename Number>
+DeclException2(ExcNonMatchingElementsSetDofValuesByInterpolation,
+               Number,
+               Number,
+               << "Called set_dof_values_by_interpolation(), but"
+               << " the element to be set, value " << std::setprecision(16)
+               << arg1 << ", does not match with the non-zero value "
+               << std::setprecision(16) << arg2 << " already set before.");
+
+namespace internal
+{
+  /**
+   * In the set_dof_values(), we need to invoke abs() also on unsigned data
+   * types, which is ill-formed on newer C++ standards. To avoid this, we use
+   * std::abs on default types, but simply return the number on unsigned types.
+   */
+  template <typename Number>
+  std::enable_if_t<!std::is_unsigned<Number>::value,
+                   typename numbers::NumberTraits<Number>::real_type>
+  get_abs(const Number a)
+  {
+    return std::abs(a);
+  }
+
+  template <typename Number>
+  std::enable_if_t<std::is_unsigned<Number>::value, Number>
+  get_abs(const Number a)
+  {
+    return a;
+  }
+
+  /**
+   * Check if a vector is a deal.II vector.
+   */
+  template <typename VectorType>
+  constexpr bool is_dealii_vector =
+    std::is_same<VectorType,
+                 dealii::Vector<typename VectorType::value_type>>::value ||
+    std::is_same<VectorType,
+                 dealii::BlockVector<typename VectorType::value_type>>::value ||
+    std::is_same<
+      VectorType,
+      dealii::LinearAlgebra::Vector<typename VectorType::value_type>>::value ||
+    std::is_same<VectorType,
+                 dealii::LinearAlgebra::distributed::Vector<
+                   typename VectorType::value_type>>::value ||
+    std::is_same<VectorType,
+                 dealii::LinearAlgebra::distributed::BlockVector<
+                   typename VectorType::value_type>>::value;
+
+  /**
+   * Helper functions that call set_ghost_state() if the vector supports this
+   * operation.
+   */
+  template <typename T>
+  using set_ghost_state_t =
+    decltype(std::declval<T const>().set_ghost_state(std::declval<bool>()));
+
+  template <typename T>
+  constexpr bool has_set_ghost_state =
+    is_supported_operation<set_ghost_state_t, T>;
+
+  template <typename VectorType,
+            typename std::enable_if<has_set_ghost_state<VectorType>,
+                                    VectorType>::type * = nullptr>
+  void
+  set_ghost_state(VectorType &vector, const bool ghosted)
+  {
+    vector.set_ghost_state(ghosted);
+  }
+
+  template <typename VectorType,
+            typename std::enable_if<!has_set_ghost_state<VectorType>,
+                                    VectorType>::type * = nullptr>
+  void
+  set_ghost_state(VectorType &, const bool)
+  {
+    // serial vector: nothing to do
+  }
+
+  /**
+   * Helper function that sets the values on a cell, but also checks if the
+   * new values are similar to the old values.
+   */
+  template <int  dim,
+            int  spacedim,
+            bool lda,
+            class OutputVector,
+            typename number>
+  void
+  set_dof_values(const DoFCellAccessor<dim, spacedim, lda> &cell,
+                 const Vector<number> &                     local_values,
+                 OutputVector &                             values,
+                 const bool                                 perform_check)
+  {
+    (void)perform_check;
+
+#ifdef DEBUG
+    if (perform_check && is_dealii_vector<OutputVector>)
+      {
+        const bool old_ghost_state = values.has_ghost_elements();
+        set_ghost_state(values, true);
+
+        Vector<number> local_values_old(cell.get_fe().n_dofs_per_cell());
+        cell.get_dof_values(values, local_values_old);
+
+        for (unsigned int i = 0; i < cell.get_fe().n_dofs_per_cell(); ++i)
+          {
+            // a check consistent with the one in
+            // Utilities::MPI::Partitioner::import_from_ghosted_array_finish()
+            Assert(local_values_old[i] == number() ||
+                     get_abs(local_values_old[i] - local_values[i]) <=
+                       get_abs(local_values_old[i] + local_values[i]) *
+                         100000. *
+                         std::numeric_limits<typename numbers::NumberTraits<
+                           number>::real_type>::epsilon(),
+                   ExcNonMatchingElementsSetDofValuesByInterpolation<number>(
+                     local_values[i], local_values_old[i]));
+          }
+
+        set_ghost_state(values, old_ghost_state);
+      }
+#endif
+
+    cell.set_dof_values(local_values, values);
+  }
+
+
+  template <int  dim,
+            int  spacedim,
+            bool lda,
+            class OutputVector,
+            typename number>
+  void
+  process_by_interpolation(
+    const DoFCellAccessor<dim, spacedim, lda> &      cell,
+    const Vector<number> &                           local_values,
+    OutputVector &                                   values,
+    const unsigned int                               fe_index_,
+    const std::function<void(const DoFCellAccessor<dim, spacedim, lda> &cell,
+                             const Vector<number> &local_values,
+                             OutputVector &        values)> &processor)
+  {
+    const unsigned int fe_index =
+      (cell.get_dof_handler().has_hp_capabilities() == false &&
+       fe_index_ == DoFHandler<dim, spacedim>::invalid_fe_index) ?
+        DoFHandler<dim, spacedim>::default_fe_index :
+        fe_index_;
+
+    if (cell.is_active() && !cell.is_artificial())
+      {
+        if ((cell.get_dof_handler().has_hp_capabilities() == false) ||
+            // for hp-DoFHandlers, we need to require that on
+            // active cells, you either don't specify an fe_index,
+            // or that you specify the correct one
+            (fe_index == cell.active_fe_index()) ||
+            (fe_index == DoFHandler<dim, spacedim>::invalid_fe_index))
+          // simply set the values on this cell
+          processor(cell, local_values, values);
+        else
+          {
+            Assert(local_values.size() ==
+                     cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell(),
+                   ExcMessage("Incorrect size of local_values vector."));
+
+            FullMatrix<double> interpolation(
+              cell.get_fe().n_dofs_per_cell(),
+              cell.get_dof_handler().get_fe(fe_index).n_dofs_per_cell());
+
+            cell.get_fe().get_interpolation_matrix(
+              cell.get_dof_handler().get_fe(fe_index), interpolation);
+
+            // do the interpolation to the target space. for historical
+            // reasons, matrices are set to size 0x0 internally even if
+            // we reinit as 4x0, so we have to treat this case specially
+            Vector<number> tmp(cell.get_fe().n_dofs_per_cell());
+            if ((tmp.size() > 0) && (local_values.size() > 0))
+              interpolation.vmult(tmp, local_values);
+
+            // now set the dof values in the global vector
+            processor(cell, tmp, values);
+          }
+      }
+    else
+      // otherwise distribute them to the children
+      {
+        Assert((cell.get_dof_handler().has_hp_capabilities() == false) ||
+                 (fe_index != DoFHandler<dim, spacedim>::invalid_fe_index),
+               ExcMessage(
+                 "You cannot call this function on non-active cells "
+                 "of DoFHandler objects unless you provide an explicit "
+                 "finite element index because they do not have naturally "
+                 "associated finite element spaces associated: degrees "
+                 "of freedom are only distributed on active cells for which "
+                 "the active FE index has been set."));
+
+        const FiniteElement<dim, spacedim> &fe =
+          cell.get_dof_handler().get_fe(fe_index);
+        const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+
+        Assert(local_values.size() == dofs_per_cell,
+               (typename DoFCellAccessor<dim, spacedim, lda>::BaseClass::
+                  ExcVectorDoesNotMatch()));
+        Assert(values.size() == cell.get_dof_handler().n_dofs(),
+               (typename DoFCellAccessor<dim, spacedim, lda>::BaseClass::
+                  ExcVectorDoesNotMatch()));
+
+        Vector<number> tmp(dofs_per_cell);
+
+        for (unsigned int child = 0; child < cell.n_children(); ++child)
+          {
+            if (tmp.size() > 0)
+              fe.get_prolongation_matrix(child, cell.refinement_case())
+                .vmult(tmp, local_values);
+            process_by_interpolation(
+              *cell.child(child), tmp, values, fe_index, processor);
+          }
+      }
+  }
+
+} // namespace internal
+
+
 
 template <int dim, int spacedim, bool lda>
 template <class OutputVector, typename number>
@@ -48,84 +271,19 @@ void
 DoFCellAccessor<dim, spacedim, lda>::set_dof_values_by_interpolation(
   const Vector<number> &local_values,
   OutputVector &        values,
-  const unsigned int    fe_index_) const
+  const unsigned int    fe_index_,
+  const bool            perform_check) const
 {
-  const unsigned int fe_index =
-    (this->dof_handler->hp_capability_enabled == false &&
-     fe_index_ == DoFHandler<dim, spacedim>::invalid_fe_index) ?
-      DoFHandler<dim, spacedim>::default_fe_index :
-      fe_index_;
-
-  if (this->is_active() && !this->is_artificial())
-    {
-      if ((this->dof_handler->hp_capability_enabled == false) ||
-          // for hp-DoFHandlers, we need to require that on
-          // active cells, you either don't specify an fe_index,
-          // or that you specify the correct one
-          (fe_index == this->active_fe_index()) ||
-          (fe_index == DoFHandler<dim, spacedim>::invalid_fe_index))
-        // simply set the values on this cell
-        this->set_dof_values(local_values, values);
-      else
-        {
-          Assert(local_values.size() ==
-                   this->dof_handler->get_fe(fe_index).n_dofs_per_cell(),
-                 ExcMessage("Incorrect size of local_values vector."));
-
-          FullMatrix<double> interpolation(
-            this->get_fe().n_dofs_per_cell(),
-            this->dof_handler->get_fe(fe_index).n_dofs_per_cell());
-
-          this->get_fe().get_interpolation_matrix(
-            this->dof_handler->get_fe(fe_index), interpolation);
-
-          // do the interpolation to the target space. for historical
-          // reasons, matrices are set to size 0x0 internally even
-          // we reinit as 4x0, so we have to treat this case specially
-          Vector<number> tmp(this->get_fe().n_dofs_per_cell());
-          if ((tmp.size() > 0) && (local_values.size() > 0))
-            interpolation.vmult(tmp, local_values);
-
-          // now set the dof values in the global vector
-          this->set_dof_values(tmp, values);
-        }
-    }
-  else
-    // otherwise distribute them to the children
-    {
-      Assert((this->dof_handler->hp_capability_enabled == false) ||
-               (fe_index != DoFHandler<dim, spacedim>::invalid_fe_index),
-             ExcMessage(
-               "You cannot call this function on non-active cells "
-               "of DoFHandler objects unless you provide an explicit "
-               "finite element index because they do not have naturally "
-               "associated finite element spaces associated: degrees "
-               "of freedom are only distributed on active cells for which "
-               "the active FE index has been set."));
-
-      const FiniteElement<dim, spacedim> &fe =
-        this->get_dof_handler().get_fe(fe_index);
-      const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
-
-      Assert(this->dof_handler != nullptr,
-             typename BaseClass::ExcInvalidObject());
-      Assert(local_values.size() == dofs_per_cell,
-             typename BaseClass::ExcVectorDoesNotMatch());
-      Assert(values.size() == this->dof_handler->n_dofs(),
-             typename BaseClass::ExcVectorDoesNotMatch());
-
-      Vector<number> tmp(dofs_per_cell);
-
-      for (unsigned int child = 0; child < this->n_children(); ++child)
-        {
-          if (tmp.size() > 0)
-            fe.get_prolongation_matrix(child, this->refinement_case())
-              .vmult(tmp, local_values);
-          this->child(child)->set_dof_values_by_interpolation(tmp,
-                                                              values,
-                                                              fe_index);
-        }
-    }
+  internal::process_by_interpolation<dim, spacedim, lda, OutputVector, number>(
+    *this,
+    local_values,
+    values,
+    fe_index_,
+    [perform_check](const DoFCellAccessor<dim, spacedim, lda> &cell,
+                    const Vector<number> &                     local_values,
+                    OutputVector &                             values) {
+      internal::set_dof_values(cell, local_values, values, perform_check);
+    });
 }
 
 
index cc02b5b864600ec08d29407ebb6329027a5538fc..6bdc91315faf0d906f9f65b4f4d81231033caff2 100644 (file)
@@ -20,7 +20,8 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL)
     template void DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>::
       set_dof_values_by_interpolation(const Vector<VEC::value_type> &,
                                       VEC &,
-                                      const unsigned int fe_index) const;
+                                      const unsigned int fe_index,
+                                      const bool) const;
 
 #if deal_II_dimension != 3
 
@@ -28,15 +29,18 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; lda : BOOL)
     DoFCellAccessor<deal_II_dimension, deal_II_dimension + 1, lda>::
       set_dof_values_by_interpolation(const Vector<VEC::value_type> &,
                                       VEC &,
-                                      const unsigned int fe_index) const;
+                                      const unsigned int fe_index,
+                                      const bool) const;
 
 #endif
 
 #if deal_II_dimension == 3
 
     template void DoFCellAccessor<1, 3, lda>::set_dof_values_by_interpolation(
-      const Vector<VEC::value_type> &, VEC &, const unsigned int fe_index)
-      const;
+      const Vector<VEC::value_type> &,
+      VEC &,
+      const unsigned int fe_index,
+      const bool) const;
 
 #endif
   }
index a7175306cf5b35aa029f693c6f168faee692121f..685c486763b7d7f58b83f97d6a2cc6e88217d700 100644 (file)
@@ -192,7 +192,8 @@ SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
               in, (*pointerstruct->second.indices_ptr)[i]);
           cell->set_dof_values_by_interpolation(local_values,
                                                 out,
-                                                this_fe_index);
+                                                this_fe_index,
+                                                true);
         }
     }
 }
@@ -520,7 +521,8 @@ SolutionTransfer<dim, VectorType, spacedim>::interpolate(
 
                   cell->set_dof_values_by_interpolation(tmp,
                                                         all_out[j],
-                                                        old_fe_index);
+                                                        old_fe_index,
+                                                        true);
                 }
             }
           else if (valuesptr)
index 9dc5e2f8ae9494107a786a08723fb644f8b56e52..01e95e7471f8d78afda6a1b447dfb4f0f59fc5f9 100644 (file)
@@ -88,7 +88,6 @@ main()
   // get the interpolated solution
   // back
   Vector<double> tmp(dh.n_dofs());
-  tmp = 2;
   soltrans.interpolate(solution, tmp);
 
   deallog << "New values:" << std::endl;

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.