// ---------------------------------------------------------------------
//
-// Copyright (C) 1998 - 2017 by the deal.II authors
+// Copyright (C) 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
#ifndef dealii_constraint_matrix_h
#define dealii_constraint_matrix_h
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/exceptions.h>
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/subscriptor.h>
-#include <deal.II/base/template_constraints.h>
-
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/vector_element_access.h>
-
-#include <boost/range/iterator_range.hpp>
-
-#include <set>
-#include <utility>
-#include <vector>
-
+#include <deal.II/lac/affine_constraints.h>
DEAL_II_NAMESPACE_OPEN
-template <int dim, class T>
-class Table;
-template <typename>
-class FullMatrix;
-class SparsityPattern;
-class DynamicSparsityPattern;
-class BlockSparsityPattern;
-class BlockDynamicSparsityPattern;
-template <typename number>
-class SparseMatrix;
-template <typename number>
-class BlockSparseMatrix;
-
-namespace internals
-{
- class GlobalRowsFromLocal;
-}
-
-
-// TODO[WB]: We should have a function of the kind
-// ConstraintMatrix::add_constraint (const size_type constrained_dof,
-// const std::vector<std::pair<size_type, double> > &entries,
-// const double inhomogeneity = 0);
-// rather than building up constraints piecemeal through add_line/add_entry
-// etc. This would also eliminate the possibility of accidentally changing
-// existing constraints into something pointless, see the discussion on the
-// mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.
-
-/**
- * This class implements dealing with linear (possibly inhomogeneous)
- * constraints on degrees of freedom. The concept and origin of such
- * constraints is extensively described in the
- * @ref constraints
- * module. The class is meant to deal with a limited number of constraints
- * relative to the total number of degrees of freedom, for example a few per
- * cent up to maybe 30 per cent; and with a linear combination of <i>M</i>
- * other degrees of freedom where <i>M</i> is also relatively small (no larger
- * than at most around the average number of entries per row of a linear
- * system). It is <em>not</em> meant to describe full rank linear systems.
- *
- * The algorithms used in the implementation of this class are described in
- * some detail in the
- * @ref hp_paper "hp paper".
- * There is also a significant amount of documentation on how to use this
- * class in the
- * @ref constraints
- * module.
- *
- *
- * <h3>Description of constraints</h3>
- *
- * Each "line" in objects of this class corresponds to one constrained degree
- * of freedom, with the number of the line being <i>i</i>, entered by using
- * add_line() or add_lines(). The entries in this line are pairs of the form
- * (<i>j</i>,<i>a<sub>ij</sub></i>), which are added by add_entry() or
- * add_entries(). The organization is essentially a SparsityPattern, but with
- * only a few lines containing nonzero elements, and therefore no data wasted
- * on the others. For each line, which has been added by the mechanism above,
- * an elimination of the constrained degree of freedom of the form
- * @f[
- * x_i = \sum_j a_{ij} x_j + b_i
- * @f]
- * is performed, where <i>b<sub>i</sub></i> is optional and set by
- * set_inhomogeneity(). Thus, if a constraint is formulated for instance as a
- * zero mean value of several degrees of freedom, one of the degrees has to be
- * chosen to be eliminated.
- *
- * Note that the constraints are linear in the <i>x<sub>i</sub></i>, and that
- * there might be a constant (non-homogeneous) term in the constraint. This is
- * exactly the form we need for hanging node constraints, where we need to
- * constrain one degree of freedom in terms of others. There are other
- * conditions of this form possible, for example for implementing mean value
- * conditions as is done in the step-11 tutorial program. The name of the
- * class stems from the fact that these constraints can be represented in
- * matrix form as <b>X</b> <i>x</i> = <i>b</i>, and this object then describes
- * the matrix <b>X</b> (and the vector <i>b</i>; originally, the
- * ConstraintMatrix class was only meant to handle homogenous constraints
- * where <i>b</i>=0, thus the name). The most frequent way to create/fill
- * objects of this type is using the DoFTools::make_hanging_node_constraints()
- * function. The use of these objects is first explained in step-6.
- *
- * Objects of the present type are organized in lines (rows), but only those
- * lines are stored where constraints are present. New constraints are added
- * by adding new lines using the add_line() function, and then populating it
- * using the add_entry() function to a given line, or add_entries() to add
- * more than one entry at a time. The right hand side element, if nonzero, can
- * be set using the set_inhomogeneity() function. After all constraints have
- * been added, you need to call close(), which compresses the storage format
- * and sorts the entries.
- *
- * @note Many of the algorithms this class implements are discussed in the
- * @ref hp_paper.
- * The algorithms are also related to those shown in <i>M. S. Shephard: Linear
- * multipoint constraints applied via transformation as part of a direct
- * stiffness assembly process. Int. J. Numer. Meth. Engrg., vol. 20 (1984),
- * pp. 2107-2112.</i>, with the difference that the algorithms shown there
- * completely eliminated constrained degrees of freedom, whereas we usually
- * keep them as part of the linear system.
- *
- * @ingroup dofs
- * @ingroup constraints
- * @author Wolfgang Bangerth, Martin Kronbichler, 1998, 2004, 2008, 2009
- */
-class ConstraintMatrix : public Subscriptor
-{
-public:
- /**
- * Declare the type for container size.
- */
- typedef types::global_dof_index size_type;
-
- /**
- * An enum that describes what should happen if the two ConstraintMatrix
- * objects involved in a call to the merge() function happen to have
- * constraints on the same degrees of freedom.
- */
- enum MergeConflictBehavior
- {
- /**
- * Throw an exception if the two objects concerned have conflicting
- * constraints on the same degree of freedom.
- */
- no_conflicts_allowed,
-
- /**
- * In an operation <code>cm1.merge(cm2)</code>, if <code>cm1</code> and
- * <code>cm2</code> have constraints on the same degree of freedom, take
- * the one from <code>cm1</code>.
- */
- left_object_wins,
-
- /**
- * In an operation <code>cm1.merge(cm2)</code>, if <code>cm1</code> and
- * <code>cm2</code> have constraints on the same degree of freedom, take
- * the one from <code>cm2</code>.
- */
- right_object_wins
- };
-
- /**
- * Constructor. The supplied IndexSet defines which indices might be
- * constrained inside this ConstraintMatrix. In a calculation with a
- * DoFHandler object based on parallel::distributed::Triangulation or
- * parallel::shared::Triangulation, one should use the set of locally
- * relevant dofs (see
- * @ref GlossLocallyRelevantDof).
- *
- * The given IndexSet allows the ConstraintMatrix to save memory by just not
- * caring about degrees of freedom that are not of importance to the current
- * processor. Alternatively, if no such IndexSet is provided, internal data
- * structures for <i>all</i> possible indices will be created, leading to
- * memory consumption on every processor that is proportional to the
- * <i>overall</i> size of the problem, not just proportional to the size of
- * the portion of the overall problem that is handled by the current
- * processor.
- */
- explicit ConstraintMatrix(const IndexSet &local_constraints = IndexSet());
-
- /**
- * Copy constructor
- */
- explicit ConstraintMatrix(const ConstraintMatrix &constraint_matrix);
-
- /**
- * Move constructor
- */
- ConstraintMatrix(ConstraintMatrix &&constraint_matrix) = default;
-
- /**
- * Copy operator. Like for many other large objects, this operator
- * is deleted to avoid its inadvertent use in places such as
- * accidentally declaring a @p ConstraintMatrix object as a
- * function argument by value, rather than by reference.
- *
- * However, you can use the copy_from() function to explicitly
- * copy ConstraintMatrix objects.
- */
- ConstraintMatrix &
- operator=(const ConstraintMatrix &) = delete;
-
- /**
- * Move assignment operator
- */
- ConstraintMatrix &
- operator=(ConstraintMatrix &&constraint_matrix) = default;
-
- /**
- * Copy the given object to the current one.
- *
- * This function exists because @p operator=() is explicitly
- * disabled.
- */
- void
- copy_from(const ConstraintMatrix &other);
-
- /**
- * clear() the ConstraintMatrix object and supply an IndexSet with lines
- * that may be constrained. This function is only relevant in the
- * distributed case to supply a different IndexSet. Otherwise this routine
- * is equivalent to calling clear(). See the constructor for details.
- */
- void
- reinit(const IndexSet &local_constraints = IndexSet());
-
- /**
- * Determines if we can store a constraint for the given @p line_index. This
- * routine only matters in the distributed case and checks if the IndexSet
- * allows storage of this line. Always returns true if not in the
- * distributed case.
- */
- bool
- can_store_line(const size_type line_index) const;
-
- /**
- * Return the index set describing locally relevant lines if any are
- * present. Note that if no local lines were given, this represents an empty
- * IndexSet, whereas otherwise it contains the global problem size and the
- * local range.
- */
- const IndexSet &
- get_local_lines() const;
-
- /**
- * This function copies the content of @p constraints_in with DoFs that are
- * element of the IndexSet @p filter. Elements that are not present in the
- * IndexSet are ignored. All DoFs will be transformed to local index space
- * of the filter, both the constrained DoFs and the other DoFs these entries
- * are constrained to. The local index space of the filter is a contiguous
- * numbering of all (global) DoFs that are elements in the filter.
- *
- * If, for example, the filter represents the range <tt>[10,20)</tt>, and
- * the constraint matrix @p constraints_in includes the global indices
- * <tt>{7,13,14}</tt>, the indices <tt>{3,4}</tt> are added to the calling
- * constraint matrix (since 13 and 14 are elements in the filter and element
- * 13 is the fourth element in the index, and 14 is the fifth).
- *
- * This function provides an easy way to create a ConstraintMatrix for
- * certain vector components in a vector-valued problem from a full
- * ConstraintMatrix, i.e. extracting a diagonal subblock from a larger
- * ConstraintMatrix. The block is specified by the IndexSet argument.
- */
- void
- add_selected_constraints(const ConstraintMatrix &constraints_in,
- const IndexSet & filter);
-
- /**
- * @name Adding constraints
- * @{
- */
-
- /**
- * Add a new line to the matrix. If the line already exists, then the
- * function simply returns without doing anything.
- */
- void
- add_line(const size_type line);
-
- /**
- * Call the first add_line() function for every index <code>i</code> for
- * which <code>lines[i]</code> is true.
- *
- * This function essentially exists to allow adding several constraints of
- * the form <i>x<sub>i</sub></i>=0 all at once, where the set of indices
- * <i>i</i> for which these constraints should be added are given by the
- * argument of this function. On the other hand, just as if the single-
- * argument add_line() function were called repeatedly, the constraints can
- * later be modified to include linear dependencies using the add_entry()
- * function as well as inhomogeneities using set_inhomogeneity().
- */
- void
- add_lines(const std::vector<bool> &lines);
-
- /**
- * Call the first add_line() function for every index <code>i</code> that
- * appears in the argument.
- *
- * This function essentially exists to allow adding several constraints of
- * the form <i>x<sub>i</sub></i>=0 all at once, where the set of indices
- * <i>i</i> for which these constraints should be added are given by the
- * argument of this function. On the other hand, just as if the single-
- * argument add_line() function were called repeatedly, the constraints can
- * later be modified to include linear dependencies using the add_entry()
- * function as well as inhomogeneities using set_inhomogeneity().
- */
- void
- add_lines(const std::set<size_type> &lines);
-
- /**
- * Call the first add_line() function for every index <code>i</code> that
- * appears in the argument.
- *
- * This function essentially exists to allow adding several constraints of
- * the form <i>x<sub>i</sub></i>=0 all at once, where the set of indices
- * <i>i</i> for which these constraints should be added are given by the
- * argument of this function. On the other hand, just as if the single-
- * argument add_line() function were called repeatedly, the constraints can
- * later be modified to include linear dependencies using the add_entry()
- * function as well as inhomogeneities using set_inhomogeneity().
- */
- void
- add_lines(const IndexSet &lines);
-
- /**
- * Add an entry to a given line. The list of lines is searched from the back
- * to the front, so clever programming would add a new line (which is pushed
- * to the back) and immediately afterwards fill the entries of that line.
- * This way, no expensive searching is needed.
- *
- * If an entry with the same indices as the one this function call denotes
- * already exists, then this function simply returns provided that the value
- * of the entry is the same. Thus, it does no harm to enter a constraint
- * twice.
- */
- void
- add_entry(const size_type line, const size_type column, const double value);
-
- /**
- * Add a whole series of entries, denoted by pairs of column indices and
- * values, to a line of constraints. This function is equivalent to calling
- * the preceding function several times, but is faster.
- */
- void
- add_entries(const size_type line,
- const std::vector<std::pair<size_type, double>> &col_val_pairs);
-
- /**
- * Set an inhomogeneity to the constraint line <i>i</i>, according to the
- * discussion in the general class description.
- *
- * @note the line needs to be added with one of the add_line() calls first.
- */
- void
- set_inhomogeneity(const size_type line, const double value);
-
- /**
- * Close the filling of entries. Since the lines of a matrix of this type
- * are usually filled in an arbitrary order and since we do not want to use
- * associative constrainers to store the lines, we need to sort the lines
- * and within the lines the columns before usage of the matrix. This is done
- * through this function.
- *
- * Also, zero entries are discarded, since they are not needed.
- *
- * After closing, no more entries are accepted. If the object was already
- * closed, then this function returns immediately.
- *
- * This function also resolves chains of constraints. For example, degree of
- * freedom 13 may be constrained to $u_{13} = \frac{u_3}{2} + \frac{u_7}{2}$
- * while degree of freedom 7 is itself constrained as $u_{7} = \frac{u_2}{2}
- * + \frac{u_4}{2}$. Then, the resolution will be that $u_{13} =
- * \frac{u_3}{2} + \frac{u_2}{4} + \frac{u_4}{4}$. Note, however, that
- * cycles in this graph of constraints are not allowed, i.e. for example
- * $u_4$ may not be constrained, directly or indirectly, to $u_{13}$ again.
- */
- void
- close();
-
- /**
- * Merge the constraints represented by the object given as argument into
- * the constraints represented by this object. Both objects may or may not
- * be closed (by having their function close() called before). If this
- * object was closed before, then it will be closed afterwards as well.
- * Note, however, that if the other argument is closed, then merging may be
- * significantly faster.
- *
- * Using the default value of the second arguments, the constraints in each
- * of the two objects (the old one represented by this object and the
- * argument) may not refer to the same degree of freedom, i.e. a degree of
- * freedom that is constrained in one object may not be constrained in the
- * second. If this is nevertheless the case, an exception is thrown.
- * However, this behavior can be changed by providing a different value for
- * the second argument.
- *
- * By default, merging two ConstraintMatrix objects that are initialized
- * with different IndexSet objects is not allowed.
- * This behavior can be altered by setting @p allow_different_local_lines
- * appropriately.
- *
- * Merging a ConstraintMatrix that is initialized with an IndexSet
- * and one that is not initialized with an IndexSet is not yet implemented.
- */
- void
- merge(
- const ConstraintMatrix & other_constraints,
- const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed,
- const bool allow_different_local_lines = false);
-
- /**
- * Shift all entries of this matrix down @p offset rows and over @p offset
- * columns. If this object is initialized with an IndexSet, local_lines are
- * shifted as well.
- *
- * This function is useful if you are building block matrices, where all
- * blocks are built by the same DoFHandler object, i.e. the matrix size is
- * larger than the number of degrees of freedom. Since several matrix rows
- * and columns correspond to the same degrees of freedom, you'd generate
- * several constraint objects, then shift them, and finally merge() them
- * together again.
- */
- void
- shift(const size_type offset);
-
- /**
- * Clear all entries of this matrix. Reset the flag determining whether new
- * entries are accepted or not.
- *
- * This function may be called also on objects which are empty or already
- * cleared.
- */
- void
- clear();
-
- /**
- * @}
- */
-
-
- /**
- * @name Querying constraints
- * @{
- */
-
- /**
- * Return number of constraints stored in this matrix.
- */
- size_type
- n_constraints() const;
-
- /**
- * Return whether the degree of freedom with number @p index is a
- * constrained one.
- *
- * Note that if close() was called before, then this function is
- * significantly faster, since then the constrained degrees of freedom are
- * sorted and we can do a binary search, while before close() was called, we
- * have to perform a linear search through all entries.
- */
- bool
- is_constrained(const size_type index) const;
-
- /**
- * Return whether the dof is constrained, and whether it is constrained to
- * only one other degree of freedom with weight one. The function therefore
- * returns whether the degree of freedom would simply be eliminated in favor
- * of exactly one other degree of freedom.
- *
- * The function returns @p false if either the degree of freedom is not
- * constrained at all, or if it is constrained to more than one other degree
- * of freedom, or if it is constrained to only one degree of freedom but
- * with a weight different from one.
- */
- bool
- is_identity_constrained(const size_type index) const;
-
- /**
- * Return whether the two given degrees of freedom are linked by an equality
- * constraint that either constrains index1 to be so that
- * <code>index1=index2</code> or constrains index2 so that
- * <code>index2=index1</code>.
- */
- bool
- are_identity_constrained(const size_type index1,
- const size_type index2) const;
-
- /**
- * Return the maximum number of other dofs that one dof is constrained to.
- * For example, in 2d a hanging node is constrained only to its two
- * neighbors, so the returned value would be 2. However, for higher order
- * elements and/or higher dimensions, or other types of constraints, this
- * number is no more obvious.
- *
- * The name indicates that within the system matrix, references to a
- * constrained node are indirected to the nodes it is constrained to.
- */
- size_type
- max_constraint_indirections() const;
-
- /**
- * Return <tt>true</tt> in case the dof is constrained and there is a non-
- * trivial inhomogeneous values set to the dof.
- */
- bool
- is_inhomogeneously_constrained(const size_type index) const;
-
- /**
- * Return <tt>false</tt> if all constraints in the ConstraintMatrix are
- * homogeneous ones, and <tt>true</tt> if there is at least one
- * inhomogeneity.
- */
- bool
- has_inhomogeneities() const;
-
- /**
- * Return a pointer to the vector of entries if a line is constrained,
- * and a zero pointer in case the dof is not constrained.
- */
- const std::vector<std::pair<size_type, double>> *
- get_constraint_entries(const size_type line) const;
-
- /**
- * Return the value of the inhomogeneity stored in the constrained dof @p
- * line. Unconstrained dofs also return a zero value.
- */
- double
- get_inhomogeneity(const size_type line) const;
-
- /**
- * Print the constraints represented by the current object to the
- * given stream.
- *
- * For each constraint of the form
- * @f[
- * x_{42} = 0.5 x_2 + 0.25 x_{14} + 2.75
- * @f]
- * this function will write a sequence of lines that look like this:
- * @code
- * 42 2 : 0.5
- * 42 14 : 0.25
- * 42 : 2.75
- * @endcode
- * The last line is only shown if the inhomogeneity (here: 2.75) is
- * nonzero.
- *
- * A block of lines such as the one above is repeated for each
- * constrained degree of freedom.
- */
- void
- print(std::ostream &out) const;
-
- /**
- * Write the graph of constraints in 'dot' format. 'dot' is a program that
- * can take a list of nodes and produce a graphical representation of the
- * graph of constrained degrees of freedom and the degrees of freedom they
- * are constrained to.
- *
- * The output of this function can be used as input to the 'dot' program
- * that can convert the graph into a graphical representation in postscript,
- * png, xfig, and a number of other formats.
- *
- * This function exists mostly for debugging purposes.
- */
- void
- write_dot(std::ostream &) const;
-
- /**
- * Determine an estimate for the memory consumption (in bytes) of this
- * object.
- */
- std::size_t
- memory_consumption() const;
-
- /**
- * Add the constraint indices associated to the indices in the given vector.
- * After a call to this function, the indices vector contains the initial
- * elements and all the associated constrained indices. This function sorts
- * the elements and suppresses duplicates.
- */
- void
- resolve_indices(std::vector<types::global_dof_index> &indices) const;
-
- /**
- * @}
- */
-
- /**
- * @name Eliminating constraints from linear systems after their creation
- * @{
- */
-
-
- /**
- * Condense a sparsity pattern. The name of the function mimics the name of
- * the function we use to condense linear systems, but it is a bit of a
- * misnomer for the current context. This is because in the context of
- * linear systems, we eliminate certain rows and columns of the linear
- * system, i.e., we "reduce" or "condense" the linear system. On the other
- * hand, in the current context, the functions does not remove nonzero
- * entries from the sparsity pattern. Rather, it adds those nonzero entry
- * locations to the sparsity pattern that will later be needed for the
- * process of condensation of constrained degrees of freedom from a linear
- * system.
- *
- * Since this function adds new nonzero entries to the sparsity pattern, the
- * given sparsity pattern must not be compressed. The constraint matrix
- * (i.e., the current object) must be closed. The sparsity pattern is
- * compressed at the end of the function.
- */
- void
- condense(SparsityPattern &sparsity) const;
-
- /**
- * Same function as above, but condenses square block sparsity patterns.
- */
- void
- condense(BlockSparsityPattern &sparsity) const;
-
- /**
- * Same function as above, but condenses square compressed sparsity
- * patterns.
- */
- void
- condense(DynamicSparsityPattern &sparsity) const;
-
- /**
- * Same function as above, but condenses square compressed sparsity
- * patterns.
- */
- void
- condense(BlockDynamicSparsityPattern &sparsity) const;
-
- /**
- * Condense a given matrix, i.e., eliminate the rows and columns of the
- * matrix that correspond to constrained degrees of freedom.
- *
- * See the general documentation of this class for more detailed
- * information.
- */
- template <typename number>
- void
- condense(SparseMatrix<number> &matrix) const;
-
- /**
- * Same function as above, but condenses square block sparse matrices.
- */
- template <typename number>
- void
- condense(BlockSparseMatrix<number> &matrix) const;
-
- /**
- * Condense the given vector in-place. The @p VectorType may be a
- * Vector<float>, Vector<double>, BlockVector<tt><...></tt>, a PETSc or
- * Trilinos vector wrapper class, or any other type having the same
- * interface. Note that this function does not take any inhomogeneity into
- * account and throws an exception in case there are any inhomogeneities.
- * Use the function using both a matrix and vector for that case.
- *
- * @note This function does not work for MPI vectors. Use condense() with
- * two vector arguments instead.
- */
- template <class VectorType>
- void
- condense(VectorType &vec) const;
-
- /**
- * The function copies and condenses values from @p vec_ghosted into @p
- * output. In a serial code it is equivalent to calling condense (vec). If
- * called in parallel, @p vec_ghosted is supposed to contain ghost elements
- * while @p output should not.
- */
- template <class VectorType>
- void
- condense(const VectorType &vec_ghosted, VectorType &output) const;
-
- /**
- * Condense a given matrix and a given vector by eliminating rows and
- * columns of the linear system that correspond to constrained degrees of
- * freedom. The sparsity pattern associated with the matrix needs to be
- * condensed and compressed. This function is the appropriate choice for
- * applying inhomogeneous constraints.
- *
- * The constraint matrix object must be closed to call this function.
- *
- * See the general documentation of this class for more detailed
- * information.
- */
- template <typename number, class VectorType>
- void
- condense(SparseMatrix<number> &matrix, VectorType &vector) const;
-
- /**
- * Same function as above, but condenses square block sparse matrices and
- * vectors.
- */
- template <typename number, class BlockVectorType>
- void
- condense(BlockSparseMatrix<number> &matrix, BlockVectorType &vector) const;
-
- /**
- * Set the values of all constrained DoFs in a vector to zero. The @p
- * VectorType may be a Vector<float>, Vector<double>,
- * BlockVector<tt><...></tt>, a PETSc or Trilinos vector wrapper class, or
- * any other type having the same interface.
- */
- template <class VectorType>
- void
- set_zero(VectorType &vec) const;
-
- /**
- * @}
- */
-
- /**
- * @name Eliminating constraints from linear systems during their creation
- * @{
- */
-
- /**
- * This function takes a vector of local contributions (@p local_vector)
- * corresponding to the degrees of freedom indices given in @p
- * local_dof_indices and distributes them to the global vector. In most
- * cases, these local contributions will be the result of an integration
- * over a cell or face of a cell. However, as long as @p local_vector and @p
- * local_dof_indices have the same number of elements, this function is
- * happy with whatever it is given.
- *
- * In contrast to the similar function in the DoFAccessor class, this
- * function also takes care of constraints, i.e. if one of the elements of
- * @p local_dof_indices belongs to a constrained node, then rather than
- * writing the corresponding element of @p local_vector into @p
- * global_vector, the element is distributed to the entries in the global
- * vector to which this particular degree of freedom is constrained.
- *
- * Thus, by using this function to distribute local contributions to the
- * global object, one saves the call to the condense function after the
- * vectors and matrices are fully assembled. On the other hand, by
- * consequence, the function does not only write into the entries enumerated
- * by the @p local_dof_indices array, but also (possibly) others as
- * necessary.
- *
- * Note that this function will apply all constraints as if they were
- * homogeneous. For correctly setting inhomogeneous constraints, use the
- * similar function with a matrix argument or the function with both matrix
- * and vector arguments.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global vector allows for
- * simultaneous access and the access is not to rows with the same global
- * index at the same time. This needs to be made sure from the caller's
- * site. There is no locking mechanism inside this method to prevent data
- * races.
- *
- * @param[in] local_vector Vector of local contributions.
- * @param[in] local_dof_indices Local degrees of freedom indices
- * corresponding to the vector of local contributions.
- * @param[out] global_vector The global vector to which all local
- * contributions will be added.
- */
- template <class InVector, class OutVector>
- void
- distribute_local_to_global(const InVector & local_vector,
- const std::vector<size_type> &local_dof_indices,
- OutVector & global_vector) const;
-
- /**
- * This function takes a vector of local contributions (@p local_vector)
- * corresponding to the degrees of freedom indices given in @p
- * local_dof_indices and distributes them to the global vector. In most
- * cases, these local contributions will be the result of an integration
- * over a cell or face of a cell. However, as long as @p local_vector and @p
- * local_dof_indices have the same number of elements, this function is
- * happy with whatever it is given.
- *
- * In contrast to the similar function in the DoFAccessor class, this
- * function also takes care of constraints, i.e. if one of the elements of
- * @p local_dof_indices belongs to a constrained node, then rather than
- * writing the corresponding element of @p local_vector into @p
- * global_vector, the element is distributed to the entries in the global
- * vector to which this particular degree of freedom is constrained.
- *
- * Thus, by using this function to distribute local contributions to the
- * global object, one saves the call to the condense function after the
- * vectors and matrices are fully assembled. On the other hand, by
- * consequence, the function does not only write into the entries enumerated
- * by the @p local_dof_indices array, but also (possibly) others as
- * necessary. This includes writing into diagonal elements of the matrix if
- * the corresponding degree of freedom is constrained.
- *
- * The fourth argument <tt>local_matrix</tt> is intended to be used in case
- * one wants to apply inhomogeneous constraints on the vector only. Such a
- * situation could be where one wants to assemble of a right hand side
- * vector on a problem with inhomogeneous constraints, but the global matrix
- * has been assembled previously. A typical example of this is a time
- * stepping algorithm where the stiffness matrix is assembled once, and the
- * right hand side updated every time step. Note that, however, the entries
- * in the columns of the local matrix have to be exactly the same as those
- * that have been written into the global matrix. Otherwise, this function
- * will not be able to correctly handle inhomogeneities.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global vector allows for
- * simultaneous access and the access is not to rows with the same global
- * index at the same time. This needs to be made sure from the caller's
- * site. There is no locking mechanism inside this method to prevent data
- * races.
- */
- template <typename VectorType, typename LocalType>
- void
- distribute_local_to_global(const Vector<LocalType> & local_vector,
- const std::vector<size_type> &local_dof_indices,
- VectorType & global_vector,
- const FullMatrix<LocalType> & local_matrix) const;
-
- /**
- * Same as the previous function, except that it uses two (possibly) different
- * index sets to correctly handle inhomogeneities when the local matrix is
- * computed from a combination of two neighboring elements, for example for an
- * edge integral term in DG. Note that in the case that these two elements
- * have different polynomial degree, the local matrix is rectangular.
- *
- * <tt>local_dof_indices_row</tt> is the set of row indices and
- * <tt>local_dof_indices_col</tt> is the set of column indices of the local
- * matrix. <tt>diagonal=false</tt> says whether the two index sets are equal
- * or not.
- *
- * If both index sets are equal, <tt>diagonal</tt> must be set to true or we
- * simply use the previous function. If both index sets are different
- * (diagonal=false) the <tt>global_vector</tt> is modified to handle
- * inhomogeneities but no entries from <tt>local_vector</tt> are added. Note
- * that the edge integrals for inner edged for DG do not contribute any values
- * to the right hand side.
- */
- template <typename VectorType, typename LocalType>
- void
- distribute_local_to_global(
- const Vector<LocalType> & local_vector,
- const std::vector<size_type> &local_dof_indices_row,
- const std::vector<size_type> &local_dof_indices_col,
- VectorType & global_vector,
- const FullMatrix<LocalType> & local_matrix,
- bool diagonal = false) const;
-
- /**
- * Enter a single value into a result vector, obeying constraints.
- */
- template <class VectorType>
- void
- distribute_local_to_global(const size_type index,
- const double value,
- VectorType & global_vector) const;
-
- /**
- * This function takes a pointer to a vector of local contributions (@p
- * local_vector) corresponding to the degrees of freedom indices given in @p
- * local_dof_indices and distributes them to the global vector. In most
- * cases, these local contributions will be the result of an integration
- * over a cell or face of a cell. However, as long as the entries in @p
- * local_dof_indices indicate reasonable global vector entries, this
- * function is happy with whatever it is given.
- *
- * If one of the elements of @p local_dof_indices belongs to a constrained
- * node, then rather than writing the corresponding element of @p
- * local_vector into @p global_vector, the element is distributed to the
- * entries in the global vector to which this particular degree of freedom
- * is constrained.
- *
- * Thus, by using this function to distribute local contributions to the
- * global object, one saves the call to the condense function after the
- * vectors and matrices are fully assembled. Note that this function
- * completely ignores inhomogeneous constraints.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global vector allows for
- * simultaneous access and the access is not to rows with the same global
- * index at the same time. This needs to be made sure from the caller's
- * site. There is no locking mechanism inside this method to prevent data
- * races.
- */
- template <typename ForwardIteratorVec,
- typename ForwardIteratorInd,
- class VectorType>
- void
- distribute_local_to_global(ForwardIteratorVec local_vector_begin,
- ForwardIteratorVec local_vector_end,
- ForwardIteratorInd local_indices_begin,
- VectorType & global_vector) const;
-
- /**
- * This function takes a matrix of local contributions (@p local_matrix)
- * corresponding to the degrees of freedom indices given in @p
- * local_dof_indices and distributes them to the global matrix. In most
- * cases, these local contributions will be the result of an integration
- * over a cell or face of a cell. However, as long as @p local_matrix and @p
- * local_dof_indices have the same number of elements, this function is
- * happy with whatever it is given.
- *
- * In contrast to the similar function in the DoFAccessor class, this
- * function also takes care of constraints, i.e. if one of the elements of
- * @p local_dof_indices belongs to a constrained node, then rather than
- * writing the corresponding element of @p local_matrix into @p
- * global_matrix, the element is distributed to the entries in the global
- * matrix to which this particular degree of freedom is constrained.
- *
- * With this scheme, we never write into rows or columns of constrained
- * degrees of freedom. In order to make sure that the resulting matrix can
- * still be inverted, we need to do something with the diagonal elements
- * corresponding to constrained nodes. Thus, if a degree of freedom in @p
- * local_dof_indices is constrained, we distribute the corresponding entries
- * in the matrix, but also add the absolute value of the diagonal entry of
- * the local matrix to the corresponding entry in the global matrix.
- * Assuming the discretized operator is positive definite, this guarantees
- * that the diagonal entry is always non-zero, positive, and of the same
- * order of magnitude as the other entries of the matrix. On the other hand,
- * when solving a source problem $Au=f$ the exact value of the diagonal
- * element is not important, since the value of the respective degree of
- * freedom will be overwritten by the distribute() call later on anyway.
- *
- * @note The procedure described above adds an unforeseeable number of
- * artificial eigenvalues to the spectrum of the matrix. Therefore, it is
- * recommended to use the equivalent function with two local index vectors
- * in such a case.
- *
- * By using this function to distribute local contributions to the global
- * object, one saves the call to the condense function after the vectors and
- * matrices are fully assembled.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global matrix allows for
- * simultaneous access and the access is not to rows with the same global
- * index at the same time. This needs to be made sure from the caller's
- * site. There is no locking mechanism inside this method to prevent data
- * races.
- */
- template <typename MatrixType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix) const;
-
- /**
- * Does almost the same as the function above but can treat general
- * rectangular matrices. The main difference to achieve this is that the
- * diagonal entries in constrained rows are left untouched instead of being
- * filled with arbitrary values.
- *
- * Since the diagonal entries corresponding to eliminated degrees of freedom
- * are not set, the result may have a zero eigenvalue, if applied to a
- * square matrix. This has to be considered when solving the resulting
- * problems. For solving a source problem $Au=f$, it is possible to set the
- * diagonal entry after building the matrix by a piece of code of the form
- *
- * @code
- * for (unsigned int i=0;i<matrix.m();++i)
- * if (constraints.is_constrained(i))
- * matrix.diag_element(i) = 1.;
- * @endcode
- *
- * The value of one which is used here is arbitrary, but in the context of
- * Krylov space methods uncritical, since it corresponds to an invariant
- * subspace. If the other matrix entries are smaller or larger by a factor
- * close to machine accuracy, it may be advisable to adjust it.
- *
- * For solving eigenvalue problems, this will only add one spurious zero
- * eigenvalue (with a multiplicity that is possibly greater than one).
- * Taking this into account, nothing else has to be changed.
- */
- template <typename MatrixType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & row_indices,
- const std::vector<size_type> & col_indices,
- MatrixType & global_matrix) const;
-
- /**
- * Does almost the same as the function above for general rectangular
- * matrices but uses different ConstraintMatrix objects on the row and
- * column indices. The convention is that row indices are constrained
- * according to the calling ConstraintMatrix <code>*this</code>, whereas
- * column indices are constrained according to the given ConstraintMatrix
- * <code>column_constraint_matrix</code>. This function allows to handle the
- * case where rows and columns of a matrix are represented by different
- * function spaces with their own enumeration of indices, as e.g. in mixed
- * finite element problems with separate DoFHandler objects or for flux
- * matrices between different levels in multigrid methods.
- *
- * Like the other method with separate slots for row and column indices,
- * this method does not add diagonal entries to eliminated degrees of
- * freedom. See there for a more elaborate description.
- */
- template <typename MatrixType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & row_indices,
- const ConstraintMatrix & column_constraint_matrix,
- const std::vector<size_type> & column_indices,
- MatrixType & global_matrix) const;
-
- /**
- * This function simultaneously writes elements into matrix and vector,
- * according to the constraints specified by the calling ConstraintMatrix.
- * This function can correctly handle inhomogeneous constraints as well. For
- * the parameter use_inhomogeneities_for_rhs see the documentation in
- * @ref constraints
- * module.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global matrix and vector allow
- * for simultaneous access and the access is not to rows with the same
- * global index at the same time. This needs to be made sure from the
- * caller's site. There is no locking mechanism inside this method to
- * prevent data races.
- */
- template <typename MatrixType, typename VectorType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs = false) const;
-
- /**
- * Do a similar operation as the distribute_local_to_global() function that
- * distributes writing entries into a matrix for constrained degrees of
- * freedom, except that here we don't write into a matrix but only allocate
- * sparsity pattern entries.
- *
- * As explained in the
- * @ref hp_paper "hp paper"
- * and in step-27, first allocating a sparsity pattern and later coming back
- * and allocating additional entries for those matrix entries that will be
- * written to due to the elimination of constrained degrees of freedom
- * (using ConstraintMatrix::condense() ), can be a very expensive procedure.
- * It is cheaper to allocate these entries right away without having to do a
- * second pass over the sparsity pattern object. This function does exactly
- * that.
- *
- * Because the function only allocates entries in a sparsity pattern, all it
- * needs to know are the degrees of freedom that couple to each other.
- * Unlike the previous function, no actual values are written, so the second
- * input argument is not necessary here.
- *
- * The third argument to this function, keep_constrained_entries determines
- * whether the function shall allocate entries in the sparsity pattern at
- * all for entries that will later be set to zero upon condensation of the
- * matrix. These entries are necessary if the matrix is built unconstrained,
- * and only later condensed. They are not necessary if the matrix is built
- * using the distribute_local_to_global() function of this class which
- * distributes entries right away when copying a local matrix into a global
- * object. The default of this argument is true, meaning to allocate the few
- * entries that may later be set to zero.
- *
- * By default, the function adds entries for all pairs of indices given in
- * the first argument to the sparsity pattern (unless
- * keep_constrained_entries is false). However, sometimes one would like to
- * only add a subset of all of these pairs. In that case, the last argument
- * can be used which specifies a boolean mask which of the pairs of indices
- * should be considered. If the mask is false for a pair of indices, then no
- * entry will be added to the sparsity pattern for this pair, irrespective
- * of whether one or both of the indices correspond to constrained degrees
- * of freedom.
- *
- * This function is not typically called from user code, but is used in the
- * DoFTools::make_sparsity_pattern() function when passed a constraint
- * matrix object.
- *
- * @note This function in itself is thread-safe, i.e., it works properly
- * also when several threads call it simultaneously. However, the function
- * call is only thread-safe if the underlying global sparsity pattern allows
- * for simultaneous access and the access is not to rows with the same
- * global index at the same time. This needs to be made sure from the
- * caller's site. There is no locking mechanism inside this method to
- * prevent data races.
- */
- template <typename SparsityPatternType>
- void
- add_entries_local_to_global(
- const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries = true,
- const Table<2, bool> & dof_mask = default_empty_table) const;
-
- /**
- * Similar to the other function, but for non-quadratic sparsity patterns.
- */
- template <typename SparsityPatternType>
- void
- add_entries_local_to_global(
- const std::vector<size_type> &row_indices,
- const std::vector<size_type> &col_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries = true,
- const Table<2, bool> & dof_mask = default_empty_table) const;
-
- /**
- * This function imports values from a global vector (@p global_vector) by
- * applying the constraints to a vector of local values, expressed in
- * iterator format. In most cases, the local values will be identified by
- * the local dof values on a cell. However, as long as the entries in @p
- * local_dof_indices indicate reasonable global vector entries, this
- * function is happy with whatever it is given.
- *
- * If one of the elements of @p local_dof_indices belongs to a constrained
- * node, then rather than writing the corresponding element of @p
- * global_vector into @p local_vector, the constraints are resolved as the
- * respective distribute function does, i.e., the local entry is constructed
- * from the global entries to which this particular degree of freedom is
- * constrained.
- *
- * In contrast to the similar function get_dof_values in the DoFAccessor
- * class, this function does not need the constrained values to be correctly
- * set (i.e., distribute to be called).
- */
- template <typename ForwardIteratorVec,
- typename ForwardIteratorInd,
- class VectorType>
- void
- get_dof_values(const VectorType & global_vector,
- ForwardIteratorInd local_indices_begin,
- ForwardIteratorVec local_vector_begin,
- ForwardIteratorVec local_vector_end) const;
-
- /**
- * @}
- */
-
- /**
- * @name Dealing with constraints after solving a linear system
- * @{
- */
-
- /**
- * Given a vector, set all constrained degrees of freedom to values so that
- * the constraints are satisfied. For example, if the current object stores
- * the constraint $x_3=\frac 12 x_1 + \frac 12 x_2$, then this function will
- * read the values of $x_1$ and $x_1$ from the given vector and set the
- * element $x_3$ according to this constraints. Similarly, if the current
- * object stores the constraint $x_{42}=208$, then this function will set
- * the 42nd element of the given vector to 208.
- *
- * @note If this function is called with a parallel vector @p vec, then the
- * vector must not contain ghost elements.
- */
- template <class VectorType>
- void
- distribute(VectorType &vec) const;
-
- /**
- * @}
- */
-
-
-
- /**
- * This class represents one line of a constraint matrix.
- */
- struct ConstraintLine
- {
- /**
- * A data type in which we store the list of entries that make up the
- * homogenous part of a constraint.
- */
- typedef std::vector<std::pair<size_type, double>> Entries;
-
- /**
- * Global DoF index of this line. Since only very few lines are stored,
- * we can not assume a specific order and have to store the index
- * explicitly.
- */
- size_type index;
-
- /**
- * Row numbers and values of the entries in this line.
- *
- * For the reason why we use a vector instead of a map and the
- * consequences thereof, the same applies as what is said for
- * ConstraintMatrix::lines.
- */
- Entries entries;
-
- /**
- * Value of the inhomogeneity.
- */
- double inhomogeneity;
-
- /**
- * This operator is a bit weird and unintuitive: it compares the line
- * numbers of two lines. We need this to sort the lines; in fact we could
- * do this using a comparison predicate. However, this way, it is easier,
- * albeit unintuitive since two lines really have no god-given order
- * relation.
- */
- bool
- operator<(const ConstraintLine &) const;
-
- /**
- * This operator is likewise weird: it checks whether the line indices of
- * the two operands are equal, irrespective of the fact that the contents
- * of the line may be different.
- */
- bool
- operator==(const ConstraintLine &) const;
-
- /**
- * Determine an estimate for the memory consumption (in bytes) of this
- * object.
- */
- std::size_t
- memory_consumption() const;
-
- /**
- * Support for boost:serialization.
- */
- template <class Archive>
- void
- serialize(Archive &ar, const unsigned int)
- {
- ar &index &entries &inhomogeneity;
- }
- };
-
-
- /**
- * Typedef for the iterator type that is used in the LineRange container.
- */
- typedef std::vector<ConstraintLine>::const_iterator const_iterator;
-
-
- /**
- * Typedef for the return type used by get_lines().
- */
- typedef boost::iterator_range<const_iterator> LineRange;
-
-
- /**
- * Return a range object containing (const) iterators to all line entries
- * stored in the ConstraintMatrix. Such a range is useful to initialize
- * range-based for loops as supported by C++11.
- *
- * @return A range object for the half open range <code>[this->begin(),
- * this->end())</code> of line entries.
- */
- const LineRange
- get_lines() const;
-
-
- /**
- * Check if the current object is consistent on all processors
- * in a distributed computation.
- *
- * This method checks if all processors agree on the constraints for their
- * local lines as given by @p locally_active_dofs. This method is a collective
- * operation and will return @p true only if all processors are consistent.
- *
- * Please supply the owned DoFs per processor as returned by
- * DoFHandler::locally_owned_dofs_per_processor() as @p locally_owned_dofs
- * and the result of DoFTools::extract_locally_active_dofs() as
- * @p locally_active_dofs. The
- * former is used to determine ownership of the specific DoF, while the latter
- * is used as the set of rows that need to be checked.
- *
- * If @p verbose is set to @p true, additional debug information is written
- * to std::cout.
- *
- * @note This method exchanges all constraint information of locally active
- * lines and is as such slow for large computations and should probably
- * only be used in debug mode. We do not check all lines returned by
- * get_local_lines() but only the locally active ones, as we allow processors
- * to not know about some locally relevant rows.
- *
- * @return Whether all ConstraintMatrix objects are consistent. Returns the
- * same value on all processors.
- */
- bool
- is_consistent_in_parallel(const std::vector<IndexSet> &locally_owned_dofs,
- const IndexSet & locally_active_dofs,
- const MPI_Comm mpi_communicator,
- const bool verbose = false) const;
-
-
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException0(ExcMatrixIsClosed);
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException0(ExcMatrixNotClosed);
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException1(ExcLineInexistant,
- size_type,
- << "The specified line " << arg1 << " does not exist.");
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException4(ExcEntryAlreadyExists,
- size_type,
- size_type,
- double,
- double,
- << "The entry for the indices " << arg1 << " and " << arg2
- << " already exists, but the values " << arg3 << " (old) and "
- << arg4 << " (new) differ "
- << "by " << (arg4 - arg3) << ".");
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException2(ExcDoFConstrainedToConstrainedDoF,
- int,
- int,
- << "You tried to constrain DoF " << arg1 << " to DoF " << arg2
- << ", but that one is also constrained. This is not allowed!");
- /**
- * Exception.
- *
- * @ingroup Exceptions
- */
- DeclException1(ExcDoFIsConstrainedFromBothObjects,
- size_type,
- << "Degree of freedom " << arg1
- << " is constrained from both object in a merge operation.");
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException1(ExcDoFIsConstrainedToConstrainedDoF,
- size_type,
- << "In the given argument a degree of freedom is constrained "
- << "to another DoF with number " << arg1
- << ", which however is constrained by this object. This is not"
- << " allowed.");
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException1(ExcRowNotStoredHere,
- size_type,
- << "The index set given to this constraint matrix indicates "
- << "constraints for degree of freedom " << arg1
- << " should not be stored by this object, but a constraint "
- << "is being added.");
-
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException2(ExcColumnNotStoredHere,
- size_type,
- size_type,
- << "The index set given to this constraint matrix indicates "
- << "constraints using degree of freedom " << arg2
- << " should not be stored by this object, but a constraint "
- << "for degree of freedom " << arg1 << " uses it.");
-
- /**
- * Exception
- *
- * @ingroup Exceptions
- */
- DeclException2(ExcIncorrectConstraint,
- int,
- int,
- << "While distributing the constraint for DoF " << arg1
- << ", it turns out that one of the processors "
- << "who own the " << arg2 << " degrees of freedom that x_"
- << arg1 << " is constrained against does not know about "
- << "the constraint on x_" << arg1
- << ". Did you not initialize the ConstraintMatrix "
- << "with the appropriate locally_relevant set so "
- << "that every processor who owns a DoF that constrains "
- << "another DoF also knows about this constraint?");
-
-private:
- /**
- * Store the lines of the matrix. Entries are usually appended in an
- * arbitrary order and insertion into a vector is done best at the end, so
- * the order is unspecified after all entries are inserted. Sorting of the
- * entries takes place when calling the <tt>close()</tt> function.
- *
- * We could, instead of using a vector, use an associative array, like a map
- * to store the lines. This, however, would mean a much more fragmented heap
- * since it allocates many small objects, and would additionally make usage
- * of this matrix much slower.
- */
- std::vector<ConstraintLine> lines;
-
- /**
- * A list of size_type that contains the position of the ConstraintLine of a
- * constrained degree of freedom, or numbers::invalid_size_type if the
- * degree of freedom is not constrained. The numbers::invalid_size_type
- * return value returns thus whether there is a constraint line for a given
- * degree of freedom index. Note that this class has no notion of how many
- * degrees of freedom there really are, so if we check whether there is a
- * constraint line for a given degree of freedom, then this vector may
- * actually be shorter than the index of the DoF we check for.
- *
- * This field exists since when adding a new constraint line we have to
- * figure out whether it already exists. Previously, we would simply walk
- * the unsorted list of constraint lines until we either hit the end or
- * found it. This algorithm is O(N) if N is the number of constraints, which
- * makes it O(N^2) when inserting all constraints. For large problems with
- * many constraints, this could easily take 5-10 per cent of the total run
- * time. With this field, we can save this time since we find any constraint
- * in O(1) time or get to know that it a certain degree of freedom is not
- * constrained.
- *
- * To make things worse, traversing the list of existing constraints
- * requires reads from many different places in memory. Thus, in large 3d
- * applications, the add_line() function showed up very prominently in the
- * overall compute time, mainly because it generated a lot of cache misses.
- * This should also be fixed by using the O(1) algorithm to access the
- * fields of this array.
- *
- * The field is useful in a number of other contexts as well, e.g. when one
- * needs random access to the constraints as in all the functions that apply
- * constraints on the fly while add cell contributions into vectors and
- * matrices.
- */
- std::vector<size_type> lines_cache;
-
- /**
- * This IndexSet is used to limit the lines to save in the ConstraintMatrix
- * to a subset. This is necessary, because the lines_cache vector would
- * become too big in a distributed calculation.
- */
- IndexSet local_lines;
-
- /**
- * Store whether the arrays are sorted. If so, no new entries can be added.
- */
- bool sorted;
-
- /**
- * Internal function to calculate the index of line @p line in the vector
- * lines_cache using local_lines.
- */
- size_type
- calculate_line_index(const size_type line) const;
-
- /**
- * Return @p true if the weight of an entry (the second element of the pair)
- * equals zero. This function is used to delete entries with zero weight.
- */
- static bool
- check_zero_weight(const std::pair<size_type, double> &p);
-
- /**
- * Dummy table that serves as default argument for function
- * <tt>add_entries_local_to_global()</tt>.
- */
- static const Table<2, bool> default_empty_table;
-
- /**
- * This function actually implements the local_to_global function for
- * standard (non-block) matrices.
- */
- template <typename MatrixType, typename VectorType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs,
- std::integral_constant<bool, false>) const;
-
- /**
- * This function actually implements the local_to_global function for block
- * matrices.
- */
- template <typename MatrixType, typename VectorType>
- void
- distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs,
- std::integral_constant<bool, true>) const;
-
- /**
- * This function actually implements the local_to_global function for
- * standard (non-block) sparsity types.
- */
- template <typename SparsityPatternType>
- void
- add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> &dof_mask,
- std::integral_constant<bool, false>) const;
-
- /**
- * This function actually implements the local_to_global function for block
- * sparsity types.
- */
- template <typename SparsityPatternType>
- void
- add_entries_local_to_global(const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> &dof_mask,
- std::integral_constant<bool, true>) const;
-
- /**
- * Internal helper function for distribute_local_to_global function.
- *
- * Creates a list of affected global rows for distribution, including the
- * local rows where the entries come from. The list is sorted according to
- * the global row indices.
- */
- void
- make_sorted_row_list(const std::vector<size_type> & local_dof_indices,
- internals::GlobalRowsFromLocal &global_rows) const;
-
- /**
- * Internal helper function for add_entries_local_to_global function.
- *
- * Creates a list of affected rows for distribution without any additional
- * information, otherwise similar to the other make_sorted_row_list()
- * function.
- */
- void
- make_sorted_row_list(const std::vector<size_type> &local_dof_indices,
- std::vector<size_type> & active_dofs) const;
-
- /**
- * Internal helper function for distribute_local_to_global function.
- */
- template <typename MatrixScalar, typename VectorScalar>
- typename ProductType<VectorScalar, MatrixScalar>::type
- resolve_vector_entry(const size_type i,
- const internals::GlobalRowsFromLocal &global_rows,
- const Vector<VectorScalar> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- const FullMatrix<MatrixScalar> &local_matrix) const;
-};
-
-
-
-/* ---------------- template and inline functions ----------------- */
-
-inline ConstraintMatrix::ConstraintMatrix(const IndexSet &local_constraints) :
- lines(),
- local_lines(local_constraints),
- sorted(false)
-{
- // make sure the IndexSet is compressed. Otherwise this can lead to crashes
- // that are hard to find (only happen in release mode).
- // see tests/mpi/constraint_matrix_crash_01
- local_lines.compress();
-}
-
-
-
-inline ConstraintMatrix::ConstraintMatrix(
- const ConstraintMatrix &constraint_matrix) :
- Subscriptor(),
- lines(constraint_matrix.lines),
- lines_cache(constraint_matrix.lines_cache),
- local_lines(constraint_matrix.local_lines),
- sorted(constraint_matrix.sorted)
-{}
-
-
-inline void
-ConstraintMatrix::add_line(const size_type line)
-{
- Assert(sorted == false, ExcMatrixIsClosed());
-
- // the following can happen when we compute with distributed meshes and dof
- // handlers and we constrain a degree of freedom whose number we don't have
- // locally. if we don't abort here the program will try to allocate several
- // terabytes of memory to resize the various arrays below :-)
- Assert(line != numbers::invalid_size_type, ExcInternalError());
- const size_type line_index = calculate_line_index(line);
-
- // check whether line already exists; it may, in which case we can just quit
- if (is_constrained(line))
- return;
-
- // if necessary enlarge vector of existing entries for cache
- if (line_index >= lines_cache.size())
- lines_cache.resize(
- std::max(2 * static_cast<size_type>(lines_cache.size()), line_index + 1),
- numbers::invalid_size_type);
-
- // push a new line to the end of the list
- lines.emplace_back();
- lines.back().index = line;
- lines.back().inhomogeneity = 0.;
- lines_cache[line_index] = lines.size() - 1;
-}
-
-
-
-inline void
-ConstraintMatrix::add_entry(const size_type line,
- const size_type column,
- const double value)
-{
- Assert(sorted == false, ExcMatrixIsClosed());
- Assert(line != column,
- ExcMessage("Can't constrain a degree of freedom to itself"));
-
- // Ensure that the current line is present in the cache:
- const size_type line_index = calculate_line_index(line);
- Assert(line_index < lines_cache.size(),
- ExcMessage("The current ConstraintMatrix does not contain the line "
- "for the current entry. Call ConstraintMatrix::add_line "
- "before calling this function."));
-
- // if in debug mode, check whether an entry for this column already exists
- // and if it's the same as the one entered at present
- //
- // in any case: exit the function if an entry for this column already
- // exists, since we don't want to enter it twice
- Assert(lines_cache[line_index] != numbers::invalid_size_type,
- ExcInternalError());
- Assert(!local_lines.size() || local_lines.is_element(column),
- ExcColumnNotStoredHere(line, column));
- ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
- Assert(line_ptr->index == line, ExcInternalError());
- for (ConstraintLine::Entries::const_iterator p = line_ptr->entries.begin();
- p != line_ptr->entries.end();
- ++p)
- if (p->first == column)
- {
- Assert(std::fabs(p->second - value) < 1.e-14,
- ExcEntryAlreadyExists(line, column, p->second, value));
- return;
- }
-
- line_ptr->entries.emplace_back(column, value);
-}
-
-
-
-inline void
-ConstraintMatrix::set_inhomogeneity(const size_type line, const double value)
-{
- const size_type line_index = calculate_line_index(line);
- Assert(line_index < lines_cache.size() &&
- lines_cache[line_index] != numbers::invalid_size_type,
- ExcMessage("call add_line() before calling set_inhomogeneity()"));
- Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
- ConstraintLine *line_ptr = &lines[lines_cache[line_index]];
- line_ptr->inhomogeneity = value;
-}
-
-
-
-inline types::global_dof_index
-ConstraintMatrix::n_constraints() const
-{
- return lines.size();
-}
-
-
-
-inline bool
-ConstraintMatrix::is_constrained(const size_type index) const
-{
- const size_type line_index = calculate_line_index(index);
- return ((line_index < lines_cache.size()) &&
- (lines_cache[line_index] != numbers::invalid_size_type));
-}
-
-
-
-inline bool
-ConstraintMatrix::is_inhomogeneously_constrained(const size_type index) const
-{
- // check whether the entry is constrained. could use is_constrained, but
- // that means computing the line index twice
- const size_type line_index = calculate_line_index(index);
- if (line_index >= lines_cache.size() ||
- lines_cache[line_index] == numbers::invalid_size_type)
- return false;
- else
- {
- Assert(lines_cache[line_index] < lines.size(), ExcInternalError());
- return !(lines[lines_cache[line_index]].inhomogeneity == 0);
- }
-}
-
-
-
-inline const std::vector<std::pair<types::global_dof_index, double>> *
-ConstraintMatrix::get_constraint_entries(const size_type line) const
-{
- // check whether the entry is constrained. could use is_constrained, but
- // that means computing the line index twice
- const size_type line_index = calculate_line_index(line);
- if (line_index >= lines_cache.size() ||
- lines_cache[line_index] == numbers::invalid_size_type)
- return nullptr;
- else
- return &lines[lines_cache[line_index]].entries;
-}
-
-
-
-inline double
-ConstraintMatrix::get_inhomogeneity(const size_type line) const
-{
- // check whether the entry is constrained. could use is_constrained, but
- // that means computing the line index twice
- const size_type line_index = calculate_line_index(line);
- if (line_index >= lines_cache.size() ||
- lines_cache[line_index] == numbers::invalid_size_type)
- return 0;
- else
- return lines[lines_cache[line_index]].inhomogeneity;
-}
-
-
-
-inline types::global_dof_index
-ConstraintMatrix::calculate_line_index(const size_type line) const
-{
- // IndexSet is unused (serial case)
- if (!local_lines.size())
- return line;
-
- Assert(local_lines.is_element(line), ExcRowNotStoredHere(line));
-
- return local_lines.index_within_set(line);
-}
-
-
-
-inline bool
-ConstraintMatrix::can_store_line(size_type line_index) const
-{
- return !local_lines.size() || local_lines.is_element(line_index);
-}
-
-
-
-inline const IndexSet &
-ConstraintMatrix::get_local_lines() const
-{
- return local_lines;
-}
-
-
-
-template <class VectorType>
-inline void
-ConstraintMatrix::distribute_local_to_global(const size_type index,
- const double value,
- VectorType &global_vector) const
-{
- Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
-
- if (is_constrained(index) == false)
- global_vector(index) += value;
- else
- {
- const ConstraintLine &position =
- lines[lines_cache[calculate_line_index(index)]];
- for (size_type j = 0; j < position.entries.size(); ++j)
- global_vector(position.entries[j].first) +=
- value * position.entries[j].second;
- }
-}
-
-
-template <typename ForwardIteratorVec,
- typename ForwardIteratorInd,
- class VectorType>
-inline void
-ConstraintMatrix::distribute_local_to_global(
- ForwardIteratorVec local_vector_begin,
- ForwardIteratorVec local_vector_end,
- ForwardIteratorInd local_indices_begin,
- VectorType & global_vector) const
-{
- Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
- for (; local_vector_begin != local_vector_end;
- ++local_vector_begin, ++local_indices_begin)
- {
- if (is_constrained(*local_indices_begin) == false)
- internal::ElementAccess<VectorType>::add(
- *local_vector_begin, *local_indices_begin, global_vector);
- else
- {
- const ConstraintLine &position =
- lines[lines_cache[calculate_line_index(*local_indices_begin)]];
- for (size_type j = 0; j < position.entries.size(); ++j)
- internal::ElementAccess<VectorType>::add(
- (*local_vector_begin) * position.entries[j].second,
- position.entries[j].first,
- global_vector);
- }
- }
-}
-
-
-template <class InVector, class OutVector>
-inline void
-ConstraintMatrix::distribute_local_to_global(
- const InVector & local_vector,
- const std::vector<size_type> &local_dof_indices,
- OutVector & global_vector) const
-{
- Assert(local_vector.size() == local_dof_indices.size(),
- ExcDimensionMismatch(local_vector.size(), local_dof_indices.size()));
- distribute_local_to_global(local_vector.begin(),
- local_vector.end(),
- local_dof_indices.begin(),
- global_vector);
-}
-
-
-
-template <typename ForwardIteratorVec,
- typename ForwardIteratorInd,
- class VectorType>
-inline void
-ConstraintMatrix::get_dof_values(const VectorType & global_vector,
- ForwardIteratorInd local_indices_begin,
- ForwardIteratorVec local_vector_begin,
- ForwardIteratorVec local_vector_end) const
-{
- Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
- for (; local_vector_begin != local_vector_end;
- ++local_vector_begin, ++local_indices_begin)
- {
- if (is_constrained(*local_indices_begin) == false)
- *local_vector_begin = global_vector(*local_indices_begin);
- else
- {
- const ConstraintLine &position =
- lines[lines_cache[calculate_line_index(*local_indices_begin)]];
- typename VectorType::value_type value = position.inhomogeneity;
- for (size_type j = 0; j < position.entries.size(); ++j)
- value += (global_vector(position.entries[j].first) *
- position.entries[j].second);
- *local_vector_begin = value;
- }
- }
-}
-
-
-template <typename MatrixType>
-class BlockMatrixBase;
-template <typename SparsityPatternType>
-class BlockSparsityPatternBase;
-template <typename number>
-class BlockSparseMatrixEZ;
-
/**
- * A class that can be used to determine whether a given type is a block
- * matrix type or not. For example,
- * @code
- * IsBlockMatrix<SparseMatrix<double> >::value
- * @endcode
- * has the value false, whereas
- * @code
- * IsBlockMatrix<BlockSparseMatrix<double> >::value
- * @endcode
- * is true. This is sometimes useful in template contexts where we may want to
- * do things differently depending on whether a template type denotes a
- * regular or a block matrix type.
+ * ConstraintMatrix has been renamed to AffineConstraints. Provide a
+ * compatibility typedef that defaults to AffineConstraints<double>.
*
- * @see
- * @ref GlossBlockLA "Block (linear algebra)"
- * @author Wolfgang Bangerth, 2009
+ * @deprecated Use AffineConstraints
*/
-template <typename MatrixType>
-struct IsBlockMatrix
-{
-private:
- struct yes_type
- {
- char c[1];
- };
- struct no_type
- {
- char c[2];
- };
-
- /**
- * Overload returning true if the class is derived from BlockMatrixBase,
- * which is what block matrices do (with the exception of
- * BlockSparseMatrixEZ).
- */
- template <typename T>
- static yes_type
- check_for_block_matrix(const BlockMatrixBase<T> *);
-
- /**
- * Overload returning true if the class is derived from
- * BlockSparsityPatternBase, which is what block sparsity patterns do.
- */
- template <typename T>
- static yes_type
- check_for_block_matrix(const BlockSparsityPatternBase<T> *);
-
- /**
- * Overload for BlockSparseMatrixEZ, which is the only block matrix not
- * derived from BlockMatrixBase at the time of writing this class.
- */
- template <typename T>
- static yes_type
- check_for_block_matrix(const BlockSparseMatrixEZ<T> *);
-
- /**
- * Catch all for all other potential matrix types that are not block
- * matrices.
- */
- static no_type
- check_for_block_matrix(...);
-
-public:
- /**
- * A statically computable value that indicates whether the template
- * argument to this class is a block matrix (in fact whether the type is
- * derived from BlockMatrixBase<T>).
- */
- static const bool value =
- (sizeof(check_for_block_matrix((MatrixType *)nullptr)) == sizeof(yes_type));
-};
-
-
-// instantiation of the static member
-template <typename MatrixType>
-const bool IsBlockMatrix<MatrixType>::value;
-
-
-template <typename MatrixType>
-inline void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix) const
-{
- // create a dummy and hand on to the function actually implementing this
- // feature in the cm.templates.h file.
- Vector<typename MatrixType::value_type> dummy(0);
- distribute_local_to_global(
- local_matrix,
- dummy,
- local_dof_indices,
- global_matrix,
- dummy,
- false,
- std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
-}
-
-
-
-template <typename MatrixType, typename VectorType>
-inline void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs) const
-{
- // enter the internal function with the respective block information set,
- // 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,
- std::integral_constant<bool, IsBlockMatrix<MatrixType>::value>());
-}
-
-
-
-template <typename SparsityPatternType>
-inline void
-ConstraintMatrix::add_entries_local_to_global(
- const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> & dof_mask) const
-{
- // enter the internal function with the respective block information set,
- // the actual implementation follows in the cm.templates.h file.
- add_entries_local_to_global(
- local_dof_indices,
- sparsity_pattern,
- keep_constrained_entries,
- dof_mask,
- std::integral_constant<bool, IsBlockMatrix<SparsityPatternType>::value>());
-}
-
+using ConstraintMatrix = AffineConstraints<double>;
DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 1999 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-
-#ifndef dealii_constraint_matrix_templates_h
-#define dealii_constraint_matrix_templates_h
-
-
-#include <deal.II/base/table.h>
-#include <deal.II/base/thread_local_storage.h>
-
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_sparsity_pattern.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/la_vector.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/sparsity_pattern.h>
-#include <deal.II/lac/trilinos_vector.h>
-
-#include <complex>
-#include <iomanip>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-template <typename number>
-void
-ConstraintMatrix::condense(SparseMatrix<number> &uncondensed) const
-{
- Vector<number> dummy(0);
- condense(uncondensed, dummy);
-}
-
-
-
-template <typename number>
-void
-ConstraintMatrix::condense(BlockSparseMatrix<number> &uncondensed) const
-{
- BlockVector<number> dummy(0);
- condense(uncondensed, dummy);
-}
-
-
-
-template <class VectorType>
-void
-ConstraintMatrix::condense(const VectorType &vec_ghosted, VectorType &vec) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
-
- // if this is called with different arguments, we need to copy the data over:
- if (&vec != &vec_ghosted)
- vec = vec_ghosted;
-
- // distribute all entries, and set them to zero. do so in
- // two loops because in the first one we need to add to elements
- // and in the second one we need to set elements to zero. for
- // parallel vectors, this can only work if we can put a compress()
- // in between, but we don't want to call compress() twice per entry
- for (std::vector<ConstraintLine>::const_iterator constraint_line =
- lines.begin();
- constraint_line != lines.end();
- ++constraint_line)
- {
- // in case the constraint is
- // inhomogeneous, this function is not
- // appropriate. Throw an exception.
- Assert(constraint_line->inhomogeneity == 0.,
- ExcMessage("Inhomogeneous constraint cannot be condensed "
- "without any matrix specified."));
-
- const typename VectorType::value_type old_value =
- vec_ghosted(constraint_line->index);
- for (size_type q = 0; q != constraint_line->entries.size(); ++q)
- if (vec.in_local_range(constraint_line->entries[q].first) == true)
- vec(constraint_line->entries[q].first) +=
- (static_cast<typename VectorType::value_type>(old_value) *
- constraint_line->entries[q].second);
- }
-
- vec.compress(VectorOperation::add);
-
- for (std::vector<ConstraintLine>::const_iterator constraint_line =
- lines.begin();
- constraint_line != lines.end();
- ++constraint_line)
- if (vec.in_local_range(constraint_line->index) == true)
- vec(constraint_line->index) = 0.;
-
- vec.compress(VectorOperation::insert);
-}
-
-
-
-template <class VectorType>
-void
-ConstraintMatrix::condense(VectorType &vec) const
-{
- condense(vec, vec);
-}
-
-
-
-template <typename number, class VectorType>
-void
-ConstraintMatrix::condense(SparseMatrix<number> &uncondensed,
- VectorType & vec) const
-{
- // check whether we work on real vectors
- // or we just used a dummy when calling
- // the other function above.
- const bool use_vectors = vec.size() == 0 ? false : true;
-
- const SparsityPattern &sparsity = uncondensed.get_sparsity_pattern();
-
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.is_compressed() == true, ExcMatrixNotClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
- if (use_vectors == true)
- AssertDimension(vec.size(), sparsity.n_rows());
-
- double average_diagonal = 0;
- for (size_type i = 0; i < uncondensed.m(); ++i)
- average_diagonal += std::abs(uncondensed.diag_element(i));
- average_diagonal /= uncondensed.m();
-
- // store for each index whether it must be
- // distributed or not. If entry is
- // invalid_size_type, no distribution is
- // necessary. otherwise, the number states
- // which line in the constraint matrix
- // handles this index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = c;
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- if (distribute[row] == numbers::invalid_size_type)
- // regular line. loop over cols
- {
- for (typename SparseMatrix<number>::iterator entry =
- uncondensed.begin(row);
- entry != uncondensed.end(row);
- ++entry)
- {
- const size_type column = entry->column();
-
- // end of row reached?
- // this should not
- // happen, since we only
- // operate on compressed
- // matrices!
- Assert(column != SparsityPattern::invalid_entry,
- ExcMatrixNotClosed());
-
- if (distribute[column] != numbers::invalid_size_type)
- // distribute entry at
- // regular row @p row
- // and irregular column
- // sparsity.get_column_numbers()[j];
- // set old entry to
- // zero
- {
- for (size_type q = 0;
- q != lines[distribute[column]].entries.size();
- ++q)
- {
- // need a temporary variable to avoid errors like
- // no known conversion from 'complex<typename
- // ProductType<float, double>::type>' to 'const
- // complex<float>' for 3rd argument
- number v = static_cast<number>(entry->value());
- v *= lines[distribute[column]].entries[q].second;
- uncondensed.add(
- row, lines[distribute[column]].entries[q].first, v);
- }
-
- // need to subtract this element from the
- // vector. this corresponds to an
- // explicit elimination in the respective
- // row of the inhomogeneous constraint in
- // the matrix with Gauss elimination
- if (use_vectors == true)
- vec(row) -= static_cast<number>(entry->value()) *
- lines[distribute[column]].inhomogeneity;
-
- // set old value to zero
- entry->value() = 0.;
- }
- }
- }
- else
- // row must be distributed
- {
- for (typename SparseMatrix<number>::iterator entry =
- uncondensed.begin(row);
- entry != uncondensed.end(row);
- ++entry)
- {
- const size_type column = entry->column();
-
- // end of row reached?
- // this should not
- // happen, since we only
- // operate on compressed
- // matrices!
- Assert(column != SparsityPattern::invalid_entry,
- ExcMatrixNotClosed());
-
- if (distribute[column] == numbers::invalid_size_type)
- // distribute entry at
- // irregular row
- // @p row and regular
- // column
- // column. set
- // old entry to zero
- {
- for (size_type q = 0;
- q != lines[distribute[row]].entries.size();
- ++q)
- {
- // need a temporary variable to avoid errors like
- // no known conversion from 'complex<typename
- // ProductType<float, double>::type>' to 'const
- // complex<float>' for 3rd argument
- number v = static_cast<number>(entry->value());
- v *= lines[distribute[row]].entries[q].second;
- uncondensed.add(
- lines[distribute[row]].entries[q].first, column, v);
- }
-
- // set old entry to zero
- entry->value() = 0.;
- }
- else
- // distribute entry at
- // irregular row @p row and
- // irregular column
- // @p column set old entry
- // to one on main
- // diagonal, zero otherwise
- {
- for (size_type p = 0;
- p != lines[distribute[row]].entries.size();
- ++p)
- {
- for (size_type q = 0;
- q != lines[distribute[column]].entries.size();
- ++q)
- {
- // need a temporary variable to avoid errors like
- // no known conversion from 'complex<typename
- // ProductType<float, double>::type>' to 'const
- // complex<float>' for 3rd argument
- number v = static_cast<number>(entry->value());
- v *= lines[distribute[row]].entries[p].second *
- lines[distribute[column]].entries[q].second;
- uncondensed.add(
- lines[distribute[row]].entries[p].first,
- lines[distribute[column]].entries[q].first,
- v);
- }
-
- if (use_vectors == true)
- vec(lines[distribute[row]].entries[p].first) -=
- static_cast<number>(entry->value()) *
- lines[distribute[row]].entries[p].second *
- lines[distribute[column]].inhomogeneity;
- }
-
- // set old entry to correct
- // value
- entry->value() = (row == column ? average_diagonal : 0.);
- }
- }
-
- // take care of vector
- if (use_vectors == true)
- {
- for (size_type q = 0; q != lines[distribute[row]].entries.size();
- ++q)
- vec(lines[distribute[row]].entries[q].first) +=
- (vec(row) * lines[distribute[row]].entries[q].second);
-
- vec(lines[distribute[row]].index) = 0.;
- }
- }
- }
-}
-
-
-
-template <typename number, class BlockVectorType>
-void
-ConstraintMatrix::condense(BlockSparseMatrix<number> &uncondensed,
- BlockVectorType & vec) const
-{
- // check whether we work on real vectors
- // or we just used a dummy when calling
- // the other function above.
- const bool use_vectors = vec.n_blocks() == 0 ? false : true;
-
- const size_type blocks = uncondensed.n_block_rows();
-
- const BlockSparsityPattern &sparsity = uncondensed.get_sparsity_pattern();
-
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.is_compressed() == true, ExcMatrixNotClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
- Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
- Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
- Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
- ExcNotQuadratic());
-
- if (use_vectors == true)
- {
- AssertDimension(vec.size(), sparsity.n_rows());
- AssertDimension(vec.n_blocks(), sparsity.n_block_rows());
- }
-
- double average_diagonal = 0;
- for (size_type b = 0; b < uncondensed.n_block_rows(); ++b)
- for (size_type i = 0; i < uncondensed.block(b, b).m(); ++i)
- average_diagonal += std::fabs(uncondensed.block(b, b).diag_element(i));
- average_diagonal /= uncondensed.m();
-
- const BlockIndices &index_mapping = sparsity.get_column_indices();
-
- // store for each index whether it must be
- // distributed or not. If entry is
- // numbers::invalid_size_type,
- // no distribution is necessary.
- // otherwise, the number states which line
- // in the constraint matrix handles this
- // index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = c;
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- // get index of this row
- // within the blocks
- const std::pair<size_type, size_type> block_index =
- index_mapping.global_to_local(row);
- const size_type block_row = block_index.first;
-
- if (distribute[row] == numbers::invalid_size_type)
- // regular line. loop over
- // all columns and see
- // whether this column must
- // be distributed
- {
- // to loop over all entries
- // in this row, we have to
- // loop over all blocks in
- // this blockrow and the
- // corresponding row
- // therein
- for (size_type block_col = 0; block_col < blocks; ++block_col)
- {
- for (typename SparseMatrix<number>::iterator entry =
- uncondensed.block(block_row, block_col)
- .begin(block_index.second);
- entry != uncondensed.block(block_row, block_col)
- .end(block_index.second);
- ++entry)
- {
- const size_type global_col =
- index_mapping.local_to_global(block_col, entry->column());
-
- if (distribute[global_col] != numbers::invalid_size_type)
- // distribute entry at
- // regular row @p row
- // and irregular column
- // global_col; set old
- // entry to zero
- {
- const double old_value = entry->value();
-
- for (size_type q = 0;
- q != lines[distribute[global_col]].entries.size();
- ++q)
- uncondensed.add(
- row,
- lines[distribute[global_col]].entries[q].first,
- old_value *
- lines[distribute[global_col]].entries[q].second);
-
- // need to subtract this element from the
- // vector. this corresponds to an
- // explicit elimination in the respective
- // row of the inhomogeneous constraint in
- // the matrix with Gauss elimination
- if (use_vectors == true)
- vec(row) -= entry->value() *
- lines[distribute[global_col]].inhomogeneity;
-
- entry->value() = 0.;
- }
- }
- }
- }
- else
- {
- // row must be
- // distributed. split the
- // whole row into the
- // chunks defined by the
- // blocks
- for (size_type block_col = 0; block_col < blocks; ++block_col)
- {
- for (typename SparseMatrix<number>::iterator entry =
- uncondensed.block(block_row, block_col)
- .begin(block_index.second);
- entry != uncondensed.block(block_row, block_col)
- .end(block_index.second);
- ++entry)
- {
- const size_type global_col =
- index_mapping.local_to_global(block_col, entry->column());
-
- if (distribute[global_col] == numbers::invalid_size_type)
- // distribute
- // entry at
- // irregular
- // row @p row
- // and regular
- // column
- // global_col. set
- // old entry to
- // zero
- {
- const double old_value = entry->value();
-
- for (size_type q = 0;
- q != lines[distribute[row]].entries.size();
- ++q)
- uncondensed.add(
- lines[distribute[row]].entries[q].first,
- global_col,
- old_value * lines[distribute[row]].entries[q].second);
-
- entry->value() = 0.;
- }
- else
- // distribute entry at
- // irregular row @p row
- // and irregular column
- // @p global_col set old
- // entry to one if on
- // main diagonal, zero
- // otherwise
- {
- const double old_value = entry->value();
-
- for (size_type p = 0;
- p != lines[distribute[row]].entries.size();
- ++p)
- {
- for (size_type q = 0;
- q !=
- lines[distribute[global_col]].entries.size();
- ++q)
- uncondensed.add(
- lines[distribute[row]].entries[p].first,
- lines[distribute[global_col]].entries[q].first,
- old_value *
- lines[distribute[row]].entries[p].second *
- lines[distribute[global_col]]
- .entries[q]
- .second);
-
- if (use_vectors == true)
- vec(lines[distribute[row]].entries[p].first) -=
- old_value *
- lines[distribute[row]].entries[p].second *
- lines[distribute[global_col]].inhomogeneity;
- }
-
- entry->value() =
- (row == global_col ? average_diagonal : 0.);
- }
- }
- }
-
- // take care of vector
- if (use_vectors == true)
- {
- for (size_type q = 0; q != lines[distribute[row]].entries.size();
- ++q)
- vec(lines[distribute[row]].entries[q].first) +=
- (vec(row) * lines[distribute[row]].entries[q].second);
-
- vec(lines[distribute[row]].index) = 0.;
- }
- }
- }
-}
-
-
-// TODO: I'm sure the following could be made more elegant by using a bit of
-// introspection using static member variables of the various vector
-// classes to dispatch between the different functions, rather than using
-// knowledge of the individual types
-
-// number of functions to select the right implementation for set_zero().
-namespace internal
-{
- namespace ConstraintMatrixImplementation
- {
- namespace
- {
- typedef types::global_dof_index size_type;
-
- template <class VectorType>
- void
- set_zero_parallel(const std::vector<size_type> &cm,
- VectorType & vec,
- size_type shift = 0)
- {
- Assert(!vec.has_ghost_elements(), ExcInternalError());
- IndexSet locally_owned = vec.locally_owned_elements();
- for (typename std::vector<size_type>::const_iterator it = cm.begin();
- it != cm.end();
- ++it)
- {
- // If shift>0 then we are working on a part of a BlockVector
- // so vec(i) is actually the global entry i+shift.
- // We first make sure the line falls into the range of vec,
- // then check if is part of the local part of the vector, before
- // finally setting it to 0.
- if ((*it) < shift)
- continue;
- size_type idx = *it - shift;
- if (idx < vec.size() && locally_owned.is_element(idx))
- internal::ElementAccess<VectorType>::set(0., idx, vec);
- }
- }
-
- template <typename Number>
- void
- set_zero_parallel(const std::vector<size_type> & cm,
- LinearAlgebra::distributed::Vector<Number> &vec,
- size_type shift = 0)
- {
- for (typename std::vector<size_type>::const_iterator it = cm.begin();
- it != cm.end();
- ++it)
- {
- // If shift>0 then we are working on a part of a BlockVector
- // so vec(i) is actually the global entry i+shift.
- // We first make sure the line falls into the range of vec,
- // then check if is part of the local part of the vector, before
- // finally setting it to 0.
- if ((*it) < shift)
- continue;
- size_type idx = *it - shift;
- if (vec.in_local_range(idx))
- vec(idx) = 0.;
- }
- vec.zero_out_ghosts();
- }
-
- template <class VectorType>
- void
- set_zero_in_parallel(const std::vector<size_type> &cm,
- VectorType & vec,
- std::integral_constant<bool, false>)
- {
- set_zero_parallel(cm, vec, 0);
- }
-
- // in parallel for BlockVectors
- template <class VectorType>
- void
- set_zero_in_parallel(const std::vector<size_type> &cm,
- VectorType & vec,
- std::integral_constant<bool, true>)
- {
- size_type start_shift = 0;
- for (size_type j = 0; j < vec.n_blocks(); ++j)
- {
- set_zero_parallel(cm, vec.block(j), start_shift);
- start_shift += vec.block(j).size();
- }
- }
-
- template <class VectorType>
- void
- set_zero_serial(const std::vector<size_type> &cm, VectorType &vec)
- {
- for (typename std::vector<size_type>::const_iterator it = cm.begin();
- it != cm.end();
- ++it)
- vec(*it) = 0.;
- }
-
- template <class VectorType>
- void
- set_zero_all(const std::vector<size_type> &cm, VectorType &vec)
- {
- set_zero_in_parallel<VectorType>(
- cm,
- vec,
- std::integral_constant<bool, IsBlockVector<VectorType>::value>());
- vec.compress(VectorOperation::insert);
- }
-
-
- template <class T>
- void
- set_zero_all(const std::vector<size_type> &cm, dealii::Vector<T> &vec)
- {
- set_zero_serial(cm, vec);
- }
-
- template <class T>
- void
- set_zero_all(const std::vector<size_type> &cm,
- dealii::BlockVector<T> & vec)
- {
- set_zero_serial(cm, vec);
- }
- } // namespace
- } // namespace ConstraintMatrixImplementation
-} // namespace internal
-
-
-template <class VectorType>
-void
-ConstraintMatrix::set_zero(VectorType &vec) const
-{
- // since we lines is a private member, we cannot pass it to the functions
- // above. therefore, copy the content which is cheap
- std::vector<size_type> constrained_lines(lines.size());
- for (unsigned int i = 0; i < lines.size(); ++i)
- constrained_lines[i] = lines[i].index;
- internal::ConstraintMatrixImplementation::set_zero_all(constrained_lines,
- vec);
-}
-
-
-
-template <typename VectorType, typename LocalType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const Vector<LocalType> & local_vector,
- const std::vector<size_type> &local_dof_indices,
- VectorType & global_vector,
- const FullMatrix<LocalType> & local_matrix) const
-{
- distribute_local_to_global(local_vector,
- local_dof_indices,
- local_dof_indices,
- global_vector,
- local_matrix,
- true);
-}
-
-
-
-template <typename VectorType, typename LocalType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const Vector<LocalType> & local_vector,
- const std::vector<size_type> &local_dof_indices_row,
- const std::vector<size_type> &local_dof_indices_col,
- VectorType & global_vector,
- const FullMatrix<LocalType> & local_matrix,
- bool diagonal) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
- AssertDimension(local_vector.size(), local_dof_indices_row.size());
- AssertDimension(local_matrix.m(), local_dof_indices_row.size());
- AssertDimension(local_matrix.n(), local_dof_indices_col.size());
-
- // diagonal checks if we have only one index set (if both are equal
- // diagonal should be set to true).
- // If true we do both, assembly of the right hand side (next lines)
- // and (see further below) modifications of the right hand side
- // according to the inhomogeneous constraints.
- // Otherwise we only modify the right hand side according to
- // local_matrix and the inhomogeneous constraints, and omit the vector add.
-
- const size_type m_local_dofs = local_dof_indices_row.size();
- const size_type n_local_dofs = local_dof_indices_col.size();
- if (lines.empty())
- {
- if (diagonal)
- global_vector.add(local_dof_indices_row, local_vector);
- }
- else
- for (size_type i = 0; i < n_local_dofs; ++i)
- {
- // check whether the current index is
- // constrained. if not, just write the entry
- // into the vector. otherwise, need to resolve
- // the constraint
- if (is_constrained(local_dof_indices_col[i]) == false)
- {
- if (diagonal)
- global_vector(local_dof_indices_row[i]) += local_vector(i);
- continue;
- }
-
- // find the constraint line to the given
- // global dof index
- const size_type line_index =
- calculate_line_index(local_dof_indices_col[i]);
- const ConstraintLine *position = lines_cache.size() <= line_index ?
- nullptr :
- &lines[lines_cache[line_index]];
-
- // Gauss elimination of the matrix columns with the inhomogeneity.
- // Go through them one by one and again check whether they are
- // constrained. If so, distribute the constraint
- const double val = position->inhomogeneity;
- if (val != 0)
- for (size_type j = 0; j < m_local_dofs; ++j)
- {
- if (is_constrained(local_dof_indices_row[j]) == false)
- {
- global_vector(local_dof_indices_row[j]) -=
- val * local_matrix(j, i);
- continue;
- }
-
- const LocalType matrix_entry = local_matrix(j, i);
-
- if (matrix_entry == LocalType())
- continue;
-
- const ConstraintLine &position_j =
- lines[lines_cache[calculate_line_index(
- local_dof_indices_row[j])]];
-
- for (size_type q = 0; q < position_j.entries.size(); ++q)
- {
- Assert(
- !(!local_lines.size() ||
- local_lines.is_element(position_j.entries[q].first)) ||
- is_constrained(position_j.entries[q].first) == false,
- ExcMessage("Tried to distribute to a fixed dof."));
- global_vector(position_j.entries[q].first) -=
- val * position_j.entries[q].second * matrix_entry;
- }
- }
-
- // now distribute the constraint,
- // but make sure we don't touch
- // the entries of fixed dofs
- if (diagonal)
- {
- for (size_type j = 0; j < position->entries.size(); ++j)
- {
- Assert(!(!local_lines.size() ||
- local_lines.is_element(position->entries[j].first)) ||
- is_constrained(position->entries[j].first) == false,
- ExcMessage("Tried to distribute to a fixed dof."));
- global_vector(position->entries[j].first) +=
- local_vector(i) * position->entries[j].second;
- }
- }
- }
-}
-
-
-
-namespace internal
-{
- namespace
- {
- // create an output vector that consists of the input vector's locally owned
- // elements plus some ghost elements that need to be imported from elsewhere
- //
- // this is an operation that is different for all vector types and so we
- // need a few overloads
-#ifdef DEAL_II_WITH_TRILINOS
- void
- import_vector_with_ghost_elements(
- const TrilinosWrappers::MPI::Vector &vec,
- const IndexSet & /*locally_owned_elements*/,
- const IndexSet & needed_elements,
- TrilinosWrappers::MPI::Vector &output,
- const std::integral_constant<bool, false> /*is_block_vector*/)
- {
- Assert(!vec.has_ghost_elements(), ExcGhostsPresent());
-# ifdef DEAL_II_WITH_MPI
- const Epetra_MpiComm *mpi_comm =
- dynamic_cast<const Epetra_MpiComm *>(&vec.trilinos_vector().Comm());
-
- Assert(mpi_comm != nullptr, ExcInternalError());
- output.reinit(needed_elements, mpi_comm->GetMpiComm());
-# else
- output.reinit(needed_elements, MPI_COMM_SELF);
-# endif
- output = vec;
- }
-#endif
-
-#ifdef DEAL_II_WITH_PETSC
- void
- import_vector_with_ghost_elements(
- const PETScWrappers::MPI::Vector &vec,
- const IndexSet & locally_owned_elements,
- const IndexSet & needed_elements,
- PETScWrappers::MPI::Vector & output,
- const std::integral_constant<bool, false> /*is_block_vector*/)
- {
- output.reinit(
- locally_owned_elements, needed_elements, vec.get_mpi_communicator());
- output = vec;
- }
-#endif
-
- template <typename number>
- void
- import_vector_with_ghost_elements(
- const LinearAlgebra::distributed::Vector<number> &vec,
- const IndexSet & locally_owned_elements,
- const IndexSet & needed_elements,
- LinearAlgebra::distributed::Vector<number> & output,
- const std::integral_constant<bool, false> /*is_block_vector*/)
- {
- // TODO: the in vector might already have all elements. need to find a
- // way to efficiently avoid the copy then
- const_cast<LinearAlgebra::distributed::Vector<number> &>(vec)
- .zero_out_ghosts();
- output.reinit(
- locally_owned_elements, needed_elements, vec.get_mpi_communicator());
- output = vec;
- output.update_ghost_values();
- }
-
-
- // all other vector non-block vector types are sequential and we should
- // not have this function called at all -- so throw an exception
- template <typename Vector>
- void
- import_vector_with_ghost_elements(
- const Vector & /*vec*/,
- const IndexSet & /*locally_owned_elements*/,
- const IndexSet & /*needed_elements*/,
- Vector & /*output*/,
- const std::integral_constant<bool, false> /*is_block_vector*/)
- {
- Assert(false, ExcMessage("We shouldn't even get here!"));
- }
-
-
- // for block vectors, simply dispatch to the individual blocks
- template <class VectorType>
- void
- import_vector_with_ghost_elements(
- const VectorType &vec,
- const IndexSet & locally_owned_elements,
- const IndexSet & needed_elements,
- VectorType & output,
- const std::integral_constant<bool, true> /*is_block_vector*/)
- {
- output.reinit(vec.n_blocks());
-
- types::global_dof_index block_start = 0;
- for (unsigned int b = 0; b < vec.n_blocks(); ++b)
- {
- import_vector_with_ghost_elements(
- vec.block(b),
- locally_owned_elements.get_view(block_start,
- block_start + vec.block(b).size()),
- needed_elements.get_view(block_start,
- block_start + vec.block(b).size()),
- output.block(b),
- std::integral_constant<bool, false>());
- block_start += vec.block(b).size();
- }
-
- output.collect_sizes();
- }
- } // namespace
-} // namespace internal
-
-
-template <class VectorType>
-void
-ConstraintMatrix::distribute(VectorType &vec) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
-
- // if the vector type supports parallel storage and if the vector actually
- // does store only part of the vector, distributing is slightly more
- // complicated. we might be able to skip the complicated part if one
- // processor owns everything and pretend that this is a sequential vector,
- // but it is difficult for the other processors to know whether they should
- // not do anything or if other processors will create a temporary vector,
- // exchange data (requiring communication, maybe even with the processors
- // that do not own anything because of that particular parallel model), and
- // call compress() finally. the first case here is for the complicated case,
- // the last else is for the simple case (sequential vector)
- const IndexSet vec_owned_elements = vec.locally_owned_elements();
-
- if (dealii::is_serial_vector<VectorType>::value == false)
- {
- // This processor owns only part of the vector. one may think that
- // every processor should be able to simply communicate those elements
- // it owns and for which it knows that they act as sources to constrained
- // DoFs to the owner of these DoFs. This would lead to a scheme where all
- // we need to do is to add some local elements to (possibly non-local)
- // ones and then call compress().
- //
- // Alas, this scheme does not work as evidenced by the disaster of bug
- // #51, see http://code.google.com/p/dealii/issues/detail?id=51 and the
- // reversion of one attempt that implements this in r29662. Rather, we
- // need to get a vector that has all the *sources* or constraints we
- // own locally, possibly as ghost vector elements, then read from them,
- // and finally throw away the ghosted vector. Implement this in the
- // following.
- IndexSet needed_elements = vec_owned_elements;
-
- typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
- for (constraint_iterator it = lines.begin(); it != lines.end(); ++it)
- if (vec_owned_elements.is_element(it->index))
- for (unsigned int i = 0; i < it->entries.size(); ++i)
- if (!vec_owned_elements.is_element(it->entries[i].first))
- needed_elements.add_index(it->entries[i].first);
-
- VectorType ghosted_vector;
- internal::import_vector_with_ghost_elements(
- vec,
- vec_owned_elements,
- needed_elements,
- ghosted_vector,
- std::integral_constant<bool, IsBlockVector<VectorType>::value>());
-
- for (constraint_iterator it = lines.begin(); it != lines.end(); ++it)
- if (vec_owned_elements.is_element(it->index))
- {
- typename VectorType::value_type new_value = it->inhomogeneity;
- for (unsigned int i = 0; i < it->entries.size(); ++i)
- new_value += (static_cast<typename VectorType::value_type>(
- internal::ElementAccess<VectorType>::get(
- ghosted_vector, it->entries[i].first)) *
- it->entries[i].second);
- AssertIsFinite(new_value);
- internal::ElementAccess<VectorType>::set(new_value, it->index, vec);
- }
-
- // now compress to communicate the entries that we added to
- // and that weren't to local processors to the owner
- //
- // this shouldn't be strictly necessary but it probably doesn't
- // hurt either
- vec.compress(VectorOperation::insert);
- }
- else
- // purely sequential vector (either because the type doesn't
- // support anything else or because it's completely stored
- // locally)
- {
- std::vector<ConstraintLine>::const_iterator next_constraint =
- lines.begin();
- for (; next_constraint != lines.end(); ++next_constraint)
- {
- // fill entry in line
- // next_constraint.index by adding the
- // different contributions
- typename VectorType::value_type new_value =
- next_constraint->inhomogeneity;
- for (unsigned int i = 0; i < next_constraint->entries.size(); ++i)
- new_value += (static_cast<typename VectorType::value_type>(
- internal::ElementAccess<VectorType>::get(
- vec, next_constraint->entries[i].first)) *
- next_constraint->entries[i].second);
- AssertIsFinite(new_value);
- internal::ElementAccess<VectorType>::set(
- new_value, next_constraint->index, vec);
- }
- }
-}
-
-
-
-// Some helper definitions for the local_to_global functions.
-namespace internals
-{
- typedef types::global_dof_index size_type;
-
- // this struct contains all the information we need to store about each of
- // the global entries (global_row): are they obtained directly by some local
- // entry (local_row) or some constraints (constraint_position). This is not
- // directly used in the user code, but accessed via the GlobalRowsFromLocal.
- //
- // The actions performed here correspond to reshaping the constraint
- // information from global degrees of freedom to local ones (i.e.,
- // cell-related DoFs), and also transforming the constraint information from
- // compressed row storage (each local dof that is constrained has a list of
- // constraint entries associated to it) into compressed column storage based
- // on the cell-related DoFs (we have a list of global degrees of freedom,
- // and to each we have a list of local rows where the entries come from). To
- // increase the speed, we additionally store whether an entry is generated
- // directly from the local degrees of freedom or whether it comes from a
- // constraint.
- struct Distributing
- {
- Distributing(const size_type global_row = numbers::invalid_size_type,
- const size_type local_row = numbers::invalid_size_type);
- Distributing(const Distributing &in);
- Distributing &
- operator=(const Distributing &in);
- bool
- operator<(const Distributing &in) const
- {
- return global_row < in.global_row;
- }
-
- size_type global_row;
- size_type local_row;
- mutable size_type constraint_position;
- };
-
- inline Distributing::Distributing(const size_type global_row,
- const size_type local_row) :
- global_row(global_row),
- local_row(local_row),
- constraint_position(numbers::invalid_size_type)
- {}
-
- inline Distributing::Distributing(const Distributing &in) :
- constraint_position(numbers::invalid_size_type)
- {
- *this = (in);
- }
-
- inline Distributing &
- Distributing::operator=(const Distributing &in)
- {
- global_row = in.global_row;
- local_row = in.local_row;
- // the constraints pointer should not contain any data here.
- Assert(constraint_position == numbers::invalid_size_type,
- ExcInternalError());
-
- if (in.constraint_position != numbers::invalid_size_type)
- {
- constraint_position = in.constraint_position;
- in.constraint_position = numbers::invalid_size_type;
- }
- return *this;
- }
-
-
-
- // this is a cache for constraints that are encountered on a local level.
- // The functionality is similar to
- // std::vector<std::vector<std::pair<uint,double> > >, but tuned so that
- // frequent memory allocation for each entry is avoided. The data is put
- // into a std::vector<std::pair<uint,double> > and the row length is kept
- // fixed at row_length. Both the number of rows and the row length can
- // change is this structure is filled. In that case, the data is
- // rearranged. This is not directly used in the user code, but accessed via
- // the GlobalRowsFromLocal.
- struct DataCache
- {
- DataCache() : row_length(8)
- {}
-
- void
- reinit()
- {
- individual_size.resize(0);
- data.resize(0);
- }
-
- size_type
- insert_new_index(const std::pair<size_type, double> &pair)
- {
- Assert(row_length > 0, ExcInternalError());
- const unsigned int index = individual_size.size();
- individual_size.push_back(1);
- data.resize(individual_size.size() * row_length);
- data[index * row_length] = pair;
- individual_size[index] = 1;
- return index;
- }
-
- void
- append_index(const size_type index,
- const std::pair<size_type, double> &pair)
- {
- AssertIndexRange(index, individual_size.size());
- const size_type my_length = individual_size[index];
- if (my_length == row_length)
- {
- AssertDimension(data.size(), individual_size.size() * row_length);
- // no space left in this row, need to double row_length and
- // rearrange the data items. Move all items to the right except the
- // first one, starting at the back. Since individual_size contains
- // at least one element when we get here, subtracting 1 works fine.
- data.resize(2 * data.size());
- for (size_type i = individual_size.size() - 1; i > 0; --i)
- {
- const auto ptr = data.data();
- std::move_backward(ptr + i * row_length,
- ptr + i * row_length + individual_size[i],
- ptr + i * 2 * row_length + individual_size[i]);
- }
- row_length *= 2;
- }
- data[index * row_length + my_length] = pair;
- individual_size[index] = my_length + 1;
- }
-
- size_type
- get_size(const size_type index) const
- {
- return individual_size[index];
- }
-
- const std::pair<size_type, double> *
- get_entry(const size_type index) const
- {
- return &data[index * row_length];
- }
-
- size_type row_length;
-
- std::vector<std::pair<size_type, double>> data;
-
- std::vector<size_type> individual_size;
- };
-
-
-
- // collects all the global rows from a local contribution (cell) and their
- // origin (direct/constraint). this is basically a vector consisting of
- // "Distributing" structs using access via the DataCache. Provides some
- // specialized sort and insert functions.
- //
- // in case there are no constraints, this is basically a list of pairs
- // <uint,unit> with the first index being the global index and the second
- // index the local index. The list is sorted with respect to the global
- // index.
- //
- // in case there are constraints, a global dof might get a contribution also
- // because it gets data from a constrained dof. This means that a global dof
- // might also have indirect contributions from a local dof via a constraint,
- // besides the direct ones.
- //
- // The actions performed here correspond to reshaping the constraint
- // information from global degrees of freedom to local ones (i.e.,
- // cell-related DoFs), and also transforming the constraint information from
- // compressed row storage (each local dof that is constrained has a list of
- // constraint entries associated to it) into compressed column storage based
- // on the cell-related DoFs (we have a list of global degrees of freedom,
- // and to each we have a list of local rows where the entries come from). To
- // increase the speed, we additionally store whether an entry is generated
- // directly from the local degrees of freedom or whether it comes from a
- // constraint.
- class GlobalRowsFromLocal
- {
- public:
- GlobalRowsFromLocal() : n_active_rows(0), n_inhomogeneous_rows(0)
- {}
-
- void
- reinit(const size_type n_local_rows)
- {
- total_row_indices.resize(n_local_rows);
- for (unsigned int i = 0; i < n_local_rows; ++i)
- total_row_indices[i].constraint_position = numbers::invalid_size_type;
- n_active_rows = n_local_rows;
- n_inhomogeneous_rows = 0;
- data_cache.reinit();
- }
-
- // implemented below
- void
- insert_index(const size_type global_row,
- const size_type local_row,
- const double constraint_value);
- void
- sort();
-
- // Print object for debugging purpose
- void
- print(std::ostream &os)
- {
- os << "Active rows " << n_active_rows << std::endl
- << "Constr rows " << n_constraints() << std::endl
- << "Inhom rows " << n_inhomogeneous_rows << std::endl
- << "Local: ";
- for (size_type i = 0; i < total_row_indices.size(); ++i)
- os << ' ' << std::setw(4) << total_row_indices[i].local_row;
- os << std::endl << "Global:";
- for (size_type i = 0; i < total_row_indices.size(); ++i)
- os << ' ' << std::setw(4) << total_row_indices[i].global_row;
- os << std::endl << "ConPos:";
- for (size_type i = 0; i < total_row_indices.size(); ++i)
- os << ' ' << std::setw(4) << total_row_indices[i].constraint_position;
- os << std::endl;
- }
-
-
- // return all kind of information on the constraints
-
- // returns the number of global indices in the struct
- size_type
- size() const
- {
- return n_active_rows;
- }
-
- // returns the number of constraints that are associated to the
- // counter_index-th entry in the list
- size_type
- size(const size_type counter_index) const
- {
- return (total_row_indices[counter_index].constraint_position ==
- numbers::invalid_size_type ?
- 0 :
- data_cache.get_size(
- total_row_indices[counter_index].constraint_position));
- }
-
- // returns the global row of the counter_index-th entry in the list
- size_type
- global_row(const size_type counter_index) const
- {
- return total_row_indices[counter_index].global_row;
- }
-
- // returns the global row of the counter_index-th entry in the list
- size_type &
- global_row(const size_type counter_index)
- {
- return total_row_indices[counter_index].global_row;
- }
-
- // returns the local row in the cell matrix associated with the
- // counter_index-th entry in the list. Returns invalid_size_type for
- // constrained rows
- size_type
- local_row(const size_type counter_index) const
- {
- return total_row_indices[counter_index].local_row;
- }
-
- // writable index
- size_type &
- local_row(const size_type counter_index)
- {
- return total_row_indices[counter_index].local_row;
- }
-
- // returns the local row in the cell matrix associated with the
- // counter_index-th entry in the list in the index_in_constraint-th
- // position of constraints
- size_type
- local_row(const size_type counter_index,
- const size_type index_in_constraint) const
- {
- return (data_cache.get_entry(total_row_indices[counter_index]
- .constraint_position)[index_in_constraint])
- .first;
- }
-
- // returns the value of the constraint in the counter_index-th entry in
- // the list in the index_in_constraint-th position of constraints
- double
- constraint_value(const size_type counter_index,
- const size_type index_in_constraint) const
- {
- return (data_cache.get_entry(total_row_indices[counter_index]
- .constraint_position)[index_in_constraint])
- .second;
- }
-
- // returns whether there is one row with indirect contributions (i.e.,
- // there has been at least one constraint with non-trivial ConstraintLine)
- bool
- have_indirect_rows() const
- {
- return data_cache.individual_size.empty() == false;
- }
-
- // append an entry that is constrained. This means that there is one less
- // nontrivial row
- void
- insert_constraint(const size_type constrained_local_dof)
- {
- --n_active_rows;
- total_row_indices[n_active_rows].local_row = constrained_local_dof;
- total_row_indices[n_active_rows].global_row = numbers::invalid_size_type;
- }
-
- // returns the number of constrained dofs in the structure. Constrained
- // dofs do not contribute directly to the matrix, but are needed in order
- // to set matrix diagonals and resolve inhomogeneities
- size_type
- n_constraints() const
- {
- return total_row_indices.size() - n_active_rows;
- }
-
- // returns the number of constrained dofs in the structure that have an
- // inhomogeneity
- size_type
- n_inhomogeneities() const
- {
- return n_inhomogeneous_rows;
- }
-
- // tells the structure that the ith constraint is
- // inhomogeneous. inhomogeneous constraints contribute to right hand
- // sides, so to have fast access to them, put them before homogeneous
- // constraints
- void
- set_ith_constraint_inhomogeneous(const size_type i)
- {
- Assert(i >= n_inhomogeneous_rows, ExcInternalError());
- std::swap(total_row_indices[n_active_rows + i],
- total_row_indices[n_active_rows + n_inhomogeneous_rows]);
- n_inhomogeneous_rows++;
- }
-
- // the local row where constraint number i was detected, to find that row
- // easily when the GlobalRowsToLocal has been set up
- size_type
- constraint_origin(size_type i) const
- {
- return total_row_indices[n_active_rows + i].local_row;
- }
-
- // a vector that contains all the global ids and the corresponding local
- // ids as well as a pointer to that data where we store how to resolve
- // constraints.
- std::vector<Distributing> total_row_indices;
-
- private:
- // holds the actual data from the constraints
- DataCache data_cache;
-
- // how many rows there are, constraints disregarded
- size_type n_active_rows;
-
- // the number of rows with inhomogeneous constraints
- size_type n_inhomogeneous_rows;
- };
-
- // a function that appends an additional row to the list of values, or
- // appends a value to an already existing row. Similar functionality as for
- // std::map<size_type,Distributing>, but here done for a
- // std::vector<Distributing>, much faster for short lists as we have them
- // here
- inline void
- GlobalRowsFromLocal::insert_index(const size_type global_row,
- const size_type local_row,
- const double constraint_value)
- {
- typedef std::vector<Distributing>::iterator index_iterator;
- index_iterator pos, pos1;
- Distributing row_value(global_row);
- std::pair<size_type, double> constraint(local_row, constraint_value);
-
- // check whether the list was really sorted before entering here
- for (size_type i = 1; i < n_active_rows; ++i)
- Assert(total_row_indices[i - 1] < total_row_indices[i],
- ExcInternalError());
-
- pos = Utilities::lower_bound(total_row_indices.begin(),
- total_row_indices.begin() + n_active_rows,
- row_value);
- if (pos->global_row == global_row)
- pos1 = pos;
- else
- {
- pos1 = total_row_indices.insert(pos, row_value);
- ++n_active_rows;
- }
-
- if (pos1->constraint_position == numbers::invalid_size_type)
- pos1->constraint_position = data_cache.insert_new_index(constraint);
- else
- data_cache.append_index(pos1->constraint_position, constraint);
- }
-
- // this sort algorithm sorts std::vector<Distributing>, but does not take
- // the constraints into account. this means that in case that constraints
- // are already inserted, this function does not work as expected. Use
- // shellsort, which is very fast in case the indices are already sorted
- // (which is the usual case with DG elements), and not too slow in other
- // cases
- inline void
- GlobalRowsFromLocal::sort()
- {
- size_type i, j, j2, temp, templ, istep;
- size_type step;
-
- // check whether the constraints are really empty.
- const size_type length = size();
-
- // make sure that we are in the range of the vector
- AssertIndexRange(length, total_row_indices.size() + 1);
- for (size_type i = 0; i < length; ++i)
- Assert(total_row_indices[i].constraint_position ==
- numbers::invalid_size_type,
- ExcInternalError());
-
- step = length / 2;
- while (step > 0)
- {
- for (i = step; i < length; i++)
- {
- istep = step;
- j = i;
- j2 = j - istep;
- temp = total_row_indices[i].global_row;
- templ = total_row_indices[i].local_row;
- if (total_row_indices[j2].global_row > temp)
- {
- while ((j >= istep) &&
- (total_row_indices[j2].global_row > temp))
- {
- total_row_indices[j].global_row =
- total_row_indices[j2].global_row;
- total_row_indices[j].local_row =
- total_row_indices[j2].local_row;
- j = j2;
- j2 -= istep;
- }
- total_row_indices[j].global_row = temp;
- total_row_indices[j].local_row = templ;
- }
- }
- step = step >> 1;
- }
- }
-
-
-
- /**
- * Scratch data that is used during calls to distribute_local_to_global and
- * add_entries_local_to_global. In order to avoid frequent memory
- * allocation, we keep the data alive from one call to the next in a static
- * variable. Since we want to allow for different number types in matrices,
- * this is a template.
- *
- * Since each thread gets its private version of scratch data out of the
- * ThreadLocalStorage, no conflicting access can occur. For this to be
- * valid, we need to make sure that no call within
- * distribute_local_to_global is made that by itself can spawn tasks.
- * Otherwise, we might end up in a situation where several threads fight for
- * the data.
- *
- * Access to the scratch data is only through the accessor class which
- * handles the access as well as marking the data as used.
- */
- template <typename MatrixScalar, typename VectorScalar = MatrixScalar>
- class ConstraintMatrixData
- {
- public:
- struct ScratchData
- {
- /**
- * Constructor, does nothing.
- */
- ScratchData() : in_use(false)
- {}
-
- /**
- * Copy constructor, does nothing
- */
- ScratchData(const ScratchData &) : in_use(false)
- {}
-
- /**
- * Stores whether the data is currently in use.
- */
- bool in_use;
-
- /**
- * Temporary array for column indices
- */
- std::vector<size_type> columns;
-
- /**
- * Temporary array for column values
- */
- std::vector<MatrixScalar> values;
-
- /**
- * Temporary array for block start indices
- */
- std::vector<size_type> block_starts;
-
- /**
- * Temporary array for vector indices
- */
- std::vector<size_type> vector_indices;
-
- /**
- * Temporary array for vector values
- */
- std::vector<VectorScalar> vector_values;
-
- /**
- * Data array for reorder row/column indices.
- */
- GlobalRowsFromLocal global_rows;
-
- /**
- * Data array for reorder row/column indices.
- */
- GlobalRowsFromLocal global_columns;
- };
-
- /**
- * Accessor class to guard access to scratch_data
- */
- class ScratchDataAccessor
- {
- public:
- /**
- * Constructor. Grabs a scratch data object on the current thread and
- * mark it as used
- */
- ScratchDataAccessor() :
- my_scratch_data(&ConstraintMatrixData::scratch_data.get())
- {
- Assert(my_scratch_data->in_use == false,
- ExcMessage(
- "Access to thread-local scratch data tried, but it is already "
- "in use"));
- my_scratch_data->in_use = true;
- }
-
- /**
- * Destructor. Mark scratch data as available again.
- */
- ~ScratchDataAccessor()
- {
- my_scratch_data->in_use = false;
- }
-
- /**
- * Dereferencing operator.
- */
- ScratchData &operator*()
- {
- return *my_scratch_data;
- }
-
- /**
- * Dereferencing operator.
- */
- ScratchData *operator->()
- {
- return my_scratch_data;
- }
-
- private:
- ScratchData *my_scratch_data;
- };
-
- private:
- /**
- * The actual data object that contains a scratch data for each thread.
- */
- static Threads::ThreadLocalStorage<ScratchData> scratch_data;
- };
-
-
-
- // function for block matrices: Find out where in the list of local dofs
- // (sorted according to global ids) the individual blocks start. Transform
- // the global indices to block-local indices in order to be able to use
- // functions like vector.block(1)(block_local_id), instead of
- // vector(global_id). This avoids transforming indices one-by-one later on.
- template <class BlockType>
- inline void
- make_block_starts(const BlockType & block_object,
- GlobalRowsFromLocal & global_rows,
- std::vector<size_type> &block_starts)
- {
- AssertDimension(block_starts.size(), block_object.n_block_rows() + 1);
-
- typedef std::vector<Distributing>::iterator row_iterator;
- row_iterator block_indices = global_rows.total_row_indices.begin();
-
- const size_type num_blocks = block_object.n_block_rows();
- const size_type n_active_rows = global_rows.size();
-
- // find end of rows.
- block_starts[0] = 0;
- for (size_type i = 1; i < num_blocks; ++i)
- {
- row_iterator first_block = Utilities::lower_bound(
- block_indices,
- global_rows.total_row_indices.begin() + n_active_rows,
- Distributing(block_object.get_row_indices().block_start(i)));
- block_starts[i] = first_block - global_rows.total_row_indices.begin();
- block_indices = first_block;
- }
- block_starts[num_blocks] = n_active_rows;
-
- // transform row indices to block-local index space
- for (size_type i = block_starts[1]; i < n_active_rows; ++i)
- global_rows.global_row(i) = block_object.get_row_indices()
- .global_to_local(global_rows.global_row(i))
- .second;
- }
-
-
-
- // same as before, but for std::vector<uint> instead of
- // GlobalRowsFromLocal. Used in functions for sparsity patterns.
- template <class BlockType>
- inline void
- make_block_starts(const BlockType & block_object,
- std::vector<size_type> &row_indices,
- std::vector<size_type> &block_starts)
- {
- AssertDimension(block_starts.size(), block_object.n_block_rows() + 1);
-
- typedef std::vector<size_type>::iterator row_iterator;
- row_iterator col_indices = row_indices.begin();
-
- const size_type num_blocks = block_object.n_block_rows();
-
- // find end of rows.
- block_starts[0] = 0;
- for (size_type i = 1; i < num_blocks; ++i)
- {
- row_iterator first_block =
- Utilities::lower_bound(col_indices,
- row_indices.end(),
- block_object.get_row_indices().block_start(i));
- block_starts[i] = first_block - row_indices.begin();
- col_indices = first_block;
- }
- block_starts[num_blocks] = row_indices.size();
-
- // transform row indices to local index space
- for (size_type i = block_starts[1]; i < row_indices.size(); ++i)
- row_indices[i] =
- block_object.get_row_indices().global_to_local(row_indices[i]).second;
- }
-
-
-
- // resolves constraints of one column at the innermost loop. goes through
- // the origin of each global entry and finds out which data we need to
- // collect.
- template <typename LocalType>
- static inline LocalType
- resolve_matrix_entry(const GlobalRowsFromLocal & global_rows,
- const GlobalRowsFromLocal & global_cols,
- const size_type i,
- const size_type j,
- const size_type loc_row,
- const FullMatrix<LocalType> &local_matrix)
- {
- const size_type loc_col = global_cols.local_row(j);
- LocalType col_val;
-
- // case 1: row has direct contribution in local matrix. decide whether col
- // has a direct contribution. if not, set the value to zero.
- if (loc_row != numbers::invalid_size_type)
- {
- col_val = ((loc_col != numbers::invalid_size_type) ?
- local_matrix(loc_row, loc_col) :
- 0);
-
- // account for indirect contributions by constraints in column
- for (size_type p = 0; p < global_cols.size(j); ++p)
- col_val += (local_matrix(loc_row, global_cols.local_row(j, p)) *
- global_cols.constraint_value(j, p));
- }
-
- // case 2: row has no direct contribution in local matrix
- else
- col_val = 0;
-
- // account for indirect contributions by constraints in row, going trough
- // the direct and indirect references in the given column.
- for (size_type q = 0; q < global_rows.size(i); ++q)
- {
- LocalType add_this =
- (loc_col != numbers::invalid_size_type) ?
- local_matrix(global_rows.local_row(i, q), loc_col) :
- 0;
-
- for (size_type p = 0; p < global_cols.size(j); ++p)
- add_this += (local_matrix(global_rows.local_row(i, q),
- global_cols.local_row(j, p)) *
- global_cols.constraint_value(j, p));
- col_val += add_this * global_rows.constraint_value(i, q);
- }
- return col_val;
- }
-
-
-
- // computes all entries that need to be written into global_rows[i]. Lists
- // the resulting values in val_ptr, and the corresponding column indices in
- // col_ptr.
- template <typename number, typename LocalType>
- inline void
- resolve_matrix_row(const GlobalRowsFromLocal & global_rows,
- const GlobalRowsFromLocal & global_cols,
- const size_type i,
- const size_type column_start,
- const size_type column_end,
- const FullMatrix<LocalType> &local_matrix,
- size_type *& col_ptr,
- number *& val_ptr)
- {
- if (column_end == column_start)
- return;
-
- AssertIndexRange(column_end - 1, global_cols.size());
- const size_type loc_row = global_rows.local_row(i);
-
- // fast function if there are no indirect references to any of the local
- // rows at all on this set of dofs (saves a lot of checks). the only check
- // we actually need to perform is whether the matrix element is zero.
- if (global_rows.have_indirect_rows() == false &&
- global_cols.have_indirect_rows() == false)
- {
- AssertIndexRange(loc_row, local_matrix.m());
- const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
-
- for (size_type j = column_start; j < column_end; ++j)
- {
- const size_type loc_col = global_cols.local_row(j);
- AssertIndexRange(loc_col, local_matrix.n());
- const LocalType col_val = matrix_ptr[loc_col];
- if (col_val != LocalType())
- {
- *val_ptr++ = static_cast<number>(col_val);
- *col_ptr++ = global_cols.global_row(j);
- }
- }
- }
-
- // more difficult part when there are indirect references and when we need
- // to do some more checks.
- else
- {
- for (size_type j = column_start; j < column_end; ++j)
- {
- LocalType col_val = resolve_matrix_entry(
- global_rows, global_cols, i, j, loc_row, local_matrix);
-
- // if we got some nontrivial value, append it to the array of
- // values.
- if (col_val != LocalType())
- {
- *val_ptr++ = static_cast<number>(col_val);
- *col_ptr++ = global_cols.global_row(j);
- }
- }
- }
- }
-
-
-
- // specialized function that can write into the row of a
- // SparseMatrix<number>.
- namespace dealiiSparseMatrix
- {
- template <typename SparseMatrixIterator, typename LocalType>
- static inline void
- add_value(const LocalType value,
- const size_type row,
- const size_type column,
- SparseMatrixIterator &matrix_values)
- {
- (void)row;
- if (value != LocalType())
- {
- while (matrix_values->column() < column)
- ++matrix_values;
- Assert(
- matrix_values->column() == column,
- typename SparseMatrix<typename SparseMatrixIterator::MatrixType::
- value_type>::ExcInvalidIndex(row, column));
- matrix_values->value() += value;
- }
- }
- } // namespace dealiiSparseMatrix
-
-
- // similar as before, now with shortcut for deal.II sparse matrices. this
- // lets us avoid using extra arrays, and does all the operations just in
- // place, i.e., in the respective matrix row
- template <typename number, typename LocalType>
- inline void
- resolve_matrix_row(const GlobalRowsFromLocal & global_rows,
- const size_type i,
- const size_type column_start,
- const size_type column_end,
- const FullMatrix<LocalType> &local_matrix,
- SparseMatrix<number> * sparse_matrix)
- {
- if (column_end == column_start)
- return;
-
- AssertIndexRange(column_end - 1, global_rows.size());
- const SparsityPattern &sparsity = sparse_matrix->get_sparsity_pattern();
-
- if (sparsity.n_nonzero_elements() == 0)
- return;
-
- const size_type row = global_rows.global_row(i);
- const size_type loc_row = global_rows.local_row(i);
-
- typename SparseMatrix<number>::iterator matrix_values =
- sparse_matrix->begin(row);
- const bool optimize_diagonal = sparsity.n_rows() == sparsity.n_cols();
-
- // distinguish three cases about what can happen for checking whether the
- // diagonal is the first element of the row. this avoids if statements at
- // the innermost loop positions
-
- if (!optimize_diagonal) // case 1: no diagonal optimization in matrix
- {
- if (global_rows.have_indirect_rows() == false)
- {
- AssertIndexRange(loc_row, local_matrix.m());
- const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
-
- for (size_type j = column_start; j < column_end; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- const LocalType col_val = matrix_ptr[loc_col];
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- else
- {
- for (size_type j = column_start; j < column_end; ++j)
- {
- LocalType col_val = resolve_matrix_entry(
- global_rows, global_rows, i, j, loc_row, local_matrix);
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- }
- else if (i >= column_start && i < column_end) // case 2: can split loop
- {
- ++matrix_values; // jump over diagonal element
- if (global_rows.have_indirect_rows() == false)
- {
- AssertIndexRange(loc_row, local_matrix.m());
- const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
-
- sparse_matrix->begin(row)->value() += matrix_ptr[loc_row];
- for (size_type j = column_start; j < i; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- const LocalType col_val = matrix_ptr[loc_col];
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- for (size_type j = i + 1; j < column_end; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- const LocalType col_val = matrix_ptr[loc_col];
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- else
- {
- sparse_matrix->begin(row)->value() += resolve_matrix_entry(
- global_rows, global_rows, i, i, loc_row, local_matrix);
- for (size_type j = column_start; j < i; ++j)
- {
- LocalType col_val = resolve_matrix_entry(
- global_rows, global_rows, i, j, loc_row, local_matrix);
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- for (size_type j = i + 1; j < column_end; ++j)
- {
- LocalType col_val = resolve_matrix_entry(
- global_rows, global_rows, i, j, loc_row, local_matrix);
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- }
- // case 3: can't say - need to check inside the loop
- else if (global_rows.have_indirect_rows() == false)
- {
- ++matrix_values; // jump over diagonal element
- AssertIndexRange(loc_row, local_matrix.m());
- const LocalType *matrix_ptr = &local_matrix(loc_row, 0);
-
- for (size_type j = column_start; j < column_end; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- const LocalType col_val = matrix_ptr[loc_col];
- if (row == global_rows.global_row(j))
- sparse_matrix->begin(row)->value() += col_val;
- else
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- else
- {
- ++matrix_values; // jump over diagonal element
- for (size_type j = column_start; j < column_end; ++j)
- {
- LocalType col_val = resolve_matrix_entry(
- global_rows, global_rows, i, j, loc_row, local_matrix);
- if (row == global_rows.global_row(j))
- sparse_matrix->begin(row)->value() += col_val;
- else
- dealiiSparseMatrix::add_value(
- col_val, row, global_rows.global_row(j), matrix_values);
- }
- }
- }
-
-
-
- // Same function to resolve all entries that will be added to the given
- // global row global_rows[i] as before, now for sparsity pattern
- inline void
- resolve_matrix_row(const GlobalRowsFromLocal & global_rows,
- const size_type i,
- const size_type column_start,
- const size_type column_end,
- const Table<2, bool> & dof_mask,
- std::vector<size_type>::iterator &col_ptr)
- {
- if (column_end == column_start)
- return;
-
- const size_type loc_row = global_rows.local_row(i);
-
- // fast function if there are no indirect references to any of the local
- // rows at all on this set of dofs
- if (global_rows.have_indirect_rows() == false)
- {
- Assert(loc_row < dof_mask.n_rows(), ExcInternalError());
-
- for (size_type j = column_start; j < column_end; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
-
- if (dof_mask(loc_row, loc_col) == true)
- *col_ptr++ = global_rows.global_row(j);
- }
- }
-
- // slower functions when there are indirect references and when we need to
- // do some more checks.
- else
- {
- for (size_type j = column_start; j < column_end; ++j)
- {
- const size_type loc_col = global_rows.local_row(j);
- if (loc_row != numbers::invalid_size_type)
- {
- Assert(loc_row < dof_mask.n_rows(), ExcInternalError());
- if (loc_col != numbers::invalid_size_type)
- {
- Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
- if (dof_mask(loc_row, loc_col) == true)
- goto add_this_index;
- }
-
- for (size_type p = 0; p < global_rows.size(j); ++p)
- if (dof_mask(loc_row, global_rows.local_row(j, p)) == true)
- goto add_this_index;
- }
-
- for (size_type q = 0; q < global_rows.size(i); ++q)
- {
- if (loc_col != numbers::invalid_size_type)
- {
- Assert(loc_col < dof_mask.n_cols(), ExcInternalError());
- if (dof_mask(global_rows.local_row(i, q), loc_col) == true)
- goto add_this_index;
- }
-
- for (size_type p = 0; p < global_rows.size(j); ++p)
- if (dof_mask(global_rows.local_row(i, q),
- global_rows.local_row(j, p)) == true)
- goto add_this_index;
- }
-
- continue;
- // if we got some nontrivial value, append it to the array of
- // values.
- add_this_index:
- *col_ptr++ = global_rows.global_row(j);
- }
- }
- }
-
-
-
- // to make sure that the global matrix remains invertible, we need to do
- // something with the diagonal elements. add the absolute value of the local
- // matrix, so the resulting entry will always be positive and furthermore be
- // in the same order of magnitude as the other elements of the matrix
- //
- // note that this also captures the special case that a dof is both
- // constrained and fixed (this can happen for hanging nodes in 3d that also
- // happen to be on the boundary). in that case, following the program flow
- // in distribute_local_to_global, it is realized that when distributing the
- // row and column no elements of the matrix are actually touched if all the
- // degrees of freedom to which this dof is constrained are also constrained
- // (the usual case with hanging nodes in 3d). however, in the line below, we
- // do actually do something with this dof
- template <typename MatrixType, typename VectorType>
- inline void
- set_matrix_diagonals(
- const internals::GlobalRowsFromLocal & global_rows,
- const std::vector<size_type> & local_dof_indices,
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const ConstraintMatrix & constraints,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs)
- {
- if (global_rows.n_constraints() > 0)
- {
- typename MatrixType::value_type average_diagonal =
- typename MatrixType::value_type();
- for (size_type i = 0; i < local_matrix.m(); ++i)
- average_diagonal += std::abs(local_matrix(i, i));
- average_diagonal /= static_cast<double>(local_matrix.m());
-
- for (size_type i = 0; i < global_rows.n_constraints(); i++)
- {
- const size_type local_row = global_rows.constraint_origin(i);
- const size_type global_row = local_dof_indices[local_row];
- const typename MatrixType::value_type new_diagonal =
- (std::abs(local_matrix(local_row, local_row)) != 0 ?
- std::abs(local_matrix(local_row, local_row)) :
- average_diagonal);
- global_matrix.add(global_row, global_row, new_diagonal);
-
- // if the use_inhomogeneities_for_rhs flag is set to true, the
- // inhomogeneities are used to create the global vector. instead
- // of fill in a zero in the ith components with an inhomogeneity,
- // we set those to: inhomogeneity(i)*global_matrix (i,i).
- if (use_inhomogeneities_for_rhs == true)
- global_vector(global_row) +=
- new_diagonal * constraints.get_inhomogeneity(global_row);
- }
- }
- }
-
-
-
- // similar function as the one above for setting matrix diagonals, but now
- // doing that for sparsity patterns when setting them up using
- // add_entries_local_to_global. In case we keep constrained entries, add all
- // the rows and columns related to the constrained dof, otherwise just add
- // the diagonal
- template <typename SparsityPatternType>
- inline void
- set_sparsity_diagonals(const internals::GlobalRowsFromLocal &global_rows,
- const std::vector<size_type> &local_dof_indices,
- const Table<2, bool> & dof_mask,
- const bool keep_constrained_entries,
- SparsityPatternType & sparsity_pattern)
- {
- // if we got constraints, need to add the diagonal element and, if the
- // user requested so, also the rest of the entries in rows and columns
- // that have been left out above
- if (global_rows.n_constraints() > 0)
- {
- for (size_type i = 0; i < global_rows.n_constraints(); i++)
- {
- const size_type local_row = global_rows.constraint_origin(i);
- const size_type global_row = local_dof_indices[local_row];
- if (keep_constrained_entries == true)
- {
- for (size_type j = 0; j < local_dof_indices.size(); ++j)
- {
- if (dof_mask(local_row, j) == true)
- sparsity_pattern.add(global_row, local_dof_indices[j]);
- if (dof_mask(j, local_row) == true)
- sparsity_pattern.add(local_dof_indices[j], global_row);
- }
- }
- else
- // don't keep constrained entries - just add the diagonal.
- sparsity_pattern.add(global_row, global_row);
- }
- }
- }
-
-} // end of namespace internals
-
-
-
-// Basic idea of setting up a list of
-// all global dofs: first find all rows and columns
-// that we are going to write touch,
-// and then go through the
-// lines and collect all the local rows that
-// are related to it.
-void
-ConstraintMatrix::make_sorted_row_list(
- const std::vector<size_type> & local_dof_indices,
- internals::GlobalRowsFromLocal &global_rows) const
-{
- const size_type n_local_dofs = local_dof_indices.size();
- AssertDimension(n_local_dofs, global_rows.size());
-
- // when distributing the local data to the global matrix, we can quite
- // cheaply sort the indices (obviously, this introduces the need for
- // allocating some memory on the way, but we need to do this only for rows,
- // whereas the distribution process itself goes over rows and columns). This
- // has the advantage that when writing into the global matrix, we can make
- // use of the sortedness.
-
- // so the first step is to create a sorted list of all row values that are
- // possible. these values are either the rows from unconstrained dofs, or
- // some indices introduced by dofs constrained to a combination of some
- // other dofs. regarding the data type, choose a <tt>std::vector</tt> of a
- // pair of unsigned ints (for global columns) and internal data (containing
- // local columns + possible jumps from constraints). Choosing
- // <tt>std::map</tt> or anything else M.K. knows of would be much more
- // expensive here!
-
- // cache whether we have to resolve any indirect rows generated from
- // resolving constrained dofs.
- size_type added_rows = 0;
-
- // first add the indices in an unsorted way and only keep track of the
- // constraints that appear. They are resolved in a second step.
- for (size_type i = 0; i < n_local_dofs; ++i)
- {
- if (is_constrained(local_dof_indices[i]) == false)
- {
- global_rows.global_row(added_rows) = local_dof_indices[i];
- global_rows.local_row(added_rows++) = i;
- }
- else
- global_rows.insert_constraint(i);
- }
- global_rows.sort();
-
- const size_type n_constrained_rows = n_local_dofs - added_rows;
- for (size_type i = 0; i < n_constrained_rows; ++i)
- {
- const size_type local_row = global_rows.constraint_origin(i);
- AssertIndexRange(local_row, n_local_dofs);
- const size_type global_row = local_dof_indices[local_row];
- Assert(is_constrained(global_row), ExcInternalError());
- const ConstraintLine &position =
- lines[lines_cache[calculate_line_index(global_row)]];
- if (position.inhomogeneity != 0)
- global_rows.set_ith_constraint_inhomogeneous(i);
- for (size_type q = 0; q < position.entries.size(); ++q)
- global_rows.insert_index(
- position.entries[q].first, local_row, position.entries[q].second);
- }
-}
-
-
-
-// Same function as before, but now do only extract the global indices that
-// come from the local ones without storing their origin. Used for sparsity
-// pattern generation.
-inline void
-ConstraintMatrix::make_sorted_row_list(
- const std::vector<size_type> &local_dof_indices,
- std::vector<size_type> & active_dofs) const
-{
- const size_type n_local_dofs = local_dof_indices.size();
- size_type added_rows = 0;
- for (size_type i = 0; i < n_local_dofs; ++i)
- {
- if (is_constrained(local_dof_indices[i]) == false)
- {
- active_dofs[added_rows++] = local_dof_indices[i];
- continue;
- }
-
- active_dofs[n_local_dofs - i + added_rows - 1] = i;
- }
- std::sort(active_dofs.begin(), active_dofs.begin() + added_rows);
-
- const size_type n_constrained_dofs = n_local_dofs - added_rows;
- for (size_type i = n_constrained_dofs; i > 0; --i)
- {
- const size_type local_row = active_dofs.back();
-
- // remove constrained entry since we are going to resolve it in place
- active_dofs.pop_back();
- const size_type global_row = local_dof_indices[local_row];
- const ConstraintLine &position =
- lines[lines_cache[calculate_line_index(global_row)]];
- for (size_type q = 0; q < position.entries.size(); ++q)
- {
- const size_type new_index = position.entries[q].first;
- if (active_dofs[active_dofs.size() - i] < new_index)
- active_dofs.insert(active_dofs.end() - i + 1, new_index);
-
- // make binary search to find where to put the new index in order to
- // keep the list sorted
- else
- {
- std::vector<size_type>::iterator it = Utilities::lower_bound(
- active_dofs.begin(), active_dofs.end() - i + 1, new_index);
- if (*it != new_index)
- active_dofs.insert(it, new_index);
- }
- }
- }
-}
-
-
-
-// Resolve the constraints from the vector and apply inhomogeneities.
-template <typename MatrixScalar, typename VectorScalar>
-inline typename ProductType<VectorScalar, MatrixScalar>::type
-ConstraintMatrix::resolve_vector_entry(
- const size_type i,
- const internals::GlobalRowsFromLocal &global_rows,
- const Vector<VectorScalar> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- const FullMatrix<MatrixScalar> & local_matrix) const
-{
- const size_type loc_row = global_rows.local_row(i);
- const size_type n_inhomogeneous_rows = global_rows.n_inhomogeneities();
- typename ProductType<VectorScalar, MatrixScalar>::type val = 0;
- // has a direct contribution from some local entry. If we have inhomogeneous
- // constraints, compute the contribution of the inhomogeneity in the current
- // row.
- if (loc_row != numbers::invalid_size_type)
- {
- val = local_vector(loc_row);
- for (size_type i = 0; i < n_inhomogeneous_rows; ++i)
- val -= (local_matrix(loc_row, global_rows.constraint_origin(i)) *
- lines[lines_cache[calculate_line_index(
- local_dof_indices[global_rows.constraint_origin(i)])]]
- .inhomogeneity);
- }
-
- // go through the indirect contributions
- for (size_type q = 0; q < global_rows.size(i); ++q)
- {
- const size_type loc_row_q = global_rows.local_row(i, q);
- typename ProductType<VectorScalar, MatrixScalar>::type add_this =
- local_vector(loc_row_q);
- for (size_type k = 0; k < n_inhomogeneous_rows; ++k)
- add_this -=
- (local_matrix(loc_row_q, global_rows.constraint_origin(k)) *
- lines[lines_cache[calculate_line_index(
- local_dof_indices[global_rows.constraint_origin(k)])]]
- .inhomogeneity);
- val += add_this * global_rows.constraint_value(i, q);
- }
- return val;
-}
-
-
-// internal implementation for distribute_local_to_global for standard
-// (non-block) matrices
-template <typename MatrixType, typename VectorType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs,
- std::integral_constant<bool, false>) const
-{
- // check whether we work on real vectors or we just used a dummy when
- // calling the other function above.
- const bool use_vectors =
- (local_vector.size() == 0 && global_vector.size() == 0) ? false : true;
- typedef typename MatrixType::value_type number;
- const bool use_dealii_matrix =
- std::is_same<MatrixType, SparseMatrix<number>>::value;
-
- AssertDimension(local_matrix.n(), local_dof_indices.size());
- AssertDimension(local_matrix.m(), local_dof_indices.size());
- Assert(global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
- if (use_vectors == true)
- {
- AssertDimension(local_matrix.m(), local_vector.size());
- AssertDimension(global_matrix.m(), global_vector.size());
- }
- Assert(lines.empty() || sorted == true, ExcMatrixNotClosed());
-
- const size_type n_local_dofs = local_dof_indices.size();
-
- typename internals::ConstraintMatrixData<
- typename MatrixType::value_type,
- typename VectorType::value_type>::ScratchDataAccessor scratch_data;
-
- internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
- global_rows.reinit(n_local_dofs);
- make_sorted_row_list(local_dof_indices, global_rows);
-
- const size_type n_actual_dofs = global_rows.size();
-
- // create arrays for the column data (indices and values) that will then be
- // written into the matrix. Shortcut for deal.II sparse matrix. We can use
- // the scratch data if we have a double matrix. Otherwise, we need to create
- // an array in any case since we cannot know about the actual data type in
- // the ConstraintMatrix class (unless we do cast). This involves a little
- // bit of logic to determine the type of the matrix value.
- std::vector<size_type> &cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
- // create arrays for writing into the vector as well
- std::vector<size_type> &vector_indices = scratch_data->vector_indices;
- std::vector<typename VectorType::value_type> &vector_values =
- scratch_data->vector_values;
- vector_indices.resize(n_actual_dofs);
- vector_values.resize(n_actual_dofs);
- SparseMatrix<number> *sparse_matrix =
- dynamic_cast<SparseMatrix<number> *>(&global_matrix);
- if (use_dealii_matrix == false)
- {
- cols.resize(n_actual_dofs);
- vals.resize(n_actual_dofs);
- }
- else
- Assert(sparse_matrix != nullptr, ExcInternalError());
-
- // now do the actual job. go through all the global rows that we will touch
- // and call resolve_matrix_row for each of those.
- size_type local_row_n = 0;
- for (size_type i = 0; i < n_actual_dofs; ++i)
- {
- const size_type row = global_rows.global_row(i);
-
- // calculate all the data that will be written into the matrix row.
- if (use_dealii_matrix == false)
- {
- size_type *col_ptr = &cols[0];
- // cast is uncritical here and only used to avoid compiler
- // warnings. We never access a non-double array
- number *val_ptr = &vals[0];
- internals::resolve_matrix_row(global_rows,
- global_rows,
- i,
- 0,
- n_actual_dofs,
- local_matrix,
- col_ptr,
- val_ptr);
- const size_type n_values = col_ptr - &cols[0];
- if (n_values > 0)
- global_matrix.add(row, n_values, &cols[0], &vals[0], false, true);
- }
- else
- internals::resolve_matrix_row(
- global_rows, i, 0, n_actual_dofs, local_matrix, sparse_matrix);
-
- // now to the vectors. besides doing the same job as we did above (i.e.,
- // distribute the content of the local vector into the global one), need
- // to account for inhomogeneities here: this corresponds to eliminating
- // the respective column in the local matrix with value on the right
- // hand side.
- if (use_vectors == true)
- {
- const typename VectorType::value_type val = resolve_vector_entry(
- i, global_rows, local_vector, local_dof_indices, local_matrix);
- AssertIsFinite(val);
-
- if (val != typename VectorType::value_type())
- {
- vector_indices[local_row_n] = row;
- vector_values[local_row_n] = val;
- ++local_row_n;
- }
- }
- }
- // Drop the elements of vector_indices and vector_values that we do not use
- // (we may always elide writing zero values to vectors)
- const size_type n_local_rows = local_row_n;
- vector_indices.resize(n_local_rows);
- vector_values.resize(n_local_rows);
-
- // While the standard case is that these types are equal, they need not be, so
- // only do a bulk update if they are. Note that the types in the arguments to
- // add must be equal if we have a Trilinos or PETSc vector but do not have to
- // be if we have a deal.II native vector: one could further optimize this for
- // Vector, LinearAlgebra::distributed::vector, etc.
- if (std::is_same<typename VectorType::value_type, number>::value)
- {
- global_vector.add(
- vector_indices,
- *reinterpret_cast<std::vector<number> *>(&vector_values));
- }
- else
- {
- for (size_type row_n = 0; row_n < n_local_rows; ++row_n)
- {
- global_vector(vector_indices[row_n]) +=
- static_cast<typename VectorType::value_type>(vector_values[row_n]);
- }
- }
-
- internals::set_matrix_diagonals(global_rows,
- local_dof_indices,
- local_matrix,
- *this,
- global_matrix,
- global_vector,
- use_inhomogeneities_for_rhs);
-}
-
-
-
-// similar function as above, but now specialized for block matrices. See the
-// other function for additional comments.
-template <typename MatrixType, typename VectorType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const Vector<typename VectorType::value_type> & local_vector,
- const std::vector<size_type> & local_dof_indices,
- MatrixType & global_matrix,
- VectorType & global_vector,
- bool use_inhomogeneities_for_rhs,
- std::integral_constant<bool, true>) const
-{
- const bool use_vectors =
- (local_vector.size() == 0 && global_vector.size() == 0) ? false : true;
- typedef typename MatrixType::value_type number;
- const bool use_dealii_matrix =
- std::is_same<MatrixType, BlockSparseMatrix<number>>::value;
-
- AssertDimension(local_matrix.n(), local_dof_indices.size());
- AssertDimension(local_matrix.m(), local_dof_indices.size());
- Assert(global_matrix.m() == global_matrix.n(), ExcNotQuadratic());
- Assert(global_matrix.n_block_rows() == global_matrix.n_block_cols(),
- ExcNotQuadratic());
- if (use_vectors == true)
- {
- AssertDimension(local_matrix.m(), local_vector.size());
- AssertDimension(global_matrix.m(), global_vector.size());
- }
- Assert(sorted == true, ExcMatrixNotClosed());
-
- typename internals::ConstraintMatrixData<
- typename MatrixType::value_type,
- typename VectorType::value_type>::ScratchDataAccessor scratch_data;
-
- const size_type n_local_dofs = local_dof_indices.size();
- internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
- global_rows.reinit(n_local_dofs);
-
- make_sorted_row_list(local_dof_indices, global_rows);
- const size_type n_actual_dofs = global_rows.size();
-
- std::vector<size_type> &global_indices = scratch_data->vector_indices;
- if (use_vectors == true)
- {
- global_indices.resize(n_actual_dofs);
- for (size_type i = 0; i < n_actual_dofs; ++i)
- global_indices[i] = global_rows.global_row(i);
- }
-
- // additional construct that also takes care of block indices.
- const size_type num_blocks = global_matrix.n_block_rows();
- std::vector<size_type> &block_starts = scratch_data->block_starts;
- block_starts.resize(num_blocks + 1);
- internals::make_block_starts(global_matrix, global_rows, block_starts);
-
- std::vector<size_type> &cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
- if (use_dealii_matrix == false)
- {
- cols.resize(n_actual_dofs);
- vals.resize(n_actual_dofs);
- }
-
- // the basic difference to the non-block variant from now onwards is that we
- // go through the blocks of the matrix separately, which allows us to set
- // the block entries individually
- for (size_type block = 0; block < num_blocks; ++block)
- {
- const size_type next_block = block_starts[block + 1];
- for (size_type i = block_starts[block]; i < next_block; ++i)
- {
- const size_type row = global_rows.global_row(i);
-
- for (size_type block_col = 0; block_col < num_blocks; ++block_col)
- {
- const size_type start_block = block_starts[block_col],
- end_block = block_starts[block_col + 1];
- if (use_dealii_matrix == false)
- {
- size_type *col_ptr = &cols[0];
- number * val_ptr = &vals[0];
- internals::resolve_matrix_row(global_rows,
- global_rows,
- i,
- start_block,
- end_block,
- local_matrix,
- col_ptr,
- val_ptr);
- const size_type n_values = col_ptr - &cols[0];
- if (n_values > 0)
- global_matrix.block(block, block_col)
- .add(row, n_values, &cols[0], &vals[0], false, true);
- }
- else
- {
- SparseMatrix<number> *sparse_matrix =
- dynamic_cast<SparseMatrix<number> *>(
- &global_matrix.block(block, block_col));
- Assert(sparse_matrix != nullptr, ExcInternalError());
- internals::resolve_matrix_row(global_rows,
- i,
- start_block,
- end_block,
- local_matrix,
- sparse_matrix);
- }
- }
-
- if (use_vectors == true)
- {
- const number val = resolve_vector_entry(
- i, global_rows, local_vector, local_dof_indices, local_matrix);
-
- if (val != number())
- global_vector(global_indices[i]) +=
- static_cast<typename VectorType::value_type>(val);
- }
- }
- }
-
- internals::set_matrix_diagonals(global_rows,
- local_dof_indices,
- local_matrix,
- *this,
- global_matrix,
- global_vector,
- use_inhomogeneities_for_rhs);
-}
-
-
-
-template <typename MatrixType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & row_indices,
- const std::vector<size_type> & col_indices,
- MatrixType & global_matrix) const
-{
- distribute_local_to_global(
- local_matrix, row_indices, *this, col_indices, global_matrix);
-}
-
-
-
-template <typename MatrixType>
-void
-ConstraintMatrix::distribute_local_to_global(
- const FullMatrix<typename MatrixType::value_type> &local_matrix,
- const std::vector<size_type> & row_indices,
- const ConstraintMatrix & col_constraint_matrix,
- const std::vector<size_type> & col_indices,
- MatrixType & global_matrix) const
-{
- typedef typename MatrixType::value_type number;
-
- AssertDimension(local_matrix.m(), row_indices.size());
- AssertDimension(local_matrix.n(), col_indices.size());
-
- const size_type n_local_row_dofs = row_indices.size();
- const size_type n_local_col_dofs = col_indices.size();
-
- typename internals::ConstraintMatrixData<
- typename MatrixType::value_type>::ScratchDataAccessor scratch_data;
- internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
- global_rows.reinit(n_local_row_dofs);
- internals::GlobalRowsFromLocal &global_cols = scratch_data->global_columns;
- global_cols.reinit(n_local_col_dofs);
- make_sorted_row_list(row_indices, global_rows);
- col_constraint_matrix.make_sorted_row_list(col_indices, global_cols);
-
- const size_type n_actual_row_dofs = global_rows.size();
- const size_type n_actual_col_dofs = global_cols.size();
-
- // create arrays for the column data (indices and values) that will then be
- // written into the matrix. Shortcut for deal.II sparse matrix
- std::vector<size_type> &cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
- cols.resize(n_actual_col_dofs);
- vals.resize(n_actual_col_dofs);
-
- // now do the actual job.
- for (size_type i = 0; i < n_actual_row_dofs; ++i)
- {
- const size_type row = global_rows.global_row(i);
-
- // calculate all the data that will be written into the matrix row.
- size_type *col_ptr = &cols[0];
- number * val_ptr = &vals[0];
- internals::resolve_matrix_row(global_rows,
- global_cols,
- i,
- 0,
- n_actual_col_dofs,
- local_matrix,
- col_ptr,
- val_ptr);
- const size_type n_values = col_ptr - &cols[0];
- if (n_values > 0)
- global_matrix.add(row, n_values, &cols[0], &vals[0], false, true);
- }
-}
-
-
-
-template <typename SparsityPatternType>
-void
-ConstraintMatrix::add_entries_local_to_global(
- const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> & dof_mask,
- std::integral_constant<bool, false>) const
-{
- Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
- ExcNotQuadratic());
-
- const size_type n_local_dofs = local_dof_indices.size();
- bool dof_mask_is_active = false;
- if (dof_mask.n_rows() == n_local_dofs)
- {
- dof_mask_is_active = true;
- AssertDimension(dof_mask.n_cols(), n_local_dofs);
- }
-
- internals::ConstraintMatrixData<double>::ScratchDataAccessor scratch_data;
-
- // if the dof mask is not active, all we have to do is to add some indices
- // in a matrix format. To do this, we first create an array of all the
- // indices that are to be added. these indices are the local dof indices
- // plus some indices that come from constraints.
- if (dof_mask_is_active == false)
- {
- std::vector<size_type> &actual_dof_indices = scratch_data->columns;
- actual_dof_indices.resize(n_local_dofs);
- make_sorted_row_list(local_dof_indices, actual_dof_indices);
- const size_type n_actual_dofs = actual_dof_indices.size();
-
- // now add the indices we collected above to the sparsity pattern. Very
- // easy here - just add the same array to all the rows...
- for (size_type i = 0; i < n_actual_dofs; ++i)
- sparsity_pattern.add_entries(actual_dof_indices[i],
- actual_dof_indices.begin(),
- actual_dof_indices.end(),
- true);
-
- // need to add the whole row and column structure in case we keep
- // constrained entries. Unfortunately, we can't use the nice matrix
- // structure we use elsewhere, so manually add those indices one by one.
- for (size_type i = 0; i < n_local_dofs; i++)
- if (is_constrained(local_dof_indices[i]))
- {
- if (keep_constrained_entries == true)
- for (size_type j = 0; j < n_local_dofs; j++)
- {
- sparsity_pattern.add(local_dof_indices[i],
- local_dof_indices[j]);
- sparsity_pattern.add(local_dof_indices[j],
- local_dof_indices[i]);
- }
- else
- sparsity_pattern.add(local_dof_indices[i], local_dof_indices[i]);
- }
-
- return;
- }
-
-
- // complicated case: we need to filter out some indices. then the function
- // gets similar to the function for distributing matrix entries, see there
- // for additional comments.
- internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
- global_rows.reinit(n_local_dofs);
- make_sorted_row_list(local_dof_indices, global_rows);
- const size_type n_actual_dofs = global_rows.size();
-
- // create arrays for the column indices that will then be written into the
- // sparsity pattern.
- std::vector<size_type> &cols = scratch_data->columns;
- cols.resize(n_actual_dofs);
-
- for (size_type i = 0; i < n_actual_dofs; ++i)
- {
- std::vector<size_type>::iterator col_ptr = cols.begin();
- const size_type row = global_rows.global_row(i);
- internals::resolve_matrix_row(
- global_rows, i, 0, n_actual_dofs, dof_mask, col_ptr);
-
- // finally, write all the information that accumulated under the given
- // process into the global matrix row and into the vector
- if (col_ptr != cols.begin())
- sparsity_pattern.add_entries(row, cols.begin(), col_ptr, true);
- }
- internals::set_sparsity_diagonals(global_rows,
- local_dof_indices,
- dof_mask,
- keep_constrained_entries,
- sparsity_pattern);
-}
-
-
-
-template <typename SparsityPatternType>
-void
-ConstraintMatrix::add_entries_local_to_global(
- const std::vector<size_type> &row_indices,
- const std::vector<size_type> &col_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> & dof_mask) const
-{
- const size_type n_local_rows = row_indices.size();
- const size_type n_local_cols = col_indices.size();
- bool dof_mask_is_active = false;
- if (dof_mask.n_rows() == n_local_rows && dof_mask.n_cols() == n_local_cols)
- dof_mask_is_active = true;
-
- // if constrained entries should be kept, need to add rows and columns of
- // those to the sparsity pattern
- if (keep_constrained_entries == true)
- {
- for (size_type i = 0; i < row_indices.size(); i++)
- if (is_constrained(row_indices[i]))
- for (size_type j = 0; j < col_indices.size(); j++)
- sparsity_pattern.add(row_indices[i], col_indices[j]);
- for (size_type i = 0; i < col_indices.size(); i++)
- if (is_constrained(col_indices[i]))
- for (size_type j = 0; j < row_indices.size(); j++)
- sparsity_pattern.add(row_indices[j], col_indices[i]);
- }
-
- // if the dof mask is not active, all we have to do is to add some indices
- // in a matrix format. To do this, we first create an array of all the
- // indices that are to be added. these indices are the local dof indices
- // plus some indices that come from constraints.
- if (dof_mask_is_active == false)
- {
- std::vector<size_type> actual_row_indices(n_local_rows);
- std::vector<size_type> actual_col_indices(n_local_cols);
- make_sorted_row_list(row_indices, actual_row_indices);
- make_sorted_row_list(col_indices, actual_col_indices);
- const size_type n_actual_rows = actual_row_indices.size();
-
- // now add the indices we collected above to the sparsity pattern. Very
- // easy here - just add the same array to all the rows...
- for (size_type i = 0; i < n_actual_rows; ++i)
- sparsity_pattern.add_entries(actual_row_indices[i],
- actual_col_indices.begin(),
- actual_col_indices.end(),
- true);
- return;
- }
-
-
- // TODO: implement this
- Assert(false, ExcNotImplemented());
-}
-
-
-
-template <typename SparsityPatternType>
-void
-ConstraintMatrix::add_entries_local_to_global(
- const std::vector<size_type> &local_dof_indices,
- SparsityPatternType & sparsity_pattern,
- const bool keep_constrained_entries,
- const Table<2, bool> & dof_mask,
- std::integral_constant<bool, true>) const
-{
- // just as the other add_entries_local_to_global function, but now
- // specialized for block matrices.
- Assert(sparsity_pattern.n_rows() == sparsity_pattern.n_cols(),
- ExcNotQuadratic());
- Assert(sparsity_pattern.n_block_rows() == sparsity_pattern.n_block_cols(),
- ExcNotQuadratic());
-
- const size_type n_local_dofs = local_dof_indices.size();
- const size_type num_blocks = sparsity_pattern.n_block_rows();
-
- internals::ConstraintMatrixData<double>::ScratchDataAccessor scratch_data;
-
- bool dof_mask_is_active = false;
- if (dof_mask.n_rows() == n_local_dofs)
- {
- dof_mask_is_active = true;
- AssertDimension(dof_mask.n_cols(), n_local_dofs);
- }
-
- if (dof_mask_is_active == false)
- {
- std::vector<size_type> &actual_dof_indices = scratch_data->columns;
- actual_dof_indices.resize(n_local_dofs);
- make_sorted_row_list(local_dof_indices, actual_dof_indices);
- const size_type n_actual_dofs = actual_dof_indices.size();
- (void)n_actual_dofs;
-
- // additional construct that also takes care of block indices.
- std::vector<size_type> &block_starts = scratch_data->block_starts;
- block_starts.resize(num_blocks + 1);
- internals::make_block_starts(
- sparsity_pattern, actual_dof_indices, block_starts);
-
- for (size_type block = 0; block < num_blocks; ++block)
- {
- const size_type next_block = block_starts[block + 1];
- for (size_type i = block_starts[block]; i < next_block; ++i)
- {
- Assert(i < n_actual_dofs, ExcInternalError());
- const size_type row = actual_dof_indices[i];
- Assert(row < sparsity_pattern.block(block, 0).n_rows(),
- ExcInternalError());
- std::vector<size_type>::iterator index_it =
- actual_dof_indices.begin();
- for (size_type block_col = 0; block_col < num_blocks; ++block_col)
- {
- const size_type next_block_col = block_starts[block_col + 1];
- sparsity_pattern.block(block, block_col)
- .add_entries(row,
- index_it,
- actual_dof_indices.begin() + next_block_col,
- true);
- index_it = actual_dof_indices.begin() + next_block_col;
- }
- }
- }
-
- for (size_type i = 0; i < n_local_dofs; i++)
- if (is_constrained(local_dof_indices[i]))
- {
- if (keep_constrained_entries == true)
- for (size_type j = 0; j < n_local_dofs; j++)
- {
- sparsity_pattern.add(local_dof_indices[i],
- local_dof_indices[j]);
- sparsity_pattern.add(local_dof_indices[j],
- local_dof_indices[i]);
- }
- else
- sparsity_pattern.add(local_dof_indices[i], local_dof_indices[i]);
- }
-
- return;
- }
-
- // difficult case with dof_mask, similar to the distribute_local_to_global
- // function for block matrices
- internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
- global_rows.reinit(n_local_dofs);
- make_sorted_row_list(local_dof_indices, global_rows);
- const size_type n_actual_dofs = global_rows.size();
-
- // additional construct that also takes care of block indices.
- std::vector<size_type> &block_starts = scratch_data->block_starts;
- block_starts.resize(num_blocks + 1);
- internals::make_block_starts(sparsity_pattern, global_rows, block_starts);
-
- std::vector<size_type> &cols = scratch_data->columns;
- cols.resize(n_actual_dofs);
-
- // the basic difference to the non-block variant from now onwards is that we
- // go through the blocks of the matrix separately.
- for (size_type block = 0; block < num_blocks; ++block)
- {
- const size_type next_block = block_starts[block + 1];
- for (size_type i = block_starts[block]; i < next_block; ++i)
- {
- const size_type row = global_rows.global_row(i);
- for (size_type block_col = 0; block_col < num_blocks; ++block_col)
- {
- const size_type begin_block = block_starts[block_col],
- end_block = block_starts[block_col + 1];
- std::vector<size_type>::iterator col_ptr = cols.begin();
- internals::resolve_matrix_row(
- global_rows, i, begin_block, end_block, dof_mask, col_ptr);
-
- sparsity_pattern.block(block, block_col)
- .add_entries(row, cols.begin(), col_ptr, true);
- }
- }
- }
-
- internals::set_sparsity_diagonals(global_rows,
- local_dof_indices,
- dof_mask,
- keep_constrained_entries,
- sparsity_pattern);
-}
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
block_vector.cc
chunk_sparse_matrix.cc
chunk_sparsity_pattern.cc
- constraint_matrix.cc
dynamic_sparsity_pattern.cc
exceptions.cc
full_matrix.cc
block_sparse_matrix.inst.in
block_vector.inst.in
chunk_sparse_matrix.inst.in
- constraint_matrix.inst.in
full_matrix.inst.in
lapack_full_matrix.inst.in
la_vector.inst.in
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 1998 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/memory_consumption.h>
-
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_sparse_matrix_ez.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/chunk_sparse_matrix.h>
-#include <deal.II/lac/constraint_matrix.h>
-#include <deal.II/lac/constraint_matrix.templates.h>
-#include <deal.II/lac/diagonal_matrix.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/la_parallel_vector.h>
-#include <deal.II/lac/la_vector.h>
-#include <deal.II/lac/matrix_block.h>
-#include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
-#include <deal.II/lac/petsc_parallel_block_vector.h>
-#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/petsc_sparse_matrix.h>
-#include <deal.II/lac/sparse_matrix_ez.h>
-#include <deal.II/lac/trilinos_block_sparse_matrix.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
-#include <deal.II/lac/trilinos_sparse_matrix.h>
-#include <deal.II/lac/trilinos_vector.h>
-
-#include <boost/serialization/utility.hpp>
-
-#include <algorithm>
-#include <numeric>
-#include <ostream>
-#include <set>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-
-// Static member variable
-const Table<2, bool> ConstraintMatrix::default_empty_table = Table<2, bool>();
-
-
-
-void
-ConstraintMatrix::copy_from(const ConstraintMatrix &other)
-{
- lines = other.lines;
- lines_cache = other.lines_cache;
- local_lines = other.local_lines;
- sorted = other.sorted;
-}
-
-
-
-bool
-ConstraintMatrix::check_zero_weight(const std::pair<size_type, double> &p)
-{
- return (p.second == 0);
-}
-
-
-
-bool
-ConstraintMatrix::ConstraintLine::operator<(const ConstraintLine &a) const
-{
- return index < a.index;
-}
-
-
-
-bool
-ConstraintMatrix::ConstraintLine::operator==(const ConstraintLine &a) const
-{
- return index == a.index;
-}
-
-
-
-std::size_t
-ConstraintMatrix::ConstraintLine::memory_consumption() const
-{
- return (MemoryConsumption::memory_consumption(index) +
- MemoryConsumption::memory_consumption(entries) +
- MemoryConsumption::memory_consumption(inhomogeneity));
-}
-
-
-
-const ConstraintMatrix::LineRange
-ConstraintMatrix::get_lines() const
-{
- return boost::make_iterator_range(lines.begin(), lines.end());
-}
-
-
-
-bool
-ConstraintMatrix::is_consistent_in_parallel(
- const std::vector<IndexSet> &locally_owned_dofs,
- const IndexSet & locally_active_dofs,
- const MPI_Comm mpi_communicator,
- const bool verbose) const
-{
- ConstraintLine empty;
- empty.inhomogeneity = 0.0;
-
- // Helper to return a reference to the ConstraintLine object that belongs to row @p row.
- // We don't want to make copies but to return a reference, we need an empty
- // object that we store above.
- auto get_line = [&](const size_type row) -> const ConstraintLine & {
- const size_type line_index = calculate_line_index(row);
- if (line_index >= lines_cache.size() ||
- lines_cache[line_index] == numbers::invalid_size_type)
- {
- empty.index = row;
- return empty;
- }
- else
- return lines[lines_cache[line_index]];
- };
-
- // identify non-owned rows and send to owner:
- std::map<unsigned int, std::vector<ConstraintLine>> to_send;
-
- const unsigned int myid =
- dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
- const unsigned int nproc =
- dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
-
- // We will send all locally active dofs that are not locally owned for
- // checking. Note that we allow constraints to differ on locally_relevant (and
- // not active) DoFs.
- IndexSet non_owned = locally_active_dofs;
- non_owned.subtract_set(locally_owned_dofs[myid]);
- for (unsigned int owner = 0; owner < nproc; ++owner)
- {
- // find all lines to send to @p owner
- IndexSet indices_to_send = non_owned & locally_owned_dofs[owner];
- for (const auto &row_idx : indices_to_send)
- {
- to_send[owner].push_back(get_line(row_idx));
- }
- }
-
- std::map<unsigned int, std::vector<ConstraintLine>> received =
- Utilities::MPI::some_to_some(mpi_communicator, to_send);
-
- unsigned int inconsistent = 0;
-
- // from each processor:
- for (const auto &kv : received)
- {
- // for each incoming line:
- for (auto &lineit : kv.second)
- {
- const ConstraintLine &reference = get_line(lineit.index);
-
- if (lineit.inhomogeneity != reference.inhomogeneity)
- {
- ++inconsistent;
-
- if (verbose)
- std::cout << "Proc " << myid << " got line " << lineit.index
- << " from " << kv.first << " inhomogeneity "
- << lineit.inhomogeneity
- << " != " << reference.inhomogeneity << std::endl;
- }
- else if (lineit.entries != reference.entries)
- {
- ++inconsistent;
- if (verbose)
- std::cout << "Proc " << myid << " got line " << lineit.index
- << " from " << kv.first << " wrong values!"
- << std::endl;
- }
- }
- }
-
- const unsigned int total =
- Utilities::MPI::sum(inconsistent, mpi_communicator);
- if (verbose && total > 0 && myid == 0)
- std::cout << total << " inconsistent lines discovered!" << std::endl;
- return total == 0;
-}
-
-
-
-void
-ConstraintMatrix::add_lines(const std::set<size_type> &lines)
-{
- for (std::set<size_type>::const_iterator i = lines.begin(); i != lines.end();
- ++i)
- add_line(*i);
-}
-
-
-
-void
-ConstraintMatrix::add_lines(const std::vector<bool> &lines)
-{
- for (size_type i = 0; i < lines.size(); ++i)
- if (lines[i] == true)
- add_line(i);
-}
-
-
-
-void
-ConstraintMatrix::add_lines(const IndexSet &lines)
-{
- for (size_type i = 0; i < lines.n_elements(); ++i)
- add_line(lines.nth_index_in_set(i));
-}
-
-
-
-void
-ConstraintMatrix::add_entries(
- const size_type line,
- const std::vector<std::pair<size_type, double>> &col_val_pairs)
-{
- Assert(sorted == false, ExcMatrixIsClosed());
- Assert(is_constrained(line), ExcLineInexistant(line));
-
- ConstraintLine *line_ptr = &lines[lines_cache[calculate_line_index(line)]];
- Assert(line_ptr->index == line, ExcInternalError());
-
- // if in debug mode, check whether an entry for this column already
- // exists and if its the same as the one entered at present
- //
- // in any case: skip this entry if an entry for this column already
- // exists, since we don't want to enter it twice
- for (std::vector<std::pair<size_type, double>>::const_iterator col_val_pair =
- col_val_pairs.begin();
- col_val_pair != col_val_pairs.end();
- ++col_val_pair)
- {
- Assert(line != col_val_pair->first,
- ExcMessage("Can't constrain a degree of freedom to itself"));
-
- for (ConstraintLine::Entries::const_iterator p =
- line_ptr->entries.begin();
- p != line_ptr->entries.end();
- ++p)
- if (p->first == col_val_pair->first)
- {
- // entry exists, break innermost loop
- Assert(
- p->second == col_val_pair->second,
- ExcEntryAlreadyExists(
- line, col_val_pair->first, p->second, col_val_pair->second));
- break;
- }
-
- line_ptr->entries.push_back(*col_val_pair);
- }
-}
-
-
-
-void
-ConstraintMatrix::add_selected_constraints(const ConstraintMatrix &constraints,
- const IndexSet & filter)
-{
- if (constraints.n_constraints() == 0)
- return;
-
- Assert(filter.size() > constraints.lines.back().index,
- ExcMessage("Filter needs to be larger than constraint matrix size."));
- for (std::vector<ConstraintLine>::const_iterator line =
- constraints.lines.begin();
- line != constraints.lines.end();
- ++line)
- if (filter.is_element(line->index))
- {
- const size_type row = filter.index_within_set(line->index);
- add_line(row);
- set_inhomogeneity(row, line->inhomogeneity);
- for (size_type i = 0; i < line->entries.size(); ++i)
- if (filter.is_element(line->entries[i].first))
- add_entry(row,
- filter.index_within_set(line->entries[i].first),
- line->entries[i].second);
- }
-}
-
-
-
-void
-ConstraintMatrix::close()
-{
- if (sorted == true)
- return;
-
- // sort the lines
- std::sort(lines.begin(), lines.end());
-
- // update list of pointers and give the vector a sharp size since we
- // won't modify the size any more after this point.
- {
- std::vector<size_type> new_lines(lines_cache.size(),
- numbers::invalid_size_type);
- size_type counter = 0;
- for (std::vector<ConstraintLine>::const_iterator line = lines.begin();
- line != lines.end();
- ++line, ++counter)
- new_lines[calculate_line_index(line->index)] = counter;
- std::swap(lines_cache, new_lines);
- }
-
- // in debug mode: check whether we really set the pointers correctly.
- for (size_type i = 0; i < lines_cache.size(); ++i)
- if (lines_cache[i] != numbers::invalid_size_type)
- Assert(i == calculate_line_index(lines[lines_cache[i]].index),
- ExcInternalError());
-
- // first, strip zero entries, as we have to do that only once
- for (std::vector<ConstraintLine>::iterator line = lines.begin();
- line != lines.end();
- ++line)
- // first remove zero entries. that would mean that in the linear
- // constraint for a node, x_i = ax_1 + bx_2 + ..., another node times 0
- // appears. obviously, 0*something can be omitted
- line->entries.erase(std::remove_if(line->entries.begin(),
- line->entries.end(),
- &check_zero_weight),
- line->entries.end());
-
-
-
-#ifdef DEBUG
- // In debug mode we are computing an estimate for the maximum number
- // of constraints so that we can bail out if there is a cycle in the
- // constraints (which is easier than searching for cycles in the graph).
- //
- // Let us figure out the largest dof index. This is an upper bound for the
- // number of constraints because it is an approximation for the number of dofs
- // in our system.
- size_type largest_idx = 0;
- for (std::vector<ConstraintLine>::iterator line = lines.begin();
- line != lines.end();
- ++line)
- {
- for (ConstraintLine::Entries::iterator it = line->entries.begin();
- it != line->entries.end();
- ++it)
- {
- largest_idx = std::max(largest_idx, it->first);
- }
- }
-#endif
-
- // replace references to dofs that are themselves constrained. note that
- // because we may replace references to other dofs that may themselves be
- // constrained to third ones, we have to iterate over all this until we
- // replace no chains of constraints any more
- //
- // the iteration replaces references to constrained degrees of freedom by
- // second-order references. for example if x3=x0/2+x2/2 and x2=x0/2+x1/2,
- // then the new list will be x3=x0/2+x0/4+x1/4. note that x0 appear
- // twice. we will throw this duplicate out in the following step, where
- // we sort the list so that throwing out duplicates becomes much more
- // efficient. also, we have to do it only once, rather than in each
- // iteration
- size_type iteration = 0;
- while (true)
- {
- bool chained_constraint_replaced = false;
-
- for (std::vector<ConstraintLine>::iterator line = lines.begin();
- line != lines.end();
- ++line)
- {
-#ifdef DEBUG
- // we need to keep track of how many replacements we do in this line,
- // because we can end up in a cycle A->B->C->A without the number of
- // entries growing.
- size_type n_replacements = 0;
-#endif
-
- // loop over all entries of this line (including ones that we
- // have appended in this go around) and see whether they are
- // further constrained. ignore elements that we don't store on
- // the current processor
- size_type entry = 0;
- while (entry < line->entries.size())
- if (((local_lines.size() == 0) ||
- (local_lines.is_element(line->entries[entry].first))) &&
- is_constrained(line->entries[entry].first))
- {
- // ok, this entry is further constrained:
- chained_constraint_replaced = true;
-
- // look up the chain of constraints for this entry
- const size_type dof_index = line->entries[entry].first;
- const double weight = line->entries[entry].second;
-
- Assert(dof_index != line->index,
- ExcMessage("Cycle in constraints detected!"));
-
- const ConstraintLine *constrained_line =
- &lines[lines_cache[calculate_line_index(dof_index)]];
- Assert(constrained_line->index == dof_index,
- ExcInternalError());
-
- // now we have to replace an entry by its expansion. we do
- // that by overwriting the entry by the first entry of the
- // expansion and adding the remaining ones to the end,
- // where we will later process them once more
- //
- // we can of course only do that if the DoF that we are
- // currently handle is constrained by a linear combination
- // of other dofs:
- if (constrained_line->entries.size() > 0)
- {
- for (size_type i = 0; i < constrained_line->entries.size();
- ++i)
- Assert(dof_index != constrained_line->entries[i].first,
- ExcMessage("Cycle in constraints detected!"));
-
- // replace first entry, then tack the rest to the end
- // of the list
- line->entries[entry] = std::make_pair(
- constrained_line->entries[0].first,
- constrained_line->entries[0].second * weight);
-
- for (size_type i = 1; i < constrained_line->entries.size();
- ++i)
- line->entries.emplace_back(
- constrained_line->entries[i].first,
- constrained_line->entries[i].second * weight);
-
-#ifdef DEBUG
- // keep track of how many entries we replace in this
- // line. If we do more than there are constraints or
- // dofs in our system, we must have a cycle.
- ++n_replacements;
- Assert(n_replacements / 2 < largest_idx,
- ExcMessage("Cycle in constraints detected!"));
- if (n_replacements / 2 >= largest_idx)
- return; // this enables us to test for this Exception.
-#endif
- }
- else
- // the DoF that we encountered is not constrained by a
- // linear combination of other dofs but is equal to just
- // the inhomogeneity (i.e. its chain of entries is
- // empty). in that case, we can't just overwrite the
- // current entry, but we have to actually eliminate it
- {
- line->entries.erase(line->entries.begin() + entry);
- }
-
- line->inhomogeneity += constrained_line->inhomogeneity * weight;
-
- // now that we're here, do not increase index by one but
- // rather make another pass for the present entry because
- // we have replaced the present entry by another one, or
- // because we have deleted it and shifted all following
- // ones one forward
- }
- else
- // entry not further constrained. just move ahead by one
- ++entry;
- }
-
- // if we didn't do anything in this round, then quit the loop
- if (chained_constraint_replaced == false)
- break;
-
- // increase iteration count. note that we should not iterate more
- // times than there are constraints, since this puts a natural upper
- // bound on the length of constraint chains
- ++iteration;
- Assert(iteration <= lines.size(), ExcInternalError());
- }
-
- // finally sort the entries and re-scale them if necessary. in this step,
- // we also throw out duplicates as mentioned above. moreover, as some
- // entries might have had zero weights, we replace them by a vector with
- // sharp sizes.
- for (std::vector<ConstraintLine>::iterator line = lines.begin();
- line != lines.end();
- ++line)
- {
- std::sort(line->entries.begin(), line->entries.end());
-
- // loop over the now sorted list and see whether any of the entries
- // references the same dofs more than once in order to find how many
- // non-duplicate entries we have. This lets us allocate the correct
- // amount of memory for the constraint entries.
- size_type duplicates = 0;
- for (size_type i = 1; i < line->entries.size(); ++i)
- if (line->entries[i].first == line->entries[i - 1].first)
- duplicates++;
-
- if (duplicates > 0 || line->entries.size() < line->entries.capacity())
- {
- ConstraintLine::Entries new_entries;
-
- // if we have no duplicates, copy verbatim the entries. this way,
- // the final size is of the vector is correct.
- if (duplicates == 0)
- new_entries = line->entries;
- else
- {
- // otherwise, we need to go through the list by and and
- // resolve the duplicates
- new_entries.reserve(line->entries.size() - duplicates);
- new_entries.push_back(line->entries[0]);
- for (size_type j = 1; j < line->entries.size(); ++j)
- if (line->entries[j].first == line->entries[j - 1].first)
- {
- Assert(new_entries.back().first == line->entries[j].first,
- ExcInternalError());
- new_entries.back().second += line->entries[j].second;
- }
- else
- new_entries.push_back(line->entries[j]);
-
- Assert(new_entries.size() == line->entries.size() - duplicates,
- ExcInternalError());
-
- // make sure there are really no duplicates left and that the
- // list is still sorted
- for (size_type j = 1; j < new_entries.size(); ++j)
- {
- Assert(new_entries[j].first != new_entries[j - 1].first,
- ExcInternalError());
- Assert(new_entries[j].first > new_entries[j - 1].first,
- ExcInternalError());
- }
- }
-
- // replace old list of constraints for this dof by the new one
- line->entries.swap(new_entries);
- }
-
- // finally do the following check: if the sum of weights for the
- // constraints is close to one, but not exactly one, then rescale all
- // the weights so that they sum up to 1. this adds a little numerical
- // stability and avoids all sorts of problems where the actual value
- // is close to, but not quite what we expected
- //
- // the case where the weights don't quite sum up happens when we
- // compute the interpolation weights "on the fly", i.e. not from
- // precomputed tables. in this case, the interpolation weights are
- // also subject to round-off
- double sum = 0;
- for (size_type i = 0; i < line->entries.size(); ++i)
- sum += line->entries[i].second;
- if ((sum != 1.0) && (std::fabs(sum - 1.) < 1.e-13))
- {
- for (size_type i = 0; i < line->entries.size(); ++i)
- line->entries[i].second /= sum;
- line->inhomogeneity /= sum;
- }
- } // end of loop over all constraint lines
-
-#ifdef DEBUG
- // if in debug mode: check that no dof is constrained to another dof that
- // is also constrained. exclude dofs from this check whose constraint
- // lines are not stored on the local processor
- for (std::vector<ConstraintLine>::const_iterator line = lines.begin();
- line != lines.end();
- ++line)
- for (ConstraintLine::Entries::const_iterator entry = line->entries.begin();
- entry != line->entries.end();
- ++entry)
- if ((local_lines.size() == 0) || (local_lines.is_element(entry->first)))
- {
- // make sure that entry->first is not the index of a line itself
- const bool is_circle = is_constrained(entry->first);
- Assert(is_circle == false,
- ExcDoFConstrainedToConstrainedDoF(line->index, entry->first));
- }
-#endif
-
- sorted = true;
-}
-
-
-
-void
-ConstraintMatrix::merge(const ConstraintMatrix & other_constraints,
- const MergeConflictBehavior merge_conflict_behavior,
- const bool allow_different_local_lines)
-{
- (void)allow_different_local_lines;
- Assert(
- allow_different_local_lines || local_lines == other_constraints.local_lines,
- ExcMessage("local_lines for this and the other objects are not the same "
- "although allow_different_local_lines is false."));
-
- // store the previous state with respect to sorting
- const bool object_was_sorted = sorted;
- sorted = false;
-
- // first action is to fold into the present object possible constraints
- // in the second object. we don't strictly need to do this any more since
- // the ConstraintMatrix has learned to deal with chains of constraints in
- // the close() function, but we have traditionally done this and it's not
- // overly hard to do.
- //
- // for this, loop over all constraints and replace the constraint lines
- // with a new one where constraints are replaced if necessary.
- ConstraintLine::Entries tmp;
- for (std::vector<ConstraintLine>::iterator line = lines.begin();
- line != lines.end();
- ++line)
- {
- tmp.clear();
- for (size_type i = 0; i < line->entries.size(); ++i)
- {
- // if the present dof is not stored, or not constrained, or if we
- // won't take the constraint from the other object, then simply copy
- // it over
- if ((other_constraints.local_lines.size() != 0 &&
- other_constraints.local_lines.is_element(
- line->entries[i].first) == false) ||
- other_constraints.is_constrained(line->entries[i].first) ==
- false ||
- ((merge_conflict_behavior != right_object_wins) &&
- other_constraints.is_constrained(line->entries[i].first) &&
- this->is_constrained(line->entries[i].first)))
- tmp.push_back(line->entries[i]);
- else
- // otherwise resolve further constraints by replacing the old
- // entry by a sequence of new entries taken from the other
- // object, but with multiplied weights
- {
- const ConstraintLine::Entries *other_line =
- other_constraints.get_constraint_entries(
- line->entries[i].first);
- Assert(other_line != nullptr, ExcInternalError());
-
- const double weight = line->entries[i].second;
-
- for (ConstraintLine::Entries::const_iterator j =
- other_line->begin();
- j != other_line->end();
- ++j)
- tmp.emplace_back(j->first, j->second * weight);
-
- line->inhomogeneity +=
- other_constraints.get_inhomogeneity(line->entries[i].first) *
- weight;
- }
- }
- // finally exchange old and newly resolved line
- line->entries.swap(tmp);
- }
-
- if (local_lines.size() != 0)
- local_lines.add_indices(other_constraints.local_lines);
-
- {
- // do not bother to resize the lines cache exactly since it is pretty
- // cheap to adjust it along the way.
- std::fill(
- lines_cache.begin(), lines_cache.end(), numbers::invalid_size_type);
-
- // reset lines_cache for our own constraints
- size_type index = 0;
- for (std::vector<ConstraintLine>::const_iterator line = lines.begin();
- line != lines.end();
- ++line)
- {
- size_type local_line_no = calculate_line_index(line->index);
- if (local_line_no >= lines_cache.size())
- lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
- lines_cache[local_line_no] = index++;
- }
-
- // Add other_constraints to lines cache and our list of constraints
- for (std::vector<ConstraintLine>::const_iterator line =
- other_constraints.lines.begin();
- line != other_constraints.lines.end();
- ++line)
- {
- const size_type local_line_no = calculate_line_index(line->index);
- if (local_line_no >= lines_cache.size())
- {
- lines_cache.resize(local_line_no + 1, numbers::invalid_size_type);
- lines.push_back(*line);
- lines_cache[local_line_no] = index++;
- }
- else if (lines_cache[local_line_no] == numbers::invalid_size_type)
- {
- // there are no constraints for that line yet
- lines.push_back(*line);
- AssertIndexRange(local_line_no, lines_cache.size());
- lines_cache[local_line_no] = index++;
- }
- else
- {
- // we already store that line
- switch (merge_conflict_behavior)
- {
- case no_conflicts_allowed:
- AssertThrow(false,
- ExcDoFIsConstrainedFromBothObjects(line->index));
- break;
-
- case left_object_wins:
- // ignore this constraint
- break;
-
- case right_object_wins:
- AssertIndexRange(local_line_no, lines_cache.size());
- lines[lines_cache[local_line_no]] = *line;
- break;
-
- default:
- Assert(false, ExcNotImplemented());
- }
- }
- }
-
- // check that we set the pointers correctly
- for (size_type i = 0; i < lines_cache.size(); ++i)
- if (lines_cache[i] != numbers::invalid_size_type)
- Assert(i == calculate_line_index(lines[lines_cache[i]].index),
- ExcInternalError());
- }
-
- // if the object was sorted before, then make sure it is so afterward as
- // well. otherwise leave everything in the unsorted state
- if (object_was_sorted == true)
- close();
-}
-
-
-
-void
-ConstraintMatrix::shift(const size_type offset)
-{
- if (local_lines.size() == 0)
- lines_cache.insert(lines_cache.begin(), offset, numbers::invalid_size_type);
- else
- {
- // shift local_lines
- IndexSet new_local_lines(local_lines.size());
- new_local_lines.add_indices(local_lines, offset);
- std::swap(local_lines, new_local_lines);
- }
-
- for (std::vector<ConstraintLine>::iterator i = lines.begin();
- i != lines.end();
- ++i)
- {
- i->index += offset;
- for (ConstraintLine::Entries::iterator j = i->entries.begin();
- j != i->entries.end();
- ++j)
- j->first += offset;
- }
-
-#ifdef DEBUG
- // make sure that lines, lines_cache and local_lines
- // are still linked correctly
- for (size_type i = 0; i < lines_cache.size(); ++i)
- Assert(lines_cache[i] == numbers::invalid_size_type ||
- calculate_line_index(lines[lines_cache[i]].index) == i,
- ExcInternalError());
-#endif
-}
-
-
-
-void
-ConstraintMatrix::clear()
-{
- {
- std::vector<ConstraintLine> tmp;
- lines.swap(tmp);
- }
-
- {
- std::vector<size_type> tmp;
- lines_cache.swap(tmp);
- }
-
- sorted = false;
-}
-
-
-
-void
-ConstraintMatrix::reinit(const IndexSet &local_constraints)
-{
- local_lines = local_constraints;
-
- // make sure the IndexSet is compressed. Otherwise this can lead to crashes
- // that are hard to find (only happen in release mode).
- // see tests/mpi/constraint_matrix_crash_01
- local_lines.compress();
-
- clear();
-}
-
-
-
-void
-ConstraintMatrix::condense(SparsityPattern &sparsity) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.is_compressed() == false, ExcMatrixIsClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
-
- // store for each index whether it must be distributed or not. If entry
- // is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
- // handles this index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = c;
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- if (distribute[row] == numbers::invalid_size_type)
- {
- // regular line. loop over cols all valid cols. note that this
- // changes the line we are presently working on: we add additional
- // entries. these are put to the end of the row. however, as
- // constrained nodes cannot be constrained to other constrained
- // nodes, nothing will happen if we run into these added nodes, as
- // they can't be distributed further. we might store the position of
- // the last old entry and stop work there, but since operating on
- // the newly added ones only takes two comparisons (column index
- // valid, distribute[column] necessarily
- // ==numbers::invalid_size_type), it is cheaper to not do so and
- // run right until the end of the line
- for (SparsityPattern::iterator entry = sparsity.begin(row);
- ((entry != sparsity.end(row)) && entry->is_valid_entry());
- ++entry)
- {
- const size_type column = entry->column();
-
- if (distribute[column] != numbers::invalid_size_type)
- {
- // distribute entry at regular row @p{row} and irregular
- // column sparsity.colnums[j]
- for (size_type q = 0;
- q != lines[distribute[column]].entries.size();
- ++q)
- sparsity.add(row,
- lines[distribute[column]].entries[q].first);
- }
- }
- }
- else
- // row must be distributed. note that here the present row is not
- // touched (unlike above)
- {
- for (SparsityPattern::iterator entry = sparsity.begin(row);
- (entry != sparsity.end(row)) && entry->is_valid_entry();
- ++entry)
- {
- const size_type column = entry->column();
- if (distribute[column] == numbers::invalid_size_type)
- // distribute entry at irregular row @p{row} and regular
- // column sparsity.colnums[j]
- for (size_type q = 0;
- q != lines[distribute[row]].entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[q].first, column);
- else
- // distribute entry at irregular row @p{row} and irregular
- // column sparsity.get_column_numbers()[j]
- for (size_type p = 0;
- p != lines[distribute[row]].entries.size();
- ++p)
- for (size_type q = 0;
- q != lines[distribute[column]].entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[p].first,
- lines[distribute[column]].entries[q].first);
- }
- }
- }
-
- sparsity.compress();
-}
-
-
-
-void
-ConstraintMatrix::condense(DynamicSparsityPattern &sparsity) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
-
- // store for each index whether it must be distributed or not. If entry
- // is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
- // handles this index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = c;
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- if (distribute[row] == numbers::invalid_size_type)
- // regular line. loop over cols. note that as we proceed to
- // distribute cols, the loop may get longer
- for (size_type j = 0; j < sparsity.row_length(row); ++j)
- {
- const size_type column = sparsity.column_number(row, j);
-
- if (distribute[column] != numbers::invalid_size_type)
- {
- // distribute entry at regular row @p{row} and irregular
- // column column. note that this changes the line we are
- // presently working on: we add additional entries. if we
- // add another entry at a column behind the present one, we
- // will encounter it later on (but since it can't be
- // further constrained, won't have to do anything about
- // it). if we add it up front of the present column, we
- // will find the present column later on again as it was
- // shifted back (again nothing happens, in particular no
- // endless loop, as when we encounter it the second time we
- // won't be able to add more entries as they all already
- // exist, but we do the same work more often than
- // necessary, and the loop gets longer), so move the cursor
- // one to the right in the case that we add an entry up
- // front that did not exist before. check whether it
- // existed before by tracking the length of this row
- size_type old_rowlength = sparsity.row_length(row);
- for (size_type q = 0;
- q != lines[distribute[column]].entries.size();
- ++q)
- {
- const size_type new_col =
- lines[distribute[column]].entries[q].first;
-
- sparsity.add(row, new_col);
-
- const size_type new_rowlength = sparsity.row_length(row);
- if ((new_col < column) && (old_rowlength != new_rowlength))
- ++j;
- old_rowlength = new_rowlength;
- }
- }
- }
- else
- // row must be distributed
- for (size_type j = 0; j < sparsity.row_length(row); ++j)
- {
- const size_type column = sparsity.column_number(row, j);
-
- if (distribute[column] == numbers::invalid_size_type)
- // distribute entry at irregular row @p{row} and regular
- // column sparsity.colnums[j]
- for (size_type q = 0; q != lines[distribute[row]].entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[q].first, column);
- else
- // distribute entry at irregular row @p{row} and irregular
- // column sparsity.get_column_numbers()[j]
- for (size_type p = 0; p != lines[distribute[row]].entries.size();
- ++p)
- for (size_type q = 0;
- q != lines[distribute[sparsity.column_number(row, j)]]
- .entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[p].first,
- lines[distribute[sparsity.column_number(row, j)]]
- .entries[q]
- .first);
- }
- }
-}
-
-
-
-void
-ConstraintMatrix::condense(BlockSparsityPattern &sparsity) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.is_compressed() == false, ExcMatrixIsClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
- Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
- Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
- ExcNotQuadratic());
-
- const BlockIndices &index_mapping = sparsity.get_column_indices();
-
- const size_type n_blocks = sparsity.n_block_rows();
-
- // store for each index whether it must be distributed or not. If entry
- // is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
- // handles this index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = c;
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- // get index of this row within the blocks
- const std::pair<size_type, size_type> block_index =
- index_mapping.global_to_local(row);
- const size_type block_row = block_index.first;
-
- if (distribute[row] == numbers::invalid_size_type)
- // regular line. loop over all columns and see whether this column
- // must be distributed
- {
- // to loop over all entries in this row, we have to loop over all
- // blocks in this blockrow and the corresponding row therein
- for (size_type block_col = 0; block_col < n_blocks; ++block_col)
- {
- const SparsityPattern &block_sparsity =
- sparsity.block(block_row, block_col);
-
- for (SparsityPattern::const_iterator entry =
- block_sparsity.begin(block_index.second);
- (entry != block_sparsity.end(block_index.second)) &&
- entry->is_valid_entry();
- ++entry)
- {
- const size_type global_col =
- index_mapping.local_to_global(block_col, entry->column());
-
- if (distribute[global_col] != numbers::invalid_size_type)
- // distribute entry at regular row @p{row} and
- // irregular column global_col
- {
- for (size_type q = 0;
- q != lines[distribute[global_col]].entries.size();
- ++q)
- sparsity.add(
- row, lines[distribute[global_col]].entries[q].first);
- }
- }
- }
- }
- else
- {
- // row must be distributed. split the whole row into the chunks
- // defined by the blocks
- for (size_type block_col = 0; block_col < n_blocks; ++block_col)
- {
- const SparsityPattern &block_sparsity =
- sparsity.block(block_row, block_col);
-
- for (SparsityPattern::const_iterator entry =
- block_sparsity.begin(block_index.second);
- (entry != block_sparsity.end(block_index.second)) &&
- entry->is_valid_entry();
- ++entry)
- {
- const size_type global_col =
- index_mapping.local_to_global(block_col, entry->column());
-
- if (distribute[global_col] == numbers::invalid_size_type)
- // distribute entry at irregular row @p{row} and
- // regular column global_col.
- {
- for (size_type q = 0;
- q != lines[distribute[row]].entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[q].first,
- global_col);
- }
- else
- // distribute entry at irregular row @p{row} and
- // irregular column @p{global_col}
- {
- for (size_type p = 0;
- p != lines[distribute[row]].entries.size();
- ++p)
- for (size_type q = 0;
- q != lines[distribute[global_col]].entries.size();
- ++q)
- sparsity.add(
- lines[distribute[row]].entries[p].first,
- lines[distribute[global_col]].entries[q].first);
- }
- }
- }
- }
- }
-
- sparsity.compress();
-}
-
-
-
-void
-ConstraintMatrix::condense(BlockDynamicSparsityPattern &sparsity) const
-{
- Assert(sorted == true, ExcMatrixNotClosed());
- Assert(sparsity.n_rows() == sparsity.n_cols(), ExcNotQuadratic());
- Assert(sparsity.n_block_rows() == sparsity.n_block_cols(), ExcNotQuadratic());
- Assert(sparsity.get_column_indices() == sparsity.get_row_indices(),
- ExcNotQuadratic());
-
- const BlockIndices &index_mapping = sparsity.get_column_indices();
-
- const size_type n_blocks = sparsity.n_block_rows();
-
- // store for each index whether it must be distributed or not. If entry
- // is numbers::invalid_unsigned_int, no distribution is necessary.
- // otherwise, the number states which line in the constraint matrix
- // handles this index
- std::vector<size_type> distribute(sparsity.n_rows(),
- numbers::invalid_size_type);
-
- for (size_type c = 0; c < lines.size(); ++c)
- distribute[lines[c].index] = static_cast<signed int>(c);
-
- const size_type n_rows = sparsity.n_rows();
- for (size_type row = 0; row < n_rows; ++row)
- {
- // get index of this row within the blocks
- const std::pair<size_type, size_type> block_index =
- index_mapping.global_to_local(row);
- const size_type block_row = block_index.first;
- const size_type local_row = block_index.second;
-
- if (distribute[row] == numbers::invalid_size_type)
- // regular line. loop over all columns and see whether this column
- // must be distributed. note that as we proceed to distribute cols,
- // the loop over cols may get longer.
- //
- // don't try to be clever here as in the algorithm for the
- // DynamicSparsityPattern, as that would be much more
- // complicated here. after all, we know that compressed patterns
- // are inefficient...
- {
- // to loop over all entries in this row, we have to loop over all
- // blocks in this blockrow and the corresponding row therein
- for (size_type block_col = 0; block_col < n_blocks; ++block_col)
- {
- const DynamicSparsityPattern &block_sparsity =
- sparsity.block(block_row, block_col);
-
- for (size_type j = 0; j < block_sparsity.row_length(local_row);
- ++j)
- {
- const size_type global_col = index_mapping.local_to_global(
- block_col, block_sparsity.column_number(local_row, j));
-
- if (distribute[global_col] != numbers::invalid_size_type)
- // distribute entry at regular row @p{row} and
- // irregular column global_col
- {
- for (size_type q = 0;
- q != lines[distribute[global_col]].entries.size();
- ++q)
- sparsity.add(
- row, lines[distribute[global_col]].entries[q].first);
- }
- }
- }
- }
- else
- {
- // row must be distributed. split the whole row into the chunks
- // defined by the blocks
- for (size_type block_col = 0; block_col < n_blocks; ++block_col)
- {
- const DynamicSparsityPattern &block_sparsity =
- sparsity.block(block_row, block_col);
-
- for (size_type j = 0; j < block_sparsity.row_length(local_row);
- ++j)
- {
- const size_type global_col = index_mapping.local_to_global(
- block_col, block_sparsity.column_number(local_row, j));
-
- if (distribute[global_col] == numbers::invalid_size_type)
- // distribute entry at irregular row @p{row} and
- // regular column global_col.
- {
- for (size_type q = 0;
- q != lines[distribute[row]].entries.size();
- ++q)
- sparsity.add(lines[distribute[row]].entries[q].first,
- global_col);
- }
- else
- // distribute entry at irregular row @p{row} and
- // irregular column @p{global_col}
- {
- for (size_type p = 0;
- p != lines[distribute[row]].entries.size();
- ++p)
- for (size_type q = 0;
- q != lines[distribute[global_col]].entries.size();
- ++q)
- sparsity.add(
- lines[distribute[row]].entries[p].first,
- lines[distribute[global_col]].entries[q].first);
- }
- }
- }
- }
- }
-}
-
-
-
-bool
-ConstraintMatrix::is_identity_constrained(const size_type index) const
-{
- if (is_constrained(index) == false)
- return false;
-
- const ConstraintLine &p = lines[lines_cache[calculate_line_index(index)]];
- Assert(p.index == index, ExcInternalError());
-
- // return if an entry for this line was found and if it has only one
- // entry equal to 1.0
- return ((p.entries.size() == 1) && (p.entries[0].second == 1.0));
-}
-
-
-bool
-ConstraintMatrix::are_identity_constrained(const size_type index1,
- const size_type index2) const
-{
- if (is_constrained(index1) == true)
- {
- const ConstraintLine &p =
- lines[lines_cache[calculate_line_index(index1)]];
- Assert(p.index == index1, ExcInternalError());
-
- // return if an entry for this line was found and if it has only one
- // entry equal to 1.0 and that one is index2
- return ((p.entries.size() == 1) && (p.entries[0].first == index2) &&
- (p.entries[0].second == 1.0));
- }
- else if (is_constrained(index2) == true)
- {
- const ConstraintLine &p =
- lines[lines_cache[calculate_line_index(index2)]];
- Assert(p.index == index2, ExcInternalError());
-
- // return if an entry for this line was found and if it has only one
- // entry equal to 1.0 and that one is index1
- return ((p.entries.size() == 1) && (p.entries[0].first == index1) &&
- (p.entries[0].second == 1.0));
- }
- else
- return false;
-}
-
-
-
-ConstraintMatrix::size_type
-ConstraintMatrix::max_constraint_indirections() const
-{
- size_type return_value = 0;
- for (std::vector<ConstraintLine>::const_iterator i = lines.begin();
- i != lines.end();
- ++i)
- // use static cast, since typeof(size)==std::size_t, which is !=
- // size_type on AIX
- return_value =
- std::max(return_value, static_cast<size_type>(i->entries.size()));
-
- return return_value;
-}
-
-
-
-bool
-ConstraintMatrix::has_inhomogeneities() const
-{
- for (std::vector<ConstraintLine>::const_iterator i = lines.begin();
- i != lines.end();
- ++i)
- if (i->inhomogeneity != 0.)
- return true;
-
- return false;
-}
-
-
-void
-ConstraintMatrix::print(std::ostream &out) const
-{
- for (size_type i = 0; i != lines.size(); ++i)
- {
- // output the list of constraints as pairs of dofs and their weights
- if (lines[i].entries.size() > 0)
- {
- for (size_type j = 0; j < lines[i].entries.size(); ++j)
- out << " " << lines[i].index << " " << lines[i].entries[j].first
- << ": " << lines[i].entries[j].second << "\n";
-
- // print out inhomogeneity.
- if (lines[i].inhomogeneity != 0)
- out << " " << lines[i].index << ": " << lines[i].inhomogeneity
- << "\n";
- }
- else
- // but also output something if the constraint simply reads
- // x[13]=0, i.e. where the right hand side is not a linear
- // combination of other dofs
- {
- if (lines[i].inhomogeneity != 0)
- out << " " << lines[i].index << " = " << lines[i].inhomogeneity
- << "\n";
- else
- out << " " << lines[i].index << " = 0\n";
- }
- }
-
- AssertThrow(out, ExcIO());
-}
-
-
-
-void
-ConstraintMatrix::write_dot(std::ostream &out) const
-{
- out << "digraph constraints {" << std::endl;
- for (size_type i = 0; i != lines.size(); ++i)
- {
- // same concept as in the previous function
- if (lines[i].entries.size() > 0)
- for (size_type j = 0; j < lines[i].entries.size(); ++j)
- out << " " << lines[i].index << "->" << lines[i].entries[j].first
- << "; // weight: " << lines[i].entries[j].second << "\n";
- else
- out << " " << lines[i].index << "\n";
- }
- out << "}" << std::endl;
-}
-
-
-
-std::size_t
-ConstraintMatrix::memory_consumption() const
-{
- return (MemoryConsumption::memory_consumption(lines) +
- MemoryConsumption::memory_consumption(lines_cache) +
- MemoryConsumption::memory_consumption(sorted) +
- MemoryConsumption::memory_consumption(local_lines));
-}
-
-
-
-void
-ConstraintMatrix::resolve_indices(
- std::vector<types::global_dof_index> &indices) const
-{
- const unsigned int indices_size = indices.size();
- const std::vector<std::pair<types::global_dof_index, double>> *line_ptr;
- for (unsigned int i = 0; i < indices_size; ++i)
- {
- line_ptr = get_constraint_entries(indices[i]);
- // if the index is constraint, the constraints indices are added to the
- // indices vector
- if (line_ptr != nullptr)
- {
- const unsigned int line_size = line_ptr->size();
- for (unsigned int j = 0; j < line_size; ++j)
- indices.push_back((*line_ptr)[j].first);
- }
- }
-
- // keep only the unique elements
- std::sort(indices.begin(), indices.end());
- std::vector<types::global_dof_index>::iterator it;
- it = std::unique(indices.begin(), indices.end());
- indices.resize(it - indices.begin());
-}
-
-
-
-// explicit instantiations
-//
-// define a list of functions for vectors and matrices, respectively, where
-// the vector/matrix can be replaced using a preprocessor variable
-// VectorType/MatrixType. note that we need a space between "VectorType" and
-// ">" to disambiguate ">>" when VectorType trails in an angle bracket
-
-// TODO: The way we define all the instantiations is probably not the very
-// best one. Try to find a better description.
-
-#define VECTOR_FUNCTIONS(VectorType) \
- template void ConstraintMatrix::condense<VectorType>( \
- const VectorType &uncondensed, VectorType &condensed) const; \
- template void ConstraintMatrix::condense<VectorType>(VectorType & vec) \
- const; \
- template void ConstraintMatrix::distribute_local_to_global<VectorType>( \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- VectorType &, \
- const FullMatrix<VectorType::value_type> &) const; \
- template void ConstraintMatrix::distribute_local_to_global<VectorType>( \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- VectorType &, \
- const FullMatrix<VectorType::value_type> &, \
- bool) const
-
-#define PARALLEL_VECTOR_FUNCTIONS(VectorType) \
- template void ConstraintMatrix::distribute_local_to_global<VectorType>( \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- VectorType &, \
- const FullMatrix<VectorType::value_type> &) const; \
- template void ConstraintMatrix::distribute_local_to_global<VectorType>( \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- VectorType &, \
- const FullMatrix<VectorType::value_type> &, \
- bool) const
-
-#ifdef DEAL_II_WITH_PETSC
-VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector);
-VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockVector);
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
-PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::Vector);
-PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector);
-#endif
-
-#define MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \
- template void \
- ConstraintMatrix::distribute_local_to_global<MatrixType, VectorType>( \
- const FullMatrix<MatrixType::value_type> &, \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &, \
- VectorType &, \
- bool, \
- std::integral_constant<bool, false>) const
-#define MATRIX_FUNCTIONS(MatrixType, VectorScalar) \
- template void \
- ConstraintMatrix::distribute_local_to_global<MatrixType, \
- Vector<VectorScalar>>( \
- const FullMatrix<MatrixType::value_type> &, \
- const Vector<VectorScalar> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &, \
- Vector<VectorScalar> &, \
- bool, \
- std::integral_constant<bool, false>) const
-#define BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \
- template void \
- ConstraintMatrix::distribute_local_to_global<MatrixType, VectorType>( \
- const FullMatrix<MatrixType::value_type> &, \
- const Vector<VectorType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &, \
- VectorType &, \
- bool, \
- std::integral_constant<bool, true>) const
-#define BLOCK_MATRIX_FUNCTIONS(MatrixType) \
- template void ConstraintMatrix:: \
- distribute_local_to_global<MatrixType, Vector<MatrixType::value_type>>( \
- const FullMatrix<MatrixType::value_type> &, \
- const Vector<MatrixType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &, \
- Vector<MatrixType::value_type> &, \
- bool, \
- std::integral_constant<bool, true>) const
-
-MATRIX_FUNCTIONS(FullMatrix<double>, double);
-MATRIX_FUNCTIONS(FullMatrix<float>, float);
-MATRIX_FUNCTIONS(FullMatrix<double>, std::complex<double>);
-MATRIX_FUNCTIONS(FullMatrix<std::complex<double>>, std::complex<double>);
-
-MATRIX_FUNCTIONS(SparseMatrix<double>, double);
-MATRIX_FUNCTIONS(SparseMatrix<float>, float);
-MATRIX_FUNCTIONS(SparseMatrix<double>, std::complex<double>);
-MATRIX_FUNCTIONS(SparseMatrix<float>, std::complex<float>);
-MATRIX_FUNCTIONS(SparseMatrix<std::complex<double>>, std::complex<double>);
-MATRIX_FUNCTIONS(SparseMatrix<std::complex<float>>, std::complex<float>);
-
-MATRIX_FUNCTIONS(SparseMatrixEZ<double>, double);
-MATRIX_FUNCTIONS(SparseMatrixEZ<float>, float);
-MATRIX_FUNCTIONS(ChunkSparseMatrix<double>, double);
-MATRIX_FUNCTIONS(ChunkSparseMatrix<float>, float);
-
-
-BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
-BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
-BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<double>, BlockVector<double>);
-BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix<float>, BlockVector<float>);
-
-// BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ<double>);
-// BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ<float>, Vector<float>);
-
-#ifdef DEAL_II_WITH_PETSC
-MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix, PetscScalar);
-MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PetscScalar);
-BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
-MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix,
- PETScWrappers::MPI::Vector);
-MATRIX_VECTOR_FUNCTIONS(PETScWrappers::SparseMatrix,
- PETScWrappers::MPI::Vector);
-BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix,
- PETScWrappers::MPI::BlockVector);
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
-MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix, double);
-BLOCK_MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix);
-MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::SparseMatrix,
- TrilinosWrappers::MPI::Vector);
-BLOCK_MATRIX_VECTOR_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix,
- TrilinosWrappers::MPI::BlockVector);
-#endif
-
-
-#define SPARSITY_FUNCTIONS(SparsityPatternType) \
- template void \
- ConstraintMatrix::add_entries_local_to_global<SparsityPatternType>( \
- const std::vector<ConstraintMatrix::size_type> &, \
- SparsityPatternType &, \
- const bool, \
- const Table<2, bool> &, \
- std::integral_constant<bool, false>) const; \
- template void \
- ConstraintMatrix::add_entries_local_to_global<SparsityPatternType>( \
- const std::vector<ConstraintMatrix::size_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- SparsityPatternType &, \
- const bool, \
- const Table<2, bool> &) const
-#define BLOCK_SPARSITY_FUNCTIONS(SparsityPatternType) \
- template void \
- ConstraintMatrix::add_entries_local_to_global<SparsityPatternType>( \
- const std::vector<ConstraintMatrix::size_type> &, \
- SparsityPatternType &, \
- const bool, \
- const Table<2, bool> &, \
- std::integral_constant<bool, true>) const; \
- template void \
- ConstraintMatrix::add_entries_local_to_global<SparsityPatternType>( \
- const std::vector<ConstraintMatrix::size_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- SparsityPatternType &, \
- const bool, \
- const Table<2, bool> &) const
-
-SPARSITY_FUNCTIONS(SparsityPattern);
-SPARSITY_FUNCTIONS(DynamicSparsityPattern);
-BLOCK_SPARSITY_FUNCTIONS(BlockSparsityPattern);
-BLOCK_SPARSITY_FUNCTIONS(BlockDynamicSparsityPattern);
-
-#ifdef DEAL_II_WITH_TRILINOS
-SPARSITY_FUNCTIONS(TrilinosWrappers::SparsityPattern);
-BLOCK_SPARSITY_FUNCTIONS(TrilinosWrappers::BlockSparsityPattern);
-#endif
-
-
-#define ONLY_MATRIX_FUNCTIONS(MatrixType) \
- template void ConstraintMatrix::distribute_local_to_global<MatrixType>( \
- const FullMatrix<MatrixType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &) const; \
- template void ConstraintMatrix::distribute_local_to_global<MatrixType>( \
- const FullMatrix<MatrixType::value_type> &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- const ConstraintMatrix &, \
- const std::vector<ConstraintMatrix::size_type> &, \
- MatrixType &) const
-
-ONLY_MATRIX_FUNCTIONS(FullMatrix<float>);
-ONLY_MATRIX_FUNCTIONS(FullMatrix<double>);
-ONLY_MATRIX_FUNCTIONS(SparseMatrix<float>);
-ONLY_MATRIX_FUNCTIONS(SparseMatrix<double>);
-ONLY_MATRIX_FUNCTIONS(MatrixBlock<SparseMatrix<float>>);
-ONLY_MATRIX_FUNCTIONS(MatrixBlock<SparseMatrix<double>>);
-ONLY_MATRIX_FUNCTIONS(BlockSparseMatrix<float>);
-ONLY_MATRIX_FUNCTIONS(BlockSparseMatrix<double>);
-
-#ifdef DEAL_II_WITH_TRILINOS
-ONLY_MATRIX_FUNCTIONS(TrilinosWrappers::SparseMatrix);
-ONLY_MATRIX_FUNCTIONS(TrilinosWrappers::BlockSparseMatrix);
-#endif
-
-#ifdef DEAL_II_WITH_PETSC
-ONLY_MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix);
-ONLY_MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix);
-ONLY_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix);
-#endif
-
-#include "constraint_matrix.inst"
-
-DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2013 - 2017 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE at
-// the top level of the deal.II distribution.
-//
-// ---------------------------------------------------------------------
-
-for (S : REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES)
- {
- template void ConstraintMatrix::condense<T<S>>(const T<S> &, T<S> &) const;
- template void ConstraintMatrix::condense<T<S>>(T<S> & vec) const;
- template void ConstraintMatrix::distribute_local_to_global<T<S>>(
- const Vector<S> &,
- const std::vector<types::global_dof_index> &,
- T<S> &,
- const FullMatrix<S> &) const;
- template void ConstraintMatrix::distribute_local_to_global<T<S>>(
- const Vector<S> &,
- const std::vector<types::global_dof_index> &,
- const std::vector<types::global_dof_index> &,
- T<S> &,
- const FullMatrix<S> &,
- bool) const;
- template void ConstraintMatrix::set_zero<T<S>>(T<S> &) const;
- }
-
-
-for (S : REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES)
- {
- template void ConstraintMatrix::condense<LinearAlgebra::distributed::T<S>>(
- const LinearAlgebra::distributed::T<S> &,
- LinearAlgebra::distributed::T<S> &) const;
- template void ConstraintMatrix::condense<LinearAlgebra::distributed::T<S>>(
- LinearAlgebra::distributed::T<S> & vec) const;
-
- template void ConstraintMatrix::distribute_local_to_global<
- LinearAlgebra::distributed::T<S>>(
- const Vector<S> &,
- const std::vector<types::global_dof_index> &,
- LinearAlgebra::distributed::T<S> &,
- const FullMatrix<S> &) const;
-
- template void ConstraintMatrix::distribute_local_to_global<
- LinearAlgebra::distributed::T<S>>(
- const Vector<S> &,
- const std::vector<types::global_dof_index> &,
- const std::vector<types::global_dof_index> &,
- LinearAlgebra::distributed::T<S> &,
- const FullMatrix<S> &,
- bool) const;
-
- template void ConstraintMatrix::distribute_local_to_global<
- DiagonalMatrix<LinearAlgebra::distributed::T<S>>>(
- const FullMatrix<S> &,
- const std::vector<size_type> &,
- DiagonalMatrix<LinearAlgebra::distributed::T<S>> &) const;
-
- template void ConstraintMatrix::distribute_local_to_global<
- DiagonalMatrix<LinearAlgebra::distributed::T<S>>,
- LinearAlgebra::distributed::T<S>>(
- const FullMatrix<S> &,
- const Vector<S> &,
- const std::vector<size_type> &,
- DiagonalMatrix<LinearAlgebra::distributed::T<S>> &,
- LinearAlgebra::distributed::T<S> &,
- bool,
- std::integral_constant<bool, false>) const;
-
- template void ConstraintMatrix::distribute_local_to_global<
- DiagonalMatrix<LinearAlgebra::distributed::T<S>>,
- T<S>>(const FullMatrix<S> &,
- const Vector<S> &,
- const std::vector<size_type> &,
- DiagonalMatrix<LinearAlgebra::distributed::T<S>> &,
- T<S> &,
- bool,
- std::integral_constant<bool, false>) const;
-
- template void ConstraintMatrix::set_zero<LinearAlgebra::distributed::T<S>>(
- LinearAlgebra::distributed::T<S> &) const;
- }
-
-
-for (V : EXTERNAL_PARALLEL_VECTORS)
- {
- template void ConstraintMatrix::set_zero<V>(V &) const;
- }
-
-
-for (S : REAL_SCALARS)
- {
- template void ConstraintMatrix::condense<S>(SparseMatrix<S> &) const;
- template void ConstraintMatrix::condense<S>(BlockSparseMatrix<S> &) const;
- }
-
-
-for (S1 : REAL_SCALARS; S2 : REAL_SCALARS)
- {
- template void ConstraintMatrix::condense<S1, Vector<S2>>(
- SparseMatrix<S1> &, Vector<S2> &) const;
- template void ConstraintMatrix::condense<S1, BlockVector<S2>>(
- BlockSparseMatrix<S1> &, BlockVector<S2> &) const;
- }
-
-for (S1 : COMPLEX_SCALARS)
- {
- template void ConstraintMatrix::condense<S1, Vector<S1>>(
- SparseMatrix<S1> &, Vector<S1> &) const;
- }
-
-
-for (Vec : VECTOR_TYPES)
- {
- template void ConstraintMatrix::distribute<Vec>(Vec &) const;
- }