From: Pasquale Africa <pasqualeclaudio.africa@polimi.it>
Date: Tue, 22 Oct 2019 15:17:15 +0000 (+0000)
Subject: Make VectorTools::interpolate working in parallel + update documentation.
X-Git-Tag: v9.2.0-rc1~947^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=322d36bc0d33fafad3d4eb232c6459c3873b4912;p=dealii.git

Make VectorTools::interpolate working in parallel + update documentation.

Co-authored-by: Ivan Fumagalli <ivan.fumagalli@polimi.it>
---

diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 52ae5023a1..341e3429c5 100644
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -639,10 +639,11 @@ namespace VectorTools
 
   /**
    * Interpolate different finite element spaces. The interpolation of vector
-   * @p data_1 is executed from the FE space represented by @p dof_1 to the
-   * vector @p data_2 on FE space @p dof_2. The interpolation on each cell is
-   * represented by the matrix @p transfer. Curved boundaries are neglected so
-   * far.
+   * @p data_1 (which is assumed to be ghosted, see @ref GlossGhostedVector)
+   * is executed from the FE space represented by @p dof_1
+   * to the vector @p data_2 (which must be non-ghosted) on FE space @p dof_2.
+   * The interpolation on each cell is represented by the matrix @p transfer.
+   * Curved boundaries are neglected so far.
    *
    * Note that you may have to call <tt>hanging_nodes.distribute(data_2)</tt>
    * with the hanging nodes from space @p dof_2 afterwards, to make the result
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 35b085efd8..91eb05be4b 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -602,53 +602,60 @@ namespace VectorTools
     Vector<number> cell_data_1(dof_1.get_fe().dofs_per_cell);
     Vector<number> cell_data_2(dof_2.get_fe().dofs_per_cell);
 
-    std::vector<short unsigned int> touch_count(
-      dof_2.n_dofs(), 0); // TODO: check on datatype... kinda strange (UK)
+    // Store how many cells share each dof (unghosted).
+    OutVector touch_count;
+    touch_count.reinit(data_2);
+
     std::vector<types::global_dof_index> local_dof_indices(
       dof_2.get_fe().dofs_per_cell);
 
-    typename DoFHandler<dim, spacedim>::active_cell_iterator h =
+    typename DoFHandler<dim, spacedim>::active_cell_iterator cell_1 =
       dof_1.begin_active();
-    typename DoFHandler<dim, spacedim>::active_cell_iterator l =
+    typename DoFHandler<dim, spacedim>::active_cell_iterator cell_2 =
       dof_2.begin_active();
-    const typename DoFHandler<dim, spacedim>::cell_iterator endh = dof_1.end();
+    const typename DoFHandler<dim, spacedim>::cell_iterator end_1 = dof_1.end();
 
-    for (; h != endh; ++h, ++l)
+    for (; cell_1 != end_1; ++cell_1, ++cell_2)
       {
-        h->get_dof_values(data_1, cell_data_1);
-        transfer.vmult(cell_data_2, cell_data_1);
+        if (cell_1->is_locally_owned())
+          {
+            Assert(cell_2->is_locally_owned(), ExcInternalError());
 
-        l->get_dof_indices(local_dof_indices);
+            // Perform dof interpolation.
+            cell_1->get_dof_values(data_1, cell_data_1);
+            transfer.vmult(cell_data_2, cell_data_1);
 
-        // distribute cell vector
-        for (unsigned int j = 0; j < dof_2.get_fe().dofs_per_cell; ++j)
-          {
-            ::dealii::internal::ElementAccess<OutVector>::add(
-              cell_data_2(j), local_dof_indices[j], data_2);
-
-            // count, how often we have
-            // added to this dof
-            Assert(
-              touch_count[local_dof_indices[j]] <
-                std::numeric_limits<decltype(touch_count)::value_type>::max(),
-              ExcInternalError());
-            ++touch_count[local_dof_indices[j]];
+            cell_2->get_dof_indices(local_dof_indices);
+
+            // Distribute cell vector.
+            for (unsigned int j = 0; j < dof_2.get_fe().dofs_per_cell; ++j)
+              {
+                ::dealii::internal::ElementAccess<OutVector>::add(
+                  cell_data_2(j), local_dof_indices[j], data_2);
+
+                // Count cells that share each dof.
+                ::dealii::internal::ElementAccess<OutVector>::add(
+                  1, local_dof_indices[j], touch_count);
+              }
           }
       }
 
-    // compute the mean value of the
-    // sum which we have placed in each
-    // entry of the output vector
-    for (unsigned int i = 0; i < dof_2.n_dofs(); ++i)
+    // Collect information from other parallel processes.
+    data_2.compress(VectorOperation::add);
+    touch_count.compress(VectorOperation::add);
+
+    // Compute the mean value of the sum which has been placed in
+    // each entry of the output vector only at locally owned elements.
+    for (const auto &i : data_2.locally_owned_elements())
       {
         Assert(touch_count[i] != 0, ExcInternalError());
-        using value_type = typename OutVector::value_type;
-        const value_type val =
-          ::dealii::internal::ElementAccess<OutVector>::get(data_2, i);
 
         ::dealii::internal::ElementAccess<OutVector>::set(
-          val / static_cast<value_type>(touch_count[i]), i, data_2);
+          data_2[i] / touch_count[i], i, data_2);
       }
+
+    // Compress data_2 to set the proper values on all the processors.
+    data_2.compress(VectorOperation::insert);
   }