dune-pdelab  2.0.0
functionutilities.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
5 #define DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
6 
7 #include <limits>
8 #include <ostream>
9 
10 #include <dune/common/debugstream.hh>
11 #include <dune/common/shared_ptr.hh>
12 
13 #include <dune/geometry/quadraturerules.hh>
14 #include <dune/geometry/type.hh>
15 
16 #include <dune/grid/common/gridenums.hh>
17 #include <dune/grid/utility/hierarchicsearch.hh>
18 
19 namespace Dune {
20  namespace PDELab {
21 
25 
27 
49  template<typename GF>
50  void integrateGridFunction(const GF& gf,
51  typename GF::Traits::RangeType& sum,
52  unsigned qorder = 1) {
53  typedef typename GF::Traits::GridViewType GV;
54  typedef typename GV::template Codim<0>::
55  template Partition<Interior_Partition>::Iterator EIterator;
56  typedef typename GV::template Codim<0>::Geometry Geometry;
57  typedef typename GF::Traits::RangeType Range;
58  typedef typename GF::Traits::DomainFieldType DF;
59  static const int dimD = GF::Traits::dimDomain;
60  typedef Dune::QuadratureRule<DF,dimD> QR;
61  typedef Dune::QuadratureRules<DF,dimD> QRs;
62  typedef typename QR::const_iterator QIterator;
63 
64  sum = 0;
65  Range val;
66  const EIterator eend = gf.getGridView().template end<0,
67  Interior_Partition>();
68  for(EIterator eit = gf.getGridView().template begin<0,
69  Interior_Partition>(); eit != eend; ++eit) {
70  const Geometry& geo = eit->geometry();
71  Dune::GeometryType gt = geo.type();
72  const QR& rule = QRs::rule(gt,qorder);
73  const QIterator qend = rule.end();
74 
75  for (QIterator qit=rule.begin(); qit != qend; ++qit)
76  {
77  // evaluate the given grid functions at integration point
78  gf.evaluate(*eit,qit->position(),val);
79 
80  // accumulate error
81  val *= qit->weight() * geo.integrationElement(qit->position());
82  sum += val;
83  }
84  }
85  }
86 
88 
97  template<typename GF>
99  typedef typename GF::Traits::GridViewType GV;
100  typedef typename GV::template Codim<0>::EntityPointer EPtr;
101  typedef typename GF::Traits::DomainType Domain;
102  typedef typename GF::Traits::RangeType Range;
103 
104  public:
106 
115  template<class GFHandle>
116  GridFunctionProbe(const GFHandle& gf, const Domain& xg)
117  {
118  setGridFunction(gf);
119  xl = 0;
120  evalRank = gfp->getGridView().comm().size();
121  int myRank = gfp->getGridView().comm().rank();
122  try {
123  e.reset(new EPtr
124  (HierarchicSearch<typename GV::Grid, typename GV::IndexSet>
125  (gfp->getGridView().grid(), gfp->getGridView().indexSet()).
126  template findEntity<Interior_Partition>(xg)));
127  // make sure only interior entities are accepted
128  if((*e)->partitionType() == InteriorEntity)
129  evalRank = myRank;
130  }
131  catch(const Dune::GridError&) { /* do nothing */ }
132  evalRank = gfp->getGridView().comm().min(evalRank);
133  if(myRank == evalRank)
134  xl = (*e)->geometry().local(xg);
135  else
136  e.reset();
137  if(myRank == 0 && evalRank == gfp->getGridView().comm().size())
138  dwarn << "Warning: GridFunctionProbe at (" << xg << ") is outside "
139  << "the grid" << std::endl;
140  }
141 
143 
148  void setGridFunction(const GF &gf) {
149  gfsp.reset();
150  gfp = &gf;
151  }
152 
154 
160  void setGridFunction(const GF *gf) {
161  gfsp.reset(gf);
162  gfp = gf;
163  }
164 
166 
174  void setGridFunction(const Dune::shared_ptr<const GF> &gf) {
175  gfsp = gf;
176  gfp = &*gf;
177  }
178 
180 
186  void eval_all(Range& val) const {
187  typedef typename GF::Traits::RangeFieldType RF;
188  if(evalRank == gfp->getGridView().comm().size())
189  val = std::numeric_limits<RF>::quiet_NaN();
190  else {
191  if(gfp->getGridView().comm().rank() == evalRank)
192  gfp->evaluate(**e, xl, val);
193  gfp->getGridView().comm().broadcast(&val,1,evalRank);
194  }
195  }
196 
198 
210  void eval(Range& val, int rank = 0) const {
211  eval_all(val);
212  }
213 
214  private:
215  Dune::shared_ptr<const GF> gfsp;
216  const GF *gfp;
217  Dune::shared_ptr<EPtr> e;
218  Domain xl;
219  int evalRank;
220  };
221 
223 
224  } // namespace PDELab
225 } // namespace Dune
226 
227 #endif // DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
void eval(Range &val, int rank=0) const
evaluate the GridFunction and communicate result to the given rank
Definition: functionutilities.hh:210
void eval_all(Range &val) const
evaluate the GridFunction and broadcast result to all ranks
Definition: functionutilities.hh:186
GridFunctionProbe(const GFHandle &gf, const Domain &xg)
Constructor.
Definition: functionutilities.hh:116
void integrateGridFunction(const GF &gf, typename GF::Traits::RangeType &sum, unsigned qorder=1)
Integrate a GridFunction.
Definition: functionutilities.hh:50
XG & xg
Definition: interpolate.hh:67
void setGridFunction(const Dune::shared_ptr< const GF > &gf)
Set a new GridFunction.
Definition: functionutilities.hh:174
void setGridFunction(const GF *gf)
Set a new GridFunction.
Definition: functionutilities.hh:160
void setGridFunction(const GF &gf)
Set a new GridFunction.
Definition: functionutilities.hh:148
Evaluate a GridFunction at a certain global coordinate.
Definition: functionutilities.hh:98