]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Generalize CM::distribute_local_to_global for rectangular case. 4343/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 4 May 2017 16:32:12 +0000 (18:32 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 5 May 2017 06:07:08 +0000 (08:07 +0200)
doc/news/changes/minor/20170504MartinKronbichler [new file with mode: 0644]
include/deal.II/lac/constraint_matrix.h
include/deal.II/lac/constraint_matrix.templates.h
tests/lac/constraints_local_to_global_rect.cc [new file with mode: 0644]
tests/lac/constraints_local_to_global_rect.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170504MartinKronbichler b/doc/news/changes/minor/20170504MartinKronbichler
new file mode 100644 (file)
index 0000000..b0f574a
--- /dev/null
@@ -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.
+<br>
+(Martin Kronbichler, 2017/04/12)
index 148640b06ac54f9986d4a72582b6e4e3ddf11532..54504012b453f05af0d3d701b91889e3fe455e2c 100644 (file)
@@ -916,6 +916,30 @@ public:
                               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.
index f809b665e214419151dcbef44492fecf9b244631..6ae5b7f2563377debdee17c6528204710fb39d60 100644 (file)
@@ -2365,60 +2365,6 @@ ConstraintMatrix::distribute_local_to_global (
 
 
 
-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>
@@ -2541,6 +2487,75 @@ distribute_local_to_global (const FullMatrix<typename MatrixType::value_type>  &
 
 
 
+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::
diff --git a/tests/lac/constraints_local_to_global_rect.cc b/tests/lac/constraints_local_to_global_rect.cc
new file mode 100644 (file)
index 0000000..091d8f1
--- /dev/null
@@ -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 <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();
+}
diff --git a/tests/lac/constraints_local_to_global_rect.output b/tests/lac/constraints_local_to_global_rect.output
new file mode 100644 (file)
index 0000000..d212302
--- /dev/null
@@ -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  

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.