3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
localassemblerhelper.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 *****************************************************************************/
26#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH
27#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH
28
29#include <algorithm>
30#include <vector>
31#include <cassert>
32#include <utility>
33#include <type_traits>
34
36
37namespace Dumux {
38
46{
47 // Helper structs to detect if matrix has resize function
48 struct HasMatrixResize
49 {
50 template<class M>
51 auto operator()(const M& m) -> decltype(std::declval<M>().resize(0, 0))
52 {}
53 };
54
55 // Helper structs to detect if vector has resize function
56 struct HasVectorResize
57 {
58 template<class V>
59 auto operator()(const V& v) -> decltype(std::declval<V>().resize(0))
60 {}
61 };
62
63 template<class Matrix>
64 static constexpr auto matrixHasResizeFunction()
65 { return decltype( isValid(HasMatrixResize())(std::declval<Matrix>()) )::value; }
66
67 template<class Vector>
68 static constexpr auto vectorHasResizeFunction()
69 { return decltype( isValid(HasVectorResize())(std::declval<Vector>()) )::value; }
70
71public:
81 template< class FVElementGeometry, class DataHandle, class IV >
82 static void solveLocalSystem(const FVElementGeometry& fvGeometry,
83 DataHandle& handle,
84 IV& iv)
85 {
86 assert(iv.numUnknowns() > 0);
87
88 // T = C*(A^-1)*B + D
89 handle.A().invert();
90 handle.CA().rightmultiply(handle.A());
91 handle.T() += multiplyMatrices(handle.CA(), handle.AB());
92 handle.AB().leftmultiply(handle.A());
93
94 // On surface grids, compute the "outside" transmissibilities
95 using GridView = typename IV::Traits::GridView;
96 static constexpr int dim = GridView::dimension;
97 static constexpr int dimWorld = GridView::dimensionworld;
98 if (dim < dimWorld)
99 {
100 // bring outside tij container to the right size
101 auto& tijOut = handle.tijOutside();
102 tijOut.resize(iv.numFaces());
103 using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
104 for (LocalIndexType fIdx = 0; fIdx < iv.numFaces(); ++fIdx)
105 {
106 const auto& curGlobalScvf = fvGeometry.scvf(iv.localScvf(fIdx).gridScvfIndex());
107 const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs();
108 // resize each face entry to the right number of outside faces
109 tijOut[fIdx].resize(numOutsideFaces);
110 std::for_each(tijOut[fIdx].begin(),
111 tijOut[fIdx].end(),
112 [&](auto& v) { resizeVector(v, iv.numKnowns()); });
113 }
114
115 // compute outside transmissibilities
116 for (const auto& localFaceData : iv.localFaceData())
117 {
118 // continue only for "outside" faces
119 if (!localFaceData.isOutsideFace())
120 continue;
121
122 const auto scvIdx = localFaceData.ivLocalInsideScvIndex();
123 const auto scvfIdx = localFaceData.ivLocalScvfIndex();
124 const auto idxInOut = localFaceData.scvfLocalOutsideScvfIndex();
125
126 const auto& wijk = handle.omegas()[scvfIdx][idxInOut+1];
127 auto& tij = tijOut[scvfIdx][idxInOut];
128
129 tij = 0.0;
130 for (unsigned int dir = 0; dir < dim; dir++)
131 {
132 // the scvf corresponding to this local direction in the scv
133 const auto& scvf = iv.localScvf(iv.localScv(scvIdx).localScvfIndex(dir));
134
135 // on interior faces the coefficients of the AB matrix come into play
136 if (!scvf.isDirichlet())
137 {
138 auto tmp = handle.AB()[scvf.localDofIndex()];
139 tmp *= wijk[dir];
140 tij -= tmp;
141 }
142 else
143 tij[scvf.localDofIndex()] -= wijk[dir];
144
145 // add entry from the scv unknown
146 tij[scvIdx] += wijk[dir];
147 }
148 }
149 }
150 }
151
159 template< class DataHandle, class IV >
160 static typename IV::Traits::MatVecTraits::FaceVector
161 assembleFaceUnkowns(const DataHandle& handle, const IV& iv)
162 {
163 typename IV::Traits::MatVecTraits::FaceVector u;
164 resizeVector(u, iv.numFaces());
165
166 handle.AB().mv(handle.uj(), u);
167
168 // maybe add gravity terms
169 if (handle.deltaG().size() == iv.numUnknowns())
170 handle.A().umv(handle.deltaG(), u);
171
172 return u;
173 }
174
183 template< class DataHandle, class IV >
184 static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate >
185 assembleScvGradients(const DataHandle& handle, const IV& iv)
186 {
187 const auto u = assembleFaceUnkowns(handle, iv);
188
189 using LocalScv = typename IV::Traits::LocalScvType;
190 using Gradient = typename LocalScv::GlobalCoordinate;
191
192 std::vector<Gradient> result; result.reserve(iv.numScvs());
193 for (unsigned int scvIdx = 0; scvIdx < iv.numScvs(); ++scvIdx)
194 {
195 const auto& scv = iv.localScv(scvIdx);
196
197 Gradient gradU(0.0);
198 for (unsigned int dir = 0; dir < LocalScv::myDimension; ++dir)
199 {
200 auto nu = scv.nu(dir);
201
202 // obtain face pressure
203 const auto& scvf = iv.localScvf( scv.localScvfIndex(dir) );
204 const auto faceU = !scvf.isDirichlet() ? u[scvf.localDofIndex()]
205 : handle.uj()[scvf.localDofIndex()];
206
207 nu *= faceU - handle.uj()[scv.localDofIndex()];
208 gradU += nu;
209 }
210
211 gradU /= scv.detX();
212 result.emplace_back( std::move(gradU) );
213 }
214
215 return result;
216 }
217
219 template< class Matrix,
220 class size_type,
221 std::enable_if_t<matrixHasResizeFunction<Matrix>(), int> = 0 >
222 static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
223 {
224 M.resize(rows, cols);
225 }
226
228 template< class Matrix,
229 class size_type,
230 std::enable_if_t<!matrixHasResizeFunction<Matrix>(), int> = 0 >
231 static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
232 {}
233
235 template< class Vector,
236 class size_type,
237 std::enable_if_t<vectorHasResizeFunction<Vector>(), int> = 0 >
238 static void resizeVector(Vector& v, size_type size)
239 {
240 v.resize(size);
241 }
242
244 template< class Vector,
245 class size_type,
246 std::enable_if_t<!vectorHasResizeFunction<Vector>(), int> = 0 >
247 static void resizeVector(Vector& v, size_type rows)
248 {}
249};
250
251} // end namespace Dumux
252
253#endif
A helper function for class member function introspection.
Dune::DynamicMatrix< Scalar > multiplyMatrices(const Dune::DynamicMatrix< Scalar > &M1, const Dune::DynamicMatrix< Scalar > &M2)
Multiply two dynamic matrices.
Definition: math.hh:720
Definition: adapt.hh:29
constexpr auto isValid(const Expression &t)
A function that creates a test functor to do class member introspection at compile time.
Definition: isvalid.hh:93
A class that contains helper functions as well as functionality which is common to different mpfa sch...
Definition: localassemblerhelper.hh:46
static void resizeVector(Vector &v, size_type rows)
resizes a vector to the given size (specialization for static vector type - do nothing)
Definition: localassemblerhelper.hh:247
static void solveLocalSystem(const FVElementGeometry &fvGeometry, DataHandle &handle, IV &iv)
Solves a previously assembled iv-local system of equations and stores the resulting transmissibilitie...
Definition: localassemblerhelper.hh:82
static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate > assembleScvGradients(const DataHandle &handle, const IV &iv)
Assembles the solution gradients in the sub-control volumes within an interaction volume.
Definition: localassemblerhelper.hh:185
static void resizeMatrix(Matrix &M, size_type rows, size_type cols)
resizes a matrix to the given sizes (specialization for dynamic matrix type)
Definition: localassemblerhelper.hh:222
static void resizeVector(Vector &v, size_type size)
resizes a vector to the given size (specialization for dynamic matrix type)
Definition: localassemblerhelper.hh:238
static IV::Traits::MatVecTraits::FaceVector assembleFaceUnkowns(const DataHandle &handle, const IV &iv)
Assembles the vector of face unknowns within an interaction volume.
Definition: localassemblerhelper.hh:161