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.
-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
-{
- typedef typename MatrixType::value_type number;
-
- AssertDimension (local_matrix.m(), row_indices.size());
- AssertDimension (local_matrix.n(), col_indices.size());
- //Assert (sorted == true, ExcMatrixNotClosed());
-
- const size_type n_local_row_dofs = row_indices.size();
- const size_type n_local_col_dofs = col_indices.size();
-
- typename internals::ConstraintMatrixData<number>::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);
- 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);
- }
-}
-
-
// similar function as above, but now specialized for block matrices. See the
// other function for additional comments.
template <typename MatrixType, typename VectorType>
+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<number>::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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of
+// ConstraintMatrix::distribute_local_to_global for row and column indices
+// with different constraints on the rows and columns.
+
+#include "../tests.h"
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/constraint_matrix.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+void test ()
+{
+ FullMatrix<double> local(2,2);
+ local(0,0) = 8.;
+ local(0,1) = 2.;
+ local(1,0) = -2.;
+ local(1,1) = 12.;
+ FullMatrix<double> global1(4,6), global2(4,6);
+ ConstraintMatrix cm1, cm2;
+ cm1.add_line(2);
+ cm1.add_entry(2, 1, 0.5);
+ cm1.add_entry(2, 3, 0.5);
+ cm1.close();
+ cm2.add_line(4);
+ cm2.add_entry(4, 2, 0.3);
+ cm2.add_entry(4, 5, 0.7);
+ cm2.close();
+
+ std::vector<unsigned int> indices1(2), indices2(2);
+ indices1[0] = 1;
+ indices1[1] = 2;
+ indices2[0] = 4;
+ indices2[1] = 5;
+ cm1.distribute_local_to_global(local, indices1, indices2, global1);
+ deallog << "Same constraints: " << std::endl;
+ global1.print_formatted(deallog.get_file_stream(), 2, true, 0, "0");
+ cm1.distribute_local_to_global(local, indices1, cm2, indices2, global2);
+ deallog << "Different constraints: " << std::endl;
+ global2.print_formatted(deallog.get_file_stream(), 2, true, 0, "0");
+}
+
+
+int main ()
+{
+ initlog();
+
+ test();
+}
projects / dealii.git / commitdiff