--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/p4est_wrappers.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_tools_interpolate.templates.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/la_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+
+#include <queue>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FETools
+{
+
+ namespace internal
+ {
+#ifndef DEAL_II_WITH_P4EST
+ // Dummy implementation in case p4est is not available.
+ template <int dim,int spacedim,class OutVector>
+ class ExtrapolateImplementation
+ {
+ public:
+
+ ExtrapolateImplementation ()
+ {
+ Assert(false, ExcNotImplemented());
+ };
+
+ template <class InVector>
+ void extrapolate_parallel (const InVector &/*u2_relevant*/,
+ const DoFHandler<dim,spacedim> &/*dof2*/,
+ OutVector &/*u2*/)
+ {
+ Assert(false, ExcNotImplemented())
+ }
+ };
+#else
+ // Implementation of the @p extrapolate function
+ // on parallel distributed grids.
+ template <int dim,int spacedim,class OutVector>
+ class ExtrapolateImplementation
+ {
+ public:
+
+ ExtrapolateImplementation ();
+
+ template <class InVector>
+ void extrapolate_parallel (const InVector &u2_relevant,
+ const DoFHandler<dim,spacedim> &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<dim>::forest forest;
+ const typename dealii::internal::p4est::types<dim>::tree tree;
+ const typename dealii::internal::p4est::types<dim>::locidx tree_index;
+ const typename DoFHandler<dim,spacedim>::cell_iterator dealii_cell;
+ const typename dealii::internal::p4est::types<dim>::quadrant p4est_cell;
+
+ WorkPackage(const typename dealii::internal::p4est::types<dim>::forest &forest_,
+ const typename dealii::internal::p4est::types<dim>::tree &tree_,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index_,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell_,
+ const typename dealii::internal::p4est::types<dim>::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<value_type> dof_values;
+
+ unsigned int tree_index;
+
+ typename dealii::internal::p4est::types<dim>::quadrant quadrant;
+
+ int receiver;
+
+ bool operator < (const CellData &rhs) const
+ {
+ if (dealii::internal::p4est::functions<dim>::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<dim>::quadrant)); // quadrant
+ }
+
+ void pack_data (std::vector<char> &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<dim>::quadrant));
+ ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
+
+ Assert (ptr == &buffer[0]+buffer.size(),
+ ExcInternalError());
+ }
+
+ void unpack_data (const std::vector<char> &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<dim>::quadrant));
+ ptr += sizeof(typename dealii::internal::p4est::types<dim>::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 <class InVector>
+ void compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
+ const InVector &u);
+
+ // traverse recursively over
+ // the cells of this tree and
+ // interpolate over patches which
+ // are part of this process
+ template <class InVector>
+ void interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::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 <class InVector>
+ void get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const InVector &u,
+ Vector<value_type> &interpolated_values,
+ std::vector<CellData> &new_needs);
+
+ // set dof values for this
+ // cell by interpolation
+ void set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const Vector<value_type> &interpolated_values,
+ OutVector &u);
+
+ // compute all cell_data
+ // needed from other processes
+ // to interpolate the part of
+ // this process
+ void compute_needs (const DoFHandler<dim,spacedim> &dof2,
+ std::vector<CellData> &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<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &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<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &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 <class InVector>
+ void compute_cells (const DoFHandler<dim,spacedim> &dof2,
+ const InVector &u,
+ std::vector<CellData> &cells_to_compute,
+ std::vector<CellData> &computed_cells,
+ std::vector<CellData> &new_needs);
+
+ // traverse over the tree
+ // and compute cells
+ template <class InVector>
+ void compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const InVector &u,
+ std::vector<CellData> &cells_to_compute,
+ std::vector<CellData> &computed_cells,
+ std::vector<CellData> &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<CellData> &cells_to_send,
+ std::vector<CellData> &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<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &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<CellData> &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<CellData> &cells_list);
+
+ MPI_Comm communicator;
+
+ // a vector of all cells this process has
+ // computed or received data
+ std::vector<CellData> available_cells;
+
+ // stores the indices of dofs on more refined ghosted cells along
+ // with the maximum level
+ std::map<types::global_dof_index, int> dofs_on_refined_neighbors;
+
+ // counts the send/receive round we are in
+ unsigned int round;
+ };
+
+ template <class OutVector>
+ class ExtrapolateImplementation<1,1,OutVector>
+ {
+ public:
+
+ ExtrapolateImplementation ()
+ {
+ AssertThrow(false, ExcNotImplemented())
+ }
+
+ template <class InVector>
+ void extrapolate_parallel (const InVector &/*u2_relevant*/,
+ const DoFHandler<1,1> &/*dof2*/,
+ OutVector &/*u2*/)
+ {}
+ };
+
+ template <class OutVector>
+ class ExtrapolateImplementation<1,2,OutVector>
+ {
+ public:
+
+ ExtrapolateImplementation ()
+ {
+ AssertThrow(false, ExcNotImplemented())
+ }
+
+ template <class InVector>
+ void extrapolate_parallel (const InVector &/*u2_relevant*/,
+ const DoFHandler<1,2> &/*dof2*/,
+ OutVector &/*u2*/)
+ {}
+ };
+
+ template <class OutVector>
+ class ExtrapolateImplementation<1,3,OutVector>
+ {
+ public:
+
+ ExtrapolateImplementation ()
+ {
+ AssertThrow(false, ExcNotImplemented())
+ }
+
+ template <class InVector>
+ void extrapolate_parallel (const InVector &/*u2_relevant*/,
+ const DoFHandler<1,3> &/*dof2*/,
+ OutVector &/*u2*/)
+ {}
+ };
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ ExtrapolateImplementation ()
+ : round(0)
+ {}
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ CellData::CellData ()
+ : tree_index (0),
+ receiver (0)
+ {}
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ CellData::CellData (const unsigned int dofs_per_cell)
+ : tree_index (0),
+ receiver (0)
+ {
+ dof_values.reinit (dofs_per_cell);
+ }
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const InVector &u1,
+ OutVector &u2)
+ {
+ // check if this cell exists in the local p4est
+ int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
+ &p4est_cell,
+ dealii::internal::p4est::functions<dim>::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<dim>::
+ quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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_n<dealii_cell->n_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<dim,spacedim> &fe = dealii_cell->get_dof_handler().get_fe();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ Vector<typename OutVector::value_type> interpolated_values(dofs_per_cell);
+
+ std::vector<CellData> 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 <int dim,int spacedim, class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const InVector &u,
+ Vector<value_type> &interpolated_values,
+ std::vector<CellData> &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<dim,spacedim> &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<value_type> tmp1(dofs_per_cell);
+ Vector<value_type> tmp2(dofs_per_cell);
+
+ interpolated_values = 0;
+ std::vector<bool> restriction_is_additive (dofs_per_cell);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ restriction_is_additive[i] = fe.restriction_is_additive(i);
+
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ bool found_child = true;
+ for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+ {
+ if (dealii::internal::p4est::functions<dim>::
+ quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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<dofs_per_cell; ++i)
+ if (tmp2(i) != value_type())
+ {
+ if (restriction_is_additive[i])
+ interpolated_values(i) += tmp2(i);
+ else
+ interpolated_values(i) = tmp2(i);
+ }
+ }
+ }
+
+ if (found_child == false)
+ interpolated_values = 0;
+ }
+ }
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const Vector<value_type> &local_values,
+ OutVector &u)
+ {
+ const FiniteElement<dim,spacedim> &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<types::global_dof_index> indices(dofs_per_cell);
+ dealii_cell->get_dof_indices(indices);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ // don't set this dof if there is a more refined ghosted neighbor setting this dof.
+ const bool on_refined_neighbor
+ = (dofs_on_refined_neighbors.find(indices[j])!=dofs_on_refined_neighbors.end());
+ if (!(on_refined_neighbor && dofs_on_refined_neighbors[indices[j]]>dealii_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<value_type> tmp(dofs_per_cell);
+
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::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 <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ compute_needs (const DoFHandler<dim,spacedim> &dof2,
+ std::vector<CellData> &new_needs)
+ {
+ const parallel::distributed::Triangulation< dim, spacedim > *tr
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&dof2.get_tria()));
+
+ Assert (tr != nullptr, ExcInternalError());
+
+ typename DoFHandler<dim,spacedim>::cell_iterator
+ cell=dof2.begin(0),
+ endc=dof2.end(0);
+ for (; cell!=endc; ++cell)
+ {
+ if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
+ tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
+ == false)
+ continue;
+
+ typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+ const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
+ typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+ dealii::internal::p4est::init_coarse_quadrant<dim>(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<sc_array_t *>(&tree->quadrants),
+ &p4est_coarse_cell,
+ dealii::internal::p4est::functions<dim>::quadrant_compare);
+
+ AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ());
+ }
+
+ traverse_tree_recursively (*tr->parallel_forest,
+ *tree,
+ tree_index,
+ cell,
+ p4est_coarse_cell,
+ new_needs);
+ }
+ }
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &new_needs)
+ {
+ // check if this cell exists in the local p4est
+ int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
+ &p4est_cell,
+ dealii::internal::p4est::functions<dim>::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<dim>::
+ quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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_n<dealii_cell->n_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<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+ {
+ traverse_tree_recursively (forest,
+ tree,
+ tree_index,
+ dealii_cell->child(c),
+ p4est_child[c],
+ new_needs);
+ }
+ }
+ }
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &new_needs)
+ {
+ if (dealii_cell->has_children ())
+ {
+ Assert (dealii_cell->has_children (), ExcInternalError ());
+
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+ {
+ // check if this child
+ // is locally available
+ // in the p4est
+ if (dealii::internal::p4est::functions<dim>::
+ quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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 <int dim,int spacedim,class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ compute_cells (const DoFHandler<dim,spacedim> &dof2,
+ const InVector &u,
+ std::vector<CellData> &cells_to_compute,
+ std::vector<CellData> &computed_cells,
+ std::vector<CellData> &new_needs)
+ {
+ const parallel::distributed::Triangulation< dim, spacedim > *tr
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&dof2.get_tria()));
+
+ Assert (tr != nullptr, ExcInternalError());
+
+ // collect in a set all trees this
+ // process has to compute cells on
+ std::set<unsigned int> trees;
+ for (typename std::vector<CellData>::const_iterator it=cells_to_compute.begin(); it!=cells_to_compute.end(); ++it)
+ trees.insert (it->tree_index);
+
+ typename DoFHandler<dim,spacedim>::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<dim> (tr->parallel_forest,
+ tree_index)
+ == false))
+ continue;
+
+ typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+ typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+ dealii::internal::p4est::init_coarse_quadrant<dim>(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 <int dim,int spacedim,class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::tree &tree,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ const InVector &u,
+ std::vector<CellData> &cells_to_compute,
+ std::vector<CellData> &computed_cells,
+ std::vector<CellData> &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<sc_array_t *>(&tree.quadrants),
+ &p4est_cell,
+ dealii::internal::p4est::functions<dim>::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<dim>::
+ quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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<CellData> 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<tmp.size (); ++i)
+ cell_data_insert (tmp[i], new_needs);
+ }
+ }
+
+ // search recursively over this tree
+ if (p4est_has_children)
+ {
+ typename dealii::internal::p4est::types<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::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 <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ send_cells (const std::vector<CellData> &cells_to_send,
+ std::vector<CellData> &received_cells) const
+ {
+ std::vector<std::vector<char> > sendbuffers (cells_to_send.size());
+ std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
+ std::vector<MPI_Request> requests (cells_to_send.size());
+ std::vector<unsigned int> destinations;
+
+ // send data
+ unsigned int idx=0;
+ for (typename std::vector<CellData>::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<unsigned int> senders
+ = Utilities::MPI::compute_point_to_point_communication_pattern(communicator, destinations);
+
+ // receive data
+ std::vector<char> receive;
+ CellData cell_data;
+ for (unsigned int index=0; index<senders.size (); ++index)
+ {
+ MPI_Status status;
+ int len;
+ int ierr = MPI_Probe(MPI_ANY_SOURCE, round, communicator, &status);
+ AssertThrowMPI(ierr);
+ ierr = MPI_Get_count(&status, MPI_BYTE, &len);
+ AssertThrowMPI(ierr);
+ receive.resize (len);
+
+ char *buf = &receive[0];
+ ierr = MPI_Recv (buf, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, communicator, &status);
+ AssertThrowMPI(ierr);
+
+ cell_data.unpack_data (receive);
+
+ // this process has to send this
+ // cell back to the sender
+ // the receiver is the old sender
+ cell_data.receiver = status.MPI_SOURCE;
+
+ received_cells.push_back (cell_data);
+ }
+
+ if (requests.size () > 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 <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ add_new_need (const typename dealii::internal::p4est::types<dim>::forest &forest,
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
+ const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
+ std::vector<CellData> &new_needs)
+ {
+ const FiniteElement<dim,spacedim> &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<dim>::
+ comm_find_owner (const_cast<typename dealii::internal::p4est::types<dim>::forest *> (&forest),
+ tree_index,
+ &p4est_cell,
+ dealii_cell->level_subdomain_id ());
+
+ cell_data_insert (cell_data, new_needs);
+ }
+
+
+
+ template <int dim,int spacedim,class OutVector>
+ int
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ cell_data_search (const CellData &cell_data,
+ const std::vector<CellData> &cells_list)
+ {
+ typename std::vector<CellData>::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 <int dim,int spacedim,class OutVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ cell_data_insert (const CellData &cell_data,
+ std::vector<CellData> &cells_list)
+ {
+ typename std::vector<CellData>::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 <int dim,int spacedim,class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
+ const InVector &u)
+ {
+ std::vector<CellData> 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<received_needs.size(); ++i)
+ cell_data_insert (received_needs[i], cells_to_compute);
+
+ compute_cells (dof2, u, cells_to_compute, computed_cells, new_needs);
+
+ // if there are no cells to compute and no new needs, stop
+ ready = Utilities::MPI::sum (new_needs.size () + cells_to_compute.size (), communicator);
+
+ for (typename std::vector<CellData>::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<CellData>::iterator recv=received_needs.begin();
+ while (recv != received_needs.end())
+ {
+ if (dealii::internal::p4est::quadrant_is_equal<dim>(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<CellData>::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 <int dim,int spacedim,class OutVector>
+ template <class InVector>
+ void
+ ExtrapolateImplementation<dim,spacedim,OutVector>::
+ extrapolate_parallel (const InVector &u2_relevant,
+ const DoFHandler<dim,spacedim> &dof2,
+ OutVector &u2)
+ {
+ const parallel::distributed::Triangulation< dim, spacedim > *tr
+ = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
+ (&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<dim,spacedim> &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<types::global_dof_index> indices (dofs_per_cell);
+ typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (!cell->at_boundary(face))
+ {
+ const typename DoFHandler<dim,spacedim>::cell_iterator neighbor = cell->neighbor(face);
+ if (neighbor->level() != cell->level())
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ {
+ const types::global_dof_index index = indices[fe.face_to_cell_index(i,face)];;
+ const bool index_stored
+ = (dofs_on_refined_neighbors.find(index)!=dofs_on_refined_neighbors.end());
+ const unsigned int level = index_stored ?
+ std::max(cell->level(), 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<WorkPackage> queue;
+ {
+ typename DoFHandler<dim,spacedim>::cell_iterator
+ cell=dof2.begin(0),
+ endc=dof2.end(0);
+ for (; cell!=endc; ++cell)
+ {
+ if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
+ tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
+ == false)
+ continue;
+
+ typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
+ const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
+ typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
+
+ dealii::internal::p4est::init_coarse_quadrant<dim>(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<dim>::forest &forest = data.forest;
+ const typename dealii::internal::p4est::types<dim>::tree &tree = data.tree;
+ const typename dealii::internal::p4est::types<dim>::locidx &tree_index= data.tree_index;
+ const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell = data.dealii_cell;
+ const typename dealii::internal::p4est::types<dim>::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<dim>::quadrant
+ p4est_child[GeometryInfo<dim>::max_children_per_cell];
+
+ dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
+ p4est_child);
+
+ for (unsigned int c=0; c<GeometryInfo<dim>::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 <class VectorType, typename dummy = void>
+ struct BlockTypeHelper
+ {
+ using type = VectorType;
+ };
+
+ template <class VectorType>
+ struct BlockTypeHelper<VectorType, typename std::enable_if<IsBlockVector<VectorType>::value>::type>
+ {
+ using type = typename VectorType::BlockType;
+ };
+
+ template <class VectorType>
+ using BlockType = typename BlockTypeHelper<VectorType>::type;
+
+ template <class VectorType, class DH>
+ void reinit_distributed(const DH &dh,
+ VectorType &vector)
+ {
+ vector.reinit(dh.n_dofs());
+ }
+
+#ifdef DEAL_II_WITH_PETSC
+ template <int dim, int spacedim>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ PETScWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ TrilinosWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim, typename Number>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ LinearAlgebra::distributed::Vector<Number> &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <class VectorType, class DH>
+ void reinit_ghosted(const DH &/*dh*/,
+ VectorType &/*vector*/)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+#ifdef DEAL_II_WITH_PETSC
+ template <int dim, int spacedim>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ PETScWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ TrilinosWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim, typename Number>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ LinearAlgebra::distributed::Vector<Number> &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate_serial(const InVector &u3,
+ const DoFHandler<dim,spacedim> &dof2,
+ OutVector &u2)
+ {
+ const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell;
+ Vector<typename OutVector::value_type> dof_values(dofs_per_cell);
+
+ // then traverse grid bottom up
+ for (unsigned int level=0; level<dof2.get_triangulation().n_levels()-1; ++level)
+ {
+ typename DoFHandler<dim,spacedim>::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_n<cell->n_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 <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &dof2,
+ OutVector &u2)
+ {
+ ConstraintMatrix dummy;
+ dummy.close();
+ extrapolate(dof1, u1, dof2, dummy, u2);
+ }
+
+
+
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<dim,spacedim>::cell_iterator cell = dof2.begin(0),
+ endc = dof2.end(0);
+ for (; cell!=endc; ++cell)
+ Assert (cell->has_children() || cell->is_artificial(), ExcGridNotRefinedAtLeastOnce());
+ }
+
+
+ internal::BlockType<OutVector> u3;
+ internal::reinit_distributed(dof2, u3);
+ if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != nullptr)
+ {
+ interpolate(dof1, u1, dof2, constraints, u3);
+
+ internal::BlockType<OutVector> u3_relevant;
+ internal::reinit_ghosted(dof2, u3_relevant);
+ u3_relevant = u3;
+
+ internal::ExtrapolateImplementation<dim,spacedim,OutVector> 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/la_vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_block_vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/hp/dof_handler.h>
+
+
+#include <deal.II/base/index_set.h>
+
+#include <iostream>
+#include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace FETools
+{
+ template <int dim, int spacedim,
+ template <int, int> class DoFHandlerType1,
+ template <int, int> class DoFHandlerType2,
+ class InVector, class OutVector>
+ void
+ interpolate(const DoFHandlerType1<dim, spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandlerType2<dim, spacedim> &dof2,
+ OutVector &u2)
+ {
+ ConstraintMatrix dummy;
+ dummy.close();
+ interpolate(dof1, u1, dof2, dummy, u2);
+ }
+
+
+
+ template <int dim, int spacedim,
+ template <int, int> class DoFHandlerType1,
+ template <int, int> class DoFHandlerType2,
+ class InVector, class OutVector>
+ void
+ interpolate (const DoFHandlerType1<dim, spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandlerType2<dim, spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
+ Vector<typename OutVector::value_type> 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<const FiniteElement<dim,spacedim> *,
+ std::map<const FiniteElement<dim,spacedim> *,
+ std::shared_ptr<FullMatrix<double> > > >
+ interpolation_matrices;
+
+ typename DoFHandlerType1<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active(),
+ endc1 = dof1.end();
+ typename DoFHandlerType2<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active(),
+ endc2 = dof2.end();
+ (void)endc2;
+
+ std::vector<types::global_dof_index> 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<GeometryInfo<dim>::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<FullMatrix<double> >
+ interpolation_matrix (new FullMatrix<double> (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<dofs_per_cell2; ++i)
+ {
+ // if dof is locally_owned
+ const types::global_dof_index gdi = dofs[i];
+ if (u2_elements.is_element(gdi))
+ {
+ u2(dofs[i])+=u2_local(i);
+ touch_count(dofs[i]) += 1;
+ }
+ }
+ }
+ // cell1 is at the end, so should
+ // be cell2
+ Assert (cell2 == endc2, ExcInternalError());
+
+ u2.compress(VectorOperation::add);
+ touch_count.compress(VectorOperation::add);
+
+ // if we work on parallel distributed
+ // vectors, we have to ensure, that we only
+ // work on dofs this processor owns.
+ IndexSet locally_owned_dofs = dof2.locally_owned_dofs();
+
+ // when a discontinuous element is
+ // interpolated to a continuous
+ // one, we take the mean values.
+ // for parallel vectors check,
+ // if this component is owned by
+ // this processor.
+ for (types::global_dof_index i=0; i<dof2.n_dofs(); ++i)
+ if (locally_owned_dofs.is_element(i))
+ {
+ Assert(static_cast<typename OutVector::value_type>(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 <int dim,
+ template <int, int> class DoFHandlerType,
+ class InVector, class OutVector, int spacedim>
+ void
+ back_interpolate(const DoFHandlerType<dim, spacedim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
+ Vector<typename OutVector::value_type> u1_int_local(DoFTools::max_dofs_per_cell(dof1));
+
+ const types::subdomain_id subdomain_id =
+ dof1.get_triangulation().locally_owned_subdomain();
+
+ typename DoFHandlerType<dim, spacedim>::active_cell_iterator cell = dof1.begin_active(),
+ endc = dof1.end();
+
+ // map from possible fe objects in
+ // dof1 to the back_interpolation
+ // matrices
+ std::map<const FiniteElement<dim> *,
+ std::shared_ptr<FullMatrix<double> > > 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<GeometryInfo<dim>::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<FullMatrix<double> >
+ (new FullMatrix<double>(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 <int dim, int spacedim, class InVector>
+ void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &dof2,
+ const ConstraintMatrix &constraints2,
+ InVector &u1_interpolated)
+ {
+ Vector<typename InVector::value_type> 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 <int dim, int spacedim>
+ void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const PETScWrappers::MPI::Vector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, int spacedim>
+ void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const TrilinosWrappers::MPI::Vector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, int spacedim, typename Number>
+ void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const LinearAlgebra::distributed::Vector<Number> &u1,
+ const DoFHandler<dim,spacedim> &dof2,
+ const ConstraintMatrix &constraints2,
+ LinearAlgebra::distributed::Vector<Number> &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<Number>
+ 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 <int dim, class InVector, class OutVector, int spacedim>
+ void back_interpolate(const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell);
+ Vector<typename OutVector::value_type> u1_diff_local(dofs_per_cell);
+
+ const types::subdomain_id subdomain_id =
+ dof1.get_triangulation().locally_owned_subdomain();
+
+ FullMatrix<double> difference_matrix(dofs_per_cell, dofs_per_cell);
+ get_interpolation_difference_matrix(dof1.get_fe(), fe2,
+ difference_matrix);
+
+ typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::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 <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, int spacedim>
+ void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const TrilinosWrappers::MPI::Vector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void project_dg(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active();
+ typename DoFHandler<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active();
+ typename DoFHandler<dim,spacedim>::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<typename OutVector::value_type> u1_local(n1);
+ Vector<typename OutVector::value_type> u2_local(n2);
+ std::vector<types::global_dof_index> dofs(n2);
+
+ FullMatrix<double> 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<n2; ++i)
+ {
+ u2(dofs[i])+=u2_local(i);
+ }
+
+ ++cell1;
+ ++cell2;
+ }
+ }
+} // end of namespace FETools
+
+DEAL_II_NAMESPACE_CLOSE
+
+/*---------------------------- fe_tools_interpolate_templates.h ---------------------------*/
+/* end of #ifndef dealii__fe_tools_interpolate_templates_H */
+#endif
//
// ---------------------------------------------------------------------
-#include <deal.II/distributed/p4est_wrappers.h>
-#include <deal.II/distributed/tria.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/la_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/petsc_block_vector.h>
-#include <deal.II/lac/petsc_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
-#include <deal.II/lac/trilinos_block_vector.h>
-
-#include <queue>
+#include <deal.II/fe/fe_tools_extrapolate.templates.h>
DEAL_II_NAMESPACE_OPEN
-namespace FETools
-{
-
- namespace internal
- {
-#ifndef DEAL_II_WITH_P4EST
- // Dummy implementation in case p4est is not available.
- template <int dim,int spacedim,class OutVector>
- class ExtrapolateImplementation
- {
- public:
-
- ExtrapolateImplementation ()
- {
- Assert(false, ExcNotImplemented());
- };
-
- template <class InVector>
- void extrapolate_parallel (const InVector &/*u2_relevant*/,
- const DoFHandler<dim,spacedim> &/*dof2*/,
- OutVector &/*u2*/)
- {
- Assert(false, ExcNotImplemented())
- }
- };
-#else
- // Implementation of the @p extrapolate function
- // on parallel distributed grids.
- template <int dim,int spacedim,class OutVector>
- class ExtrapolateImplementation
- {
- public:
-
- ExtrapolateImplementation ();
-
- template <class InVector>
- void extrapolate_parallel (const InVector &u2_relevant,
- const DoFHandler<dim,spacedim> &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<dim>::forest forest;
- const typename dealii::internal::p4est::types<dim>::tree tree;
- const typename dealii::internal::p4est::types<dim>::locidx tree_index;
- const typename DoFHandler<dim,spacedim>::cell_iterator dealii_cell;
- const typename dealii::internal::p4est::types<dim>::quadrant p4est_cell;
-
- WorkPackage(const typename dealii::internal::p4est::types<dim>::forest &forest_,
- const typename dealii::internal::p4est::types<dim>::tree &tree_,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index_,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell_,
- const typename dealii::internal::p4est::types<dim>::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<value_type> dof_values;
-
- unsigned int tree_index;
-
- typename dealii::internal::p4est::types<dim>::quadrant quadrant;
-
- int receiver;
-
- bool operator < (const CellData &rhs) const
- {
- if (dealii::internal::p4est::functions<dim>::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<dim>::quadrant)); // quadrant
- }
-
- void pack_data (std::vector<char> &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<dim>::quadrant));
- ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant);
-
- Assert (ptr == &buffer[0]+buffer.size(),
- ExcInternalError());
- }
-
- void unpack_data (const std::vector<char> &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<dim>::quadrant));
- ptr += sizeof(typename dealii::internal::p4est::types<dim>::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 <class InVector>
- void compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
- const InVector &u);
-
- // traverse recursively over
- // the cells of this tree and
- // interpolate over patches which
- // are part of this process
- template <class InVector>
- void interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::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 <class InVector>
- void get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const InVector &u,
- Vector<value_type> &interpolated_values,
- std::vector<CellData> &new_needs);
-
- // set dof values for this
- // cell by interpolation
- void set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const Vector<value_type> &interpolated_values,
- OutVector &u);
-
- // compute all cell_data
- // needed from other processes
- // to interpolate the part of
- // this process
- void compute_needs (const DoFHandler<dim,spacedim> &dof2,
- std::vector<CellData> &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<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &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<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &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 <class InVector>
- void compute_cells (const DoFHandler<dim,spacedim> &dof2,
- const InVector &u,
- std::vector<CellData> &cells_to_compute,
- std::vector<CellData> &computed_cells,
- std::vector<CellData> &new_needs);
-
- // traverse over the tree
- // and compute cells
- template <class InVector>
- void compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const InVector &u,
- std::vector<CellData> &cells_to_compute,
- std::vector<CellData> &computed_cells,
- std::vector<CellData> &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<CellData> &cells_to_send,
- std::vector<CellData> &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<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &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<CellData> &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<CellData> &cells_list);
-
- MPI_Comm communicator;
-
- // a vector of all cells this process has
- // computed or received data
- std::vector<CellData> available_cells;
-
- // stores the indices of dofs on more refined ghosted cells along
- // with the maximum level
- std::map<types::global_dof_index, int> dofs_on_refined_neighbors;
-
- // counts the send/receive round we are in
- unsigned int round;
- };
-
- template <class OutVector>
- class ExtrapolateImplementation<1,1,OutVector>
- {
- public:
-
- ExtrapolateImplementation ()
- {
- AssertThrow(false, ExcNotImplemented())
- }
-
- template <class InVector>
- void extrapolate_parallel (const InVector &/*u2_relevant*/,
- const DoFHandler<1,1> &/*dof2*/,
- OutVector &/*u2*/)
- {}
- };
-
- template <class OutVector>
- class ExtrapolateImplementation<1,2,OutVector>
- {
- public:
-
- ExtrapolateImplementation ()
- {
- AssertThrow(false, ExcNotImplemented())
- }
-
- template <class InVector>
- void extrapolate_parallel (const InVector &/*u2_relevant*/,
- const DoFHandler<1,2> &/*dof2*/,
- OutVector &/*u2*/)
- {}
- };
-
- template <class OutVector>
- class ExtrapolateImplementation<1,3,OutVector>
- {
- public:
-
- ExtrapolateImplementation ()
- {
- AssertThrow(false, ExcNotImplemented())
- }
-
- template <class InVector>
- void extrapolate_parallel (const InVector &/*u2_relevant*/,
- const DoFHandler<1,3> &/*dof2*/,
- OutVector &/*u2*/)
- {}
- };
-
-
-
- template <int dim,int spacedim,class OutVector>
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- ExtrapolateImplementation ()
- : round(0)
- {}
-
-
-
- template <int dim,int spacedim,class OutVector>
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- CellData::CellData ()
- : tree_index (0),
- receiver (0)
- {}
-
-
-
- template <int dim,int spacedim,class OutVector>
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- CellData::CellData (const unsigned int dofs_per_cell)
- : tree_index (0),
- receiver (0)
- {
- dof_values.reinit (dofs_per_cell);
- }
-
-
-
- template <int dim,int spacedim,class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- interpolate_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const InVector &u1,
- OutVector &u2)
- {
- // check if this cell exists in the local p4est
- int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
- &p4est_cell,
- dealii::internal::p4est::functions<dim>::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<dim>::
- quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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_n<dealii_cell->n_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<dim,spacedim> &fe = dealii_cell->get_dof_handler().get_fe();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
- Vector<typename OutVector::value_type> interpolated_values(dofs_per_cell);
-
- std::vector<CellData> 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 <int dim,int spacedim, class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- get_interpolated_dof_values (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const InVector &u,
- Vector<value_type> &interpolated_values,
- std::vector<CellData> &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<dim,spacedim> &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<value_type> tmp1(dofs_per_cell);
- Vector<value_type> tmp2(dofs_per_cell);
-
- interpolated_values = 0;
- std::vector<bool> restriction_is_additive (dofs_per_cell);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- restriction_is_additive[i] = fe.restriction_is_additive(i);
-
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- bool found_child = true;
- for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
- {
- if (dealii::internal::p4est::functions<dim>::
- quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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<dofs_per_cell; ++i)
- if (tmp2(i) != value_type())
- {
- if (restriction_is_additive[i])
- interpolated_values(i) += tmp2(i);
- else
- interpolated_values(i) = tmp2(i);
- }
- }
- }
-
- if (found_child == false)
- interpolated_values = 0;
- }
- }
-
-
-
- template <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- set_dof_values_by_interpolation (const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const Vector<value_type> &local_values,
- OutVector &u)
- {
- const FiniteElement<dim,spacedim> &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<types::global_dof_index> indices(dofs_per_cell);
- dealii_cell->get_dof_indices(indices);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- // don't set this dof if there is a more refined ghosted neighbor setting this dof.
- const bool on_refined_neighbor
- = (dofs_on_refined_neighbors.find(indices[j])!=dofs_on_refined_neighbors.end());
- if (!(on_refined_neighbor && dofs_on_refined_neighbors[indices[j]]>dealii_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<value_type> tmp(dofs_per_cell);
-
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::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 <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- compute_needs (const DoFHandler<dim,spacedim> &dof2,
- std::vector<CellData> &new_needs)
- {
- const parallel::distributed::Triangulation< dim, spacedim > *tr
- = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
- (&dof2.get_tria()));
-
- Assert (tr != nullptr, ExcInternalError());
-
- typename DoFHandler<dim,spacedim>::cell_iterator
- cell=dof2.begin(0),
- endc=dof2.end(0);
- for (; cell!=endc; ++cell)
- {
- if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
- tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
- == false)
- continue;
-
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
- typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
- dealii::internal::p4est::init_coarse_quadrant<dim>(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<sc_array_t *>(&tree->quadrants),
- &p4est_coarse_cell,
- dealii::internal::p4est::functions<dim>::quadrant_compare);
-
- AssertThrow (idx == -1, ExcGridNotRefinedAtLeastOnce ());
- }
-
- traverse_tree_recursively (*tr->parallel_forest,
- *tree,
- tree_index,
- cell,
- p4est_coarse_cell,
- new_needs);
- }
- }
-
-
-
- template <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- traverse_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &new_needs)
- {
- // check if this cell exists in the local p4est
- int idx = sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
- &p4est_cell,
- dealii::internal::p4est::functions<dim>::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<dim>::
- quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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_n<dealii_cell->n_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<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
- {
- traverse_tree_recursively (forest,
- tree,
- tree_index,
- dealii_cell->child(c),
- p4est_child[c],
- new_needs);
- }
- }
- }
-
-
-
- template <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- traverse_patch_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &new_needs)
- {
- if (dealii_cell->has_children ())
- {
- Assert (dealii_cell->has_children (), ExcInternalError ());
-
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
- {
- // check if this child
- // is locally available
- // in the p4est
- if (dealii::internal::p4est::functions<dim>::
- quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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 <int dim,int spacedim,class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- compute_cells (const DoFHandler<dim,spacedim> &dof2,
- const InVector &u,
- std::vector<CellData> &cells_to_compute,
- std::vector<CellData> &computed_cells,
- std::vector<CellData> &new_needs)
- {
- const parallel::distributed::Triangulation< dim, spacedim > *tr
- = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
- (&dof2.get_tria()));
-
- Assert (tr != nullptr, ExcInternalError());
-
- // collect in a set all trees this
- // process has to compute cells on
- std::set<unsigned int> trees;
- for (typename std::vector<CellData>::const_iterator it=cells_to_compute.begin(); it!=cells_to_compute.end(); ++it)
- trees.insert (it->tree_index);
-
- typename DoFHandler<dim,spacedim>::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<dim> (tr->parallel_forest,
- tree_index)
- == false))
- continue;
-
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
- dealii::internal::p4est::init_coarse_quadrant<dim>(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 <int dim,int spacedim,class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- compute_cells_in_tree_recursively (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::tree &tree,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- const InVector &u,
- std::vector<CellData> &cells_to_compute,
- std::vector<CellData> &computed_cells,
- std::vector<CellData> &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<sc_array_t *>(&tree.quadrants),
- &p4est_cell,
- dealii::internal::p4est::functions<dim>::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<dim>::
- quadrant_overlaps_tree (const_cast<typename dealii::internal::p4est::types<dim>::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<CellData> 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<tmp.size (); ++i)
- cell_data_insert (tmp[i], new_needs);
- }
- }
-
- // search recursively over this tree
- if (p4est_has_children)
- {
- typename dealii::internal::p4est::types<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::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 <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- send_cells (const std::vector<CellData> &cells_to_send,
- std::vector<CellData> &received_cells) const
- {
- std::vector<std::vector<char> > sendbuffers (cells_to_send.size());
- std::vector<std::vector<char> >::iterator buffer = sendbuffers.begin();
- std::vector<MPI_Request> requests (cells_to_send.size());
- std::vector<unsigned int> destinations;
-
- // send data
- unsigned int idx=0;
- for (typename std::vector<CellData>::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<unsigned int> senders
- = Utilities::MPI::compute_point_to_point_communication_pattern(communicator, destinations);
-
- // receive data
- std::vector<char> receive;
- CellData cell_data;
- for (unsigned int index=0; index<senders.size (); ++index)
- {
- MPI_Status status;
- int len;
- int ierr = MPI_Probe(MPI_ANY_SOURCE, round, communicator, &status);
- AssertThrowMPI(ierr);
- ierr = MPI_Get_count(&status, MPI_BYTE, &len);
- AssertThrowMPI(ierr);
- receive.resize (len);
-
- char *buf = &receive[0];
- ierr = MPI_Recv (buf, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, communicator, &status);
- AssertThrowMPI(ierr);
-
- cell_data.unpack_data (receive);
-
- // this process has to send this
- // cell back to the sender
- // the receiver is the old sender
- cell_data.receiver = status.MPI_SOURCE;
-
- received_cells.push_back (cell_data);
- }
-
- if (requests.size () > 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 <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- add_new_need (const typename dealii::internal::p4est::types<dim>::forest &forest,
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index,
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell,
- const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell,
- std::vector<CellData> &new_needs)
- {
- const FiniteElement<dim,spacedim> &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<dim>::
- comm_find_owner (const_cast<typename dealii::internal::p4est::types<dim>::forest *> (&forest),
- tree_index,
- &p4est_cell,
- dealii_cell->level_subdomain_id ());
-
- cell_data_insert (cell_data, new_needs);
- }
-
-
-
- template <int dim,int spacedim,class OutVector>
- int
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- cell_data_search (const CellData &cell_data,
- const std::vector<CellData> &cells_list)
- {
- typename std::vector<CellData>::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 <int dim,int spacedim,class OutVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- cell_data_insert (const CellData &cell_data,
- std::vector<CellData> &cells_list)
- {
- typename std::vector<CellData>::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 <int dim,int spacedim,class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- compute_all_non_local_data (const DoFHandler<dim,spacedim> &dof2,
- const InVector &u)
- {
- std::vector<CellData> 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<received_needs.size(); ++i)
- cell_data_insert (received_needs[i], cells_to_compute);
-
- compute_cells (dof2, u, cells_to_compute, computed_cells, new_needs);
-
- // if there are no cells to compute and no new needs, stop
- ready = Utilities::MPI::sum (new_needs.size () + cells_to_compute.size (), communicator);
-
- for (typename std::vector<CellData>::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<CellData>::iterator recv=received_needs.begin();
- while (recv != received_needs.end())
- {
- if (dealii::internal::p4est::quadrant_is_equal<dim>(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<CellData>::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 <int dim,int spacedim,class OutVector>
- template <class InVector>
- void
- ExtrapolateImplementation<dim,spacedim,OutVector>::
- extrapolate_parallel (const InVector &u2_relevant,
- const DoFHandler<dim,spacedim> &dof2,
- OutVector &u2)
- {
- const parallel::distributed::Triangulation< dim, spacedim > *tr
- = (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>
- (&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<dim,spacedim> &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<types::global_dof_index> indices (dofs_per_cell);
- typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::faces_per_cell; ++face)
- if (!cell->at_boundary(face))
- {
- const typename DoFHandler<dim,spacedim>::cell_iterator neighbor = cell->neighbor(face);
- if (neighbor->level() != cell->level())
- for (unsigned int i=0; i<dofs_per_face; ++i)
- {
- const types::global_dof_index index = indices[fe.face_to_cell_index(i,face)];;
- const bool index_stored
- = (dofs_on_refined_neighbors.find(index)!=dofs_on_refined_neighbors.end());
- const unsigned int level = index_stored ?
- std::max(cell->level(), 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<WorkPackage> queue;
- {
- typename DoFHandler<dim,spacedim>::cell_iterator
- cell=dof2.begin(0),
- endc=dof2.end(0);
- for (; cell!=endc; ++cell)
- {
- if (dealii::internal::p4est::tree_exists_locally<dim> (tr->parallel_forest,
- tr->coarse_cell_to_p4est_tree_permutation[cell->index()])
- == false)
- continue;
-
- typename dealii::internal::p4est::types<dim>::quadrant p4est_coarse_cell;
- const unsigned int tree_index = tr->coarse_cell_to_p4est_tree_permutation[cell->index()];
- typename dealii::internal::p4est::types<dim>::tree *tree = tr->init_tree(cell->index());
-
- dealii::internal::p4est::init_coarse_quadrant<dim>(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<dim>::forest &forest = data.forest;
- const typename dealii::internal::p4est::types<dim>::tree &tree = data.tree;
- const typename dealii::internal::p4est::types<dim>::locidx &tree_index= data.tree_index;
- const typename DoFHandler<dim,spacedim>::cell_iterator &dealii_cell = data.dealii_cell;
- const typename dealii::internal::p4est::types<dim>::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<dim>::quadrant
- p4est_child[GeometryInfo<dim>::max_children_per_cell];
-
- dealii::internal::p4est::init_quadrant_children<dim> (p4est_cell,
- p4est_child);
-
- for (unsigned int c=0; c<GeometryInfo<dim>::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 <class VectorType, class DH>
- void reinit_distributed(const DH &dh,
- VectorType &vector)
- {
- vector.reinit(dh.n_dofs());
- }
-
-#ifdef DEAL_II_WITH_PETSC
- template <int dim, int spacedim>
- void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
- PETScWrappers::MPI::Vector &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim>
- void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
- TrilinosWrappers::MPI::Vector &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim, typename Number>
- void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
- LinearAlgebra::distributed::Vector<Number> &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <class VectorType, class DH>
- void reinit_ghosted(const DH &/*dh*/,
- VectorType &/*vector*/)
- {
- Assert(false, ExcNotImplemented());
- }
-
-#ifdef DEAL_II_WITH_PETSC
- template <int dim, int spacedim>
- void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
- PETScWrappers::MPI::Vector &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim>
- void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
- TrilinosWrappers::MPI::Vector &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, int spacedim, typename Number>
- void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
- LinearAlgebra::distributed::Vector<Number> &vector)
- {
- const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
- dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&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 <int dim, class InVector, class OutVector, int spacedim>
- void extrapolate_serial(const InVector &u3,
- const DoFHandler<dim,spacedim> &dof2,
- OutVector &u2)
- {
- const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell;
- Vector<typename OutVector::value_type> dof_values(dofs_per_cell);
-
- // then traverse grid bottom up
- for (unsigned int level=0; level<dof2.get_triangulation().n_levels()-1; ++level)
- {
- typename DoFHandler<dim,spacedim>::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_n<cell->n_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 <int dim, class InVector, class OutVector, int spacedim>
- void extrapolate(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &dof2,
- OutVector &u2)
- {
- ConstraintMatrix dummy;
- dummy.close();
- extrapolate(dof1, u1, dof2, dummy, u2);
- }
-
-
-
- template <int dim, class InVector, class OutVector, int spacedim>
- void extrapolate(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<dim,spacedim>::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<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != nullptr)
- {
- interpolate(dof1, u1, dof2, constraints, u3);
-
- OutVector u3_relevant;
- internal::reinit_ghosted(dof2, u3_relevant);
- u3_relevant = u3;
-
- internal::ExtrapolateImplementation<dim,spacedim,OutVector> 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"
// ---------------------------------------------------------------------
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/qprojector.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/la_vector.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/petsc_block_vector.h>
-#include <deal.II/lac/petsc_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
-#include <deal.II/lac/trilinos_block_vector.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/hp/dof_handler.h>
-
-
-#include <deal.II/base/index_set.h>
-
-#include <iostream>
-#include <memory>
-
+#include <deal.II/fe/fe_tools_interpolate.templates.h>
DEAL_II_NAMESPACE_OPEN
-namespace FETools
-{
- template <int dim, int spacedim,
- template <int, int> class DoFHandlerType1,
- template <int, int> class DoFHandlerType2,
- class InVector, class OutVector>
- void
- interpolate(const DoFHandlerType1<dim, spacedim> &dof1,
- const InVector &u1,
- const DoFHandlerType2<dim, spacedim> &dof2,
- OutVector &u2)
- {
- ConstraintMatrix dummy;
- dummy.close();
- interpolate(dof1, u1, dof2, dummy, u2);
- }
-
-
-
- template <int dim, int spacedim,
- template <int, int> class DoFHandlerType1,
- template <int, int> class DoFHandlerType2,
- class InVector, class OutVector>
- void
- interpolate (const DoFHandlerType1<dim, spacedim> &dof1,
- const InVector &u1,
- const DoFHandlerType2<dim, spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
- Vector<typename OutVector::value_type> 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<const FiniteElement<dim,spacedim> *,
- std::map<const FiniteElement<dim,spacedim> *,
- std::shared_ptr<FullMatrix<double> > > >
- interpolation_matrices;
-
- typename DoFHandlerType1<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active(),
- endc1 = dof1.end();
- typename DoFHandlerType2<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active(),
- endc2 = dof2.end();
- (void)endc2;
-
- std::vector<types::global_dof_index> 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<GeometryInfo<dim>::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<FullMatrix<double> >
- interpolation_matrix (new FullMatrix<double> (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<dofs_per_cell2; ++i)
- {
- // if dof is locally_owned
- const types::global_dof_index gdi = dofs[i];
- if (u2_elements.is_element(gdi))
- {
- u2(dofs[i])+=u2_local(i);
- touch_count(dofs[i]) += 1;
- }
- }
- }
- // cell1 is at the end, so should
- // be cell2
- Assert (cell2 == endc2, ExcInternalError());
-
- u2.compress(VectorOperation::add);
- touch_count.compress(VectorOperation::add);
-
- // if we work on parallel distributed
- // vectors, we have to ensure, that we only
- // work on dofs this processor owns.
- IndexSet locally_owned_dofs = dof2.locally_owned_dofs();
-
- // when a discontinuous element is
- // interpolated to a continuous
- // one, we take the mean values.
- // for parallel vectors check,
- // if this component is owned by
- // this processor.
- for (types::global_dof_index i=0; i<dof2.n_dofs(); ++i)
- if (locally_owned_dofs.is_element(i))
- {
- Assert(static_cast<typename OutVector::value_type>(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 <int dim,
- template <int, int> class DoFHandlerType,
- class InVector, class OutVector, int spacedim>
- void
- back_interpolate(const DoFHandlerType<dim, spacedim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
- Vector<typename OutVector::value_type> u1_int_local(DoFTools::max_dofs_per_cell(dof1));
-
- const types::subdomain_id subdomain_id =
- dof1.get_triangulation().locally_owned_subdomain();
-
- typename DoFHandlerType<dim, spacedim>::active_cell_iterator cell = dof1.begin_active(),
- endc = dof1.end();
-
- // map from possible fe objects in
- // dof1 to the back_interpolation
- // matrices
- std::map<const FiniteElement<dim> *,
- std::shared_ptr<FullMatrix<double> > > 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<GeometryInfo<dim>::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<FullMatrix<double> >
- (new FullMatrix<double>(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 <int dim, int spacedim, class InVector>
- void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &dof2,
- const ConstraintMatrix &constraints2,
- InVector &u1_interpolated)
- {
- Vector<typename InVector::value_type> 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 <int dim, int spacedim>
- void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const PETScWrappers::MPI::Vector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, int spacedim>
- void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const TrilinosWrappers::MPI::Vector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, int spacedim, typename Number>
- void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const LinearAlgebra::distributed::Vector<Number> &u1,
- const DoFHandler<dim,spacedim> &dof2,
- const ConstraintMatrix &constraints2,
- LinearAlgebra::distributed::Vector<Number> &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<Number>
- 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 <int dim, class InVector, class OutVector, int spacedim>
- void back_interpolate(const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell);
- Vector<typename OutVector::value_type> u1_diff_local(dofs_per_cell);
-
- const types::subdomain_id subdomain_id =
- dof1.get_triangulation().locally_owned_subdomain();
-
- FullMatrix<double> difference_matrix(dofs_per_cell, dofs_per_cell);
- get_interpolation_difference_matrix(dof1.get_fe(), fe2,
- difference_matrix);
-
- typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::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 <int dim, class InVector, class OutVector, int spacedim>
- void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, int spacedim>
- void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const TrilinosWrappers::MPI::Vector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- void project_dg(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active();
- typename DoFHandler<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active();
- typename DoFHandler<dim,spacedim>::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<typename OutVector::value_type> u1_local(n1);
- Vector<typename OutVector::value_type> u2_local(n2);
- std::vector<types::global_dof_index> dofs(n2);
-
- FullMatrix<double> 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<n2; ++i)
- {
- u2(dofs[i])+=u2_local(i);
- }
-
- ++cell1;
- ++cell2;
- }
- }
-} // end of namespace FETools
-
-
-
/*-------------- Explicit Instantiations -------------------------------*/
#include "fe_tools_interpolate.inst"
-
-/*---------------------------- fe_tools.cc ---------------------------*/
-
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/trilinos_parallel_block_vector.h>
+#include "fe_tools_extrapolate_common.h"
+
+// check FETools::extrapolate on distributed triangulations
+// for TrilinosWrappers::MPI::BlockVector
+
+template <int dim>
+void
+check (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &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<dim, TrilinosWrappers::MPI::BlockVector> (fe1, fe2);
+}
+
+#define CHECK(EL1,deg1,EL2,deg2,dim)\
+ { FE_ ## EL1<dim> fe1(deg1); \
+ FE_ ## EL2<dim> 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);
+}
+
--- /dev/null
+
+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
--- /dev/null
+
+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
+
+
+
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/petsc_parallel_block_vector.h>
+#include "fe_tools_extrapolate_common.h"
+
+// check FETools::extrapolate on distributed triangulations
+// for PETScWrappers::MPI::BlockVector
+
+template <int dim>
+void
+check (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &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<dim, PETScWrappers::MPI::BlockVector> (fe1, fe2);
+}
+
+#define CHECK(EL1,deg1,EL2,deg2,dim)\
+ { FE_ ## EL1<dim> fe1(deg1); \
+ FE_ ## EL2<dim> 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);
+}
+
--- /dev/null
+
+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
--- /dev/null
+
+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
+
+
+
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/la_parallel_block_vector.h>
+#include "fe_tools_extrapolate_common.h"
+
+// check FETools::extrapolate on distributed triangulations
+// for LinearAlgebra::distributed::BlockVector<double>
+
+template <int dim>
+void
+check (const FiniteElement<dim> &fe1,
+ const FiniteElement<dim> &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<dim, LinearAlgebra::distributed::BlockVector<double> > (fe1, fe2);
+}
+
+#define CHECK(EL1,deg1,EL2,deg2,dim)\
+ { FE_ ## EL1<dim> fe1(deg1); \
+ FE_ ## EL2<dim> 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);
+}
+
--- /dev/null
+
+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
--- /dev/null
+
+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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector_base.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
+template <typename VectorType>
+typename std::enable_if<!IsBlockVector<VectorType>::value, VectorType>::type
+build_ghosted (const IndexSet &owned_indices,
+ const IndexSet &ghosted_indices)
+{
+ return VectorType(owned_indices, ghosted_indices, MPI_COMM_WORLD);
+}
+
+template <typename VectorType>
+typename std::enable_if<IsBlockVector<VectorType>::value, VectorType>::type
+build_ghosted (const IndexSet &owned_indices,
+ const IndexSet &ghosted_indices)
+{
+ std::vector<IndexSet> owned_indices_vector(1, owned_indices);
+ std::vector<IndexSet> ghosted_indices_vector(1, ghosted_indices);
+ return VectorType(owned_indices_vector, ghosted_indices_vector, MPI_COMM_WORLD);
+}
+
+
+
+template <typename VectorType>
+typename std::enable_if<!IsBlockVector<VectorType>::value, VectorType>::type
+build_distributed (const IndexSet &owned_indices)
+{
+ return VectorType(owned_indices, MPI_COMM_WORLD);
+}
+
+template <typename VectorType>
+typename std::enable_if<IsBlockVector<VectorType>::value, VectorType>::type
+build_distributed (const IndexSet &owned_indices)
+{
+ std::vector<IndexSet> owned_indices_vector(1, owned_indices);
+ return VectorType(owned_indices_vector, MPI_COMM_WORLD);
+}
+
+
+
template <int dim, typename VectorType>
void
check_this (const FiniteElement<dim> &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<VectorType> (locally_owned_dofs1, locally_relevant_dofs1);
+ VectorType in_distributed = build_distributed<VectorType> (locally_owned_dofs1);
+ VectorType out_distributed = build_distributed<VectorType> (locally_owned_dofs2);
+ VectorType out_ghosted = build_ghosted<VectorType> (locally_owned_dofs2, locally_relevant_dofs2);
+ VectorType out_reference = build_distributed<VectorType> (locally_owned_dofs2);
// Choose some reference function of sufficiently high polynomial degree.
TestFunction<dim> function (std::max(fe1.degree,fe2.degree));
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<VectorType> (locally_owned_dofs1, locally_relevant_dofs1);
+ VectorType out_ghosted = build_ghosted<VectorType> (locally_owned_dofs2, locally_relevant_dofs2);
+ VectorType out_reference = build_ghosted<VectorType> (locally_owned_dofs2, locally_relevant_dofs2);
// Choose some reference function of sufficiently high polynomial degree.
TestFunction<dim> function (std::max(fe1.degree,fe2.degree));
projects / dealii.git / commitdiff