dune-pdelab  2.7-git
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 #include <memory>
10 
11 #include <dune/common/shared_ptr.hh>
12 #include <dune/common/debugstream.hh>
13 
14 #include <dune/geometry/quadraturerules.hh>
15 #include <dune/geometry/type.hh>
16 
17 #include <dune/grid/common/gridenums.hh>
18 #include <dune/grid/utility/hierarchicsearch.hh>
19 
20 namespace Dune {
21  namespace PDELab {
22 
26 
27 
29 
50  template<typename GF>
51  typename GF::Traits::RangeType integrateGridFunction(const GF& gf,
52  unsigned qorder = 1) {
53  typename GF::Traits::RangeType sum;
54  integrateGridFunction(gf, sum, qorder);
55  return sum;
56  }
57 
59 
81  template<typename GF>
82  void integrateGridFunction(const GF& gf,
83  typename GF::Traits::RangeType& sum,
84  unsigned qorder = 1) {
85  typedef typename GF::Traits::GridViewType GV;
86  typedef typename GV::template Codim<0>::Geometry Geometry;
87  typedef typename GF::Traits::RangeType Range;
88  typedef typename GF::Traits::DomainFieldType DF;
89  static const int dimD = GF::Traits::dimDomain;
90  typedef Dune::QuadratureRule<DF,dimD> QR;
91  typedef Dune::QuadratureRules<DF,dimD> QRs;
92 
93  sum = 0;
94  Range val;
95  for(const auto& cell : elements(gf.getGridView(),Dune::Partitions::interior)) {
96  const Geometry& geo = cell.geometry();
97  Dune::GeometryType gt = geo.type();
98  const QR& rule = QRs::rule(gt,qorder);
99  for (const auto& qip : rule) {
100  // evaluate the given grid functions at integration point
101  gf.evaluate(cell,qip.position(),val);
102 
103  // accumulate error
104  val *= qip.weight() * geo.integrationElement(qip.position());
105  sum += val;
106  }
107  }
108  }
109 
111 
120  template<typename GF>
122  typedef typename GF::Traits::GridViewType GV;
123  typedef typename GV::template Codim<0>::Entity Entity;
124  typedef typename GF::Traits::DomainType Domain;
125  typedef typename GF::Traits::RangeType Range;
126 
127  public:
129 
138  GridFunctionProbe(const GF& gf, const Domain& xg)
139  : gfp(Dune::stackobject_to_shared_ptr(gf))
140  {
141  hierarchic_search(gfp->getGridView(), xg);
142  }
143 
144  GridFunctionProbe(std::shared_ptr<const GF> gf, const Domain& xg)
145  : gfp(gf)
146  {
147  hierarchic_search(gfp->getGridView(), xg);
148  }
149 
157  GridFunctionProbe(const GV& gv, const Domain& xg)
158  {
159  hierarchic_search(gv, xg);
160  }
161 
163 
168  void setGridFunction(const GF &gf) {
169  gfp = Dune::stackobject_to_shared_ptr(gf);
170  }
171 
173 
179  void setGridFunction(const GF *gf) {
180  gfp.reset(gf);
181  }
182 
184 
192  void setGridFunction(const std::shared_ptr<const GF> &gf) {
193  gfp = gf;
194  }
195 
197 
203  void eval_all(Range& val) const {
204  typedef typename GF::Traits::RangeFieldType RF;
205  if(evalRank == gfp->getGridView().comm().size())
206  val = std::numeric_limits<RF>::quiet_NaN();
207  else {
208  if(gfp->getGridView().comm().rank() == evalRank)
209  gfp->evaluate(*e, xl, val);
210  gfp->getGridView().comm().broadcast(&val,1,evalRank);
211  }
212  }
213 
215 
227  void eval(Range& val, int rank = 0) const {
228  eval_all(val);
229  }
230 
232 
235  int getEvalRank() {
236  return evalRank;
237  }
238 
239  private:
240  void hierarchic_search(const GV& gv, const Domain& xg)
241  {
242  xl = 0;
243  evalRank = gv.comm().size();
244  int myRank = gv.comm().rank();
245  try {
246  e.reset(new Entity
247  (HierarchicSearch<typename GV::Grid, typename GV::IndexSet>
248  (gv.grid(), gv.indexSet()).
249  template findEntity<Interior_Partition>(xg)));
250  // make sure only interior entities are accepted
251  if(e->partitionType() == InteriorEntity)
252  evalRank = myRank;
253  }
254  catch(const Dune::GridError&) { /* do nothing */ }
255  evalRank = gv.comm().min(evalRank);
256  if(myRank == evalRank)
257  xl = e->geometry().local(xg);
258  else
259  e.reset();
260  if(myRank == 0 && evalRank == gv.comm().size())
261  dwarn << "Warning: GridFunctionProbe at (" << xg << ") is outside "
262  << "the grid" << std::endl;
263  }
264 
265  std::shared_ptr<const GF> gfp;
266  std::shared_ptr<Entity> e;
267  Domain xl;
268  int evalRank;
269  };
270 
272 
273  } // namespace PDELab
274 } // namespace Dune
275 
276 #endif // DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
XG & xg
Definition: interpolate.hh:55
GF::Traits::RangeType integrateGridFunction(const GF &gf, unsigned qorder=1)
Integrate a GridFunction.
Definition: functionutilities.hh:51
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Evaluate a GridFunction at a certain global coordinate.
Definition: functionutilities.hh:121
GridFunctionProbe(const GV &gv, const Domain &xg)
Definition: functionutilities.hh:157
void setGridFunction(const std::shared_ptr< const GF > &gf)
Set a new GridFunction.
Definition: functionutilities.hh:192
void eval(Range &val, int rank=0) const
evaluate the GridFunction and communicate result to the given rank
Definition: functionutilities.hh:227
GridFunctionProbe(std::shared_ptr< const GF > gf, const Domain &xg)
Definition: functionutilities.hh:144
void setGridFunction(const GF &gf)
Set a new GridFunction.
Definition: functionutilities.hh:168
void setGridFunction(const GF *gf)
Set a new GridFunction.
Definition: functionutilities.hh:179
GridFunctionProbe(const GF &gf, const Domain &xg)
Constructor.
Definition: functionutilities.hh:138
void eval_all(Range &val) const
evaluate the GridFunction and broadcast result to all ranks
Definition: functionutilities.hh:203
int getEvalRank()
MPI rank of coordinate.
Definition: functionutilities.hh:235