From: Daniel Arndt Date: Tue, 27 Sep 2016 15:33:34 +0000 (+0200) Subject: Add dealii comment blob and split extrapolate into an own file X-Git-Tag: v8.5.0-rc1~615^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e3a9962a22dfac87f3d513cd545051dc2692562;p=dealii.git Add dealii comment blob and split extrapolate into an own file --- diff --git a/include/deal.II/distributed/p4est_wrappers.h b/include/deal.II/distributed/p4est_wrappers.h index 1bbad8f114..67181ec79e 100644 --- a/include/deal.II/distributed/p4est_wrappers.h +++ b/include/deal.II/distributed/p4est_wrappers.h @@ -1,5 +1,20 @@ -#ifndef P4EST_WRAPPERS_H -#define P4EST_WRAPPERS_H +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__p4est_wrappers_h +#define dealii__p4est_wrappers_h #include @@ -531,4 +546,4 @@ DEAL_II_NAMESPACE_CLOSE #endif //DEAL_II_WITH_P4EST -#endif //P4EST_WRAPPERS_H +#endif //dealii__p4est_wrappers_h diff --git a/source/distributed/p4est_wrappers.inst.in b/source/distributed/p4est_wrappers.inst.in index ec5547c9e4..129cc271aa 100644 --- a/source/distributed/p4est_wrappers.inst.in +++ b/source/distributed/p4est_wrappers.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2010 - 2014 by the deal.II authors +// Copyright (C) 2016 - 2016 by the deal.II authors // // This file is part of the deal.II library. // diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index e71ac66a0e..1c744274b4 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -48,6 +48,7 @@ SET(_src fe_enriched.cc fe_tools.cc fe_tools_interpolate.cc + fe_tools_extrapolate.cc fe_trace.cc fe_values.cc fe_values_inst2.cc @@ -91,6 +92,7 @@ SET(_inst fe_enriched.inst.in fe_tools.inst.in fe_tools_interpolate.inst.in + fe_tools_extrapolate.inst.in fe_trace.inst.in fe_values.decl.1.inst.in fe_values.decl.2.inst.in diff --git a/source/fe/fe_tools_extrapolate.cc b/source/fe/fe_tools_extrapolate.cc new file mode 100755 index 0000000000..6c82be1bfb --- /dev/null +++ b/source/fe/fe_tools_extrapolate.cc @@ -0,0 +1,1567 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace FETools +{ + + namespace internal + { +#ifdef DEAL_II_WITH_P4EST + // 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 to + // set dofs recursively on cells of arbitrary level + 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_) + {} + }; + + + // 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 (dealii::internal::p4est::functions::quadrant_compare (&quadrant, &rhs.quadrant) < 0) + return true; + + return false; + } + + 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; + + // stores the indices of dofs on more refined ghosted cells along + // with the maximum level + std::map dofs_on_refined_neighbors; + }; + + template <> + class ExtrapolateImplementation<1,1> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + ~ExtrapolateImplementation () + {} + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,1> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + template <> + class ExtrapolateImplementation<1,2> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + ~ExtrapolateImplementation () + {} + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,2> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + template <> + class ExtrapolateImplementation<1,3> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + ~ExtrapolateImplementation () + {} + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,3> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + + + 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, + dealii::internal::p4est::functions::quadrant_compare); + + // 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) + == 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); + } + + + + } + + + + 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) + { + 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 (dealii::internal::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 is + // 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 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) + { + 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 ()) + { + 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 + { + 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) + { + 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); + + // 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, + dealii::internal::p4est::functions::quadrant_compare); + + AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ()); + } + + traverse_tree_recursively (*tr->parallel_forest, + *tree, + tree_index, + cell, + p4est_coarse_cell, + new_needs); + } + } + + + + 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, + dealii::internal::p4est::functions::quadrant_compare); + + // 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) + == 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) + { + 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 (dealii::internal::p4est::functions:: + quadrant_overlaps_tree (const_cast::tree *>(&tree), + &p4est_child[c]) + == false) + { + // 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, + dealii_cell->child(c), + p4est_child[c], + new_needs); + } + else + { + // 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 ()) + || + (dealii::internal::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, + dealii::internal::p4est::functions::quadrant_compare); + + // 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) + == 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 = dealii::internal::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 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 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. + // 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 + // 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=received_needs.begin(); + while (recv != received_needs.end()) + { + 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 = received_needs.begin(); + } + 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 ("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; + + 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()) + { + const WorkPackage &data = queue.front(); + + 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; + + interpolate_recursively (forest, tree, tree_index, dealii_cell, p4est_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); + } +#endif //DEAL_II_WITH_P4EST + + 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); + } + } + } + } + + + +#ifdef DEAL_II_WITH_P4EST + template + void extrapolate_parallel(const InVector &u2_relevant, + const DoFHandler &dof2, + OutVector &u2) + { + if (dynamic_cast*>(&dof2.get_tria()) != 0) + { + ExtrapolateImplementation implementation; + implementation.extrapolate_parallel (u2_relevant, dof2, u2); + } + else + { + extrapolate_serial (u2_relevant, dof2, u2); + } + } +#endif //DEAL_II_WITH_P4EST + } + } + + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + OutVector &u2) + { + ConstraintMatrix dummy; + dummy.close(); + extrapolate(dof1, u1, dof2, dummy, 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_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())); + + // 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); + + internal::extrapolate_serial (u3, dof2, u2); + + constraints.distribute(u2); + } + + +#ifdef DEAL_II_WITH_P4EST + 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); + + VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ()); + interpolate (dof1, u1, dof2, constraints2, u3); + 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); + + constraints2.distribute(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, + u1.get_mpi_communicator ()); + interpolate (dof1, u1, dof2, constraints2, u3); + u3.update_ghost_values(); + + internal::extrapolate_parallel (u3, dof2, u2); + + constraints2.distribute(u2); + } +#endif //DEAL_II_WITH_P4EST + +} // end of namespace FETools + + +/*-------------- Explicit Instantiations -------------------------------*/ +#include "fe_tools_extrapolate.inst" + + +/*---------------------------- fe_tools.cc ---------------------------*/ + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_tools_extrapolate.inst.in b/source/fe/fe_tools_extrapolate.inst.in new file mode 100755 index 0000000000..dbcb9375d3 --- /dev/null +++ b/source/fe/fe_tools_extrapolate.inst.in @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; VEC : SERIAL_VECTORS) +{ + namespace FETools + \{ +#if deal_II_dimension == deal_II_space_dimension + template + void extrapolate + (const DoFHandler &, const VEC &, + const DoFHandler &, VEC &); + + template + void extrapolate + (const DoFHandler &, const VEC &, + const DoFHandler &, const ConstraintMatrix &, + VEC &); +#endif + \} +} + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#ifdef DEAL_II_WITH_P4EST + 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 + \} +#endif //DEAL_II_WITH_P4EST +} + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; scalar : REAL_SCALARS) +{ +#ifdef DEAL_II_WITH_P4EST + 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 + \} +#endif //DEAL_II_WITH_P4EST +} + diff --git a/source/fe/fe_tools_interpolate.cc b/source/fe/fe_tools_interpolate.cc old mode 100644 new mode 100755 index 0b1d5cff97..b82ece6a35 --- a/source/fe/fe_tools_interpolate.cc +++ b/source/fe/fe_tools_interpolate.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2000 - 2015 by the deal.II authors +// Copyright (C) 2000 - 2016 by the deal.II authors // // This file is part of the deal.II library. // @@ -44,9 +44,7 @@ #include -#include #include -#include DEAL_II_NAMESPACE_OPEN @@ -718,1527 +716,6 @@ namespace FETools ++cell2; } } - - - namespace internal - { -#ifdef DEAL_II_WITH_P4EST - // 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 to - // set dofs recursively on cells of arbitrary level - 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_) - {} - }; - - - // 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 (dealii::internal::p4est::functions::quadrant_compare (&quadrant, &rhs.quadrant) < 0) - return true; - - return false; - } - - 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; - - // stores the indices of dofs on more refined ghosted cells along - // with the maximum level - std::map dofs_on_refined_neighbors; - }; - - template <> - class ExtrapolateImplementation<1,1> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - ~ExtrapolateImplementation () - {} - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,1> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - template <> - class ExtrapolateImplementation<1,2> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - ~ExtrapolateImplementation () - {} - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,2> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - template <> - class ExtrapolateImplementation<1,3> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - ~ExtrapolateImplementation () - {} - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,3> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - - - 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, - dealii::internal::p4est::functions::quadrant_compare); - - // 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) - == 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); - } - - - - } - - - - 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) - { - 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 (dealii::internal::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 is - // 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 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) - { - 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 ()) - { - 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 - { - 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) - { - 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); - - // 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, - dealii::internal::p4est::functions::quadrant_compare); - - AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ()); - } - - traverse_tree_recursively (*tr->parallel_forest, - *tree, - tree_index, - cell, - p4est_coarse_cell, - new_needs); - } - } - - - - 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, - dealii::internal::p4est::functions::quadrant_compare); - - // 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) - == 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) - { - 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 (dealii::internal::p4est::functions:: - quadrant_overlaps_tree (const_cast::tree *>(&tree), - &p4est_child[c]) - == false) - { - // 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, - dealii_cell->child(c), - p4est_child[c], - new_needs); - } - else - { - // 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 ()) - || - (dealii::internal::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, - dealii::internal::p4est::functions::quadrant_compare); - - // 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) - == 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 = dealii::internal::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 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 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. - // 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 - // 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=received_needs.begin(); - while (recv != received_needs.end()) - { - 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 = received_needs.begin(); - } - 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 ("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; - - 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()) - { - const WorkPackage &data = queue.front(); - - 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; - - interpolate_recursively (forest, tree, tree_index, dealii_cell, p4est_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); - } -#endif //DEAL_II_WITH_P4EST - - 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); - } - } - } - } - - - -#ifdef DEAL_II_WITH_P4EST - template - void extrapolate_parallel(const InVector &u2_relevant, - const DoFHandler &dof2, - OutVector &u2) - { - if (dynamic_cast*>(&dof2.get_tria()) != 0) - { - ExtrapolateImplementation implementation; - implementation.extrapolate_parallel (u2_relevant, dof2, u2); - } - else - { - extrapolate_serial (u2_relevant, dof2, u2); - } - } -#endif //DEAL_II_WITH_P4EST - } - } - - template - void extrapolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - OutVector &u2) - { - ConstraintMatrix dummy; - dummy.close(); - extrapolate(dof1, u1, dof2, dummy, 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_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())); - - // 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); - - internal::extrapolate_serial (u3, dof2, u2); - - constraints.distribute(u2); - } - - -#ifdef DEAL_II_WITH_P4EST - 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); - - VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ()); - interpolate (dof1, u1, dof2, constraints2, u3); - 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); - - constraints2.distribute(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, - u1.get_mpi_communicator ()); - interpolate (dof1, u1, dof2, constraints2, u3); - u3.update_ghost_values(); - - internal::extrapolate_parallel (u3, dof2, u2); - - constraints2.distribute(u2); - } -#endif //DEAL_II_WITH_P4EST - } // end of namespace FETools diff --git a/source/fe/fe_tools_interpolate.inst.in b/source/fe/fe_tools_interpolate.inst.in old mode 100644 new mode 100755 index 24fcd7dd13..d051fdfbe6 --- a/source/fe/fe_tools_interpolate.inst.in +++ b/source/fe/fe_tools_interpolate.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 1998 - 2014 by the deal.II authors +// Copyright (C) 1998 - 2016 by the deal.II authors // // This file is part of the deal.II library. // @@ -97,67 +97,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS 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 &, - const DoFHandler &, const ConstraintMatrix &, - VEC &); -#endif - \} -} - -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) -{ -#ifdef DEAL_II_WITH_P4EST - 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 - \} -#endif //DEAL_II_WITH_P4EST -} - -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; scalar : REAL_SCALARS) -{ -#ifdef DEAL_II_WITH_P4EST - 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 \} -#endif //DEAL_II_WITH_P4EST } diff --git a/tests/mpi/fe_tools_extrapolate_03.mpirun=10.output b/tests/mpi/fe_tools_extrapolate_03.mpirun=10.output new file mode 100644 index 0000000000..f406bbaa3d --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_03.mpirun=10.output @@ -0,0 +1,115 @@ + +DEAL:0::Checking Q1 against Q1 in 2d: +DEAL:0::23.00000000 3.944933460 1.000000000 +DEAL:0::23.00000000 3.944933460 1.000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking Q1 against Q1 in 3d: +DEAL:0::162.3750000 12.71625436 1.500000000 +DEAL:0::162.3750000 12.71625436 1.500000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking Q1 against Q2 in 2d: +DEAL:0::9.125000000 1.682283604 0.5000000000 +DEAL:0::26.88281250 2.700093063 0.5000000000 +DEAL:0::global_error_before: 0.01934481330 +DEAL:0::global_error_after: 0.001208379635 +DEAL:0::OK +DEAL:0::Checking Q1 against Q2 in 3d: +DEAL:0::66.28125000 5.415995551 0.7500000000 +DEAL:0::342.1987305 11.49092805 0.7500000000 +DEAL:0::global_error_before: 0.03152756518 +DEAL:0::global_error_after: 0.001607487444 +DEAL:0::OK +DEAL:0::Checking Q2 against Q1 in 2d: +DEAL:0::26.84375000 2.699035020 0.5000000000 +DEAL:0::10.25000000 1.831793861 0.5000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0.07284222390 +DEAL:0::OK +DEAL:0::Checking Q2 against Q1 in 3d: +DEAL:0::341.5859375 11.47797042 0.7500000000 +DEAL:0::75.09375000 5.992915740 0.7500000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0.1190491594 +DEAL:0::OK +DEAL:0::Checking Q2 against Q2 in 2d: +DEAL:0::26.84375000 2.699035020 0.5000000000 +DEAL:0::26.84375000 2.699035020 0.5000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking Q2 against Q2 in 3d: +DEAL:0::341.5859375 11.47797042 0.7500000000 +DEAL:0::341.5859375 11.47797042 0.7500000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking DGQ0 against DGQ0 in 2d: +DEAL:0::56.00000000 10.58300524 2.000000000 +DEAL:0::56.00000000 10.58300524 2.000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking DGQ0 against DGQ0 in 3d: +DEAL:0::276.0000000 28.77498914 3.000000000 +DEAL:0::276.0000000 28.77498914 3.000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking DGQ0 against DGQ1 in 2d: +DEAL:0::14.00000000 2.828427125 0.8750000000 +DEAL:0::56.00000000 5.678908346 0.9375000000 +DEAL:0::global_error_before: 0.09199750902 +DEAL:0::global_error_after: 0.04599875451 +DEAL:0::OK +DEAL:0::Checking DGQ0 against DGQ1 in 3d: +DEAL:0::72.50000000 7.854139036 1.312500000 +DEAL:0::580.6250000 22.31268054 1.406250000 +DEAL:0::global_error_before: 0.1220351512 +DEAL:0::global_error_after: 0.06040494959 +DEAL:0::OK +DEAL:0::Checking DGQ1 against DGQ0 in 2d: +DEAL:0::56.00000000 5.830951895 1.000000000 +DEAL:0::14.00000000 2.738612788 0.7500000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0.1839950180 +DEAL:0::OK +DEAL:0::Checking DGQ1 against DGQ0 in 3d: +DEAL:0::580.0000000 22.78157150 1.500000000 +DEAL:0::73.00000000 7.697402159 1.125000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0.2440703024 +DEAL:0::OK +DEAL:0::Checking DGQ1 against DGQ1 in 2d: +DEAL:0::56.00000000 5.830951895 1.000000000 +DEAL:0::56.00000000 5.830951895 1.000000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK +DEAL:0::Checking DGQ1 against DGQ1 in 3d: +DEAL:0::580.0000000 22.78157150 1.500000000 +DEAL:0::580.0000000 22.78157150 1.500000000 +DEAL:0::global_error_before: 0 +DEAL:0::global_error_after: 0 +DEAL:0::OK + + + + + + + + + + + + + + + + + + diff --git a/tests/mpi/fe_tools_extrapolate_common.h b/tests/mpi/fe_tools_extrapolate_common.h old mode 100644 new mode 100755 index 353fea4ac3..5cd255fbf3 --- a/tests/mpi/fe_tools_extrapolate_common.h +++ b/tests/mpi/fe_tools_extrapolate_common.h @@ -73,7 +73,8 @@ parallel::distributed::Triangulation *make_tria () cell = tria->begin_active(); cell->set_refine_flag(); ++cell; - cell->set_refine_flag(); + if (cell != tria->end()) + cell->set_refine_flag(); tria->execute_coarsening_and_refinement (); } @@ -98,8 +99,6 @@ template void output_vector (const VectorType &v, const std::string &output_name, const DoFHandler &dof_handler) { - //v.print(deallog.get_file_stream()); - DataOut data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (v, output_name, DataOut::type_dof_data);