From a2a829daea41795dceb59d0aeb90f451c941b585 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 7 Sep 2022 11:55:57 +0200 Subject: [PATCH] TensorProductMatrixSymmetricSum for matrices with rows/columns filled with zero --- include/deal.II/lac/tensor_product_matrix.h | 46 +++- tests/lac/tensor_product_matrix_07.cc | 210 ++++++++++++++++++ ..._product_matrix_07.with_lapack=true.output | 77 +++++++ 3 files changed, 322 insertions(+), 11 deletions(-) create mode 100644 tests/lac/tensor_product_matrix_07.cc create mode 100644 tests/lac/tensor_product_matrix_07.with_lapack=true.output diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 0d2ef2a199..774cf9d8f7 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -304,7 +304,8 @@ public: * and TensorProductMatrixSymmetricSumBase::eigenvectors, respectively. * Note that the current implementation requires each $M_{d}$ to be symmetric * and positive definite and every $A_{d}$ to be symmetric and invertible but - * not necessarily positive definite. + * not necessarily positive definite. Columns and rows filled with zero are + * ignored. */ void reinit(const std::array, dim> &mass_matrix, @@ -450,13 +451,35 @@ namespace internal { Assert(n_rows == n_cols, ExcNotImplemented()); - auto &&transpose_fill_nm = [](Number * out, - const Number * in, - const unsigned int n, - const unsigned int m) { - for (unsigned int mm = 0; mm < m; ++mm) - for (unsigned int nn = 0; nn < n; ++nn) - out[mm + nn * m] = *(in++); + std::vector constrained_dofs(n_rows, false); + + for (unsigned int i = 0; i < n_rows; ++i) + { + if (mass_matrix[i + i * n_rows] == 0.0) + { + Assert(derivative_matrix[i + i * n_rows] == 0.0, + ExcInternalError()); + + for (unsigned int j = 0; j < n_rows; ++j) + { + Assert(derivative_matrix[i + j * n_rows] == 0, + ExcInternalError()); + Assert(derivative_matrix[j + i * n_rows] == 0, + ExcInternalError()); + } + + constrained_dofs[i] = true; + } + } + + const auto transpose_fill_nm = [&constrained_dofs](Number * out, + const Number * in, + const unsigned int n, + const unsigned int m) { + for (unsigned int mm = 0, c = 0; mm < m; ++mm) + for (unsigned int nn = 0; nn < n; ++nn, ++c) + out[mm + nn * m] = + (mm == nn && constrained_dofs[mm]) ? Number(1.0) : in[c]; }; std::vector> eigenvecs(n_rows); @@ -469,9 +492,10 @@ namespace internal deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy, eigenvecs); AssertDimension(eigenvecs.size(), n_rows); - for (unsigned int i = 0; i < n_rows; ++i) - for (unsigned int j = 0; j < n_cols; ++j, ++eigenvectors) - *eigenvectors = eigenvecs[j][i]; + for (unsigned int i = 0, c = 0; i < n_rows; ++i) + for (unsigned int j = 0; j < n_cols; ++j, ++c) + if (constrained_dofs[i] == false) + eigenvectors[c] = eigenvecs[j][i]; for (unsigned int i = 0; i < n_rows; ++i, ++eigenvalues) *eigenvalues = deriv_copy.eigenvalue(i).real(); diff --git a/tests/lac/tensor_product_matrix_07.cc b/tests/lac/tensor_product_matrix_07.cc new file mode 100644 index 0000000000..b73b64c034 --- /dev/null +++ b/tests/lac/tensor_product_matrix_07.cc @@ -0,0 +1,210 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test TensorProductMatrixSymmetricSum for zero (constrained) rows and colums. +// We consider a single cell with DBC applied to face 2*(dim-1). + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +#include "../testmatrix.h" + +template +std::tuple, FullMatrix> +compute_reference_matrices(const unsigned int fe_degree) +{ + MappingQ1 mapping; + QGauss quadrature(fe_degree + 1); + FE_DGQ fe(fe_degree); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + const unsigned int n_dofs = fe.n_dofs_per_cell(); + + FullMatrix mass_matrix_reference(n_dofs, n_dofs); + FullMatrix derivative_matrix_reference(n_dofs, n_dofs); + + FEValues fe_values(mapping, + fe, + quadrature, + update_values | update_gradients | update_JxW_values); + + fe_values.reinit(dof_handler.begin()); + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + { + mass_matrix_reference(i, j) += + (fe_values.shape_value(i, q_index) * + fe_values.shape_value(j, q_index) * fe_values.JxW(q_index)); + + derivative_matrix_reference(i, j) += + (fe_values.shape_grad(i, q_index) * + fe_values.shape_grad(j, q_index) * fe_values.JxW(q_index)); + } + + + + return std::tuple, FullMatrix>{ + mass_matrix_reference, derivative_matrix_reference}; +} + +template +void +print(const FullMatrix &matrix, const std::string label) +{ + deallog << label << std::endl; + matrix.print_formatted(deallog.get_file_stream(), 10, true, 15); + deallog << std::endl << std::endl; +} + + +void +do_test(const bool zero_out_constraints) +{ + using Number = double; + const unsigned int dim = 2; + const unsigned int n_rows_1d = 4; + + // compute 2D stiffness matrix + const auto reference_matrices_2D = + compute_reference_matrices(n_rows_1d - 1); + auto K_2D = std::get<1>(reference_matrices_2D); + + // ... and apply DBC on face 2*(dim-1) + for (unsigned int j = 0; j < Utilities::pow(n_rows_1d, dim); ++j) + for (unsigned int i = 0; i < Utilities::pow(n_rows_1d, dim - 1); ++i) + { + K_2D[i][j] = 0.0; + K_2D[j][i] = 0.0; + } + + if (zero_out_constraints == false) + for (unsigned int i = 0; i < Utilities::pow(n_rows_1d, dim - 1); ++i) + K_2D[i][i] = 1.0; + + // ... print + print(K_2D, "K_2D:"); + + if (zero_out_constraints) + for (unsigned int i = 0; i < Utilities::pow(n_rows_1d, dim - 1); ++i) + K_2D[i][i] = 1.0; + + K_2D.gauss_jordan(); + + if (zero_out_constraints) + for (unsigned int i = 0; i < Utilities::pow(n_rows_1d, dim - 1); ++i) + K_2D[i][i] = 0.0; + + // ... print + print(K_2D, "K_2D^-1:"); + + // compute 1D stiffness and mass matrix + const auto reference_matrices_1D = + compute_reference_matrices<1, Number>(n_rows_1d - 1); + + // ... setup FDM + std::array, dim> mass_matrix; + std::array, dim> derivative_matrix; + + for (unsigned int d = 0; d < dim; ++d) + { + mass_matrix[d] = std::get<0>(reference_matrices_1D); + derivative_matrix[d] = std::get<1>(reference_matrices_1D); + + // apply constraints + if ((d + 1) == dim) + { + for (unsigned int i = 0; i < n_rows_1d; ++i) + { + mass_matrix[d][i][0] = 0.0; + mass_matrix[d][0][i] = 0.0; + derivative_matrix[d][i][0] = 0.0; + derivative_matrix[d][0][i] = 0.0; + } + + if (zero_out_constraints == false) + { + mass_matrix[d][0][0] = 1.0; + derivative_matrix[d][0][0] = 1.0; + } + } + } + + // ... print matrix + TensorProductMatrixSymmetricSum fdm; + fdm.reinit(mass_matrix, derivative_matrix); + + FullMatrix matrix(fdm.m(), fdm.m()); + + for (unsigned int i = 0; i < fdm.m(); ++i) + { + Vector dst(fdm.m()); + Vector src(fdm.n()); + + src[i] = 1.0; + + fdm.vmult(make_array_view(dst), make_array_view(src)); + + for (unsigned int j = 0; j < fdm.m(); ++j) + matrix[j][i] = dst[j]; + } + + print(matrix, "K_FDM:"); + + // ... print inverse matrix + for (unsigned int i = 0; i < fdm.m(); ++i) + { + Vector dst(fdm.m()); + Vector src(fdm.n()); + + src[i] = 1.0; + + fdm.apply_inverse(make_array_view(dst), make_array_view(src)); + + for (unsigned int j = 0; j < fdm.m(); ++j) + matrix[j][i] = dst[j]; + } + + print(matrix, "K_FDM^-1:"); +} + + +int +main() +{ + initlog(); + + do_test(true); + + return 0; +} diff --git a/tests/lac/tensor_product_matrix_07.with_lapack=true.output b/tests/lac/tensor_product_matrix_07.with_lapack=true.output new file mode 100644 index 0000000000..ec02d62911 --- /dev/null +++ b/tests/lac/tensor_product_matrix_07.with_lapack=true.output @@ -0,0 +1,77 @@ + +DEAL::K_2D: + + + + + 2.1428571429e+00 -1.5204601112e+00 3.2364873140e-02 3.9682539683e-02 -3.9682539683e-02 -4.0129811229e-01 1.5328223927e-01 -5.9523809524e-02 -6.4513310372e-02 1.4880952381e-01 -3.7893453497e-02 1.2909876605e-02 + -1.5204601112e+00 5.9523809524e+00 -9.9206349206e-01 3.2364873140e-02 -4.0129811229e-01 -9.9206349206e-01 -4.9603174603e-01 1.5328223927e-01 1.4880952381e-01 3.2364873140e-02 1.5328223927e-01 -3.7893453497e-02 + 3.2364873140e-02 -9.9206349206e-01 5.9523809524e+00 -1.5204601112e+00 1.5328223927e-01 -4.9603174603e-01 -9.9206349206e-01 -4.0129811229e-01 -3.7893453497e-02 1.5328223927e-01 3.2364873140e-02 1.4880952381e-01 + 3.9682539683e-02 3.2364873140e-02 -1.5204601112e+00 2.1428571429e+00 -5.9523809524e-02 1.5328223927e-01 -4.0129811229e-01 -3.9682539683e-02 1.2909876605e-02 -3.7893453497e-02 1.4880952381e-01 -6.4513310372e-02 + -3.9682539683e-02 -4.0129811229e-01 1.5328223927e-01 -5.9523809524e-02 2.1428571429e+00 -1.5204601112e+00 3.2364873140e-02 3.9682539683e-02 -2.3310573725e-01 -2.5972559412e-01 1.4880952381e-01 -6.2513051208e-02 + -4.0129811229e-01 -9.9206349206e-01 -4.9603174603e-01 1.5328223927e-01 -1.5204601112e+00 5.9523809524e+00 -9.9206349206e-01 3.2364873140e-02 -2.5972559412e-01 -1.5204601112e+00 -4.0129811229e-01 1.4880952381e-01 + 1.5328223927e-01 -4.9603174603e-01 -9.9206349206e-01 -4.0129811229e-01 3.2364873140e-02 -9.9206349206e-01 5.9523809524e+00 -1.5204601112e+00 1.4880952381e-01 -4.0129811229e-01 -1.5204601112e+00 -2.5972559412e-01 + -5.9523809524e-02 1.5328223927e-01 -4.0129811229e-01 -3.9682539683e-02 3.9682539683e-02 3.2364873140e-02 -1.5204601112e+00 2.1428571429e+00 -6.2513051208e-02 1.4880952381e-01 -2.5972559412e-01 -2.3310573725e-01 + -6.4513310372e-02 1.4880952381e-01 -3.7893453497e-02 1.2909876605e-02 -2.3310573725e-01 -2.5972559412e-01 1.4880952381e-01 -6.2513051208e-02 6.1904761905e-01 -2.3310573725e-01 -6.4513310372e-02 3.9682539683e-02 + 1.4880952381e-01 3.2364873140e-02 1.5328223927e-01 -3.7893453497e-02 -2.5972559412e-01 -1.5204601112e+00 -4.0129811229e-01 1.4880952381e-01 -2.3310573725e-01 2.1428571429e+00 -3.9682539683e-02 -6.4513310372e-02 + -3.7893453497e-02 1.5328223927e-01 3.2364873140e-02 1.4880952381e-01 1.4880952381e-01 -4.0129811229e-01 -1.5204601112e+00 -2.5972559412e-01 -6.4513310372e-02 -3.9682539683e-02 2.1428571429e+00 -2.3310573725e-01 + 1.2909876605e-02 -3.7893453497e-02 1.4880952381e-01 -6.4513310372e-02 -6.2513051208e-02 1.4880952381e-01 -2.5972559412e-01 -2.3310573725e-01 3.9682539683e-02 -6.4513310372e-02 -2.3310573725e-01 6.1904761905e-01 +DEAL:: +DEAL:: +DEAL::K_2D^-1: + + + + + 7.1398913857e-01 2.7936057528e-01 1.2595048784e-01 9.6173972846e-02 3.5466634519e-01 3.3341836231e-01 2.2437530470e-01 2.1636531977e-01 3.7849889433e-01 2.9761957985e-01 2.4739510402e-01 2.1314611331e-01 + 2.7936057528e-01 3.3864576440e-01 1.4763570837e-01 1.2595048784e-01 3.3341836231e-01 3.1114491295e-01 2.4929635364e-01 2.2437530470e-01 2.9761957985e-01 3.0390259267e-01 2.5043815595e-01 2.4739510402e-01 + 1.2595048784e-01 1.4763570837e-01 3.3864576440e-01 2.7936057528e-01 2.2437530470e-01 2.4929635364e-01 3.1114491295e-01 3.3341836231e-01 2.4739510402e-01 2.5043815595e-01 3.0390259267e-01 2.9761957985e-01 + 9.6173972846e-02 1.2595048784e-01 2.7936057528e-01 7.1398913857e-01 2.1636531977e-01 2.2437530470e-01 3.3341836231e-01 3.5466634519e-01 2.1314611331e-01 2.4739510402e-01 2.9761957985e-01 3.7849889433e-01 + 3.5466634519e-01 3.3341836231e-01 2.2437530470e-01 2.1636531977e-01 1.2712646737e+00 7.9275803884e-01 5.0423937608e-01 4.4702982470e-01 1.0548034198e+00 8.7872567784e-01 5.4153032266e-01 5.2719815068e-01 + 3.3341836231e-01 3.1114491295e-01 2.4929635364e-01 2.2437530470e-01 7.9275803884e-01 8.2875711684e-01 5.5249971478e-01 5.0423937608e-01 8.7872567784e-01 8.3420787522e-01 6.1839723928e-01 5.4153032266e-01 + 2.2437530470e-01 2.4929635364e-01 3.1114491295e-01 3.3341836231e-01 5.0423937608e-01 5.5249971478e-01 8.2875711684e-01 7.9275803884e-01 5.4153032266e-01 6.1839723928e-01 8.3420787522e-01 8.7872567784e-01 + 2.1636531977e-01 2.2437530470e-01 3.3341836231e-01 3.5466634519e-01 4.4702982470e-01 5.0423937608e-01 7.9275803884e-01 1.2712646737e+00 5.2719815068e-01 5.4153032266e-01 8.7872567784e-01 1.0548034198e+00 + 3.7849889433e-01 2.9761957985e-01 2.4739510402e-01 2.1314611331e-01 1.0548034198e+00 8.7872567784e-01 5.4153032266e-01 5.2719815068e-01 2.7334715515e+00 1.0955555964e+00 6.6133919416e-01 4.8205449545e-01 + 2.9761957985e-01 3.0390259267e-01 2.5043815595e-01 2.4739510402e-01 8.7872567784e-01 8.3420787522e-01 6.1839723928e-01 5.4153032266e-01 1.0955555964e+00 1.3563953245e+00 6.9222571734e-01 6.6133919416e-01 + 2.4739510402e-01 2.5043815595e-01 3.0390259267e-01 2.9761957985e-01 5.4153032266e-01 6.1839723928e-01 8.3420787522e-01 8.7872567784e-01 6.6133919416e-01 6.9222571734e-01 1.3563953245e+00 1.0955555964e+00 + 2.1314611331e-01 2.4739510402e-01 2.9761957985e-01 3.7849889433e-01 5.2719815068e-01 5.4153032266e-01 8.7872567784e-01 1.0548034198e+00 4.8205449545e-01 6.6133919416e-01 1.0955555964e+00 2.7334715515e+00 +DEAL:: +DEAL:: +DEAL::K_FDM: + + + + + 2.1428571429e+00 -1.5204601112e+00 3.2364873140e-02 3.9682539683e-02 -3.9682539683e-02 -4.0129811229e-01 1.5328223927e-01 -5.9523809524e-02 -6.4513310372e-02 1.4880952381e-01 -3.7893453497e-02 1.2909876605e-02 + -1.5204601112e+00 5.9523809524e+00 -9.9206349206e-01 3.2364873140e-02 -4.0129811229e-01 -9.9206349206e-01 -4.9603174603e-01 1.5328223927e-01 1.4880952381e-01 3.2364873140e-02 1.5328223927e-01 -3.7893453497e-02 + 3.2364873140e-02 -9.9206349206e-01 5.9523809524e+00 -1.5204601112e+00 1.5328223927e-01 -4.9603174603e-01 -9.9206349206e-01 -4.0129811229e-01 -3.7893453497e-02 1.5328223927e-01 3.2364873140e-02 1.4880952381e-01 + 3.9682539683e-02 3.2364873140e-02 -1.5204601112e+00 2.1428571429e+00 -5.9523809524e-02 1.5328223927e-01 -4.0129811229e-01 -3.9682539683e-02 1.2909876605e-02 -3.7893453497e-02 1.4880952381e-01 -6.4513310372e-02 + -3.9682539683e-02 -4.0129811229e-01 1.5328223927e-01 -5.9523809524e-02 2.1428571429e+00 -1.5204601112e+00 3.2364873140e-02 3.9682539683e-02 -2.3310573725e-01 -2.5972559412e-01 1.4880952381e-01 -6.2513051208e-02 + -4.0129811229e-01 -9.9206349206e-01 -4.9603174603e-01 1.5328223927e-01 -1.5204601112e+00 5.9523809524e+00 -9.9206349206e-01 3.2364873140e-02 -2.5972559412e-01 -1.5204601112e+00 -4.0129811229e-01 1.4880952381e-01 + 1.5328223927e-01 -4.9603174603e-01 -9.9206349206e-01 -4.0129811229e-01 3.2364873140e-02 -9.9206349206e-01 5.9523809524e+00 -1.5204601112e+00 1.4880952381e-01 -4.0129811229e-01 -1.5204601112e+00 -2.5972559412e-01 + -5.9523809524e-02 1.5328223927e-01 -4.0129811229e-01 -3.9682539683e-02 3.9682539683e-02 3.2364873140e-02 -1.5204601112e+00 2.1428571429e+00 -6.2513051208e-02 1.4880952381e-01 -2.5972559412e-01 -2.3310573725e-01 + -6.4513310372e-02 1.4880952381e-01 -3.7893453497e-02 1.2909876605e-02 -2.3310573725e-01 -2.5972559412e-01 1.4880952381e-01 -6.2513051208e-02 6.1904761905e-01 -2.3310573725e-01 -6.4513310372e-02 3.9682539683e-02 + 1.4880952381e-01 3.2364873140e-02 1.5328223927e-01 -3.7893453497e-02 -2.5972559412e-01 -1.5204601112e+00 -4.0129811229e-01 1.4880952381e-01 -2.3310573725e-01 2.1428571429e+00 -3.9682539683e-02 -6.4513310372e-02 + -3.7893453497e-02 1.5328223927e-01 3.2364873140e-02 1.4880952381e-01 1.4880952381e-01 -4.0129811229e-01 -1.5204601112e+00 -2.5972559412e-01 -6.4513310372e-02 -3.9682539683e-02 2.1428571429e+00 -2.3310573725e-01 + 1.2909876605e-02 -3.7893453497e-02 1.4880952381e-01 -6.4513310372e-02 -6.2513051208e-02 1.4880952381e-01 -2.5972559412e-01 -2.3310573725e-01 3.9682539683e-02 -6.4513310372e-02 -2.3310573725e-01 6.1904761905e-01 +DEAL:: +DEAL:: +DEAL::K_FDM^-1: + + + + + 7.1398913857e-01 2.7936057528e-01 1.2595048784e-01 9.6173972846e-02 3.5466634519e-01 3.3341836231e-01 2.2437530470e-01 2.1636531977e-01 3.7849889433e-01 2.9761957985e-01 2.4739510402e-01 2.1314611331e-01 + 2.7936057528e-01 3.3864576440e-01 1.4763570837e-01 1.2595048784e-01 3.3341836231e-01 3.1114491295e-01 2.4929635364e-01 2.2437530470e-01 2.9761957985e-01 3.0390259267e-01 2.5043815595e-01 2.4739510402e-01 + 1.2595048784e-01 1.4763570837e-01 3.3864576440e-01 2.7936057528e-01 2.2437530470e-01 2.4929635364e-01 3.1114491295e-01 3.3341836231e-01 2.4739510402e-01 2.5043815595e-01 3.0390259267e-01 2.9761957985e-01 + 9.6173972846e-02 1.2595048784e-01 2.7936057528e-01 7.1398913857e-01 2.1636531977e-01 2.2437530470e-01 3.3341836231e-01 3.5466634519e-01 2.1314611331e-01 2.4739510402e-01 2.9761957985e-01 3.7849889433e-01 + 3.5466634519e-01 3.3341836231e-01 2.2437530470e-01 2.1636531977e-01 1.2712646737e+00 7.9275803884e-01 5.0423937608e-01 4.4702982470e-01 1.0548034198e+00 8.7872567784e-01 5.4153032266e-01 5.2719815068e-01 + 3.3341836231e-01 3.1114491295e-01 2.4929635364e-01 2.2437530470e-01 7.9275803884e-01 8.2875711684e-01 5.5249971478e-01 5.0423937608e-01 8.7872567784e-01 8.3420787522e-01 6.1839723928e-01 5.4153032266e-01 + 2.2437530470e-01 2.4929635364e-01 3.1114491295e-01 3.3341836231e-01 5.0423937608e-01 5.5249971478e-01 8.2875711684e-01 7.9275803884e-01 5.4153032266e-01 6.1839723928e-01 8.3420787522e-01 8.7872567784e-01 + 2.1636531977e-01 2.2437530470e-01 3.3341836231e-01 3.5466634519e-01 4.4702982470e-01 5.0423937608e-01 7.9275803884e-01 1.2712646737e+00 5.2719815068e-01 5.4153032266e-01 8.7872567784e-01 1.0548034198e+00 + 3.7849889433e-01 2.9761957985e-01 2.4739510402e-01 2.1314611331e-01 1.0548034198e+00 8.7872567784e-01 5.4153032266e-01 5.2719815068e-01 2.7334715515e+00 1.0955555964e+00 6.6133919416e-01 4.8205449545e-01 + 2.9761957985e-01 3.0390259267e-01 2.5043815595e-01 2.4739510402e-01 8.7872567784e-01 8.3420787522e-01 6.1839723928e-01 5.4153032266e-01 1.0955555964e+00 1.3563953245e+00 6.9222571734e-01 6.6133919416e-01 + 2.4739510402e-01 2.5043815595e-01 3.0390259267e-01 2.9761957985e-01 5.4153032266e-01 6.1839723928e-01 8.3420787522e-01 8.7872567784e-01 6.6133919416e-01 6.9222571734e-01 1.3563953245e+00 1.0955555964e+00 + 2.1314611331e-01 2.4739510402e-01 2.9761957985e-01 3.7849889433e-01 5.2719815068e-01 5.4153032266e-01 8.7872567784e-01 1.0548034198e+00 4.8205449545e-01 6.6133919416e-01 1.0955555964e+00 2.7334715515e+00 +DEAL:: +DEAL:: -- 2.39.5