From 08fa5fa6b0fd623584818b259514b08f7aff839a Mon Sep 17 00:00:00 2001 From: mbreit Date: Thu, 27 Jun 2019 11:10:20 +0200 Subject: [PATCH 1/2] Corrected hanging node constraints In cases where the hanging dof i had a connection to dof j in the sparse matrix, but not j to i (asymmetric matrix) the constraint did not work properly. --- .../p1_continuity_constraints_impl.h | 606 +++++++++--------- 1 file changed, 318 insertions(+), 288 deletions(-) diff --git a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h index c743a60e0..0a57e481f 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h +++ b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h @@ -38,44 +38,32 @@ namespace ug { -/// sets a matrix row corresponding to averaging the constrained index template -void SetInterpolation(TMatrix& A, - std::vector & constrainedIndex, - std::vector >& vConstrainingIndex, - bool assembleLinearProblem = true) +static void setConstrainedRow +( + TMatrix& J, + size_t constrdInd, + const std::vector& vConstrgInd, + bool assembleLinearProblem +) { - //typedef typename TMatrix::row_iterator row_iterator; - //typedef typename TMatrix::value_type block_type; + const size_t nConstrg = vConstrgInd.size(); - // check number of indices passed - for(size_t i = 0; i < vConstrainingIndex.size(); ++i) - UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(), - "Wrong number of indices."); + // set a unity row + SetRow(J, constrdInd, 0.0); + J(constrdInd, constrdInd) = 1.0; -// loop all constrained dofs - for(size_t i = 0; i < constrainedIndex.size(); ++i) + // set an interpolation row + // this is only required if assembling for a linear problem. + if (assembleLinearProblem) { - const size_t ai = constrainedIndex[i]; - - // remove all couplings - SetRow(A, ai, 0.0); - - // set diag of row to identity - A(ai, ai) = 1.0; - - // set coupling to all constraining dofs the inverse of the - // number of constraining dofs - // This is only required if assembling for a linear problem. - if (assembleLinearProblem) - { - const number frac = -1.0/(vConstrainingIndex.size()); - for(size_t j=0; j < vConstrainingIndex.size(); ++j) - A(ai, vConstrainingIndex[j][i]) = frac; - } + const number frac = 1.0 / nConstrg; + for (size_t k = 0; k < nConstrg; ++k) + J(constrdInd, vConstrgInd[k]) = -frac; } } + template void InterpolateValues(TVector& u, std::vector& constrainedIndex, @@ -109,163 +97,6 @@ void InterpolateValues(TVector& u, -template -void SplitAddRow_Symmetric(TMatrix& A, - std::vector& constrainedIndex, - std::vector >& vConstrainingIndex) -{ - typedef typename TMatrix::value_type block_type; - typedef typename TMatrix::row_iterator row_iterator; - - size_t nConstrg = vConstrainingIndex.size(); - UG_ASSERT(nConstrg, "There have to be constraining indices!"); - - // check number of indices passed - for(size_t i = 0; i < nConstrg; ++i) - UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(), - "Wrong number of indices."); - -// scaling factor - const number frac = 1.0 / nConstrg; - -// handle each constrained index - for(size_t i = 0; i < constrainedIndex.size(); ++i) - { - const size_t algInd = constrainedIndex[i]; - - // add coupling constrained index -> constrained index - // we don't have to adjust the block itself, since the row of - // constraints will be set to interpolation afterwards - block_type diagEntry = A(algInd, algInd); - - // scale by weight - diagEntry *= frac*frac; - - // add coupling constrained dof -> constrained dof - for(size_t k = 0; k < vConstrainingIndex.size(); ++k) - for(size_t m = 0; m < vConstrainingIndex.size(); ++m) - A(vConstrainingIndex[k][i], vConstrainingIndex[m][i]) += diagEntry; - -#if 0 - - // (handled separately as it would otherwise trigger an assert - - // manipulation of matrix row while iterating over it) - const block_type block = A(algInd, algInd); - size_t nBlockCols = GetCols(block); - UG_ASSERT(nBlockCols == constrainedIndex.size(), - "Number of block cols and number of constrained DoFs in hanging vertex not equal."); - for (size_t blockIndConn = 0; blockIndConn < nBlockCols; ++blockIndConn) - { - const number val = BlockRef(block, blockInd, blockIndConn) * frac*frac; - for (size_t k = 0; k < nConstrg; ++k) - for (size_t m = 0; m < nConstrg; ++m) - DoFRef(A, vConstrainingIndex[k][i], vConstrainingIndex[m][blockIndConn]) += val; - } -#endif - - // loop coupling between constrained dof -> any other dof - for(row_iterator conn = A.begin_row(algInd); conn != A.end_row(algInd); ++conn) - { - const size_t algIndConn = conn.index(); - - // diagonal entry already handled - if (algIndConn == algInd) - continue; - - // warning: do NOT use references here! - // they might become invalid when A is accessed at an entry that does not exist yet - // FIXME: This will only work properly if there is an entry A(i,j) for any entry - // A(j,i) in the matrix! Is this always the case!? - block_type block = conn.value(); - block_type blockT = A(algIndConn, algInd); - - // multiply the cpl value by the inverse number of constraining - block *= frac; - blockT *= frac; - - // add the coupling to the constraining indices rows - for (size_t k = 0; k < nConstrg; ++k) - { - UG_ASSERT(vConstrainingIndex[k][i] != constrainedIndex[i], - "Modifying 'this' (=conn referenced) matrix row is invalid!" << constrainedIndex[i]); - A(vConstrainingIndex[k][i], algIndConn) += block; - A(algIndConn, vConstrainingIndex[k][i]) += blockT; - } - - // set the split coupling to zero - // this needs to be done only in columns, since the row associated to - // the constrained index will be set to unity (or interpolation). - A(algIndConn, algInd) = 0.0; - } - } -} - -template -void SplitAddRow_OneSide(TMatrix& A, - std::vector& constrainedIndex, - std::vector >& vConstrainingIndex) -{ - typedef typename TMatrix::value_type block_type; - typedef typename TMatrix::row_iterator row_iterator; - - UG_ASSERT(!vConstrainingIndex.empty(), "There have to be constraining indices!"); - - // check number of indices passed - for(size_t i = 0; i < vConstrainingIndex.size(); ++i) - UG_ASSERT(vConstrainingIndex[i].size() == constrainedIndex.size(), - "Wrong number of indices."); - -// scaling factor - const number frac = 1./(vConstrainingIndex.size()); - - for(size_t i = 0; i < constrainedIndex.size(); ++i) - { - const size_t algInd = constrainedIndex[i]; - - // choose randomly the first dof to add whole row - const size_t addTo = vConstrainingIndex[0][i]; - - block_type diagEntry = A(algInd, algInd); - - // scale by weight - diagEntry *= frac; - - // add coupling - for(size_t k = 0; k < vConstrainingIndex.size(); ++k) - A(addTo, vConstrainingIndex[k][i]) += diagEntry; - - // loop coupling between constrained dof -> any other dof - for(row_iterator conn = A.begin_row(algInd); conn != A.end_row(algInd); ++conn) - { - const size_t algIndConn = conn.index(); - - // skip self-coupling (already handled) - if (algIndConn == algInd) continue; - - // warning: do NOT use references here! - // they might become invalid when A is accessed at an entry that does not exist yet - // FIXME: This will only work properly if there is an entry A(i,j) for any entry - // A(j,i) in the matrix! Is this always the case!? - const block_type block = conn.value(); - block_type blockT = A(algIndConn, algInd); - - blockT *= frac; - - // add the coupling to the constraining indices rows - for(size_t k = 0; k < vConstrainingIndex.size(); ++k) - A(algIndConn, vConstrainingIndex[k][i]) += blockT; - - // coupling due to one side adding - A(addTo, algIndConn) += block; - - // set the split coupling to zero - // this must only be done in columns, since the row associated to - // the constrained index will be set to an interpolation. - A(algIndConn, algInd) = 0.0; - } - } -} - template void SplitAddRhs_Symmetric(TVector& rhs, std::vector & constrainedIndex, @@ -465,45 +296,185 @@ adjust_rhs(vector_type& rhs, const vector_type& u, } } -template -void -SymP1Constraints:: -adjust_jacobian(matrix_type& J, const vector_type& u, - ConstSmartPtr dd, int type, number time, - ConstSmartPtr > vSol, - const number s_a0) +static void fillConstraintMapSymmetric +( + std::map >& constraintMap, + ConstSmartPtr dd +) { - if(this->m_spAssTuner->single_index_assembling_enabled()) - UG_THROW("index-wise assemble routine is not " - "implemented for SymP1Constraints \n"); - -// storage for indices and vertices std::vector > vConstrainingInd; std::vector constrainedInd; + std::vector constrainers; std::vector vConstrainingVrt; -// get begin end of hanging vertices + // get begin end of hanging vertices DoFDistribution::traits::const_iterator iter, iterEnd; iter = dd->begin(); iterEnd = dd->end(); -// loop constrained vertices - for(; iter != iterEnd; ++iter) + // loop constrained vertices + for (; iter != iterEnd; ++iter) { - // get hanging vert + // get hanging vert ConstrainedVertex* hgVrt = *iter; - // get algebra indices for constrained and constraining vertices + // get algebra indices for constrained and constraining vertices get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd); - // Split using indices - SplitAddRow_Symmetric(J, constrainedInd, vConstrainingInd); + // save in constraint map + const size_t nInd = constrainedInd.size(); + const size_t nConstrainers = vConstrainingInd.size(); + constrainers.resize(nConstrainers); + for (size_t i = 0; i < nInd; ++i) + { + for (size_t j = 0; j < nConstrainers; ++j) + constrainers[j] = vConstrainingInd[j][i]; + constraintMap[constrainedInd[i]] = constrainers; + } + } +} + +template +static void splitConstrainedRowSymmetric +( + TMatrix& J, + size_t constrdInd, + const std::vector& vConstrgInd +) +{ + typedef typename TMatrix::value_type block_type; + typedef typename TMatrix::row_iterator row_iterator; + + const size_t nConstrg = vConstrgInd.size(); + const number frac = 1.0 / nConstrg; + + // distribute row equally to constrainers + for (row_iterator conn = J.begin_row(constrdInd); conn != J.end_row(constrdInd); ++conn) + { + const size_t connInd = conn.index(); + block_type block = conn.value(); + block *= frac; + for (size_t k = 0; k < nConstrg; ++k) + J(vConstrgInd[k], connInd) += block; + } +} + +template +static void splitConstrainedRhsEntrySymmetric +( + TVector& rhs, + size_t constrdInd, + const std::vector& vConstrgInd +) +{ + typedef typename TVector::value_type block_type; + + const size_t nConstrg = vConstrgInd.size(); + const number frac = 1.0 / nConstrg; + + // get constrained rhs + // modify block directly since set to zero afterwards + block_type& val = rhs[constrdInd]; + val *= frac; + + // split equally on all constraining indices + for (size_t i = 0; i < nConstrg; ++i) + rhs[vConstrgInd[i]] += val; + + // set rhs to zero for constrained index + val = 0.0; +} + +template +static void splitConstrainedCols +( + TMatrix& J, + const std::map >& constraintMap +) +{ + typedef typename TMatrix::value_type block_type; + typedef typename TMatrix::row_iterator row_iterator; + typedef typename std::map >::const_iterator constraint_iter; + + const size_t nRows = J.num_rows(); + std::vector tmpConstraints; + for (size_t r = 0; r < nRows; ++r) + { + // leave out constrained rows (already final) + if (constraintMap.find(r) != constraintMap.end()) + continue; + + // loop row + tmpConstraints.clear(); + for (row_iterator conn = J.begin_row(r); conn != J.end_row(r); ++conn) + { + const size_t connInd = conn.index(); + + // remember constrained columns + constraint_iter it = constraintMap.find(connInd); + if (it != constraintMap.end()) + tmpConstraints.push_back(it); + } + + // split constrained columns to constrainers (separate iteration to not throw off iterators) + const size_t nConstraints = tmpConstraints.size(); + for (size_t i = 0; i < nConstraints; ++i) + { + const size_t cdi = tmpConstraints[i]->first; + const std::vector& vConstrgInd = tmpConstraints[i]->second; + const size_t nConstrg = vConstrgInd.size(); + const number frac = 1.0 / nConstrg; + + // add to constraining column entries + block_type block = J(r, cdi); + block *= frac; + for (size_t k = 0; k < nConstrg; ++k) + J(r, vConstrgInd[k]) += block; - // set interpolation - SetInterpolation(J, constrainedInd, vConstrainingInd, m_bAssembleLinearProblem); + // set constrained column entry zero + J(r, cdi) = 0.0; + } } } + + +template +void +SymP1Constraints:: +adjust_jacobian(matrix_type& J, const vector_type& u, + ConstSmartPtr dd, int type, number time, + ConstSmartPtr > vSol, + const number s_a0) +{ + if (this->m_spAssTuner->single_index_assembling_enabled()) + UG_THROW("index-wise assemble routine is not " + "implemented for SymP1Constraints \n"); + + std::map > constraintMap; + typedef typename std::map >::const_iterator constraint_iter; + fillConstraintMapSymmetric(constraintMap, dd); + + // split constrained rows to constrainers and set interpolation row (or identity) + constraint_iter it = constraintMap.begin(); + constraint_iter itEnd = constraintMap.end(); + for (; it != itEnd; ++it) + { + size_t constrdInd = it->first; + const std::vector& vConstrgInd = it->second; + + // split original row + splitConstrainedRowSymmetric(J, constrdInd, vConstrgInd); + + // set identity row (or interpolating row in linear case) + setConstrainedRow(J, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + } + + // apply chain rule, i.e., also split constrained columns + // (unfortunately, we have to iterate the whole matrix for this) + splitConstrainedCols(J, constraintMap); +} + template void SymP1Constraints:: @@ -512,38 +483,35 @@ adjust_linear(matrix_type& mat, vector_type& rhs, { m_bAssembleLinearProblem = true; - if(this->m_spAssTuner->single_index_assembling_enabled()) + if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " "implemented for SymP1Constraints \n"); -// storage for indices and vertices - std::vector > vConstrainingInd; - std::vector constrainedInd; - std::vector vConstrainingVrt; - -// get begin end of hanging vertices - DoFDistribution::traits::const_iterator iter, iterEnd; - iter = dd->begin(); - iterEnd = dd->end(); + std::map > constraintMap; + typedef typename std::map >::const_iterator constraint_iter; + fillConstraintMapSymmetric(constraintMap, dd); -// loop constrained vertices - for(; iter != iterEnd; ++iter) + // split constrained rows to constrainers and set interpolation row (or identity) + constraint_iter it = constraintMap.begin(); + constraint_iter itEnd = constraintMap.end(); + for (; it != itEnd; ++it) { - // get hanging vert - ConstrainedVertex* hgVrt = *iter; + size_t constrdInd = it->first; + const std::vector& vConstrgInd = it->second; - // get algebra indices for constrained and constraining vertices - get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd); + // split original row + splitConstrainedRowSymmetric(mat, constrdInd, vConstrgInd); - // Split using indices - SplitAddRow_Symmetric(mat, constrainedInd, vConstrainingInd); + // split original rhs entry + splitConstrainedRhsEntrySymmetric(rhs, constrdInd, vConstrgInd); - // set interpolation - SetInterpolation(mat, constrainedInd, vConstrainingInd, true); - - // adapt rhs - SplitAddRhs_Symmetric(rhs, constrainedInd, vConstrainingInd); + // set identity row (or interpolating row in linear case) + setConstrainedRow(mat, constrdInd, vConstrgInd, m_bAssembleLinearProblem); } + + // apply chain rule, i.e., also split constrained columns + // (unfortunately, we have to iterate the whole matrix for this) + splitConstrainedCols(mat, constraintMap); } template @@ -694,7 +662,7 @@ adjust_correction template struct SortVertexPos { - SortVertexPos(SmartPtr > spDomain) + SortVertexPos(ConstSmartPtr > spDomain) : m_aaPos(spDomain->position_accessor()) {} @@ -702,7 +670,7 @@ struct SortVertexPos { {UG_THROW(dim <<" not implemented.");} protected: - typename Domain::position_accessor_type& m_aaPos; + const typename Domain::position_accessor_type& m_aaPos; }; template<> @@ -911,53 +879,126 @@ adjust_rhs(vector_type& rhs, const vector_type& u, } } -template -void -OneSideP1Constraints:: -adjust_jacobian(matrix_type& J, const vector_type& u, - ConstSmartPtr dd, int type, number time, - ConstSmartPtr > vSol, - const number s_a0) -{ - if(this->m_spAssTuner->single_index_assembling_enabled()) - UG_THROW("index-wise assemble routine is not " - "implemented for OneSideP1Constraints \n"); -// storage for indices and vertices +template +static void fillConstraintMapOneSide +( + std::map >& constraintMap, + ConstSmartPtr dd, + ConstSmartPtr dom +) +{ std::vector > vConstrainingInd; - std::vector constrainedInd; + std::vector constrainedInd; + std::vector constrainers; std::vector vConstrainingVrt; #ifdef UG_PARALLEL - SortVertexPos sortVertexPos(this->approximation_space()->domain()); + SortVertexPos sortVertexPos(dom); #endif -// get begin end of hanging vertices + // get begin end of hanging vertices DoFDistribution::traits::const_iterator iter, iterEnd; iter = dd->begin(); iterEnd = dd->end(); -// loop constrained vertices - for(; iter != iterEnd; ++iter) + // loop constrained vertices + for (; iter != iterEnd; ++iter) { - // get hanging vert + // get hanging vert ConstrainedVertex* hgVrt = *iter; - // get algebra indices for constrained and constraining vertices + // get algebra indices for constrained and constraining vertices #ifdef UG_PARALLEL get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos); #else get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd); #endif - // Split using indices - SplitAddRow_OneSide(J, constrainedInd, vConstrainingInd); + // save in constraint map + const size_t nInd = constrainedInd.size(); + const size_t nConstrainers = vConstrainingInd.size(); + constrainers.resize(nConstrainers); + for (size_t i = 0; i < nInd; ++i) + { + for (size_t j = 0; j < nConstrainers; ++j) + constrainers[j] = vConstrainingInd[j][i]; + constraintMap[constrainedInd[i]] = constrainers; + } + } +} - // set interpolation - SetInterpolation(J, constrainedInd, vConstrainingInd, m_bAssembleLinearProblem); +template +static void splitConstrainedRowOneSide +( + TMatrix& J, + size_t constrdInd, + const std::vector& vConstrgInd +) +{ + typedef typename TMatrix::row_iterator row_iterator; + + // add complete row to first constrainer + for (row_iterator conn = J.begin_row(constrdInd); conn != J.end_row(constrdInd); ++conn) + { + const size_t connInd = conn.index(); + J(vConstrgInd[0], connInd) += conn.value(); } } +template +static void splitConstrainedRhsEntryOneSide +( + TVector& rhs, + size_t constrdInd, + const std::vector& vConstrgInd +) +{ + // add complete entry to first constrainer index + typename TVector::value_type& val = rhs[constrdInd]; + rhs[vConstrgInd[0]] += val; + + // set rhs to zero for constrained index + val = 0.0; +} + + +template +void +OneSideP1Constraints:: +adjust_jacobian(matrix_type& J, const vector_type& u, + ConstSmartPtr dd, int type, number time, + ConstSmartPtr > vSol, + const number s_a0) +{ + if (this->m_spAssTuner->single_index_assembling_enabled()) + UG_THROW("index-wise assemble routine is not " + "implemented for OneSideP1Constraints."); + + std::map > constraintMap; + typedef typename std::map >::const_iterator constraint_iter; + fillConstraintMapOneSide(constraintMap, dd, this->approximation_space()->domain()); + + // split constrained rows to constrainers and set interpolation row (or identity) + constraint_iter it = constraintMap.begin(); + constraint_iter itEnd = constraintMap.end(); + for (; it != itEnd; ++it) + { + size_t constrdInd = it->first; + const std::vector& vConstrgInd = it->second; + + // split original row + splitConstrainedRowOneSide(J, constrdInd, vConstrgInd); + + // set identity row (or interpolating row in linear case) + setConstrainedRow(J, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + } + + // apply chain rule, i.e., also split constrained columns + // (unfortunately, we have to iterate the whole matrix for this) + splitConstrainedCols(J, constraintMap); +} + template void OneSideP1Constraints:: @@ -966,46 +1007,35 @@ adjust_linear(matrix_type& mat, vector_type& rhs, { m_bAssembleLinearProblem = true; - if(this->m_spAssTuner->single_index_assembling_enabled()) + if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " - "implemented for OneSideP1Constraints \n"); - -// storage for indices and vertices - std::vector > vConstrainingInd; - std::vector constrainedInd; - std::vector vConstrainingVrt; - -#ifdef UG_PARALLEL - SortVertexPos sortVertexPos(this->approximation_space()->domain()); -#endif + "implemented for OneSideP1Constraints"); -// get begin end of hanging vertices - DoFDistribution::traits::const_iterator iter, iterEnd; - iter = dd->begin(); - iterEnd = dd->end(); + std::map > constraintMap; + typedef typename std::map >::const_iterator constraint_iter; + fillConstraintMapOneSide(constraintMap, dd, this->approximation_space()->domain()); -// loop constraining edges - for(; iter != iterEnd; ++iter) + // split constrained rows to constrainers and set interpolation row (or identity) + constraint_iter it = constraintMap.begin(); + constraint_iter itEnd = constraintMap.end(); + for (; it != itEnd; ++it) { - // get hanging vert - ConstrainedVertex* hgVrt = *iter; - - // get algebra indices for constrained and constraining vertices -#ifdef UG_PARALLEL - get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd, sortVertexPos); -#else - get_algebra_indices(dd, hgVrt, vConstrainingVrt, constrainedInd, vConstrainingInd); -#endif + size_t constrdInd = it->first; + const std::vector& vConstrgInd = it->second; - // Split using indices - SplitAddRow_OneSide(mat, constrainedInd, vConstrainingInd); + // split original row + splitConstrainedRowOneSide(mat, constrdInd, vConstrgInd); - // Set interpolation - SetInterpolation(mat, constrainedInd, vConstrainingInd, true); + // split original rhs entry + splitConstrainedRhsEntryOneSide(rhs, constrdInd, vConstrgInd); - // adapt rhs - SplitAddRhs_OneSide(rhs, constrainedInd, vConstrainingInd); + // set identity row (or interpolating row in linear case) + setConstrainedRow(mat, constrdInd, vConstrgInd, m_bAssembleLinearProblem); } + + // apply chain rule, i.e., also split constrained columns + // (unfortunately, we have to iterate the whole matrix for this) + splitConstrainedCols(mat, constraintMap); } template From 3f5d44d0dee5bc9f89356465ab1a37016a15109a Mon Sep 17 00:00:00 2001 From: mbreit Date: Sun, 10 Jan 2021 15:10:11 +0100 Subject: [PATCH 2/2] Fix for bad implementation of hanging node constraints for linear problems The implementation no longer distinguishes the linear from the non-linear case. Both handle the constrained nodes by effectively taking them out of the computation by making sure that the defect and corrections of constrained DoFs are always zero and their value in the solution does not change. To make this possible in the linear assembling case, the initial solution had to be introduced as an additional argument to the assembling function. ATTENTION: Unlike before, it is now necessary to call adjust_solution() once after the solver is finished to enforce the hanging node constraints in linear problems. Otherwise, the DoFs of hanging nodes will have value 0 in the solution. --- ugbase/bridge/disc_bridges/algebra_bridge.cpp | 5 ++- .../disc_bridges/domain_disc_bridge.cpp | 1 + .../obstacles/obstacle_constraint_interface.h | 2 +- ugbase/lib_disc/assemble_interface.h | 14 +++++- .../assembled_linear_operator.h | 3 ++ .../assembled_linear_operator_impl.h | 35 ++++++++++++--- .../nested_iteration/nested_iteration_impl.h | 2 +- .../constraints/constraint_interface.h | 2 +- .../p1_continuity_constraints.h | 14 ++---- .../p1_continuity_constraints_impl.h | 44 +++++++------------ .../lagrange_dirichlet_boundary.h | 2 +- .../lagrange_dirichlet_boundary_impl.h | 2 +- ugbase/lib_disc/spatial_disc/domain_disc.h | 15 +++++-- .../lib_disc/spatial_disc/domain_disc_impl.h | 6 +-- .../spatial_disc/domain_disc_interface.h | 4 +- .../lib_disc/time_disc/composite_time_disc.h | 2 +- .../time_disc/composite_time_disc_impl.h | 5 ++- ugbase/lib_disc/time_disc/theta_time_step.h | 4 +- .../lib_disc/time_disc/theta_time_step_impl.h | 4 +- 19 files changed, 98 insertions(+), 68 deletions(-) diff --git a/ugbase/bridge/disc_bridges/algebra_bridge.cpp b/ugbase/bridge/disc_bridges/algebra_bridge.cpp index aa7decf9d..319adc2f5 100644 --- a/ugbase/bridge/disc_bridges/algebra_bridge.cpp +++ b/ugbase/bridge/disc_bridges/algebra_bridge.cpp @@ -123,6 +123,8 @@ static void Algebra(Registry& reg, string parentGroup) .add_method("assemble_defect", static_cast(&T::assemble_defect),"", "d(u)#u#GridLevel", "Assembles Defect on grid level") .add_method("assemble_linear", static_cast(&T::assemble_linear),"", "A#b", "Assembles Matrix and rhs on surface grid.") .add_method("assemble_linear", static_cast(&T::assemble_linear),"", "A#b#GridLevel", "Assembles Matrix and rhs on grid level.") + .add_method("assemble_linear", static_cast(&T::assemble_linear),"", "A#b#u", "Assembles Matrix and rhs on surface grid.") + .add_method("assemble_linear", static_cast(&T::assemble_linear),"", "A#b#u#GridLevel", "Assembles Matrix and rhs on grid level.") .add_method("assemble_rhs", static_cast(&T::assemble_rhs),"", "rhs#u", "assembles right-hand side on surface grid") .add_method("assemble_rhs", static_cast(&T::assemble_rhs),"", "rhs#u#GridLevel", "assembles right-hand side on grid level") .add_method("assemble_rhs", static_cast(&T::assemble_rhs),"", "rhs#GridLevel", "assembles right-hand side on surface grid for linear case") @@ -256,7 +258,8 @@ static void Algebra(Registry& reg, string parentGroup) .add_method("set_discretization", &T::set_discretization) .add_method("set_level", &T::set_level) .add_method("set_dirichlet_values", &T::set_dirichlet_values) - .add_method("init_op_and_rhs", &T::init_op_and_rhs) + .add_method("init_op_and_rhs", static_cast(&T::init_op_and_rhs)) + .add_method("init_op_and_rhs", static_cast(&T::init_op_and_rhs)) .add_method("level", &T::level) .set_construct_as_smart_pointer(true); reg.add_class_to_group(name, "AssembledLinearOperator", tag); diff --git a/ugbase/bridge/disc_bridges/domain_disc_bridge.cpp b/ugbase/bridge/disc_bridges/domain_disc_bridge.cpp index d3f206122..724f54ef5 100644 --- a/ugbase/bridge/disc_bridges/domain_disc_bridge.cpp +++ b/ugbase/bridge/disc_bridges/domain_disc_bridge.cpp @@ -98,6 +98,7 @@ static void DomainAlgebra(Registry& reg, string grp) .add_method("remove", static_cast >)>(&T::remove), "", "Element Discretization") .add_method("add", static_cast >)>(&T::add), "", "DiscItem") .add_method("assemble_linear", static_cast&)>(&T::assemble_linear)) + .add_method("assemble_linear", static_cast&, const GridFunction&)>(&T::assemble_linear)) .add_method("assemble_rhs", static_cast&)>(&T::assemble_rhs)) .add_method("assemble_rhs", static_cast&)>(&T::assemble_rhs)) .add_method("adjust_solution", static_cast&)>(&T::adjust_solution)) diff --git a/ugbase/lib_algebra/operator/preconditioner/projected_gauss_seidel/obstacles/obstacle_constraint_interface.h b/ugbase/lib_algebra/operator/preconditioner/projected_gauss_seidel/obstacles/obstacle_constraint_interface.h index 18ee9de18..3a8d69a88 100644 --- a/ugbase/lib_algebra/operator/preconditioner/projected_gauss_seidel/obstacles/obstacle_constraint_interface.h +++ b/ugbase/lib_algebra/operator/preconditioner/projected_gauss_seidel/obstacles/obstacle_constraint_interface.h @@ -200,7 +200,7 @@ class IObstacleConstraint: }; /// sets unity rows in A and dirichlet values in right-hand side b - void adjust_linear(matrix_type& A, vector_type& b, + void adjust_linear(matrix_type& A, vector_type& b, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0){ UG_LOG("IObstacleConstraint::adjust_linear() \n"); }; diff --git a/ugbase/lib_disc/assemble_interface.h b/ugbase/lib_disc/assemble_interface.h index 1b67ec1a8..f481482a3 100644 --- a/ugbase/lib_disc/assemble_interface.h +++ b/ugbase/lib_disc/assemble_interface.h @@ -151,9 +151,19 @@ class IAssemble * \param[out] b Right-Hand-Side * \param[in] gl Grid Level */ - virtual void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) = 0; + virtual void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl) = 0; + void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) + { + UG_LOGN("HINT: assemble_linear without given solution is deprecated.\n" + " Supply the solution as it will enter the solver as a third argument."); + // Use the rhs as solution. This will reproduce previous behavior + // when entries of b are assigned the corresponding entries of u. + assemble_linear(A, b, b, gl); + } + void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u) + {assemble_linear(A, b, u, GridLevel());} void assemble_linear(matrix_type& A, vector_type& b) - {assemble_linear(A,b, GridLevel());} + {assemble_linear(A, b, GridLevel());} /// assembles rhs virtual void assemble_rhs(vector_type& rhs, const vector_type& u, const GridLevel& gl) = 0; diff --git a/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator.h b/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator.h index a381c361c..a3f8535ca 100644 --- a/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator.h +++ b/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator.h @@ -103,6 +103,9 @@ class AssembledLinearOperator : /// initializes the operator and assembles the passed rhs vector void init_op_and_rhs(vector_type& b); + /// initializes the operator and assembles the passed rhs vector + void init_op_and_rhs(vector_type& b, const vector_type& u); + /// compute d = J(u)*c (here, J(u) is a Matrix) virtual void apply(vector_type& d, const vector_type& c); diff --git a/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator_impl.h b/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator_impl.h index 721fb81c6..f0809fd60 100644 --- a/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator_impl.h +++ b/ugbase/lib_disc/operator/linear_operator/assembled_linear_operator_impl.h @@ -88,6 +88,9 @@ AssembledLinearOperator::init_op_and_rhs(vector_type& b) { // todo: check that assembling is linear + UG_LOGN("Hint: init_op_and_rhs without solution vector is deprecated.\n" + "Supply the solution vector (as it will enter the solver) as second argument."); + if(m_spAss.invalid()) UG_THROW("AssembledLinearOperator: Assembling routine not set."); @@ -99,6 +102,24 @@ AssembledLinearOperator::init_op_and_rhs(vector_type& b) " Cannot assemble Matrix and Rhs."); } +// Initialize the operator +template +void +AssembledLinearOperator::init_op_and_rhs(vector_type& b, const vector_type& u) +{ +// todo: check that assembling is linear + + if(m_spAss.invalid()) + UG_THROW("AssembledLinearOperator: Assembling routine not set."); + +// assemble matrix and rhs in one loop + try{ + m_spAss->assemble_linear(*this, b, u, m_gridLevel); + } + UG_CATCH_THROW("AssembledLinearOperator::init_op_and_rhs:" + " Cannot assemble Matrix and Rhs."); +} + template void AssembledLinearOperator::apply(vector_type& d, const vector_type& c) @@ -170,13 +191,6 @@ void AssembleLinearOperatorRhsAndSolution { ASS_PROFILE_BEGIN(ASS_AssembleLinearOperatorRhsAndSolution); -// initialize operator - ASS_PROFILE_BEGIN(ASS_InitOperatorAndRhs); - try{ - op.init_op_and_rhs(b); - }UG_CATCH_THROW("Cannot init the operator (assembling failed)."); - ASS_PROFILE_END(); - // sets the dirichlet values in the solution ASS_PROFILE_BEGIN(ASS_SetDirValues); try{ @@ -184,6 +198,13 @@ void AssembleLinearOperatorRhsAndSolution } UG_CATCH_THROW("Cannot set the dirichlet values in the solution."); ASS_PROFILE_END(); +// initialize operator + ASS_PROFILE_BEGIN(ASS_InitOperatorAndRhs); + try{ + op.init_op_and_rhs(b, u); + }UG_CATCH_THROW("Cannot init the operator (assembling failed)."); + ASS_PROFILE_END(); + ASS_PROFILE_END(); } diff --git a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h index 42ac68923..d823029ed 100644 --- a/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h +++ b/ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h @@ -262,7 +262,7 @@ bool NestedIterationSolver::apply(vector_type& u) // Assemble. NESTED_ITER_PROFILE_BEGIN(NestedIterationAssemble); - m_spAss->assemble_linear(*m_J, *spD, surfGridLevel); // todo: replace for non-linear problems + m_spAss->assemble_linear(*m_J, *spD, u, surfGridLevel); // todo: replace for non-linear problems m_spAss->adjust_solution(u, surfGridLevel); NESTED_ITER_PROFILE_END(); diff --git a/ugbase/lib_disc/spatial_disc/constraints/constraint_interface.h b/ugbase/lib_disc/spatial_disc/constraints/constraint_interface.h index eb4f925c3..03dffaf2a 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/constraint_interface.h +++ b/ugbase/lib_disc/spatial_disc/constraints/constraint_interface.h @@ -86,7 +86,7 @@ class IConstraint const std::vector* vScaleStiff = NULL) = 0; /// adapts matrix and rhs (linear case) to enforce constraints - virtual void adjust_linear(matrix_type& mat, vector_type& rhs, + virtual void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0) = 0; /// adapts a rhs to enforce constraints diff --git a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h index 4702e6dd1..1f9708960 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h +++ b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h @@ -61,7 +61,7 @@ class SymP1Constraints typedef typename algebra_type::vector_type vector_type; public: - SymP1Constraints() : IDomainConstraint(), m_bAssembleLinearProblem(false) {} + SymP1Constraints() : IDomainConstraint() {} virtual ~SymP1Constraints() {} virtual int type() const {return CT_HANGING;} @@ -80,7 +80,7 @@ class SymP1Constraints void adjust_rhs(vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0); - void adjust_linear(matrix_type& mat, vector_type& rhs, + void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0); void adjust_solution(vector_type& u, ConstSmartPtr dd, @@ -104,9 +104,6 @@ class SymP1Constraints int type, number time = 0.0 ); - - protected: - bool m_bAssembleLinearProblem; }; @@ -129,7 +126,7 @@ class OneSideP1Constraints typedef IDomainConstraint base_type; public: - OneSideP1Constraints() : IDomainConstraint(), m_bAssembleLinearProblem(false) {} + OneSideP1Constraints() : IDomainConstraint() {} virtual ~OneSideP1Constraints() {} virtual int type() const {return CT_HANGING;} @@ -148,7 +145,7 @@ class OneSideP1Constraints void adjust_rhs(vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0); - void adjust_linear(matrix_type& mat, vector_type& rhs, + void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0); void adjust_solution(vector_type& u, ConstSmartPtr dd, @@ -172,9 +169,6 @@ class OneSideP1Constraints int type, number time = 0.0 ); - - protected: - bool m_bAssembleLinearProblem; }; }; // namespace ug diff --git a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h index 0a57e481f..a4bcdd9d5 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h +++ b/ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints_impl.h @@ -43,24 +43,12 @@ static void setConstrainedRow ( TMatrix& J, size_t constrdInd, - const std::vector& vConstrgInd, - bool assembleLinearProblem + const std::vector& vConstrgInd ) { - const size_t nConstrg = vConstrgInd.size(); - // set a unity row SetRow(J, constrdInd, 0.0); J(constrdInd, constrdInd) = 1.0; - - // set an interpolation row - // this is only required if assembling for a linear problem. - if (assembleLinearProblem) - { - const number frac = 1.0 / nConstrg; - for (size_t k = 0; k < nConstrg; ++k) - J(constrdInd, vConstrgInd[k]) = -frac; - } } @@ -293,6 +281,10 @@ adjust_rhs(vector_type& rhs, const vector_type& u, // adapt rhs SplitAddRhs_Symmetric(rhs, constrainedInd, vConstrainingInd); + + const size_t nCstrd = constrainedInd.size(); + for (size_t i = 0; i < nCstrd; ++i) + rhs[constrainedInd[i]] = u[constrainedInd[i]]; } } @@ -467,7 +459,7 @@ adjust_jacobian(matrix_type& J, const vector_type& u, splitConstrainedRowSymmetric(J, constrdInd, vConstrgInd); // set identity row (or interpolating row in linear case) - setConstrainedRow(J, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + setConstrainedRow(J, constrdInd, vConstrgInd); } // apply chain rule, i.e., also split constrained columns @@ -478,11 +470,9 @@ adjust_jacobian(matrix_type& J, const vector_type& u, template void SymP1Constraints:: -adjust_linear(matrix_type& mat, vector_type& rhs, +adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time) { - m_bAssembleLinearProblem = true; - if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " "implemented for SymP1Constraints \n"); @@ -504,9 +494,10 @@ adjust_linear(matrix_type& mat, vector_type& rhs, // split original rhs entry splitConstrainedRhsEntrySymmetric(rhs, constrdInd, vConstrgInd); + rhs[constrdInd] = u[constrdInd]; // set identity row (or interpolating row in linear case) - setConstrainedRow(mat, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + setConstrainedRow(mat, constrdInd, vConstrgInd); } // apply chain rule, i.e., also split constrained columns @@ -561,8 +552,6 @@ adjust_prolongation number time ) { - if (m_bAssembleLinearProblem) return; - if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " "implemented for SymP1Constraints \n"); @@ -876,6 +865,10 @@ adjust_rhs(vector_type& rhs, const vector_type& u, // adapt rhs SplitAddRhs_OneSide(rhs, constrainedInd, vConstrainingInd); + + const size_t nCstrd = constrainedInd.size(); + for (size_t i = 0; i < nCstrd; ++i) + rhs[constrainedInd[i]] = u[constrainedInd[i]]; } } @@ -991,7 +984,7 @@ adjust_jacobian(matrix_type& J, const vector_type& u, splitConstrainedRowOneSide(J, constrdInd, vConstrgInd); // set identity row (or interpolating row in linear case) - setConstrainedRow(J, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + setConstrainedRow(J, constrdInd, vConstrgInd); } // apply chain rule, i.e., also split constrained columns @@ -1002,11 +995,9 @@ adjust_jacobian(matrix_type& J, const vector_type& u, template void OneSideP1Constraints:: -adjust_linear(matrix_type& mat, vector_type& rhs, +adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd, int type, number time) { - m_bAssembleLinearProblem = true; - if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " "implemented for OneSideP1Constraints"); @@ -1028,9 +1019,10 @@ adjust_linear(matrix_type& mat, vector_type& rhs, // split original rhs entry splitConstrainedRhsEntryOneSide(rhs, constrdInd, vConstrgInd); + rhs[constrdInd] = u[constrdInd]; // set identity row (or interpolating row in linear case) - setConstrainedRow(mat, constrdInd, vConstrgInd, m_bAssembleLinearProblem); + setConstrainedRow(mat, constrdInd, vConstrgInd); } // apply chain rule, i.e., also split constrained columns @@ -1086,8 +1078,6 @@ adjust_prolongation number time ) { - if (m_bAssembleLinearProblem) return; - if (this->m_spAssTuner->single_index_assembling_enabled()) UG_THROW("index-wise assemble routine is not " "implemented for SymP1Constraints \n"); diff --git a/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h b/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h index e362a5af2..6580fc5c3 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h +++ b/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h @@ -203,7 +203,7 @@ class DirichletBoundary number time = 0.0); /// sets unity rows in A and dirichlet values in right-hand side b - void adjust_linear(matrix_type& A, vector_type& b, + void adjust_linear(matrix_type& A, vector_type& b, const vector_type& u, ConstSmartPtr dd, int type, number time = 0.0); /// sets the dirichlet value in the right-hand side diff --git a/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary_impl.h b/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary_impl.h index 5b10bffdf..6cb943944 100644 --- a/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary_impl.h +++ b/ugbase/lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary_impl.h @@ -1272,7 +1272,7 @@ adjust_correction(const std::vector& vUserData, int si, template void DirichletBoundary:: -adjust_linear(matrix_type& A, vector_type& b, +adjust_linear(matrix_type& A, vector_type& b, const vector_type& u, ConstSmartPtr dd, int type, number time) { extract_data(); diff --git a/ugbase/lib_disc/spatial_disc/domain_disc.h b/ugbase/lib_disc/spatial_disc/domain_disc.h index 01adacb98..b2f2cad82 100644 --- a/ugbase/lib_disc/spatial_disc/domain_disc.h +++ b/ugbase/lib_disc/spatial_disc/domain_disc.h @@ -129,9 +129,9 @@ class DomainDiscretizationBase {assemble_defect(d, u, dd(gl));} /// \copydoc IAssemble::assemble_linear() - virtual void assemble_linear(matrix_type& A, vector_type& b, ConstSmartPtr dd); - virtual void assemble_linear(matrix_type& mat, vector_type& rhs, const GridLevel& gl) - {assemble_linear(mat, rhs, dd(gl));} + virtual void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, ConstSmartPtr dd); + virtual void assemble_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, const GridLevel& gl) + {assemble_linear(mat, rhs, u, dd(gl));} /// \copydoc IAssemble::assemble_rhs() virtual void assemble_rhs(vector_type& rhs, const vector_type& u, ConstSmartPtr dd); @@ -158,7 +158,14 @@ class DomainDiscretizationBase {assemble_defect(d, u, u.dof_distribution());} void assemble_linear(matrix_type& A, GridFunction& rhs) - {assemble_linear(A, rhs, rhs.dof_distribution());} + { + UG_LOGN("HINT: assemble_linear without passing a solution is deprecated.\n" + " Please provide the solution as it will enter the solver as third argument."); + assemble_linear(A, rhs, rhs, rhs.dof_distribution()); + } + + void assemble_linear(matrix_type& A, GridFunction& rhs, const GridFunction& u) + {assemble_linear(A, rhs, u, rhs.dof_distribution());} void assemble_rhs(vector_type& rhs, GridFunction& u) {assemble_rhs(rhs, u, u.dof_distribution());} diff --git a/ugbase/lib_disc/spatial_disc/domain_disc_impl.h b/ugbase/lib_disc/spatial_disc/domain_disc_impl.h index ebd497789..13dd39771 100644 --- a/ugbase/lib_disc/spatial_disc/domain_disc_impl.h +++ b/ugbase/lib_disc/spatial_disc/domain_disc_impl.h @@ -834,7 +834,7 @@ AssembleDefect( const std::vector*>& vElemDisc, /////////////////////////////////////////////////////////////////////////////// template void DomainDiscretizationBase:: -assemble_linear(matrix_type& mat, vector_type& rhs, +assemble_linear(matrix_type& mat, vector_type& rhs, const vector_type& u, ConstSmartPtr dd) { PROFILE_FUNC_GROUP("discretization"); @@ -933,7 +933,7 @@ assemble_linear(matrix_type& mat, vector_type& rhs, if(m_vConstraint[i]->type() & type) { m_vConstraint[i]->set_ass_tuner(m_spAssTuner); - m_vConstraint[i]->adjust_linear(mat, rhs, dd, type); + m_vConstraint[i]->adjust_linear(mat, rhs, u, dd, type); } } post_assemble_loop(m_vElemDisc); @@ -1839,7 +1839,7 @@ assemble_linear(matrix_type& mat, vector_type& rhs, if(m_vConstraint[i]->type() & type) { m_vConstraint[i]->set_ass_tuner(m_spAssTuner); - m_vConstraint[i]->adjust_linear(mat, rhs, dd, type, vSol->time(0)); + m_vConstraint[i]->adjust_linear(mat, rhs, *vSol->latest(), dd, type, vSol->time(0)); } } } UG_CATCH_THROW("Cannot adjust linear."); diff --git a/ugbase/lib_disc/spatial_disc/domain_disc_interface.h b/ugbase/lib_disc/spatial_disc/domain_disc_interface.h index f944863ed..b0d646fff 100644 --- a/ugbase/lib_disc/spatial_disc/domain_disc_interface.h +++ b/ugbase/lib_disc/spatial_disc/domain_disc_interface.h @@ -178,8 +178,8 @@ class IDomainDiscretization : public IAssemble, public IDomainErrorInd * \param[out] b Right-Hand-Side * \param[in] dd DoF Distribution */ - virtual void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) = 0; - virtual void assemble_linear(matrix_type& A, vector_type& b, ConstSmartPtr dd) = 0; + virtual void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl) = 0; + virtual void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, ConstSmartPtr dd) = 0; /// assembles the rhs virtual void assemble_rhs(vector_type& rhs, const vector_type& u, const GridLevel& gl) = 0; diff --git a/ugbase/lib_disc/time_disc/composite_time_disc.h b/ugbase/lib_disc/time_disc/composite_time_disc.h index 20608e515..ef0948c5f 100644 --- a/ugbase/lib_disc/time_disc/composite_time_disc.h +++ b/ugbase/lib_disc/time_disc/composite_time_disc.h @@ -126,7 +126,7 @@ class CompositeTimeDiscretization void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl); /// @copydoc IAssemble::assemble_linear - void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl); + void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl); /// @copydoc IAssemble::assemble_rhs(vector_type&, vector_type&, GridLevel&) void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl); diff --git a/ugbase/lib_disc/time_disc/composite_time_disc_impl.h b/ugbase/lib_disc/time_disc/composite_time_disc_impl.h index c4f001e40..ee9036d1e 100644 --- a/ugbase/lib_disc/time_disc/composite_time_disc_impl.h +++ b/ugbase/lib_disc/time_disc/composite_time_disc_impl.h @@ -157,19 +157,20 @@ void CompositeTimeDiscretization::assemble_linear ( matrix_type& A, vector_type& b, + const vector_type& u, const GridLevel& gl ) { UG_COND_THROW(!m_vTimeDisc.size(), "At least one time disc must be added to CompositeTimeDiscretization.") - m_vTimeDisc[0]->assemble_linear(A, b, gl); + m_vTimeDisc[0]->assemble_linear(A, b, u, gl); for (size_t i = 1; i < m_vTimeDisc.size(); ++i) { // avoid clearing of the matrix before assembling m_vTimeDisc[i]->domain_disc()->ass_tuner()->disable_clear_on_resize(); - m_vTimeDisc[i]->assemble_linear(A, b, gl); + m_vTimeDisc[i]->assemble_linear(A, b, u, gl); } } diff --git a/ugbase/lib_disc/time_disc/theta_time_step.h b/ugbase/lib_disc/time_disc/theta_time_step.h index ee2be2bd8..567626b45 100644 --- a/ugbase/lib_disc/time_disc/theta_time_step.h +++ b/ugbase/lib_disc/time_disc/theta_time_step.h @@ -105,7 +105,7 @@ class MultiStepTimeDiscretization void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl); - void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl); + void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl); void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl); @@ -478,7 +478,7 @@ class SDIRK void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl); - void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl); + void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl); void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl); diff --git a/ugbase/lib_disc/time_disc/theta_time_step_impl.h b/ugbase/lib_disc/time_disc/theta_time_step_impl.h index c5b3e2000..68788356a 100644 --- a/ugbase/lib_disc/time_disc/theta_time_step_impl.h +++ b/ugbase/lib_disc/time_disc/theta_time_step_impl.h @@ -181,7 +181,7 @@ adjust_solution(vector_type& u, const GridLevel& gl) template void MultiStepTimeDiscretization:: -assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) +assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl) { PROFILE_BEGIN_GROUP(MultiStepTimeDiscretization_assemble_linear, "discretization MultiStepTimeDiscretization"); // perform checks @@ -705,7 +705,7 @@ adjust_solution(vector_type& u, const GridLevel& gl) template void SDIRK:: -assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) +assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl) { UG_THROW("Not implemented") }