From: Pasquale Africa 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 --- 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 hanging_nodes.distribute(data_2) * 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 cell_data_1(dof_1.get_fe().dofs_per_cell); Vector cell_data_2(dof_2.get_fe().dofs_per_cell); - std::vector 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 local_dof_indices( dof_2.get_fe().dofs_per_cell); - typename DoFHandler::active_cell_iterator h = + typename DoFHandler::active_cell_iterator cell_1 = dof_1.begin_active(); - typename DoFHandler::active_cell_iterator l = + typename DoFHandler::active_cell_iterator cell_2 = dof_2.begin_active(); - const typename DoFHandler::cell_iterator endh = dof_1.end(); + const typename DoFHandler::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::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::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::add( + cell_data_2(j), local_dof_indices[j], data_2); + + // Count cells that share each dof. + ::dealii::internal::ElementAccess::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::get(data_2, i); ::dealii::internal::ElementAccess::set( - val / static_cast(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); }