From 1bc426369f5ffa5abd87a7c5b92c0c14a5741151 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 4 May 2017 18:32:12 +0200 Subject: [PATCH] Generalize CM::distribute_local_to_global for rectangular case. --- .../changes/minor/20170504MartinKronbichler | 5 + include/deal.II/lac/constraint_matrix.h | 24 ++++ .../deal.II/lac/constraint_matrix.templates.h | 123 ++++++++++-------- tests/lac/constraints_local_to_global_rect.cc | 70 ++++++++++ .../constraints_local_to_global_rect.output | 11 ++ 5 files changed, 179 insertions(+), 54 deletions(-) create mode 100644 doc/news/changes/minor/20170504MartinKronbichler create mode 100644 tests/lac/constraints_local_to_global_rect.cc create mode 100644 tests/lac/constraints_local_to_global_rect.output diff --git a/doc/news/changes/minor/20170504MartinKronbichler b/doc/news/changes/minor/20170504MartinKronbichler new file mode 100644 index 0000000000..b0f574af9c --- /dev/null +++ b/doc/news/changes/minor/20170504MartinKronbichler @@ -0,0 +1,5 @@ +New: ConstraintMatrix::distribute_local_to_global can now also assemble to +rectangular matrices where rows and columns are described by different +ConstraintMatrix objects. +
+(Martin Kronbichler, 2017/04/12) diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h index 148640b06a..54504012b4 100644 --- a/include/deal.II/lac/constraint_matrix.h +++ b/include/deal.II/lac/constraint_matrix.h @@ -916,6 +916,30 @@ public: const std::vector &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 *this, whereas + * column indices are constrained according to the given ConstraintMatrix + * column_constraint_matrix. 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 + void + distribute_local_to_global (const FullMatrix &local_matrix, + const std::vector &row_indices, + const ConstraintMatrix &column_constraint_matrix, + const std::vector &column_indices, + MatrixType &global_matrix) const; + /** * This function simultaneously writes elements into matrix and vector, * according to the constraints specified by the calling ConstraintMatrix. diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h index f809b665e2..6ae5b7f256 100644 --- a/include/deal.II/lac/constraint_matrix.templates.h +++ b/include/deal.II/lac/constraint_matrix.templates.h @@ -2365,60 +2365,6 @@ ConstraintMatrix::distribute_local_to_global ( -template -void -ConstraintMatrix::distribute_local_to_global ( - const FullMatrix &local_matrix, - const std::vector &row_indices, - const std::vector &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::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 &cols = scratch_data->columns; - std::vector &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 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 @@ -2541,6 +2487,75 @@ distribute_local_to_global (const FullMatrix & +template +void +ConstraintMatrix::distribute_local_to_global ( + const FullMatrix &local_matrix, + const std::vector &row_indices, + const std::vector &col_indices, + MatrixType &global_matrix) const +{ + distribute_local_to_global(local_matrix, row_indices, *this, + col_indices, global_matrix); +} + + + +template +void +ConstraintMatrix::distribute_local_to_global ( + const FullMatrix &local_matrix, + const std::vector &row_indices, + const ConstraintMatrix &col_constraint_matrix, + const std::vector &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::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 &cols = scratch_data->columns; + std::vector &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 0) + global_matrix.add(row, n_values, &cols[0], &vals[0], false, true); + } +} + + + template void ConstraintMatrix:: diff --git a/tests/lac/constraints_local_to_global_rect.cc b/tests/lac/constraints_local_to_global_rect.cc new file mode 100644 index 0000000000..091d8f1ff7 --- /dev/null +++ b/tests/lac/constraints_local_to_global_rect.cc @@ -0,0 +1,70 @@ +// --------------------------------------------------------------------- +// +// 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 +#include +#include +#include + +#include +#include +#include + +void test () +{ + FullMatrix local(2,2); + local(0,0) = 8.; + local(0,1) = 2.; + local(1,0) = -2.; + local(1,1) = 12.; + FullMatrix 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 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(); +} diff --git a/tests/lac/constraints_local_to_global_rect.output b/tests/lac/constraints_local_to_global_rect.output new file mode 100644 index 0000000000..d212302281 --- /dev/null +++ b/tests/lac/constraints_local_to_global_rect.output @@ -0,0 +1,11 @@ + +DEAL::Same constraints: +0 0 0 0 0 0 +0 0 0 0 7.00e+00 8.00e+00 +0 0 0 0 0 0 +0 0 0 0 -1.00e+00 6.00e+00 +DEAL::Different constraints: +0 0 0 0 0 0 +0 0 2.10e+00 0 0 1.29e+01 +0 0 0 0 0 0 +0 0 -3.00e-01 0 0 5.30e+00 -- 2.39.5