You are here: Concepts>Concepts Web>Class documentation

Class documentation of Concepts

Loading...
Searching...
No Matches
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
29namespace 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
conceptsException
#define conceptsException(exc)
Definition exceptions.hh:344
concepts::CellEdgeIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::CellEdgeIntegral::operator()
virtual void operator()(const Element< Real > &elm)
concepts::CellEdgeIntegral::CellEdgeIntegral
CellEdgeIntegral(F &val, const concepts::Formula< F > &formula, const std::set< uint > *attributes=0)
concepts::CellFaceIntegral::operator()
virtual void operator()(const Element< Real > &elm)
concepts::CellFaceIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::CellFaceIntegral::CellFaceIntegral
CellFaceIntegral(F &val, const concepts::Formula< F > &formula, const std::set< uint > *attributes=0)
concepts::CellIntegral::info
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
concepts::CellIntegral::attributes_
const std::set< uint > * attributes_
Attributes, which contribute.
Definition integral.hh:332
concepts::CellIntegral::val_
F & val_
Integration value.
Definition integral.hh:328
concepts::CellIntegral::formula_
const concepts::Formula< F > & formula_
Formula to integrate.
Definition integral.hh:330
concepts::CellIntegral::CellIntegral
CellIntegral(F &val, const Formula< F > &formula, const std::set< uint > *attributes=0)
concepts::Cell::connector
virtual Connector & connector() const =0
Returns the connector.
concepts::ElementWithCell::cell
virtual const Cell & cell() const =0
Returns the cell on which the element is built.
concepts::IntegrationCell::quadraturePoint
virtual bool quadraturePoint(uint i, intPoint &p, intFormType form=ZERO, bool localCoord=false) const =0
concepts::SpaceOnCells::scan
virtual Scanner * scan() const =0
Returns a scanner to iterate over the elements of the space.
DEBUGL
#define DEBUGL(doit, msg)
Definition debug.hh:40
concepts::Real
double Real
Definition typedefs.hh:39
concepts::integrate
Real integrate(const Element< G > &elm)
Definition integral.hh:43
concepts::makeSet
Set< F > makeSet(uint n, const F &first,...)
Definition set.hh:320
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)
Definition integral.hh:214
Generated on Wed Sep 13 2023 21:06:23 for Concepts by doxygen 1.9.8