Quadrature.h 16.2 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
// ==                                                                        ==
// ============================================================================

/** \file Quadrature.h */

#ifndef AMDIS_QUADRATURE_H
#define AMDIS_QUADRATURE_H

Thomas Witkowski's avatar
Thomas Witkowski committed
25
26
#include <list>
#include "AMDiS_fwd.h"
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include "BasisFunction.h"
#include "Flag.h"
#include "FixVec.h"

namespace AMDiS {

  /** 
   * \ingroup Assembler
   *
   * \brief
   * For the assemblage of the system matrix and right hand side vector of the 
   * linear system, we have to compute integrals, for example: 
   * \f[ \int_{\Omega} f(x)\varphi_i(x) dx \f]
   * For general data A, b, c, and f, these integrals can not be calculated 
   * exactly. Quadrature formulas have to be used in order to calculate the 
   * integrals approximately. Numerical integration in finite element methods is
   * done by looping over all grid elements and using a quadrature formula on 
   * each element.
   */
  class Quadrature
  {
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
49
    /// Avoids call of default constructor
50
51
52
53
54
55
56
57
58
    Quadrature();

    /** \brief
     * Constructs a Quadrature with name name_ of degree degree_ for dim dim_.
     * The Quadrature has n_points_ quadrature points with barycentric 
     * coordinates lambda_ and weights w_. The constructor is protected because
     * access to a Quadrature should be done via \ref provideQuadrature.
     */
    Quadrature(const char* name_,
59
60
61
	       int degree_,
	       int dim_,
	       int n_points_,
62
	       VectorOfFixVecs<DimVec<double> > *lambda_,
63
	       double* w_)
64
65
66
67
68
69
      : name(name_),
	degree(degree_),
	dim(dim_),
	n_points(n_points_),
	lambda(lambda_),
	w(w_) 
70
    {}
71
72

  public:
73
    /// Copy constructor
74
75
    Quadrature(const Quadrature&);

76
    /// Destructor
77
78
    ~Quadrature();

79
    /// Returns a Quadrature for dimension dim exact for degree degree.
Thomas Witkowski's avatar
Thomas Witkowski committed
80
    static Quadrature *provideQuadrature(int dim, int degree);
81
82
83
84
85
86
87
88
89
90
91
92
93

    /** \brief
     * Approximates an integral by the numerical quadrature described by quad;
     * f is a pointer to an AbstractFunction to be integrated, evaluated in 
     * barycentric coordinates; the return value is 
     * \f[ \sum_{k = 0}^{n_points-1} w[k] * (*f)(lambda[k]) \f]
     * For the approximation of \f$ \int_S f\f$ we have to multiply this value 
     * with d!|S| for a simplex S; for a parametric simplex f should be a pointer
     * to a function which calculates 
     * \f$ f(\lambda)|det DF_S(\hat{x}(\lambda))| \f$.
     */
    double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f);
  
Thomas Witkowski's avatar
Thomas Witkowski committed
94
    /// Returns \ref name
Thomas Witkowski's avatar
Thomas Witkowski committed
95
    inline std::string getName() 
Thomas Witkowski's avatar
Thomas Witkowski committed
96
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
97
      return name; 
98
    }
99

100
    /// Returns \ref n_points
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    inline int getNumPoints() const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
103
      return n_points;
104
    }
105

Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
108
    /// Returns \ref w[p]
    inline double getWeight(int p) const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
      return w[p];
110
    }
111

Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
    /// Returns \ref w.
    inline double* getWeight() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
115
      return w; 
116
    }
117

Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
120
    /// Returns \ref dim
    inline int getDim() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
121
      return dim; 
122
    }
123

Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
    /// Returns \ref degree
    inline int getDegree() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
127
      return degree; 
