From: Daniel Arndt Date: Thu, 13 Apr 2017 14:22:27 +0000 (+0200) Subject: Fix FETools::extrapolate for block vectors X-Git-Tag: v9.0.0-rc1~1690^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=01bc19fd8f74639a8b13ed9289ffebcd5219886b;p=dealii.git Fix FETools::extrapolate for block vectors --- diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h new file mode 100644 index 0000000000..a3671bfe5b --- /dev/null +++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h @@ -0,0 +1,1634 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2017 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__fe_tools_extrapolate_templates_H +#define dealii__fe_tools_extrapolate_templates_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace FETools +{ + + namespace internal + { +#ifndef DEAL_II_WITH_P4EST + // Dummy implementation in case p4est is not available. + template + class ExtrapolateImplementation + { + public: + + ExtrapolateImplementation () + { + Assert(false, ExcNotImplemented()); + }; + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler &/*dof2*/, + OutVector &/*u2*/) + { + Assert(false, ExcNotImplemented()) + } + }; +#else + // Implementation of the @p extrapolate function + // on parallel distributed grids. + template + class ExtrapolateImplementation + { + public: + + ExtrapolateImplementation (); + + template + void extrapolate_parallel (const InVector &u2_relevant, + const DoFHandler &dof2, + OutVector &u2); + + private: + // A shortcut for the type of the OutVector + typedef typename OutVector::value_type value_type; + + // 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(value_type) + // 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(value_type)); + ptr += n_dofs*sizeof(value_type); + + 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(value_type)); + ptr += n_dofs * sizeof(value_type); + + 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 + 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 + static 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); + + // 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; + + // counts the send/receive round we are in + unsigned int round; + }; + + template + class ExtrapolateImplementation<1,1,OutVector> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,1> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + template + class ExtrapolateImplementation<1,2,OutVector> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,2> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + template + class ExtrapolateImplementation<1,3,OutVector> + { + public: + + ExtrapolateImplementation () + { + AssertThrow(false, ExcNotImplemented()) + } + + template + void extrapolate_parallel (const InVector &/*u2_relevant*/, + const DoFHandler<1,3> &/*dof2*/, + OutVector &/*u2*/) + {} + }; + + + + template + ExtrapolateImplementation:: + ExtrapolateImplementation () + : round(0) + {} + + + + 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 (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 + 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 (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 != nullptr, 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 != nullptr, 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); + const int ierr = MPI_Isend (&(*buffer)[0], buffer->size(), + MPI_BYTE, + it->receiver, + round, + communicator, + &requests[idx]); + AssertThrowMPI(ierr); + } + + 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) + { + const int ierr = MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE); + AssertThrowMPI(ierr); + } + + // 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 FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + 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; + + // reset the round count we will use in send_cells + round = 0; + + // Compute all the cells needed from other processes. + compute_needs (dof2, cells_we_need); + + // Send the cells needed to their 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; + } + } + + // increase the round counter, such that we are sure to only send + // and receive data from the correct call + ++round; + + send_cells (cells_to_send, received_cells); + + // store received cell_data + for (typename std::vector::const_iterator recv=received_cells.begin (); + recv != received_cells.end (); + ++recv) + { + cell_data_insert (*recv, available_cells); + } + + // increase the round counter, such that we are sure to only send + // and receive data from the correct call + ++round; + + // 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 != nullptr, + 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 + struct BlockTypeHelper + { + using type = VectorType; + }; + + template + struct BlockTypeHelper::value>::type> + { + using type = typename VectorType::BlockType; + }; + + template + using BlockType = typename BlockTypeHelper::type; + + template + void reinit_distributed(const DH &dh, + VectorType &vector) + { + vector.reinit(dh.n_dofs()); + } + +#ifdef DEAL_II_WITH_PETSC + template + void reinit_distributed(const DoFHandler &dh, + PETScWrappers::MPI::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert (parallel_tria != nullptr, ExcNotImplemented()); + + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); + } +#endif //DEAL_II_WITH_PETSC + +#ifdef DEAL_II_WITH_TRILINOS + template + void reinit_distributed(const DoFHandler &dh, + TrilinosWrappers::MPI::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert (parallel_tria != nullptr, ExcNotImplemented()); + + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); + } +#endif //DEAL_II_WITH_TRILINOS + + template + void reinit_distributed(const DoFHandler &dh, + LinearAlgebra::distributed::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert (parallel_tria != nullptr, ExcNotImplemented()); + + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); + } + + + + template + void reinit_ghosted(const DH &/*dh*/, + VectorType &/*vector*/) + { + Assert(false, ExcNotImplemented()); + } + +#ifdef DEAL_II_WITH_PETSC + template + void reinit_ghosted(const DoFHandler &dh, + PETScWrappers::MPI::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert(parallel_tria != nullptr, ExcNotImplemented()); + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); + vector.reinit (locally_owned_dofs, locally_relevant_dofs, + parallel_tria->get_communicator()); + } +#endif //DEAL_II_WITH_PETSC + +#ifdef DEAL_II_WITH_TRILINOS + template + void reinit_ghosted(const DoFHandler &dh, + TrilinosWrappers::MPI::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert (parallel_tria != nullptr, ExcNotImplemented()); + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); + vector.reinit (locally_owned_dofs, locally_relevant_dofs, + parallel_tria->get_communicator()); + } +#endif //DEAL_II_WITH_TRILINOS + + template + void reinit_ghosted(const DoFHandler &dh, + LinearAlgebra::distributed::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast*>(&dh.get_tria()); + Assert (parallel_tria != nullptr, ExcNotImplemented()); + const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); + vector.reinit (locally_owned_dofs, locally_relevant_dofs, + parallel_tria->get_communicator()); + } + + + + template + void extrapolate_serial(const InVector &u3, + const DoFHandler &dof2, + OutVector &u2) + { + const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; + Vector dof_values(dofs_per_cell); + + // then traverse grid bottom up + for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), + endc=dof2.end(level); + + for (; cell!=endc; ++cell) + if (!cell->active()) + { + // check whether this + // cell has active + // children + bool active_children=false; + for (unsigned int child_n=0; child_nn_children(); ++child_n) + if (cell->child(child_n)->active()) + { + active_children=true; + break; + } + + // if there are active + // children, this process + // has to work on this + // cell. get the data + // from the one vector + // and set it on the + // other + if (active_children) + { + cell->get_interpolated_dof_values(u3, dof_values); + cell->set_dof_values_by_interpolation(dof_values, u2); + } + } + } + } + } + } + + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + 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() || cell->is_artificial(), ExcGridNotRefinedAtLeastOnce()); + } + + + internal::BlockType u3; + internal::reinit_distributed(dof2, u3); + if (dynamic_cast*>(&dof2.get_tria()) != nullptr) + { + interpolate(dof1, u1, dof2, constraints, u3); + + internal::BlockType u3_relevant; + internal::reinit_ghosted(dof2, u3_relevant); + u3_relevant = u3; + + internal::ExtrapolateImplementation implementation; + implementation.extrapolate_parallel (u3_relevant, dof2, u2); + } + else + { + interpolate(dof1, u1, dof2, constraints, u3); + + internal::extrapolate_serial (u3, dof2, u2); + } + + constraints.distribute(u2); + } +} // end of namespace FETools + +DEAL_II_NAMESPACE_CLOSE + +/*---------------------------- fe_tools_extrapolate_templates.h ---------------------------*/ +/* end of #ifndef dealii__fe_tools_extrapolate_templates_H */ +#endif diff --git a/include/deal.II/fe/fe_tools_interpolate.templates.h b/include/deal.II/fe/fe_tools_interpolate.templates.h new file mode 100644 index 0000000000..fc987af8d0 --- /dev/null +++ b/include/deal.II/fe/fe_tools_interpolate.templates.h @@ -0,0 +1,655 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2017 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__fe_tools_interpolate_templates_H +#define dealii__fe_tools_interpolate_templates_H + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace FETools +{ + template class DoFHandlerType1, + template class DoFHandlerType2, + class InVector, class OutVector> + void + interpolate(const DoFHandlerType1 &dof1, + const InVector &u1, + const DoFHandlerType2 &dof2, + OutVector &u2) + { + ConstraintMatrix dummy; + dummy.close(); + interpolate(dof1, u1, dof2, dummy, u2); + } + + + + template class DoFHandlerType1, + template class DoFHandlerType2, + class InVector, class OutVector> + void + interpolate (const DoFHandlerType1 &dof1, + const InVector &u1, + const DoFHandlerType2 &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) + { + 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())); + + + const IndexSet u2_elements = u2.locally_owned_elements(); +#ifdef DEBUG + const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); + const IndexSet &dof2_local_dofs = dof2.locally_owned_dofs(); + const IndexSet u1_elements = u1.locally_owned_elements(); + Assert(u1_elements == dof1_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); + Assert(u2_elements == dof2_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); +#endif + + // allocate vectors at maximal + // size. will be reinited in inner + // cell, but Vector makes sure that + // this does not lead to + // reallocation of memory + Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); + Vector u2_local(DoFTools::max_dofs_per_cell(dof2)); + + // have a map for interpolation + // matrices. shared_ptr make sure + // that memory is released again + std::map *, + std::map *, + std::shared_ptr > > > + interpolation_matrices; + + typename DoFHandlerType1::active_cell_iterator cell1 = dof1.begin_active(), + endc1 = dof1.end(); + typename DoFHandlerType2::active_cell_iterator cell2 = dof2.begin_active(), + endc2 = dof2.end(); + (void)endc2; + + std::vector dofs; + dofs.reserve (DoFTools::max_dofs_per_cell (dof2)); + + u2 = 0; + OutVector touch_count(u2); + touch_count = 0; + + // for distributed triangulations, + // we can only interpolate u1 on + // a cell, which this processor owns, + // so we have to know the subdomain_id + const types::subdomain_id subdomain_id = + dof1.get_triangulation().locally_owned_subdomain(); + + for (; cell1!=endc1; ++cell1, ++cell2) + if ((cell1->subdomain_id() == subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id)) + { + Assert(cell1->get_fe().n_components() == cell2->get_fe().n_components(), + ExcDimensionMismatch (cell1->get_fe().n_components(), + cell2->get_fe().n_components())); + + // for continuous elements on + // grids with hanging nodes we + // need hanging node + // constraints. Consequently, + // if there are no constraints + // then hanging nodes are not + // allowed. + const bool hanging_nodes_not_allowed + = ((cell2->get_fe().dofs_per_vertex != 0) && + (constraints.n_constraints() == 0)); + + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell1->at_boundary(face) || + cell1->neighbor(face)->level() == cell1->level(), + ExcHangingNodesNotAllowed(0)); + + + const unsigned int dofs_per_cell1 = cell1->get_fe().dofs_per_cell; + const unsigned int dofs_per_cell2 = cell2->get_fe().dofs_per_cell; + u1_local.reinit (dofs_per_cell1); + u2_local.reinit (dofs_per_cell2); + + // check if interpolation + // matrix for this particular + // pair of elements is already + // there + if (interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()].get() == nullptr) + { + std::shared_ptr > + interpolation_matrix (new FullMatrix (dofs_per_cell2, + dofs_per_cell1)); + interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] + = interpolation_matrix; + + get_interpolation_matrix(cell1->get_fe(), + cell2->get_fe(), + *interpolation_matrix); + } + + cell1->get_dof_values(u1, u1_local); + interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] + ->vmult(u2_local, u1_local); + + dofs.resize (dofs_per_cell2); + cell2->get_dof_indices(dofs); + + for (unsigned int i=0; i(touch_count(i)) != typename OutVector::value_type(0), + ExcInternalError()); + u2(i) /= touch_count(i); + } + + // finish the work on parallel vectors + u2.compress(VectorOperation::insert); + // Apply hanging node constraints. + constraints.distribute(u2); + } + + + + template class DoFHandlerType, + class InVector, class OutVector, int spacedim> + void + back_interpolate(const DoFHandlerType &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated) + { + Assert(dof1.get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); + Assert(u1.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + +#ifdef DEBUG + const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); + const IndexSet u1_elements = u1.locally_owned_elements(); + const IndexSet u1_interpolated_elements = u1_interpolated.locally_owned_elements(); + Assert(u1_elements == dof1_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); + Assert(u1_interpolated_elements == dof1_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); +#endif + + Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); + Vector u1_int_local(DoFTools::max_dofs_per_cell(dof1)); + + const types::subdomain_id subdomain_id = + dof1.get_triangulation().locally_owned_subdomain(); + + typename DoFHandlerType::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + + // map from possible fe objects in + // dof1 to the back_interpolation + // matrices + std::map *, + std::shared_ptr > > interpolation_matrices; + + for (; cell!=endc; ++cell) + if ((cell->subdomain_id() == subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id)) + { + // For continuous elements on + // grids with hanging nodes we + // need hanging node + // constraints. Consequently, + // when the elements are + // continuous no hanging node + // constraints are allowed. + const bool hanging_nodes_not_allowed= + (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell; + + // make sure back_interpolation + // matrix is available + if (interpolation_matrices[&cell->get_fe()] == nullptr) + { + interpolation_matrices[&cell->get_fe()] = + std::shared_ptr > + (new FullMatrix(dofs_per_cell1, dofs_per_cell1)); + get_back_interpolation_matrix(cell->get_fe(), fe2, + *interpolation_matrices[&cell->get_fe()]); + } + + u1_local.reinit (dofs_per_cell1); + u1_int_local.reinit (dofs_per_cell1); + + cell->get_dof_values(u1, u1_local); + interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local); + cell->set_dof_values(u1_int_local, u1_interpolated); + }; + + // if we work on a parallel vector, we have to finish the work + u1_interpolated.compress(VectorOperation::insert); + } + + + + namespace internal + { + namespace + { + template + void back_interpolate (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + InVector &u1_interpolated) + { + Vector u2(dof2.n_dofs()); + interpolate(dof1, u1, dof2, constraints2, u2); + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } + + // special version for PETSc +#ifdef DEAL_II_WITH_PETSC + template + void back_interpolate (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const PETScWrappers::MPI::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + PETScWrappers::MPI::Vector &u1_interpolated) + { + // if u1 is a parallel distributed PETSc vector, we create a + // vector u2 with based on the sets of locally owned and relevant + // dofs of dof2 + IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); + IndexSet dof2_locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof2, + dof2_locally_relevant_dofs); + + PETScWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs, + u1.get_mpi_communicator()); + interpolate(dof1, u1, dof2, constraints2, u2_out); + PETScWrappers::MPI::Vector u2 (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator()); + u2 = u2_out; + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } +#endif + + // special version for Trilinos +#ifdef DEAL_II_WITH_TRILINOS + template + void back_interpolate (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const TrilinosWrappers::MPI::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + TrilinosWrappers::MPI::Vector &u1_interpolated) + { + // if u1 is a parallel distributed Trilinos vector, we create a + // vector u2 with based on the sets of locally owned and relevant + // dofs of dof2 + IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); + IndexSet dof2_locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs (dof2, + dof2_locally_relevant_dofs); + + TrilinosWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs, + u1.get_mpi_communicator()); + interpolate(dof1, u1, dof2, constraints2, u2_out); + TrilinosWrappers::MPI::Vector u2 (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator()); + u2 = u2_out; + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } +#endif + + // special version for LinearAlgebra::distributed::Vector + template + void back_interpolate (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const LinearAlgebra::distributed::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + LinearAlgebra::distributed::Vector &u1_interpolated) + { + 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 + u2 (dof2_locally_owned_dofs, + dof2_locally_relevant_dofs, + u1.get_mpi_communicator()); + + interpolate(dof1, u1, dof2, constraints2, u2); + u2.update_ghost_values (); + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } + } + } + + + + template + void back_interpolate(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + OutVector &u1_interpolated) + { + // For discontinuous elements without constraints take the simpler version + // of the back_interpolate function. + if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 + && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) + back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated); + else + { + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + + // For continuous elements first interpolate to dof2, taking into + // account constraints2, and then interpolate back to dof1 taking into + // account constraints1 + internal::back_interpolate(dof1, constraints1, u1, dof2, constraints2, + u1_interpolated); + } + } + + + + template + void interpolation_difference (const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_difference) + { + Assert(dof1.get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_difference.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs())); + +#ifdef DEBUG + const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); + const IndexSet u1_elements = u1.locally_owned_elements(); + const IndexSet u1_difference_elements = u1_difference.locally_owned_elements(); + Assert(u1_elements == dof1_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); + Assert(u1_difference_elements == dof1_local_dofs, + ExcMessage("The provided vector and DoF handler should have the same" + " index sets.")); +#endif + + // For continuous elements on grids + // with hanging nodes we need + // hanging node + // constraints. Consequently, when + // the elements are continuous no + // hanging node constraints are + // allowed. + const bool hanging_nodes_not_allowed= + (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell; + + Vector u1_local(dofs_per_cell); + Vector u1_diff_local(dofs_per_cell); + + const types::subdomain_id subdomain_id = + dof1.get_triangulation().locally_owned_subdomain(); + + FullMatrix difference_matrix(dofs_per_cell, dofs_per_cell); + get_interpolation_difference_matrix(dof1.get_fe(), fe2, + difference_matrix); + + typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + + for (; cell!=endc; ++cell) + if ((cell->subdomain_id() == subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id)) + { + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + cell->get_dof_values(u1, u1_local); + difference_matrix.vmult(u1_diff_local, u1_local); + cell->set_dof_values(u1_diff_local, u1_difference); + } + + // if we work on a parallel PETSc vector + // we have to finish the work and + // update ghost values + u1_difference.compress(VectorOperation::insert); + } + + + + namespace internal + { + namespace + { + template + void interpolation_difference (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + OutVector &u1_difference) + { + back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); + u1_difference.sadd(-1., 1., u1); + } + + // special version for Trilinos +#ifdef DEAL_II_WITH_TRILINOS + template + void interpolation_difference (const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const TrilinosWrappers::MPI::Vector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + TrilinosWrappers::MPI::Vector &u1_difference) + { + back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); + + // Trilinos vectors with and without ghost entries are very different + // and we cannot use the sadd function directly, so we have to create + // a completely distributed vector first and copy the local entries + // from the vector with ghost entries + TrilinosWrappers::MPI::Vector u1_completely_distributed (u1_difference.vector_partitioner ()); + + u1_completely_distributed = u1; + + u1_difference.sadd(-1, u1_completely_distributed); + } +#endif + } + } + + + + template + void interpolation_difference(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + OutVector &u1_difference) + { + // For discontinuous elements + // without constraints take the + // cheaper version of the + // interpolation_difference function. + if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 + && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) + interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference); + else + { + internal::interpolation_difference(dof1, constraints1, u1, dof2, constraints2, u1_difference); + } + } + + + + template + void project_dg(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + OutVector &u2) + { + Assert(&dof1.get_triangulation()==&dof2.get_triangulation(), ExcTriangulationMismatch()); + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + + typename DoFHandler::active_cell_iterator cell1 = dof1.begin_active(); + typename DoFHandler::active_cell_iterator cell2 = dof2.begin_active(); + typename DoFHandler::active_cell_iterator end = dof2.end(); + + const unsigned int n1 = dof1.get_fe().dofs_per_cell; + const unsigned int n2 = dof2.get_fe().dofs_per_cell; + + Vector u1_local(n1); + Vector u2_local(n2); + std::vector dofs(n2); + + FullMatrix matrix(n2,n1); + get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix); + + u2 = 0; + while (cell2 != end) + { + cell1->get_dof_values(u1, u1_local); + matrix.vmult(u2_local, u1_local); + cell2->get_dof_indices(dofs); + for (unsigned int i=0; i -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include +#include DEAL_II_NAMESPACE_OPEN -namespace FETools -{ - - namespace internal - { -#ifndef DEAL_II_WITH_P4EST - // Dummy implementation in case p4est is not available. - template - class ExtrapolateImplementation - { - public: - - ExtrapolateImplementation () - { - Assert(false, ExcNotImplemented()); - }; - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler &/*dof2*/, - OutVector &/*u2*/) - { - Assert(false, ExcNotImplemented()) - } - }; -#else - // Implementation of the @p extrapolate function - // on parallel distributed grids. - template - class ExtrapolateImplementation - { - public: - - ExtrapolateImplementation (); - - template - void extrapolate_parallel (const InVector &u2_relevant, - const DoFHandler &dof2, - OutVector &u2); - - private: - // A shortcut for the type of the OutVector - typedef typename OutVector::value_type value_type; - - // 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(value_type) + // 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(value_type)); - ptr += n_dofs*sizeof(value_type); - - 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(value_type)); - ptr += n_dofs * sizeof(value_type); - - 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 - 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 - static 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); - - // 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; - - // counts the send/receive round we are in - unsigned int round; - }; - - template - class ExtrapolateImplementation<1,1,OutVector> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,1> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - template - class ExtrapolateImplementation<1,2,OutVector> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,2> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - template - class ExtrapolateImplementation<1,3,OutVector> - { - public: - - ExtrapolateImplementation () - { - AssertThrow(false, ExcNotImplemented()) - } - - template - void extrapolate_parallel (const InVector &/*u2_relevant*/, - const DoFHandler<1,3> &/*dof2*/, - OutVector &/*u2*/) - {} - }; - - - - template - ExtrapolateImplementation:: - ExtrapolateImplementation () - : round(0) - {} - - - - 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 (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 - 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 (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 != nullptr, 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 != nullptr, 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); - const int ierr = MPI_Isend (&(*buffer)[0], buffer->size(), - MPI_BYTE, - it->receiver, - round, - communicator, - &requests[idx]); - AssertThrowMPI(ierr); - } - - 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) - { - const int ierr = MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE); - AssertThrowMPI(ierr); - } - - // 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 FiniteElement &fe = dealii_cell->get_dof_handler().get_fe(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - 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; - - // reset the round count we will use in send_cells - round = 0; - - // Compute all the cells needed from other processes. - compute_needs (dof2, cells_we_need); - - // Send the cells needed to their 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; - } - } - - // increase the round counter, such that we are sure to only send - // and receive data from the correct call - ++round; - - send_cells (cells_to_send, received_cells); - - // store received cell_data - for (typename std::vector::const_iterator recv=received_cells.begin (); - recv != received_cells.end (); - ++recv) - { - cell_data_insert (*recv, available_cells); - } - - // increase the round counter, such that we are sure to only send - // and receive data from the correct call - ++round; - - // 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 != nullptr, - 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 reinit_distributed(const DH &dh, - VectorType &vector) - { - vector.reinit(dh.n_dofs()); - } - -#ifdef DEAL_II_WITH_PETSC - template - void reinit_distributed(const DoFHandler &dh, - PETScWrappers::MPI::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert (parallel_tria != nullptr, ExcNotImplemented()); - - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); - } -#endif //DEAL_II_WITH_PETSC - -#ifdef DEAL_II_WITH_TRILINOS - template - void reinit_distributed(const DoFHandler &dh, - TrilinosWrappers::MPI::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert (parallel_tria != nullptr, ExcNotImplemented()); - - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); - } -#endif //DEAL_II_WITH_TRILINOS - - template - void reinit_distributed(const DoFHandler &dh, - LinearAlgebra::distributed::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert (parallel_tria != nullptr, ExcNotImplemented()); - - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); - } - - - - template - void reinit_ghosted(const DH &/*dh*/, - VectorType &/*vector*/) - { - Assert(false, ExcNotImplemented()); - } - -#ifdef DEAL_II_WITH_PETSC - template - void reinit_ghosted(const DoFHandler &dh, - PETScWrappers::MPI::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert(parallel_tria != nullptr, ExcNotImplemented()); - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); - vector.reinit (locally_owned_dofs, locally_relevant_dofs, - parallel_tria->get_communicator()); - } -#endif //DEAL_II_WITH_PETSC - -#ifdef DEAL_II_WITH_TRILINOS - template - void reinit_ghosted(const DoFHandler &dh, - TrilinosWrappers::MPI::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert (parallel_tria != nullptr, ExcNotImplemented()); - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); - vector.reinit (locally_owned_dofs, locally_relevant_dofs, - parallel_tria->get_communicator()); - } -#endif //DEAL_II_WITH_TRILINOS - - template - void reinit_ghosted(const DoFHandler &dh, - LinearAlgebra::distributed::Vector &vector) - { - const parallel::distributed::Triangulation *parallel_tria = - dynamic_cast*>(&dh.get_tria()); - Assert (parallel_tria != nullptr, ExcNotImplemented()); - const IndexSet locally_owned_dofs = dh.locally_owned_dofs(); - IndexSet locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs); - vector.reinit (locally_owned_dofs, locally_relevant_dofs, - parallel_tria->get_communicator()); - } - - - - template - void extrapolate_serial(const InVector &u3, - const DoFHandler &dof2, - OutVector &u2) - { - const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; - Vector dof_values(dofs_per_cell); - - // then traverse grid bottom up - for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), - endc=dof2.end(level); - - for (; cell!=endc; ++cell) - if (!cell->active()) - { - // check whether this - // cell has active - // children - bool active_children=false; - for (unsigned int child_n=0; child_nn_children(); ++child_n) - if (cell->child(child_n)->active()) - { - active_children=true; - break; - } - - // if there are active - // children, this process - // has to work on this - // cell. get the data - // from the one vector - // and set it on the - // other - if (active_children) - { - cell->get_interpolated_dof_values(u3, dof_values); - cell->set_dof_values_by_interpolation(dof_values, u2); - } - } - } - } - } - } - - template - void extrapolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - 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() || cell->is_artificial(), ExcGridNotRefinedAtLeastOnce()); - } - - - OutVector u3; - internal::reinit_distributed(dof2, u3); - if (dynamic_cast*>(&dof2.get_tria()) != nullptr) - { - interpolate(dof1, u1, dof2, constraints, u3); - - OutVector u3_relevant; - internal::reinit_ghosted(dof2, u3_relevant); - u3_relevant = u3; - - internal::ExtrapolateImplementation implementation; - implementation.extrapolate_parallel (u3_relevant, dof2, u2); - } - else - { - interpolate(dof1, u1, dof2, constraints, u3); - - internal::extrapolate_serial (u3, dof2, u2); - } - - constraints.distribute(u2); - } -} // end of namespace FETools - - /*-------------- Explicit Instantiations -------------------------------*/ #include "fe_tools_extrapolate.inst" diff --git a/source/fe/fe_tools_interpolate.cc b/source/fe/fe_tools_interpolate.cc index 8ea418fa52..97a014bbcc 100644 --- a/source/fe/fe_tools_interpolate.cc +++ b/source/fe/fe_tools_interpolate.cc @@ -14,643 +14,11 @@ // --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include - -#include -#include - +#include DEAL_II_NAMESPACE_OPEN -namespace FETools -{ - template class DoFHandlerType1, - template class DoFHandlerType2, - class InVector, class OutVector> - void - interpolate(const DoFHandlerType1 &dof1, - const InVector &u1, - const DoFHandlerType2 &dof2, - OutVector &u2) - { - ConstraintMatrix dummy; - dummy.close(); - interpolate(dof1, u1, dof2, dummy, u2); - } - - - - template class DoFHandlerType1, - template class DoFHandlerType2, - class InVector, class OutVector> - void - interpolate (const DoFHandlerType1 &dof1, - const InVector &u1, - const DoFHandlerType2 &dof2, - const ConstraintMatrix &constraints, - OutVector &u2) - { - 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())); - - - const IndexSet u2_elements = u2.locally_owned_elements(); -#ifdef DEBUG - const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); - const IndexSet &dof2_local_dofs = dof2.locally_owned_dofs(); - const IndexSet u1_elements = u1.locally_owned_elements(); - Assert(u1_elements == dof1_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); - Assert(u2_elements == dof2_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); -#endif - - // allocate vectors at maximal - // size. will be reinited in inner - // cell, but Vector makes sure that - // this does not lead to - // reallocation of memory - Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); - Vector u2_local(DoFTools::max_dofs_per_cell(dof2)); - - // have a map for interpolation - // matrices. shared_ptr make sure - // that memory is released again - std::map *, - std::map *, - std::shared_ptr > > > - interpolation_matrices; - - typename DoFHandlerType1::active_cell_iterator cell1 = dof1.begin_active(), - endc1 = dof1.end(); - typename DoFHandlerType2::active_cell_iterator cell2 = dof2.begin_active(), - endc2 = dof2.end(); - (void)endc2; - - std::vector dofs; - dofs.reserve (DoFTools::max_dofs_per_cell (dof2)); - - u2 = 0; - OutVector touch_count(u2); - touch_count = 0; - - // for distributed triangulations, - // we can only interpolate u1 on - // a cell, which this processor owns, - // so we have to know the subdomain_id - const types::subdomain_id subdomain_id = - dof1.get_triangulation().locally_owned_subdomain(); - - for (; cell1!=endc1; ++cell1, ++cell2) - if ((cell1->subdomain_id() == subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id)) - { - Assert(cell1->get_fe().n_components() == cell2->get_fe().n_components(), - ExcDimensionMismatch (cell1->get_fe().n_components(), - cell2->get_fe().n_components())); - - // for continuous elements on - // grids with hanging nodes we - // need hanging node - // constraints. Consequently, - // if there are no constraints - // then hanging nodes are not - // allowed. - const bool hanging_nodes_not_allowed - = ((cell2->get_fe().dofs_per_vertex != 0) && - (constraints.n_constraints() == 0)); - - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell1->at_boundary(face) || - cell1->neighbor(face)->level() == cell1->level(), - ExcHangingNodesNotAllowed(0)); - - - const unsigned int dofs_per_cell1 = cell1->get_fe().dofs_per_cell; - const unsigned int dofs_per_cell2 = cell2->get_fe().dofs_per_cell; - u1_local.reinit (dofs_per_cell1); - u2_local.reinit (dofs_per_cell2); - - // check if interpolation - // matrix for this particular - // pair of elements is already - // there - if (interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()].get() == nullptr) - { - std::shared_ptr > - interpolation_matrix (new FullMatrix (dofs_per_cell2, - dofs_per_cell1)); - interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] - = interpolation_matrix; - - get_interpolation_matrix(cell1->get_fe(), - cell2->get_fe(), - *interpolation_matrix); - } - - cell1->get_dof_values(u1, u1_local); - interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] - ->vmult(u2_local, u1_local); - - dofs.resize (dofs_per_cell2); - cell2->get_dof_indices(dofs); - - for (unsigned int i=0; i(touch_count(i)) != typename OutVector::value_type(0), - ExcInternalError()); - u2(i) /= touch_count(i); - } - - // finish the work on parallel vectors - u2.compress(VectorOperation::insert); - // Apply hanging node constraints. - constraints.distribute(u2); - } - - - - template class DoFHandlerType, - class InVector, class OutVector, int spacedim> - void - back_interpolate(const DoFHandlerType &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated) - { - Assert(dof1.get_fe().n_components() == fe2.n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); - Assert(u1.size() == dof1.n_dofs(), - ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_interpolated.size() == dof1.n_dofs(), - ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); - -#ifdef DEBUG - const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); - const IndexSet u1_elements = u1.locally_owned_elements(); - const IndexSet u1_interpolated_elements = u1_interpolated.locally_owned_elements(); - Assert(u1_elements == dof1_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); - Assert(u1_interpolated_elements == dof1_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); -#endif - - Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); - Vector u1_int_local(DoFTools::max_dofs_per_cell(dof1)); - - const types::subdomain_id subdomain_id = - dof1.get_triangulation().locally_owned_subdomain(); - - typename DoFHandlerType::active_cell_iterator cell = dof1.begin_active(), - endc = dof1.end(); - - // map from possible fe objects in - // dof1 to the back_interpolation - // matrices - std::map *, - std::shared_ptr > > interpolation_matrices; - - for (; cell!=endc; ++cell) - if ((cell->subdomain_id() == subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id)) - { - // For continuous elements on - // grids with hanging nodes we - // need hanging node - // constraints. Consequently, - // when the elements are - // continuous no hanging node - // constraints are allowed. - const bool hanging_nodes_not_allowed= - (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); - - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell->at_boundary(face) || - cell->neighbor(face)->level() == cell->level(), - ExcHangingNodesNotAllowed(0)); - - const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell; - - // make sure back_interpolation - // matrix is available - if (interpolation_matrices[&cell->get_fe()] == nullptr) - { - interpolation_matrices[&cell->get_fe()] = - std::shared_ptr > - (new FullMatrix(dofs_per_cell1, dofs_per_cell1)); - get_back_interpolation_matrix(cell->get_fe(), fe2, - *interpolation_matrices[&cell->get_fe()]); - } - - u1_local.reinit (dofs_per_cell1); - u1_int_local.reinit (dofs_per_cell1); - - cell->get_dof_values(u1, u1_local); - interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local); - cell->set_dof_values(u1_int_local, u1_interpolated); - }; - - // if we work on a parallel vector, we have to finish the work - u1_interpolated.compress(VectorOperation::insert); - } - - - - namespace internal - { - namespace - { - template - void back_interpolate (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - InVector &u1_interpolated) - { - Vector u2(dof2.n_dofs()); - interpolate(dof1, u1, dof2, constraints2, u2); - interpolate(dof2, u2, dof1, constraints1, u1_interpolated); - } - - // special version for PETSc -#ifdef DEAL_II_WITH_PETSC - template - void back_interpolate (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const PETScWrappers::MPI::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - PETScWrappers::MPI::Vector &u1_interpolated) - { - // if u1 is a parallel distributed PETSc vector, we create a - // vector u2 with based on the sets of locally owned and relevant - // dofs of dof2 - IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); - IndexSet dof2_locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dof2, - dof2_locally_relevant_dofs); - - PETScWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs, - u1.get_mpi_communicator()); - interpolate(dof1, u1, dof2, constraints2, u2_out); - PETScWrappers::MPI::Vector u2 (dof2_locally_owned_dofs, - dof2_locally_relevant_dofs, - u1.get_mpi_communicator()); - u2 = u2_out; - interpolate(dof2, u2, dof1, constraints1, u1_interpolated); - } -#endif - - // special version for Trilinos -#ifdef DEAL_II_WITH_TRILINOS - template - void back_interpolate (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const TrilinosWrappers::MPI::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - TrilinosWrappers::MPI::Vector &u1_interpolated) - { - // if u1 is a parallel distributed Trilinos vector, we create a - // vector u2 with based on the sets of locally owned and relevant - // dofs of dof2 - IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs(); - IndexSet dof2_locally_relevant_dofs; - DoFTools::extract_locally_relevant_dofs (dof2, - dof2_locally_relevant_dofs); - - TrilinosWrappers::MPI::Vector u2_out (dof2_locally_owned_dofs, - u1.get_mpi_communicator()); - interpolate(dof1, u1, dof2, constraints2, u2_out); - TrilinosWrappers::MPI::Vector u2 (dof2_locally_owned_dofs, - dof2_locally_relevant_dofs, - u1.get_mpi_communicator()); - u2 = u2_out; - interpolate(dof2, u2, dof1, constraints1, u1_interpolated); - } -#endif - - // special version for LinearAlgebra::distributed::Vector - template - void back_interpolate (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const LinearAlgebra::distributed::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - LinearAlgebra::distributed::Vector &u1_interpolated) - { - 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 - u2 (dof2_locally_owned_dofs, - dof2_locally_relevant_dofs, - u1.get_mpi_communicator()); - - interpolate(dof1, u1, dof2, constraints2, u2); - u2.update_ghost_values (); - interpolate(dof2, u2, dof1, constraints1, u1_interpolated); - } - } - } - - - - template - void back_interpolate(const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - OutVector &u1_interpolated) - { - // For discontinuous elements without constraints take the simpler version - // of the back_interpolate function. - if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 - && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) - back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated); - else - { - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_interpolated.size()==dof1.n_dofs(), - ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); - - // For continuous elements first interpolate to dof2, taking into - // account constraints2, and then interpolate back to dof1 taking into - // account constraints1 - internal::back_interpolate(dof1, constraints1, u1, dof2, constraints2, - u1_interpolated); - } - } - - - - template - void interpolation_difference (const DoFHandler &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_difference) - { - Assert(dof1.get_fe().n_components() == fe2.n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_difference.size()==dof1.n_dofs(), - ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs())); - -#ifdef DEBUG - const IndexSet &dof1_local_dofs = dof1.locally_owned_dofs(); - const IndexSet u1_elements = u1.locally_owned_elements(); - const IndexSet u1_difference_elements = u1_difference.locally_owned_elements(); - Assert(u1_elements == dof1_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); - Assert(u1_difference_elements == dof1_local_dofs, - ExcMessage("The provided vector and DoF handler should have the same" - " index sets.")); -#endif - - // For continuous elements on grids - // with hanging nodes we need - // hanging node - // constraints. Consequently, when - // the elements are continuous no - // hanging node constraints are - // allowed. - const bool hanging_nodes_not_allowed= - (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); - - const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell; - - Vector u1_local(dofs_per_cell); - Vector u1_diff_local(dofs_per_cell); - - const types::subdomain_id subdomain_id = - dof1.get_triangulation().locally_owned_subdomain(); - - FullMatrix difference_matrix(dofs_per_cell, dofs_per_cell); - get_interpolation_difference_matrix(dof1.get_fe(), fe2, - difference_matrix); - - typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), - endc = dof1.end(); - - for (; cell!=endc; ++cell) - if ((cell->subdomain_id() == subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id)) - { - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell->at_boundary(face) || - cell->neighbor(face)->level() == cell->level(), - ExcHangingNodesNotAllowed(0)); - - cell->get_dof_values(u1, u1_local); - difference_matrix.vmult(u1_diff_local, u1_local); - cell->set_dof_values(u1_diff_local, u1_difference); - } - - // if we work on a parallel PETSc vector - // we have to finish the work and - // update ghost values - u1_difference.compress(VectorOperation::insert); - } - - - - namespace internal - { - namespace - { - template - void interpolation_difference (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - OutVector &u1_difference) - { - back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); - u1_difference.sadd(-1., 1., u1); - } - - // special version for Trilinos -#ifdef DEAL_II_WITH_TRILINOS - template - void interpolation_difference (const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const TrilinosWrappers::MPI::Vector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - TrilinosWrappers::MPI::Vector &u1_difference) - { - back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); - - // Trilinos vectors with and without ghost entries are very different - // and we cannot use the sadd function directly, so we have to create - // a completely distributed vector first and copy the local entries - // from the vector with ghost entries - TrilinosWrappers::MPI::Vector u1_completely_distributed (u1_difference.vector_partitioner ()); - - u1_completely_distributed = u1; - - u1_difference.sadd(-1, u1_completely_distributed); - } -#endif - } - } - - - - template - void interpolation_difference(const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - OutVector &u1_difference) - { - // For discontinuous elements - // without constraints take the - // cheaper version of the - // interpolation_difference function. - if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 - && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) - interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference); - else - { - internal::interpolation_difference(dof1, constraints1, u1, dof2, constraints2, u1_difference); - } - } - - - - template - void project_dg(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - OutVector &u2) - { - Assert(&dof1.get_triangulation()==&dof2.get_triangulation(), ExcTriangulationMismatch()); - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); - - typename DoFHandler::active_cell_iterator cell1 = dof1.begin_active(); - typename DoFHandler::active_cell_iterator cell2 = dof2.begin_active(); - typename DoFHandler::active_cell_iterator end = dof2.end(); - - const unsigned int n1 = dof1.get_fe().dofs_per_cell; - const unsigned int n2 = dof2.get_fe().dofs_per_cell; - - Vector u1_local(n1); - Vector u2_local(n2); - std::vector dofs(n2); - - FullMatrix matrix(n2,n1); - get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix); - - u2 = 0; - while (cell2 != end) - { - cell1->get_dof_values(u1, u1_local); - matrix.vmult(u2_local, u1_local); - cell2->get_dof_indices(dofs); - for (unsigned int i=0; i +#include "fe_tools_extrapolate_common.h" + +// check FETools::extrapolate on distributed triangulations +// for TrilinosWrappers::MPI::BlockVector + +template +void +check (const FiniteElement &fe1, + const FiniteElement &fe2, + const std::string &name) +{ + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) + deallog << "Checking " << name + << " in " << dim << "d:" + << std::endl; + + // call main function in .cc files + check_this (fe1, fe2); +} + +#define CHECK(EL1,deg1,EL2,deg2,dim)\ + { FE_ ## EL1 fe1(deg1); \ + FE_ ## EL2 fe2(deg2); \ + check(fe1, fe2, #EL1 #deg1 " against " #EL2 #deg2); \ + } + +#define CHECK_ALL(EL1,deg1,EL2,deg2)\ + { CHECK(EL1,deg1,EL2,deg2,2); \ + CHECK(EL1,deg1,EL2,deg2,3); \ + } + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + MPILogInitAll log; + deallog.threshold_double(1.e-10); + + CHECK_ALL(Q,1,Q,1); + CHECK_ALL(Q,1,Q,2); + CHECK_ALL(Q,2,Q,1); + CHECK_ALL(Q,2,Q,2); + + CHECK_ALL(DGQ,0,DGQ,0); + CHECK_ALL(DGQ,0,DGQ,1); + CHECK_ALL(DGQ,1,DGQ,0); + CHECK_ALL(DGQ,1,DGQ,1); +} + diff --git a/tests/mpi/fe_tools_extrapolate_04.with_trilinos=true.mpirun=1.output b/tests/mpi/fe_tools_extrapolate_04.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..39f85f69f6 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_04.with_trilinos=true.mpirun=1.output @@ -0,0 +1,97 @@ + +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_04.with_trilinos=true.mpirun=3.output b/tests/mpi/fe_tools_extrapolate_04.with_trilinos=true.mpirun=3.output new file mode 100644 index 0000000000..8f5cfcce45 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_04.with_trilinos=true.mpirun=3.output @@ -0,0 +1,101 @@ + +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_05.cc b/tests/mpi/fe_tools_extrapolate_05.cc new file mode 100644 index 0000000000..7972a66e89 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_05.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 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 "fe_tools_extrapolate_common.h" + +// check FETools::extrapolate on distributed triangulations +// for PETScWrappers::MPI::BlockVector + +template +void +check (const FiniteElement &fe1, + const FiniteElement &fe2, + const std::string &name) +{ + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) + deallog << "Checking " << name + << " in " << dim << "d:" + << std::endl; + + // call main function in .cc files + check_this (fe1, fe2); +} + +#define CHECK(EL1,deg1,EL2,deg2,dim)\ + { FE_ ## EL1 fe1(deg1); \ + FE_ ## EL2 fe2(deg2); \ + check(fe1, fe2, #EL1 #deg1 " against " #EL2 #deg2); \ + } + +#define CHECK_ALL(EL1,deg1,EL2,deg2)\ + { CHECK(EL1,deg1,EL2,deg2,2); \ + CHECK(EL1,deg1,EL2,deg2,3); \ + } + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + MPILogInitAll log; + deallog.threshold_double(1.e-10); + + CHECK_ALL(Q,1,Q,1); + CHECK_ALL(Q,1,Q,2); + CHECK_ALL(Q,2,Q,1); + CHECK_ALL(Q,2,Q,2); + + CHECK_ALL(DGQ,0,DGQ,0); + CHECK_ALL(DGQ,0,DGQ,1); + CHECK_ALL(DGQ,1,DGQ,0); + CHECK_ALL(DGQ,1,DGQ,1); +} + diff --git a/tests/mpi/fe_tools_extrapolate_05.with_petsc=true.with_petsc_with_complex=false.mpirun=1.output b/tests/mpi/fe_tools_extrapolate_05.with_petsc=true.with_petsc_with_complex=false.mpirun=1.output new file mode 100644 index 0000000000..39f85f69f6 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_05.with_petsc=true.with_petsc_with_complex=false.mpirun=1.output @@ -0,0 +1,97 @@ + +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_05.with_petsc=true.with_petsc_with_complex=false.mpirun=3.output b/tests/mpi/fe_tools_extrapolate_05.with_petsc=true.with_petsc_with_complex=false.mpirun=3.output new file mode 100644 index 0000000000..8f5cfcce45 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_05.with_petsc=true.with_petsc_with_complex=false.mpirun=3.output @@ -0,0 +1,101 @@ + +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_06.cc b/tests/mpi/fe_tools_extrapolate_06.cc new file mode 100644 index 0000000000..2c3befabad --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_06.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2003 - 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 "fe_tools_extrapolate_common.h" + +// check FETools::extrapolate on distributed triangulations +// for LinearAlgebra::distributed::BlockVector + +template +void +check (const FiniteElement &fe1, + const FiniteElement &fe2, + const std::string &name) +{ + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0) + deallog << "Checking " << name + << " in " << dim << "d:" + << std::endl; + + // call main function in .cc files + check_this_dealii > (fe1, fe2); +} + +#define CHECK(EL1,deg1,EL2,deg2,dim)\ + { FE_ ## EL1 fe1(deg1); \ + FE_ ## EL2 fe2(deg2); \ + check(fe1, fe2, #EL1 #deg1 " against " #EL2 #deg2); \ + } + +#define CHECK_ALL(EL1,deg1,EL2,deg2)\ + { CHECK(EL1,deg1,EL2,deg2,2); \ + CHECK(EL1,deg1,EL2,deg2,3); \ + } + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + MPILogInitAll log; + deallog.threshold_double(1.e-10); + + CHECK_ALL(Q,1,Q,1); + CHECK_ALL(Q,1,Q,2); + CHECK_ALL(Q,2,Q,1); + CHECK_ALL(Q,2,Q,2); + + CHECK_ALL(DGQ,0,DGQ,0); + CHECK_ALL(DGQ,0,DGQ,1); + CHECK_ALL(DGQ,1,DGQ,0); + CHECK_ALL(DGQ,1,DGQ,1); +} + diff --git a/tests/mpi/fe_tools_extrapolate_06.mpirun=1.output b/tests/mpi/fe_tools_extrapolate_06.mpirun=1.output new file mode 100644 index 0000000000..39f85f69f6 --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_06.mpirun=1.output @@ -0,0 +1,97 @@ + +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_06.mpirun=10.output b/tests/mpi/fe_tools_extrapolate_06.mpirun=10.output new file mode 100644 index 0000000000..f406bbaa3d --- /dev/null +++ b/tests/mpi/fe_tools_extrapolate_06.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 index b0eebff448..8b649dc341 100755 --- a/tests/mpi/fe_tools_extrapolate_common.h +++ b/tests/mpi/fe_tools_extrapolate_common.h @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -117,6 +118,43 @@ output_vector (const VectorType &v, const std::string &output_name, const DoFHan +template +typename std::enable_if::value, VectorType>::type +build_ghosted (const IndexSet &owned_indices, + const IndexSet &ghosted_indices) +{ + return VectorType(owned_indices, ghosted_indices, MPI_COMM_WORLD); +} + +template +typename std::enable_if::value, VectorType>::type +build_ghosted (const IndexSet &owned_indices, + const IndexSet &ghosted_indices) +{ + std::vector owned_indices_vector(1, owned_indices); + std::vector ghosted_indices_vector(1, ghosted_indices); + return VectorType(owned_indices_vector, ghosted_indices_vector, MPI_COMM_WORLD); +} + + + +template +typename std::enable_if::value, VectorType>::type +build_distributed (const IndexSet &owned_indices) +{ + return VectorType(owned_indices, MPI_COMM_WORLD); +} + +template +typename std::enable_if::value, VectorType>::type +build_distributed (const IndexSet &owned_indices) +{ + std::vector owned_indices_vector(1, owned_indices); + return VectorType(owned_indices_vector, MPI_COMM_WORLD); +} + + + template void check_this (const FiniteElement &fe1, @@ -157,11 +195,11 @@ check_this (const FiniteElement &fe1, IndexSet locally_relevant_dofs2; DoFTools::extract_locally_relevant_dofs (*dof2, locally_relevant_dofs2); - VectorType in_ghosted (locally_owned_dofs1, locally_relevant_dofs1, MPI_COMM_WORLD); - VectorType in_distributed (locally_owned_dofs1, MPI_COMM_WORLD); - VectorType out_distributed (locally_owned_dofs2, MPI_COMM_WORLD); - VectorType out_ghosted (locally_owned_dofs2, locally_relevant_dofs2, MPI_COMM_WORLD); - VectorType out_reference (locally_owned_dofs2, MPI_COMM_WORLD); + VectorType in_ghosted = build_ghosted (locally_owned_dofs1, locally_relevant_dofs1); + VectorType in_distributed = build_distributed (locally_owned_dofs1); + VectorType out_distributed = build_distributed (locally_owned_dofs2); + VectorType out_ghosted = build_ghosted (locally_owned_dofs2, locally_relevant_dofs2); + VectorType out_reference = build_distributed (locally_owned_dofs2); // Choose some reference function of sufficiently high polynomial degree. TestFunction function (std::max(fe1.degree,fe2.degree)); @@ -281,9 +319,9 @@ check_this_dealii (const FiniteElement &fe1, IndexSet locally_relevant_dofs2; DoFTools::extract_locally_relevant_dofs (*dof2, locally_relevant_dofs2); - VectorType in_ghosted (locally_owned_dofs1, locally_relevant_dofs1, MPI_COMM_WORLD); - VectorType out_ghosted (locally_owned_dofs2, locally_relevant_dofs2, MPI_COMM_WORLD); - VectorType out_reference (locally_owned_dofs2, locally_relevant_dofs2, MPI_COMM_WORLD); + VectorType in_ghosted = build_ghosted (locally_owned_dofs1, locally_relevant_dofs1); + VectorType out_ghosted = build_ghosted (locally_owned_dofs2, locally_relevant_dofs2); + VectorType out_reference = build_ghosted (locally_owned_dofs2, locally_relevant_dofs2); // Choose some reference function of sufficiently high polynomial degree. TestFunction function (std::max(fe1.degree,fe2.degree));