From: msteigemann Date: Mon, 15 Jun 2015 13:59:08 +0000 (+0200) Subject: FETools::extrapolate: Working on the extrapolate algorithm on distributed parallel... X-Git-Tag: v8.5.0-rc1~615^2~14 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=af2ce9928fc2013bb3021159168fe5ea28efa408;p=dealii.git FETools::extrapolate: Working on the extrapolate algorithm on distributed parallel grids --- diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 0db12becd7..cde26ae7cc 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -967,6 +967,8 @@ namespace parallel mark_locally_active_vertices_on_level(const int level) const; template friend class dealii::internal::DoFHandler::Policy::ParallelDistributed; + + template friend class dealii::FETools::internal::ExtrapolateImplementation; }; diff --git a/source/fe/fe_tools_interpolate.cc b/source/fe/fe_tools_interpolate.cc index 170cd029d8..9485121a39 100644 --- a/source/fe/fe_tools_interpolate.cc +++ b/source/fe/fe_tools_interpolate.cc @@ -44,6 +44,22 @@ #include +#ifdef DEAL_II_WITH_P4EST +# include +# include +# include +# include +# include +# include + +# include +# include +# include +# include +# include +# include +#endif + #include @@ -499,6 +515,7 @@ namespace FETools } + template void back_interpolate(const DoFHandler &dof1, const ConstraintMatrix &constraints1, @@ -713,6 +730,7 @@ namespace FETools } + template void extrapolate(const DoFHandler &dof1, const InVector &u1, @@ -726,73 +744,1493 @@ namespace FETools - template - void extrapolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints, - OutVector &u2) + namespace internal { - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(&dof1.get_triangulation()==&dof2.get_triangulation(), ExcTriangulationMismatch()); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + namespace p4est + { + /** + * A structure whose explicit specializations contain pointers to the + * relevant p4est_* and p8est_* functions. Using this structure, for + * example by saying functions::quadrant_compare, we can write code + * in a dimension independent way, either calling p4est_quadrant_compare + * or p8est_quadrant_compare, depending on template argument. + */ + template struct functions; + + template <> + struct functions<2> + { + static + int (&quadrant_compare) (const void *v1, const void *v2); + + static + int (&quadrant_overlaps_tree) (dealii::internal::p4est::types<2>::tree *tree, + const dealii::internal::p4est::types<2>::quadrant *q); + + static + int (&comm_find_owner) (dealii::internal::p4est::types<2>::forest *p4est, + const dealii::internal::p4est::types<2>::locidx which_tree, + const dealii::internal::p4est::types<2>::quadrant *q, + const int guess); + }; + + int (&functions<2>::quadrant_compare) (const void *v1, const void *v2) + = p4est_quadrant_compare; + + int (&functions<2>::quadrant_overlaps_tree) (dealii::internal::p4est::types<2>::tree *tree, + const dealii::internal::p4est::types<2>::quadrant *q) + = p4est_quadrant_overlaps_tree; + + int (&functions<2>::comm_find_owner) (dealii::internal::p4est::types<2>::forest *p4est, + const dealii::internal::p4est::types<2>::locidx which_tree, + const dealii::internal::p4est::types<2>::quadrant *q, + const int guess) + = p4est_comm_find_owner; + + template <> + struct functions<3> + { + static + int (&quadrant_compare) (const void *v1, const void *v2); + + static + int (&quadrant_overlaps_tree) (dealii::internal::p4est::types<3>::tree *tree, + const dealii::internal::p4est::types<3>::quadrant *q); + + static + int (&comm_find_owner) (dealii::internal::p4est::types<3>::forest *p4est, + const dealii::internal::p4est::types<3>::locidx which_tree, + const dealii::internal::p4est::types<3>::quadrant *q, + const int guess); + }; + + int (&functions<3>::quadrant_compare) (const void *v1, const void *v2) + = p8est_quadrant_compare; + + int (&functions<3>::quadrant_overlaps_tree) (dealii::internal::p4est::types<3>::tree *tree, + const dealii::internal::p4est::types<3>::quadrant *q) + = p8est_quadrant_overlaps_tree; + + int (&functions<3>::comm_find_owner) (dealii::internal::p4est::types<3>::forest *p4est, + const dealii::internal::p4est::types<3>::locidx which_tree, + const dealii::internal::p4est::types<3>::quadrant *q, + const int guess) + = p8est_comm_find_owner; + + template + bool + tree_exists_locally (const typename dealii::internal::p4est::types::forest *parallel_forest, + const typename dealii::internal::p4est::types::topidx coarse_grid_cell) + { + Assert (coarse_grid_cell < parallel_forest->connectivity->num_trees, + ExcInternalError()); + return ((coarse_grid_cell >= parallel_forest->first_local_tree) + && + (coarse_grid_cell <= parallel_forest->last_local_tree)); + } + } + + /** + * Implementation of the @p extrapolate function + * on parallel distributed grids. + */ + template + class ExtrapolateImplementation + { + public: + + ExtrapolateImplementation (); + + ~ExtrapolateImplementation (); + + template + void extrapolate_parallel (const InVector &u2_relevant, + const DoFHandler &dof2, + OutVector &u2); + + private: + + /** + * A structure holding all data + * of cells needed from other processes + * for the extrapolate algorithm. + */ + struct CellData + { + CellData (); + + CellData (const unsigned int dofs_per_cell); + + Vector dof_values; + + unsigned int tree_index; + + typename dealii::internal::p4est::types::quadrant quadrant; + + int receiver; + + bool operator < (const CellData &rhs) const + { + if (p4est::functions::quadrant_compare (&quadrant, &rhs.quadrant) < 0) + return true; + + return false; + } - OutVector u3; - u3.reinit(u2); - interpolate(dof1, u1, dof2, constraints, u3); + unsigned int bytes_for_buffer () const + { + return (sizeof(unsigned int) + // dofs_per_cell + dof_values.size() * sizeof(double) + // dof_values + sizeof(unsigned int) + // tree_index + sizeof(typename dealii::internal::p4est::types::quadrant)); // quadrant + } + + void pack_data (std::vector &buffer) const + { + buffer.resize(bytes_for_buffer()); + + char *ptr = &buffer[0]; + + unsigned int n_dofs = dof_values.size (); + std::memcpy(ptr, &n_dofs, sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + std::memcpy(ptr,dof_values.begin(),n_dofs*sizeof(double)); + ptr += n_dofs*sizeof(double); + + std::memcpy(ptr,&tree_index,sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + std::memcpy(ptr,&quadrant,sizeof(typename dealii::internal::p4est::types::quadrant)); + ptr += sizeof(typename dealii::internal::p4est::types::quadrant); + + Assert (ptr == &buffer[0]+buffer.size(), + ExcInternalError()); + } + + void unpack_data (const std::vector &buffer) + { + const char *ptr = &buffer[0]; + unsigned int n_dofs; + memcpy(&n_dofs, ptr, sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + dof_values.reinit(n_dofs); + std::memcpy(dof_values.begin(),ptr,n_dofs * sizeof(double)); + ptr += n_dofs * sizeof(double); + + std::memcpy(&tree_index,ptr,sizeof(unsigned int)); + ptr += sizeof(unsigned int); + + std::memcpy(&quadrant,ptr,sizeof(typename dealii::internal::p4est::types::quadrant)); + ptr += sizeof(typename dealii::internal::p4est::types::quadrant); + + Assert (ptr == &buffer[0]+buffer.size(), + ExcInternalError()); + } + }; + + // Problem: The function extrapolates a polynomial + // function from a finer mesh of size $h$ to a polynmial + // function of higher degree but on a coarser mesh of + // size $2h$. Therefor the mesh has to consist of patches + // of four (2d) or eight (3d) cells and the function requires + // that the mesh is refined globally at least once. + // The algorithm starts on the coarsest level of the grid, + // loops over all cells and if a cell has at least one active child, + // dof values are set via + // cell->get_interpolated_dof_values(u_input, dof_values) + // cell->set_dof_values_by_interpolation(dof_values, u_output) + // both *_interpolation_* functions traverse recursively + // over all children down to the active cells + // and get/set dof values by interpolation. + // + // On distributed parallel grids the problem arises, that + // if a cell has at least one active child, there is no guarantee + // that all children of this cell belong to this process. + // There might be children which are owned by other processes + // and the algorithm needs to find and has to get the + // dof values from these processes and so on. + // + // Algorithm: + // 1) Loop over all coarse cells + // From each coarse cell traverse down the tree + // of refined cells and search for active children + // If there is an active child, check, if all + // other children down to the finest level are part + // of this process, if not, add the cell to the list + // of needs. + // 2) Send the list of needs to other processes + // Each process has a list of needs and after + // the loop over all coarse cells (all trees) + // is finished, this list can be send to other processes. + // 3) Compute the needs required by other processes + // While computing the values required from other + // processes there can arise new needs and they are + // stored in a list again. + // 4) Send the computed values and the list of new needs around + // + // This procedure has to be repeated until no process needs any + // new need and all needs are computed, but there are at most the + // "number of grid levels" rounds of sending/receiving cell data. + // + // Then each process has all data needed for extrapolation. + + // driver function sending/receiving all values from + // different processes + template + void compute_all_non_local_data (const DoFHandler &dof2, + const InVector &u); + + // traverse recursively over + // the cells of this tree and + // interpolate over patches which + // are part of this process + template + void interpolate_recursively (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, + const InVector &u1, + OutVector &u2); + + // get dof values for this + // cell by interpolation + // if a child is reached which + // is not part of this process + // a new need is created + template + void get_interpolated_dof_values (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, + const InVector &u, + Vector &interpolated_values, + std::vector &new_needs); + + // set dof values for this + // cell by interpolation + template + void set_dof_values_by_interpolation (const typename DoFHandler::cell_iterator &dealii_cell, + const typename dealii::internal::p4est::types::quadrant &p4est_cell, + const Vector &interpolated_values, + OutVector &u); + + // compute all cell_data + // needed from other processes + // to interpolate the part of + // this process + void compute_needs (const DoFHandler &dof2, + std::vector &new_needs); + + // traverse over the tree + // and look for patches this + // process has to work on + void traverse_tree_recursively (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, + std::vector &new_needs); + + // traverse recursively + // over a patch and look + // for cells needed from + // other processes for interpolation + void traverse_patch_recursively (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, + std::vector &new_needs); + + // compute dof values of all + // cells collected in cells_to_compute + // computed cells are deleted + // from cells_to_compute and + // stored in computed_cells + // during computation there can + // be new cells needed from + // other processes, these cells + // are stored in new_needs + template + void compute_cells (const DoFHandler &dof2, + const InVector &u, + std::vector &cells_to_compute, + std::vector &computed_cells, + std::vector &new_needs); + + // traverse over the tree + // and compute cells + template + void compute_cells_in_tree_recursively (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, + const InVector &u, + std::vector &cells_to_compute, + std::vector &computed_cells, + std::vector &new_needs); + + // send all cell_data in the vector + // cells_to_send to their receivers + // and receives a vector of cell_data + void send_cells (const std::vector &cells_to_send, + std::vector &received_cells) const; + + // add new cell_data to + // the ordered list new_needs + // uses cell_data_insert + void add_new_need (const typename dealii::internal::p4est::types::forest &forest, + 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, + std::vector &new_needs) const; + + // binary search in cells_list + // assume that cells_list + // is ordered + static int cell_data_search (const CellData &cell_data, + const std::vector &cells_list); + + // insert cell_data into a sorted + // vector cells_list at the correct + // position if cell_data + // not exists already in cells_list + static void cell_data_insert (const CellData &cell_data, + std::vector &cells_list); + + MPI_Comm communicator; + + // a vector of all cells this process has + // computed or received data + std::vector available_cells; + }; + + template + ExtrapolateImplementation:: + ExtrapolateImplementation () + {} + + + + template + ExtrapolateImplementation:: + ~ExtrapolateImplementation () + {} + + + + template + ExtrapolateImplementation:: + CellData::CellData () + : tree_index (0), + receiver (0) + {} + + + + template + ExtrapolateImplementation:: + CellData::CellData (const unsigned int dofs_per_cell) + : tree_index (0), + receiver (0) + { + dof_values.reinit (dofs_per_cell); + } + + + + template + template + void + ExtrapolateImplementation:: + interpolate_recursively (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, + const InVector &u1, + OutVector &u2) + { + // check if this cell exists in the local p4est + int idx = sc_array_bsearch(const_cast(&tree.quadrants), + &p4est_cell, + p4est::functions::quadrant_compare); + + // this cell and none of it's children belongs to us + if (idx == -1 && (p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_cell) + == false)) + return; + + bool p4est_has_children = (idx == -1); + + bool locally_owned_children = false; + if (p4est_has_children) + { + Assert (dealii_cell->has_children (), ExcInternalError ()); + + // 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()) + { + locally_owned_children=true; + break; + } + } + + if (locally_owned_children) + { + const FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + Vector interpolated_values(dofs_per_cell); + + std::vector new_needs; + get_interpolated_dof_values (forest, + tree, + tree_index, + dealii_cell, + p4est_cell, + u1, + interpolated_values, + new_needs); + + // at this point of + // the procedure no new + // needs should come up + Assert (new_needs.size () == 0, + ExcInternalError ()); + + set_dof_values_by_interpolation (dealii_cell, + p4est_cell, + interpolated_values, + u2); + } + + // traverse recursively over this tree + if (p4est_has_children) + { + typename dealii::internal::p4est::types::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; - const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; - Vector dof_values(dofs_per_cell); + dealii::internal::p4est::init_quadrant_children (p4est_cell, + p4est_child); - // make sure that each cell on the - // coarsest level is at least once - // refined. otherwise, we can't - // treat these cells and would - // generate a bogus result + 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); + } + } + } + + + + template + template + void + ExtrapolateImplementation:: + get_interpolated_dof_values (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, + const InVector &u, + Vector &interpolated_values, + std::vector &new_needs) { - typename DoFHandler::cell_iterator cell = dof2.begin(0), - endc = dof2.end(0); + if (!dealii_cell->has_children ()) + { + if (dealii_cell->is_locally_owned ()) + { + // if this is one of our cells, + // get dof values from input vector + dealii_cell->get_dof_values (u, interpolated_values); + } + else + { + add_new_need (forest, + tree_index, + dealii_cell, + p4est_cell, + new_needs); + } + } + 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 (interpolated_values.size() == dofs_per_cell, + ExcDimensionMismatch(interpolated_values.size(), dofs_per_cell)); + Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(), + ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs())); + + Vector tmp1(dofs_per_cell); + Vector tmp2(dofs_per_cell); + + interpolated_values = 0; + std::vector restriction_is_additive (dofs_per_cell); + for (unsigned int i=0; i::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; + + dealii::internal::p4est::init_quadrant_children (p4est_cell, + p4est_child); + + bool found_child = true; + for (unsigned int c=0; c::max_children_per_cell; ++c) + { + if (p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_child[c]) + == false) + { + // this is a cell this process needs + // data from another process + + // check if this cell + // available from other + // computed patches + CellData cell_data; + cell_data.quadrant = p4est_child[c]; + int pos = cell_data_search (cell_data, available_cells); + + if (pos == -1) + { + // data is not available + // create a new need + found_child = false; + + add_new_need (forest, + tree_index, + dealii_cell->child(c), + p4est_child[c], + new_needs); + } + else + { + Assert (available_cells[pos].dof_values.size() == dofs_per_cell, + ExcDimensionMismatch(available_cells[pos].dof_values.size(), dofs_per_cell)); + + tmp1 = available_cells[pos].dof_values; + } + } + else + { + // get the values from the present child, if necessary by + // interpolation itself 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, + tree, + tree_index, + dealii_cell->child(c), + p4est_child[c], + u, + tmp1, + new_needs); + } + + if (found_child) + { + // interpolate these to the mother cell + fe.get_restriction_matrix(c, dealii_cell->refinement_case()).vmult (tmp2, tmp1); + + // and add up or set them in the output vector + for (unsigned int i=0; i + template + void + ExtrapolateImplementation:: + set_dof_values_by_interpolation (const typename DoFHandler::cell_iterator &dealii_cell, + const typename dealii::internal::p4est::types::quadrant &p4est_cell, + const Vector &local_values, + OutVector &u) + { + 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); + } + } + 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, + ExcDimensionMismatch(local_values.size(), dofs_per_cell)); + Assert (u.size() == dealii_cell->get_dof_handler().n_dofs(), + ExcDimensionMismatch(u.size(), dealii_cell->get_dof_handler().n_dofs())); + + Vector tmp(dofs_per_cell); + + 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) + { + if (tmp.size() > 0) + fe.get_prolongation_matrix(c, dealii_cell->refinement_case()) + .vmult (tmp, local_values); + + set_dof_values_by_interpolation (dealii_cell->child(c), + p4est_child[c], + tmp, + u); + } + } + } + + + + template + void + ExtrapolateImplementation:: + compute_needs (const DoFHandler &dof2, + std::vector &new_needs) + { + const parallel::distributed::Triangulation< dim, spacedim > *tr + = (dynamic_cast*> + (&dof2.get_tria())); + + Assert (tr != 0, ExcInternalError()); + + typename DoFHandler::cell_iterator + cell=dof2.begin(0), + endc=dof2.end(0); for (; cell!=endc; ++cell) - Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce()); + { + if (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); + + // make sure that each cell on the + // coarsest level is at least once + // refined, otherwise, these cells + // can't be treated and would + // generate a bogus result + { + int idx = sc_array_bsearch(const_cast(&tree->quadrants), + &p4est_coarse_cell, + p4est::functions::quadrant_compare); + + AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ()); + } + + traverse_tree_recursively (*tr->parallel_forest, + *tree, + tree_index, + cell, + p4est_coarse_cell, + new_needs); + } } - // then traverse grid bottom up - for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), - endc=dof2.end(level); - for (; cell!=endc; ++cell) - if (!cell->active()) + + template + void + ExtrapolateImplementation:: + traverse_tree_recursively (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, + std::vector &new_needs) + { + // check if this cell exists in the local p4est + int idx = sc_array_bsearch(const_cast(&tree.quadrants), + &p4est_cell, + p4est::functions::quadrant_compare); + + // this cell and none of it's children belongs to us + if (idx == -1 && (p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_cell) + == false)) + return; + + bool p4est_has_children = (idx == -1); + + // this cell is part of a patch + // this process has to interpolate on + // if there is at least one locally + // owned child + bool locally_owned_children = false; + if (p4est_has_children) + { + Assert (dealii_cell->has_children (), ExcInternalError ()); + + 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()) + { + locally_owned_children=true; + break; + } + } + + // traverse the patch recursively + // to find cells needed from + // other processes + if (locally_owned_children) + { + traverse_patch_recursively (forest, + tree, + tree_index, + dealii_cell, + p4est_cell, + new_needs); + } + + // 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) { - // check whether this - // cell has active - // children - bool active_children=false; - for (unsigned int child_n=0; child_nn_children(); ++child_n) - if (cell->child(child_n)->active()) - { - active_children=true; - break; - } - - // if there are active - // children, the we have - // to work on this - // cell. get the data - // from the one vector - // and set it on the - // other - if (active_children) + traverse_tree_recursively (forest, + tree, + tree_index, + dealii_cell->child(c), + p4est_child[c], + new_needs); + } + } + } + + + + template + void + ExtrapolateImplementation:: + traverse_patch_recursively (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, + std::vector &new_needs) + { + 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) + { + // check if this child + // is locally available + // in the p4est + if (p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_child[c]) + == false) + { + // this is a cell this process needs + // data from another process + // so add the cell to the list + add_new_need (forest, + tree_index, + dealii_cell->child(c), + p4est_child[c], + new_needs); + } + else { - cell->get_interpolated_dof_values(u3, dof_values); - cell->set_dof_values_by_interpolation(dof_values, u2); + // at least some part of + // the tree rooted in this + // child is locally available + traverse_patch_recursively (forest, + tree, + tree_index, + dealii_cell->child(c), + p4est_child[c], + new_needs); } } + } + } + + + + + template + template + void + ExtrapolateImplementation:: + compute_cells (const DoFHandler &dof2, + const InVector &u, + std::vector &cells_to_compute, + std::vector &computed_cells, + std::vector &new_needs) + { + const parallel::distributed::Triangulation< dim, spacedim > *tr + = (dynamic_cast*> + (&dof2.get_tria())); + + Assert (tr != 0, ExcInternalError()); + + // collect in a set all trees this + // process has to compute cells on + std::set trees; + for (typename std::vector::const_iterator it=cells_to_compute.begin(); it!=cells_to_compute.end(); ++it) + trees.insert (it->tree_index); + + typename DoFHandler::cell_iterator + cell=dof2.begin(0), + endc=dof2.end(0); + for (; cell!=endc; ++cell) + { + // check if this is a tree this process has to + // work on and that this tree is in the p4est + const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; + + if ((trees.find (tree_index) == trees.end ()) + || + (p4est::tree_exists_locally (tr->parallel_forest, + tree_index) + == false)) + continue; + + typename dealii::internal::p4est::types::quadrant p4est_coarse_cell; + typename dealii::internal::p4est::types::tree *tree = tr->init_tree(cell->index()); + + dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + + compute_cells_in_tree_recursively (*tr->parallel_forest, + *tree, + tree_index, + cell, + p4est_coarse_cell, + u, + cells_to_compute, + computed_cells, + new_needs); + } + } + + + + template + template + void + ExtrapolateImplementation:: + compute_cells_in_tree_recursively (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, + const InVector &u, + std::vector &cells_to_compute, + std::vector &computed_cells, + std::vector &new_needs) + { + if (cells_to_compute.size () == 0) + return; + + // check if this cell exists in the local p4est + int idx = sc_array_bsearch(const_cast(&tree.quadrants), + &p4est_cell, + p4est::functions::quadrant_compare); + + // this cell and none of it's children belongs to us + if (idx == -1 && (p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_cell) + == false)) + return; + + bool p4est_has_children = (idx == -1); + + // check if this quadrant is in the list + CellData cell_data; + cell_data.quadrant = p4est_cell; + int pos = cell_data_search (cell_data, cells_to_compute); + if (pos != -1) + { + std::vector tmp; + // compute dof values + get_interpolated_dof_values (forest, + tree, + tree_index, + dealii_cell, + p4est_cell, + u, + cells_to_compute[pos].dof_values, + tmp); + // there is no new cell_data + // this process needs for computing this cell, + // store cell_data in the list of + // computed cells and erase this cell + // from the list of cells to compute + if (tmp.size () == 0) + { + cell_data_insert (cells_to_compute[pos], computed_cells); + cells_to_compute.erase (cells_to_compute.begin () + pos); + } + else + { + for (unsigned int i=0; i::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) + { + compute_cells_in_tree_recursively (forest, + tree, + tree_index, + dealii_cell->child(c), + p4est_child[c], + u, + cells_to_compute, + computed_cells, + new_needs); + } + } + } + + + + template + void + ExtrapolateImplementation:: + send_cells (const std::vector &cells_to_send, + std::vector &received_cells) const + { + std::vector > sendbuffers (cells_to_send.size()); + std::vector >::iterator buffer = sendbuffers.begin(); + std::vector requests (cells_to_send.size()); + std::vector destinations; + + // send data + unsigned int idx=0; + for (typename std::vector::const_iterator it=cells_to_send.begin(); + it!=cells_to_send.end(); + ++it, ++idx) + { + destinations.push_back (it->receiver); + + it->pack_data (*buffer); + MPI_Isend (&(*buffer)[0], buffer->size(), + MPI_BYTE, + it->receiver, + 123, + communicator, + &requests[idx]); + } + + Assert(destinations.size()==cells_to_send.size(), ExcInternalError()); + + std::vector senders + = Utilities::MPI::compute_point_to_point_communication_pattern(communicator, destinations); + + // receive data + std::vector receive; + CellData cell_data; + for (unsigned int index=0; index 0) + MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE); + + // finally sort the list of cells + std::sort (received_cells.begin (), received_cells.end ()); + } + + + + template + void + ExtrapolateImplementation:: + add_new_need (const typename dealii::internal::p4est::types::forest &forest, + 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, + std::vector &new_needs) const + { + const FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + Assert (&fe != 0, ExcNotInitialized()); + + CellData cell_data (dofs_per_cell); + cell_data.quadrant = p4est_cell; + cell_data.tree_index = tree_index; + cell_data.receiver = p4est::functions:: + comm_find_owner (const_cast::forest *> (&forest), + tree_index, + &p4est_cell, + dealii_cell->level_subdomain_id ()); + + cell_data_insert (cell_data, new_needs); + } + + + + template + int + ExtrapolateImplementation:: + cell_data_search (const CellData &cell_data, + const std::vector &cells_list) + { + typename std::vector::const_iterator bound + = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data); + + if ((bound != cells_list.end ()) + && + !(cell_data < *bound)) + return (int)(bound - cells_list.begin ()); + + return (-1); + } + + + + template + void + ExtrapolateImplementation:: + cell_data_insert (const CellData &cell_data, + std::vector &cells_list) + { + typename std::vector::iterator bound + = std::lower_bound (cells_list.begin(), cells_list.end(), cell_data); + + if ((bound == cells_list.end ()) + || + (cell_data < *bound)) + cells_list.insert (bound, 1, cell_data); + } + + + + template + template + void + ExtrapolateImplementation:: + compute_all_non_local_data (const DoFHandler &dof2, + const InVector &u) + { + std::vector cells_we_need, + cells_to_compute, + received_cells, + received_needs, + new_needs, + computed_cells, + cells_to_send; + + // 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_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 + // 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 + // 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 + // a cell needed from another process, + // itself needs some data from us + unsigned int ready = 0; + do + { + for (unsigned int i=0; i::const_iterator comp=computed_cells.begin (); + comp != computed_cells.end (); + ++comp) + { + // store computed cells + cell_data_insert (*comp, available_cells); + + // and generate a vector + // of computed cells with + // correct receivers + // then delete this received + // need from the list + typename std::vector::iterator recv=std::begin (received_needs); + while (recv != std::end (received_needs)) + { + if (dealii::internal::p4est::quadrant_is_equal(recv->quadrant, comp->quadrant)) + { + recv->dof_values = comp->dof_values; + cells_to_send.push_back (*recv); + received_needs.erase (recv); + recv = std::begin (received_needs); + } + else + ++recv; + } + } + + send_cells (cells_to_send, received_cells); + + // strore received cell_data + for (typename std::vector::const_iterator recv=received_cells.begin (); + recv != received_cells.end (); + ++recv) + { + cell_data_insert (*recv, available_cells); + } + + // finally send and receive new + // needs and start a new round + send_cells (new_needs, received_needs); + } + while (ready != 0); + } + + + + template + template + void + ExtrapolateImplementation:: + extrapolate_parallel (const InVector &u2_relevant, + const DoFHandler &dof2, + OutVector &u2) + { + const parallel::distributed::Triangulation< dim, spacedim > *tr + = (dynamic_cast*> + (&dof2.get_tria())); + + Assert (tr != 0, ExcMessage ("Exrapolate in parallel only works for parallel distributed triangulations!")); + + communicator = tr->get_communicator (); + + compute_all_non_local_data (dof2, u2_relevant); + + // 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) + { + if (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); + + interpolate_recursively (*tr->parallel_forest, + *tree, + tree_index, + cell, + p4est_coarse_cell, + u2_relevant, + u2); + } + + u2.compress(VectorOperation::insert); + } + + + + namespace + { + template + void extrapolate_serial(const InVector &u3, + const DoFHandler &dof2, + OutVector &u2) + { + const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; + Vector dof_values(dofs_per_cell); + + // then traverse grid bottom up + for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), + endc=dof2.end(level); + + for (; cell!=endc; ++cell) + if (!cell->active()) + { + // check whether this + // cell has active + // children + bool active_children=false; + for (unsigned int child_n=0; child_nn_children(); ++child_n) + if (cell->child(child_n)->active()) + { + active_children=true; + break; + } + + // if there are active + // children, this process + // has to work on this + // cell. get the data + // from the one vector + // and set it on the + // other + if (active_children) + { + cell->get_interpolated_dof_values(u3, dof_values); + cell->set_dof_values_by_interpolation(dof_values, u2); + } + } + } + } + + + + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) + { + // make sure that each cell on the + // coarsest level is at least once + // refined, otherwise, these cells + // can't be treated and would + // generate a bogus result + { + typename DoFHandler::cell_iterator cell = dof2.begin(0), + endc = dof2.end(0); + for (; cell!=endc; ++cell) + Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce()); + } + + OutVector u3; + u3.reinit(u2); + interpolate(dof1, u1, dof2, constraints, u3); + + extrapolate_serial (u3, dof2, u2); + } + + + + template + void extrapolate_parallel(const InVector &u2_relevant, + const DoFHandler &dof2, + OutVector &u2) + { + if (dynamic_cast*>(&dof2.get_tria()) != 0) + { + internal::ExtrapolateImplementation implementation; + implementation.extrapolate_parallel (u2_relevant, dof2, u2); + } + else + { + extrapolate_serial (u2_relevant, dof2, u2); + } } + + + // 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) + { + 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 ()); + interpolate (dof1, u1, dof2, constraints2, u3); + PETScWrappers::MPI::Vector u3_relevant (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator ()); + u3_relevant = u3; + + 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; + + extrapolate_parallel (u3_relevant, dof2, u2); + } +#endif + + + // special version for parallel::distributed::Vector + template + void extrapolate(const DoFHandler &dof1, + const parallel::distributed::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + parallel::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); + + parallel::distributed::Vector u3 (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u2.get_mpi_communicator ()); + + interpolate (dof1, u1, dof2, constraints2, u3); + u3.update_ghost_values (); + + extrapolate_parallel (u3, dof2, u2); + } + } + } + + + + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) + { + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + + internal::extrapolate (dof1, u1, dof2, constraints, u2); + // Apply hanging node constraints. constraints.distribute(u2); } diff --git a/source/fe/fe_tools_interpolate.inst.in b/source/fe/fe_tools_interpolate.inst.in index 08b6654304..af93c4cc92 100644 --- a/source/fe/fe_tools_interpolate.inst.in +++ b/source/fe/fe_tools_interpolate.inst.in @@ -92,6 +92,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS void project_dg (const DoFHandler &, const VEC &, const DoFHandler &, VEC &); +# if deal_II_dimension > 1 template void extrapolate (const DoFHandler &, const VEC &, @@ -101,6 +102,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (const DoFHandler &, const VEC &, const DoFHandler &, const ConstraintMatrix &, VEC &); +# endif #endif \} }