3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
croperatoradaptive.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_CROPERATOR2PADAPTIVE_HH
25#define DUMUX_CROPERATOR2PADAPTIVE_HH
26
27#include<iostream>
28#include<vector>
29#include<cassert>
30#include<stdio.h>
31#include<stdlib.h>
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
41#include<dune/istl/bvector.hh>
42#include<dune/istl/operators.hh>
43#include<dune/istl/bcrsmatrix.hh>
44
47#include "localstiffness.hh"
49
50namespace Dumux {
51
67template<class TypeTag>
69{
72 enum {dim=GridView::dimension};
73 using Element = typename GridView::template Codim<0>::Entity;
74 using IS = typename GridView::IndexSet;
75 using BlockType = Dune::FieldMatrix<Scalar, 1, 1>;
77 using MBlockType = typename MatrixType::block_type;
78 using rowiterator = typename MatrixType::RowIterator;
79 using coliterator = typename MatrixType::ColIterator;
80 using BCBlockType = std::array<BoundaryConditions::Flags, 1>; // componentwise boundary conditions
81 using SatType = Dune::BlockVector< Dune::FieldVector<double, 1> >;
83
85
86 enum
87 {
88 pressureEqIdx = Indices::pressureEqIdx,
89 };
90
91 // return number of rows/columns
92 int size () const
93 {
95 }
96
97public:
99
100 CROperatorAssemblerTwoPAdaptive (const GridView& gridview)
101 : gridView_(gridview), is_(gridView_.indexSet()), intersectionMapper_(gridView_)
102 {
103 A_.setBuildMode(MatrixType::random);
104 }
105
107 {
108 adapt();
109 }
110
111 void adapt()
112 {
114 A_.setSize(size(), size());
115 updateMatrix();
116 }
117
120 {
121 return A_;
122 }
123
126 {
127 return A_;
128 }
129
131 {
132 return intersectionMapper_;
133 }
134
136 {
137 return intersectionMapper_;
138 }
139
157 template <class LocalStiffness, class Vector>
158 void assemble (LocalStiffness& loc, Vector& u, Vector& f);
159
160 void updateMatrix();
161
162protected:
163 const GridView gridView_;
164 const IS& is_;
167};
168
169template<class TypeTag>
171{
172 // set size of all rows to zero
173 for (int i = 0; i < size(); i++)
174 A_.setrowsize(i, 0);
175
176 // build needs a flag for all entities of all codims
177 std::vector<bool> visited(size(), false);
178
179 int numElem = gridView_.size(0);
180
181 for (int elemIdx = 0; elemIdx < numElem; elemIdx++)
182 {
183 int numFaces = intersectionMapper_.size(elemIdx);
184 for (int fIdx = 0; fIdx < numFaces; fIdx++)
185 {
186 int faceIdxGlobal = intersectionMapper_.subIndex(elemIdx, fIdx);
187 if (!visited[faceIdxGlobal])
188 {
189 A_.incrementrowsize(faceIdxGlobal);
190 visited[faceIdxGlobal] = true;
191 }
192 for (int k = 0; k < numFaces-1; k++) {
193 A_.incrementrowsize(faceIdxGlobal);
194 }
195
196 }
197
198 }
199
200 A_.endrowsizes();
201
202 // clear the flags for the next round, actually that is not necessary because addindex takes care of this
203 for (int i = 0; i < size(); i++)
204 visited[i] = false;
205
206 for (int elemIdx = 0; elemIdx < numElem; elemIdx++)
207 {
208 int numFaces = intersectionMapper_.size(elemIdx);
209 for (int fIdx = 0; fIdx < numFaces; fIdx++)
210 {
211 int faceIdxGlobalI = intersectionMapper_.subIndex(elemIdx, fIdx);
212 if (!visited[faceIdxGlobalI])
213 {
214 A_.addindex(faceIdxGlobalI,faceIdxGlobalI);
215 visited[faceIdxGlobalI] = true;
216 }
217 for (int k = 0; k < numFaces; k++)
218 if (k != fIdx) {
219 int faceIdxGlobalJ = intersectionMapper_.subIndex(elemIdx, k);
220 A_.addindex(faceIdxGlobalI, faceIdxGlobalJ);
221 }
222
223 }
224
225 }
226
227 A_.endindices();
228}
229
230template<class TypeTag>
231template <class LocalStiffness, class Vector>
233{
234
235 // check size
236 if (u.N()!=A_.M() || f.N()!=A_.N())
237 DUNE_THROW(Dune::MathError,"CROperatorAssemblerTwoPAdaptive::assemble(): size mismatch");
238 // clear global stiffness matrix and right hand side
239 A_ = 0;
240 f = 0;
241
242 // allocate flag vector to hold flags for essential boundary conditions
243 std::vector<BCBlockType> essential(intersectionMapper_.size());
244 for (typename std::vector<BCBlockType>::size_type i=0; i<essential.size(); i++)
245 essential[i][0] = BoundaryConditions::neumann;
246
247 // local to global id mapping (do not ask vertex mapper repeatedly
248 std::vector<int> local2Global(2*dim);
249
250 // run over all leaf elements
251 for (const auto& element : elements(gridView_))
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 int eIdxGlobal = intersectionMapper_.index(element);
258
259 unsigned int numFaces = intersectionMapper_.size(eIdxGlobal);
260 local2Global.resize(numFaces);
261
262 for (unsigned int i = 0; i < numFaces; i++)
263 {
264 int idx = intersectionMapper_.subIndex(element, i);
265 local2Global[i] = idx;
266 }
267
268 // accumulate local matrix into global matrix for non-hanging nodes
269 for (unsigned int i=0; i<numFaces; i++) // loop over rows, i.e. test functions
270 {
271 // accumulate matrix
272 for (unsigned int j=0; j<numFaces; j++)
273 {
274 // the standard entry
275 A_[local2Global[i]][local2Global[j]] += loc.mat(i,j);
276 }
277
278 // essential boundary condition and rhs
279 if (loc.bc(i).isDirichlet(pressureEqIdx))
280 {
281 essential[local2Global[i]][0] = BoundaryConditions::dirichlet;
282 f[local2Global[i]][0] = loc.rhs(i)[0];
283 }
284 else
285 f[local2Global[i]][0] += loc.rhs(i)[0];
286 }
287
288 }
289 // run over all leaf elements
290 for (const auto& element : elements(gridView_))
291 {
292 int eIdxGlobal = intersectionMapper_.index(element);
293
294 unsigned int numFaces = intersectionMapper_.size(eIdxGlobal);
295 local2Global.resize(numFaces);
296
297 for (unsigned int i = 0; i < numFaces; i++)
298 {
299 int idx = intersectionMapper_.subIndex(element, i);
300 local2Global[i] = idx;
301 }
302
303 loc.completeRHS(element, local2Global, f);
304 }
305
306 // put in essential boundary conditions
307 rowiterator endi=A_.end();
308 for (rowiterator i=A_.begin(); i!=endi; ++i)
309 {
310 // muck up extra rows
311 if ((int) i.index() >= (int) size())
312 {
313 coliterator endj=(*i).end();
314 for (coliterator j=(*i).begin(); j!=endj; ++j)
315 {
316 (*j) = 0;
317 if (j.index()==i.index())
318 (*j)[0][0] = 1;
319 }
320 f[i.index()] = 0;
321 continue;
322 }
323
324 // insert Dirichlet ans processor boundary conditions
325 if (essential[i.index()][0]!=BoundaryConditions::neumann)
326 {
327 coliterator endj=(*i).end();
328 for (coliterator j=(*i).begin(); j!=endj; ++j)
329 if (j.index()==i.index())
330 {
331 (*j)[0][0] = 1;
332 }
333 else
334 {
335 (*j)[0][0] = 0;
336 }
337 u[i.index()][0] = f[i.index()][0];
338 }
339 }
340}
341} // end namespace Dumux
342
343#endif
Definition of boundary condition types, extend if necessary.
defines intersection mappers.
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
@ dirichlet
Dirichlet boundary.
Definition: boundaryconditions.hh:41
@ neumann
Neumann boundary.
Definition: boundaryconditions.hh:39
bool isDirichlet(unsigned eqIdx) const
Returns true if an equation is used to specify a Dirichlet condition.
Definition: common/boundarytypes.hh:236
defines an intersection mapper for mapping of global DOFs assigned to faces which also works for adap...
Definition: intersectionmapper.hh:185
unsigned int size() const
Definition: intersectionmapper.hh:292
void update()
Definition: intersectionmapper.hh:307
Definition: matrix.hh:32
Extends CROperatorBase by a generic methods to assemble global stiffness matrix from local stiffness ...
Definition: croperatoradaptive.hh:69
void updateMatrix()
Definition: croperatoradaptive.hh:170
const IntersectionMapper & intersectionMapper()
Definition: croperatoradaptive.hh:130
CROperatorAssemblerTwoPAdaptive(const GridView &gridview)
Definition: croperatoradaptive.hh:100
void assemble(LocalStiffness &loc, Vector &u, Vector &f)
Assembles global stiffness matrix.
Definition: croperatoradaptive.hh:232
void adapt()
Definition: croperatoradaptive.hh:111
void initialize()
Definition: croperatoradaptive.hh:106
const IS & is_
Definition: croperatoradaptive.hh:164
const RepresentationType & operator*() const
Returns const reference to operator matrix.
Definition: croperatoradaptive.hh:119
const IntersectionMapper & intersectionMapper() const
Definition: croperatoradaptive.hh:135
const GridView gridView_
Definition: croperatoradaptive.hh:163
IntersectionMapper intersectionMapper_
Definition: croperatoradaptive.hh:165
RepresentationType A_
Definition: croperatoradaptive.hh:166
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.