128
    }
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

    /** \brief
     * Returns a pointer to a vector storing the values of a doubled valued 
     * function at all quadrature points; f is that AbstractFunction
     * , evaluated in barycentric coordinates; if vec is not NULL, the values are
     * stored in this vector, otherwise the values are stored in some static 
     * local vector, which is overwritten on the next call
     */
    const double *fAtQp(const AbstractFunction<double, DimVec<double> >& f,
			double *vec) const ;

    /** \brief
     * Returns a pointer to a vector storing the gradient (with respect to world
     * coordinates) of a double valued function at all quadrature points;
     * grdF is a pointer to a AbstractFunction, evaluated in barycentric
     * coordinates and returning a pointer to a WorldVector storing the gradient;
     * if vec is not NULL, the values are stored in this vector, otherwise the 
     * values are stored in some local static vector, which is overwritten on the
     * next call
     */
    const WorldVector<double> *grdFAtQp(const AbstractFunction<WorldVector<double>, 
					DimVec<double> >& grdF, 
					WorldVector<double>* vec) const;
  


    /** \brief
     * Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th
     * quadrature point
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
    inline double getLambda(int a, int b) const 
    {
161
      return (lambda ? (*lambda)[a][b] : 0.0);
162
    }
163
164
165
166
167

    /** \brief
     * Returns \ref lambda[a] which is a DimVec<double> containing the 
     * coordiantes of the a-th quadrature point
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
168
169
    inline const DimVec<double>& getLambda(int a) const 
    {
170
      return (*lambda)[a]; 
171
    }
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
175
    /// Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
    VectorOfFixVecs<DimVec<double> > *getLambda() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
176
      return lambda; 
177
    }
178
179
180


  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    /// Maximal number of quadrature points for the different dimensions
182
183
184
    static const int maxNQuadPoints[4];

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
185
    /// Name of this Quadrature
186
    std::string name;
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
    /// Quadrature is exact of this degree
189
190
    int degree;

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    /// Quadrature for dimension dim
192
193
    int dim;

194
    /// Number of quadrature points
195
196
    int n_points;

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    /// Vector of quadrature points given in barycentric coordinates
198
199
    VectorOfFixVecs<DimVec<double> > *lambda; 

Thomas Witkowski's avatar
Thomas Witkowski committed
200
    /// Vector of quadrature weights
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    double *w;

  protected:
    /** \brief
     * Initialisation of all static Quadrature objects which will be returned
     * by \ref provideQuadrature()
     */
    static void initStaticQuadratures();

    /** \name static quadratures, used weights, and barycentric coords
     * \{
     */
    static Quadrature **quad_nd[4];
    static Quadrature *quad_0d[1];
    static Quadrature *quad_1d[20];
    static Quadrature *quad_2d[18];
    static Quadrature *quad_3d[8];

    static VectorOfFixVecs<DimVec<double> > *x_0d;
    static double *w_0d;

    static VectorOfFixVecs<DimVec<double> > *x0_1d;
    static VectorOfFixVecs<DimVec<double> > *x1_1d;
    static VectorOfFixVecs<DimVec<double> > *x2_1d;
    static VectorOfFixVecs<DimVec<double> > *x3_1d;
    static VectorOfFixVecs<DimVec<double> > *x4_1d;
    static VectorOfFixVecs<DimVec<double> > *x5_1d;
    static VectorOfFixVecs<DimVec<double> > *x6_1d;
    static VectorOfFixVecs<DimVec<double> > *x7_1d;
    static VectorOfFixVecs<DimVec<double> > *x8_1d;
    static VectorOfFixVecs<DimVec<double> > *x9_1d;
    static double *w0_1d;
    static double *w1_1d;
    static double *w2_1d;
    static double *w3_1d;
    static double *w4_1d;
    static double *w5_1d;
    static double *w6_1d;
    static double *w7_1d;
    static double *w8_1d;
    static double *w9_1d;

    static VectorOfFixVecs<DimVec<double> > *x1_2d;
    static VectorOfFixVecs<DimVec<double> > *x2_2d;
    static VectorOfFixVecs<DimVec<double> > *x3_2d;
    static VectorOfFixVecs<DimVec<double> > *x4_2d;
    static VectorOfFixVecs<DimVec<double> > *x5_2d;
    static VectorOfFixVecs<DimVec<double> > *x7_2d;
    static VectorOfFixVecs<DimVec<double> > *x8_2d;
    static VectorOfFixVecs<DimVec<double> > *x9_2d;
    static VectorOfFixVecs<DimVec<double> > *x10_2d;
    static VectorOfFixVecs<DimVec<double> > *x11_2d;
    static VectorOfFixVecs<DimVec<double> > *x12_2d;
    static VectorOfFixVecs<DimVec<double> > *x17_2d;
    static double *w1_2d;
    static double *w2_2d;
    static double *w3_2d;
    static double *w4_2d;
    static double *w5_2d;
    static double *w7_2d;
    static double *w8_2d;
    static double *w9_2d;
    static double *w10_2d;
    static double *w11_2d;
    static double *w12_2d;
    static double *w17_2d;

    static VectorOfFixVecs<DimVec<double> > *x1_3d;
    static VectorOfFixVecs<DimVec<double> > *x2_3d;
    static VectorOfFixVecs<DimVec<double> > *x3_3d;
    static VectorOfFixVecs<DimVec<double> > *x4_3d;
    static VectorOfFixVecs<DimVec<double> > *x5_3d;
    static VectorOfFixVecs<DimVec<double> > *x7_3d;
    static double *w1_3d;
    static double *w2_3d;
    static double *w3_3d;
    static double *w4_3d;
    static double *w5_3d;
    static double *w7_3d;

    /** \} */
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
284

