# Generator specific values:
IF(CMAKE_GENERATOR MATCHES "Ninja")
#
- # Ninja doesn't like a redifinition of the "help" target, so use "info"
+ # Ninja doesn't like a redefinition of the "help" target, so use "info"
# in this case...
#
SET(_make_command "$ ninja")
<h3>Specific improvements</h3>
<ol>
+ <li> Fixed: The method DoFTools::extract_constant_modes only worked for
+ elements where the constant function 1 is represented by all ones. This
+ is now fixed by querying the element for its constant modes on each cell.
+ <br>
+ (Martin Kronbichler, 2014/01/19)
+ </li>
+
+ <li> Fixed: PETScWrappers::MPI::Vector::all_zero() was broken with more than
+ one processor (illegal memory access) and did not communicate between all
+ processors. Documentation for many vector types of all_zero() has been extended.
+ <br>
+ (Timo Heister, 2014/01/17)
+ </li>
+
+ <li> Fixed/new: DoFCellAccessor::set_dof_values_by_interpolation and
+ DoFCellAccessor::get_interpolated_dof_values could previously be
+ called for hp::DoFHandler objects on cells that are non-active. This
+ makes no sense since these cells have no associated finite element
+ space. Doing so now raises an exception.
+ <br>
+ However, there are legitimate cases where one may want to interpolate
+ from children to a parent's finite element space. Since in the hp
+ case no finite element space is naturally associated with an inactive
+ cell, it is now possible to pass an explicit finite element index
+ argument to the function specifying which element of an hp::FECollection
+ object describes the space onto which you want to interpolate.
+ <br>
+ (Wolfgang Bangerth, 2014/01/18)
+ </li>
+
+ <li> Fixed: The methods IndexSet::do_compress() and
+ IndexSet::add_indices(IndexSet&) had quadratic complexity in the number of
+ ranges. The algorithms have been changed into linear complexity ones.
+ <br>
+ (Martin Kronbichler, 2014/01/15)
+ </li>
<li> Fixed: There were several bugs in functions like
FEValues::get_function_values() where the code did not properly handle the
// with a gotcha that is easily forgotten: The indices of global
// degrees of freedom we want here are the ones for current level, not
// for the global matrix. We therefore need the function
- // MGDoFAccessorLLget_mg_dof_indices, not
+ // MGDoFAccessor::get_mg_dof_indices, not
// MGDoFAccessor::get_dof_indices as used in the assembly of the
// global system:
cell->get_mg_dof_indices (local_dof_indices);
/* ---------------------------------------------------------------------
* $Id$
*
- * Copyright (C) 2008 - 2013 by the deal.II authors
+ * Copyright (C) 2008 - 2014 by the deal.II authors
*
* This file is part of the deal.II library.
*
// that have been broken out of the ones listed above. Specifically, there
// are first three functions that we call from <code>setup_dofs</code> and
// then the ones that do the assembling of linear systems:
- void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
- void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
- void setup_temperature_matrices (const IndexSet &temperature_partitioning);
+ void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning,
+ const std::vector<IndexSet> &stokes_relevant_partitioning);
+ void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning,
+ const std::vector<IndexSet> &stokes_relevant_partitioning);
+ void setup_temperature_matrices (const IndexSet &temperature_partitioning,
+ const IndexSet &temperature_relevant_partitioning);
// Following the @ref MTWorkStream "task-based parallelization" paradigm,
// the globally assembled sparsity pattern available with which the global
// matrix can be initialized.
//
+ // There is one important aspect when initializing Trilinos sparsity
+ // patterns in parallel: In addition to specifying the locally owned rows
+ // and columns of the matrices via the @p stokes_partitioning index set, we
+ // also supply information about all the rows we are possibly going to write
+ // into when assembling on a certain processor. The set of locally relevant
+ // rows contains all such rows (possibly also a few unnecessary ones, but it
+ // is difficult to find the exact row indices before actually getting
+ // indices on all cells and resolving constraints). This additional
+ // information allows to exactly determine the structure for the
+ // off-processor data found during assembly. While Trilinos matrices are
+ // able to collect this information on the fly as well (when initializing
+ // them from some other reinit method), it is less efficient and leads to
+ // problems when assembling matrices with multiple threads. In this program,
+ // we pessimistically assume that only one processor at a time can write
+ // into the matrix while assembly (whereas the computation is parallel),
+ // which is fine for Trilinos matrices. In practice, one can do better by
+ // hinting WorkStream at cells that do not share vertices, allowing for
+ // parallelism among those cells (see the graph coloring algorithms and
+ // WorkStream with colored iterators argument). However, that only works
+ // when only one MPI processor is present because Trilinos' internal data
+ // structures for accumulating off-processor data on the fly are not thread
+ // safe. With the initialization presented here, there is no such problem
+ // and one could safely introduce graph coloring for this algorithm.
+ //
// The only other change we need to make is to tell the
// DoFTools::make_sparsity_pattern() function that it is only supposed to
// work on a subset of cells, namely the ones whose
// once the matrix has been given the sparsity structure.
template <int dim>
void BoussinesqFlowProblem<dim>::
- setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning)
+ setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning,
+ const std::vector<IndexSet> &stokes_relevant_partitioning)
{
stokes_matrix.clear ();
- TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
- MPI_COMM_WORLD);
+ TrilinosWrappers::BlockSparsityPattern sp(stokes_partitioning, stokes_partitioning,
+ stokes_relevant_partitioning,
+ MPI_COMM_WORLD);
Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
-
for (unsigned int c=0; c<dim+1; ++c)
for (unsigned int d=0; d<dim+1; ++d)
if (! ((c==dim) && (d==dim)))
template <int dim>
void BoussinesqFlowProblem<dim>::
- setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning)
+ setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning,
+ const std::vector<IndexSet> &stokes_relevant_partitioning)
{
Amg_preconditioner.reset ();
Mp_preconditioner.reset ();
stokes_preconditioner_matrix.clear ();
- TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
- MPI_COMM_WORLD);
+ TrilinosWrappers::BlockSparsityPattern sp(stokes_partitioning, stokes_partitioning,
+ stokes_relevant_partitioning,
+ MPI_COMM_WORLD);
Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
for (unsigned int c=0; c<dim+1; ++c)
template <int dim>
void BoussinesqFlowProblem<dim>::
- setup_temperature_matrices (const IndexSet &temperature_partitioner)
+ setup_temperature_matrices (const IndexSet &temperature_partitioner,
+ const IndexSet &temperature_relevant_partitioner)
{
T_preconditioner.reset ();
temperature_mass_matrix.clear ();
temperature_stiffness_matrix.clear ();
temperature_matrix.clear ();
- TrilinosWrappers::SparsityPattern sp (temperature_partitioner,
- MPI_COMM_WORLD);
+ TrilinosWrappers::SparsityPattern sp(temperature_partitioner,
+ temperature_partitioner,
+ temperature_relevant_partitioner,
+ MPI_COMM_WORLD);
DoFTools::make_sparsity_pattern (temperature_dof_handler, sp,
temperature_constraints, false,
Utilities::MPI::
// All this done, we can then initialize the various matrix and vector
// objects to their proper sizes. At the end, we also record that all
// matrices and preconditioners have to be re-computed at the beginning of
- // the next time step.
- setup_stokes_matrix (stokes_partitioning);
- setup_stokes_preconditioner (stokes_partitioning);
- setup_temperature_matrices (temperature_partitioning);
-
- stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD);
+ // the next time step. Note how we initialize the vectors for the Stokes
+ // and temperature right hand sides: These are writable vectors (last
+ // boolean argument set to @p true) that have the correct one-to-one
+ // partitioning of locally owned elements but are still given the relevant
+ // partitioning for means of figuring out the vector entries that are
+ // going to be set right away. As for matrices, this allows for writing
+ // local contributions into the vector with multiple threads (always
+ // assuming that the same vector entry is not accessed by multiple threads
+ // at the same time). The other vectors only allow for read access of
+ // individual elements, including ghosts, but are not suitable for
+ // solvers.
+ setup_stokes_matrix (stokes_partitioning, stokes_relevant_partitioning);
+ setup_stokes_preconditioner (stokes_partitioning,
+ stokes_relevant_partitioning);
+ setup_temperature_matrices (temperature_partitioning,
+ temperature_relevant_partitioning);
+
+ stokes_rhs.reinit (stokes_partitioning, stokes_relevant_partitioning,
+ MPI_COMM_WORLD, true);
stokes_solution.reinit (stokes_relevant_partitioning, MPI_COMM_WORLD);
old_stokes_solution.reinit (stokes_solution);
- temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD);
+ temperature_rhs.reinit (temperature_partitioning,
+ temperature_relevant_partitioning,
+ MPI_COMM_WORLD, true);
temperature_solution.reinit (temperature_relevant_partitioning, MPI_COMM_WORLD);
old_temperature_solution.reinit (temperature_solution);
old_old_temperature_solution.reinit (temperature_solution);
Assembly::CopyData::
StokesSystem<dim> (stokes_fe));
- stokes_matrix.compress(VectorOperation::add);
+ if (rebuild_stokes_matrix == true)
+ stokes_matrix.compress(VectorOperation::add);
stokes_rhs.compress(VectorOperation::add);
rebuild_stokes_matrix = false;
using namespace Step32;
using namespace dealii;
- Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
+ numbers::invalid_unsigned_int);
try
{
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1999 - 2013 by the deal.II authors
+// Copyright (C) 1999 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
* of an STL <tt>std::map</tt> type with a certain number of
* elements?), this is only an estimate. however often quite close
* to the true value.
- */
+ */
std::size_t memory_consumption () const;
-
- private:
-
};
* understand.
*/
static const unsigned int format_version = 3;
+
private:
/**
* Dummy entry to suppress compiler warnings when copying an empty
using DataOutBase::get_output_format_names;
using DataOutBase::determine_intermediate_format_dimensions;
+ using DataOutBase::Patch;
+
/**
* Constructor.
*/
IndexSet ();
/**
- * Constructor that also sets the
- * overall size of the index
- * range.
+ * Constructor that also sets the overall size of the index range.
*/
explicit IndexSet (const types::global_dof_index size);
/**
- * Remove all indices from this
- * index set. The index set retains
- * its size, however.
+ * Remove all indices from this index set. The index set retains its size,
+ * however.
*/
void clear ();
/**
- * Set the maximal size of the
- * indices upon which this object
- * operates.
+ * Set the maximal size of the indices upon which this object operates.
*
- * This function can only be
- * called if the index set does
- * not yet contain any elements.
- * This can be achieved by calling
- * clear(), for example.
+ * This function can only be called if the index set does not yet contain
+ * any elements. This can be achieved by calling clear(), for example.
*/
void set_size (const types::global_dof_index size);
/**
- * Return the size of the index
- * space of which this index set
- * is a subset of.
+ * Return the size of the index space of which this index set is a subset
+ * of.
*
- * Note that the result is not equal to
- * the number of indices within this
- * set. The latter information is
- * returned by n_elements().
+ * Note that the result is not equal to the number of indices within this
+ * set. The latter information is returned by n_elements().
*/
types::global_dof_index size () const;
/**
- * Add the half-open range
- * $[\text{begin},\text{end})$ to
- * the set of indices represented
- * by this class.
+ * Add the half-open range $[\text{begin},\text{end})$ to the set of indices
+ * represented by this class.
*/
void add_range (const types::global_dof_index begin,
const types::global_dof_index end);
/**
- * Add an individual index to the
- * set of indices.
+ * Add an individual index to the set of indices.
*/
void add_index (const types::global_dof_index index);
/**
- * Add a whole set of indices
- * described by dereferencing
- * every element of the the
- * iterator range
- * <code>[begin,end)</code>.
+ * Add a whole set of indices described by dereferencing every element of
+ * the the iterator range <code>[begin,end)</code>.
*/
template <typename ForwardIterator>
void add_indices (const ForwardIterator &begin,
const ForwardIterator &end);
/**
- * Add the given IndexSet @p other to the
- * current one, constructing the union of
- * *this and @p other.
+ * Add the given IndexSet @p other to the current one, constructing the
+ * union of *this and @p other.
*
- * If the @p offset argument is nonzero, then every
- * index in @p other is shifted by @p offset before being
- * added to the current index set. This allows to construct,
- * for example, one index set from several others that are
- * supposed to represent index sets corresponding to
- * different ranges (e.g., when constructing the set of
- * nonzero entries of a block vector from the sets of nonzero
- * elements of the individual blocks of a vector).
+ * If the @p offset argument is nonzero, then every index in @p other is
+ * shifted by @p offset before being added to the current index set. This
+ * allows to construct, for example, one index set from several others that
+ * are supposed to represent index sets corresponding to different ranges
+ * (e.g., when constructing the set of nonzero entries of a block vector
+ * from the sets of nonzero elements of the individual blocks of a vector).
*
- * This function will generate an exception if any of the
- * (possibly shifted) indices of the @p other index set
- * lie outside the range <code>[0,size())</code> represented
- * by the current object.
+ * This function will generate an exception if any of the (possibly shifted)
+ * indices of the @p other index set lie outside the range
+ * <code>[0,size())</code> represented by the current object.
*/
void add_indices(const IndexSet &other,
const unsigned int offset = 0);
/**
- * Return whether the specified
- * index is an element of the
- * index set.
+ * Return whether the specified index is an element of the index set.
*/
bool is_element (const types::global_dof_index index) const;
/**
- * Return whether the index set
- * stored by this object defines
- * a contiguous range. This is
- * true also if no indices are
- * stored at all.
+ * Return whether the index set stored by this object defines a contiguous
+ * range. This is true also if no indices are stored at all.
*/
bool is_contiguous () const;
/**
- * Return the number of elements
- * stored in this index set.
+ * Return the number of elements stored in this index set.
*/
types::global_dof_index n_elements () const;
/**
- * Return the global index of the local
- * index with number @p local_index
- * stored in this index set. @p
- * local_index obviously needs to be less
- * than n_elements().
+ * Return the global index of the local index with number @p local_index
+ * stored in this index set. @p local_index obviously needs to be less than
+ * n_elements().
*/
types::global_dof_index nth_index_in_set (const unsigned int local_index) const;
/**
- * Return the how-manyth element of this
- * set (counted in ascending order) @p
- * global_index is. @p global_index needs
- * to be less than the size(). This
- * function throws an exception if the
- * index @p global_index is not actually
- * a member of this index set, i.e. if
- * is_element(global_index) is false.
+ * Return the how-manyth element of this set (counted in ascending order) @p
+ * global_index is. @p global_index needs to be less than the size(). This
+ * function throws an exception if the index @p global_index is not actually
+ * a member of this index set, i.e. if is_element(global_index) is false.
*/
types::global_dof_index index_within_set (const types::global_dof_index global_index) const;
/**
- * Each index set can be
- * represented as the union of a
- * number of contiguous intervals
- * of indices, where if necessary
- * intervals may only consist of
- * individual elements to
- * represent isolated members of
- * the index set.
+ * Each index set can be represented as the union of a number of contiguous
+ * intervals of indices, where if necessary intervals may only consist of
+ * individual elements to represent isolated members of the index set.
*
- * This function returns the
- * minimal number of such
- * intervals that are needed to
- * represent the index set under
- * consideration.
+ * This function returns the minimal number of such intervals that are
+ * needed to represent the index set under consideration.
*/
unsigned int n_intervals () const;
/**
- * Compress the internal
- * representation by merging
- * individual elements with
- * contiguous ranges, etc. This
- * function does not have any
- * external effect.
+ * Compress the internal representation by merging individual elements with
+ * contiguous ranges, etc. This function does not have any external effect.
*/
void compress () const;
/**
- * Comparison for equality of
- * index sets. This operation is
- * only allowed if the size of
- * the two sets is the same
- * (though of course they do not
- * have to have the same number
- * of indices).
+ * Comparison for equality of index sets. This operation is only allowed if
+ * the size of the two sets is the same (though of course they do not have
+ * to have the same number of indices).
*/
bool operator == (const IndexSet &is) const;
/**
- * Comparison for inequality of
- * index sets. This operation is
- * only allowed if the size of
- * the two sets is the same
- * (though of course they do not
- * have to have the same number
- * of indices).
+ * Comparison for inequality of index sets. This operation is only allowed
+ * if the size of the two sets is the same (though of course they do not
+ * have to have the same number of indices).
*/
bool operator != (const IndexSet &is) const;
/**
- * Return the intersection of the
- * current index set and the
- * argument given, i.e. a set of
- * indices that are elements of
- * both index sets. The two index
- * sets must have the same size
- * (though of course they do not
- * have to have the same number
- * of indices).
+ * Return the intersection of the current index set and the argument given,
+ * i.e. a set of indices that are elements of both index sets. The two index
+ * sets must have the same size (though of course they do not have to have
+ * the same number of indices).
*/
IndexSet operator & (const IndexSet &is) const;
/**
- * This command takes an interval
- * <tt>[begin, end)</tt> and returns
- * the intersection of the current
- * index set with the interval, shifted
- * to the range <tt>[0,
- * end-begin)</tt>.
+ * This command takes an interval <tt>[begin, end)</tt> and returns the
+ * intersection of the current index set with the interval, shifted to the
+ * range <tt>[0, end-begin)</tt>.
*
- * In other words, the result of this operation is the
- * intersection of the set represented by the current object
- * and the interval <tt>[begin, end)</tt>, as seen
- * <i>within the interval <tt>[begin, end)</tt></i> by
- * shifting the result of the intersection operation to
- * the left by <tt>begin</tt>. This corresponds to the notion
- * of a <i>view</i>: The interval <tt>[begin, end)</tt> is
- * a <i>window</i> through which we see the set represented
- * by the current object.
+ * In other words, the result of this operation is the intersection of the
+ * set represented by the current object and the interval <tt>[begin,
+ * end)</tt>, as seen <i>within the interval <tt>[begin, end)</tt></i> by
+ * shifting the result of the intersection operation to the left by
+ * <tt>begin</tt>. This corresponds to the notion of a <i>view</i>: The
+ * interval <tt>[begin, end)</tt> is a <i>window</i> through which we see
+ * the set represented by the current object.
*/
IndexSet get_view (const types::global_dof_index begin,
const types::global_dof_index end) const;
/**
- * Removes all elements contained in @p
- * other from this set. In other words,
- * if $x$ is the current object and $o$
- * the argument, then we compute $x
+ * Removes all elements contained in @p other from this set. In other words,
+ * if $x$ is the current object and $o$ the argument, then we compute $x
* \leftarrow x \backslash o$.
*/
void subtract_set (const IndexSet &other);
/**
- * Fills the given vector with all
- * indices contained in this IndexSet.
+ * Fills the given vector with all indices contained in this IndexSet.
*/
void fill_index_vector(std::vector<types::global_dof_index> &indices) const;
/**
- * Fill the given vector with either
- * zero or one elements, providing
- * a binary representation of this
- * index set. The given vector is
- * assumed to already have the correct
- * size.
+ * Fill the given vector with either zero or one elements, providing a
+ * binary representation of this index set. The given vector is assumed to
+ * already have the correct size.
*
- * The given argument is filled with
- * integer values zero and one, using
- * <code>vector.operator[]</code>. Thus,
- * any object that has such an operator
- * can be used as long as it allows
- * conversion of integers zero and one to
- * elements of the vector. Specifically,
- * this is the case for classes Vector,
- * BlockVector, but also
- * std::vector@<bool@>,
- * std::vector@<int@>, and
- * std::vector@<double@>.
+ * The given argument is filled with integer values zero and one, using
+ * <code>vector.operator[]</code>. Thus, any object that has such an
+ * operator can be used as long as it allows conversion of integers zero and
+ * one to elements of the vector. Specifically, this is the case for classes
+ * Vector, BlockVector, but also std::vector@<bool@>, std::vector@<int@>,
+ * and std::vector@<double@>.
*/
template <typename Vector>
void fill_binary_vector (Vector &vector) const;
/**
- * Outputs a text representation of this
- * IndexSet to the given stream. Used for
- * testing.
+ * Outputs a text representation of this IndexSet to the given stream. Used
+ * for testing.
*/
template <class STREAM>
void print(STREAM &out) const;
/**
- * Writes the IndexSet into a text based
- * file format, that can be read in again
- * using the read() function.
+ * Writes the IndexSet into a text based file format, that can be read in
+ * again using the read() function.
*/
void write(std::ostream &out) const;
/**
- * Constructs the IndexSet from a text
- * based representation given by the
- * stream @param in written by the
- * write() function.
+ * Constructs the IndexSet from a text based representation given by the
+ * stream @param in written by the write() function.
*/
void read(std::istream &in);
/**
- * Writes the IndexSet into a binary,
- * compact representation, that can be
- * read in again using the block_read()
- * function.
+ * Writes the IndexSet into a binary, compact representation, that can be
+ * read in again using the block_read() function.
*/
void block_write(std::ostream &out) const;
/**
- * Constructs the IndexSet from a binary
- * representation given by the stream
- * @param in written by the write_block()
- * function.
+ * Constructs the IndexSet from a binary representation given by the stream
+ * @param in written by the write_block() function.
*/
void block_read(std::istream &in);
#ifdef DEAL_II_WITH_TRILINOS
/**
- * Given an MPI communicator,
- * create a Trilinos map object
- * that represents a distribution
- * of vector elements or matrix
- * rows in which we will locally
- * store those elements or rows
- * for which we store the index
- * in the current index set, and
- * all the other elements/rows
- * elsewhere on one of the other
+ * Given an MPI communicator, create a Trilinos map object that represents a
+ * distribution of vector elements or matrix rows in which we will locally
+ * store those elements or rows for which we store the index in the current
+ * index set, and all the other elements/rows elsewhere on one of the other
* MPI processes.
*
- * The last argument only plays a
- * role if the communicator is a
- * parallel one, distributing
- * computations across multiple
- * processors. In that case, if
- * the last argument is false,
- * then it is assumed that the
- * index sets this function is
- * called on on all processors
- * are mutually exclusive but
- * together enumerate each index
- * exactly once. In other words,
- * if you call this function on
- * two processors, then the index
- * sets this function is called
- * with must together have all
- * possible indices from zero to
- * size()-1, and no index must
- * appear in both index
- * sets. This corresponds, for
- * example, to the case where we
- * want to split the elements of
- * vectors into unique subsets to
- * be stored on different
- * processors -- no element
- * should be owned by more than
- * one processor, but each
- * element must be owned by one.
+ * The last argument only plays a role if the communicator is a parallel
+ * one, distributing computations across multiple processors. In that case,
+ * if the last argument is false, then it is assumed that the index sets
+ * this function is called on on all processors are mutually exclusive but
+ * together enumerate each index exactly once. In other words, if you call
+ * this function on two processors, then the index sets this function is
+ * called with must together have all possible indices from zero to
+ * size()-1, and no index must appear in both index sets. This corresponds,
+ * for example, to the case where we want to split the elements of vectors
+ * into unique subsets to be stored on different processors -- no element
+ * should be owned by more than one processor, but each element must be
+ * owned by one.
*
- * On the other hand, if the
- * second argument is true, then
- * the index sets can be
- * overlapping, though they still
- * need to contain each index
- * exactly once on all processors
- * taken together. This is a
- * useful operation if we want to
- * create vectors that not only
- * contain the locally owned
- * indices, but for example also
- * the elements that correspond
- * to degrees of freedom located
- * on ghost cells.
+ * On the other hand, if the second argument is true, then the index sets
+ * can be overlapping, though they still need to contain each index exactly
+ * once on all processors taken together. This is a useful operation if we
+ * want to create vectors that not only contain the locally owned indices,
+ * but for example also the elements that correspond to degrees of freedom
+ * located on ghost cells.
*/
Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD,
const bool overlapping = false) const;
/**
- * Determine an estimate for the memory
- * consumption (in bytes) of this
+ * Determine an estimate for the memory consumption (in bytes) of this
* object.
*/
std::size_t memory_consumption () const;
private:
/**
- * A type that denotes the half
- * open index range
- * <code>[begin,end)</code>.
+ * A type that denotes the half open index range <code>[begin,end)</code>.
*
- * The nth_index_in_set denotes
- * the how many-th index within
- * this IndexSet the first
- * element of the current range
- * is. This information is only
- * accurate if
- * IndexSet::compress() has been
- * called after the last
+ * The nth_index_in_set denotes the how many-th index within this IndexSet
+ * the first element of the current range is. This information is only
+ * accurate if IndexSet::compress() has been called after the last
* insertion.
*/
struct Range
}
/**
- * Write or read the data of this object to or
- * from a stream for the purpose of serialization
+ * Write or read the data of this object to or from a stream for the
+ * purpose of serialization
*/
template <class Archive>
void serialize (Archive &ar, const unsigned int version);
};
/**
- * A set of contiguous ranges of
- * indices that make up (part of)
- * this index set. This variable
- * is always kept sorted.
+ * A set of contiguous ranges of indices that make up (part of) this index
+ * set. This variable is always kept sorted.
*
- * The variable is marked
- * "mutable" so that it can be
- * changed by compress(), though
- * this of course doesn't change
- * anything about the external
- * representation of this index
- * set.
+ * The variable is marked "mutable" so that it can be changed by compress(),
+ * though this of course doesn't change anything about the external
+ * representation of this index set.
*/
mutable std::vector<Range> ranges;
/**
- * True if compress() has been
- * called after the last change
- * in the set of indices.
+ * True if compress() has been called after the last change in the set of
+ * indices.
*
- * The variable is marked
- * "mutable" so that it can be
- * changed by compress(), though
- * this of course doesn't change
- * anything about the external
- * representation of this index
- * set.
+ * The variable is marked "mutable" so that it can be changed by compress(),
+ * though this of course doesn't change anything about the external
+ * representation of this index set.
*/
mutable bool is_compressed;
/**
- * The overall size of the index
- * range. Elements of this index
- * set have to have a smaller
- * number than this value.
+ * The overall size of the index range. Elements of this index set have to
+ * have a smaller number than this value.
*/
types::global_dof_index index_space_size;
/**
- * This integer caches the index of the
- * largest range in @p ranges. This gives
- * <tt>O(1)</tt> access to the range with
- * most elements, while general access
- * costs <tt>O(log(n_ranges))</tt>. The
- * largest range is needed for the
- * methods @p is_element(), @p
- * index_within_set(), @p
- * nth_index_in_set. In many
- * applications, the largest range
- * contains most elements (the locally
- * owned range), whereas there are only a
- * few other elements (ghosts).
+ * This integer caches the index of the largest range in @p ranges. This
+ * gives <tt>O(1)</tt> access to the range with most elements, while general
+ * access costs <tt>O(log(n_ranges))</tt>. The largest range is needed for
+ * the methods @p is_element(), @p index_within_set(), @p
+ * nth_index_in_set. In many applications, the largest range contains most
+ * elements (the locally owned range), whereas there are only a few other
+ * elements (ghosts).
*/
mutable types::global_dof_index largest_range;
/**
- * Actually perform the compress()
- * operation.
+ * Actually perform the compress() operation.
*/
void do_compress() const;
};
-inline
-void
-IndexSet::add_range (const types::global_dof_index begin,
- const types::global_dof_index end)
-{
- Assert ((begin < index_space_size)
- ||
- ((begin == index_space_size) && (end == index_space_size)),
- ExcIndexRangeType<types::global_dof_index> (begin, 0, index_space_size));
- Assert (end <= index_space_size,
- ExcIndexRangeType<types::global_dof_index> (end, 0, index_space_size+1));
- Assert (begin <= end,
- ExcIndexRangeType<types::global_dof_index> (begin, 0, end));
-
- if (begin != end)
- {
- const Range new_range(begin,end);
-
- // the new index might be larger than the last
- // index present in the ranges. Then we can
- // skip the binary search
- if (ranges.size() == 0 || begin > ranges.back().end)
- ranges.push_back(new_range);
- else
- ranges.insert (Utilities::lower_bound (ranges.begin(),
- ranges.end(),
- new_range),
- new_range);
- is_compressed = false;
- }
-}
-
-
-
inline
void
IndexSet::add_index (const types::global_dof_index index)
IndexSet::add_indices (const ForwardIterator &begin,
const ForwardIterator &end)
{
- // insert each element of the
- // range. if some of them happen to
- // be consecutive, merge them to a
- // range
+ // insert each element of the range. if some of them happen to be
+ // consecutive, merge them to a range
for (ForwardIterator p=begin; p!=end;)
{
const types::global_dof_index begin_index = *p;
-inline
-void
-IndexSet::add_indices(const IndexSet &other,
- const unsigned int offset)
-{
- if ((this == &other) && (offset == 0))
- return;
-
- for (std::vector<Range>::iterator range = other.ranges.begin();
- range != other.ranges.end();
- ++range)
- {
- add_range(range->begin+offset, range->end+offset);
- }
-
- compress();
-}
-
-
-
inline
bool
IndexSet::is_element (const types::global_dof_index index) const
{
compress ();
- // fast check whether the index is in the
- // largest range
+ // fast check whether the index is in the largest range
Assert (largest_range < ranges.size(), ExcInternalError());
if (index >= ranges[largest_range].begin &&
index < ranges[largest_range].end)
return true;
- // get the element after which
- // we would have to insert a
- // range that consists of all
- // elements from this element
- // to the end of the index
- // range plus one. after this
- // call we know that if
- // p!=end() then
- // p->begin<=index unless there
- // is no such range at all
+ // get the element after which we would have to insert a range that
+ // consists of all elements from this element to the end of the index
+ // range plus one. after this call we know that if p!=end() then
+ // p->begin<=index unless there is no such range at all
//
- // if the searched for element
- // is an element of this range,
- // then we're done. otherwise,
- // the element can't be in one
- // of the following ranges
- // because otherwise p would be
- // a different iterator
+ // if the searched for element is an element of this range, then we're
+ // done. otherwise, the element can't be in one of the following ranges
+ // because otherwise p would be a different iterator
//
- // since we already know the position
- // relative to the largest range (we
- // called compress!), we can perform
- // the binary search on ranges with
- // lower/higher number compared to the
- // largest range
+ // since we already know the position relative to the largest range (we
+ // called compress!), we can perform the binary search on ranges with
+ // lower/higher number compared to the largest range
std::vector<Range>::const_iterator
p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
0 : largest_range+1),
Assert ((p == ranges.end()) || (p->begin > index),
ExcInternalError());
- // now move to that previous
- // range
+ // now move to that previous range
--p;
Assert (p->begin <= index, ExcInternalError());
return (p->end > index);
}
- // didn't find this index, so it's
- // not in the set
+ // didn't find this index, so it's not in the set
return false;
}
types::global_dof_index
IndexSet::n_elements () const
{
- // make sure we have
- // non-overlapping ranges
+ // make sure we have non-overlapping ranges
compress ();
types::global_dof_index v = 0;
+inline
+unsigned int
+IndexSet::n_intervals () const
+{
+ compress ();
+ return ranges.size();
+}
+
+
+
inline
types::global_dof_index
IndexSet::nth_index_in_set (const unsigned int n) const
{
- // to make this call thread-safe, compress()
- // must not be called through this function
+ // to make this call thread-safe, compress() must not be called through this
+ // function
Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
Assert (n < n_elements(), ExcIndexRangeType<types::global_dof_index> (n, 0, n_elements()));
- // first check whether the index is in the
- // largest range
+ // first check whether the index is in the largest range
Assert (largest_range < ranges.size(), ExcInternalError());
std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
if (n>=main_range->nth_index_in_set &&
n<main_range->nth_index_in_set+(main_range->end-main_range->begin))
return main_range->begin + (n-main_range->nth_index_in_set);
- // find out which chunk the local index n
- // belongs to by using a binary search. the
- // comparator is based on the end of the
- // ranges. Use the position relative to main_range to
- // subdivide the ranges
+ // find out which chunk the local index n belongs to by using a binary
+ // search. the comparator is based on the end of the ranges. Use the
+ // position relative to main_range to subdivide the ranges
Range r (n,n+1);
r.nth_index_in_set = n;
std::vector<Range>::const_iterator range_begin, range_end;
types::global_dof_index
IndexSet::index_within_set (const types::global_dof_index n) const
{
- // to make this call thread-safe, compress()
- // must not be called through this function
+ // to make this call thread-safe, compress() must not be called through this
+ // function
Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
Assert (is_element(n) == true, ExcIndexNotPresent (n));
Assert (n < size(), ExcIndexRangeType<types::global_dof_index> (n, 0, size()));
- // check whether the index is in the largest
- // range. use the result to perform a
- // one-sided binary search afterward
+ // check whether the index is in the largest range. use the result to
+ // perform a one-sided binary search afterward
Assert (largest_range < ranges.size(), ExcInternalError());
std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
if (n >= main_range->begin && n < main_range->end)
ExcDimensionMismatch (vector.size(), size()));
compress();
- // first fill all elements of the vector
- // with zeroes.
+ // first fill all elements of the vector with zeroes.
std::fill (vector.begin(), vector.end(), 0);
- // then write ones into the elements whose
- // indices are contained in the index set
+ // then write ones into the elements whose indices are contained in the
+ // index set
for (std::vector<Range>::iterator it = ranges.begin();
it != ranges.end();
++it)
/**
* Return the interpolation of the given finite element function to
* the present cell. In the simplest case, the cell is a terminal
- * one, i.e. has no children; then, the returned value is the vector
- * of nodal values on that cell. You could then as well get the
+ * one, i.e., it has no children; then, the returned value is the vector
+ * of nodal values on that cell. You could as well get the
* desired values through the @p get_dof_values function. In the
* other case, when the cell has children, we use the restriction
* matrices provided by the finite element class to compute the
* interpolation from the children to the present cell.
*
- * It is assumed that both vectors already have the right size
+ * If the cell is part of a hp::DoFHandler object, cells only have an
+ * associated finite element space if they are active. However, this
+ * function is supposed to also provide information on inactive cells
+ * with children. Consequently, it carries a third argument that can be
+ * used in the hp context that denotes the finite element space we are
+ * supposed to interpolate onto. If the cell is active, this function
+ * then obtains the finite element function from the <code>values</code>
+ * vector on this cell and interpolates it onto the space described
+ * by the <code>fe_index</code>th element of the hp::FECollection associated
+ * with the hp::DoFHandler of which this cell is a part of. If the cell
+ * is not active, then we first perform this interpolation on all of its
+ * terminal children and then interpolate this function down to the
+ * cell requested keeping the function space the same.
+ *
+ * It is assumed that both input vectors already have the right size
* beforehand.
*
- * Unlike the get_dof_values() function, this function works on
- * cells rather than to lines, quads, and hexes, since interpolation
+ * @note Unlike the get_dof_values() function, this function is only available
+ * on cells, rather than on lines, quads, and hexes, since interpolation
* is presently only provided for cells by the finite element
* classes.
*/
template <class InputVector, typename number>
void get_interpolated_dof_values (const InputVector &values,
- Vector<number> &interpolated_values) const;
+ Vector<number> &interpolated_values,
+ const unsigned int fe_index = DH::default_fe_index) const;
/**
- * This, again, is the counterpart to get_interpolated_dof_values():
+ * This function is the counterpart to get_interpolated_dof_values():
* you specify the dof values on a cell and these are interpolated
* to the children of the present cell and set on the terminal
* cells.
*/
template <class OutputVector, typename number>
void set_dof_values_by_interpolation (const Vector<number> &local_values,
- OutputVector &values) const;
+ OutputVector &values,
+ const unsigned int fe_index = DH::default_fe_index) const;
/**
* Distribute a local (cell based) vector to a global one by mapping
: BaseClass(other),
dof_handler(0)
{
- Assert (false, ExcInvalidObject());
+ Assert (false, ExcMessage("You are trying to assign iterators that are incompatible. "
+ "Reasons for incompatibility are that they point to different "
+ "types of DoFHandlers (e.g., dealii::DoFHandler and "
+ "dealii::hp::DoFHandler) or that the refer to objects of "
+ "different dimensionality (e.g., assigning a line iterator "
+ "to a quad iterator)."));
}
BlockMask
block_mask (const ComponentMask &component_mask) const;
+ /**
+ * Returns a list of constant modes of the element. The returns table has as
+ * many rows as there are components in the element and dofs_per_cell
+ * columns. To each component of the finite element, the row in the returned
+ * table contains a basis representation of the constant function 1 on the
+ * element.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
//@}
/**
static const unsigned int n_projection_matrices;
};
+ /**
+ * Returns a list of constant modes of the element. For this element, the
+ * first entry is true, all other are false.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
protected:
/**
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
+ /**
+ * Returns a list of constant modes of the element. For this element, it
+ * simply returns one row with all entries set to true.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
/**
* Determine an estimate for the
* memory consumption (in bytes)
FiniteElementDomination::Domination
compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+
+ /**
+ * Returns a list of constant modes of the element. For this element, it
+ * simply returns one row with all entries set to true.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
private:
/**
* Return vector with dofs per vertex, line, quad, hex.
FiniteElementDomination::Domination
compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+ /**
+ * Returns a list of constant modes of the element. For this element, the
+ * first entry on each face is true, all other are false (as the constant
+ * function is represented by the first base function of Legendre
+ * polynomials.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
private:
/**
* Return vector with dofs per vertex, line, quad, hex.
virtual void interpolate (std::vector<double> &local_dofs,
const VectorSlice<const std::vector<std::vector<double> > > &values)
const;
+
+ /**
+ * Returns a list of constant modes of the element.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
virtual std::size_t memory_consumption () const;
virtual FiniteElement<dim> *clone() const;
virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
const Point<dim> &p,
const unsigned int component) const;
+
protected:
virtual
const bool face_flip = false,
const bool face_rotation = false) const;
+ /**
+ * Returns a list of constant modes of the element. For this element, the
+ * list consists of true arguments for all components.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
/**
* @name Functions to support hp
* @{
*/
std::vector<unsigned int> get_embedding_dofs (const unsigned int sub_degree) const;
+ /**
+ * Returns a list of constant modes of the element. For this element, the
+ * list consists of true arguments for the first vertex shape functions and
+ * false for the remaining ones.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
protected:
/**
* @p clone function instead of
virtual void interpolate(
std::vector<double> &local_dofs,
const VectorSlice<const std::vector<std::vector<double> > > &values) const;
+
+ /**
+ * Returns a list of constant modes of the element. For this element, the
+ * list consists of true arguments for all components.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
virtual std::size_t memory_consumption () const;
virtual FiniteElement<dim> *clone() const;
Point<dim-1>
unit_face_support_point (const unsigned int index) const;
+ /**
+ * Returns a list of constant modes of the element. The returns table has as
+ * many rows as there are components in the element and dofs_per_cell
+ * columns. To each component of the finite element, the row in the returned
+ * table contains a basis representation of the constant function 1 on the
+ * element. Concatenates the constant modes of each base element.
+ */
+ virtual Table<2,bool> get_constant_modes () const;
+
/**
* @name Functions to support hp
* @{
const double half_length = 1.0);
/**
- * Initialize the given triangulation with a hyper-L consisting of exactly
+ * Initialize the given triangulation with a hyper-L (in 2d or 3d)
+ * consisting of exactly
* <tt>2^dim-1</tt> cells. It produces the hypercube with the interval
* [<i>left,right</i>] without the hypercube made out of the interval
- * [<i>(a+b)/2,b</i>].
+ * [<i>(a+b)/2,b</i>]. This will result in the classical L-shape in 2d.
+ * The shape will look like the following in 3d:
*
* @image html hyper_l.png
*
* indicator 1, while the inner boundary has id zero. If the flag is @p
* false, both have indicator zero.
*
- * In 2D, the number <tt>n_cells</tt> of elements for this initial
+ * In 2d, the number <tt>n_cells</tt> of elements for this initial
* triangulation can be chosen arbitrarily. If the number of initial cells
* is zero (as is the default), then it is computed adaptively such that the
* resulting elements have the least aspect ratio.
*
- * In 3D, only two different numbers are meaningful, 6 for a surface based
+ * In 3d, only two different numbers are meaningful, 6 for a surface based
* on a hexahedron (i.e. 6 panels on the inner sphere extruded in radial
* direction to form 6 cells) and 12 for the rhombic dodecahedron. These
* give rise to the following meshes upon one refinement:
const bool colorize = false);
/**
- * Produce a ring of cells in 3D that is cut open, twisted and glued
+ * Produce a ring of cells in 3d that is cut open, twisted and glued
* together again. This results in a kind of moebius-loop.
*
* @param tria The triangulation to be worked on.
typedef types::global_dof_index size_type;
/**
- * An iterator that can be used to
- * iterate over the elements of a single
- * row. The result of dereferencing such
- * an iterator is a column index.
+ * An iterator that can be used to iterate over the elements of a single
+ * row. The result of dereferencing such an iterator is a column index.
*/
typedef std::vector<size_type>::const_iterator row_iterator;
/**
- * Initialize the matrix empty,
- * that is with no memory
- * allocated. This is useful if
- * you want such objects as
- * member variables in other
- * classes. You can make the
- * structure usable by calling
- * the reinit() function.
+ * Initialize the matrix empty, that is with no memory allocated. This is
+ * useful if you want such objects as member variables in other classes. You
+ * can make the structure usable by calling the reinit() function.
*/
CompressedSimpleSparsityPattern ();
/**
- * Copy constructor. This constructor is
- * only allowed to be called if the
- * matrix structure to be copied is
- * empty. This is so in order to prevent
- * involuntary copies of objects for
- * temporaries, which can use large
- * amounts of computing time. However,
- * copy constructors are needed if you
- * want to use the STL data types on
- * classes like this, e.g. to write such
- * statements like <tt>v.push_back
- * (CompressedSparsityPattern());</tt>,
- * with @p v a vector of @p
- * CompressedSparsityPattern objects.
+ * Copy constructor. This constructor is only allowed to be called if the
+ * matrix structure to be copied is empty. This is so in order to prevent
+ * involuntary copies of objects for temporaries, which can use large
+ * amounts of computing time. However, copy constructors are needed if you
+ * want to use the STL data types on classes like this, e.g. to write such
+ * statements like <tt>v.push_back (CompressedSparsityPattern());</tt>, with
+ * @p v a vector of @p CompressedSparsityPattern objects.
*/
CompressedSimpleSparsityPattern (const CompressedSimpleSparsityPattern &);
/**
- * Initialize a rectangular
- * matrix with @p m rows and
- * @p n columns. The @p rowset
- * restricts the storage to
- * elements in rows of this set.
- * Adding elements outside of
- * this set has no effect. The
- * default argument keeps all
- * entries.
+ * Initialize a rectangular matrix with @p m rows and @p n columns. The @p
+ * rowset restricts the storage to elements in rows of this set. Adding
+ * elements outside of this set has no effect. The default argument keeps
+ * all entries.
*/
CompressedSimpleSparsityPattern (const size_type m,
const size_type n,
const IndexSet &rowset = IndexSet());
/**
- * Create a square SparsityPattern using
- * the index set.
+ * Create a square SparsityPattern using the index set.
*/
CompressedSimpleSparsityPattern (const IndexSet &indexset);
/**
- * Initialize a square matrix of
- * dimension @p n.
+ * Initialize a square matrix of dimension @p n.
*/
CompressedSimpleSparsityPattern (const size_type n);
/**
- * Copy operator. For this the
- * same holds as for the copy
- * constructor: it is declared,
- * defined and fine to be called,
- * but the latter only for empty
+ * Copy operator. For this the same holds as for the copy constructor: it is
+ * declared, defined and fine to be called, but the latter only for empty
* objects.
*/
CompressedSimpleSparsityPattern &operator = (const CompressedSimpleSparsityPattern &);
/**
- * Reallocate memory and set up
- * data structures for a new
- * matrix with @p m rows and
- * @p n columns, with at most
- * max_entries_per_row() nonzero
- * entries per row. The @p rowset
- * restricts the storage to
- * elements in rows of this set.
- * Adding elements outside of
- * this set has no effect. The
- * default argument keeps all
- * entries.
+ * Reallocate memory and set up data structures for a new matrix with @p m
+ * rows and @p n columns, with at most max_entries_per_row() nonzero entries
+ * per row. The @p rowset restricts the storage to elements in rows of this
+ * set. Adding elements outside of this set has no effect. The default
+ * argument keeps all entries.
*/
void reinit (const size_type m,
const size_type n,
const IndexSet &rowset = IndexSet());
/**
- * Since this object is kept
- * compressed at all times anway,
- * this function does nothing,
- * but is declared to make the
- * interface of this class as
- * much alike as that of the
- * SparsityPattern class.
+ * Since this object is kept compressed at all times anway, this function
+ * does nothing, but is declared to make the interface of this class as much
+ * alike as that of the SparsityPattern class.
*/
void compress ();
/**
- * Return whether the object is
- * empty. It is empty if no
- * memory is allocated, which is
- * the same as that both
- * dimensions are zero.
+ * Return whether the object is empty. It is empty if no memory is
+ * allocated, which is the same as that both dimensions are zero.
*/
bool empty () const;
/**
- * Return the maximum number of
- * entries per row. Note that
- * this number may change as
- * entries are added.
+ * Return the maximum number of entries per row. Note that this number may
+ * change as entries are added.
*/
size_type max_entries_per_row () const;
/**
- * Add a nonzero entry to the
- * matrix. If the entry already
- * exists, nothing bad happens.
+ * Add a nonzero entry to the matrix. If the entry already exists, nothing
+ * bad happens.
*/
void add (const size_type i,
const size_type j);
/**
- * Add several nonzero entries to the
- * specified row of the matrix. If the
- * entries already exist, nothing bad
- * happens.
+ * Add several nonzero entries to the specified row of the matrix. If the
+ * entries already exist, nothing bad happens.
*/
template <typename ForwardIterator>
void add_entries (const size_type row,
const bool indices_are_unique_and_sorted = false);
/**
- * Check if a value at a certain
- * position may be non-zero.
+ * Check if a value at a certain position may be non-zero.
*/
bool exists (const size_type i,
const size_type j) const;
/**
- * Make the sparsity pattern
- * symmetric by adding the
- * sparsity pattern of the
+ * Make the sparsity pattern symmetric by adding the sparsity pattern of the
* transpose object.
*
- * This function throws an
- * exception if the sparsity
- * pattern does not represent a
- * square matrix.
+ * This function throws an exception if the sparsity pattern does not
+ * represent a square matrix.
*/
void symmetrize ();
/**
- * Print the sparsity of the
- * matrix. The output consists of
- * one line per row of the format
- * <tt>[i,j1,j2,j3,...]</tt>. <i>i</i>
- * is the row number and
- * <i>jn</i> are the allocated
- * columns in this row.
+ * Print the sparsity of the matrix. The output consists of one line per row
+ * of the format <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and
+ * <i>jn</i> are the allocated columns in this row.
*/
void print (std::ostream &out) const;
/**
- * Print the sparsity of the matrix in a
- * format that @p gnuplot understands and
- * which can be used to plot the sparsity
- * pattern in a graphical way. The format
- * consists of pairs <tt>i j</tt> of
- * nonzero elements, each representing
- * one entry of this matrix, one per line
- * of the output file. Indices are
- * counted from zero on, as usual. Since
- * sparsity patterns are printed in the
- * same way as matrices are displayed, we
- * print the negative of the column
- * index, which means that the
- * <tt>(0,0)</tt> element is in the top
- * left rather than in the bottom left
- * corner.
+ * Print the sparsity of the matrix in a format that @p gnuplot understands
+ * and which can be used to plot the sparsity pattern in a graphical
+ * way. The format consists of pairs <tt>i j</tt> of nonzero elements, each
+ * representing one entry of this matrix, one per line of the output
+ * file. Indices are counted from zero on, as usual. Since sparsity patterns
+ * are printed in the same way as matrices are displayed, we print the
+ * negative of the column index, which means that the <tt>(0,0)</tt> element
+ * is in the top left rather than in the bottom left corner.
*
- * Print the sparsity pattern in
- * gnuplot by setting the data style
- * to dots or points and use the
- * @p plot command.
+ * Print the sparsity pattern in gnuplot by setting the data style to dots
+ * or points and use the @p plot command.
*/
void print_gnuplot (std::ostream &out) const;
/**
- * Return number of rows of this
- * matrix, which equals the dimension
- * of the image space.
+ * Return number of rows of this matrix, which equals the dimension of the
+ * image space.
*/
size_type n_rows () const;
/**
- * Return number of columns of this
- * matrix, which equals the dimension
- * of the range space.
+ * Return number of columns of this matrix, which equals the dimension of
+ * the range space.
*/
size_type n_cols () const;
/**
- * Number of entries in a
- * specific row. This function
- * can only be called if the
- * given row is a member of the
- * index set of rows that we want
- * to store.
+ * Number of entries in a specific row. This function can only be called if
+ * the given row is a member of the index set of rows that we want to store.
*/
size_type row_length (const size_type row) const;
/**
- * Access to column number field.
- * Return the column number of
- * the @p indexth entry in @p row.
+ * Access to column number field. Return the column number of the @p
+ * indexth entry in @p row.
*/
size_type column_number (const size_type row,
const size_type index) const;
/**
- * Return an iterator that can loop over
- * all entries in the given
- * row. Dereferencing the iterator yields
- * a column index.
+ * Return an iterator that can loop over all entries in the given
+ * row. Dereferencing the iterator yields a column index.
*/
row_iterator row_begin (const size_type row) const;
*/
row_iterator row_end (const size_type row) const;
/**
- * Compute the bandwidth of the matrix
- * represented by this structure. The
- * bandwidth is the maximum of
- * $|i-j|$ for which the index pair
- * $(i,j)$ represents a nonzero entry
- * of the matrix.
+ * Compute the bandwidth of the matrix represented by this structure. The
+ * bandwidth is the maximum of $|i-j|$ for which the index pair $(i,j)$
+ * represents a nonzero entry of the matrix.
*/
size_type bandwidth () const;
/**
- * Return the number of nonzero elements
- * allocated through this sparsity pattern.
+ * Return the number of nonzero elements allocated through this sparsity
+ * pattern.
*/
size_type n_nonzero_elements () const;
/**
- * Return the IndexSet that sets which
- * rows are active on the current
- * processor. It corresponds to the
- * IndexSet given to this class in the
+ * Return the IndexSet that sets which rows are active on the current
+ * processor. It corresponds to the IndexSet given to this class in the
* constructor or in the reinit function.
*/
const IndexSet &row_index_set () const;
/**
- * return whether this object stores only
- * those entries that have been added
- * explicitly, or if the sparsity pattern
- * contains elements that have been added
- * through other means (implicitly) while
- * building it. For the current class,
- * the result is always true.
+ * return whether this object stores only those entries that have been added
+ * explicitly, or if the sparsity pattern contains elements that have been
+ * added through other means (implicitly) while building it. For the current
+ * class, the result is always true.
*
- * This function mainly serves the
- * purpose of describing the current
- * class in cases where several kinds of
- * sparsity patterns can be passed as
+ * This function mainly serves the purpose of describing the current class
+ * in cases where several kinds of sparsity patterns can be passed as
* template arguments.
*/
static
bool stores_only_added_elements ();
/**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object.
+ * Determine an estimate for the memory consumption (in bytes) of this
+ * object.
*/
size_type memory_consumption () const;
private:
/**
- * Number of rows that this sparsity
- * structure shall represent.
+ * Number of rows that this sparsity structure shall represent.
*/
size_type rows;
/**
- * Number of columns that this sparsity
- * structure shall represent.
+ * Number of columns that this sparsity structure shall represent.
*/
size_type cols;
/**
- * Store some data for each row
- * describing which entries of this row
- * are nonzero. Data is stored sorted in
- * the @p entries std::vector.
- * The vector per row is dynamically
- * growing upon insertion doubling its
- * memory each time.
+ * Store some data for each row describing which entries of this row are
+ * nonzero. Data is stored sorted in the @p entries std::vector. The vector
+ * per row is dynamically growing upon insertion doubling its memory each
+ * time.
*/
struct Line
{
public:
/**
- * Storage for the column indices of
- * this row. This array is always
- * kept sorted.
+ * Storage for the column indices of this row. This array is always kept
+ * sorted.
*/
std::vector<size_type> entries;
Line ();
/**
- * Add the given column number to
- * this line.
+ * Add the given column number to this line.
*/
void add (const size_type col_num);
/**
- * Add the columns specified by the
- * iterator range to this line.
+ * Add the columns specified by the iterator range to this line.
*/
template <typename ForwardIterator>
void add_entries (ForwardIterator begin,
/**
- * Actual data: store for each
- * row the set of nonzero
- * entries.
+ * Actual data: store for each row the set of nonzero entries.
*/
std::vector<Line> lines;
};
void
CompressedSimpleSparsityPattern::Line::add (const size_type j)
{
- // first check the last element (or if line
- // is still empty)
+ // first check the last element (or if line is still empty)
if ( (entries.size()==0) || ( entries.back() < j) )
{
entries.push_back(j);
return;
}
- // do a binary search to find the place
- // where to insert:
+ // do a binary search to find the place where to insert:
std::vector<size_type>::iterator
it = Utilities::lower_bound(entries.begin(),
entries.end(),
j);
- // If this entry is a duplicate, exit
- // immediately
+ // If this entry is a duplicate, exit immediately
if (*it == j)
return;
- // Insert at the right place in the
- // vector. Vector grows automatically to
+ // Insert at the right place in the vector. Vector grows automatically to
// fit elements. Always doubles its size.
entries.insert(it, j);
}
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
// the actual implementation follows in the cm.templates.h file.
distribute_local_to_global (local_matrix, local_vector, local_dof_indices,
global_matrix, global_vector, use_inhomogeneities_for_rhs,
- internal::bool2type<IsBlockMatrix<MatrixType>::value>());
+ dealii::internal::bool2type<IsBlockMatrix<MatrixType>::value>());
}
* Functions related to parallel functionality:
* - The function <code>compress()</code> goes through the data associated
* with ghost indices and communicates it to the owner process, which can
- * then add/set it to the correct position. This can be used e.g. after
+ * then add it to the correct position. This can be used e.g. after
* having run an assembly routine involving ghosts that fill this vector.
+ * Note that the @p insert mode of @p compress() does not set the
+ * elements included in ghost entries but simply discards them, assuming
+ * that the owning processor has set them to the desired value already.
* - The <code>update_ghost_values()</code> function imports the data from
* the owning processor to the ghost indices in order to provide read
* access to the data associated with ghosts.
* writing into ghost elements but only reading from them. This is in
* order to avoid undesired ghost data artifacts when calling compress()
* after modifying some vector entries.
- * The current statues of the ghost entries (read mode or write mode) can
+ * The current status of the ghost entries (read mode or write mode) can
* be queried by the method has_ghost_elements(), which returns
* <code>true</code> exactly when ghost elements have been updated and
* <code>false</code> otherwise, irrespective of the actual number of
* both the element on the ghost site as well as the owning site), this
* operation makes the assumption that all data is set correctly on the
* owning processor. Upon call of compress(VectorOperation::insert), all
- * ghost entries are therefore simply zeroed out (using
- * zero_ghost_values()). In debug mode, a check is performed that makes
- * sure that the data set is actually consistent between processors,
+ * ghost entries are thus simply zeroed out (using
+ * zero_ghost_values()). In debug mode, a check is performed for whether
+ * the data set is actually consistent between processors,
* i.e., whenever a non-zero ghost element is found, it is compared to
* the value on the owning processor and an exception is thrown if these
* elements do not agree.
/**
* Return whether the vector contains only elements with value
- * zero. This function is mainly for internal consistency checks and
- * should seldom be used when not in debug mode since it uses quite some
- * time.
+ * zero. This is a collective operation. This function is expensive, because
+ * potentially all elements have to be checked.
*/
bool all_zero () const;
const bool scientific = true,
const bool across = true) const;
+ /**
+ * @copydoc PETScWrappers::VectorBase::all_zero()
+ *
+ * @note This function overloads the one in the base class
+ * to make this a collective operation.
+ */
+ bool all_zero () const;
+
protected:
/**
* Create a vector of length
const VectorBase &v);
/**
- * Return whether the vector contains
- * only elements with value zero. This
- * function is mainly for internal
- * consistency checks and should
- * seldom be used when not in debug
- * mode since it uses quite some time.
+ * Return whether the vector contains only elements with value
+ * zero. This is a collective operation. This function is expensive, because
+ * potentially all elements have to be checked.
*/
bool all_zero () const;
typedef types::global_dof_index size_type;
/**
- * Use the METIS partitioner to generate
- * a partitioning of the degrees of
- * freedom represented by this sparsity
- * pattern. In effect, we view this
- * sparsity pattern as a graph of
- * connections between various degrees of
- * freedom, where each nonzero entry in
- * the sparsity pattern corresponds to an
- * edge between two nodes in the
- * connection graph. The goal is then to
- * decompose this graph into groups of
- * nodes so that a minimal number of
- * edges are cut by the boundaries
- * between node groups. This partitioning
- * is done by METIS. Note that METIS can
- * only partition symmetric sparsity
- * patterns, and that of course the
- * sparsity pattern has to be square. We
- * do not check for symmetry of the
- * sparsity pattern, since this is an
- * expensive operation, but rather leave
- * this as the responsibility of caller
- * of this function.
+ * Use the METIS partitioner to generate a partitioning of the degrees of
+ * freedom represented by this sparsity pattern. In effect, we view this
+ * sparsity pattern as a graph of connections between various degrees of
+ * freedom, where each nonzero entry in the sparsity pattern corresponds to
+ * an edge between two nodes in the connection graph. The goal is then to
+ * decompose this graph into groups of nodes so that a minimal number of
+ * edges are cut by the boundaries between node groups. This partitioning is
+ * done by METIS. Note that METIS can only partition symmetric sparsity
+ * patterns, and that of course the sparsity pattern has to be square. We do
+ * not check for symmetry of the sparsity pattern, since this is an
+ * expensive operation, but rather leave this as the responsibility of
+ * caller of this function.
*
- * After calling this function, the
- * output array will have values between
- * zero and @p n_partitions-1 for each
- * node (i.e. row or column of the
+ * After calling this function, the output array will have values between
+ * zero and @p n_partitions-1 for each node (i.e. row or column of the
* matrix).
*
- * This function will generate an error
- * if METIS is not installed unless
- * @p n_partitions is one. I.e., you can
- * write a program so that it runs in the
- * single-processor single-partition case
- * without METIS installed, and only
- * requires METIS when multiple
- * partitions are required.
+ * This function will generate an error if METIS is not installed unless @p
+ * n_partitions is one. I.e., you can write a program so that it runs in the
+ * single-processor single-partition case without METIS installed, and only
+ * requires METIS when multiple partitions are required.
*
- * Note that the sparsity pattern itself
- * is not changed by calling this
- * function. However, you will likely use
- * the information generated by calling
- * this function to renumber degrees of
- * freedom, after which you will of
- * course have to regenerate the sparsity
- * pattern.
+ * Note that the sparsity pattern itself is not changed by calling this
+ * function. However, you will likely use the information generated by
+ * calling this function to renumber degrees of freedom, after which you
+ * will of course have to regenerate the sparsity pattern.
*
- * This function will rarely be called
- * separately, since in finite element
- * methods you will want to partition the
- * mesh, not the matrix. This can be done
- * by calling
- * @p GridTools::partition_triangulation.
+ * This function will rarely be called separately, since in finite element
+ * methods you will want to partition the mesh, not the matrix. This can be
+ * done by calling @p GridTools::partition_triangulation.
*/
void partition (const SparsityPattern &sparsity_pattern,
const unsigned int n_partitions,
std::vector<unsigned int> &partition_indices);
/**
- * For a given sparsity pattern, compute a
- * re-enumeration of row/column indices
- * based on the algorithm by Cuthill-McKee.
+ * For a given sparsity pattern, compute a re-enumeration of row/column
+ * indices based on the algorithm by Cuthill-McKee.
*
- * This algorithm is a graph renumbering
- * algorithm in which we attempt to find a
- * new numbering of all nodes of a graph
- * based on their connectivity to other
- * nodes (i.e. the edges that connect
- * nodes). This connectivity is here
- * represented by the sparsity pattern. In
- * many cases within the library, the nodes
- * represent degrees of freedom and edges
- * are nonzero entries in a matrix,
- * i.e. pairs of degrees of freedom that
- * couple through the action of a bilinear
- * form.
+ * This algorithm is a graph renumbering algorithm in which we attempt to
+ * find a new numbering of all nodes of a graph based on their connectivity
+ * to other nodes (i.e. the edges that connect nodes). This connectivity is
+ * here represented by the sparsity pattern. In many cases within the
+ * library, the nodes represent degrees of freedom and edges are nonzero
+ * entries in a matrix, i.e. pairs of degrees of freedom that couple through
+ * the action of a bilinear form.
*
- * The algorithms starts at a node,
- * searches the other nodes for
- * those which are coupled with the one we
- * started with and numbers these in a
- * certain way. It then finds the second
- * level of nodes, namely those that couple
- * with those of the previous level (which
- * were those that coupled with the initial
- * node) and numbers these. And so on. For
- * the details of the algorithm, especially
- * the numbering within each level, we
- * refer the reader to the book of Schwarz
- * (H. R. Schwarz: Methode der finiten
+ * The algorithms starts at a node, searches the other nodes for those which
+ * are coupled with the one we started with and numbers these in a certain
+ * way. It then finds the second level of nodes, namely those that couple
+ * with those of the previous level (which were those that coupled with the
+ * initial node) and numbers these. And so on. For the details of the
+ * algorithm, especially the numbering within each level, we refer the
+ * reader to the book of Schwarz (H. R. Schwarz: Methode der finiten
* Elemente).
*
- * These algorithms have one major
- * drawback: they require a good starting
- * node, i.e. node that will have number
- * zero in the output array. A starting
- * node forming the initial level of nodes
- * can thus be given by the user, e.g. by
- * exploiting knowledge of the actual
- * topology of the domain. It is also
- * possible to give several starting
- * indices, which may be used to simulate a
- * simple upstream numbering (by giving the
- * inflow nodes as starting values) or to
- * make preconditioning faster (by letting
- * the Dirichlet boundary indices be
- * starting points).
+ * These algorithms have one major drawback: they require a good starting
+ * node, i.e. node that will have number zero in the output array. A
+ * starting node forming the initial level of nodes can thus be given by the
+ * user, e.g. by exploiting knowledge of the actual topology of the
+ * domain. It is also possible to give several starting indices, which may
+ * be used to simulate a simple upstream numbering (by giving the inflow
+ * nodes as starting values) or to make preconditioning faster (by letting
+ * the Dirichlet boundary indices be starting points).
*
- * If no starting index is given, one is
- * chosen automatically, namely one with
- * the smallest coordination number (the
- * coordination number is the number of
- * other nodes this node couples
- * with). This node is usually located on
- * the boundary of the domain. There is,
- * however, large ambiguity in this when
- * using the hierarchical meshes used in
- * this library, since in most cases the
- * computational domain is not approximated
- * by tilting and deforming elements and by
- * plugging together variable numbers of
- * elements at vertices, but rather by
- * hierarchical refinement. There is
- * therefore a large number of nodes with
- * equal coordination numbers. The
- * renumbering algorithms will therefore
- * not give optimal results.
+ * If no starting index is given, one is chosen automatically, namely one
+ * with the smallest coordination number (the coordination number is the
+ * number of other nodes this node couples with). This node is usually
+ * located on the boundary of the domain. There is, however, large ambiguity
+ * in this when using the hierarchical meshes used in this library, since in
+ * most cases the computational domain is not approximated by tilting and
+ * deforming elements and by plugging together variable numbers of elements
+ * at vertices, but rather by hierarchical refinement. There is therefore a
+ * large number of nodes with equal coordination numbers. The renumbering
+ * algorithms will therefore not give optimal results.
*
- * If the graph has two or more
- * unconnected components and if no
- * starting indices are given, the
- * algorithm will number each
- * component
- * consecutively. However, this
- * requires the determination of a
- * starting index for each
- * component; as a consequence, the
- * algorithm will produce an
- * exception if starting indices
- * are given, taking the latter as
- * an indication that the caller of
- * the function would like to
- * override the part of the
- * algorithm that chooses starting
- * indices.
+ * If the graph has two or more unconnected components and if no starting
+ * indices are given, the algorithm will number each component
+ * consecutively. However, this requires the determination of a starting
+ * index for each component; as a consequence, the algorithm will produce an
+ * exception if starting indices are given, taking the latter as an
+ * indication that the caller of the function would like to override the
+ * part of the algorithm that chooses starting indices.
*/
void
reorder_Cuthill_McKee (const SparsityPattern &sparsity,
#ifdef DEAL_II_WITH_MPI
/**
- * Communciate rows in a compressed
- * sparsity pattern over MPI.
+ * Communciate rows in a compressed sparsity pattern over MPI.
*
- * @param csp is the sparsity
- * pattern that has been built
- * locally and for which we need to
- * exchange entries with other
- * processors to make sure that
- * each processor knows all the
- * elements of the rows of a matrix
- * it stores and that may
- * eventually be written to. This
- * sparsity pattern will be changed
- * as a result of this function:
- * All entries in rows that belong
- * to a different processor are
- * sent to them and added there.
+ * @param csp is the sparsity pattern that has been built locally and for
+ * which we need to exchange entries with other processors to make sure that
+ * each processor knows all the elements of the rows of a matrix it stores
+ * and that may eventually be written to. This sparsity pattern will be
+ * changed as a result of this function: All entries in rows that belong to
+ * a different processor are sent to them and added there.
*
* @param rows_per_cpu determines ownership of rows.
*
- * @param mpi_comm is the MPI
- * communicator that is shared
- * between the processors that all
- * participate in this operation.
+ * @param mpi_comm is the MPI communicator that is shared between the
+ * processors that all participate in this operation.
*
- * @param myrange indicates the
- * range of elements stored locally
- * and should be the one used in
- * the constructor of the
- * CompressedSimpleSparsityPattern.
- * This should be the locally relevant set.
- * Only
- * rows contained in myrange are
- * checked in csp for transfer.
- * This function needs to be used
- * with
- * PETScWrappers::MPI::SparseMatrix
- * for it to work correctly in a
- * parallel computation.
+ * @param myrange indicates the range of elements stored locally and should
+ * be the one used in the constructor of the
+ * CompressedSimpleSparsityPattern. This should be the locally relevant
+ * set. Only rows contained in myrange are checked in csp for transfer.
+ * This function needs to be used with PETScWrappers::MPI::SparseMatrix for
+ * it to work correctly in a parallel computation.
*/
template <class CSP_t>
void distribute_sparsity_pattern(CSP_t &csp,
*/
template <typename BlockSparsityType>
void reinit (const std::vector<Epetra_Map> &input_maps,
- const BlockSparsityType &block_sparsity_pattern);
+ const BlockSparsityType &block_sparsity_pattern,
+ const bool exchange_data = false);
/**
* Resize the matrix, by using an
template <typename BlockSparsityType>
void reinit (const std::vector<IndexSet> &input_maps,
const BlockSparsityType &block_sparsity_pattern,
- const MPI_Comm &communicator = MPI_COMM_WORLD);
+ const MPI_Comm &communicator = MPI_COMM_WORLD,
+ const bool exchange_data = false);
/**
* Resize the matrix and initialize it
# include "Epetra_SerialComm.h"
# endif
+class Epetra_Export;
+
DEAL_II_NAMESPACE_OPEN
// forward declarations
template <typename number> class SparseMatrix;
class SparsityPattern;
+
namespace TrilinosWrappers
{
// forward declarations
* matrix row at the same time can lead to data races and must be explicitly
* avoided by the user. However, it is possible to access <b>different</b>
* rows of the matrix from several threads simultaneously under the
- * following two conditions:
+ * following three conditions:
* <ul>
* <li> The matrix uses only one MPI process.
+ * <li> The matrix has been initialized with the reinit() method
+ * with a CompressedSimpleSparsityPattern (that includes the set of locally
+ * relevant rows, i.e., the rows that an assembly routine will possibly
+ * write into).
* <li> The matrix has been initialized from a
* TrilinosWrappers::SparsityPattern object that in turn has been
* initialized with the reinit function specifying three index sets, one
* for the rows, one for the columns and for the larger set of @p
- * writeable_rows, and the operation is an addition. Note that all other
- * reinit methods and constructors of TrilinosWrappers::SparsityPattern
- * will result in a matrix that needs to allocate off-processor entries on
- * demand, which breaks thread-safety. Of course, using the respective
- * reinit method for the block Trilinos sparsity pattern and block matrix
- * also results in thread-safety.
+ * writeable_rows, and the operation is an addition.
* </ul>
*
+ * Note that all other reinit methods and constructors of
+ * TrilinosWrappers::SparsityPattern will result in a matrix that needs to
+ * allocate off-processor entries on demand, which breaks thread-safety. Of
+ * course, using the respective reinit method for the block Trilinos
+ * sparsity pattern and block matrix also results in thread-safety.
+ *
* @ingroup TrilinosWrappers
* @ingroup Matrix1
* @author Martin Kronbichler, Wolfgang Bangerth, 2008, 2009
* set, each processor just sets the elements in the sparsity pattern that
* belong to its rows.
*
+ * If the sparsity pattern given to this function is of type
+ * CompressedSimpleSparsity pattern, then a matrix will be created that
+ * allows several threads to write into different rows of the matrix at
+ * the same also with MPI, as opposed to most other reinit() methods.
+ *
* This is a collective operation that needs to be called on all
* processors in order to avoid a dead lock.
*/
*/
std_cxx1x::shared_ptr<Epetra_CrsMatrix> nonlocal_matrix;
+ /**
+ * An export object used to communicate the nonlocal matrix.
+ */
+ std_cxx1x::shared_ptr<Epetra_Export> nonlocal_matrix_exporter;
+
/**
* Trilinos doesn't allow to mix additions to matrix entries and
* overwriting them (to make synchronisation of %parallel computations
for (TrilinosWrappers::types::int_type i=0; i<n_columns; ++i)
std::cout << col_index_ptr[i] << " ";
std::cout << std::endl << std::endl;
- std::cout << "Matrix row has the following indices:" << std::endl;
- int n_indices, *indices;
- trilinos_sparsity_pattern().ExtractMyRowView(row_partitioner().LID(static_cast<TrilinosWrappers::types::int_type>(row)),
- n_indices,
- indices);
+ std::cout << "Matrix row "
+ << (row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) == false ? "(nonlocal part)" : "")
+ << " has the following indices:" << std::endl;
+ std::vector<TrilinosWrappers::types::int_type> indices;
+ const Epetra_CrsGraph* graph =
+ (nonlocal_matrix.get() != 0 &&
+ row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) == false) ?
+ &nonlocal_matrix->Graph() : &matrix->Graph();
+
+ indices.resize(graph->NumGlobalIndices(static_cast<TrilinosWrappers::types::int_type>(row)));
+ int n_indices = 0;
+ graph->ExtractGlobalRowCopy(static_cast<TrilinosWrappers::types::int_type>(row),
+ indices.size(), n_indices, &indices[0]);
+ AssertDimension(static_cast<unsigned int>(n_indices), indices.size());
+
for (TrilinosWrappers::types::int_type i=0; i<n_indices; ++i)
std::cout << indices[i] << " ";
std::cout << std::endl << std::endl;
real_type linfty_norm () const;
/**
- * Return whether the vector contains only elements with value zero. This
- * function is mainly for internal consistency checks and should seldom be
- * used when not in debug mode since it uses quite some time.
+ * Return whether the vector contains only elements with value
+ * zero. This is a collective operation. This function is expensive, because
+ * potentially all elements have to be checked.
*/
bool all_zero () const;
&& (fe.base_element (base_indices.first).face_to_cell_index (line * fe.degree, face)
<= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
&& (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
- < fe.base_element (base_indices.first).face_to_cell_index
- ((line + 1) * fe.degree, face)))
+ <= fe.base_element (base_indices.first).face_to_cell_index
+ ((line + 1) * fe.degree - 1, face)))
|| ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (line * fe.degree <= i)
&& (i < (line + 1) * fe.degree)))
{
DEAL_II_NAMESPACE_OPEN
+
+void
+IndexSet::add_range (const types::global_dof_index begin,
+ const types::global_dof_index end)
+{
+ Assert ((begin < index_space_size)
+ ||
+ ((begin == index_space_size) && (end == index_space_size)),
+ ExcIndexRangeType<types::global_dof_index> (begin, 0, index_space_size));
+ Assert (end <= index_space_size,
+ ExcIndexRangeType<types::global_dof_index> (end, 0, index_space_size+1));
+ Assert (begin <= end,
+ ExcIndexRangeType<types::global_dof_index> (begin, 0, end));
+
+ if (begin != end)
+ {
+ const Range new_range(begin,end);
+
+ // the new index might be larger than the last index present in the
+ // ranges. Then we can skip the binary search
+ if (ranges.size() == 0 || begin > ranges.back().end)
+ ranges.push_back(new_range);
+ else
+ ranges.insert (Utilities::lower_bound (ranges.begin(),
+ ranges.end(),
+ new_range),
+ new_range);
+ is_compressed = false;
+ }
+}
+
+
+
void
IndexSet::do_compress () const
{
- // see if any of the
- // contiguous ranges can be
- // merged. since they are sorted by
- // their first index, determining
+ // see if any of the contiguous ranges can be merged. do not use
+ // std::vector::erase in-place as it is quadratic in the number of
+ // ranges. since the ranges are sorted by their first index, determining
// overlap isn't all that hard
- for (std::vector<Range>::iterator
- i = ranges.begin();
+ std::vector<Range>::iterator store = ranges.begin();
+ for (std::vector<Range>::iterator i = ranges.begin();
i != ranges.end(); )
{
std::vector<Range>::iterator
types::global_dof_index first_index = i->begin;
types::global_dof_index last_index = i->end;
- // see if we can merge any of
- // the following ranges
- bool can_merge = false;
+ // see if we can merge any of the following ranges
while (next != ranges.end() &&
(next->begin <= last_index))
{
last_index = std::max (last_index, next->end);
++next;
- can_merge = true;
}
+ i = next;
- if (can_merge == true)
- {
- // delete the old ranges
- // and insert the new range
- // in place of the previous
- // one
- *i = Range(first_index, last_index);
- i = ranges.erase (i+1, next);
- }
- else
- ++i;
+ // store the new range in the slot we last occupied
+ *store = Range(first_index, last_index);
+ ++store;
+ }
+ // use a compact array with exactly the right amount of storage
+ if (store != ranges.end())
+ {
+ std::vector<Range> new_ranges(ranges.begin(), store);
+ ranges.swap(new_ranges);
}
-
- // now compute indices within set and the
- // range with most elements
+ // now compute indices within set and the range with most elements
types::global_dof_index next_index = 0, largest_range_size = 0;
- for (std::vector<Range>::iterator
- i = ranges.begin();
- i != ranges.end();
+ for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end();
++i)
{
Assert(i->begin < i->end, ExcInternalError());
}
is_compressed = true;
- // check that next_index is
- // correct. needs to be after the
- // previous statement because we
- // otherwise will get into an
- // endless loop
+ // check that next_index is correct. needs to be after the previous
+ // statement because we otherwise will get into an endless loop
Assert (next_index == n_elements(), ExcInternalError());
}
&&
(r2 != is.ranges.end()))
{
- // if r1 and r2 do not overlap
- // at all, then move the
- // pointer that sits to the
- // left of the other up by one
+ // if r1 and r2 do not overlap at all, then move the pointer that sits
+ // to the left of the other up by one
if (r1->end <= r2->begin)
++r1;
else if (r2->end <= r1->begin)
++r2;
else
{
- // the ranges must overlap
- // somehow
+ // the ranges must overlap somehow
Assert (((r1->begin <= r2->begin) &&
(r1->end > r2->begin))
||
(r2->end > r1->begin)),
ExcInternalError());
- // add the overlapping
- // range to the result
+ // add the overlapping range to the result
result.add_range (std::max (r1->begin,
r2->begin),
std::min (r1->end,
r2->end));
- // now move that iterator
- // that ends earlier one
- // up. note that it has to
- // be this one because a
- // subsequent range may
- // still have a chance of
- // overlapping with the
- // range that ends later
+ // now move that iterator that ends earlier one up. note that it has
+ // to be this one because a subsequent range may still have a chance
+ // of overlapping with the range that ends later
if (r1->end <= r2->end)
++r1;
else
-unsigned int
-IndexSet::n_intervals () const
-{
- compress ();
- return ranges.size();
-}
-
-
-
IndexSet
IndexSet::get_view (const types::global_dof_index begin,
const types::global_dof_index end) const
+void
+IndexSet::subtract_set (const IndexSet &other)
+{
+ compress();
+ other.compress();
+ is_compressed = false;
+
+
+ // we save new ranges to be added to our IndexSet in an temporary list and
+ // add all of them in one go at the end. This is necessary because a growing
+ // ranges vector invalidates iterators.
+ std::list<Range> temp_list;
+
+ std::vector<Range>::iterator own_it = ranges.begin();
+ std::vector<Range>::iterator other_it = other.ranges.begin();
+
+ while (own_it != ranges.end() && other_it != other.ranges.end())
+ {
+ //advance own iterator until we get an overlap
+ if (own_it->end <= other_it->begin)
+ {
+ ++own_it;
+ continue;
+ }
+ //we are done with other_it, so advance
+ if (own_it->begin >= other_it->end)
+ {
+ ++other_it;
+ continue;
+ }
+
+ //Now own_it and other_it overlap. First save the part of own_it that
+ //is before other_it (if not empty).
+ if (own_it->begin < other_it->begin)
+ {
+ Range r(own_it->begin, other_it->begin);
+ r.nth_index_in_set = 0; //fix warning of unused variable
+ temp_list.push_back(r);
+ }
+ // change own_it to the sub range behind other_it. Do not delete own_it
+ // in any case. As removal would invalidate iterators, we just shrink
+ // the range to an empty one.
+ own_it->begin = other_it->end;
+ if (own_it->begin > own_it->end)
+ {
+ own_it->begin = own_it->end;
+ ++own_it;
+ }
+
+ // continue without advancing iterators, the right one will be advanced
+ // next.
+ }
+
+ // Now delete all empty ranges we might
+ // have created.
+ for (std::vector<Range>::iterator it = ranges.begin();
+ it != ranges.end(); )
+ {
+ if (it->begin >= it->end)
+ it = ranges.erase(it);
+ else
+ ++it;
+ }
+
+ // done, now add the temporary ranges
+ for (std::list<Range>::iterator it = temp_list.begin();
+ it != temp_list.end();
+ ++it)
+ add_range(it->begin, it->end);
+
+ compress();
+}
+
+
+
+void
+IndexSet::add_indices(const IndexSet &other,
+ const unsigned int offset)
+{
+ if ((this == &other) && (offset == 0))
+ return;
+
+ compress();
+ other.compress();
+
+ std::vector<Range>::const_iterator r1 = ranges.begin(),
+ r2 = other.ranges.begin();
+
+ std::vector<Range> new_ranges;
+ // just get the start and end of the ranges right in this method, everything
+ // else will be done in compress()
+ while (r1 != ranges.end() || r2 != other.ranges.end())
+ {
+ // the two ranges do not overlap or we are at the end of one of the
+ // ranges
+ if (r2 == other.ranges.end() ||
+ (r1 != ranges.end() && r1->end < (r2->begin+offset)))
+ {
+ new_ranges.push_back(*r1);
+ ++r1;
+ }
+ else if (r1 == ranges.end() || (r2->end+offset) < r1->begin)
+ {
+ new_ranges.push_back(Range(r2->begin+offset,r2->end+offset));
+ ++r2;
+ }
+ else
+ {
+ // ok, we do overlap, so just take the combination of the current
+ // range (do not bother to merge with subsequent ranges)
+ Range next(std::min(r1->begin, r2->begin+offset),
+ std::max(r1->end, r2->end+offset));
+ new_ranges.push_back(next);
+ ++r1;
+ ++r2;
+ }
+ }
+ ranges.swap(new_ranges);
+
+ is_compressed = false;
+ compress();
+}
+
+
+
void
IndexSet::write(std::ostream &out) const
{
-void
-IndexSet::subtract_set (const IndexSet &other)
-{
- compress();
- other.compress();
- is_compressed = false;
-
-
- // we save new ranges to be added to our
- // IndexSet in an temporary list and add
- // all of them in one go at the end. This
- // is necessary because a growing ranges
- // vector invalidates iterators.
- std::list<Range> temp_list;
-
- std::vector<Range>::iterator own_it = ranges.begin();
- std::vector<Range>::iterator other_it = other.ranges.begin();
-
- while (own_it != ranges.end() && other_it != other.ranges.end())
- {
- //advance own iterator until we get an
- //overlap
- if (own_it->end <= other_it->begin)
- {
- ++own_it;
- continue;
- }
- //we are done with other_it, so advance
- if (own_it->begin >= other_it->end)
- {
- ++other_it;
- continue;
- }
-
- //Now own_it and other_it overlap.
- //First save the part of own_it that is
- //before other_it (if not empty).
- if (own_it->begin < other_it->begin)
- {
- Range r(own_it->begin, other_it->begin);
- r.nth_index_in_set = 0; //fix warning of unused variable
- temp_list.push_back(r);
- }
- // change own_it to the sub range
- // behind other_it. Do not delete
- // own_it in any case. As removal would
- // invalidate iterators, we just shrink
- // the range to an empty one.
- own_it->begin = other_it->end;
- if (own_it->begin > own_it->end)
- {
- own_it->begin = own_it->end;
- ++own_it;
- }
-
- // continue without advancing
- // iterators, the right one will be
- // advanced next.
- }
-
- // Now delete all empty ranges we might
- // have created.
- for (std::vector<Range>::iterator it = ranges.begin();
- it != ranges.end(); )
- {
- if (it->begin >= it->end)
- it = ranges.erase(it);
- else
- ++it;
- }
-
- // done, now add the temporary ranges
- for (std::list<Range>::iterator it = temp_list.begin();
- it != temp_list.end();
- ++it)
- add_range(it->begin, it->end);
-
- compress();
-}
-
-
void IndexSet::fill_index_vector(std::vector<types::global_dof_index> &indices) const
{
compress();
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
void
DoFCellAccessor<DH,lda>::
get_interpolated_dof_values (const InputVector &values,
- Vector<number> &interpolated_values) const
+ Vector<number> &interpolated_values,
+ const unsigned int fe_index) const
{
- const FiniteElement<dim,spacedim> &fe = this->get_fe();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
- Assert (this->dof_handler != 0,
- typename BaseClass::ExcInvalidObject());
- Assert (&fe != 0,
- typename BaseClass::ExcInvalidObject());
- Assert (interpolated_values.size() == dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- Assert (values.size() == this->dof_handler->n_dofs(),
- typename BaseClass::ExcVectorDoesNotMatch());
-
if (!this->has_children())
// if this cell has no children: simply return the exact values on this
- // cell
- this->get_dof_values (values, interpolated_values);
+ // cell unless the finite element we need to interpolate to is different than
+ // the one we have on the current cell
+ {
+ if ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
+ (this->dof_handler)
+ != 0)
+ ||
+ (fe_index == this->active_fe_index())
+ ||
+ (fe_index == DH::default_fe_index))
+ this->get_dof_values (values, interpolated_values);
+ else
+ {
+ // well, here we need to first get the values from the current
+ // cell and then interpolate it to the element requested. this
+ // can clearly only happen for hp::DoFHandler objects
+ Vector<number> tmp (this->get_fe().dofs_per_cell);
+ this->get_dof_values (values, tmp);
+
+ FullMatrix<double> interpolation (this->dof_handler->get_fe()[fe_index].dofs_per_cell,
+ this->get_fe().dofs_per_cell);
+ this->dof_handler->get_fe()[fe_index].get_interpolation_matrix (this->get_fe(),
+ interpolation);
+ interpolation.vmult (interpolated_values, tmp);
+ }
+ }
else
- // otherwise clobber them from the children
+ // otherwise obtain them from the children
{
+ // we are on a non-active cell. these do not have any finite
+ // element associated with them in the hp context (in the non-hp
+ // context, we can simply assume that the FE space to which we
+ // want to interpolate is the same as for all elements in the
+ // mesh). consequently, we cannot interpolate from children's FE
+ // space to this cell's (unknown) FE space unless an explicit
+ // fe_index is given
+ Assert ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
+ (this->dof_handler)
+ != 0)
+ ||
+ (fe_index != DH::default_fe_index),
+ ExcMessage ("You cannot call this function on non-active cells "
+ "of hp::DoFHandler objects unless you provide an explicit "
+ "finite element index because they do not have naturally "
+ "associated finite element spaces associated: degrees "
+ "of freedom are only distributed on active cells for which "
+ "the active_fe_index has been set."));
+
+ const FiniteElement<dim,spacedim> &fe = this->get_dof_handler().get_fe()[fe_index];
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ Assert (this->dof_handler != 0,
+ typename BaseClass::ExcInvalidObject());
+ Assert (&fe != 0,
+ typename BaseClass::ExcInvalidObject());
+ Assert (interpolated_values.size() == dofs_per_cell,
+ typename BaseClass::ExcVectorDoesNotMatch());
+ Assert (values.size() == this->dof_handler->n_dofs(),
+ typename BaseClass::ExcVectorDoesNotMatch());
+
+
Vector<number> tmp1(dofs_per_cell);
Vector<number> tmp2(dofs_per_cell);
for (unsigned int child=0; child<this->n_children(); ++child)
{
// get the values from the present child, if necessary by
- // interpolation itself
+ // interpolation itself either from its own children or
+ // by interpolating from the finite element on an active
+ // child to the finite element space requested here
this->child(child)->get_interpolated_dof_values (values,
- tmp1);
+ tmp1,
+ fe_index);
// interpolate these to the mother cell
fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1);
template
void
DoFCellAccessor<DoFHandler<deal_II_dimension>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#if deal_II_dimension != 3
template
void
DoFCellAccessor<DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<DoFHandler<1,3>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<hp::DoFHandler<deal_II_dimension>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#if deal_II_dimension != 3
template
void
DoFCellAccessor<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<hp::DoFHandler<1,3>, lda>::get_interpolated_dof_values
- (const VEC&, Vector<SCALAR>&) const;
+ (const VEC&, Vector<SCALAR>&,
+ const unsigned int fe_index) const;
#endif
}
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
void
DoFCellAccessor<DH,lda>::
set_dof_values_by_interpolation (const Vector<number> &local_values,
- OutputVector &values) const
+ OutputVector &values,
+ const unsigned int fe_index) const
{
- const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell;
-
- Assert (this->dof_handler != 0,
- typename BaseClass::ExcInvalidObject());
- Assert (&this->get_fe() != 0,
- typename BaseClass::ExcInvalidObject());
- Assert (local_values.size() == dofs_per_cell,
- typename BaseClass::ExcVectorDoesNotMatch());
- Assert (values.size() == this->dof_handler->n_dofs(),
- typename BaseClass::ExcVectorDoesNotMatch());
-
if (!this->has_children())
// if this cell has no children: simply set the values on this cell
this->set_dof_values (local_values, values);
else
// otherwise distribute them to the children
{
+ // we are on a non-active cell. these do not have any finite
+ // element associated with them in the hp context (in the non-hp
+ // context, we can simply assume that the FE space to from which we
+ // want to interpolate is the same as for all elements in the
+ // mesh). consequently, we cannot interpolate to children's FE
+ // space from this cell's (unknown) FE space unless an explicit
+ // fe_index is given
+ Assert ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
+ (this->dof_handler)
+ != 0)
+ ||
+ (fe_index != DH::default_fe_index),
+ ExcMessage ("You cannot call this function on non-active cells "
+ "of hp::DoFHandler objects unless you provide an explicit "
+ "finite element index because they do not have naturally "
+ "associated finite element spaces associated: degrees "
+ "of freedom are only distributed on active cells for which "
+ "the active_fe_index has been set."));
+
+ const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell;
+
+ Assert (this->dof_handler != 0,
+ typename BaseClass::ExcInvalidObject());
+ Assert (&this->get_fe() != 0,
+ typename BaseClass::ExcInvalidObject());
+ Assert (local_values.size() == dofs_per_cell,
+ typename BaseClass::ExcVectorDoesNotMatch());
+ Assert (values.size() == this->dof_handler->n_dofs(),
+ typename BaseClass::ExcVectorDoesNotMatch());
+
Vector<number> tmp(dofs_per_cell);
for (unsigned int child=0; child<this->n_children(); ++child)
template
void
DoFCellAccessor<DoFHandler<deal_II_dimension>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#if deal_II_dimension != 3
template
void
DoFCellAccessor<DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<DoFHandler<1,3>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<hp::DoFHandler<deal_II_dimension>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#if deal_II_dimension != 3
template
void
DoFCellAccessor<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#endif
template
void
DoFCellAccessor<hp::DoFHandler<1,3>, lda>::set_dof_values_by_interpolation
- (const Vector<SCALAR>&, VEC&) const;
+ (const Vector<SCALAR>&, VEC&,
+ const unsigned int fe_index) const;
#endif
}
unsigned int sent=needs_to_get_cells.size();
unsigned int recv=senders.size();
- MPI_Reduce(&sent, &sum_send, 1, MPI_UNSIGNED, MPI_SUM, 56345, tr->get_communicator());
- MPI_Reduce(&recv, &sum_recv, 1, MPI_UNSIGNED, MPI_SUM, 56346, tr->get_communicator());
+ MPI_Allreduce(&sent, &sum_send, 1, MPI_UNSIGNED, MPI_SUM, tr->get_communicator());
+ MPI_Allreduce(&recv, &sum_recv, 1, MPI_UNSIGNED, MPI_SUM, tr->get_communicator());
Assert(sum_send==sum_recv, ExcInternalError());
}
#endif
n_selected_dofs += std::count (dofs_by_component.begin(),
dofs_by_component.end(), i);
+ // Find local numbering within the selected components
+ const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
+ std::vector<unsigned int> component_numbering(locally_owned_dofs.n_elements(),
+ numbers::invalid_unsigned_int);
+ for (unsigned int i=0, count=0; i<locally_owned_dofs.n_elements(); ++i)
+ if (component_mask[dofs_by_component[i]])
+ component_numbering[i] = count++;
+
// First count the number of dofs in the current component.
constant_modes.resize (n_components_selected, std::vector<bool>(n_selected_dofs,
false));
- std::vector<unsigned int> component_list (n_components, 0);
- for (unsigned int d=0; d<n_components; ++d)
- component_list[d] = component_mask[d];
- unsigned int counter = 0;
- for (unsigned int i=0; i<dof_handler.n_locally_owned_dofs(); ++i)
- if (component_mask[dofs_by_component[i]])
+ // Loop over all owned cells and ask the element for the constant modes
+ const dealii::hp::FECollection<DH::dimension,DH::space_dimension>
+ fe_collection (dof_handler.get_fe());
+ std::vector<Table<2,bool> > element_constant_modes;
+ for (unsigned int f=0; f<fe_collection.size(); ++f)
+ element_constant_modes.push_back(fe_collection[f].get_constant_modes());
+
+ typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ std::vector<types::global_dof_index> dof_indices;
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
{
- constant_modes[localized_component[dofs_by_component[i]]][counter] = true;
- ++counter;
+ dof_indices.resize(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(dof_indices);
+
+ for (unsigned int i=0; i<dof_indices.size(); ++i)
+ if (locally_owned_dofs.is_element(dof_indices[i]))
+ {
+ const unsigned int loc_index =
+ locally_owned_dofs.index_within_set(dof_indices[i]);
+ const unsigned int comp = dofs_by_component[loc_index];
+ if (component_mask[comp])
+ constant_modes[localized_component[comp]][component_numbering[loc_index]] =
+ element_constant_modes[cell->active_fe_index()](comp,i);
+ }
}
}
constraints.add_entries_local_to_global
(dofs_on_other_cell, dofs_on_this_cell,
sparsity, keep_constrained_dofs);
+ // only need to add this when the neighbor is not
+ // owned by the current processor, otherwise we add
+ // the entries for the neighbor there
+ if (sub_neighbor->subdomain_id() != cell->subdomain_id())
+ constraints.add_entries_local_to_global
+ (dofs_on_other_cell, sparsity, keep_constrained_dofs);
}
}
else
{
// Refinement edges are taken care of by coarser
// cells
-
- // TODO: in the distributed case, we miss out the
- // constraints when the neighbor cell is coarser, but
- // only the current cell is owned locally!
- if (cell->neighbor_is_coarser(face))
+ if (cell->neighbor_is_coarser(face) &&
+ neighbor->subdomain_id() == cell->subdomain_id())
continue;
const unsigned int n_dofs_on_neighbor
// around
if (!cell->neighbor(face)->active()
||
- (cell->neighbor(face)->subdomain_id() !=
- cell->subdomain_id()))
- constraints.add_entries_local_to_global
- (dofs_on_other_cell, dofs_on_this_cell,
- sparsity, keep_constrained_dofs);
+ (neighbor->subdomain_id() != cell->subdomain_id()))
+ {
+ constraints.add_entries_local_to_global
+ (dofs_on_other_cell, dofs_on_this_cell,
+ sparsity, keep_constrained_dofs);
+ if (neighbor->subdomain_id() != cell->subdomain_id())
+ constraints.add_entries_local_to_global
+ (dofs_on_other_cell, sparsity, keep_constrained_dofs);
+ }
}
}
}
}
+
template <int dim, int spacedim>
bool
FiniteElement<dim,spacedim>::has_support_on_face (
}
+
+template <int dim, int spacedim>
+Table<2,bool>
+FiniteElement<dim,spacedim>::get_constant_modes () const
+{
+ Assert (false, ExcNotImplemented());
+ return Table<2,bool>(this->n_components(), this->dofs_per_cell);
+}
+
+
+
template <int dim, int spacedim>
void
FiniteElement<dim,spacedim>::interpolate(
+template <int dim, int spacedim>
+Table<2,bool>
+FE_DGP<dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ constant_modes(0,0) = true;
+ return constant_modes;
+}
+
+
+
template <int dim, int spacedim>
std::size_t
FE_DGP<dim,spacedim>::memory_consumption () const
+template <int dim, int spacedim>
+Table<2,bool>
+FE_DGQ<dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ constant_modes(0,i) = true;
+ return constant_modes;
+}
+
+
+
+
template <int dim, int spacedim>
std::size_t
FE_DGQ<dim, spacedim>::memory_consumption () const
+template <int dim, int spacedim>
+Table<2,bool>
+FE_FaceQ<dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ constant_modes(0,i) = true;
+ return constant_modes;
+}
+
+
+
// --------------------------------------- FE_FaceP --------------------------
template <int dim, int spacedim>
}
+
+template <int dim, int spacedim>
+Table<2,bool>
+FE_FaceP<dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ constant_modes(0, face*this->dofs_per_face) = true;
+ return constant_modes;
+}
+
+
+
// explicit instantiations
#include "fe_face.inst"
}
+
+template <int dim>
+Table<2,bool>
+FE_Nedelec<dim>::get_constant_modes() const
+{
+ Table<2,bool> constant_modes(dim, this->dofs_per_cell);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ constant_modes(d,i) = true;
+ return constant_modes;
+}
+
+
template <int dim>
std::size_t
FE_Nedelec<dim>::memory_consumption () const
+template <typename POLY, int dim, int spacedim>
+Table<2,bool>
+FE_Q_Base<POLY,dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ // leave out last component
+ for (unsigned int i=0; i<Utilities::fixed_power<dim>(this->degree+1); ++i)
+ constant_modes(0, i) = true;
+ return constant_modes;
+}
+
+
// explicit instantiations
#include "fe_q_base.inst"
+template <int dim>
+Table<2,bool>
+FE_Q_Hierarchical<dim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(1, this->dofs_per_cell);
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ constant_modes(0,i) = true;
+ for (unsigned int i=GeometryInfo<dim>::vertices_per_cell; i<this->dofs_per_cell; ++i)
+ constant_modes(0,i) = false;
+ return constant_modes;
+}
+
+
+
template <int dim>
std::size_t
FE_Q_Hierarchical<dim>::memory_consumption () const
+template <int dim>
+Table<2,bool>
+FE_RaviartThomas<dim>::get_constant_modes() const
+{
+ Table<2,bool> constant_modes(dim, this->dofs_per_cell);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ constant_modes(d,i) = true;
+ return constant_modes;
+}
+
+
+
//---------------------------------------------------------------------------
// Data field initialization
//---------------------------------------------------------------------------
}
+
template <int dim>
std::size_t
FE_RaviartThomas<dim>::memory_consumption () const
+template <int dim, int spacedim>
+Table<2,bool>
+FESystem<dim,spacedim>::get_constant_modes () const
+{
+ Table<2,bool> constant_modes(this->n_components(), this->dofs_per_cell);
+ unsigned int comp=0;
+ for (unsigned int i=0; i<base_elements.size(); ++i)
+ {
+ Table<2,bool> base_table = base_elements[i].first->get_constant_modes();
+ const unsigned int n_base_components = base_elements[i].first->n_components();
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ {
+ std::pair<std::pair<unsigned int,unsigned int>, unsigned int> ind
+ = this->system_to_base_index(k);
+ if (ind.first.first == i)
+ for (unsigned int c=0; c<base_table.n_rows(); ++c)
+ constant_modes(comp+ind.first.second*n_base_components+c,k)
+ = base_table(c,ind.second);
+ }
+ comp += n_base_components * base_elements[i].second;
+ }
+ AssertDimension(comp, this->n_components());
+ return constant_modes;
+}
+
+
+
template <int dim, int spacedim>
std::size_t
FESystem<dim,spacedim>::memory_consumption () const
+ bool
+ Vector::all_zero() const
+ {
+ unsigned int has_nonzero = VectorBase::all_zero()?0:1;
+#ifdef DEAL_II_WITH_MPI
+ // in parallel, check that the vector
+ // is zero on _all_ processors.
+ unsigned int num_nonzero = Utilities::MPI::sum(has_nonzero, communicator);
+ return num_nonzero == 0;
+#else
+ return has_nonzero == 0;
+#endif
+ }
+
+
void
Vector::print (std::ostream &out,
const unsigned int precision,
AssertThrow (ierr == 0, ExcPETScError(ierr));
const PetscScalar *ptr = start_ptr,
- *eptr = start_ptr + size();
+ *eptr = start_ptr + local_size();
bool flag = true;
while (ptr != eptr)
{
AssertThrow (ierr == 0, ExcPETScError(ierr));
const PetscScalar *ptr = start_ptr,
- *eptr = start_ptr + size();
+ *eptr = start_ptr + local_size();
bool flag = true;
while (ptr != eptr)
{
Assert (status.MPI_TAG==124, ExcInternalError());
MPI_Get_count(&status, MPI_BYTE, &len);
- Assert( len%sizeof(unsigned int)==0, ExcInternalError());
+ Assert( len%sizeof(size_type)==0, ExcInternalError());
recv_buf.resize(len/sizeof(size_type));
}
}
- // complete all sends, so that we can
- // safely destroy the buffers.
+ // complete all sends, so that we can safely destroy the buffers.
MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
}
}
}
- // complete all sends, so that we can
- // safely destroy the buffers.
+ // complete all sends, so that we can safely destroy the buffers.
MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
}
void
BlockSparseMatrix::
reinit (const std::vector<Epetra_Map> ¶llel_partitioning,
- const BlockSparsityType &block_sparsity_pattern)
+ const BlockSparsityType &block_sparsity_pattern,
+ const bool exchange_data)
{
Assert (parallel_partitioning.size() == block_sparsity_pattern.n_block_rows(),
ExcDimensionMismatch (parallel_partitioning.size(),
{
this->sub_objects[r][c]->reinit (parallel_partitioning[r],
parallel_partitioning[c],
- block_sparsity_pattern.block(r,c));
+ block_sparsity_pattern.block(r,c),
+ exchange_data);
}
}
BlockSparseMatrix::
reinit (const std::vector<IndexSet> ¶llel_partitioning,
const BlockSparsityType &block_sparsity_pattern,
- const MPI_Comm &communicator)
+ const MPI_Comm &communicator,
+ const bool exchange_data)
{
std::vector<Epetra_Map> epetra_maps;
for (size_type i=0; i<block_sparsity_pattern.n_block_rows(); ++i)
epetra_maps.push_back
(parallel_partitioning[i].make_trilinos_map(communicator, false));
- reinit (epetra_maps, block_sparsity_pattern);
+ reinit (epetra_maps, block_sparsity_pattern, exchange_data);
}
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
- const dealii::BlockSparsityPattern &);
+ const dealii::BlockSparsityPattern &,
+ const bool);
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
- const dealii::BlockCompressedSparsityPattern &);
+ const dealii::BlockCompressedSparsityPattern &,
+ const bool);
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
- const dealii::BlockCompressedSetSparsityPattern &);
+ const dealii::BlockCompressedSetSparsityPattern &,
+ const bool);
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
- const dealii::BlockCompressedSimpleSparsityPattern &);
+ const dealii::BlockCompressedSimpleSparsityPattern &,
+ const bool);
template void
BlockSparseMatrix::reinit (const std::vector<IndexSet> &,
const dealii::BlockCompressedSimpleSparsityPattern &,
- const MPI_Comm &);
+ const MPI_Comm &,
+ const bool);
}
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2008 - 2013 by the deal.II authors
+// Copyright (C) 2008 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
# include <deal.II/lac/compressed_sparsity_pattern.h>
# include <deal.II/lac/compressed_set_sparsity_pattern.h>
# include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+# include <deal.II/lac/sparsity_tools.h>
# include <Epetra_Export.h>
# include <ml_epetra_utils.h>
SparseMatrix::copy_from (const SparseMatrix &m)
{
nonlocal_matrix.reset();
+ nonlocal_matrix_exporter.reset();
// check whether we need to update the
// partitioner or can just copy the data:
// release memory before reallocation
matrix.reset();
nonlocal_matrix.reset();
+ nonlocal_matrix_exporter.reset();
// if we want to exchange data, build a usual Trilinos sparsity pattern
// and let that handle the exchange. otherwise, manually create a
+ // specialization for CompressedSimpleSparsityPattern which can provide us
+ // with more information about the non-locally owned rows
+ template <>
+ void
+ SparseMatrix::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const CompressedSimpleSparsityPattern &sparsity_pattern,
+ const bool exchange_data)
+ {
+ matrix.reset();
+ nonlocal_matrix.reset();
+ nonlocal_matrix_exporter.reset();
+
+ AssertDimension (sparsity_pattern.n_rows(),
+ static_cast<size_type>(n_global_elements(input_row_map)));
+
+ // exchange data if requested with deal.II's own function
+#ifdef DEAL_II_WITH_MPI
+ const Epetra_MpiComm* communicator =
+ dynamic_cast<const Epetra_MpiComm*>(&input_row_map.Comm());
+ if (exchange_data && communicator != 0)
+ {
+ std::vector<size_type>
+ owned_per_proc(communicator->NumProc(), -1);
+ size_type my_elements = input_row_map.NumMyElements();
+ MPI_Allgather(&my_elements, 1,
+ Utilities::MPI::internal::mpi_type_id(&my_elements),
+ &owned_per_proc[0], 1,
+ Utilities::MPI::internal::mpi_type_id(&my_elements),
+ communicator->Comm());
+
+ AssertDimension(std::accumulate(owned_per_proc.begin(),
+ owned_per_proc.end(), size_type()), sparsity_pattern.n_rows());
+
+ SparsityTools::distribute_sparsity_pattern
+ (const_cast<CompressedSimpleSparsityPattern&>(sparsity_pattern),
+ owned_per_proc, communicator->Comm(), sparsity_pattern.row_index_set());
+ }
+#endif
+
+ if (input_row_map.Comm().MyPID() == 0)
+ {
+ AssertDimension (sparsity_pattern.n_rows(),
+ static_cast<size_type>(n_global_elements(input_row_map)));
+ AssertDimension (sparsity_pattern.n_cols(),
+ static_cast<size_type>(n_global_elements(input_col_map)));
+ }
+
+ column_space_map.reset (new Epetra_Map (input_col_map));
+
+ IndexSet relevant_rows (sparsity_pattern.row_index_set());
+ // serial case
+ if (relevant_rows.size() == 0)
+ {
+ relevant_rows.set_size(n_global_elements(input_row_map));
+ relevant_rows.add_range(0, n_global_elements(input_row_map));
+ }
+ relevant_rows.compress();
+ Assert(relevant_rows.n_elements() >= static_cast<unsigned int>(input_row_map.NumMyElements()),
+ ExcMessage("Locally relevant rows of sparsity pattern must contain "
+ "all locally owned rows"));
+
+ // check whether the relevant rows correspond to exactly the same map as
+ // the owned rows. In that case, do not create the nonlocal graph and fill
+ // the columns by demand
+ bool have_ghost_rows = false;
+ {
+ std::vector<dealii::types::global_dof_index> indices;
+ relevant_rows.fill_index_vector(indices);
+ Epetra_Map relevant_map (TrilinosWrappers::types::int_type(-1),
+ TrilinosWrappers::types::int_type(relevant_rows.n_elements()),
+ (indices.empty() ? 0 :
+ reinterpret_cast<TrilinosWrappers::types::int_type *>(&indices[0])),
+ 0, input_row_map.Comm());
+ if (relevant_map.SameAs(input_row_map))
+ have_ghost_rows = false;
+ else
+ have_ghost_rows = true;
+ }
+
+ const unsigned int n_rows = relevant_rows.n_elements();
+ std::vector<TrilinosWrappers::types::int_type> ghost_rows;
+ std::vector<int> n_entries_per_row(input_row_map.NumMyElements());
+ std::vector<int> n_entries_per_ghost_row;
+ for (unsigned int i=0, own=0; i<n_rows; ++i)
+ {
+ const TrilinosWrappers::types::int_type global_row =
+ relevant_rows.nth_index_in_set(i);
+ if (input_row_map.MyGID(global_row))
+ n_entries_per_row[own++] = sparsity_pattern.row_length(global_row);
+ else if (sparsity_pattern.row_length(global_row) > 0)
+ {
+ ghost_rows.push_back(global_row);
+ n_entries_per_ghost_row.push_back(sparsity_pattern.row_length(global_row));
+ }
+ }
+
+ // make sure all processors create an off-processor matrix with at least
+ // one entry
+ if (have_ghost_rows == true && ghost_rows.empty() == true)
+ {
+ ghost_rows.push_back(0);
+ n_entries_per_ghost_row.push_back(1);
+ }
+
+ Epetra_Map off_processor_map(-1, ghost_rows.size(), &ghost_rows[0],
+ 0, input_row_map.Comm());
+
+ std_cxx1x::shared_ptr<Epetra_CrsGraph> graph, nonlocal_graph;
+ if (input_row_map.Comm().NumProc() > 1)
+ {
+ graph.reset (new Epetra_CrsGraph (Copy, input_row_map,
+ &n_entries_per_row[0], true));
+ if (have_ghost_rows == true)
+ nonlocal_graph.reset (new Epetra_CrsGraph (Copy, off_processor_map,
+ &n_entries_per_ghost_row[0],
+ true));
+ }
+ else
+ graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
+ &n_entries_per_row[0], true));
+
+ // now insert the indices, select between the right matrix
+ std::vector<TrilinosWrappers::types::int_type> row_indices;
+
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ const TrilinosWrappers::types::int_type global_row =
+ relevant_rows.nth_index_in_set(i);
+ const int row_length = sparsity_pattern.row_length(global_row);
+ if (row_length == 0)
+ continue;
+
+ row_indices.resize (row_length, -1);
+ internal::copy_row(sparsity_pattern, global_row, row_indices);
+
+ if (input_row_map.MyGID(global_row))
+ graph->InsertGlobalIndices (global_row, row_length, &row_indices[0]);
+ else
+ {
+ Assert(nonlocal_graph.get() != 0, ExcInternalError());
+ nonlocal_graph->InsertGlobalIndices (global_row, row_length,
+ &row_indices[0]);
+ }
+ }
+
+ graph->FillComplete(input_col_map, input_row_map);
+ graph->OptimizeStorage();
+
+ AssertDimension (sparsity_pattern.n_cols(),static_cast<size_type>(
+ n_global_cols(*graph)));
+
+ matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
+
+ // and now the same operations for the nonlocal graph and matrix
+ if (nonlocal_graph.get() != 0)
+ {
+ if (nonlocal_graph->IndicesAreGlobal() == false &&
+ nonlocal_graph->RowMap().NumMyElements() > 0)
+ {
+ // insert dummy element
+ TrilinosWrappers::types::int_type row =
+ nonlocal_graph->RowMap().MyGID(TrilinosWrappers::types::int_type(0));
+ nonlocal_graph->InsertGlobalIndices(row, 1, &row);
+ }
+ Assert(nonlocal_graph->IndicesAreGlobal() == true,
+ ExcInternalError());
+ nonlocal_graph->FillComplete(input_col_map, input_row_map);
+ nonlocal_graph->OptimizeStorage();
+
+ nonlocal_matrix.reset (new Epetra_CrsMatrix(Copy, *nonlocal_graph));
+ }
+ last_action = Zero;
+
+ // In the end, the matrix needs to be compressed in order to be really
+ // ready.
+ compress();
+ }
+
+
+
+
void
SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
{
matrix.reset ();
+ nonlocal_matrix_exporter.reset();
// reinit with a (parallel) Trilinos sparsity pattern.
column_space_map.reset (new Epetra_Map
{
column_space_map.reset (new Epetra_Map (sparse_matrix.domain_partitioner()));
matrix.reset ();
+ nonlocal_matrix_exporter.reset();
matrix.reset (new Epetra_FECrsMatrix
(Copy, sparse_matrix.trilinos_sparsity_pattern(), false));
if (sparse_matrix.nonlocal_matrix != 0)
const Epetra_CrsGraph *graph = &input_matrix.Graph();
nonlocal_matrix.reset();
+ nonlocal_matrix_exporter.reset();
matrix.reset ();
matrix.reset (new Epetra_FECrsMatrix(Copy, *graph, false));
// do only export in case of an add() operation, otherwise the owning
// processor must have set the correct entry
nonlocal_matrix->FillComplete(*column_space_map, matrix->RowMap());
- Epetra_Export exporter(nonlocal_matrix->RowMap(), matrix->RowMap());
- ierr = matrix->Export(*nonlocal_matrix, exporter, mode);
+ if (nonlocal_matrix_exporter.get() == 0)
+ nonlocal_matrix_exporter.reset
+ (new Epetra_Export(nonlocal_matrix->RowMap(), matrix->RowMap()));
+ ierr = matrix->Export(*nonlocal_matrix, *nonlocal_matrix_exporter, mode);
AssertThrow(ierr == 0, ExcTrilinosError(ierr));
ierr = matrix->FillComplete(*column_space_map, matrix->RowMap());
nonlocal_matrix->PutScalar(0);
Utilities::Trilinos::comm_self()));
matrix.reset (new Epetra_FECrsMatrix(View, *column_space_map, 0));
nonlocal_matrix.reset();
+ nonlocal_matrix_exporter.reset();
matrix->FillComplete();
const CompressedSparsityPattern &,
const bool);
template void
- SparseMatrix::reinit (const Epetra_Map &,
- const Epetra_Map &,
- const CompressedSimpleSparsityPattern &,
- const bool);
- template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
const CompressedSetSparsityPattern &,
typename DoFHandler<dim,spacedim>::cell_iterator
neighbor = cell->neighbor(face);
neighbor->get_mg_dof_indices (dofs_on_other_cell);
- // only add one direction
- // The other is taken care of
- // by neighbor.
+ // only add one direction The other is taken care of by
+ // neighbor (except when the neighbor is not owned by the same
+ // processor)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
dofs_on_other_cell[j]);
}
}
+ if (neighbor->is_locally_owned_on_level() == false)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ }
}
}
}
ExcDimensionMismatch (boundary_interface_dofs.size(),
mg_dof_handler.get_tria().n_global_levels()));
- for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
- {
- interface_dofs[l].clear();
- boundary_interface_dofs[l].clear();
- }
+ std::vector<std::vector<types::global_dof_index> >
+ tmp_interface_dofs(interface_dofs.size());
+ std::vector<std::vector<types::global_dof_index> >
+ tmp_boundary_interface_dofs(interface_dofs.size());
const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
}
}
- if (has_coarser_neighbor == true)
- for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
- if (cell->at_boundary(face_nr))
- for (unsigned int j=0; j<dofs_per_face; ++j)
+ if (has_coarser_neighbor == false)
+ continue;
+
+ for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
+ if (cell->at_boundary(face_nr))
+ for (unsigned int j=0; j<dofs_per_face; ++j)
// if (cell_dofs[fe.face_to_cell_index(j,face_nr)] == true) //is this necessary?
- boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+ boundary_cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
const unsigned int level = cell->level();
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
if (cell_dofs[i])
- interface_dofs[level].add_index(local_dof_indices[i]);
+ tmp_interface_dofs[level].push_back(local_dof_indices[i]);
if (boundary_cell_dofs[i])
- boundary_interface_dofs[level].add_index(local_dof_indices[i]);
+ tmp_boundary_interface_dofs[level].push_back(local_dof_indices[i]);
}
}
+
+
+ for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
+ {
+ interface_dofs[l].clear();
+ std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
+ interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
+ std::unique(tmp_interface_dofs[l].begin(),
+ tmp_interface_dofs[l].end()));
+ interface_dofs[l].compress();
+ boundary_interface_dofs[l].clear();
+ std::sort(tmp_boundary_interface_dofs[l].begin(),
+ tmp_boundary_interface_dofs[l].end());
+ boundary_interface_dofs[l].add_indices(tmp_boundary_interface_dofs[l].begin(),
+ std::unique(tmp_boundary_interface_dofs[l].begin(),
+ tmp_boundary_interface_dofs[l].end()));
+ boundary_interface_dofs[l].compress();
+ }
+
}
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1998 - 2013 by the deal.II authors
+// Copyright (C) 1998 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
const std::vector<const InputVector *> &/*solutions*/,
std::vector<Vector<float>*> &/*errors*/,
const ComponentMask &/*component_mask_*/,
- const Function<spacedim> */*coefficient*/,
+ const Function<spacedim> * /*coefficient*/,
const unsigned int,
const types::subdomain_id /*subdomain_id*/,
const types::material_id /*material_id*/)
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include "../tests.h"
+#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/fe/fe_system.h>
+#include <fstream>
+#include <string>
+
+
+template<int dim>
+void print_constant_modes(const FiniteElement<dim> &fe)
+{
+ Table<2,bool> constant_modes = fe.get_constant_modes();
+ deallog << "Testing " << fe.get_name() << std::endl;
+
+ for (unsigned int r=0; r<constant_modes.n_rows(); ++r)
+ {
+ for (unsigned int c=0; c<constant_modes.n_cols(); ++c)
+ deallog << constant_modes(r,c) << " ";
+ deallog << std::endl;
+ }
+ deallog << std::endl;
+}
+
+
+template<int dim>
+void test()
+{
+ print_constant_modes(FE_Q<dim>(1));
+ print_constant_modes(FE_Q<dim>(2));
+ print_constant_modes(FE_DGQ<dim>(1));
+ print_constant_modes(FE_DGP<dim>(2));
+ print_constant_modes(FE_Q_Hierarchical<dim>(1));
+ print_constant_modes(FE_Q_Hierarchical<dim>(2));
+ print_constant_modes(FE_FaceQ<dim>(1));
+ print_constant_modes(FE_FaceP<dim>(1));
+ print_constant_modes(FESystem<dim>(FE_Q<dim>(1), 2, FE_Q<dim>(2), 1));
+ print_constant_modes(FESystem<dim>(FE_DGP<dim>(1), 1, FE_Q_iso_Q1<dim>(2), 1));
+}
+
+template <>
+void test<1>()
+{
+ print_constant_modes(FE_Q<1>(1));
+ print_constant_modes(FESystem<1>(FE_Q<1>(1), 2, FE_Q<1>(2), 1));
+ print_constant_modes(FESystem<1>(FE_DGP<1>(1), 1, FE_Q_iso_Q1<1>(2), 1));
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::Testing FE_Q<1>(1)
+DEAL::1 1
+DEAL::
+DEAL::Testing FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)]
+DEAL::1 0 0 1 0 0 0
+DEAL::0 1 0 0 1 0 0
+DEAL::0 0 1 0 0 1 1
+DEAL::
+DEAL::Testing FESystem<1>[FE_DGP<1>(1)-FE_Q_iso_Q1<1>(2)]
+DEAL::0 0 1 0 0
+DEAL::1 1 0 0 1
+DEAL::
+DEAL::Testing FE_Q<2>(1)
+DEAL::1 1 1 1
+DEAL::
+DEAL::Testing FE_Q<2>(2)
+DEAL::1 1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_DGQ<2>(1)
+DEAL::1 1 1 1
+DEAL::
+DEAL::Testing FE_DGP<2>(2)
+DEAL::1 0 0 0 0 0
+DEAL::
+DEAL::Testing FE_Q_Hierarchical<2>(1)
+DEAL::1 1 1 1
+DEAL::
+DEAL::Testing FE_Q_Hierarchical<2>(2)
+DEAL::1 1 1 1 0 0 0 0 0
+DEAL::
+DEAL::Testing FE_FaceQ<2>(1)
+DEAL::1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_FaceP<2>(1)
+DEAL::1 0 1 0 1 0 1 0
+DEAL::
+DEAL::Testing FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)]
+DEAL::1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0
+DEAL::0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0
+DEAL::0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FESystem<2>[FE_DGP<2>(1)-FE_Q_iso_Q1<2>(2)]
+DEAL::0 0 0 0 0 0 0 0 1 0 0 0
+DEAL::1 1 1 1 1 1 1 1 0 0 0 1
+DEAL::
+DEAL::Testing FE_Q<3>(1)
+DEAL::1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_Q<3>(2)
+DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_DGQ<3>(1)
+DEAL::1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_DGP<3>(2)
+DEAL::1 0 0 0 0 0 0 0 0 0
+DEAL::
+DEAL::Testing FE_Q_Hierarchical<3>(1)
+DEAL::1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_Q_Hierarchical<3>(2)
+DEAL::1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::
+DEAL::Testing FE_FaceQ<3>(1)
+DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FE_FaceP<3>(1)
+DEAL::1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0
+DEAL::
+DEAL::Testing FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)]
+DEAL::1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+DEAL::
+DEAL::Testing FESystem<3>[FE_DGP<3>(1)-FE_Q_iso_Q1<3>(2)]
+DEAL::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
+DEAL::1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1
+DEAL::
void compress(VectorOperation::values op)
{}
+ bool all_zero()
+ {
+ return false;
+ }
+
bool has_ghost_elements()
{
return false;
IndexSet local_active(numproc*2);
local_active.add_range(myid*2,myid*2+2);
IndexSet local_relevant= local_active;
- local_relevant.add_range(1,2);
+ local_relevant.add_range(0,1);
CompressedSimpleSparsityPattern csp (local_relevant);
if (local_relevant.is_element(i))
csp.add(i,i);
- if (myid==0)
- csp.add(0,1);
+ csp.add(0,1);
typename LA::MPI::SparseMatrix mat;
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: vec_00.cc 31681 2013-11-16 08:04:28Z maier $
+//
+// Copyright (C) 2004 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check that all_zero is working correctly (as a collective operation)
+
+#include "../tests.h"
+#include <deal.II/lac/generic_linear_algebra.h>
+#include <deal.II/base/index_set.h>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+
+#include "gla.h"
+
+template <class LA>
+void test ()
+{
+ unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+ unsigned int numproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD);
+
+ if (myid==0)
+ deallog << "numproc=" << numproc << std::endl;
+
+ // each processor owns 2 indices and all
+ // are ghosting Element 1 (the second)
+
+ IndexSet local_active(numproc*2);
+ local_active.add_range(myid*2,myid*2+2);
+ IndexSet local_relevant(numproc*2);
+ local_relevant.add_range(1,2);
+
+
+ typename LA::MPI::Vector x;
+ x.reinit(local_active, MPI_COMM_WORLD);
+ x=0;
+ typename LA::MPI::Vector g(local_active, local_relevant, MPI_COMM_WORLD);
+ g=x;
+ deallog << "all_zero? " << g.all_zero() << " (should be true)" << std::endl;
+ if (myid==0)
+ x(0)=1.0;
+ x.compress(VectorOperation::insert);
+ g=x;
+ deallog << "all_zero? " << g.all_zero() << " (should be false)" << std::endl;
+
+ // done
+ if (myid==0)
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+ MPILogInitAll log;
+ unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+ {
+ deallog.push("PETSc");
+ test<LA_PETSc>();
+ deallog.pop();
+ deallog.push("Trilinos");
+ test<LA_Trilinos>();
+ deallog.pop();
+ }
+
+ if (myid==9999)
+ test<LA_Dummy>(); // don't execute
+
+}
--- /dev/null
+
+DEAL:0:PETSc::numproc=3
+DEAL:0:PETSc::all_zero? 1 (should be true)
+DEAL:0:PETSc::all_zero? 0 (should be false)
+DEAL:0:PETSc::OK
+DEAL:0:Trilinos::numproc=3
+DEAL:0:Trilinos::all_zero? 1 (should be true)
+DEAL:0:Trilinos::all_zero? 0 (should be false)
+DEAL:0:Trilinos::OK
+
+DEAL:1:PETSc::all_zero? 1 (should be true)
+DEAL:1:PETSc::all_zero? 0 (should be false)
+DEAL:1:Trilinos::all_zero? 1 (should be true)
+DEAL:1:Trilinos::all_zero? 0 (should be false)
+
+
+DEAL:2:PETSc::all_zero? 1 (should be true)
+DEAL:2:PETSc::all_zero? 0 (should be false)
+DEAL:2:Trilinos::all_zero? 1 (should be true)
+DEAL:2:Trilinos::all_zero? 0 (should be false)
+
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// cell->get_interpolated_dof_values can not work properly in the hp
+// context when called on non-active cells because these do not have a
+// finite element associated with them
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 2
+
+
+
+
+template <int dim>
+void test ()
+{
+ // create a hp::DoFHandler with different finite elements on the
+ // cells. note that we assign active_fe_indices also to inactive
+ // elements
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr, 0., 1.);
+ tr.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ for (unsigned int i=1; i<5; ++i)
+ fe.push_back (FE_Q<dim>(i));
+
+ hp::DoFHandler<dim> dof_handler(tr);
+ for (typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+ cell!=dof_handler.end(); ++cell)
+ cell->set_active_fe_index (cell->index() % fe.size());
+
+ dof_handler.distribute_dofs (fe);
+
+ // create a mostly arbitrary FE field
+ Vector<double> solution(dof_handler.n_dofs());
+ for (unsigned int i=0; i<solution.size(); ++i)
+ solution(i) = i;
+
+ // try to interpolate from the active cell onto the coarsest cell,
+ // which is definitely not active. this can't work, so expect an
+ // exception
+ typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin(0);
+ Vector<double> local (cell->get_fe().dofs_per_cell);
+
+ try
+ {
+ cell->get_interpolated_dof_values (solution, local);
+ }
+ catch (const ExceptionBase &e)
+ {
+ deallog << "Yes, exception!" << std::endl;
+ deallog << e.get_exc_name() << std::endl;
+ }
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("output");
+ logfile.precision (PRECISION);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ deal_II_exceptions::disable_abort_on_exception();
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects unless you provide an explicit " "finite element index because they do not have naturally " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects unless you provide an explicit " "finite element index because they do not have naturally " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects unless you provide an explicit " "finite element index because they do not have naturally " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// cell->get_interpolated_dof_values can not work properly in the hp
+// context when called on non-active cells because these do not have a
+// finite element associated with them
+//
+// this test verifies that if we call the function on active cells with no
+// explicitly given fe_index that we get the same result as from
+// cell->get_dof_values, and that if we call it with an fe_index for a Q1
+// element that we simply get the vertex dof values
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 2
+
+
+
+
+template <int dim>
+void test ()
+{
+ // create a hp::DoFHandler with different finite elements on the
+ // cells. note that we assign active_fe_indices also to inactive
+ // elements
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr, 0., 1.);
+ tr.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ for (unsigned int i=1; i<5; ++i)
+ fe.push_back (FE_Q<dim>(i));
+
+ hp::DoFHandler<dim> dof_handler(tr);
+ for (typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+ cell!=dof_handler.end(); ++cell)
+ cell->set_active_fe_index (cell->index() % fe.size());
+
+ dof_handler.distribute_dofs (fe);
+
+ // create a mostly arbitrary FE field
+ Vector<double> solution(dof_handler.n_dofs());
+ for (unsigned int i=0; i<solution.size(); ++i)
+ solution(i) = i;
+
+ // do the test
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ {
+ // get values without specifying an explicit fe_index
+ Vector<double> local1 (cell->get_fe().dofs_per_cell);
+ cell->get_interpolated_dof_values (solution, local1);
+
+ // then do the same with the "correct", local fe_index
+ Vector<double> local2 (cell->get_fe().dofs_per_cell);
+ cell->get_interpolated_dof_values (solution, local2,
+ cell->active_fe_index());
+
+ // and do it a third time with the fe_index for a Q1 element
+ Vector<double> local3 (fe[0].dofs_per_cell);
+ cell->get_interpolated_dof_values (solution, local3,
+ 0);
+
+ // now verify correctness
+ Assert (local1 == local2, ExcInternalError());
+
+ // also for the second test. note that vertex dofs always come first in
+ // local1, so we can easily compare
+ for (unsigned int i=0; i<fe[0].dofs_per_cell; ++i)
+ Assert (local1[i] == local3[i],
+ ExcInternalError());
+ }
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("output");
+ logfile.precision (PRECISION);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// cell->get_interpolated_dof_values can not work properly in the hp
+// context when called on non-active cells because these do not have a
+// finite element associated with them
+//
+// this test verifies that if we call the function on active cells with an
+// explicitly given fe_index that we can interpolate from children to the
+// so-specified FE space
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+
+template <int dim>
+void test ()
+{
+ deallog << dim << "D" << std::endl;
+
+ // create a hp::DoFHandler with different finite elements on the
+ // cells. note that we assign active_fe_indices also to inactive
+ // elements
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr, 0., 1.);
+ tr.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ for (unsigned int i=1; i<5; ++i)
+ fe.push_back (FE_Q<dim>(i));
+
+ hp::DoFHandler<dim> dof_handler(tr);
+ for (typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+ cell!=dof_handler.end(); ++cell)
+ cell->set_active_fe_index (cell->index() % fe.size());
+
+ dof_handler.distribute_dofs (fe);
+
+ // create a mostly arbitrary FE field
+ Vector<double> solution(dof_handler.n_dofs());
+ for (unsigned int i=0; i<solution.size(); ++i)
+ solution(i) = i;
+
+ // do the test where we request interpolation onto the coarsest cell with an
+ // explicit Q1 space
+ typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin(0);
+ Vector<double> local (fe[0].dofs_per_cell);
+ cell->get_interpolated_dof_values (solution, local,
+ 0);
+ for (unsigned int i=0; i<local.size(); ++i)
+ deallog << local[i] << ' ';
+ deallog << std::endl;
+
+ // for comparison purposes, also output the values of DoFs on all cells
+ for (typename hp::DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ {
+ Vector<double> x (cell->get_fe().dofs_per_cell);
+ cell->get_dof_values (solution, x);
+ deallog << "cell =" << cell << ": ";
+ for (unsigned int i=0; i<x.size(); ++i)
+ deallog << x[i] << ' ';
+ deallog << std::endl;
+ }
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("output");
+ logfile.precision (1);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::1D
+DEAL::0 7.00000
+DEAL::cell =2.0: 0 1.00000
+DEAL::cell =2.1: 1.00000 2.00000 3.00000
+DEAL::cell =2.2: 2.00000 4.00000 5.00000 6.00000
+DEAL::cell =2.3: 4.00000 7.00000 8.00000 9.00000 10.0000
+DEAL::2D
+DEAL::0 44.0000 91.0000 150.000
+DEAL::cell =2.0: 0 1.00000 2.00000 3.00000
+DEAL::cell =2.1: 1.00000 40.0000 3.00000 42.0000 4.00000 5.00000 6.00000 7.00000 8.00000
+DEAL::cell =2.2: 2.00000 3.00000 82.0000 83.0000 9.00000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
+DEAL::cell =2.3: 3.00000 42.0000 83.0000 125.000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 7.00000 28.0000 29.0000 88.0000 30.0000 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000
+DEAL::cell =2.4: 40.0000 41.0000 42.0000 43.0000
+DEAL::cell =2.5: 41.0000 44.0000 43.0000 45.0000 46.0000 47.0000 48.0000 49.0000 50.0000
+DEAL::cell =2.6: 42.0000 43.0000 125.000 126.000 51.0000 52.0000 53.0000 54.0000 55.0000 56.0000 57.0000 58.0000 59.0000 60.0000 61.0000 62.0000
+DEAL::cell =2.7: 43.0000 45.0000 126.000 129.000 63.0000 64.0000 65.0000 66.0000 67.0000 68.0000 69.0000 49.0000 70.0000 71.0000 133.000 72.0000 73.0000 74.0000 75.0000 76.0000 77.0000 78.0000 79.0000 80.0000 81.0000
+DEAL::cell =2.8: 82.0000 83.0000 84.0000 85.0000
+DEAL::cell =2.9: 83.0000 125.000 85.0000 127.000 86.0000 87.0000 88.0000 89.0000 90.0000
+DEAL::cell =2.10: 84.0000 85.0000 91.0000 92.0000 93.0000 94.0000 95.0000 96.0000 97.0000 98.0000 99.0000 100.000 101.000 102.000 103.000 104.000
+DEAL::cell =2.11: 85.0000 127.000 92.0000 136.000 105.000 106.000 107.000 108.000 109.000 110.000 111.000 89.0000 112.000 113.000 114.000 115.000 116.000 117.000 118.000 119.000 120.000 121.000 122.000 123.000 124.000
+DEAL::cell =2.12: 125.000 126.000 127.000 128.000
+DEAL::cell =2.13: 126.000 129.000 128.000 130.000 131.000 132.000 133.000 134.000 135.000
+DEAL::cell =2.14: 127.000 128.000 136.000 137.000 138.000 139.000 140.000 141.000 142.000 143.000 144.000 145.000 146.000 147.000 148.000 149.000
+DEAL::cell =2.15: 128.000 130.000 137.000 150.000 151.000 152.000 153.000 154.000 155.000 156.000 157.000 134.000 158.000 159.000 160.000 161.000 162.000 163.000 164.000 165.000 166.000 167.000 168.000 169.000 170.000
+DEAL::3D
+DEAL::0 362.000 741.000 1167.00 1603.00 1921.00 2252.00 2621.00
+DEAL::cell =2.0: 0 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000
+DEAL::cell =2.1: 1.00000 354.000 3.00000 356.000 5.00000 358.000 7.00000 360.000 8.00000 9.00000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000
+DEAL::cell =2.2: 2.00000 3.00000 714.000 715.000 6.00000 7.00000 718.000 719.000 27.0000 28.0000 29.0000 30.0000 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000 48.0000 49.0000 50.0000 51.0000 52.0000 53.0000 54.0000 55.0000 56.0000 57.0000 58.0000 59.0000 60.0000 61.0000 62.0000 63.0000 64.0000 65.0000 66.0000 67.0000 68.0000 69.0000 70.0000 71.0000 72.0000 73.0000 74.0000 75.0000 76.0000 77.0000 78.0000 79.0000 80.0000 81.0000 82.0000
+DEAL::cell =2.3: 3.00000 356.000 715.000 1076.00 7.00000 360.000 719.000 1080.00 83.0000 84.0000 85.0000 86.0000 87.0000 88.0000 89.0000 90.0000 91.0000 92.0000 93.0000 94.0000 95.0000 96.0000 97.0000 98.0000 99.0000 100.000 101.000 102.000 103.000 104.000 105.000 106.000 107.000 108.000 109.000 110.000 111.000 112.000 113.000 114.000 115.000 116.000 117.000 118.000 119.000 120.000 121.000 122.000 123.000 124.000 125.000 126.000 127.000 128.000 129.000 130.000 131.000 132.000 133.000 134.000 135.000 136.000 137.000 138.000 139.000 140.000 23.0000 141.000 142.000 143.000 144.000 145.000 146.000 147.000 148.000 736.000 149.000 150.000 151.000 152.000 153.000 154.000 155.000 156.000 157.000 158.000 159.000 160.000 161.000 162.000 163.000 164.000 165.000 166.000 167.000 168.000 169.000 170.000 171.000 172.000 173.000 174.000 175.000 176.000 177.000 178.000 179.000 180.000 181.000 182.000 183.000 184.000 185.000 186.000 187.000 188.000 189.000 190.000 191.000 192.000 193.000 194.000 195.000 196.000 197.000
+DEAL::cell =2.4: 4.00000 5.00000 6.00000 7.00000 198.000 199.000 200.000 201.000
+DEAL::cell =2.5: 5.00000 358.000 7.00000 360.000 199.000 556.000 201.000 558.000 12.0000 13.0000 14.0000 15.0000 202.000 203.000 204.000 205.000 206.000 207.000 208.000 209.000 210.000 211.000 212.000 213.000 25.0000 214.000 215.000
+DEAL::cell =2.6: 6.00000 7.00000 718.000 719.000 200.000 201.000 917.000 918.000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 216.000 217.000 218.000 219.000 220.000 221.000 222.000 223.000 224.000 225.000 226.000 227.000 228.000 229.000 230.000 231.000 232.000 233.000 234.000 235.000 236.000 237.000 238.000 239.000 240.000 241.000 242.000 243.000 244.000 245.000 246.000 247.000 71.0000 72.0000 73.0000 74.0000 248.000 249.000 250.000 251.000 252.000 253.000 254.000 255.000 256.000 257.000 258.000 259.000
+DEAL::cell =2.7: 7.00000 360.000 719.000 1080.00 201.000 558.000 918.000 1285.00 95.0000 96.0000 97.0000 98.0000 99.0000 100.000 101.000 102.000 103.000 104.000 105.000 106.000 260.000 261.000 262.000 263.000 264.000 265.000 266.000 267.000 268.000 269.000 270.000 271.000 272.000 273.000 274.000 275.000 276.000 277.000 278.000 279.000 280.000 281.000 282.000 283.000 284.000 285.000 286.000 287.000 288.000 289.000 290.000 291.000 292.000 293.000 294.000 295.000 296.000 297.000 298.000 299.000 300.000 301.000 302.000 303.000 304.000 305.000 213.000 306.000 307.000 308.000 309.000 310.000 311.000 312.000 313.000 931.000 314.000 315.000 316.000 317.000 162.000 163.000 164.000 165.000 166.000 167.000 168.000 169.000 170.000 318.000 319.000 320.000 321.000 322.000 323.000 324.000 325.000 326.000 327.000 328.000 329.000 330.000 331.000 332.000 333.000 334.000 335.000 336.000 337.000 338.000 339.000 340.000 341.000 342.000 343.000 344.000 345.000 346.000 347.000 348.000 349.000 350.000 351.000 352.000 353.000
+DEAL::cell =2.8: 354.000 355.000 356.000 357.000 358.000 359.000 360.000 361.000
+DEAL::cell =2.9: 355.000 362.000 357.000 363.000 359.000 364.000 361.000 365.000 366.000 367.000 368.000 369.000 370.000 371.000 372.000 373.000 374.000 375.000 376.000 377.000 378.000 379.000 380.000 381.000 382.000 383.000 384.000
+DEAL::cell =2.10: 356.000 357.000 1076.00 1077.00 360.000 361.000 1080.00 1081.00 385.000 386.000 387.000 388.000 389.000 390.000 391.000 392.000 393.000 394.000 395.000 396.000 397.000 398.000 399.000 400.000 401.000 402.000 403.000 404.000 405.000 406.000 407.000 408.000 409.000 410.000 411.000 412.000 413.000 414.000 415.000 416.000 417.000 418.000 419.000 420.000 421.000 422.000 423.000 424.000 425.000 426.000 427.000 428.000 429.000 430.000 431.000 432.000 433.000 434.000 435.000 436.000 437.000 438.000 439.000 440.000
+DEAL::cell =2.11: 357.000 363.000 1077.00 1084.00 361.000 365.000 1081.00 1086.00 441.000 442.000 443.000 444.000 445.000 446.000 447.000 448.000 449.000 450.000 451.000 452.000 453.000 454.000 455.000 456.000 457.000 458.000 459.000 460.000 461.000 462.000 463.000 464.000 465.000 466.000 467.000 468.000 469.000 470.000 471.000 472.000 473.000 474.000 475.000 476.000 477.000 478.000 479.000 480.000 481.000 482.000 483.000 484.000 485.000 486.000 487.000 488.000 489.000 490.000 491.000 492.000 493.000 494.000 495.000 496.000 497.000 498.000 381.000 499.000 500.000 501.000 502.000 503.000 504.000 505.000 506.000 1102.00 507.000 508.000 509.000 510.000 511.000 512.000 513.000 514.000 515.000 516.000 517.000 518.000 519.000 520.000 521.000 522.000 523.000 524.000 525.000 526.000 527.000 528.000 529.000 530.000 531.000 532.000 533.000 534.000 535.000 536.000 537.000 538.000 539.000 540.000 541.000 542.000 543.000 544.000 545.000 546.000 547.000 548.000 549.000 550.000 551.000 552.000 553.000 554.000 555.000
+DEAL::cell =2.12: 358.000 359.000 360.000 361.000 556.000 557.000 558.000 559.000
+DEAL::cell =2.13: 359.000 364.000 361.000 365.000 557.000 560.000 559.000 561.000 370.000 371.000 372.000 373.000 562.000 563.000 564.000 565.000 566.000 567.000 568.000 569.000 570.000 571.000 572.000 573.000 383.000 574.000 575.000
+DEAL::cell =2.14: 360.000 361.000 1080.00 1081.00 558.000 559.000 1285.00 1286.00 393.000 394.000 395.000 396.000 397.000 398.000 399.000 400.000 576.000 577.000 578.000 579.000 580.000 581.000 582.000 583.000 584.000 585.000 586.000 587.000 588.000 589.000 590.000 591.000 592.000 593.000 594.000 595.000 596.000 597.000 598.000 599.000 600.000 601.000 602.000 603.000 604.000 605.000 606.000 607.000 429.000 430.000 431.000 432.000 608.000 609.000 610.000 611.000 612.000 613.000 614.000 615.000 616.000 617.000 618.000 619.000
+DEAL::cell =2.15: 361.000 365.000 1081.00 1086.00 559.000 561.000 1286.00 1289.00 453.000 454.000 455.000 456.000 457.000 458.000 459.000 460.000 461.000 462.000 463.000 464.000 620.000 621.000 622.000 623.000 624.000 625.000 626.000 627.000 628.000 629.000 630.000 631.000 632.000 633.000 634.000 635.000 636.000 637.000 638.000 639.000 640.000 641.000 642.000 643.000 644.000 645.000 646.000 647.000 648.000 649.000 650.000 651.000 652.000 653.000 654.000 655.000 656.000 657.000 658.000 659.000 660.000 661.000 662.000 663.000 664.000 665.000 573.000 666.000 667.000 668.000 669.000 670.000 671.000 672.000 673.000 1301.00 674.000 675.000 676.000 677.000 520.000 521.000 522.000 523.000 524.000 525.000 526.000 527.000 528.000 678.000 679.000 680.000 681.000 682.000 683.000 684.000 685.000 686.000 687.000 688.000 689.000 690.000 691.000 692.000 693.000 694.000 695.000 696.000 697.000 698.000 699.000 700.000 701.000 702.000 703.000 704.000 705.000 706.000 707.000 708.000 709.000 710.000 711.000 712.000 713.000
+DEAL::cell =2.16: 714.000 715.000 716.000 717.000 718.000 719.000 720.000 721.000
+DEAL::cell =2.17: 715.000 1076.00 717.000 1078.00 719.000 1080.00 721.000 1082.00 722.000 723.000 724.000 725.000 726.000 727.000 728.000 729.000 730.000 731.000 732.000 733.000 734.000 735.000 736.000 737.000 738.000 739.000 740.000
+DEAL::cell =2.18: 716.000 717.000 741.000 742.000 720.000 721.000 743.000 744.000 745.000 746.000 747.000 748.000 749.000 750.000 751.000 752.000 753.000 754.000 755.000 756.000 757.000 758.000 759.000 760.000 761.000 762.000 763.000 764.000 765.000 766.000 767.000 768.000 769.000 770.000 771.000 772.000 773.000 774.000 775.000 776.000 777.000 778.000 779.000 780.000 781.000 782.000 783.000 784.000 785.000 786.000 787.000 788.000 789.000 790.000 791.000 792.000 793.000 794.000 795.000 796.000 797.000 798.000 799.000 800.000
+DEAL::cell =2.19: 717.000 1078.00 742.000 1107.00 721.000 1082.00 744.000 1109.00 801.000 802.000 803.000 804.000 805.000 806.000 807.000 808.000 809.000 810.000 811.000 812.000 813.000 814.000 815.000 816.000 817.000 818.000 819.000 820.000 821.000 822.000 823.000 824.000 825.000 826.000 827.000 828.000 829.000 830.000 831.000 832.000 833.000 834.000 835.000 836.000 837.000 838.000 839.000 840.000 841.000 842.000 843.000 844.000 845.000 846.000 847.000 848.000 849.000 850.000 851.000 852.000 853.000 854.000 855.000 856.000 857.000 858.000 737.000 859.000 860.000 861.000 862.000 863.000 864.000 865.000 866.000 867.000 868.000 869.000 870.000 871.000 872.000 873.000 874.000 875.000 876.000 877.000 878.000 879.000 880.000 881.000 882.000 883.000 884.000 885.000 886.000 887.000 888.000 889.000 890.000 891.000 892.000 893.000 894.000 895.000 896.000 897.000 898.000 899.000 900.000 901.000 902.000 903.000 904.000 905.000 906.000 907.000 908.000 909.000 910.000 911.000 912.000 913.000 914.000 915.000 916.000
+DEAL::cell =2.20: 718.000 719.000 720.000 721.000 917.000 918.000 919.000 920.000
+DEAL::cell =2.21: 719.000 1080.00 721.000 1082.00 918.000 1285.00 920.000 1287.00 726.000 727.000 728.000 729.000 921.000 922.000 923.000 924.000 925.000 926.000 927.000 928.000 929.000 930.000 931.000 932.000 739.000 933.000 934.000
+DEAL::cell =2.22: 720.000 721.000 743.000 744.000 919.000 920.000 935.000 936.000 753.000 754.000 755.000 756.000 757.000 758.000 759.000 760.000 937.000 938.000 939.000 940.000 941.000 942.000 943.000 944.000 945.000 946.000 947.000 948.000 949.000 950.000 951.000 952.000 953.000 954.000 955.000 956.000 957.000 958.000 959.000 960.000 961.000 962.000 963.000 964.000 965.000 966.000 967.000 968.000 789.000 790.000 791.000 792.000 969.000 970.000 971.000 972.000 973.000 974.000 975.000 976.000 977.000 978.000 979.000 980.000
+DEAL::cell =2.23: 721.000 1082.00 744.000 1109.00 920.000 1287.00 936.000 1305.00 813.000 814.000 815.000 816.000 817.000 818.000 819.000 820.000 821.000 822.000 823.000 824.000 981.000 982.000 983.000 984.000 985.000 986.000 987.000 988.000 989.000 990.000 991.000 992.000 993.000 994.000 995.000 996.000 997.000 998.000 999.000 1000.00 1001.00 1002.00 1003.00 1004.00 1005.00 1006.00 1007.00 1008.00 1009.00 1010.00 1011.00 1012.00 1013.00 1014.00 1015.00 1016.00 1017.00 1018.00 1019.00 1020.00 1021.00 1022.00 1023.00 1024.00 1025.00 1026.00 932.000 1027.00 1028.00 1029.00 1030.00 1031.00 1032.00 1033.00 1034.00 1035.00 1036.00 1037.00 1038.00 1039.00 881.000 882.000 883.000 884.000 885.000 886.000 887.000 888.000 889.000 1040.00 1041.00 1042.00 1043.00 1044.00 1045.00 1046.00 1047.00 1048.00 1049.00 1050.00 1051.00 1052.00 1053.00 1054.00 1055.00 1056.00 1057.00 1058.00 1059.00 1060.00 1061.00 1062.00 1063.00 1064.00 1065.00 1066.00 1067.00 1068.00 1069.00 1070.00 1071.00 1072.00 1073.00 1074.00 1075.00
+DEAL::cell =2.24: 1076.00 1077.00 1078.00 1079.00 1080.00 1081.00 1082.00 1083.00
+DEAL::cell =2.25: 1077.00 1084.00 1079.00 1085.00 1081.00 1086.00 1083.00 1087.00 1088.00 1089.00 1090.00 1091.00 1092.00 1093.00 1094.00 1095.00 1096.00 1097.00 1098.00 1099.00 1100.00 1101.00 1102.00 1103.00 1104.00 1105.00 1106.00
+DEAL::cell =2.26: 1078.00 1079.00 1107.00 1108.00 1082.00 1083.00 1109.00 1110.00 1111.00 1112.00 1113.00 1114.00 1115.00 1116.00 1117.00 1118.00 1119.00 1120.00 1121.00 1122.00 1123.00 1124.00 1125.00 1126.00 1127.00 1128.00 1129.00 1130.00 1131.00 1132.00 1133.00 1134.00 1135.00 1136.00 1137.00 1138.00 1139.00 1140.00 1141.00 1142.00 1143.00 1144.00 1145.00 1146.00 1147.00 1148.00 1149.00 1150.00 1151.00 1152.00 1153.00 1154.00 1155.00 1156.00 1157.00 1158.00 1159.00 1160.00 1161.00 1162.00 1163.00 1164.00 1165.00 1166.00
+DEAL::cell =2.27: 1079.00 1085.00 1108.00 1167.00 1083.00 1087.00 1110.00 1168.00 1169.00 1170.00 1171.00 1172.00 1173.00 1174.00 1175.00 1176.00 1177.00 1178.00 1179.00 1180.00 1181.00 1182.00 1183.00 1184.00 1185.00 1186.00 1187.00 1188.00 1189.00 1190.00 1191.00 1192.00 1193.00 1194.00 1195.00 1196.00 1197.00 1198.00 1199.00 1200.00 1201.00 1202.00 1203.00 1204.00 1205.00 1206.00 1207.00 1208.00 1209.00 1210.00 1211.00 1212.00 1213.00 1214.00 1215.00 1216.00 1217.00 1218.00 1219.00 1220.00 1221.00 1222.00 1223.00 1224.00 1225.00 1226.00 1103.00 1227.00 1228.00 1229.00 1230.00 1231.00 1232.00 1233.00 1234.00 1235.00 1236.00 1237.00 1238.00 1239.00 1240.00 1241.00 1242.00 1243.00 1244.00 1245.00 1246.00 1247.00 1248.00 1249.00 1250.00 1251.00 1252.00 1253.00 1254.00 1255.00 1256.00 1257.00 1258.00 1259.00 1260.00 1261.00 1262.00 1263.00 1264.00 1265.00 1266.00 1267.00 1268.00 1269.00 1270.00 1271.00 1272.00 1273.00 1274.00 1275.00 1276.00 1277.00 1278.00 1279.00 1280.00 1281.00 1282.00 1283.00 1284.00
+DEAL::cell =2.28: 1080.00 1081.00 1082.00 1083.00 1285.00 1286.00 1287.00 1288.00
+DEAL::cell =2.29: 1081.00 1086.00 1083.00 1087.00 1286.00 1289.00 1288.00 1290.00 1092.00 1093.00 1094.00 1095.00 1291.00 1292.00 1293.00 1294.00 1295.00 1296.00 1297.00 1298.00 1299.00 1300.00 1301.00 1302.00 1105.00 1303.00 1304.00
+DEAL::cell =2.30: 1082.00 1083.00 1109.00 1110.00 1287.00 1288.00 1305.00 1306.00 1119.00 1120.00 1121.00 1122.00 1123.00 1124.00 1125.00 1126.00 1307.00 1308.00 1309.00 1310.00 1311.00 1312.00 1313.00 1314.00 1315.00 1316.00 1317.00 1318.00 1319.00 1320.00 1321.00 1322.00 1323.00 1324.00 1325.00 1326.00 1327.00 1328.00 1329.00 1330.00 1331.00 1332.00 1333.00 1334.00 1335.00 1336.00 1337.00 1338.00 1155.00 1156.00 1157.00 1158.00 1339.00 1340.00 1341.00 1342.00 1343.00 1344.00 1345.00 1346.00 1347.00 1348.00 1349.00 1350.00
+DEAL::cell =2.31: 1083.00 1087.00 1110.00 1168.00 1288.00 1290.00 1306.00 1351.00 1181.00 1182.00 1183.00 1184.00 1185.00 1186.00 1187.00 1188.00 1189.00 1190.00 1191.00 1192.00 1352.00 1353.00 1354.00 1355.00 1356.00 1357.00 1358.00 1359.00 1360.00 1361.00 1362.00 1363.00 1364.00 1365.00 1366.00 1367.00 1368.00 1369.00 1370.00 1371.00 1372.00 1373.00 1374.00 1375.00 1376.00 1377.00 1378.00 1379.00 1380.00 1381.00 1382.00 1383.00 1384.00 1385.00 1386.00 1387.00 1388.00 1389.00 1390.00 1391.00 1392.00 1393.00 1394.00 1395.00 1396.00 1397.00 1302.00 1398.00 1399.00 1400.00 1401.00 1402.00 1403.00 1404.00 1405.00 1406.00 1407.00 1408.00 1409.00 1410.00 1249.00 1250.00 1251.00 1252.00 1253.00 1254.00 1255.00 1256.00 1257.00 1411.00 1412.00 1413.00 1414.00 1415.00 1416.00 1417.00 1418.00 1419.00 1420.00 1421.00 1422.00 1423.00 1424.00 1425.00 1426.00 1427.00 1428.00 1429.00 1430.00 1431.00 1432.00 1433.00 1434.00 1435.00 1436.00 1437.00 1438.00 1439.00 1440.00 1441.00 1442.00 1443.00 1444.00 1445.00 1446.00
+DEAL::cell =2.32: 198.000 199.000 200.000 201.000 1447.00 1448.00 1449.00 1450.00
+DEAL::cell =2.33: 199.000 556.000 201.000 558.000 1448.00 1759.00 1450.00 1761.00 202.000 203.000 204.000 205.000 1451.00 1452.00 1453.00 1454.00 1455.00 1456.00 1457.00 1458.00 1459.00 1460.00 1461.00 1462.00 214.000 1463.00 1464.00
+DEAL::cell =2.34: 200.000 201.000 917.000 918.000 1449.00 1450.00 2075.00 2076.00 216.000 217.000 218.000 219.000 220.000 221.000 222.000 223.000 1465.00 1466.00 1467.00 1468.00 1469.00 1470.00 1471.00 1472.00 1473.00 1474.00 1475.00 1476.00 1477.00 1478.00 1479.00 1480.00 1481.00 1482.00 1483.00 1484.00 1485.00 1486.00 1487.00 1488.00 1489.00 1490.00 1491.00 1492.00 1493.00 1494.00 1495.00 1496.00 248.000 249.000 250.000 251.000 1497.00 1498.00 1499.00 1500.00 1501.00 1502.00 1503.00 1504.00 1505.00 1506.00 1507.00 1508.00
+DEAL::cell =2.35: 201.000 558.000 918.000 1285.00 1450.00 1761.00 2076.00 2393.00 260.000 261.000 262.000 263.000 264.000 265.000 266.000 267.000 268.000 269.000 270.000 271.000 1509.00 1510.00 1511.00 1512.00 1513.00 1514.00 1515.00 1516.00 1517.00 1518.00 1519.00 1520.00 1521.00 1522.00 1523.00 1524.00 1525.00 1526.00 1527.00 1528.00 1529.00 1530.00 1531.00 1532.00 1533.00 1534.00 1535.00 1536.00 1537.00 1538.00 1539.00 1540.00 1541.00 1542.00 1543.00 1544.00 1545.00 1546.00 1547.00 1548.00 1549.00 1550.00 1551.00 1552.00 1553.00 1554.00 1462.00 1555.00 1556.00 1557.00 1558.00 1559.00 1560.00 1561.00 1562.00 2089.00 1563.00 1564.00 1565.00 1566.00 318.000 319.000 320.000 321.000 322.000 323.000 324.000 325.000 326.000 1567.00 1568.00 1569.00 1570.00 1571.00 1572.00 1573.00 1574.00 1575.00 1576.00 1577.00 1578.00 1579.00 1580.00 1581.00 1582.00 1583.00 1584.00 1585.00 1586.00 1587.00 1588.00 1589.00 1590.00 1591.00 1592.00 1593.00 1594.00 1595.00 1596.00 1597.00 1598.00 1599.00 1600.00 1601.00 1602.00
+DEAL::cell =2.36: 1447.00 1448.00 1449.00 1450.00 1603.00 1604.00 1605.00 1606.00
+DEAL::cell =2.37: 1448.00 1759.00 1450.00 1761.00 1604.00 1917.00 1606.00 1919.00 1451.00 1452.00 1453.00 1454.00 1607.00 1608.00 1609.00 1610.00 1611.00 1612.00 1613.00 1614.00 1615.00 1616.00 1617.00 1618.00 1463.00 1619.00 1620.00
+DEAL::cell =2.38: 1449.00 1450.00 2075.00 2076.00 1605.00 1606.00 2234.00 2235.00 1465.00 1466.00 1467.00 1468.00 1469.00 1470.00 1471.00 1472.00 1621.00 1622.00 1623.00 1624.00 1625.00 1626.00 1627.00 1628.00 1629.00 1630.00 1631.00 1632.00 1633.00 1634.00 1635.00 1636.00 1637.00 1638.00 1639.00 1640.00 1641.00 1642.00 1643.00 1644.00 1645.00 1646.00 1647.00 1648.00 1649.00 1650.00 1651.00 1652.00 1497.00 1498.00 1499.00 1500.00 1653.00 1654.00 1655.00 1656.00 1657.00 1658.00 1659.00 1660.00 1661.00 1662.00 1663.00 1664.00
+DEAL::cell =2.39: 1450.00 1761.00 2076.00 2393.00 1606.00 1919.00 2235.00 2555.00 1509.00 1510.00 1511.00 1512.00 1513.00 1514.00 1515.00 1516.00 1517.00 1518.00 1519.00 1520.00 1665.00 1666.00 1667.00 1668.00 1669.00 1670.00 1671.00 1672.00 1673.00 1674.00 1675.00 1676.00 1677.00 1678.00 1679.00 1680.00 1681.00 1682.00 1683.00 1684.00 1685.00 1686.00 1687.00 1688.00 1689.00 1690.00 1691.00 1692.00 1693.00 1694.00 1695.00 1696.00 1697.00 1698.00 1699.00 1700.00 1701.00 1702.00 1703.00 1704.00 1705.00 1706.00 1707.00 1708.00 1709.00 1710.00 1618.00 1711.00 1712.00 1713.00 1714.00 1715.00 1716.00 1717.00 1718.00 2248.00 1719.00 1720.00 1721.00 1722.00 1567.00 1568.00 1569.00 1570.00 1571.00 1572.00 1573.00 1574.00 1575.00 1723.00 1724.00 1725.00 1726.00 1727.00 1728.00 1729.00 1730.00 1731.00 1732.00 1733.00 1734.00 1735.00 1736.00 1737.00 1738.00 1739.00 1740.00 1741.00 1742.00 1743.00 1744.00 1745.00 1746.00 1747.00 1748.00 1749.00 1750.00 1751.00 1752.00 1753.00 1754.00 1755.00 1756.00 1757.00 1758.00
+DEAL::cell =2.40: 556.000 557.000 558.000 559.000 1759.00 1760.00 1761.00 1762.00
+DEAL::cell =2.41: 557.000 560.000 559.000 561.000 1760.00 1763.00 1762.00 1764.00 562.000 563.000 564.000 565.000 1765.00 1766.00 1767.00 1768.00 1769.00 1770.00 1771.00 1772.00 1773.00 1774.00 1775.00 1776.00 574.000 1777.00 1778.00
+DEAL::cell =2.42: 558.000 559.000 1285.00 1286.00 1761.00 1762.00 2393.00 2394.00 576.000 577.000 578.000 579.000 580.000 581.000 582.000 583.000 1779.00 1780.00 1781.00 1782.00 1783.00 1784.00 1785.00 1786.00 1787.00 1788.00 1789.00 1790.00 1791.00 1792.00 1793.00 1794.00 1795.00 1796.00 1797.00 1798.00 1799.00 1800.00 1801.00 1802.00 1803.00 1804.00 1805.00 1806.00 1807.00 1808.00 1809.00 1810.00 608.000 609.000 610.000 611.000 1811.00 1812.00 1813.00 1814.00 1815.00 1816.00 1817.00 1818.00 1819.00 1820.00 1821.00 1822.00
+DEAL::cell =2.43: 559.000 561.000 1286.00 1289.00 1762.00 1764.00 2394.00 2397.00 620.000 621.000 622.000 623.000 624.000 625.000 626.000 627.000 628.000 629.000 630.000 631.000 1823.00 1824.00 1825.00 1826.00 1827.00 1828.00 1829.00 1830.00 1831.00 1832.00 1833.00 1834.00 1835.00 1836.00 1837.00 1838.00 1839.00 1840.00 1841.00 1842.00 1843.00 1844.00 1845.00 1846.00 1847.00 1848.00 1849.00 1850.00 1851.00 1852.00 1853.00 1854.00 1855.00 1856.00 1857.00 1858.00 1859.00 1860.00 1861.00 1862.00 1863.00 1864.00 1865.00 1866.00 1867.00 1868.00 1776.00 1869.00 1870.00 1871.00 1872.00 1873.00 1874.00 1875.00 1876.00 2409.00 1877.00 1878.00 1879.00 1880.00 678.000 679.000 680.000 681.000 682.000 683.000 684.000 685.000 686.000 1881.00 1882.00 1883.00 1884.00 1885.00 1886.00 1887.00 1888.00 1889.00 1890.00 1891.00 1892.00 1893.00 1894.00 1895.00 1896.00 1897.00 1898.00 1899.00 1900.00 1901.00 1902.00 1903.00 1904.00 1905.00 1906.00 1907.00 1908.00 1909.00 1910.00 1911.00 1912.00 1913.00 1914.00 1915.00 1916.00
+DEAL::cell =2.44: 1759.00 1760.00 1761.00 1762.00 1917.00 1918.00 1919.00 1920.00
+DEAL::cell =2.45: 1760.00 1763.00 1762.00 1764.00 1918.00 1921.00 1920.00 1922.00 1765.00 1766.00 1767.00 1768.00 1923.00 1924.00 1925.00 1926.00 1927.00 1928.00 1929.00 1930.00 1931.00 1932.00 1933.00 1934.00 1777.00 1935.00 1936.00
+DEAL::cell =2.46: 1761.00 1762.00 2393.00 2394.00 1919.00 1920.00 2555.00 2556.00 1779.00 1780.00 1781.00 1782.00 1783.00 1784.00 1785.00 1786.00 1937.00 1938.00 1939.00 1940.00 1941.00 1942.00 1943.00 1944.00 1945.00 1946.00 1947.00 1948.00 1949.00 1950.00 1951.00 1952.00 1953.00 1954.00 1955.00 1956.00 1957.00 1958.00 1959.00 1960.00 1961.00 1962.00 1963.00 1964.00 1965.00 1966.00 1967.00 1968.00 1811.00 1812.00 1813.00 1814.00 1969.00 1970.00 1971.00 1972.00 1973.00 1974.00 1975.00 1976.00 1977.00 1978.00 1979.00 1980.00
+DEAL::cell =2.47: 1762.00 1764.00 2394.00 2397.00 1920.00 1922.00 2556.00 2559.00 1823.00 1824.00 1825.00 1826.00 1827.00 1828.00 1829.00 1830.00 1831.00 1832.00 1833.00 1834.00 1981.00 1982.00 1983.00 1984.00 1985.00 1986.00 1987.00 1988.00 1989.00 1990.00 1991.00 1992.00 1993.00 1994.00 1995.00 1996.00 1997.00 1998.00 1999.00 2000.00 2001.00 2002.00 2003.00 2004.00 2005.00 2006.00 2007.00 2008.00 2009.00 2010.00 2011.00 2012.00 2013.00 2014.00 2015.00 2016.00 2017.00 2018.00 2019.00 2020.00 2021.00 2022.00 2023.00 2024.00 2025.00 2026.00 1934.00 2027.00 2028.00 2029.00 2030.00 2031.00 2032.00 2033.00 2034.00 2571.00 2035.00 2036.00 2037.00 2038.00 1881.00 1882.00 1883.00 1884.00 1885.00 1886.00 1887.00 1888.00 1889.00 2039.00 2040.00 2041.00 2042.00 2043.00 2044.00 2045.00 2046.00 2047.00 2048.00 2049.00 2050.00 2051.00 2052.00 2053.00 2054.00 2055.00 2056.00 2057.00 2058.00 2059.00 2060.00 2061.00 2062.00 2063.00 2064.00 2065.00 2066.00 2067.00 2068.00 2069.00 2070.00 2071.00 2072.00 2073.00 2074.00
+DEAL::cell =2.48: 917.000 918.000 919.000 920.000 2075.00 2076.00 2077.00 2078.00
+DEAL::cell =2.49: 918.000 1285.00 920.000 1287.00 2076.00 2393.00 2078.00 2395.00 921.000 922.000 923.000 924.000 2079.00 2080.00 2081.00 2082.00 2083.00 2084.00 2085.00 2086.00 2087.00 2088.00 2089.00 2090.00 933.000 2091.00 2092.00
+DEAL::cell =2.50: 919.000 920.000 935.000 936.000 2077.00 2078.00 2093.00 2094.00 937.000 938.000 939.000 940.000 941.000 942.000 943.000 944.000 2095.00 2096.00 2097.00 2098.00 2099.00 2100.00 2101.00 2102.00 2103.00 2104.00 2105.00 2106.00 2107.00 2108.00 2109.00 2110.00 2111.00 2112.00 2113.00 2114.00 2115.00 2116.00 2117.00 2118.00 2119.00 2120.00 2121.00 2122.00 2123.00 2124.00 2125.00 2126.00 969.000 970.000 971.000 972.000 2127.00 2128.00 2129.00 2130.00 2131.00 2132.00 2133.00 2134.00 2135.00 2136.00 2137.00 2138.00
+DEAL::cell =2.51: 920.000 1287.00 936.000 1305.00 2078.00 2395.00 2094.00 2413.00 981.000 982.000 983.000 984.000 985.000 986.000 987.000 988.000 989.000 990.000 991.000 992.000 2139.00 2140.00 2141.00 2142.00 2143.00 2144.00 2145.00 2146.00 2147.00 2148.00 2149.00 2150.00 2151.00 2152.00 2153.00 2154.00 2155.00 2156.00 2157.00 2158.00 2159.00 2160.00 2161.00 2162.00 2163.00 2164.00 2165.00 2166.00 2167.00 2168.00 2169.00 2170.00 2171.00 2172.00 2173.00 2174.00 2175.00 2176.00 2177.00 2178.00 2179.00 2180.00 2181.00 2182.00 2183.00 2184.00 2090.00 2185.00 2186.00 2187.00 2188.00 2189.00 2190.00 2191.00 2192.00 2193.00 2194.00 2195.00 2196.00 2197.00 1040.00 1041.00 1042.00 1043.00 1044.00 1045.00 1046.00 1047.00 1048.00 2198.00 2199.00 2200.00 2201.00 2202.00 2203.00 2204.00 2205.00 2206.00 2207.00 2208.00 2209.00 2210.00 2211.00 2212.00 2213.00 2214.00 2215.00 2216.00 2217.00 2218.00 2219.00 2220.00 2221.00 2222.00 2223.00 2224.00 2225.00 2226.00 2227.00 2228.00 2229.00 2230.00 2231.00 2232.00 2233.00
+DEAL::cell =2.52: 2075.00 2076.00 2077.00 2078.00 2234.00 2235.00 2236.00 2237.00
+DEAL::cell =2.53: 2076.00 2393.00 2078.00 2395.00 2235.00 2555.00 2237.00 2557.00 2079.00 2080.00 2081.00 2082.00 2238.00 2239.00 2240.00 2241.00 2242.00 2243.00 2244.00 2245.00 2246.00 2247.00 2248.00 2249.00 2091.00 2250.00 2251.00
+DEAL::cell =2.54: 2077.00 2078.00 2093.00 2094.00 2236.00 2237.00 2252.00 2253.00 2095.00 2096.00 2097.00 2098.00 2099.00 2100.00 2101.00 2102.00 2254.00 2255.00 2256.00 2257.00 2258.00 2259.00 2260.00 2261.00 2262.00 2263.00 2264.00 2265.00 2266.00 2267.00 2268.00 2269.00 2270.00 2271.00 2272.00 2273.00 2274.00 2275.00 2276.00 2277.00 2278.00 2279.00 2280.00 2281.00 2282.00 2283.00 2284.00 2285.00 2127.00 2128.00 2129.00 2130.00 2286.00 2287.00 2288.00 2289.00 2290.00 2291.00 2292.00 2293.00 2294.00 2295.00 2296.00 2297.00
+DEAL::cell =2.55: 2078.00 2395.00 2094.00 2413.00 2237.00 2557.00 2253.00 2575.00 2139.00 2140.00 2141.00 2142.00 2143.00 2144.00 2145.00 2146.00 2147.00 2148.00 2149.00 2150.00 2298.00 2299.00 2300.00 2301.00 2302.00 2303.00 2304.00 2305.00 2306.00 2307.00 2308.00 2309.00 2310.00 2311.00 2312.00 2313.00 2314.00 2315.00 2316.00 2317.00 2318.00 2319.00 2320.00 2321.00 2322.00 2323.00 2324.00 2325.00 2326.00 2327.00 2328.00 2329.00 2330.00 2331.00 2332.00 2333.00 2334.00 2335.00 2336.00 2337.00 2338.00 2339.00 2340.00 2341.00 2342.00 2343.00 2249.00 2344.00 2345.00 2346.00 2347.00 2348.00 2349.00 2350.00 2351.00 2352.00 2353.00 2354.00 2355.00 2356.00 2198.00 2199.00 2200.00 2201.00 2202.00 2203.00 2204.00 2205.00 2206.00 2357.00 2358.00 2359.00 2360.00 2361.00 2362.00 2363.00 2364.00 2365.00 2366.00 2367.00 2368.00 2369.00 2370.00 2371.00 2372.00 2373.00 2374.00 2375.00 2376.00 2377.00 2378.00 2379.00 2380.00 2381.00 2382.00 2383.00 2384.00 2385.00 2386.00 2387.00 2388.00 2389.00 2390.00 2391.00 2392.00
+DEAL::cell =2.56: 1285.00 1286.00 1287.00 1288.00 2393.00 2394.00 2395.00 2396.00
+DEAL::cell =2.57: 1286.00 1289.00 1288.00 1290.00 2394.00 2397.00 2396.00 2398.00 1291.00 1292.00 1293.00 1294.00 2399.00 2400.00 2401.00 2402.00 2403.00 2404.00 2405.00 2406.00 2407.00 2408.00 2409.00 2410.00 1303.00 2411.00 2412.00
+DEAL::cell =2.58: 1287.00 1288.00 1305.00 1306.00 2395.00 2396.00 2413.00 2414.00 1307.00 1308.00 1309.00 1310.00 1311.00 1312.00 1313.00 1314.00 2415.00 2416.00 2417.00 2418.00 2419.00 2420.00 2421.00 2422.00 2423.00 2424.00 2425.00 2426.00 2427.00 2428.00 2429.00 2430.00 2431.00 2432.00 2433.00 2434.00 2435.00 2436.00 2437.00 2438.00 2439.00 2440.00 2441.00 2442.00 2443.00 2444.00 2445.00 2446.00 1339.00 1340.00 1341.00 1342.00 2447.00 2448.00 2449.00 2450.00 2451.00 2452.00 2453.00 2454.00 2455.00 2456.00 2457.00 2458.00
+DEAL::cell =2.59: 1288.00 1290.00 1306.00 1351.00 2396.00 2398.00 2414.00 2459.00 1352.00 1353.00 1354.00 1355.00 1356.00 1357.00 1358.00 1359.00 1360.00 1361.00 1362.00 1363.00 2460.00 2461.00 2462.00 2463.00 2464.00 2465.00 2466.00 2467.00 2468.00 2469.00 2470.00 2471.00 2472.00 2473.00 2474.00 2475.00 2476.00 2477.00 2478.00 2479.00 2480.00 2481.00 2482.00 2483.00 2484.00 2485.00 2486.00 2487.00 2488.00 2489.00 2490.00 2491.00 2492.00 2493.00 2494.00 2495.00 2496.00 2497.00 2498.00 2499.00 2500.00 2501.00 2502.00 2503.00 2504.00 2505.00 2410.00 2506.00 2507.00 2508.00 2509.00 2510.00 2511.00 2512.00 2513.00 2514.00 2515.00 2516.00 2517.00 2518.00 1411.00 1412.00 1413.00 1414.00 1415.00 1416.00 1417.00 1418.00 1419.00 2519.00 2520.00 2521.00 2522.00 2523.00 2524.00 2525.00 2526.00 2527.00 2528.00 2529.00 2530.00 2531.00 2532.00 2533.00 2534.00 2535.00 2536.00 2537.00 2538.00 2539.00 2540.00 2541.00 2542.00 2543.00 2544.00 2545.00 2546.00 2547.00 2548.00 2549.00 2550.00 2551.00 2552.00 2553.00 2554.00
+DEAL::cell =2.60: 2393.00 2394.00 2395.00 2396.00 2555.00 2556.00 2557.00 2558.00
+DEAL::cell =2.61: 2394.00 2397.00 2396.00 2398.00 2556.00 2559.00 2558.00 2560.00 2399.00 2400.00 2401.00 2402.00 2561.00 2562.00 2563.00 2564.00 2565.00 2566.00 2567.00 2568.00 2569.00 2570.00 2571.00 2572.00 2411.00 2573.00 2574.00
+DEAL::cell =2.62: 2395.00 2396.00 2413.00 2414.00 2557.00 2558.00 2575.00 2576.00 2415.00 2416.00 2417.00 2418.00 2419.00 2420.00 2421.00 2422.00 2577.00 2578.00 2579.00 2580.00 2581.00 2582.00 2583.00 2584.00 2585.00 2586.00 2587.00 2588.00 2589.00 2590.00 2591.00 2592.00 2593.00 2594.00 2595.00 2596.00 2597.00 2598.00 2599.00 2600.00 2601.00 2602.00 2603.00 2604.00 2605.00 2606.00 2607.00 2608.00 2447.00 2448.00 2449.00 2450.00 2609.00 2610.00 2611.00 2612.00 2613.00 2614.00 2615.00 2616.00 2617.00 2618.00 2619.00 2620.00
+DEAL::cell =2.63: 2396.00 2398.00 2414.00 2459.00 2558.00 2560.00 2576.00 2621.00 2460.00 2461.00 2462.00 2463.00 2464.00 2465.00 2466.00 2467.00 2468.00 2469.00 2470.00 2471.00 2622.00 2623.00 2624.00 2625.00 2626.00 2627.00 2628.00 2629.00 2630.00 2631.00 2632.00 2633.00 2634.00 2635.00 2636.00 2637.00 2638.00 2639.00 2640.00 2641.00 2642.00 2643.00 2644.00 2645.00 2646.00 2647.00 2648.00 2649.00 2650.00 2651.00 2652.00 2653.00 2654.00 2655.00 2656.00 2657.00 2658.00 2659.00 2660.00 2661.00 2662.00 2663.00 2664.00 2665.00 2666.00 2667.00 2572.00 2668.00 2669.00 2670.00 2671.00 2672.00 2673.00 2674.00 2675.00 2676.00 2677.00 2678.00 2679.00 2680.00 2519.00 2520.00 2521.00 2522.00 2523.00 2524.00 2525.00 2526.00 2527.00 2681.00 2682.00 2683.00 2684.00 2685.00 2686.00 2687.00 2688.00 2689.00 2690.00 2691.00 2692.00 2693.00 2694.00 2695.00 2696.00 2697.00 2698.00 2699.00 2700.00 2701.00 2702.00 2703.00 2704.00 2705.00 2706.00 2707.00 2708.00 2709.00 2710.00 2711.00 2712.00 2713.00 2714.00 2715.00 2716.00
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// cell->set_dof_values_by_interpolation can not work properly in the hp
+// context when called on non-active cells because these do not have a
+// finite element associated with them
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#define PRECISION 2
+
+
+
+
+template <int dim>
+void test ()
+{
+ // create a hp::DoFHandler with different finite elements on the
+ // cells. note that we assign active_fe_indices also to inactive
+ // elements
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr, 0., 1.);
+ tr.refine_global (2);
+
+ hp::FECollection<dim> fe;
+ for (unsigned int i=1; i<5; ++i)
+ fe.push_back (FE_Q<dim>(i));
+
+ hp::DoFHandler<dim> dof_handler(tr);
+ for (typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin();
+ cell!=dof_handler.end(); ++cell)
+ cell->set_active_fe_index (cell->index() % fe.size());
+
+ dof_handler.distribute_dofs (fe);
+
+ // create a mostly arbitrary FE field
+ Vector<double> solution(dof_handler.n_dofs());
+ for (unsigned int i=0; i<solution.size(); ++i)
+ solution(i) = i;
+
+ // try to interpolate from the active cell onto the coarsest cell,
+ // which is definitely not active. this can't work, so expect an
+ // exception
+ typename hp::DoFHandler<dim>::cell_iterator cell=dof_handler.begin(0);
+ Vector<double> local (cell->get_fe().dofs_per_cell);
+
+ try
+ {
+ cell->get_interpolated_dof_values (solution, local);
+ }
+ catch (const ExceptionBase &e)
+ {
+ deallog << "Yes, exception!" << std::endl;
+ deallog << e.get_exc_name() << std::endl;
+ }
+}
+
+
+int
+main()
+{
+ std::ofstream logfile ("output");
+ logfile.precision (PRECISION);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ deal_II_exceptions::disable_abort_on_exception();
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects because they do not have " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects because they do not have " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
+DEAL::Yes, exception!
+DEAL::ExcMessage ("You cannot call this function on non-active cells " "of hp::DoFHandler objects because they do not have " "associated finite element spaces associated: degrees " "of freedom are only distributed on active cells for which " "the active_fe_index has been set.")
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests thread safety of parallel Trilinos matrices. Same test as
+// parallel_matrix_assemble_02 but initializing the matrix from
+// CompressedSimpleSparsityPattern instead of a Trilinos sparsity pattern.
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+ namespace Scratch
+ {
+ template <int dim>
+ struct Data
+ {
+ Data (const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature)
+ :
+ fe_values(fe,
+ quadrature,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values)
+ {}
+
+ Data (const Data &data)
+ :
+ fe_values(data.fe_values.get_mapping(),
+ data.fe_values.get_fe(),
+ data.fe_values.get_quadrature(),
+ data.fe_values.get_update_flags())
+ {}
+
+ FEValues<dim> fe_values;
+ };
+ }
+
+ namespace Copy
+ {
+ struct Data
+ {
+ Data(const bool assemble_reference)
+ :
+ assemble_reference(assemble_reference)
+ {}
+ std::vector<types::global_dof_index> local_dof_indices;
+ FullMatrix<double> local_matrix;
+ Vector<double> local_rhs;
+ const bool assemble_reference;
+ };
+ }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+
+ void run ();
+
+private:
+ void setup_system ();
+ void test_equality ();
+ void assemble_reference ();
+ void assemble_test ();
+ void solve ();
+ void create_coarse_grid ();
+ void postprocess ();
+
+ void local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+ Assembly::Scratch::Data<dim> &scratch,
+ Assembly::Copy::Data &data);
+ void copy_local_to_global (const Assembly::Copy::Data &data);
+
+ std::vector<types::global_dof_index>
+ get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const;
+
+ parallel::distributed::Triangulation<dim> triangulation;
+
+ DoFHandler<dim> dof_handler;
+ FE_Q<dim> fe;
+ QGauss<dim> quadrature;
+
+ ConstraintMatrix constraints;
+
+ TrilinosWrappers::SparseMatrix reference_matrix;
+ TrilinosWrappers::SparseMatrix test_matrix;
+
+ TrilinosWrappers::MPI::Vector reference_rhs;
+ TrilinosWrappers::MPI::Vector test_rhs;
+
+ std::vector<std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > graph;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ BoundaryValues () : Function<dim> () {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double sum = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ sum += std::sin(deal_II_numbers::PI*p[d]);
+ return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide () : Function<dim> () {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double product = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ product *= (p[d]+1);
+ return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+ :
+ triangulation (MPI_COMM_WORLD),
+ dof_handler (triangulation),
+ fe (1),
+ quadrature(fe.degree+1)
+{
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const
+{
+ std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+
+ constraints.resolve_indices(local_dof_indices);
+ return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+ reference_matrix.clear();
+ test_matrix.clear();
+ dof_handler.distribute_dofs (fe);
+
+ constraints.clear ();
+
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+ // add boundary conditions as inhomogeneous constraints here, do it after
+ // having added the hanging node constraints in order to be consistent and
+ // skip dofs that are already constrained (i.e., are hanging nodes on the
+ // boundary in 3D). In contrast to step-27, we choose a sine function.
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ constraints);
+ constraints.close ();
+
+ typedef FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> CellFilter;
+ CellFilter begin(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active());
+ CellFilter end(IteratorFilters::LocallyOwnedCell(),dof_handler.end());
+ graph = GraphColoring::make_graph_coloring(begin,end,
+ static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+ (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &)> >
+ (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+ IndexSet locally_owned = dof_handler.locally_owned_dofs();
+ {
+ TrilinosWrappers::SparsityPattern csp;
+ csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD);
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ csp.compress();
+ reference_matrix.reinit (csp);
+ reference_rhs.reinit (locally_owned, MPI_COMM_WORLD);
+ }
+ {
+ IndexSet relevant_set;
+ DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_set);
+ CompressedSimpleSparsityPattern csp(dof_handler.n_dofs(), dof_handler.n_dofs(),
+ relevant_set);
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ test_matrix.reinit (locally_owned, csp, MPI_COMM_WORLD, true);
+ test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true);
+ }
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+ Assembly::Scratch::Data<dim> &scratch,
+ Assembly::Copy::Data &data)
+{
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+ data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ data.local_matrix = 0;
+
+ data.local_rhs.reinit (dofs_per_cell);
+ data.local_rhs = 0;
+
+ scratch.fe_values.reinit (cell);
+
+ const FEValues<dim> &fe_values = scratch.fe_values;
+
+ const RightHandSide<dim> rhs_function;
+
+ for (unsigned int q_point=0;
+ q_point<fe_values.n_quadrature_points;
+ ++q_point)
+ {
+ const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_value *
+ fe_values.JxW(q_point));
+ }
+ }
+
+ data.local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
+{
+ if (data.assemble_reference)
+ constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+ data.local_dof_indices,
+ reference_matrix, reference_rhs);
+ else
+ constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+ data.local_dof_indices,
+ test_matrix, test_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+ reference_matrix = 0;
+ reference_rhs = 0;
+
+ Assembly::Copy::Data copy_data(true);
+ Assembly::Scratch::Data<dim> assembly_data(fe, quadrature);
+
+ for (unsigned int color=0; color<graph.size(); ++color)
+ for (typename std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >::const_iterator p = graph[color].begin();
+ p != graph[color].end(); ++p)
+ {
+ local_assemble(*p, assembly_data, copy_data);
+ copy_local_to_global(copy_data);
+ }
+ reference_matrix.compress(VectorOperation::add);
+ reference_rhs.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+ test_matrix = 0;
+ test_rhs = 0;
+
+ WorkStream::
+ run (graph,
+ std_cxx1x::bind (&LaplaceProblem<dim>::
+ local_assemble,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&LaplaceProblem<dim>::
+ copy_local_to_global,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::Data<dim>(fe, quadrature),
+ Assembly::Copy::Data (false),
+ 2*multithread_info.n_threads(),
+ 1);
+ test_matrix.compress(VectorOperation::add);
+ test_rhs.compress(VectorOperation::add);
+
+ test_matrix.add(-1, reference_matrix);
+
+ // there should not even be roundoff difference between matrices
+ deallog.threshold_double(1.e-30);
+ deallog << "error in matrix: " << test_matrix.frobenius_norm() << std::endl;
+ test_rhs.add(-1., reference_rhs);
+ deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+ estimated_error_per_cell(i) = i;
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<3; ++cycle)
+ {
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube(triangulation, 0, 1);
+ triangulation.refine_global(6);
+ }
+
+ setup_system ();
+
+ assemble_reference ();
+ assemble_test ();
+
+ if (cycle < 2)
+ postprocess ();
+ }
+}
+
+
+
+int main (int argc, char **argv)
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+ {
+ deallog.push("2d");
+ LaplaceProblem<2> laplace_problem;
+ laplace_problem.run ();
+ deallog.pop();
+ }
+}
+
--- /dev/null
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
--- /dev/null
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests thread safety of parallel Trilinos block matrices. Same test as
+// parallel_matrix_assemble_04 but initializing the matrix from
+// BlockCompressedSimpleSparsityPattern instead of a Trilinos sparsity pattern.
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/graph_coloring.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("output");
+
+using namespace dealii;
+
+
+namespace Assembly
+{
+ namespace Scratch
+ {
+ template <int dim>
+ struct Data
+ {
+ Data (const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature)
+ :
+ fe_values(fe,
+ quadrature,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values)
+ {}
+
+ Data (const Data &data)
+ :
+ fe_values(data.fe_values.get_mapping(),
+ data.fe_values.get_fe(),
+ data.fe_values.get_quadrature(),
+ data.fe_values.get_update_flags())
+ {}
+
+ FEValues<dim> fe_values;
+ };
+ }
+
+ namespace Copy
+ {
+ struct Data
+ {
+ Data(const bool assemble_reference)
+ :
+ assemble_reference(assemble_reference)
+ {}
+ std::vector<types::global_dof_index> local_dof_indices;
+ FullMatrix<double> local_matrix;
+ Vector<double> local_rhs;
+ const bool assemble_reference;
+ };
+ }
+}
+
+template <int dim>
+class LaplaceProblem
+{
+public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+
+ void run ();
+
+private:
+ void setup_system ();
+ void test_equality ();
+ void assemble_reference ();
+ void assemble_test ();
+ void solve ();
+ void create_coarse_grid ();
+ void postprocess ();
+
+ void local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+ Assembly::Scratch::Data<dim> &scratch,
+ Assembly::Copy::Data &data);
+ void copy_local_to_global (const Assembly::Copy::Data &data);
+
+ std::vector<types::global_dof_index>
+ get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const;
+
+ parallel::distributed::Triangulation<dim> triangulation;
+
+ DoFHandler<dim> dof_handler;
+ FESystem<dim> fe;
+ QGauss<dim> quadrature;
+
+ ConstraintMatrix constraints;
+
+ TrilinosWrappers::BlockSparseMatrix reference_matrix;
+ TrilinosWrappers::BlockSparseMatrix test_matrix;
+
+ TrilinosWrappers::MPI::BlockVector reference_rhs;
+ TrilinosWrappers::MPI::BlockVector test_rhs;
+
+ std::vector<std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > graph;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ BoundaryValues () : Function<dim> (2) {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double sum = 0;
+ for (unsigned int d=0; d<dim; ++d)
+ sum += std::sin(deal_II_numbers::PI*p[d]);
+ return sum;
+}
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide () : Function<dim> () {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double product = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ product *= (p[d]+1);
+ return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+ :
+ triangulation (MPI_COMM_WORLD),
+ dof_handler (triangulation),
+ fe (FE_Q<dim>(1),1, FE_Q<dim>(2),1),
+ quadrature(3)
+{
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+std::vector<types::global_dof_index>
+LaplaceProblem<dim>::
+get_conflict_indices (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &cell) const
+{
+ std::vector<types::global_dof_index> local_dof_indices(cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices(local_dof_indices);
+
+ constraints.resolve_indices(local_dof_indices);
+ return local_dof_indices;
+}
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+ reference_matrix.clear();
+ test_matrix.clear();
+ dof_handler.distribute_dofs (fe);
+ std::vector<unsigned int> blocks(2,0);
+ blocks[1] = 1;
+ DoFRenumbering::component_wise(dof_handler, blocks);
+
+ constraints.clear ();
+
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+
+ // add boundary conditions as inhomogeneous constraints here, do it after
+ // having added the hanging node constraints in order to be consistent and
+ // skip dofs that are already constrained (i.e., are hanging nodes on the
+ // boundary in 3D). In contrast to step-27, we choose a sine function.
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ constraints);
+ constraints.close ();
+
+ typedef FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> CellFilter;
+ CellFilter begin(IteratorFilters::LocallyOwnedCell(),dof_handler.begin_active());
+ CellFilter end(IteratorFilters::LocallyOwnedCell(),dof_handler.end());
+ graph = GraphColoring::make_graph_coloring(begin,end,
+ static_cast<std_cxx1x::function<std::vector<types::global_dof_index>
+ (FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> const &)> >
+ (std_cxx1x::bind(&LaplaceProblem<dim>::get_conflict_indices, this,std_cxx1x::_1)));
+
+ TrilinosWrappers::BlockSparsityPattern csp(2,2);
+ std::vector<IndexSet> locally_owned(2), relevant_set(2);
+ IndexSet locally_owned_total = dof_handler.locally_owned_dofs(), relevant_total;
+ DoFTools::extract_locally_relevant_dofs (dof_handler, relevant_total);
+
+ std::vector<types::global_dof_index> dofs_per_block (2);
+ DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, blocks);
+ locally_owned[0] = locally_owned_total.get_view(0, dofs_per_block[0]);
+ locally_owned[1] = locally_owned_total.get_view(dofs_per_block[0],
+ dof_handler.n_dofs());
+ relevant_set[0] = relevant_total.get_view(0, dofs_per_block[0]);
+ relevant_set[1] = relevant_total.get_view(dofs_per_block[0],
+ dof_handler.n_dofs());
+
+ {
+ csp.reinit(locally_owned, locally_owned, MPI_COMM_WORLD);
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ csp.compress();
+ reference_matrix.reinit (csp);
+ reference_rhs.reinit (locally_owned, MPI_COMM_WORLD);
+ }
+ {
+ BlockCompressedSimpleSparsityPattern csp(relevant_set);
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ test_matrix.reinit (locally_owned, csp, MPI_COMM_WORLD, true);
+ test_rhs.reinit (locally_owned, relevant_set, MPI_COMM_WORLD, true);
+ }
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::local_assemble (const FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> &cell,
+ Assembly::Scratch::Data<dim> &scratch,
+ Assembly::Copy::Data &data)
+{
+ const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+ data.local_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ data.local_matrix = 0;
+
+ data.local_rhs.reinit (dofs_per_cell);
+ data.local_rhs = 0;
+
+ scratch.fe_values.reinit (cell);
+
+ const FEValues<dim> &fe_values = scratch.fe_values;
+
+ const RightHandSide<dim> rhs_function;
+
+ // this does not make a lot of sense physically but it serves the purpose of
+ // the test well
+ for (unsigned int q_point=0;
+ q_point<fe_values.n_quadrature_points;
+ ++q_point)
+ {
+ const double rhs_value = rhs_function.value(fe_values.quadrature_point(q_point),0);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ data.local_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ data.local_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_value *
+ fe_values.JxW(q_point));
+ }
+ }
+
+ data.local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (data.local_dof_indices);
+}
+
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::copy_local_to_global (const Assembly::Copy::Data &data)
+{
+ if (data.assemble_reference)
+ constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+ data.local_dof_indices,
+ reference_matrix, reference_rhs);
+ else
+ constraints.distribute_local_to_global(data.local_matrix, data.local_rhs,
+ data.local_dof_indices,
+ test_matrix, test_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference ()
+{
+ reference_matrix = 0;
+ reference_rhs = 0;
+
+ Assembly::Copy::Data copy_data(true);
+ Assembly::Scratch::Data<dim> assembly_data(fe, quadrature);
+
+ for (unsigned int color=0; color<graph.size(); ++color)
+ for (typename std::vector<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >::const_iterator p = graph[color].begin();
+ p != graph[color].end(); ++p)
+ {
+ local_assemble(*p, assembly_data, copy_data);
+ copy_local_to_global(copy_data);
+ }
+ reference_matrix.compress(VectorOperation::add);
+ reference_rhs.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test ()
+{
+ test_matrix = 0;
+ test_rhs = 0;
+
+ WorkStream::
+ run (graph,
+ std_cxx1x::bind (&LaplaceProblem<dim>::
+ local_assemble,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&LaplaceProblem<dim>::
+ copy_local_to_global,
+ this,
+ std_cxx1x::_1),
+ Assembly::Scratch::Data<dim>(fe, quadrature),
+ Assembly::Copy::Data (false),
+ 2*multithread_info.n_threads(),
+ 1);
+ test_matrix.compress(VectorOperation::add);
+ test_rhs.compress(VectorOperation::add);
+
+ test_matrix.add(-1, reference_matrix);
+
+ // there should not even be roundoff difference between matrices
+ deallog.threshold_double(1.e-30);
+ double frobenius_norm = 0;
+ for (unsigned int i=0; i<2; ++i)
+ for (unsigned int j=0; j<2; ++j)
+ frobenius_norm += numbers::NumberTraits<double>::abs_square(test_matrix.block(i,j).frobenius_norm());
+ deallog << "error in matrix: " << std::sqrt(frobenius_norm) << std::endl;
+ test_rhs.add(-1., reference_rhs);
+ deallog << "error in vector: " << test_rhs.l2_norm() << std::endl;
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+ estimated_error_per_cell(i) = i;
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<3; ++cycle)
+ {
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_shell(triangulation,
+ Point<dim>(),
+ 0.5, 1., (dim==3) ? 96 : 12, false);
+#ifdef DEBUG
+ triangulation.refine_global(3);
+#else
+ triangulation.refine_global(5);
+#endif
+ }
+
+ setup_system ();
+
+ assemble_reference ();
+ assemble_test ();
+
+ if (cycle < 2)
+ postprocess ();
+ }
+}
+
+
+
+int main (int argc, char **argv)
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ Utilities::MPI::MPI_InitFinalize init(argc, argv, numbers::invalid_unsigned_int);
+
+ {
+ deallog.push("2d");
+ LaplaceProblem<2> laplace_problem;
+ laplace_problem.run ();
+ deallog.pop();
+ }
+}
+
--- /dev/null
+
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
+DEAL:2d::error in matrix: 0
+DEAL:2d::error in vector: 0
projects / dealii-svn.git / commitdiff