3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
croperator.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_CROPERATOR2P_HH
25#define DUMUX_CROPERATOR2P_HH
26
27#include<iostream>
28#include<vector>
29#include<set>
30#include<map>
31#include<cassert>
32
33#include<dune/common/timer.hh>
34#include<dune/common/fvector.hh>
35#include<dune/common/fmatrix.hh>
36#include<dune/common/exceptions.hh>
37#include<dune/geometry/type.hh>
38#include<dune/grid/common/grid.hh>
39#include<dune/grid/common/mcmgmapper.hh>
40#include<dune/istl/bvector.hh>
41#include<dune/istl/operators.hh>
42#include<dune/istl/bcrsmatrix.hh>
43
46#include "localstiffness.hh"
47
48namespace Dumux {
49
64template<class TypeTag>
66{
67 template<int dim>
68 struct FaceLayout
69 {
70 bool contains (Dune::GeometryType gt)
71 {
72 return gt.dim() == dim-1;
73 }
74 };
77 enum {dim=GridView::dimension};
78 using IS = typename GridView::IndexSet;
79 using BlockType = Dune::FieldMatrix<Scalar, 1, 1>;
81 using MBlockType = typename MatrixType::block_type;
82 using rowiterator = typename MatrixType::RowIterator;
83 using coliterator = typename MatrixType::ColIterator;
84 using BCBlockType = std::array<BoundaryConditions::Flags, 1>; // componentwise boundary conditions
85 using SatType = Dune::BlockVector< Dune::FieldVector<double, 1> >;
86 using FaceMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
87
89 enum
90 {
91 pressureEqIdx = Indices::pressureEqIdx,
92 };
93
95 int nnz()
96 {
97 return (4*dim - 1)*size_;
98 }
99
100public:
102
103 CROperatorAssemblerTwoP (const GridView& gridview)
104 : gridView_(gridview)
105 , faceMapper_(gridView_, Dune::mcmgLayout(Dune::Codim<1>()))
106 , size_(faceMapper_.size())
107 , A_(size_, size_, nnz()
108 , RepresentationType::random)
109 {}
110
113 {
114 faceMapper_.update();
115
116 size_ = faceMapper_.size();
117 A_.setSize(size_, size_, nnz());
118
119 assert(nnz() != 0);
120
121 // set size of all rows to zero
122 for (unsigned int i = 0; i < size_; i++)
123 A_.setrowsize(i,0);
124
125 // build needs a flag for all entities of all codims
126 std::vector<bool> visited(size_, false);
127
128 // LOOP 1 : Compute row sizes
129 for (const auto& element : elements(gridView_))
130 {
131 int numFaces = element.subEntities(1);
132
133 for (int i = 0; i < numFaces; i++)
134 {
135 int index = faceMapper_.subIndex(element, i,1);
136
137 if (!visited[index])
138 {
139 A_.incrementrowsize(index);
140 visited[index] = true;
141 // std::cout << "increment row " << index << std::endl;
142 }
143 A_.incrementrowsize(index, numFaces - 1);
144 // std::cout << "increment row " << index
145 // << " by " << numFaces - 1 << std::endl;
146 }
147 }
148
149 // now the row sizes have been set
150 A_.endrowsizes();
151
152 // clear the flags for the next round, actually that is not necessary because addindex takes care of this
153 visited.assign(size_, false);
154
155 // LOOP 2 : insert the nonzeros
156 for (const auto& element : elements(gridView_))
157 {
158 int numFaces = element.subEntities(1);
159
160 for (int i = 0; i < numFaces; i++)
161 {
162 int indexI = faceMapper_.subIndex(element, i, 1);
163
164 if (!visited[indexI])
165 {
166 A_.addindex(indexI,indexI);
167 visited[indexI] = true;
168 }
169 for (int k = 0; k < numFaces; k++)
170 if (k != i) {
171 int indexJ = faceMapper_.subIndex(element, k, 1);
172
173 A_.addindex(indexI, indexJ);
174 //std::cout << "indexI = " << indexI << ", added indexJ = " << indexJ << std::endl;
175 }
176 }
177 }
178
179 // now the matrix is ready for use
180 A_.endindices();
181 }
182
185 {
186 return A_;
187 }
188
191 {
192 return A_;
193 }
194
195 const FaceMapper& faceMapper()
196 {
197 return faceMapper_;
198 }
199
200 const FaceMapper& faceMapper() const
201 {
202 return faceMapper_;
203 }
204
222 template <class LocalStiffness, class Vector>
223 void assemble (LocalStiffness& loc, Vector& u, Vector& f)
224 {
225
226 // check size
227 if (u.N()!=A_.M() || f.N()!=A_.N())
228 DUNE_THROW(Dune::MathError,"CROperatorAssemblerTwoP::assemble(): size mismatch");
229 // clear global stiffness matrix and right hand side
230 A_ = 0;
231 f = 0;
232
233 // allocate flag vector to hold flags for essential boundary conditions
234 std::vector<BCBlockType> essential(faceMapper_.size());
235 for (typename std::vector<BCBlockType>::size_type i=0; i<essential.size(); i++)
236 essential[i][0] = BoundaryConditions::neumann;
237
238 // local to global id mapping (do not ask vertex mapper repeatedly
239 Dune::FieldVector<int, 2*dim> local2Global(0);
240
241 // run over all leaf elements
242 for (const auto& element : elements(gridView_))
243 {
244 unsigned int numFaces = element.subEntities(1);
245
246 // get local to global id map
247 for (unsigned int k = 0; k < numFaces; k++)
248 {
249 int alpha = faceMapper_.subIndex(element, k, 1);
250 local2Global[k] = alpha;
251 }
252
253 // build local stiffness matrix for CR elements
254 // inludes rhs and boundary condition information
255 loc.assemble(element, 1); // assemble local stiffness matrix
256
257
258 // accumulate local matrix into global matrix for non-hanging nodes
259 for (unsigned int i=0; i<numFaces; i++) // loop over rows, i.e. test functions
260 {
261 // accumulate matrix
262 for (unsigned int j=0; j<numFaces; j++)
263 {
264 // the standard entry
265 A_[local2Global[i]][local2Global[j]] += loc.mat(i,j);
266 }
267
268 // essential boundary condition and rhs
269 if (loc.bc(i).isDirichlet(pressureEqIdx))
270 {
271 essential[local2Global[i]][0] = BoundaryConditions::dirichlet;
272 f[local2Global[i]][0] = loc.rhs(i)[0];
273 }
274 else
275 f[local2Global[i]][0] += loc.rhs(i)[0];
276 }
277 }
278 // run over all leaf elements
279 for (const auto& element : elements(gridView_))
280 {
281 unsigned int numFaces = element.subEntities(1);
282
283 // get local to global id map
284 for (unsigned int k = 0; k < numFaces; k++)
285 {
286 int alpha = faceMapper_.subIndex(element, k, 1);
287 local2Global[k] = alpha;
288 }
289 loc.completeRHS(element, local2Global, f);
290 }
291
292 // put in essential boundary conditions
293 rowiterator endi=A_.end();
294 for (rowiterator i=A_.begin(); i!=endi; ++i)
295 {
296 // muck up extra rows
297 if ((int) i.index() >= (int) faceMapper_.size())
298 {
299 coliterator endj=(*i).end();
300 for (coliterator j=(*i).begin(); j!=endj; ++j)
301 {
302 (*j) = 0;
303 if (j.index()==i.index())
304 (*j)[0][0] = 1;
305 }
306 f[i.index()] = 0;
307 }
308 // insert dirichlet and processor boundary conditions
309 else if (essential[i.index()][0]!=BoundaryConditions::neumann)
310 {
311 coliterator endj=(*i).end();
312 for (coliterator j=(*i).begin(); j!=endj; ++j)
313 if (j.index()==i.index())
314 {
315 (*j)[0][0] = 1;
316 }
317 else
318 {
319 (*j)[0][0] = 0;
320 }
321 u[i.index()][0] = f[i.index()][0];
322 }
323 }
324 }
325
326protected:
327 const GridView gridView_;
328 FaceMapper faceMapper_;
329 unsigned int size_;
331};
332
333} // end namespace Dumux
334
335#endif
Definition of boundary condition types, extend if necessary.
Base class for assembling local stiffness matrices.
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Definition: common/pdesolver.hh:35
@ dirichlet
Dirichlet boundary.
Definition: boundaryconditions.hh:41
@ neumann
Neumann boundary.
Definition: boundaryconditions.hh:39
Definition: matrix.hh:32
Extends CROperatorBase by a generic methods to assemble global stiffness matrix from local stiffness ...
Definition: croperator.hh:66
const GridView gridView_
Definition: croperator.hh:327
CROperatorAssemblerTwoP(const GridView &gridview)
Definition: croperator.hh:103
FaceMapper faceMapper_
Definition: croperator.hh:328
unsigned int size_
Definition: croperator.hh:329
const FaceMapper & faceMapper() const
Definition: croperator.hh:200
const RepresentationType & operator*() const
Returns const reference to operator matrix.
Definition: croperator.hh:184
const FaceMapper & faceMapper()
Definition: croperator.hh:195
void assemble(LocalStiffness &loc, Vector &u, Vector &f)
Assembles global stiffness matrix.
Definition: croperator.hh:223
void initialize()
Initialize the CR operator assembler.
Definition: croperator.hh:112
RepresentationType A_
Definition: croperator.hh:330
Base class for local assemblers.
Definition: localstiffness.hh:60
const VBlockType & rhs(int i) const
Accesses right hand side.
Definition: localstiffness.hh:210
const MBlockType & mat(int i, int j) const
Accesses the local stiffness matrix.
Definition: localstiffness.hh:197
const BoundaryTypes & bc(int i) const
Accesses boundary condition for each dof.
Definition: localstiffness.hh:223
virtual void assemble(const Entity &e, int k=1)=0
Assembles local stiffness matrix including boundary conditions for given element and order.
Base file for properties related to sequential IMPET algorithms.