From d07820c8c61050e1ac459e213c2a8efc5f83e5d2 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 14 Jul 2017 14:10:50 -0500 Subject: [PATCH] Bugfix: Do not assume a fully functional matrix interface for mg matrices --- include/deal.II/multigrid/mg_block_smoother.h | 33 ++++------------- include/deal.II/multigrid/mg_coarse.h | 19 +++++++--- include/deal.II/multigrid/mg_matrix.h | 19 +++------- include/deal.II/multigrid/mg_smoother.h | 36 +++++++++++++++---- 4 files changed, 55 insertions(+), 52 deletions(-) diff --git a/include/deal.II/multigrid/mg_block_smoother.h b/include/deal.II/multigrid/mg_block_smoother.h index 15c780af93..3d5ccd7032 100644 --- a/include/deal.II/multigrid/mg_block_smoother.h +++ b/include/deal.II/multigrid/mg_block_smoother.h @@ -194,33 +194,12 @@ MGSmootherBlock::initialize (const MGMatrixT for (unsigned int i=min; i<=max; ++i) { - // Workaround: The matrix objects supplied to this class do not - // expose information how to reinitialize a vector. Therefore, - // populate vmult and Tvmult by hand... - - matrices[i] = LinearOperator>(); - matrices[i].vmult = [&m, i](BlockVector &v, - const BlockVector &u) - { - m[i].vmult(v, u); - }; - matrices[i].Tvmult = [&m, i](BlockVector &v, - const BlockVector &u) - { - m[i].Tvmult(v, u); - }; - - smoothers[i] = LinearOperator>(); - smoothers[i].vmult = [&s, i](BlockVector &v, - const BlockVector &u) - { - s[i].vmult(v, u); - }; - smoothers[i].Tvmult = [&s, i](BlockVector &v, - const BlockVector &u) - { - s[i].Tvmult(v, u); - }; + // Workaround: Unfortunately, not every "m[i]" object has a + // rich enough interface to populate reinit_(domain|range)_vector. + // Thus, apply an empty LinearOperator exemplar. + matrices[i] = linear_operator >( + LinearOperator >(), m[i]); + smoothers[i] = linear_operator >(matrices[i], s[i]); } } diff --git a/include/deal.II/multigrid/mg_coarse.h b/include/deal.II/multigrid/mg_coarse.h index 04b28f05f2..f086cecc74 100644 --- a/include/deal.II/multigrid/mg_coarse.h +++ b/include/deal.II/multigrid/mg_coarse.h @@ -376,8 +376,11 @@ MGCoarseGridLACIteration : solver(&s, typeid(*this).name()) { - matrix = linear_operator(m); - precondition = linear_operator(m, p); + // Workaround: Unfortunately, not every "m" object has a rich enough + // interface to populate reinit_(domain|range)_vector. Thus, supply an + // empty LinearOperator exemplar. + matrix = linear_operator(LinearOperator(), m); + precondition = linear_operator(matrix, p); } @@ -398,8 +401,11 @@ MGCoarseGridLACIteration const PreconditionerType &p) { solver = &s; - matrix = linear_operator(m); - precondition = linear_operator(m, p); + // Workaround: Unfortunately, not every "m" object has a rich enough + // interface to populate reinit_(domain|range)_vector. Thus, supply an + // empty LinearOperator exemplar. + matrix = linear_operator(LinearOperator(), m); + precondition = linear_operator(matrix, p); } @@ -433,7 +439,10 @@ void MGCoarseGridLACIteration ::set_matrix(const MatrixType &m) { - matrix = linear_operator(m); + // Workaround: Unfortunately, not every "m" object has a rich enough + // interface to populate reinit_(domain|range)_vector. Thus, supply an + // empty LinearOperator exemplar. + matrix = linear_operator(LinearOperator(), m); } DEAL_II_ENABLE_EXTRA_DIAGNOSTICS diff --git a/include/deal.II/multigrid/mg_matrix.h b/include/deal.II/multigrid/mg_matrix.h index 8333cf21e8..3f354be286 100644 --- a/include/deal.II/multigrid/mg_matrix.h +++ b/include/deal.II/multigrid/mg_matrix.h @@ -185,20 +185,11 @@ namespace mg matrices.resize(p.min_level(), p.max_level()); for (unsigned int level = p.min_level(); level <= p.max_level(); ++level) { - // Workaround: The matrix objects supplied to this class do not - // expose information how to reinitialize a vector. Therefore, - // populate vmult and Tvmult by hand... - matrices[level] = LinearOperator(); - matrices[level].vmult = [&p, level]( - VectorType &v, const VectorType &u) - { - p[level].vmult(v, u); - }; - matrices[level].Tvmult = [&p, level]( - VectorType &v, const VectorType &u) - { - p[level].Tvmult(v, u); - }; + // Workaround: Unfortunately, not every "p[level]" object has a + // rich enough interface to populate reinit_(domain|range)_vector. + // Thus, apply an empty LinearOperator exemplar. + matrices[level] = + linear_operator(LinearOperator(), p[level]); } } diff --git a/include/deal.II/multigrid/mg_smoother.h b/include/deal.II/multigrid/mg_smoother.h index 29726dc887..3c577b66ad 100644 --- a/include/deal.II/multigrid/mg_smoother.h +++ b/include/deal.II/multigrid/mg_smoother.h @@ -805,7 +805,11 @@ MGSmootherRelaxation::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = linear_operator(m[i]); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = + linear_operator(LinearOperator(), m[i]); smoothers[i].initialize(m[i], data); } } @@ -830,7 +834,11 @@ MGSmootherRelaxation::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = linear_operator(m[i]); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = + linear_operator(LinearOperator(), m[i]); smoothers[i].initialize(m[i], data[i]); } } @@ -852,7 +860,11 @@ MGSmootherRelaxation::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = &(m[i].block(row, col)); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = linear_operator(LinearOperator(), + m[i].block(row, col)); smoothers[i].initialize(m[i].block(row, col), data); } } @@ -879,7 +891,11 @@ MGSmootherRelaxation::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = &(m[i].block(row, col)); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = linear_operator(LinearOperator(), + m[i].block(row, col)); smoothers[i].initialize(m[i].block(row, col), data[i]); } } @@ -1009,7 +1025,11 @@ MGSmootherPrecondition::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = linear_operator(m[i]); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = + linear_operator(LinearOperator(), m[i]); smoothers[i].initialize(m[i], data); } } @@ -1036,7 +1056,11 @@ MGSmootherPrecondition::initialize for (unsigned int i=min; i<=max; ++i) { - matrices[i] = linear_operator(m[i]); + // Workaround: Unfortunately, not every "m[i]" object has a rich + // enough interface to populate reinit_(domain|range)_vector. Thus, + // apply an empty LinearOperator exemplar. + matrices[i] = + linear_operator(LinearOperator(), m[i]); smoothers[i].initialize(m[i], data[i]); } } -- 2.39.5