3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
jacobianpattern.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_JACOBIAN_PATTERN_HH
25#define DUMUX_JACOBIAN_PATTERN_HH
26
27#include <type_traits>
28#include <dune/istl/matrixindexset.hh>
30
31namespace Dumux {
32
37template<bool isImplicit, class GridGeometry,
38 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethod::box), int> = 0>
39Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
40{
41 const auto numDofs = gridGeometry.numDofs();
42 Dune::MatrixIndexSet pattern;
43 pattern.resize(numDofs, numDofs);
44
45 // matrix pattern for implicit Jacobians
46 if (isImplicit)
47 {
48 static constexpr int dim = std::decay_t<decltype(gridGeometry.gridView())>::dimension;
49 for (const auto& element : elements(gridGeometry.gridView()))
50 {
51 for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
52 {
53 const auto globalI = gridGeometry.vertexMapper().subIndex(element, vIdx, dim);
54 for (unsigned int vIdx2 = 0; vIdx2 < element.subEntities(dim); ++vIdx2)
55 {
56 const auto globalJ = gridGeometry.vertexMapper().subIndex(element, vIdx2, dim);
57 pattern.add(globalI, globalJ);
58
59 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
60 {
61 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
62 pattern.add(globalIP, globalI);
63 pattern.add(globalI, globalIP);
64 if (globalI > globalIP)
65 pattern.add(globalIP, globalJ);
66 }
67 }
68 }
69 }
70 }
71
72 // matrix pattern for explicit Jacobians -> diagonal matrix
73 else
74 {
75 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
76 pattern.add(globalI, globalI);
77 }
78
79 return pattern;
80}
81
86template<bool isImplicit, class GridGeometry,
87 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethod::cctpfa)
88 || (GridGeometry::discMethod == DiscretizationMethod::ccmpfa) ), int> = 0>
89Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
90{
91 const auto numDofs = gridGeometry.numDofs();
92 Dune::MatrixIndexSet pattern;
93 pattern.resize(numDofs, numDofs);
94
95 // matrix pattern for implicit Jacobians
96 if (isImplicit)
97 {
98 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
99 {
100 pattern.add(globalI, globalI);
101 for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
102 pattern.add(dataJ.globalJ, globalI);
103 }
104 }
105
106 // matrix pattern for explicit Jacobians -> diagonal matrix
107 else
108 {
109 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
110 pattern.add(globalI, globalI);
111 }
112
113 return pattern;
114}
115
120template<bool isImplicit, class GridGeometry,
121 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethod::staggered) ), int> = 0>
122auto getJacobianPattern(const GridGeometry& gridGeometry)
123{
124 // resize the jacobian and the residual
125 const auto numDofs = gridGeometry.numDofs();
126 Dune::MatrixIndexSet pattern(numDofs, numDofs);
127
128 const auto& connectivityMap = gridGeometry.connectivityMap();
129
130 // evaluate the acutal pattern
131 for (const auto& element : elements(gridGeometry.gridView()))
132 {
133 if(gridGeometry.isCellCenter())
134 {
135 // the global index of the element at hand
136 static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
137 const auto ccGlobalI = gridGeometry.elementMapper().index(element);
138 pattern.add(ccGlobalI, ccGlobalI);
139 for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
140 pattern.add(ccGlobalI, ccGlobalJ);
141 }
142 else
143 {
144 static constexpr auto faceIdx = GridGeometry::faceIdx();
145 auto fvGeometry = localView(gridGeometry);
146 fvGeometry.bindElement(element);
147
148 // loop over sub control faces
149 for (auto&& scvf : scvfs(fvGeometry))
150 {
151 const auto faceGlobalI = scvf.dofIndex();
152 pattern.add(faceGlobalI, faceGlobalI);
153 for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
154 pattern.add(faceGlobalI, faceGlobalJ);
155 }
156 }
157 }
158
159 return pattern;
160}
161
166template<class FEBasis>
167Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis& feBasis)
168{
169 const auto numDofs = feBasis.size();
170
171 Dune::MatrixIndexSet pattern;
172 pattern.resize(numDofs, numDofs);
173
174 // matrix pattern for implicit Jacobians
175 for (const auto& element : elements(feBasis.gridView()))
176 {
177 auto localView = feBasis.localView();
178 localView.bind(element);
179
180 const auto& finiteElement = localView.tree().finiteElement();
181 const auto numLocalDofs = finiteElement.localBasis().size();
182 for (size_t i = 0; i < numLocalDofs; i++)
183 {
184 const auto dofIdxI = localView.index(i);
185 for (size_t j = 0; j < numLocalDofs; j++)
186 {
187 const auto dofIdxJ = localView.index(j);
188 pattern.add(dofIdxI, dofIdxJ);
189 }
190 }
191 }
192
193 return pattern;
194}
195
202template<bool isImplicit, class GridGeometry,
203 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethod::fem), int> = 0>
204Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
205{ return getFEJacobianPattern(gridGeometry.feBasis()); }
206
207} // namespace Dumux
208
209#endif
The available discretization methods in Dumux.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Dune::MatrixIndexSet getJacobianPattern(const GridGeometry &gridGeometry)
Helper function to generate Jacobian pattern for the box method.
Definition: jacobianpattern.hh:39
Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis &feBasis)
Helper function to generate Jacobian pattern for finite element scheme.
Definition: jacobianpattern.hh:167
Definition: adapt.hh:29