quadRule.hh

Go to the documentation of this file.
1 
6 #ifndef quadRule_hh
7 #define quadRule_hh
8 
9 #include <map>
10 
11 #include "basics/typedefs.hh"
12 #include "basics/defines.hh"
13 #include "basics/outputOperator.hh"
15 #include "quadrature.hh"
16 #include "geometry/cell2D.hh"
17 #include "toolbox/hashMap.hh"
18 
19 //Precision libs small abcissas
20 //#include <gmp.h>
21 //#include <gmpxx.h>
22 
23 namespace concepts {
24 
25  // ******************************************************** QuadratureRule **
26 
30  class QuadratureRule1d : public virtual OutputOperator {
31  private:
34  public:
36  virtual ~QuadratureRule1d() { }
37 
39  virtual const Real* abscissas() const = 0;
40 
42  virtual const Real* weights() const = 0;
43 
45  virtual uint n() const = 0;
46 
48  virtual void printRule();
49  };
50 
51  // ************************************************* QuadratureRule1dDynamic **
52 
59  public:
62  virtual const Real* abscissas() const { return abscissas_; }
63  virtual const Real* weights() const { return weights_; }
64  protected:
66  const Real* abscissas_;
68  const Real* weights_;
69  };
70 
71  // ******************************************** QuadratureRule1dGaussLobatto **
72 
95  public:
102  weights_ = rule_.weights();
103  }
105  virtual uint n() const { return rule_.n(); }
106  protected:
107  virtual std::ostream& info(std::ostream& os) const;
108  private:
110  };
111 
112  // ********************************************* QuadratureRule1dGaussJacobi **
113 
136  public:
142  weights_ = rule_.weights();
143  }
145  virtual uint n() const { return rule_.n(); }
146  protected:
147  virtual std::ostream& info(std::ostream& os) const;
148  private:
150  };
151 
152  // ************************************************** QuadratureRule1dTrapeze **
153 
154  /* Trapeze quadrature rule including both endpoints -1 and 1.
155  \f[ \int_{-1}^1 f(x) \, dx \approx \sum_{i=0}^(n-1) w_i f(x_i) \f]
156  is only exact for \f$f \in P_{1}\f$.
157  n must be greater or equal to 2.
158  \n
159  The abscissas \f$x_i = -1 + 2i/(n-1)\f$ are equidistant
160  and the weights are
161  \f[w_i = \frac{1}{n-1}\left\{\begin{array}{ll}1 & i \in \{0,n\}\\2 & \mbox{otherwise}\end{array}\right.\f]
162 
163  The computations and the storage of the values are done by the
164  class Quadrature with template parameter 4. The difference
165  between this class and Quadrature is that it is in a class
166  hierarchy of quadrature rules. This has advantages when
167  dynamically switching quadrature rules is needed. On the other
168  hand, this class returns the values via a virtual function call
169  \c abscissas() and \c weights() should therefore not be called
170  to often (inside loops etc.).
171 
172  The abscissas are good for graphical output.
173 
174  @see Quadrature
175  @author Kersten Schmidt, 2004
176  */
177 
179  public:
184  : rule_(n), shAbscissas_(0), shWeights_(0) {
186  weights_ = rule_.weights();
187  }
200 
201  virtual const Real* abscissas() const {
203  }
204  virtual const Real* weights() const {
205  return weights_ ? weights_ : shWeights_;
206  }
207 
208  virtual uint n() const { return rule_.n(); }
209  protected:
210  virtual std::ostream& info(std::ostream& os) const;
211  private:
215  };
216 
217  // ******************************************************* QuadRuleFactory **
218 
224  class QuadRuleFactory : public virtual OutputOperator {
225  public:
226  QuadRuleFactory(enum intRule type = GAUSS_JACOBI, uint constPoints = 10,
227  uint addPoints = 2, bool constant = false);
228 
235  QuadratureRule1d* operator()(const ushort p = 1) const;
236 
243  void set(enum concepts::intRule rule, bool constant, uint points);
244 
246  void reset();
247 
249  inline const uint count() const { return cnt_; }
250 
252  const std::string integrationRule() const;
253 
255  enum intRule type() const { return integrationType_; }
256  protected:
257  virtual std::ostream& info(std::ostream& os) const;
258  private:
264  // TODO: adapt documentation: it should be set(), not setIntegrationRule()
278 
280  uint cnt_;
281  };
282 
283 
284 
285  // ********************************************************* QuadratureRule **
286 
296  class QuadratureRule : public virtual OutputOperator {
297  private:
300  public:
302  virtual ~QuadratureRule() { }
303 
308  virtual const Real* abscissas(uint i) const = 0;
309 
314  virtual const Real* weights(uint i) const = 0;
315 
320  virtual const uint n(uint i=0) const = 0 ;
321 
322  //delivers ith quadraturepoint q on reference cell [-1,1]^n and ith weight, must be implemented in child classes
327  virtual bool quadratureData(const uint i, Real3d& q, Real& w) const = 0;
328  };
329 
330  // ******************************************************* QuadratureRule2d **
331 
348  private:
351 
352  public:
354  virtual ~QuadratureRule2d() { }
355 
360  virtual const Real* abscissas(uint i) const{
362  return i ? intY_ : intX_;
363  }
364 
365  virtual const Real* weights(uint i = 0) const = 0;
366  virtual const uint n(uint i = 0) const = 0;
367 
373  virtual const bool domain() const = 0 ;
374 
378  virtual const bool tensor() const = 0 ;
379 
380  protected:
382  const Real* intX_;
383  const Real* intY_;
384  };
385 
386  // ********************************************* QuadratureRule2dQuadTensor **
387 
398  public:
406 
411  virtual const Real* weights(uint i) const ;
412 
417  virtual const uint n(uint i) const;
418 
419  virtual bool quadratureData(const uint i, Real3d& q, Real& w) const;
420 
421  //tensored rules are supposed to be computed on [-1,1]^2
422  virtual const bool domain() const { return true; }
423 
424  // returns true as class is tensor structured
425  virtual const bool tensor() const { return true; }
426 
427  //returns 1d quadrature rule in x direction
428  //@pre quadrature rule exists by constructor
429  const concepts::QuadratureRule1d* xRule() const{return xRule_.get();}
430 
431  //returns 1d quadrature rule in y direction
432  //@pre quadrature rule exists by constructor
433  const concepts::QuadratureRule1d* yRule() const{return yRule_.get();}
434 
435 
436  protected:
437  virtual std::ostream& info(std::ostream& os) const;
438  private:
440  std::unique_ptr<const concepts::QuadratureRule1d> xRule_;
441  std::unique_ptr<const concepts::QuadratureRule1d> yRule_;
442  };
443 
444  // ********************************************** QuadratureRule2dQuadDuffy **
445 
482  public:
496  QuadratureRule2dQuadDuffy(Real beta, uint nx, uint ny, uint noVtx = 0);
497 
499 
503  virtual const Real* weights(uint i = 0) const{
505  return weights_;
506  }
507 
511  const uint n(uint i = 0) const{
513  return n_;
514  }
515 
516  // returns i-th quadrature point \c q = (x_i,y_i) and i-th weight \c w
517  virtual bool quadratureData(const uint i, Real3d& q, Real& w) const;
518 
519  // abscissas and weights are computed in [0,1]^2
520  virtual const bool domain() const { return false; }
521 
522  // non-tensor structure
523  virtual const bool tensor() const { return false; }
524 
525  // prints the quadrature points
526  void printData() const;
527 
528  protected:
529  virtual std::ostream& info(std::ostream& os) const;
530 
531  private :
545  void compute_(const Real2d S,const Real2d B,const Real2d C,
546  const Real beta, const uint nx, const uint ny,
547  Real*& abscXP, Real*& abscYP, Real*& wP);
548 
549  // weights for each integration Point
550  const Real* weights_;
551  //duffy transformation parameter
553  //number of quadrature points n_ = 2*(nx+ny)
554  uint n_;
555  };
556 
557  // ************************************************** QuadRuleFactoryBase2d **
558 
567  class QuadRuleFactoryBase2d: public virtual OutputOperator {
568 
569  private :
572  public :
581  const ushort* p) const = 0;
582 
583  // print of the quadrature rule
584  virtual const std::string integrationRule() const = 0;
585  };
586 
587 
588  // ************************************************ QuadRuleFactoryTensor2d **
589 
593  //this class is a direct copy of the QuadRuleFactory
595  public:
596  QuadRuleFactoryTensor2d(enum intRule type = GAUSS_JACOBI, uint constPoints = 10,
597  uint addPoints = 2, bool constant = false);
598 
599  virtual QuadratureRule2d* operator()(const concepts::QuadNd& cell, const ushort* p) const;
600 
601 
608  virtual void setTensor(enum concepts::intRule rule, bool constant,
609  uint points);
610 
612  void reset();
613 
615  inline const uint count() const {return cnt_; }
616 
618  enum intRule type() const {return integrationType_; }
619 
621  virtual const std::string integrationRule() const;
622 
623  protected:
624  virtual std::ostream& info(std::ostream& os) const;
625 
626  // Tensor information datas
627 
633  // TODO: adapt documentation: it should be set(), not setIntegrationRule()
647 
649  uint cnt_;
650  };
651 
652  // ************************************************ QuadRuleFactoryTensorDuffy2d **
653 
671  private:
673  struct DuffyData : public OutputOperator {
674  // duffy parameter
676  // number of integration points for underlying tensor integration rule.
677  uint npx;
678  uint npy;
679 
680  DuffyData(const Real beta = 1 ,const uint npx = 20, const uint npy = 20) : beta(beta), npx(npx), npy(npy) {}
681 
683  if (this != &dd){
684  this->beta = dd.beta;
685  this->npx = dd.npx;
686  this->npy = dd.npy;
687  }
688  return *this;
689  }
690 
691  protected:
692  virtual std::ostream& info(std::ostream& os) const {
693  return os << "DuffyData{beta = "<< beta <<" , (npx,npy) = ("<< npx << " , " << npy << ")}";
694  }
695  };
696 
697  public:
702  QuadRuleFactoryTensorDuffy2d(enum intRule tensorType = GAUSS_JACOBI, uint constPointsT = 10,
703  uint addPointsT = 2, bool constantT = false);
704 
705 
714  inline void addSingularPoint(const uint vAttr,const Real beta, uint npx = 20 , uint npy = 20){
715  vtxData_[vAttr]= DuffyData(beta,npx,npy);
716  }
717 
722  virtual QuadratureRule2d* operator()(const concepts::QuadNd& cell, const ushort* p) const;
723 
725  virtual const std::string integrationRule() const;
726  protected:
727  virtual std::ostream& info(std::ostream& os) const;
728 
729  // data from each attribute to its duffy quadrature points necessarc data, i.e. beta, npx, npy
731  };
732 
733 
734 } // namespace concepts
735 
736 #endif // quadRule_hh
virtual bool quadratureData(const uint i, Real3d &q, Real &w) const
Method delivers the i-th quadrature point.
Tensor quadrature rule in 2d.
Definition: quadRule.hh:397
QuadRuleFactoryBase2d(const QuadRuleFactoryBase2d &other)
QuadratureRule(const QuadratureRule1d &other)
virtual const Real * abscissas() const
Returns a pointer into the array of the abscissas.
Definition: quadRule.hh:62
virtual const Real * abscissas() const
Returns a pointer into the array of the abscissas.
Definition: quadRule.hh:201
@ GAUSS_JACOBI
Definition: defines.hh:13
QuadratureRule1dTrapeze(uint n)
Constructor.
Definition: quadRule.hh:183
virtual bool quadratureData(const uint i, Real3d &q, Real &w) const =0
Method delivers the i-th quadrature point.
virtual const Real * weights(uint i=0) const =0
Weight return function.
virtual QuadratureRule2d * operator()(const concepts::QuadNd &cell, const ushort *p) const =0
Returns the quadrature rule, i.e.
uint n() const
Returns the number of quadrature points.
Definition: quadrature.hh:109
const uint count() const
Returns counter of changes.
Definition: quadRule.hh:249
enum intRule type() const
Returns the integration type.
Definition: quadRule.hh:618
virtual QuadratureRule2d * operator()(const concepts::QuadNd &cell, const ushort *p) const
Application operator that returns setted integration rule on a given cell for requested polynomial de...
QuadratureRule1d & operator=(const QuadratureRule1d &other)
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
DuffyData(const Real beta=1, const uint npx=20, const uint npy=20)
Definition: quadRule.hh:680
Class for creation of a quadrature rule.
Definition: quadRule.hh:224
const uint count() const
Returns counter of changes.
Definition: quadRule.hh:615
uint constNumerOfPoints_
Number of integration points to use when constant number is requested.
Definition: quadRule.hh:638
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
virtual const bool tensor() const
Method delivers information about the quadrature rule structure.
Definition: quadRule.hh:425
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
const Real * abscissas_
Abscissas.
Definition: quadRule.hh:66
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
concepts::HashMap< DuffyData > vtxData_
Definition: quadRule.hh:730
virtual const Real * abscissas() const =0
Returns a pointer into the array of the abscissas.
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
virtual const Real * weights(uint i=0) const
Returns the i-th weight, belonging to (x_i,y_i) as nontensored.
Definition: quadRule.hh:503
const concepts::QuadratureRule1d * yRule() const
Definition: quadRule.hh:433
virtual const Real * weights() const
Returns a pointer into the array of the weights.
Definition: quadRule.hh:204
virtual QuadratureRule2d * operator()(const concepts::QuadNd &cell, const ushort *p) const
Returns the quadrature rule, i.e.
virtual const Real * weights(uint i) const =0
Weight return function.
bool useConstantNumberOfPoints_
Use constant number of integration points (true) or not (false).
Definition: quadRule.hh:646
const Real * weights() const
Returns a pointer into the array of the weights.
Definition: quadrature.hh:107
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
std::unique_ptr< const concepts::QuadratureRule1d > yRule_
Definition: quadRule.hh:441
const Real * intX_
Abscissas.
Definition: quadRule.hh:382
QuadratureRule1d(const QuadratureRule1d &other)
QuadRuleFactory(enum intRule type=GAUSS_JACOBI, uint constPoints=10, uint addPoints=2, bool constant=false)
QuadRuleFactoryTensorDuffy2d(enum intRule tensorType=GAUSS_JACOBI, uint constPointsT=10, uint addPointsT=2, bool constantT=false)
Constructor setting the information for the tensor quadrature.
Gauss Jacobi quadrature rule not including both endpoints.
Definition: quadRule.hh:135
Abstract class for quadrature rule factories in 2D.
Definition: quadRule.hh:567
QuadratureRule2dQuadTensor(const concepts::QuadratureRule1d *xRule, const concepts::QuadratureRule1d *yRule)
Constructor.
virtual const uint n(uint i) const
Returns the number of quadrature points in i-th direction.
virtual const Real * abscissas(uint i) const
Returns the quadrature abcissas in the i-th direction.
Definition: quadRule.hh:360
const uint n(uint i=0) const
Returns the number of quadrature points.
Definition: quadRule.hh:511
virtual const std::string integrationRule() const
Returns information on the settings of the quadrature rule.
const Real * weights_
Weights.
Definition: quadRule.hh:68
QuadratureRule2d(const QuadratureRule2d &other)
Exception class for assertions.
Definition: exceptions.hh:258
virtual uint n() const
Returns the number of points.
Definition: quadRule.hh:145
const Real * abscissas() const
Returns a pointer into the array of the abscissas.
Definition: quadrature.hh:105
Base class for a quadrilateral in any dimension.
Definition: cell2D.hh:187
QuadratureRule1dTrapeze(uint n, Real x0, Real xend)
Constructor.
uint addNumberOfPoints_
Number of integration points to add to approximation order when varying number of integration points ...
Definition: quadRule.hh:273
uint cnt_
Counter for changes.
Definition: quadRule.hh:649
std::unique_ptr< const concepts::QuadratureRule1d > xRule_
The integration rules along each direction.
Definition: quadRule.hh:440
const std::string integrationRule() const
Returns information on the settings of the quadrature rule.
virtual const uint n(uint i=0) const =0
Returns the number of quadrature Points.
virtual const bool domain() const
Method delivers if the integration points and weights are computed in [-1,1]^2 or in [0,...
Definition: quadRule.hh:422
QuadratureRule1dGaussLobatto(uint n)
Constructor.
Definition: quadRule.hh:100
virtual const bool domain() const =0
Method delivers if the integration points and weights are computed in [-1,1]^2 or in [0,...
enum intRule integrationType_
Default behaviour: integration rule Gauss Jacobi (highest order).
Definition: quadRule.hh:634
Real * shAbscissas_
Abscissas and weights on the shifted interval.
Definition: quadRule.hh:214
uint cnt_
Counter for changes.
Definition: quadRule.hh:280
Gauss Lobatto quadrature rule including both endpoints.
Definition: quadRule.hh:94
QuadRuleFactoryBase2d & operator=(const QuadRuleFactoryBase2d &other)
virtual uint n() const
Returns the number of points.
Definition: quadRule.hh:105
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: quadRule.hh:692
enum intRule type() const
Returns the integration type.
Definition: quadRule.hh:255
virtual const std::string integrationRule() const =0
const concepts::QuadratureRule1d * xRule() const
Definition: quadRule.hh:429
intRule
Types of integration rules to choose from.
Definition: defines.hh:13
virtual const bool domain() const
Method delivers if the integration points and weights are computed in [-1,1]^2 or in [0,...
Definition: quadRule.hh:520
virtual const Real * weights() const
Returns a pointer into the array of the weights.
Definition: quadRule.hh:63
QuadratureRule2d & operator=(const QuadratureRule2d &other)
double beta(const double a, const double b)
QuadratureRule1d * operator()(const ushort p=1) const
Returns the quadrature rule, i.e.
uint addNumberOfPoints_
Number of integration points to add to approximation order when varying number of integration points ...
Definition: quadRule.hh:642
QuadratureRule1dDynamic(const Real *abscissas=0, const Real *weights=0)
Definition: quadRule.hh:60
virtual void printRule()
print weights and abscissas to stdout
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
unsigned short ushort
Abbreviation for unsigned short.
Definition: typedefs.hh:48
void reset()
Set the standard type of integration.
virtual const bool tensor() const
Method delivers information about the quadrature rule structure.
Definition: quadRule.hh:523
Duffy information data, i.e. the duffy parameter and the number of points.
Definition: quadRule.hh:673
virtual uint n() const =0
Returns the number of points.
Class representing the Generalized Duffy quadrature rule in 2d, , that is for the standard duffy int...
Definition: quadRule.hh:481
void addSingularPoint(const uint vAttr, const Real beta, uint npx=20, uint npy=20)
Adds information for a given singularity vertex attribute, and specified duffy quadrature on it.
Definition: quadRule.hh:714
QuadRuleFactoryTensor2d(enum intRule type=GAUSS_JACOBI, uint constPoints=10, uint addPoints=2, bool constant=false)
Quadrature rule for numerical integration.
Definition: quadRule.hh:30
virtual uint n() const
Returns the number of points.
Definition: quadRule.hh:208
Class representing a quadrature factory, that holds information about cells on which generalized duff...
Definition: quadRule.hh:670
void reset()
Set the standard type of integration.
virtual void setTensor(enum concepts::intRule rule, bool constant, uint points)
Sets the integration rule for cells with tensor rule.
QuadratureRule & operator=(const QuadratureRule1d &other)
void set(enum concepts::intRule rule, bool constant, uint points)
Sets the integration rule.
bool useConstantNumberOfPoints_
Use constant number of integration points (true) or not (false).
Definition: quadRule.hh:277
enum intRule integrationType_
Default behaviour: integration rule Gauss Jacobi (highest order).
Definition: quadRule.hh:265
uint constNumerOfPoints_
Number of integration points to use when constant number is requested.
Definition: quadRule.hh:269
virtual const std::string integrationRule() const
Returns information on the settings of the quadrature rule.
Abstract class for quadrature rules in.
Definition: quadRule.hh:347
Abtract class for quadrature rules in .
Definition: quadRule.hh:296
Class providing an output operator.
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
virtual const bool tensor() const =0
Method delivers information about the quadrature rule structure.
void compute_(const Real2d S, const Real2d B, const Real2d C, const Real beta, const uint nx, const uint ny, Real *&abscXP, Real *&abscYP, Real *&wP)
Routine computes generalized duffy quadrature points on the triangle conv(S,B,C), where S is the coor...
virtual const uint n(uint i=0) const =0
Returns the number of quadrature Points.
virtual bool quadratureData(const uint i, Real3d &q, Real &w) const
Method delivers the i-th quadrature point.
This class is the same as QuadRuleFactory, but returning integration rules in 2d.
Definition: quadRule.hh:594
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
DuffyData & operator=(const DuffyData &dd)
Definition: quadRule.hh:682
Basic namespace for Concepts-2.
Definition: pml_formula.h:16
virtual const Real * abscissas(uint i) const =0
Returns the quadrature abcissas in the i-th direction.
virtual const Real * weights() const =0
Returns a pointer into the array of the weights.
QuadratureRule2dQuadDuffy(Real beta, uint nx, uint ny, uint noVtx=0)
Constructor.
virtual const Real * weights(uint i) const
Returns the weights in i-th direction.
Base class for quadrature rules with dynamically allocated storage for the weights and abscissas.
Definition: quadRule.hh:58
QuadratureRule1dGaussJacobi(uint n)
Constructor.
Definition: quadRule.hh:140
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich