integral.hh

Go to the documentation of this file.
1 
7 #ifndef spcIntegral_hh
8 #define spcIntegral_hh
9 
10 #include <set>
11 #include "basics/typedefs.hh"
12 #include "basics/vectorsMatrices.hh"
13 #include "basics/vectorsMatricesForward.hh"
14 #include "basics/operations.hh"
15 #include "toolbox/sequence.hh"
16 #include "geometry/integral.hh"
17 #include "space/element.hh"
18 #include "space/formula.hh"
19 #include "space/postProcess.hh"
20 #include "space/space.hh"
21 
22 #define ElemFrmIntegrate_D 0
23 #define PWintegrate_D 0
24 #define SpaceIntegrate_D 0
25 #define ElemFrmL2Product_D 0
26 #define SpcFrmL2Product_D 0
27 #define SeqElemFrmL2Product_D 0
28 
29 namespace concepts {
30 
31  // forward declaration
32  template<typename F, typename G>
33  class ElementFormula;
34 
35 
36  // ************************************************************* integrate **
37 
42  template<typename G>
43  Real integrate(const Element<G>& elm)
44  {
45  const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
46  if (!cell)
47  throw conceptsException(MissingFeature("not a integration cell"));
48 
49  Real val(0);
50  uint i = 0;
51  IntegrationCell::intPoint p;
52 
53  while(cell->quadraturePoint(i++, p)) val += p.intermediate;
54 
55  return val;
56  }
57 
63  template<typename F, typename G>
64  F integrate(const ElementWithCell<G>& elm, const ElementFormula<F,G>& frm,
65  const Real t = 0.0,
66  IntegrationCell::intFormType form = IntegrationCell::ZERO)
67  {
68  DEBUGL(ElemFrmIntegrate_D, "elm = " << elm);
69  const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
70  if (!cell)
71  throw conceptsException(MissingFeature("not a integration cell"));
72 
73  F val(0);
74  uint i = 0;
75  IntegrationCell::intPoint p;
76  DEBUGL(ElemFrmIntegrate_D, "elm = " << elm);
77 
78  while(cell->quadraturePoint(i++, p, form, true)) {
79  DEBUGL(ElemFrmIntegrate_D, "coord = " << p.coord << ", frm(coord) = " <<
80  frm(elm, p.coord, t) << ", interme = " << p.intermediate);
81  val += frm(elm, p.coord, t) * p.intermediate;
82  }
83  return val;
84  }
85 
95  template<class F, typename G>
96  F integrate(const SpaceOnCells<G>& spc, const ElementFormula<F,G>& frm,
97  const Real t = 0.0,
98  IntegrationCell::intFormType form = IntegrationCell::ZERO) {
99  F val(0);
100 #if SpaceIntegrate_D
101  uint i = 0;
102 #endif
103  // loop over elements
104  Scan<ElementWithCell<G> >* sc(spc.scan());
105  while (*sc) {
106  const ElementWithCell<G>& elm = (*sc)++;
107 #if SpaceIntegrate_D
108  DEBUGL(1, ++i << ".Element: " << elm);
109  F tmp = integrate(elm, frm, t, form);
110  val += tmp;
111  DEBUGL(1, " Integral: " << tmp);
112 #else
113  val += integrate(elm, frm, t, form);
114 #endif
115  }
116  return val;
117  }
118 
130  template<class F, typename G>
131  F integrate(const Sequence<ElementWithCell<G> * >& elm_seq, const ElementFormula<F,G>& frm,
132  const Real t = 0.0,
133  IntegrationCell::intFormType form = IntegrationCell::ZERO) {
134 
135  F val(0);
136 
137  for (typename Sequence<ElementWithCell<G> *>::iterator it = elm_seq.begin(); it != elm_seq.end(); ++it)
138  {
139  val+= integrate(**it, frm, t, form);
140  }
141 
142  return val;
143 
144 
145  }
153  template<class F, class G>
154  F integrate(const ElementWithCell<G> &elm1, const ElementWithCell<G> &elm2,
155  const ElementFormula<F,G> &frm1, const ElementFormula<F,G> &frm2,
156  const Real t = 0,
157  IntegrationCell::intFormType form = IntegrationCell::ZERO)
158  {
159  const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm1);
160  if (!cell)
161  throw conceptsException(MissingFeature("not a integration cell"));
162 
163  F val(0);
164  uint i = 0;
165  IntegrationCell::intPoint p;
166 
167  while(cell->quadraturePoint(i++, p, form, true)) {
168  val += frm1(elm1, p.coord, t) * frm2(elm2, p.coord,t) * p.intermediate;
169  }
170  return val;
171  }
172 
181  template<class F, class G>
182  F integrate(const SpaceOnCells<G>& spc1, const SpaceOnCells<G>& spc2,
183  const ElementFormula<F,G>& frm1, const ElementFormula<F,G>& frm2,
184  const Real t = 0,
185  IntegrationCell::intFormType form = IntegrationCell::ZERO)
186  {
187  F val(0);
188  Scan<ElementWithCell<G> >* sc1(spc1.scan());
189  while( *sc1 ) {
190  const ElementWithCell<G> &elm1 = (*sc1)++;
191  Scan<ElementWithCell<G> >* sc2(spc2.scan());
192  while( *sc2 ) {
193  const ElementWithCell<G> &elm2 = (*sc2)++;
194  if(elm1.cell().connector() == elm2.cell().connector()) {
195  val += integrate(elm1, elm2, frm1, frm2, t, form);
196  break;
197  }
198  }
199  }
200  return val;
201  }
202 
203  // ************************************************************* L2product **
204 
213  template<typename F, typename G>
214  Real L2product(const ElementWithCell<G>& elm, const ElementFormula<F,G>& u,
215  const ElementFormula<Real>* c = 0, const Real t = 0.0,
216  IntegrationCell::intFormType form = IntegrationCell::ZERO) {
217  const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
218  if (!cell)
219  throw conceptsException(MissingFeature("not a integration cell"));
220 
221  DEBUGL(ElemFrmL2Product_D, "Element = " << elm);
222  DEBUGL(ElemFrmL2Product_D, "u = " << u);
223  Real val(0);
224  uint i = 0;
225  IntegrationCell::intPoint p;
226 
227  while(cell->quadraturePoint(i++, p, form, true)) {
228 #if ElemFrmL2Product_D
229  if (!c) {
230  DEBUGL(1, "coord = " << p.coord << ", u(coord) = " <<
231  u(elm, p.coord, t) << ", interme = " << p.intermediate);
232  } else {
233  DEBUGL(1, "coord = " << p.coord << ", u(coord) = " <<
234  u(elm, p.coord, t) << ", interme = " << p.intermediate <<
235  ", formula = " << (*c)(elm, p.coord, t) << ", contribution = " <<
236  u(elm, p.coord, t)*u(elm, p.coord, t)*
237  p.intermediate*(*c)(elm, p.coord, t));
238  }
239 #endif
240  F f = u(elm, p.coord, t);
241  Real add = std::norm(f) * p.intermediate;
242  if (c) add *= (*c)(elm, p.coord, t);
243  val += add;
244  }
245  // one could introduce special handling for piecewise constant
246  // formulas which might be zero on some elements or formulas whose
247  // support are just some elements
248 
249  DEBUGL(ElemFrmL2Product_D, "val = " << val);
250  return val;
251  }
252 
266  template<class F, typename G>
267  Real L2product(SpaceOnCells<F>& spc, const G& u,
268  const ElementFormula<Real>* c = 0, const Real t = 0.0,
269  IntegrationCell::intFormType form = IntegrationCell::ZERO) {
270  Real val(0);
271  DEBUGL(SpcFrmL2Product_D, "spc = " << spc);
272  DEBUGL(SpcFrmL2Product_D, "u = " << u);
273  // loop over elements
274  std::unique_ptr<Scan<ElementWithCell<F> > > sc(spc.scan());
275  uint j = 0;
276  while (*sc) {
277  const ElementWithCell<F>& elm = (*sc)++;
278  DEBUGL(SpcFrmL2Product_D, ++j << "th element = " << elm.cell());
279  val += L2product(elm, u, c, t, form);
280  DEBUGL(SpcFrmL2Product_D, "val = " << val);
281  }
282  DEBUGL(SpcFrmL2Product_D, "done.");
283  return val;
284  }
285 
295  template<class F, typename G>
296  Real L2product(const Sequence<ElementWithCell<G> * >& elm_seq, const ElementFormula<F,G>& u,
297  const ElementFormula<Real>* c = 0, const Real t = 0.0,
298  IntegrationCell::intFormType form = IntegrationCell::ZERO) {
299 
300  Real val(0);
301  DEBUGL(SeqElemFrmL2Product_D, "Looking over iterators");
302  for (typename Sequence<ElementWithCell<G> * >::const_iterator it = elm_seq.begin(); it != elm_seq.end(); ++it)
303  {
304  DEBUGL(SeqElemFrmL2Product_D, **it);
305  val+= L2product(**it, u, c, t, form);
306  }
307 
308  return val;
309  }
310 
311  // ********************************************************** CellIntegral **
312 
315  template<class F>
316  class CellIntegral : public CellPostprocess<Real> {
317  public:
323  CellIntegral(F& val, const Formula<F>& formula,
324  const std::set<uint>* attributes = 0);
325  protected:
326  virtual std::ostream& info(std::ostream& os) const;
328  F& val_;
330  const concepts::Formula<F>& formula_;
332  const std::set<uint>* attributes_;
333  };
334 
335  // ****************************************************** CellFaceIntegral **
336 
339  template<class F>
340  class CellFaceIntegral : public CellIntegral<F> {
341  public:
347  CellFaceIntegral(F& val, const concepts::Formula<F>& formula,
348  const std::set<uint>* attributes = 0);
349  virtual void operator() (const Element<Real>& elm);
350  virtual void operator() (const Cell& cell);
351 
352  protected:
353  virtual std::ostream& info(std::ostream& os) const;
354  private:
355  template<class G>
356  F compute_(const G& elm);
357  };
358 
359  // ****************************************************** CellEdgeIntegral **
360 
363  template<class F>
364  class CellEdgeIntegral : public CellIntegral<F> {
365  public:
371  CellEdgeIntegral(F& val, const concepts::Formula<F>& formula,
372  const std::set<uint>* attributes = 0);
373  virtual void operator() (const Element<Real>& elm);
374  virtual void operator() (const Cell& cell);
375  protected:
376  virtual std::ostream& info(std::ostream& os) const;
377  private:
378  template<class G>
379  F compute_(const G& elm);
380  };
381 
382 } // namespace concepts
383 
384 #endif // spcIntegral_hh
concepts::integrate
Real integrate(const Element< G > &elm)
Returns the area of the cell belonging to the element elm.
Definition: integral.hh:43
concepts::Cell::connector
virtual Connector & connector() const =0
Returns the connector.
concepts::CellIntegral::val_
F & val_
Integration value.
Definition: integral.hh:328
concepts::IntegrationCell::intPoint
Integration point consisting of coordinates and intermediate data.
Definition: integral.hh:34
conceptsException
#define conceptsException(exc)
Prepares an exception for throwing.
Definition: exceptions.hh:344
concepts::CellEdgeIntegral::CellEdgeIntegral
CellEdgeIntegral(F &val, const concepts::Formula< F > &formula, const std::set< uint > *attributes=0)
Constructor.
concepts::CellEdgeIntegral
Integral over a edge, evaluated on a cell.
Definition: integral.hh:364
concepts::CellIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::Cell
A cell in a mesh consist of topological information (neighbours, connectivity, orientation) and geome...
Definition: cell.hh:39
concepts::IntegrationCell::intFormType
intFormType
Integration form, which determines terms coming from integration over reference element.
Definition: integral.hh:29
concepts::L2product
Real L2product(const ElementWithCell< G > &elm, const ElementFormula< F, G > &u, const ElementFormula< Real > *c=0, const Real t=0.0, IntegrationCell::intFormType form=IntegrationCell::ZERO)
Returns the L2 product or with c weighted L2 product of an element formula u over the cell belonging ...
Definition: integral.hh:214
std::norm
const concepts::Real norm(const concepts::Real &v)
Returns the square of a real value.
Definition: operations.hh:90
concepts::Formula
Interface for a formula.
Definition: lform.hh:18
concepts::CellIntegral::CellIntegral
CellIntegral(F &val, const Formula< F > &formula, const std::set< uint > *attributes=0)
Constructor.
concepts::CellEdgeIntegral::operator()
virtual void operator()(const Element< Real > &elm)
Application operator.
concepts::ElementWithCell
Element with cell.
Definition: elementFormula.hh:19
SeqElemFrmL2Product_D
#define SeqElemFrmL2Product_D
Definition: integral.hh:27
DEBUGL
#define DEBUGL(doit, msg)
ElemFrmIntegrate_D
#define ElemFrmIntegrate_D
Definition: integral.hh:22
concepts::ElementFormula
Interface for a formula defined element by element.
Definition: elementFormula.hh:35
concepts::CellFaceIntegral
Integral over a face, evaluated on a cell.
Definition: integral.hh:340
concepts::CellEdgeIntegral::compute_
F compute_(const G &elm)
concepts::Scan
An abstract class for scanning a mesh (a set of cells) or a space (a set of elements).
Definition: scannerConnectors.hh:29
concepts::CellFaceIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::Sequence
Sequence with operations, output operator, and method of the particular element types.
Definition: sequence.hh:39
ElemFrmL2Product_D
#define ElemFrmL2Product_D
Definition: integral.hh:25
concepts::CellFaceIntegral::compute_
F compute_(const G &elm)
concepts::CellFaceIntegral::CellFaceIntegral
CellFaceIntegral(F &val, const concepts::Formula< F > &formula, const std::set< uint > *attributes=0)
Constructor.
concepts::IntegrationCell
Cell over which can be integrated.
Definition: integral.hh:24
SpcFrmL2Product_D
#define SpcFrmL2Product_D
Definition: integral.hh:26
concepts::CellIntegral
Integral, evaluated on a cell.
Definition: integral.hh:316
concepts::CellIntegral::formula_
const concepts::Formula< F > & formula_
Formula to integrate.
Definition: integral.hh:330
concepts::ElementWithCell::cell
virtual const Cell & cell() const =0
Returns the cell on which the element is built.
concepts::SpaceOnCells
Abstract class for a space.
Definition: space.hh:81
concepts::MissingFeature
Exception class to express a missing feature.
Definition: exceptions.hh:206
concepts::SpaceOnCells::scan
virtual Scanner * scan() const =0
Returns a scanner to iterate over the elements of the space.
concepts::CellIntegral::attributes_
const std::set< uint > * attributes_
Attributes, which contribute.
Definition: integral.hh:332
concepts::CellFaceIntegral::operator()
virtual void operator()(const Element< Real > &elm)
Application operator.
concepts::CellEdgeIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::IntegrationCell::quadraturePoint
virtual bool quadraturePoint(uint i, intPoint &p, intFormType form=ZERO, bool localCoord=false) const =0
Delivers a quadrature point.
concepts::CellPostprocess
Abstract class for per cell postprocessing.
Definition: postProcess.hh:38
concepts::Real
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
concepts
Basic namespace for Concepts-2.
Definition: pml_formula.h:16
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich