* the classes that are used to wrap
* PETSc objects.
*
- * For a replacement function,
- * see the documentation of the
- * @ref{FilteredMatrix} class in
- * the @p{LAC} sublibrary.
+ * Note that this function is not very
+ * efficient: it needs to alternatingly
+ * read and write into the matrix, a
+ * situation that PETSc does not handle
+ * too well. In addition, we only get rid
+ * of rows corresponding to boundary
+ * nodes, but the corresponding case of
+ * deleting the respective columns
+ * (i.e. if @arg eliminate_columns is @p
+ * true) is not presently implemented,
+ * and probably will never because it is
+ * too expensive without direct access to
+ * the PETSc data structures. A third
+ * reason against this function is that
+ * it doesn't handle the case where the
+ * matrix is distributed across an MPI
+ * system.
+ *
+ * In order to still be able to eliminate
+ * boundary values, it is better to get
+ * rid of them before the local matrices
+ * and vectors are distributed to the
+ * global ones, because then we don't
+ * have to mess with the sparse data
+ * structures. The
+ * local_apply_boundary_values() function
+ * does that, and is recommended for use
+ * instead of the global one for PETSc
+ * matrices and vectors.
*/
#ifdef DEAL_II_USE_PETSC
static void
PETScWrappers::VectorBase &right_hand_side,
const bool eliminate_columns = true);
#endif
+
+ /**
+ * Rather than applying boundary values
+ * to the global matrix and vector after
+ * assembly, this function does so @em
+ * before assembling, by modifying the
+ * local matrix and vector
+ * contributions. If you call this
+ * function on all local contributions,
+ * the resulting matrix will have the
+ * same entries, and the final call to
+ * apply_boundary_values() on the global
+ * system will not be necessary.
+ *
+ * Since this function does not have to
+ * work on the complicated data
+ * structures of sparse matrices, it is
+ * relatively cheap. It may therefore be
+ * a win if you have many fixed degrees
+ * of freedom (e.g. boundary nodes), or
+ * if access to the sparse matrix is
+ * expensive (e.g. for block sparse
+ * matrices, or for PETSc matrices).
+ */
+ static void
+ local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ const std::vector<unsigned int> &local_dof_indices,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_rhs,
+ const bool eliminate_columns);
/**
* Exception
#include <numerics/matrices.h>
+#include <lac/full_matrix.h>
#include <lac/vector.h>
#include <lac/block_vector.h>
#include <lac/sparse_matrix.h>
#endif
+void
+MatrixTools::
+local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+ const std::vector<unsigned int> &local_dof_indices,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_rhs,
+ const bool eliminate_columns)
+{
+ Assert (local_dof_indices.size() == local_matrix.m(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_matrix.m()));
+ Assert (local_dof_indices.size() == local_matrix.n(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_matrix.n()));
+ Assert (local_dof_indices.size() == local_rhs.size(),
+ ExcDimensionMismatch(local_dof_indices.size(),
+ local_rhs.size()));
+
+ // if there is nothing to do, then exit
+ // right away
+ if (boundary_values.size() == 0)
+ return;
+
+ // otherwise traverse all the dofs used in
+ // the local matrices and vectors and see
+ // what's there to do
+ const unsigned int n_local_dofs = local_dof_indices.size();
+ for (unsigned int i=0; i<n_local_dofs; ++i)
+ {
+ const std::map<unsigned int, double>::const_iterator
+ boundary_value = boundary_values.find (local_dof_indices[i]);
+ if (boundary_value != boundary_values.end())
+ {
+ // remove this row, except for the
+ // diagonal element
+ for (unsigned j=0; j<n_local_dofs; ++j)
+ if (i != j)
+ local_matrix(i,j) = 0;
+
+ // replace diagonal entry by its
+ // absolute value to make sure that
+ // everything remains positive, or
+ // by one if zero
+ if (local_matrix(i,i) == 0.)
+ local_matrix(i,i) = 1.;
+ else
+ local_matrix(i,i) = std::fabs(local_matrix(i,i));
+
+ // and replace rhs entry by correct
+ // value
+ local_rhs(i) = local_matrix(i,i) * boundary_value->second;
+
+ // finally do the elimination step
+ // if requested
+ if (eliminate_columns == true)
+ {
+ for (unsigned int row=0; row<n_local_dofs; ++row)
+ if (row != i)
+ {
+ local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
+ local_matrix(row,i) = 0;
+ }
+ }
+ }
+ }
+}
+
+