dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
20namespace 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
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