Commit 266db842 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[freeflow][md] Add coupling jacobian pattern for fcstaggered

parent 1bc99753
......@@ -169,6 +169,51 @@ Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingM
return pattern;
* \ingroup MultiDomain
* \brief Helper function to generate coupling Jacobian pattern (off-diagonal blocks)
* for the staggered scheme (degrees of freedom on cell centers)
template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
typename std::enable_if_t<(GridGeometryI::discMethod == DiscretizationMethod::fcstaggered), int> = 0>
Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
Dune::index_constant<i> domainI,
const GridGeometryI& gridGeometryI,
Dune::index_constant<j> domainJ,
const GridGeometryJ& gridGeometryJ)
Dune::MatrixIndexSet pattern(gridGeometryI.numDofs(), gridGeometryJ.numDofs());
auto fvGeometry = localView(gridGeometryI);
for (const auto& elementI : elements(gridGeometryI.gridView()))
for (const auto& scv : scvs(fvGeometry))
const auto globalI = scv.dofIndex();
const auto& stencil = couplingManager.couplingStencil(domainI, elementI, scv, domainJ);
for (const auto globalJ : stencil)
assert(globalJ < gridGeometryJ.numDofs());
pattern.add(globalI, globalJ);
if (gridGeometryI.isPeriodic())
if (gridGeometryI.dofOnPeriodicBoundary(globalI))
const auto globalIP = gridGeometryI.periodicallyMappedDof(globalI);
if (globalI > globalIP)
pattern.add(globalIP, globalJ);
return pattern;
} // end namespace Dumux
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment