diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh index f2340daee9c997fc9e651e5c6ff04e466dcd6cba..61e7dee1271d8d39729f82abb655c97645da12bc 100644 --- a/dumux/assembly/fvassembler.hh +++ b/dumux/assembly/fvassembler.hh @@ -474,8 +474,27 @@ private: // enforce constraint in second row res[m.second] = curSol[m.second] - curSol[m.first]; + + // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first + auto setMatrixBlock = [] (auto& matrixBlock, double diagValue) + { + for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx) + matrixBlock[eIdx][eIdx] = diagValue; + }; + for (auto it = jac[m.second].begin(); it != end; ++it) - (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0; + { + auto& matrixBlock = *it; + matrixBlock = 0.0; + + assert(matrixBlock.N() == matrixBlock.M()); + if(it.index() == m.second) + setMatrixBlock(matrixBlock, 1.0); + + if(it.index() == m.first) + setMatrixBlock(matrixBlock, -1.0); + + } } } } diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh index d619eba11fb582f843cf271c32a461b85de1afb7..f7664cddab59882c95510bdabd71f1eb1f4ea3e2 100644 --- a/dumux/multidomain/fvassembler.hh +++ b/dumux/multidomain/fvassembler.hh @@ -611,16 +611,34 @@ private: // add the second row to the first res[m.first] += res[m.second]; - // enforce the solution of the first periodic DOF to the second one - res[m.second] = curSol[m.second] - curSol[m.first]; - const auto end = jac[m.second].end(); for (auto it = jac[m.second].begin(); it != end; ++it) jac[m.first][it.index()] += (*it); - // enforce constraint in second row + + // enforce the solution of the first periodic DOF to the second one + res[m.second] = curSol[m.second] - curSol[m.first]; + + // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first + auto setMatrixBlock = [] (auto& matrixBlock, double diagValue) + { + for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx) + matrixBlock[eIdx][eIdx] = diagValue; + }; + for (auto it = jac[m.second].begin(); it != end; ++it) - (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0; + { + auto& matrixBlock = *it; + matrixBlock = 0.0; + + assert(matrixBlock.N() == matrixBlock.M()); + if(it.index() == m.second) + setMatrixBlock(matrixBlock, 1.0); + + if(it.index() == m.first) + setMatrixBlock(matrixBlock, -1.0); + + } using namespace Dune::Hybrid; forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)