version 3.10-dev
couplingjacobianpattern.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
12#ifndef DUMUX_MUTLIDOMAIN_COUPLING_JACOBIAN_PATTERN_HH
13#define DUMUX_MUTLIDOMAIN_COUPLING_JACOBIAN_PATTERN_HH
14
15#include <type_traits>
16#include <dune/common/indices.hh>
17#include <dune/istl/matrixindexset.hh>
19
20namespace Dumux {
21
27template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
28 typename std::enable_if_t<( (GridGeometryI::discMethod == DiscretizationMethods::cctpfa)
29 || (GridGeometryI::discMethod == DiscretizationMethods::ccmpfa) ), int> = 0>
30Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
31 Dune::index_constant<i> domainI,
32 const GridGeometryI& gridGeometryI,
33 Dune::index_constant<j> domainJ,
34 const GridGeometryJ& gridGeometryJ)
35{
36 const auto numDofsI = gridGeometryI.numDofs();
37 const auto numDofsJ = gridGeometryJ.numDofs();
38 Dune::MatrixIndexSet pattern;
39 pattern.resize(numDofsI, numDofsJ);
40
41 // matrix pattern for implicit Jacobians
42 if (isImplicit)
43 {
44 for (const auto& elementI : elements(gridGeometryI.gridView()))
45 {
46 const auto& stencil = couplingManager.couplingStencil(domainI, elementI, domainJ);
47 const auto globalI = gridGeometryI.elementMapper().index(elementI);
48 for (const auto globalJ : stencil)
49 pattern.add(globalI, globalJ);
50 }
51 }
52
53 // matrix pattern for explicit Jacobians
54 // -> diagonal matrix, so coupling block is empty
55 // just return the empty pattern
56
57 return pattern;
58}
59
65template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
66 typename std::enable_if_t<(GridGeometryI::discMethod == DiscretizationMethods::staggered &&
67 GridGeometryI::isCellCenter()), int> = 0>
68Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
69 Dune::index_constant<i> domainI,
70 const GridGeometryI& gridGeometryI,
71 Dune::index_constant<j> domainJ,
72 const GridGeometryJ& gridGeometryJ)
73{
74 Dune::MatrixIndexSet pattern(gridGeometryI.numDofs(), gridGeometryJ.numDofs());
75
76 for (const auto& elementI : elements(gridGeometryI.gridView()))
77 {
78 const auto ccGlobalI = gridGeometryI.elementMapper().index(elementI);
79 for (auto&& faceGlobalJ : couplingManager.couplingStencil(domainI, elementI, domainJ))
80 pattern.add(ccGlobalI, faceGlobalJ);
81 }
82
83 return pattern;
84}
85
91template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
92 typename std::enable_if_t<(GridGeometryI::discMethod == DiscretizationMethods::staggered &&
93 GridGeometryI::isFace()), int> = 0>
94Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
95 Dune::index_constant<i> domainI,
96 const GridGeometryI& gridGeometryI,
97 Dune::index_constant<j> domainJ,
98 const GridGeometryJ& gridGeometryJ)
99{
100 Dune::MatrixIndexSet pattern(gridGeometryI.numDofs(), gridGeometryJ.numDofs());
101
102 auto fvGeometry = localView(gridGeometryI);
103 for (const auto& elementI : elements(gridGeometryI.gridView()))
104 {
105 fvGeometry.bindElement(elementI);
106
107 // loop over sub control faces
108 for (auto&& scvf : scvfs(fvGeometry))
109 {
110 const auto faceGlobalI = scvf.dofIndex();
111 for (auto&& globalJ : couplingManager.couplingStencil(domainI, scvf, domainJ))
112 pattern.add(faceGlobalI, globalJ);
113 }
114 }
115
116 return pattern;
117}
118
124template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
125 typename std::enable_if_t<(GridGeometryI::discMethod == DiscretizationMethods::fcstaggered), int> = 0>
126Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
127 Dune::index_constant<i> domainI,
128 const GridGeometryI& gridGeometryI,
129 Dune::index_constant<j> domainJ,
130 const GridGeometryJ& gridGeometryJ)
131{
132 Dune::MatrixIndexSet pattern(gridGeometryI.numDofs(), gridGeometryJ.numDofs());
133
134 auto fvGeometry = localView(gridGeometryI);
135 for (const auto& elementI : elements(gridGeometryI.gridView()))
136 {
137 fvGeometry.bindElement(elementI);
138 for (const auto& scv : scvs(fvGeometry))
139 {
140 const auto globalI = scv.dofIndex();
141 const auto& stencil = couplingManager.couplingStencil(domainI, elementI, scv, domainJ);
142 for (const auto globalJ : stencil)
143 {
144 assert(globalJ < gridGeometryJ.numDofs());
145 pattern.add(globalI, globalJ);
146
147 if (gridGeometryI.isPeriodic())
148 {
149 if (gridGeometryI.dofOnPeriodicBoundary(globalI))
150 {
151 const auto globalIP = gridGeometryI.periodicallyMappedDof(globalI);
152
153 if (globalI > globalIP)
154 pattern.add(globalIP, globalJ);
155 }
156 }
157 }
158 }
159 }
160
161 return pattern;
162}
163
169template<bool isImplicit, class CouplingManager, class GridGeometryI, class GridGeometryJ, std::size_t i, std::size_t j,
170 typename std::enable_if_t<DiscretizationMethods::isCVFE<typename GridGeometryI::DiscretizationMethod>, int> = 0>
171Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager& couplingManager,
172 Dune::index_constant<i> domainI,
173 const GridGeometryI& gridGeometryI,
174 Dune::index_constant<j> domainJ,
175 const GridGeometryJ& gridGeometryJ)
176{
177 Dune::MatrixIndexSet pattern;
178
179 // matrix pattern for implicit Jacobians
180 if constexpr (isImplicit)
181 {
182 pattern.resize(gridGeometryI.numDofs(), gridGeometryJ.numDofs());
183 auto fvGeometry = localView(gridGeometryI);
184 for (const auto& elementI : elements(gridGeometryI.gridView()))
185 {
186 fvGeometry.bindElement(elementI);
187 const auto& stencil = couplingManager.couplingStencil(domainI, elementI, domainJ);
188 for (const auto& scv : scvs(fvGeometry))
189 {
190 for (const auto globalJ : stencil)
191 pattern.add(scv.dofIndex(), globalJ);
192
193 }
194 }
195 }
196
197 // matrix pattern for explicit Jacobians
198 // -> diagonal matrix, so coupling block is empty
199 // just return the empty pattern
200
201 return pattern;
202}
203
204} // end namespace Dumux
205
206#endif
The interface of the coupling manager for multi domain problems.
Definition: multidomain/couplingmanager.hh:37
const CouplingStencilType< i, j > & couplingStencil(Dune::index_constant< i > domainI, const Element< i > &elementI, Dune::index_constant< j > domainJ) const
returns an iterable container of all indices of degrees of freedom of domain j that couple with / inf...
Definition: multidomain/couplingmanager.hh:103
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:26
Dune::MatrixIndexSet getCouplingJacobianPattern(const CouplingManager &couplingManager, Dune::index_constant< i > domainI, const GridGeometryI &gridGeometryI, Dune::index_constant< j > domainJ, const GridGeometryJ &gridGeometryJ)
Helper function to generate coupling Jacobian pattern (off-diagonal blocks) for cell-centered schemes...
Definition: couplingjacobianpattern.hh:30
The available discretization methods in Dumux.
constexpr CCMpfa ccmpfa
Definition: method.hh:146
constexpr CCTpfa cctpfa
Definition: method.hh:145
constexpr Staggered staggered
Definition: method.hh:149
constexpr FCStaggered fcstaggered
Definition: method.hh:151
Definition: adapt.hh:17