From: Daniel Arndt Date: Fri, 2 Sep 2016 21:17:14 +0000 (+0200) Subject: Fix implementation X-Git-Tag: v8.5.0-rc1~615^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73d55ae1f3314f2a6136d2fd34fdac699f6c39ef;p=dealii.git Fix implementation --- diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index b609e6c6ec..274d2daf51 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -714,28 +714,29 @@ namespace FETools * corresponding cell of `dof2` using the interpolation matrix of the finite * element spaces used on these cells and provided by the finite element * objects involved. This step is done using the FETools::interpolate() - * function. - It then performs a loop over all non-active cells of `dof2`. + * function. + * - It then performs a loop over all non-active cells of `dof2`. * If such a non-active cell has at least one active child, then we call the * children of this cell a "patch". We then interpolate from the children of * this patch to the patch, using the finite element space associated with * `dof2` and immediately interpolate back to the children. In essence, this * information throws away all information in the solution vector that lives - * on a scale smaller than the patch cell. - Since we traverse non-active - * cells from the coarsest to the finest levels, we may find patches that - * correspond to child cells of previously treated patches if the mesh had - * been refined adaptively (this cannot happen if the mesh has been refined - * globally because there the children of a patch are all active). We also - * perform the operation described above on these patches, but it is easy to - * see that on patches that are children of previously treated patches, the - * operation is now the identity operation (since it interpolates from the - * children of the current patch a function that had previously been - * interpolated to these children from an even coarser patch). Consequently, - * this does not alter the solution vector any more. + * on a scale smaller than the patch cell. + * - Since we traverse non-active cells from the coarsest to the finest + * levels, we may find patches that correspond to child cells of previously + * treated patches if the mesh had been refined adaptively (this cannot + * happen if the mesh has been refined globally because there the children + * of a patch are all active). We also perform the operation described above + * on these patches, but it is easy to see that on patches that are children + * of previously treated patches, the operation is now the identity operation + * (since it interpolates from the children of the current patch a function + * that had previously been interpolated to these children from an even coarser + * patch). Consequently, this does not alter the solution vector any more. * * The name of the function originates from the fact that it can be used to * construct a representation of a function of higher polynomial degree on a * once coarser mesh. For example, if you imagine that you start with a - * $Q_1$ function on globally refined mesh, and that @p dof2 is associated + * $Q_1$ function on a globally refined mesh, and that @p dof2 is associated * with a $Q_2$ element, then this function computes the equivalent of the * operator $I_{2h}^{(2)}$ interpolating the original piecewise linear * function onto a quadratic function on a once coarser mesh with mesh size @@ -781,6 +782,29 @@ namespace FETools const DoFHandler &dof2, const ConstraintMatrix &constraints, OutVector &z2); + + /* + * Same as above for external parallel VectorTypes + * on a parallel::distributed::Triangulation + */ + template + void extrapolate_parallel(const DoFHandler &dof1, + const VectorType &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + VectorType &u2); + + /* + * Same as above for deal.II parallel Vectors + * on a parallel::distributed::Triangulation + */ + template + void extrapolate_parallel(const DoFHandler &dof1, + const LinearAlgebra::distributed::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + LinearAlgebra::distributed::Vector &u2); + //@} /** * The numbering of the degrees of freedom in continuous finite elements is diff --git a/source/fe/fe_tools_interpolate.cc b/source/fe/fe_tools_interpolate.cc index c1752b8fd9..78aed6702c 100644 --- a/source/fe/fe_tools_interpolate.cc +++ b/source/fe/fe_tools_interpolate.cc @@ -44,12 +44,9 @@ #include -namespace -{ #include -} - #include +#include DEAL_II_NAMESPACE_OPEN @@ -91,11 +88,12 @@ namespace FETools Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + + const IndexSet u2_elements = u2.locally_owned_elements(); #ifdef DEBUG const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); const IndexSet &dof2_local_dofs = dof2.locally_owned_dofs(); const IndexSet u1_elements = u1.locally_owned_elements(); - const IndexSet u2_elements = u2.locally_owned_elements(); Assert(u1_elements == dof1_local_dofs, ExcMessage("The provided vector and DoF handler should have the same" " index sets.")); @@ -198,8 +196,13 @@ namespace FETools for (unsigned int i=0; iset_dof_values(u1_int_local, u1_interpolated); } - // if we work on a parallel PETSc vector - // we have to finish the work u1_interpolated.compress(VectorOperation::insert); } @@ -1000,6 +1001,10 @@ namespace FETools // a vector of all cells this process has // computed or received data std::vector available_cells; + + // stores the indices of dofs on more refined ghosted cells along + // with the maximum level + std::map dofs_on_refined_neighbors; }; template <> @@ -1016,9 +1021,9 @@ namespace FETools {} template - void extrapolate_parallel (const InVector &u2_relevant, - const DoFHandler<1,1> &dof2, - OutVector &u2) + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,1> &/*dof2*/, + OutVector &/*u2*/) {} }; @@ -1036,9 +1041,9 @@ namespace FETools {} template - void extrapolate_parallel (const InVector &u2_relevant, - const DoFHandler<1,2> &dof2, - OutVector &u2) + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,2> &/*dof2*/, + OutVector &/*u2*/) {} }; @@ -1056,9 +1061,9 @@ namespace FETools {} template - void extrapolate_parallel (const InVector &u2_relevant, - const DoFHandler<1,3> &dof2, - OutVector &u2) + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,3> &/*dof2*/, + OutVector &/*u2*/) {} }; @@ -1115,7 +1120,7 @@ namespace FETools &p4est_cell, dealii::internal::p4est::functions::quadrant_compare); - // this cell and none of it's children belongs to us + // if neither this cell nor one of it's children belongs to us, don't do anything if (idx == -1 && (dealii::internal::p4est::functions:: quadrant_overlaps_tree (const_cast::tree *>(&tree), &p4est_cell) @@ -1129,9 +1134,7 @@ namespace FETools { Assert (dealii_cell->has_children (), ExcInternalError ()); - // check if at least one - // child is locally owned - // on our process + // check if at least one child is locally owned on our process for (unsigned int child_n=0; child_nn_children(); ++child_n) if (dealii_cell->child(child_n)->active()) if (dealii_cell->child(child_n)->is_locally_owned()) @@ -1170,26 +1173,8 @@ namespace FETools u2); } - // traverse recursively over this tree - if (p4est_has_children) - { - typename dealii::internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - dealii::internal::p4est::init_quadrant_children (p4est_cell, - p4est_child); - for (unsigned int c=0; c::max_children_per_cell; ++c) - { - interpolate_recursively (forest, - tree, - tree_index, - dealii_cell->child(c), - p4est_child[c], - u1, - u2); - } - } } @@ -1261,7 +1246,7 @@ namespace FETools // this is a cell this process needs // data from another process - // check if this cell + // check if this cell is // available from other // computed patches CellData cell_data; @@ -1291,7 +1276,7 @@ namespace FETools else { // get the values from the present child, if necessary by - // interpolation itself either from its own children or + // interpolation either from its own children or // by interpolating from the finite element on an active // child to the finite element space requested here get_interpolated_dof_values (forest, @@ -1311,10 +1296,13 @@ namespace FETools // and add up or set them in the output vector for (unsigned int i=0; i &local_values, OutVector &u) { + const FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + if (!dealii_cell->has_children ()) { if (dealii_cell->is_locally_owned ()) { - // if this is one of our cells, - // set dof values in output vector - dealii_cell->set_dof_values (local_values, u); + std::vector indices(dofs_per_cell); + dealii_cell->get_dof_indices(indices); + for (unsigned int j=0; jdealii_cell->level())) + u(indices[j]) = local_values(j); + } } } else { - const FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - Assert (&fe != 0, ExcNotInitialized()); Assert (local_values.size() == dofs_per_cell, @@ -1446,7 +1441,7 @@ namespace FETools &p4est_cell, dealii::internal::p4est::functions::quadrant_compare); - // this cell and none of it's children belongs to us + // if neither this cell nor one of it's children belongs to us, don't do anything if (idx == -1 && (dealii::internal::p4est::functions:: quadrant_overlaps_tree (const_cast::tree *>(&tree), &p4est_cell) @@ -1539,8 +1534,8 @@ namespace FETools &p4est_child[c]) == false) { - // this is a cell this process needs - // data from another process + // this is a cell for which this process + // needs data from another process // so add the cell to the list add_new_need (forest, tree_index, @@ -1646,7 +1641,7 @@ namespace FETools &p4est_cell, dealii::internal::p4est::functions::quadrant_compare); - // this cell and none of it's children belongs to us + // if neither this cell nor one one of it's children belongs to us, don't do anything if (idx == -1 && (dealii::internal::p4est::functions:: quadrant_overlaps_tree (const_cast::tree *>(&tree), &p4est_cell) @@ -1859,31 +1854,31 @@ namespace FETools computed_cells, cells_to_send; - // compute all the cells needed - // from other processes + // Compute all the cells needed + // from other processes. compute_needs (dof2, cells_we_need); - // send the cells needed to there - // owners and receive a list other - // processes need from us + // Send the cells needed to there + // owners and receive a list of cells other + // processes need from us. send_cells (cells_we_need, received_needs); - // the list of received needs can contain - // some cells more than ones because different - // processes may need data from the same cell - // to compute data only ones, generate a vector - // with unique entries and distribute computed + // The list of received needs can contain + // some cells more than once because different + // processes may need data from the same cell. + // To compute data only once, generate a vector + // with unique entries and distribute the computed // data afterwards back to a vector with correct - // receivers to send the data back - // computing cell_data can cause some new cells - // needed for this ones - // if a cell is computed send it back to + // receivers. + // Computing cell_data can cause a need for + // data from some new cells. + // If a cell is computed, send it back to // their senders, maybe receive new needs and // compute again, do not wait that all cells - // are computed or all needs are collected, - // otherwise we can run into a deadlock if + // are computed or all needs are collected. + // Otherwise we can run into a deadlock if // a cell needed from another process, - // itself needs some data from us + // itself needs some data from us. unsigned int ready = 0; do { @@ -1953,43 +1948,127 @@ namespace FETools = (dynamic_cast*> (&dof2.get_tria())); - Assert (tr != 0, ExcMessage ("Exrapolate in parallel only works for parallel distributed triangulations!")); + Assert (tr != 0, + ExcMessage ("Extrapolate in parallel only works for parallel distributed triangulations!")); communicator = tr->get_communicator (); compute_all_non_local_data (dof2, u2_relevant); + // exclude dofs on more refined ghosted cells + const FiniteElement &fe = dof2.get_fe(); + const unsigned int dofs_per_face = fe.dofs_per_face; + if (dofs_per_face > 0) + { + const unsigned int dofs_per_cell = fe.dofs_per_cell; + std::vector indices (dofs_per_cell); + typename DoFHandler::active_cell_iterator + cell=dof2.begin_active(), + endc=dof2.end(); + for (; cell!=endc; ++cell) + if (cell->is_ghost()) + { + cell->get_dof_indices(indices); + for (unsigned int face=0; face::faces_per_cell; ++face) + if (!cell->at_boundary(face)) + { + const typename DoFHandler::cell_iterator neighbor = cell->neighbor(face); + if (neighbor->level() != cell->level()) + for (unsigned int i=0; ilevel(), dofs_on_refined_neighbors[index]) : + cell->level(); + dofs_on_refined_neighbors[index] = level; + } + } + } + } + // after all dof values on // not owned patch cells // are computed, start // the interpolation u2 = 0; - typename DoFHandler::cell_iterator - cell=dof2.begin(0), - endc=dof2.end(0); - for (; cell!=endc; ++cell) + struct WorkPackage + { + const typename dealii::internal::p4est::types::forest forest; + const typename dealii::internal::p4est::types::tree tree; + const typename dealii::internal::p4est::types::locidx tree_index; + const typename DoFHandler::cell_iterator dealii_cell; + const typename dealii::internal::p4est::types::quadrant p4est_cell; + + WorkPackage(const typename dealii::internal::p4est::types::forest &forest_, + const typename dealii::internal::p4est::types::tree &tree_, + const typename dealii::internal::p4est::types::locidx &tree_index_, + const typename DoFHandler::cell_iterator &dealii_cell_, + const typename dealii::internal::p4est::types::quadrant &p4est_cell_) + : + forest(forest_), + tree(tree_), + tree_index(tree_index_), + dealii_cell(dealii_cell_), + p4est_cell(p4est_cell_) + {} + }; + + std::queue queue; + { + typename DoFHandler::cell_iterator + cell=dof2.begin(0), + endc=dof2.end(0); + for (; cell!=endc; ++cell) + { + if (dealii::internal::p4est::tree_exists_locally (tr->parallel_forest, + tr->coarse_cell_to_p4est_tree_permutation[cell->index()]) + == false) + continue; + + typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; + const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; + typename dealii::internal::p4est::types::tree *tree = tr->init_tree(cell->index()); + + dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + + const WorkPackage data (*tr->parallel_forest, *tree,tree_index,cell,p4est_coarse_cell); + + queue.push(data); + } + } + + while (!queue.empty()) { - if (dealii::internal::p4est::tree_exists_locally (tr->parallel_forest, - tr->coarse_cell_to_p4est_tree_permutation[cell->index()]) - == false) - continue; + const WorkPackage &data = queue.front(); - typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; - const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; - typename dealii::internal::p4est::types::tree *tree = tr->init_tree(cell->index()); + const typename dealii::internal::p4est::types::forest &forest = data.forest; + const typename dealii::internal::p4est::types::tree &tree = data.tree; + const typename dealii::internal::p4est::types::locidx &tree_index= data.tree_index; + const typename DoFHandler::cell_iterator &dealii_cell = data.dealii_cell; + const typename dealii::internal::p4est::types::quadrant &p4est_cell = data.p4est_cell; - dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + interpolate_recursively (forest, tree, tree_index, dealii_cell, p4est_cell, u2_relevant, u2); - interpolate_recursively (*tr->parallel_forest, - *tree, - tree_index, - cell, - p4est_coarse_cell, - u2_relevant, - u2); + // traverse recursively over this tree + if (dealii_cell->has_children()) + { + Assert(dealii_cell->has_children(), ExcInternalError()); + typename dealii::internal::p4est::types::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; + + dealii::internal::p4est::init_quadrant_children (p4est_cell, + p4est_child); + + for (unsigned int c=0; c::max_children_per_cell; ++c) + queue.push(WorkPackage(forest, tree, tree_index, dealii_cell->child(c), p4est_child[c])); + } + queue.pop(); } + u2.compress(VectorOperation::insert); } @@ -2112,84 +2191,56 @@ namespace FETools - // special version for PETSc -#ifdef DEAL_II_WITH_PETSC - template - void extrapolate(const DoFHandler &dof1, - const PETScWrappers::MPI::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - PETScWrappers::MPI::Vector &u2) + template + void extrapolate_parallel(const DoFHandler &dof1, + const VectorType &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + VectorType &u2) { IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); IndexSet dof2_locally_relevant_dofs; DoFTools::extract_locally_relevant_dofs (dof2, dof2_locally_relevant_dofs); - PETScWrappers::MPI::Vector u3 (dof2_locally_owned_dofs, - u1.get_mpi_communicator ()); + VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ()); interpolate (dof1, u1, dof2, constraints2, u3); - PETScWrappers::MPI::Vector u3_relevant (dof2_locally_owned_dofs, - dof2_locally_relevant_dofs, - u1.get_mpi_communicator ()); + VectorType u3_relevant (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator ()); u3_relevant = u3; internal::extrapolate_parallel (u3_relevant, dof2, u2); - } -#endif - - - // special version for Trilinos -#ifdef DEAL_II_WITH_TRILINOS - template - void extrapolate(const DoFHandler &dof1, - const TrilinosWrappers::MPI::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - TrilinosWrappers::MPI::Vector &u2) - { - IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); - IndexSet dof2_locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dof2, - dof2_locally_relevant_dofs); - - TrilinosWrappers::MPI::Vector u3 (dof2_locally_owned_dofs, - u1.get_mpi_communicator ()); - interpolate (dof1, u1, dof2, constraints2, u3); - TrilinosWrappers::MPI::Vector u3_relevant (dof2_locally_relevant_dofs, - u1.get_mpi_communicator ()); - u3_relevant = u3; - - internal::extrapolate_parallel (u3_relevant, dof2, u2); + constraints2.distribute(u2); } -#endif - // special version for LinearAlgebra::distributed::Vector - template - void extrapolate(const DoFHandler &dof1, - const parallel::distributed::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - parallel::distributed::Vector &u2) + template + void extrapolate_parallel(const DoFHandler &dof1, + const LinearAlgebra::distributed::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + LinearAlgebra::distributed::Vector &u2) { IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); IndexSet dof2_locally_relevant_dofs; DoFTools::extract_locally_relevant_dofs (dof2, dof2_locally_relevant_dofs); - LinearAlgebra::distributed::Vector u3 (dof2_locally_owned_dofs, - dof2_locally_relevant_dofs, - u2.get_mpi_communicator ()); - + LinearAlgebra::distributed::Vector u3(dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator ()); interpolate (dof1, u1, dof2, constraints2, u3); - u3.update_ghost_values (); + u3.update_ghost_values(); internal::extrapolate_parallel (u3, dof2, u2); + + constraints2.distribute(u2); } + } // end of namespace FETools diff --git a/source/fe/fe_tools_interpolate.inst.in b/source/fe/fe_tools_interpolate.inst.in index 08b6654304..58933f1833 100644 --- a/source/fe/fe_tools_interpolate.inst.in +++ b/source/fe/fe_tools_interpolate.inst.in @@ -43,15 +43,18 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS void interpolate (const hp::DoFHandler &, const Vector &, const hp::DoFHandler &, Vector &); + template void interpolate (const hp::DoFHandler &, const Vector &, const hp::DoFHandler &, const ConstraintMatrix &, Vector &); + template void interpolate (const hp::DoFHandler &, const Vector &, const hp::DoFHandler &, Vector &); + template void interpolate (const hp::DoFHandler &, const Vector &, @@ -61,8 +64,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS \} } - - for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; VEC : SERIAL_VECTORS) { namespace FETools @@ -72,30 +73,36 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS void back_interpolate (const DoFHandler &, const VEC &, const FiniteElement &, VEC &); + template void back_interpolate (const DoFHandler &, const ConstraintMatrix &, const VEC &, const DoFHandler &, const ConstraintMatrix &, VEC &); + template void interpolation_difference (const DoFHandler &, const VEC &, const FiniteElement &, VEC &); + template void interpolation_difference (const DoFHandler &, const ConstraintMatrix &, const VEC &, const DoFHandler &, const ConstraintMatrix &, VEC &); + template void project_dg (const DoFHandler &, const VEC &, const DoFHandler &, VEC &); + template void extrapolate (const DoFHandler &, const VEC &, const DoFHandler &, VEC &); + template void extrapolate (const DoFHandler &, const VEC &, @@ -104,3 +111,49 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #endif \} } + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +{ + namespace FETools + \{ +#if deal_II_dimension == deal_II_space_dimension + +#ifdef DEAL_II_WITH_PETSC + template + void extrapolate_parallel + (const DoFHandler &, + const PETScWrappers::MPI::Vector &, + const DoFHandler &, + const ConstraintMatrix &, + PETScWrappers::MPI::Vector &); +#endif + +#ifdef DEAL_II_WITH_TRILINOS + template + void extrapolate_parallel + (const DoFHandler &, + const TrilinosWrappers::MPI::Vector &, + const DoFHandler &, + const ConstraintMatrix &, + TrilinosWrappers::MPI::Vector &); +#endif + +#endif + \} +} + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; scalar : REAL_SCALARS) +{ + namespace FETools + \{ +#if deal_II_dimension == deal_II_space_dimension + template + void extrapolate_parallel + (const DoFHandler &, + const LinearAlgebra::distributed::Vector &, + const DoFHandler &, + const ConstraintMatrix &, + LinearAlgebra::distributed::Vector &); +#endif + \} +}