285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335

  /** \brief
   * Pre-compute the values of all basis functions at all quadrature nodes;  
   */
  const Flag INIT_PHI=1;

  /** \brief
   * Pre-compute the gradients (with respect to the barycentric coordinates) of
   * all basis functions at all quadrature nodes
   */
  const Flag INIT_GRD_PHI=2; 

  /** \brief
   * pre-compute all 2nd derivatives (with respect to the barycentric 
   * coordinates) of all basis functions at all quadrature nodes;
   */
  const Flag INIT_D2_PHI=4;


  /** 
   * \ingroup Integration
   *
   *\brief
   * Often numerical integration involves basis functions, such as the assembling
   * of the system matrix and right hand side, or the integration of finite 
   * element functions. Since numerical quadrature involves only the values at 
   * the quadrature points and the values of basis functions and its derivatives
   * are the same at these points for all elements of the grid, such routines can
   * be much more efficient, if they can use pre-computed values of the basis
   * functions at the quadrature points. In this case the basis functions do not 
   * have to be evaluated for each quadrature point on every element newly.
   * Information that should be pre-computed can be specified by the following 
   * symbolic constants:
   * \ref INIT_PHI, \ref INIT_GRD_PHI, \ref INIT_D2_PHI
   */
  class FastQuadrature
  {
  protected:
    /** \brief
     * Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
     * flag.
     */
    FastQuadrature(BasisFunction* basFcts, 
		   Quadrature* quad, 
		   Flag flag)
      : init_flag(flag), 
	phi(NULL), 
	grdPhi(NULL), 
	D2Phi(NULL), 
	quadrature(quad), 
	basisFunctions(basFcts) 
336
    {}
337

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    /// Copy constructor
339
340
    FastQuadrature(const FastQuadrature&);

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// Extended copy constructor
342
    FastQuadrature(const FastQuadrature&, const Flag);
343

Thomas Witkowski's avatar
Thomas Witkowski committed
344
    /// Destructor
345
346
347
    ~FastQuadrature();

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    /// Returns a FastQuadrature for the given BasisFunction, Quadrature, and flags.
349
350
351
352
    static FastQuadrature* provideFastQuadrature(const BasisFunction*,
						 const Quadrature&,
						 Flag);

Thomas Witkowski's avatar
Thomas Witkowski committed
353
    /// inits FastQuadrature like speciefied in flag
354
355
    void init(Flag init_flag);

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
    inline bool initialized(Flag flag) 
    {
      if (flag == INIT_PHI)
359
360
	return (phi != NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
361
      if (flag == INIT_GRD_PHI)
362
363
	return (grdPhi != NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
364
      if (flag == INIT_D2_PHI)
365
366
367
368
	return (D2Phi != NULL);

      ERROR_EXIT("invalid flag\n");
      return false;
369
    }
370

Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
    /// Returns \ref quadrature
    inline const Quadrature* getQuadrature() const 
    { 
374
      return quadrature; 
375
    }
376

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
    /// Returns \ref max_points
    inline int getMaxQuadPoints() 
    { 
380
      return max_points; 
381
    }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    /// Returns (*\ref D2Phi)[q][i][j][m]
384
    const double getSecDer(int q, int i, int j, int m) const;
385

Thomas Witkowski's avatar
Thomas Witkowski committed
386
    /// Returns (*\ref D2Phi)[q]
387
388
    const VectorOfFixVecs<DimMat<double> > *getSecDer(int q) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
391
    /// Returns (*\ref grdPhi)[q][i][j]
    inline const double getGradient(int q, int i ,int j) const 
    {
392
      return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
393
    }
394

Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
397
    /// Returns (*\ref grdPhi)[q]
    inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const 
    {
398
      return (grdPhi) ? &((*grdPhi)[q]) : NULL;
399
    }
400

401
402
403
404
405
    inline DimVec<double>& getGradient(int q, int i) const
    {
      return (*grdPhi)[q][i];
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
406
407
408
    /// Returns \ref phi[q][i]
    inline const double getPhi(int q, int i) const 
    {
409
      return (phi) ? phi[q][i] : 0.0;
410
    }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
414
    /// Returns \ref phi[q]
    inline const double *getPhi(int q) const 
    {
415
      return (phi) ? phi[q] : NULL;
416
    }
417

Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
420
    /// Returns \ref quadrature ->integrateStdSimplex(f)
    inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) 
    {
421
      return quadrature->integrateStdSimplex(f); 
422
    }
423
  
Thomas Witkowski's avatar
Thomas Witkowski committed
424
425
426
    /// Returns \ref quadrature ->getNumPoints()
    inline int getNumPoints() const 
    { 
427
      return quadrature->getNumPoints();
428
    }
429

Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
432
    /// Returns \ref quadrature ->getWeight(p)
    inline double getWeight(int p) const 
    {
433
      return quadrature->getWeight(p);
434
    }
435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
    /// Returns \ref quadrature ->getDim()
    inline int getDim() const 
    { 
439
      return quadrature->getDim(); 
440
    }
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
444
    /// Returns \ref quadrature ->getDegree()
    inline int getDegree() const 
    { 
445
      return quadrature->getDegree();
446
    }
447

Thomas Witkowski's avatar
Thomas Witkowski committed
448
    /// Returns \ref quadrature ->grdFAtQp(f, vec)
449
450
451
452
453
454
    inline const WorldVector<double> 
    *grdFAtQp(const AbstractFunction<WorldVector<double>, 
	      DimVec<double> >& f,
	      WorldVector<double>* vec) const 
    { 
      return quadrature->grdFAtQp(f, vec);
455
    }
456
  
Thomas Witkowski's avatar
Thomas Witkowski committed
457
    /// Returns \ref quadrature ->fAtQp(f, vec)
458
    inline const double *fAtQp(const AbstractFunction<double, 
459
			       DimVec<double> >& f, double *vec) const 
460
461
    {
      return quadrature->fAtQp(f, vec);    
462
    }
463

Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
466
    /// Returns \ref quadrature ->getLambda(a,b)
    inline double getLambda(int a,int b) const 
    {
467
      return quadrature->getLambda(a,b);
468
    }
469

Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
    /// Returns \ref quadrature ->getLambda(a)
    inline const DimVec<double>& getLambda(int a) const 
    { 
473
      return quadrature->getLambda(a); 
474
    }
475

Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
    /// Returns \ref basisFunctions
    inline BasisFunction* getBasisFunctions() const 
    { 
479
      return basisFunctions; 
480
    }
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517

  protected:
    /** \brief
     * Specifies which information should be pre-computed. Can be \ref INIT_PHI, 
     * \ref INIT_GRD_PHI, or \ref INIT_D2_PHI
     */
    Flag init_flag;

    /** \brief
     * Matrix storing function values if the flag \ref INIT_PHI is set;
     * phi[i][j] stores the value \ref basisFunctions->phi[j]
     * (quadrature->lambda[i]), 0 <= j < basisFunctions->getNumber()  and 
     * 0 <= i < n_points
     */
    double **phi;

    /** \brief
     * Matrix storing all gradients (with respect to the barycentric coordinates)
     * if the flag \ref INIT_GRD_PHI is set; grdPhi[i][j][k] stores the value 
     * basisFunctions->grdPhi[j](quadrature->lambda[i])[k] 
     * for 0 <= j < basisFunctions->getNumber(),
     * 0 <= i < . . . , n_points, and 0 <= k < DIM
     */
    MatrixOfFixVecs<DimVec<double> > *grdPhi;

    /** \brief
     * Matrix storing all second derivatives (with respect to the barycentric
     * coordinates) if the flag \ref INIT_D2_PHI is set; D2Phi[i][j][k][l] stores
     * the value basisFunctions->D2Phi[j](quadrature->lambda[i])[k][l] 
     * for 0 <= j < basisFunctions->getNumber(), 
     * 0 <= i < n_points, and 0 <= k,l < DIM
     */
    MatrixOfFixVecs<DimMat<double> > *D2Phi;

    /** \brief
     * List of all used FastQuadratures
     */
518
    static std::list<FastQuadrature*> fastQuadList;
519
520
521
522
523
524
525
526

    /** \brief
     * Maximal number of quadrature points for all yet initialised FastQuadrature
     * objects. This value may change after a new initialisation of a
     * FastQuadrature
     */
    static int max_points;

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    /// This FastQuadrature stores values for Quadrature quadrature
528
529
    Quadrature* quadrature;

Thomas Witkowski's avatar
Thomas Witkowski committed
530
    /// Values stored for basis functions basisFunctions
531
532
533
534
535
536
537
    BasisFunction* basisFunctions;

  };

}

#endif  // AMDIS_QUADRATURE_H