13#ifndef ERROR_DISTRIBUTION_HH
14#define ERROR_DISTRIBUTION_HH
18#include <boost/utility.hpp>
19#include <boost/fusion/include/as_vector.hpp>
37 using namespace AssemblyDetail;
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)
46 if( f.considerControlVariable_ )
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));
53 if( f.considerStateVariable_ )
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));
60 if( f.considerAdjointVariable_ )
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));
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)
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;
80 template <
class ArgU,
class ArgY,
class ArgP>
83 return u.value*u.value + y.value*y.value + p.value*p.value;
86 template <
class ArgU,
class ArgY,
class ArgP>
90 return sp(u.gradient,u.gradient) + sp(y.gradient,y.gradient) + sp(p.gradient,p.gradient);
94 template <
class Functional,
class ExtendedAnsatzVars >
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;
109 typedef typename Functional::Scalar
Scalar;
115 typedef typename AnsatzVars::template CoefficientVectorRepresentation<0,1>::type
ErrorVector;
119 static int const dim = Grid::dimension;
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;
140 template <
class Entity>
144 domainCache.moveTo(entity);
147 template <
class Position,
class Evaluators>
159 size_t nQuadPos = qr.size();
160 for (
size_t g=0; g<nQuadPos; ++g)
169 domainCache.evaluateAt(qr[g].position(),originalEvaluators);
171 evaluateData(originalEvaluators, extendedEvaluators, f, yl, yh, ul, uh, pl, ph, qr[g].weight());
177 template<
int row,
int dim>
188 u.gradient += ul.gradient;
190 y.gradient += yl.gradient;
192 p.gradient += pl.gradient;
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.;
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.;
214 result += Functional::template
D2<pIdx,pIdx>::present ? domainCache.template d2_impl<yIdx,yIdx>(p,p) : 0.;
216 return result*arg.
value[0];
219 template<
int row,
int col,
int dim>
224 typename Functional::DomainCache domainCache;
225 typename AnsatzVars::Grid::template Codim<0>::Entity
const* e;
243 template <
class FaceIterator>
247 boundaryCache.moveTo(entity);
250 template <
class Evaluators>
264 size_t nQuadPos = qr.size();
265 for (
size_t g=0; g<nQuadPos; ++g)
274 boundaryCache.evaluateAt(qr[g].position(),originalEvaluators);
276 evaluateData(originalEvaluators, extendedEvaluators, f, yl, yh, ul, uh, pl, ph, qr[g].weight());
282 template<
int row,
int dim>
293 u.gradient += ul.gradient;
295 y.gradient += yl.gradient;
297 p.gradient += pl.gradient;
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.;
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.;
316 result += Functional::template
D2<pIdx,pIdx>::present ? boundaryCache.template d2_impl<yIdx,yIdx>(p,p) : 0.;
318 return result*arg.
value[0];
321 template<
int row,
int col,
int dim>
326 typename Functional::BoundaryCache boundaryCache;
327 typename AnsatzVars::Grid::LeafGridView::IntersectionIterator
const* e;
338 ErrorDistribution(Functional
const& functional_,
typename OriginalAnsatzVars::VariableSet
const& iterate_,
339 typename OriginalAnsatzVars::VariableSet
const& errorEstimateL_,
typename ExtendedAnsatzVars::VariableSet
const& errorEstimateH_)
352 template <
int row,
int col>
360 template <
class Cell>
378 typename OriginalAnsatzVars::VariableSet
const&
iterate;
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)
void moveTo(Entity const &entity)
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)
OriginalAnsatzVars OriginVars
static ProblemType const type
bool considerAdjointVariable_
VariableSetDescription< AnsatzSpaces, VariableDescriptions > AnsatzVars
static constexpr int pSHIdx
AnsatzVars::template CoefficientVectorRepresentation< 0, 1 >::type ErrorVector
static constexpr int yIdx
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)
bool considerControlVariable_
bool considerStateVariable_
bool useStateNormForAdjoint_
void considerAdjointVariable(bool consider)
Functional const & functional
void considerControlVariable(bool consider)
Functional::Scalar Scalar
AnsatzSpaces ansatzSpaces
OriginalAnsatzVars::VariableSet const & iterate
FEFunctionSpace< DiscontinuousLagrangeMapper< Scalar, typename Grid::LeafGridView > > AnsatzSpace
Variable< SpaceIndex< 0 >, Components< 1 >, VariableId< 0 > > AnsatzVariableInformation
static constexpr int uIdx
static constexpr int uSLIdx
int integrationOrder(Cell const &, int, bool) const
static constexpr int ySLIdx
void ignoreLowerOrderError(bool ignore)
boost::fusion::vector< AnsatzVariableInformation > VariableDescriptions
static constexpr int pIdx
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.
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.
static bool const constant
static bool const present
static bool const present
static bool const symmetric
Euclidean scalar product.
A cache that stores suitable quadrature rules for quick retrieval.
static Rule const & rule(const Dune::GeometryType &t, int p)
A class defining elementary information about a single variable.
A class that stores values, gradients, and Hessians of evaluated shape functions.
ValueType value
The shape function's value, a vector of dimension components