dune-pdelab  2.4.1
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/debugstream.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>::Geometry Geometry;
55  typedef typename GF::Traits::RangeType Range;
56  typedef typename GF::Traits::DomainFieldType DF;
57  static const int dimD = GF::Traits::dimDomain;
58  typedef Dune::QuadratureRule<DF,dimD> QR;
59  typedef Dune::QuadratureRules<DF,dimD> QRs;
60 
61  sum = 0;
62  Range val;
63  for(const auto& cell : elements(gf.getGridView(),Dune::Partitions::interior)) {
64  const Geometry& geo = cell.geometry();
65  Dune::GeometryType gt = geo.type();
66  const QR& rule = QRs::rule(gt,qorder);
67  for (const auto& qip : rule) {
68  // evaluate the given grid functions at integration point
69  gf.evaluate(cell,qip.position(),val);
70 
71  // accumulate error
72  val *= qip.weight() * geo.integrationElement(qip.position());
73  sum += val;
74  }
75  }
76  }
77 
79 
88  template<typename GF>
90  typedef typename GF::Traits::GridViewType GV;
91  typedef typename GV::template Codim<0>::Entity Entity;
92  typedef typename GF::Traits::DomainType Domain;
93  typedef typename GF::Traits::RangeType Range;
94 
95  public:
97 
106  template<class GFHandle>
107  GridFunctionProbe(const GFHandle& gf, const Domain& xg)
108  {
109  setGridFunction(gf);
110  xl = 0;
111  evalRank = gfp->getGridView().comm().size();
112  int myRank = gfp->getGridView().comm().rank();
113  try {
114  e.reset(new Entity
115  (HierarchicSearch<typename GV::Grid, typename GV::IndexSet>
116  (gfp->getGridView().grid(), gfp->getGridView().indexSet()).
117  template findEntity<Interior_Partition>(xg)));
118  // make sure only interior entities are accepted
119  if(e->partitionType() == InteriorEntity)
120  evalRank = myRank;
121  }
122  catch(const Dune::GridError&) { /* do nothing */ }
123  evalRank = gfp->getGridView().comm().min(evalRank);
124  if(myRank == evalRank)
125  xl = e->geometry().local(xg);
126  else
127  e.reset();
128  if(myRank == 0 && evalRank == gfp->getGridView().comm().size())
129  dwarn << "Warning: GridFunctionProbe at (" << xg << ") is outside "
130  << "the grid" << std::endl;
131  }
132 
134 
139  void setGridFunction(const GF &gf) {
140  gfsp.reset();
141  gfp = &gf;
142  }
143 
145 
151  void setGridFunction(const GF *gf) {
152  gfsp.reset(gf);
153  gfp = gf;
154  }
155 
157 
165  void setGridFunction(const std::shared_ptr<const GF> &gf) {
166  gfsp = gf;
167  gfp = &*gf;
168  }
169 
171 
177  void eval_all(Range& val) const {
178  typedef typename GF::Traits::RangeFieldType RF;
179  if(evalRank == gfp->getGridView().comm().size())
180  val = std::numeric_limits<RF>::quiet_NaN();
181  else {
182  if(gfp->getGridView().comm().rank() == evalRank)
183  gfp->evaluate(*e, xl, val);
184  gfp->getGridView().comm().broadcast(&val,1,evalRank);
185  }
186  }
187 
189 
201  void eval(Range& val, int rank = 0) const {
202  eval_all(val);
203  }
204 
205  private:
206  std::shared_ptr<const GF> gfsp;
207  const GF *gfp;
208  std::shared_ptr<Entity> e;
209  Domain xl;
210  int evalRank;
211  };
212 
214 
215  } // namespace PDELab
216 } // namespace Dune
217 
218 #endif // DUNE_PDELAB_COMMON_FUNCTIONUTILITIES_HH
void setGridFunction(const GF *gf)
Set a new GridFunction.
Definition: functionutilities.hh:151
void eval(Range &val, int rank=0) const
evaluate the GridFunction and communicate result to the given rank
Definition: functionutilities.hh:201
GridFunctionProbe(const GFHandle &gf, const Domain &xg)
Constructor.
Definition: functionutilities.hh:107
Definition: adaptivity.hh:27
Evaluate a GridFunction at a certain global coordinate.
Definition: functionutilities.hh:89
void eval_all(Range &val) const
evaluate the GridFunction and broadcast result to all ranks
Definition: functionutilities.hh:177
void setGridFunction(const std::shared_ptr< const GF > &gf)
Set a new GridFunction.
Definition: functionutilities.hh:165
void integrateGridFunction(const GF &gf, typename GF::Traits::RangeType &sum, unsigned qorder=1)
Integrate a GridFunction.
Definition: functionutilities.hh:50
void setGridFunction(const GF &gf)
Set a new GridFunction.
Definition: functionutilities.hh:139
XG & xg
Definition: interpolate.hh:67