3.6-git
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>
29#include <dumux/discretization/method.hh>
30
31namespace Dumux {
32
37template<bool isImplicit, class GridGeometry,
38 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::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 constexpr (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 == DiscretizationMethods::pq1bubble), int> = 0>
88Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
89{
90 const auto numDofs = gridGeometry.numDofs();
91 Dune::MatrixIndexSet pattern;
92 pattern.resize(numDofs, numDofs);
93
94 // matrix pattern for implicit Jacobians
95 if constexpr (isImplicit)
96 {
97 static constexpr int dim = std::decay_t<decltype(gridGeometry.gridView())>::dimension;
98 for (const auto& element : elements(gridGeometry.gridView()))
99 {
100 const auto eIdx = gridGeometry.dofMapper().index(element);
101 pattern.add(eIdx, eIdx);
102 for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
103 {
104 const auto globalI = gridGeometry.dofMapper().subIndex(element, vIdx, dim);
105 pattern.add(eIdx, globalI);
106 pattern.add(globalI, eIdx);
107 for (unsigned int vIdx2 = 0; vIdx2 < element.subEntities(dim); ++vIdx2)
108 {
109 const auto globalJ = gridGeometry.dofMapper().subIndex(element, vIdx2, dim);
110 pattern.add(globalI, globalJ);
111
112 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
113 {
114 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
115 pattern.add(globalIP, globalI);
116 pattern.add(globalI, globalIP);
117 if (globalI > globalIP)
118 pattern.add(globalIP, globalJ);
119 }
120 }
121 }
122 }
123 }
124
125 // matrix pattern for explicit Jacobians -> diagonal matrix
126 else
127 {
128 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
129 pattern.add(globalI, globalI);
130 }
131
132 return pattern;
133}
134
139template<bool isImplicit, class GridGeometry,
140 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethods::cctpfa)
141 || (GridGeometry::discMethod == DiscretizationMethods::ccmpfa) ), int> = 0>
142Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
143{
144 const auto numDofs = gridGeometry.numDofs();
145 Dune::MatrixIndexSet pattern;
146 pattern.resize(numDofs, numDofs);
147
148 // matrix pattern for implicit Jacobians
149 if (isImplicit)
150 {
151 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
152 {
153 pattern.add(globalI, globalI);
154 for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
155 pattern.add(dataJ.globalJ, globalI);
156 }
157 }
158
159 // matrix pattern for explicit Jacobians -> diagonal matrix
160 else
161 {
162 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
163 pattern.add(globalI, globalI);
164 }
165
166 return pattern;
167}
168
173template<bool isImplicit, class GridGeometry,
174 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethods::staggered) ), int> = 0>
175auto getJacobianPattern(const GridGeometry& gridGeometry)
176{
177 // resize the jacobian and the residual
178 const auto numDofs = gridGeometry.numDofs();
179 Dune::MatrixIndexSet pattern(numDofs, numDofs);
180
181 const auto& connectivityMap = gridGeometry.connectivityMap();
182
183 auto fvGeometry = localView(gridGeometry);
184 // evaluate the actual pattern
185 for (const auto& element : elements(gridGeometry.gridView()))
186 {
187 if(gridGeometry.isCellCenter())
188 {
189 // the global index of the element at hand
190 static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
191 const auto ccGlobalI = gridGeometry.elementMapper().index(element);
192 pattern.add(ccGlobalI, ccGlobalI);
193 for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
194 pattern.add(ccGlobalI, ccGlobalJ);
195 }
196 else
197 {
198 static constexpr auto faceIdx = GridGeometry::faceIdx();
199 fvGeometry.bindElement(element);
200
201 // loop over sub control faces
202 for (auto&& scvf : scvfs(fvGeometry))
203 {
204 const auto faceGlobalI = scvf.dofIndex();
205 pattern.add(faceGlobalI, faceGlobalI);
206 for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
207 pattern.add(faceGlobalI, faceGlobalJ);
208 }
209 }
210 }
211
212 return pattern;
213}
214
219template<class FEBasis>
220Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis& feBasis)
221{
222 const auto numDofs = feBasis.size();
223
224 Dune::MatrixIndexSet pattern;
225 pattern.resize(numDofs, numDofs);
226
227 auto localView = feBasis.localView();
228 // matrix pattern for implicit Jacobians
229 for (const auto& element : elements(feBasis.gridView()))
230 {
231 localView.bind(element);
232
233 const auto& finiteElement = localView.tree().finiteElement();
234 const auto numLocalDofs = finiteElement.localBasis().size();
235 for (std::size_t i = 0; i < numLocalDofs; i++)
236 {
237 const auto dofIdxI = localView.index(i);
238 for (std::size_t j = 0; j < numLocalDofs; j++)
239 {
240 const auto dofIdxJ = localView.index(j);
241 pattern.add(dofIdxI, dofIdxJ);
242 }
243 }
244 }
245
246 return pattern;
247}
248
255template<bool isImplicit, class GridGeometry,
256 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::fem), int> = 0>
257Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
258{ return getFEJacobianPattern(gridGeometry.feBasis()); }
259
264template<bool isImplicit, class GridGeometry,
265 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::fcstaggered), int> = 0>
266Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
267{
268 // resize the jacobian and the residual
269 const auto numDofs = gridGeometry.numDofs();
270 Dune::MatrixIndexSet pattern(numDofs, numDofs);
271
272 const auto& connectivityMap = gridGeometry.connectivityMap();
273 auto fvGeometry = localView(gridGeometry);
274
275 // set the pattern
276 for (const auto& element : elements(gridGeometry.gridView()))
277 {
278 fvGeometry.bind(element);
279 for (const auto& scv : scvs(fvGeometry))
280 {
281 const auto globalI = scv.dofIndex();
282 pattern.add(globalI, globalI);
283
284 for (const auto& scvIdxJ : connectivityMap[scv.index()])
285 {
286 const auto globalJ = fvGeometry.scv(scvIdxJ).dofIndex();
287 pattern.add(globalI, globalJ);
288
289 if (gridGeometry.isPeriodic())
290 {
291 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
292 {
293 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
294 pattern.add(globalIP, globalI);
295 pattern.add(globalI, globalIP);
296
297 if (globalI > globalIP)
298 pattern.add(globalIP, globalJ);
299 }
300 }
301 }
302 }
303 }
304
305 return pattern;
306}
307
312template<bool isImplicit, class GridGeometry,
313 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::fcdiamond), int> = 0>
314Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
315{
316 // resize the jacobian and the residual
317 const auto numDofs = gridGeometry.numDofs();
318 Dune::MatrixIndexSet pattern(numDofs, numDofs);
319
320 auto fvGeometry = localView(gridGeometry);
321 for (const auto& element : elements(gridGeometry.gridView()))
322 {
323 fvGeometry.bind(element);
324 for (const auto& scv : scvs(fvGeometry))
325 {
326 const auto globalI = scv.dofIndex();
327 for (const auto& scvJ : scvs(fvGeometry))
328 {
329 const auto globalJ = scvJ.dofIndex();
330 pattern.add(globalI, globalJ);
331
332 if (gridGeometry.isPeriodic())
333 {
334 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
335 {
336 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
337 pattern.add(globalIP, globalI);
338 pattern.add(globalI, globalIP);
339
340 if (globalI > globalIP)
341 pattern.add(globalIP, globalJ);
342 }
343 }
344 }
345 }
346 }
347
348 return pattern;
349}
350
351
352} // namespace Dumux
353
354#endif
method.hh
The available discretization methods in Dumux.
Dumux::localView
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
Dumux::getJacobianPattern
Dune::MatrixIndexSet getJacobianPattern(const GridGeometry &gridGeometry)
Helper function to generate Jacobian pattern for the box method.
Definition: jacobianpattern.hh:39
Dumux::getFEJacobianPattern
Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis &feBasis)
Helper function to generate Jacobian pattern for finite element scheme.
Definition: jacobianpattern.hh:220
Dumux
Adaption of the non-isothermal two-phase two-component flow model to problems with CO2.
Definition: adapt.hh:29
Dumux::DiscretizationMethods::ccmpfa
constexpr CCMpfa ccmpfa
Definition: method.hh:135
Dumux::DiscretizationMethods::fcdiamond
constexpr FCDiamond fcdiamond
Definition: method.hh:141
Dumux::DiscretizationMethods::cctpfa
constexpr CCTpfa cctpfa
Definition: method.hh:134
Dumux::DiscretizationMethods::box
constexpr Box box
Definition: method.hh:136
Dumux::DiscretizationMethods::staggered
constexpr Staggered staggered
Definition: method.hh:138
Dumux::DiscretizationMethods::fem
constexpr FEM fem
Definition: method.hh:139
Dumux::DiscretizationMethods::pq1bubble
constexpr PQ1Bubble pq1bubble
Definition: method.hh:137
Dumux::DiscretizationMethods::fcstaggered
constexpr FCStaggered fcstaggered
Definition: method.hh:140
Generated by doxygen 1.9.3