KASKADE 7 development version
errorDistribution.hh
Go to the documentation of this file.
1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2/* */
3/* This file is part of the library KASKADE 7 */
4/* see http://www.zib.de/Numerik/numsoft/kaskade7/ */
5/* */
6/* Copyright (C) 2002-2015 Zuse Institute Berlin */
7/* */
8/* KASKADE 7 is distributed under the terms of the ZIB Academic License. */
9/* see $KASKADE/academic.txt */
10/* */
11/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12
13#ifndef ERROR_DISTRIBUTION_HH
14#define ERROR_DISTRIBUTION_HH
15
16#include <string>
17
18#include <boost/utility.hpp>
19#include <boost/fusion/include/as_vector.hpp>
20
21#include "fem/assemble.hh"
22#include "fem/functionspace.hh"
23#include "fem/lagrangespace.hh"
25
26// forward declarations
27namespace Dune
28{
29 template <class,int> class FieldVector;
30 template <class,int,int> class FieldMatrix;
31 template <class, int> class QuadratureRule;
32}
33
34namespace Kaskade
35{
36 using namespace boost::fusion;
37 using namespace AssemblyDetail;
38
39 enum class ErrorNorm { Energy, L2, H1, H1_half };
40
41 template <class OriginalEvaluators, class ExtensionEvaluators, class Functional,
42 class ArgYH, class ArgYE, class ArgUH, class ArgUE, class ArgPH, class ArgPE>
43 void evaluateData(OriginalEvaluators& originalEvaluators, ExtensionEvaluators& extendedEvaluators, Functional const& f,
44 ArgYH& yl, ArgYE& yh, ArgUH& ul, ArgUE& uh, ArgPH& pl, ArgPE& ph, double w)
45 {
46 if( f.considerControlVariable_ )
47 {
48 uh.value += w * at_c<Functional::uIdx>(f.errorEstimateH.data).value(at_c<Functional::uSHIdx>(extendedEvaluators));
49 ul.value += w * at_c<Functional::uIdx>(f.errorEstimateL.data).value(at_c<Functional::uSLIdx>(originalEvaluators));
50 uh.gradient += w * at_c<Functional::uIdx>(f.errorEstimateH.data).gradient(at_c<Functional::uSHIdx>(extendedEvaluators));
51 ul.gradient += w * at_c<Functional::uIdx>(f.errorEstimateL.data).gradient(at_c<Functional::uSLIdx>(originalEvaluators));
52 }
53 if( f.considerStateVariable_ )
54 {
55 yh.value += w * at_c<Functional::yIdx>(f.errorEstimateH.data).value(at_c<Functional::ySHIdx>(extendedEvaluators));
56 yl.value += w * at_c<Functional::yIdx>(f.errorEstimateL.data).value(at_c<Functional::ySLIdx>(originalEvaluators));
57 yh.gradient += w * at_c<Functional::yIdx>(f.errorEstimateH.data).gradient(at_c<Functional::ySHIdx>(extendedEvaluators));
58 yl.gradient += w * at_c<Functional::yIdx>(f.errorEstimateL.data).gradient(at_c<Functional::ySLIdx>(originalEvaluators));
59 }
60 if( f.considerAdjointVariable_ )
61 {
62 ph.value += w * at_c<Functional::pIdx>(f.errorEstimateH.data).value(at_c<Functional::pSHIdx>(extendedEvaluators));
63 pl.value += w * at_c<Functional::pIdx>(f.errorEstimateL.data).value(at_c<Functional::pSLIdx>(originalEvaluators));
64 ph.gradient += w * at_c<Functional::pIdx>(f.errorEstimateH.data).gradient(at_c<Functional::pSHIdx>(extendedEvaluators));
65 pl.gradient += w * at_c<Functional::pIdx>(f.errorEstimateL.data).gradient(at_c<Functional::pSLIdx>(extendedEvaluators));
66 }
67 }
68
69 template <class ArgYL, class ArgYH, class ArgUL, class ArgUH, class ArgPL, class ArgPH>
70 void clearVarArgs(ArgYL& yl, ArgYH& yh, ArgUL& ul, ArgUH& uh, ArgPL& pl, ArgPH& ph)
71 {
72 ul.value = 0; ul.gradient = 0;
73 yl.value = 0; yl.gradient = 0;
74 pl.value = 0; pl.gradient = 0;
75 uh.value = 0; uh.gradient = 0;
76 yh.value = 0; yh.gradient = 0;
77 ph.value = 0; ph.gradient = 0;
78 }
79
80 template <class ArgU, class ArgY, class ArgP>
81 double computeL2Error(ArgU const& u, ArgY const& y, ArgP const& p)
82 {
83 return u.value*u.value + y.value*y.value + p.value*p.value;
84 }
85
86 template <class ArgU, class ArgY, class ArgP>
87 double computeH1HalfError(ArgU const& u, ArgY const& y, ArgP const& p)
88 {
90 return sp(u.gradient,u.gradient) + sp(y.gradient,y.gradient) + sp(p.gradient,p.gradient);
91 }
92
93
94 template <class Functional, class ExtendedAnsatzVars >
96 {
97 typedef typename Functional::AnsatzVars OriginalAnsatzVars;
98 typedef typename OriginalAnsatzVars::Grid Grid;
101 typedef typename OriginalAnsatzVars::Spaces OriginalSpaces;
102 typedef typename ExtendedAnsatzVars::Spaces ExtendedSpaces;
103 typedef typename result_of::as_vector<typename result_of::transform<OriginalSpaces, GetEvaluators<SfCache> >::type>::type OriginalEvaluators;
104 typedef typename result_of::as_vector<typename result_of::transform<ExtendedSpaces, GetEvaluators<SfCache2> >::type>::type ExtendedEvaluators;
105 typedef typename Grid::ctype CoordType;
107 typedef Dune::QuadratureRule<typename Functional::AnsatzVars::Grid::ctype, Functional::AnsatzVars::Grid::dimension-1> FaceQuadRule;
108 public:
109 typedef typename Functional::Scalar Scalar;
111 typedef boost::fusion::vector<AnsatzSpace const*> AnsatzSpaces;
113 typedef boost::fusion::vector<AnsatzVariableInformation> VariableDescriptions;
115 typedef typename AnsatzVars::template CoefficientVectorRepresentation<0,1>::type ErrorVector;
117 typedef OriginalAnsatzVars OriginVars;
118 //typedef AnsatzVars OriginVars;
119 static int const dim = Grid::dimension;
120 static ProblemType const type = Functional::type;
121
122 static constexpr int yIdx = getStateId<Functional>();
123 static constexpr int uIdx = getControlId<Functional>();
124 static constexpr int pIdx = getAdjointId<Functional>();
125 static constexpr int uSLIdx = result_of::value_at_c<typename OriginalAnsatzVars::Variables, uIdx>::type::spaceIndex;
126 static constexpr int ySLIdx = result_of::value_at_c<typename OriginalAnsatzVars::Variables, yIdx>::type::spaceIndex;
127 static constexpr int pSLIdx = result_of::value_at_c<typename OriginalAnsatzVars::Variables, pIdx>::type::spaceIndex;
128 static constexpr int uSHIdx = result_of::value_at_c<typename ExtendedAnsatzVars::Variables, uIdx>::type::spaceIndex;
129 static constexpr int ySHIdx = result_of::value_at_c<typename ExtendedAnsatzVars::Variables, yIdx>::type::spaceIndex;
130 static constexpr int pSHIdx = result_of::value_at_c<typename ExtendedAnsatzVars::Variables, pIdx>::type::spaceIndex;
131
133 {
134 public:
135 DomainCache(ErrorDistribution const& f_, typename OriginVars::VariableSet const& /*vars*/, int flags=7)
136 : f(f_), domainCache(f.functional,f.iterate,flags)
137 {
138 }
139
140 template <class Entity>
141 void moveTo(Entity const &entity)
142 {
143 e = &entity;
144 domainCache.moveTo(entity);
145 }
146
147 template <class Position, class Evaluators>
148 void evaluateAt(Position const&, Evaluators const& evaluators)
149 {
150 clearVarArgs(yl,yh,ul,uh,pl,ph);
151 OriginalEvaluators originalEvaluators(transform(f.iterate.descriptions.spaces,GetEvaluators<SfCache>(&sfCache)));
152 ExtendedEvaluators extendedEvaluators(transform(f.errorEstimateH.descriptions.spaces,GetEvaluators<SfCache>(&extendedSFCache)));
154 moveEvaluatorsToCell(originalEvaluators,*e);
155 moveEvaluatorsToCell(extendedEvaluators,*e);
156 useQuadratureRuleInEvaluators(originalEvaluators,qr,0);
157 useQuadratureRuleInEvaluators(extendedEvaluators,qr,0);
158
159 size_t nQuadPos = qr.size();
160 for (size_t g=0; g<nQuadPos; ++g)
161 {
162 // pos of integration point
163 Dune::FieldVector<CoordType,dim> quadPos = qr[g].position();
164 // for all spaces involved, update the evaluators associated
165 // to this quadrature point
166 moveEvaluatorsToIntegrationPoint(originalEvaluators,quadPos);
167 moveEvaluatorsToIntegrationPoint(extendedEvaluators,quadPos);
168 // prepare evaluation of functional
169 domainCache.evaluateAt(qr[g].position(),originalEvaluators);
170
171 evaluateData(originalEvaluators, extendedEvaluators, f, yl, yh, ul, uh, pl, ph, qr[g].weight());
172 }
173 }
174
175 Scalar d0() const { return 0; }
176
177 template<int row, int dim>
180 {
184
185 if( !f.onlyH_ )
186 {
187 u.value += ul.value;
188 u.gradient += ul.gradient;
189 y.value += yl.value;
190 y.gradient += yl.gradient;
191 p.value += pl.value;
192 p.gradient += pl.gradient;
193 }
194
195 if( f.errorNorm == ErrorNorm::L2 ) return computeL2Error(u,y,p);
196 if( f.errorNorm == ErrorNorm::H1_half ) return computeH1HalfError(u,y,p);
197 if( f.errorNorm == ErrorNorm::H1 ) return computeL2Error(u,y,p) + computeH1HalfError(u,y,p);
198
199 Scalar result = 0;
200 result += Functional::template D2<yIdx,yIdx>::present ? domainCache.template d2_impl<yIdx,yIdx>(y,y) : 0.;
201 result += Functional::template D2<yIdx,uIdx>::present ? domainCache.template d2_impl<yIdx,uIdx>(y,u) : 0.;
202 result += Functional::template D2<uIdx,yIdx>::present ? domainCache.template d2_impl<uIdx,yIdx>(u,y) : 0.;
203 result += Functional::template D2<uIdx,uIdx>::present ? domainCache.template d2_impl<uIdx,uIdx>(u,u) : 0.;
205 {
206 result += Functional::template D2<yIdx,pIdx>::present ? domainCache.template d2_impl<yIdx,pIdx>(y,p) : 0.;
207 result += Functional::template D2<pIdx,yIdx>::present ? domainCache.template d2_impl<pIdx,yIdx>(p,y) : 0.;
208 result += Functional::template D2<pIdx,uIdx>::present ? domainCache.template d2_impl<pIdx,uIdx>(p,u) : 0.;
209 result += Functional::template D2<uIdx,pIdx>::present ? domainCache.template d2_impl<uIdx,pIdx>(u,p) : 0.;
210 result += Functional::template D2<pIdx,pIdx>::present ? domainCache.template d2_impl<pIdx,pIdx>(p,p) : 0.;
211 }
212 else
213 {
214 result += Functional::template D2<pIdx,pIdx>::present ? domainCache.template d2_impl<yIdx,yIdx>(p,p) : 0.;
215 }
216 return result*arg.value[0];
217 }
218
219 template<int row, int col, int dim>
221
222 private:
223 ErrorDistribution const& f;
224 typename Functional::DomainCache domainCache;
225 typename AnsatzVars::Grid::template Codim<0>::Entity const* e;
232 SfCache sfCache;
233 SfCache2 extendedSFCache;
234 };
235
237 {
238 public:
239 BoundaryCache(ErrorDistribution const& f_, typename OriginVars::VariableSet const& vars, int flags=7)
240 : f(f_), boundaryCache(f.functional,f.iterate,flags)
241 {}
242
243 template <class FaceIterator>
244 void moveTo(FaceIterator const& entity)
245 {
246 e = &entity;
247 boundaryCache.moveTo(entity);
248 }
249
250 template <class Evaluators>
252 {
253 clearVarArgs(yl,yh,ul,uh,pl,ph);
254 OriginalEvaluators originalEvaluators(transform(f.iterate.descriptions.spaces,GetEvaluators<SfCache>(&sfCache)));
255 ExtendedEvaluators extendedEvaluators(transform(f.errorEstimateH.descriptions.spaces,GetEvaluators<SfCache>(&extendedSFCache)));
256 FaceQuadRule qr = QuadratureTraits<FaceQuadRule>().rule((*e)->geometryInInside().type(),f.qOrder);
257 moveEvaluatorsToCell(originalEvaluators,*at_c<ySLIdx>(evaluators).cell_);
258 moveEvaluatorsToCell(extendedEvaluators,*at_c<ySLIdx>(evaluators).cell_);
259 moveEvaluatorsToCell(originalEvaluators,*e);
260 moveEvaluatorsToCell(extendedEvaluators,*e);
261 useQuadratureRuleInEvaluators(originalEvaluators,qr,(*e)->indexInInside());
262 useQuadratureRuleInEvaluators(extendedEvaluators,qr,(*e)->indexInInside());
263
264 size_t nQuadPos = qr.size();
265 for (size_t g=0; g<nQuadPos; ++g)
266 {
267 // pos of integration point
268 Dune::FieldVector<CoordType,dim> quadPos = (*e)->geometryInInside().global(qr[g].position());
269 // for all spaces involved, update the evaluators associated
270 // to this quadrature point
271 moveEvaluatorsToIntegrationPoint(originalEvaluators,quadPos);
272 moveEvaluatorsToIntegrationPoint(extendedEvaluators,quadPos);
273 // prepare evaluation of functional
274 boundaryCache.evaluateAt(qr[g].position(),originalEvaluators);
275
276 evaluateData(originalEvaluators, extendedEvaluators, f, yl, yh, ul, uh, pl, ph, qr[g].weight());
277 }
278 }
279
280 Scalar d0() const { return 0; }
281
282 template<int row, int dim>
284 {
289
290 if( !f.onlyH_ )
291 {
292 u.value += ul.value;
293 u.gradient += ul.gradient;
294 y.value += yl.value;
295 y.gradient += yl.gradient;
296 p.value += pl.value;
297 p.gradient += pl.gradient;
298 }
299
300 if( f.errorNorm != ErrorNorm::Energy ) return computeL2Error(u,y,p);
301
302 Scalar result = Functional::template D2<yIdx,yIdx>::present ? boundaryCache.template d2_impl<yIdx,yIdx>(y,y) : 0.;
303 result += Functional::template D2<yIdx,uIdx>::present ? boundaryCache.template d2_impl<yIdx,uIdx>(y,u) : 0.;
304 result += Functional::template D2<uIdx,yIdx>::present ? boundaryCache.template d2_impl<uIdx,yIdx>(u,y) : 0.;
305 result += Functional::template D2<uIdx,uIdx>::present ? boundaryCache.template d2_impl<uIdx,uIdx>(u,u) : 0.;
307 {
308 result += Functional::template D2<yIdx,pIdx>::present ? boundaryCache.template d2_impl<yIdx,pIdx>(y,p) : 0.;
309 result += Functional::template D2<pIdx,yIdx>::present ? boundaryCache.template d2_impl<pIdx,yIdx>(p,y) : 0.;
310 result += Functional::template D2<pIdx,uIdx>::present ? boundaryCache.template d2_impl<pIdx,uIdx>(p,u) : 0.;
311 result += Functional::template D2<uIdx,pIdx>::present ? boundaryCache.template d2_impl<uIdx,pIdx>(u,p) : 0.;
312 result += Functional::template D2<pIdx,pIdx>::present ? boundaryCache.template d2_impl<pIdx,pIdx>(p,p) : 0.;
313 }
314 else
315 {
316 result += Functional::template D2<pIdx,pIdx>::present ? boundaryCache.template d2_impl<yIdx,yIdx>(p,p) : 0.;
317 }
318 return result*arg.value[0];
319 }
320
321 template<int row, int col, int dim>
323
324 private:
325 ErrorDistribution const& f;
326 typename Functional::BoundaryCache boundaryCache;
327 typename AnsatzVars::Grid::LeafGridView::IntersectionIterator const* e;
334 SfCache sfCache;
335 SfCache2 extendedSFCache;
336 };
337
338 ErrorDistribution(Functional const& functional_, typename OriginalAnsatzVars::VariableSet const& iterate_,
339 typename OriginalAnsatzVars::VariableSet const& errorEstimateL_, typename ExtendedAnsatzVars::VariableSet const& errorEstimateH_)
340 : functional(functional_), iterate(iterate_), errorEstimateL(errorEstimateL_), errorEstimateH(errorEstimateH_), ansatzSpace(boost::fusion::at_c<0>(iterate.descriptions.spaces)->gridManager(), iterate.descriptions.gridView,0),
342
343 {}
344
345 template <int row>
346 struct D1
347 {
348 static bool const present = true;
349 static bool const constant = false;
350 };
351
352 template <int row, int col>
353 struct D2
354 {
355 static bool const present = false;
356 static bool const symmetric = false;
357 static bool const lumped = false;
358 };
359
360 template <class Cell>
361 int integrationOrder(Cell const&, int, bool) const { return 0; }
362
363 AnsatzSpaces const& getSpaces() const { return ansatzSpaces; }
364
366
367 void ignoreLowerOrderError(bool ignore) { onlyH_ = ignore; }
368 void considerStateVariable(bool consider) { considerStateVariable_ = consider; }
369 void considerControlVariable(bool consider) { considerControlVariable_ = consider; }
370 void considerAdjointVariable(bool consider) { considerAdjointVariable_ = consider; }
371 void setErrorNorm(ErrorNorm errNorm) { errorNorm = errNorm; }
372 void useStateNormForAdjoint(bool useStateNorm) { useStateNormForAdjoint_ = useStateNorm; }
373
374 friend class DomainCache;
375 friend class BoundaryCache;
376
377 Functional const& functional;
378 typename OriginalAnsatzVars::VariableSet const& iterate;
379 typename OriginalAnsatzVars::VariableSet const& errorEstimateL;
380 typename ExtendedAnsatzVars::VariableSet const& errorEstimateH;
383 std::string varName[1];
387 int qOrder = 6;
388 };
389}
390#endif
Dune::FieldMatrix< Scalar, 1, 1 > d2(VariationalArg< Scalar, dim > const &, VariationalArg< Scalar, dim > const &) const
void moveTo(FaceIterator const &entity)
void evaluateAt(Dune::FieldVector< typename AnsatzVars::Grid::ctype, AnsatzVars::Grid::dimension-1 > const &xlocal, Evaluators const &evaluators)
BoundaryCache(ErrorDistribution const &f_, typename OriginVars::VariableSet const &vars, int flags=7)
Dune::FieldVector< Scalar, 1 > d1(VariationalArg< Scalar, dim > const &arg) const
DomainCache(ErrorDistribution const &f_, typename OriginVars::VariableSet const &, int flags=7)
Dune::FieldMatrix< Scalar, 1, 1 > d2(VariationalArg< Scalar, dim > const &, VariationalArg< Scalar, dim > const &) const
void evaluateAt(Position const &, Evaluators const &evaluators)
Dune::FieldVector< Scalar, TestVars::template Components< row >::m > d1(VariationalArg< Scalar, dim > const &arg) const
static constexpr int ySHIdx
static constexpr int uSHIdx
void setErrorNorm(ErrorNorm errNorm)
static ProblemType const type
VariableSetDescription< AnsatzSpaces, VariableDescriptions > AnsatzVars
static constexpr int pSHIdx
AnsatzVars::template CoefficientVectorRepresentation< 0, 1 >::type ErrorVector
ExtendedAnsatzVars::VariableSet const & errorEstimateH
AnsatzSpaces const & getSpaces() const
boost::fusion::vector< AnsatzSpace const * > AnsatzSpaces
AnsatzVars const & getVariableSetDescription() const
OriginalAnsatzVars::VariableSet const & errorEstimateL
ErrorDistribution(Functional const &functional_, typename OriginalAnsatzVars::VariableSet const &iterate_, typename OriginalAnsatzVars::VariableSet const &errorEstimateL_, typename ExtendedAnsatzVars::VariableSet const &errorEstimateH_)
void considerStateVariable(bool consider)
static constexpr int pSLIdx
void useStateNormForAdjoint(bool useStateNorm)
void considerAdjointVariable(bool consider)
void considerControlVariable(bool consider)
OriginalAnsatzVars::VariableSet const & iterate
FEFunctionSpace< DiscontinuousLagrangeMapper< Scalar, typename Grid::LeafGridView > > AnsatzSpace
Variable< SpaceIndex< 0 >, Components< 1 >, VariableId< 0 > > AnsatzVariableInformation
static constexpr int uSLIdx
int integrationOrder(Cell const &, int, bool) const
static constexpr int ySLIdx
void ignoreLowerOrderError(bool ignore)
boost::fusion::vector< AnsatzVariableInformation > VariableDescriptions
A representation of a finite element function space defined over a domain covered by a grid.
A boost::fusion functor for creating an evaluator from the space, using a shape function cache that i...
FEFunctionSpace and FunctionSpaceElement and the like.
void useQuadratureRuleInEvaluators(Evaluator &evals, QuadratureRule const &qr, int subId)
Tells all evaluators to use the given quadrature rule on the given subentity.
void moveEvaluatorsToIntegrationPoint(Evaluators &evals, Dune::FieldVector< CoordType, dim > const &x, Dune::QuadratureRule< CoordType, subDim > const &qr, int ip, int subId)
Moves all provided evaluators to the given integration point, evaluating the shape functions there.
typename boost::fusion::result_of::as_vector< typename boost::fusion::result_of::transform< Spaces, GetEvaluatorTypes >::type >::type Evaluators
the heterogeneous sequence type of Evaluators for the given spaces
ProblemType
A type for describing the nature of a PDE problem.
typename GridView::template Codim< 0 >::Entity Cell
The type of cells (entities of codimension 0) in the grid view.
Definition: gridBasics.hh:35
Lagrange Finite Elements.
double computeH1HalfError(ArgU const &u, ArgY const &y, ArgP const &p)
void clearVarArgs(ArgYL &yl, ArgYH &yh, ArgUL &ul, ArgUH &uh, ArgPL &pl, ArgPH &ph)
void evaluateData(OriginalEvaluators &originalEvaluators, ExtensionEvaluators &extendedEvaluators, Functional const &f, ArgYH &yl, ArgYE &yh, ArgUH &ul, ArgUE &uh, ArgPH &pl, ArgPE &ph, double w)
void moveEvaluatorsToCell(Evaluators &evals, Cell const &cell)
Moves all provided evaluators to the given cell.
double computeL2Error(ArgU const &u, ArgY const &y, ArgP const &p)
Helper class for specifying the number of components of a variable.
Definition: variables.hh:100
A cache that stores suitable quadrature rules for quick retrieval.
Definition: quadrature.hh:35
static Rule const & rule(const Dune::GeometryType &t, int p)
Definition: quadrature.hh:38
A class defining elementary information about a single variable.
Definition: variables.hh:128
A class that stores values, gradients, and Hessians of evaluated shape functions.
ValueType value
The shape function's value, a vector of dimension components