Quadrature.h 15.4 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file Quadrature.h */

#ifndef AMDIS_QUADRATURE_H
#define AMDIS_QUADRATURE_H

Thomas Witkowski's avatar
Thomas Witkowski committed
26
#include <list>
27
28
#include <vector>
#include <boost/numeric/mtl/mtl.hpp>
Thomas Witkowski's avatar
Thomas Witkowski committed
29
#include "AMDiS_fwd.h"
30
31
32
33
34
35
#include "BasisFunction.h"
#include "Flag.h"
#include "FixVec.h"

namespace AMDiS {

36
37
  using namespace std;

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  /** 
   * \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
54
    /// Avoids call of default constructor
55
56
57
58
59
60
61
62
63
    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_,
64
65
66
	       int degree_,
	       int dim_,
	       int n_points_,
67
	       VectorOfFixVecs<DimVec<double> > *lambda_,
68
	       double* w_)
69
70
71
72
73
74
      : name(name_),
	degree(degree_),
	dim(dim_),
	n_points(n_points_),
	lambda(lambda_),
	w(w_) 
75
    {}
76
77

  public:
78
    /// Copy constructor
79
80
    Quadrature(const Quadrature&);

81
    /// Destructor
82
83
    ~Quadrature();

84
    /// Returns a Quadrature for dimension dim exact for degree degree.
Thomas Witkowski's avatar
Thomas Witkowski committed
85
    static Quadrature *provideQuadrature(int dim, int degree);
86
87
88
89
90
91
92
93
94
95
96
97
98

    /** \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
99
    /// Returns \ref name
100
    inline string getName() 
Thomas Witkowski's avatar
Thomas Witkowski committed
101
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
102
      return name; 
103
    }
104

105
    /// Returns \ref n_points
Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
    inline int getNumPoints() const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      return n_points;
109
    }
110

Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
113
    /// Returns \ref w[p]
    inline double getWeight(int p) const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
114
      return w[p];
115
    }
116

Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
119
    /// Returns \ref w.
    inline double* getWeight() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      return w; 
121
    }
122

Thomas Witkowski's avatar
Thomas Witkowski committed
123
124
125
    /// Returns \ref dim
    inline int getDim() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
126
      return dim; 
127
    }
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
130
131
    /// Returns \ref degree
    inline int getDegree() const 
    { 
Thomas Witkowski's avatar
Thomas Witkowski committed
132
      return degree; 
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
159

    /** \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;
  


160
161
    /// 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
162
163
    inline double getLambda(int a, int b) const 
    {
164
      return (lambda ? (*lambda)[a][b] : 0.0);
165
    }
166

167
168
    /// 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
169
170
    inline const DimVec<double>& getLambda(int a) const 
    {
171
      return (*lambda)[a]; 
172
    }
173

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
201
    /// Vector of quadrature weights
202
203
204
    double *w;

  protected:
205
206
    /// Initialisation of all static Quadrature objects which will be returned
    /// by \ref provideQuadrature()
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
    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
283

284

285
  /// Pre-compute the values of all basis functions at all quadrature nodes;  
286
287
  const Flag INIT_PHI=1;

288
289
  /// Pre-compute the gradients (with respect to the barycentric coordinates) of
  /// all basis functions at all quadrature nodes
290
291
  const Flag INIT_GRD_PHI=2; 

292
293
  /// pre-compute all 2nd derivatives (with respect to the barycentric 
  /// coordinates) of all basis functions at all quadrature nodes;
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
  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:
316
317
    /// Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
    /// flag.
318
    FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag)
319
320
      : init_flag(flag),
	phi(0, 0),
321
322
323
	D2Phi(NULL), 
	quadrature(quad), 
	basisFunctions(basFcts) 
324
    {}
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
    /// Copy constructor
327
328
    FastQuadrature(const FastQuadrature&);

Thomas Witkowski's avatar
Thomas Witkowski committed
329
    /// Extended copy constructor
330
    FastQuadrature(const FastQuadrature&, const Flag);
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
    /// Destructor
333
334
335
    ~FastQuadrature();

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    /// Returns a FastQuadrature for the given BasisFunction, Quadrature, and flags.
337
338
339
340
    static FastQuadrature* provideFastQuadrature(const BasisFunction*,
						 const Quadrature&,
						 Flag);

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// inits FastQuadrature like speciefied in flag
342
343
    void init(Flag init_flag);

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
    inline bool initialized(Flag flag) 
    {
      if (flag == INIT_PHI)
347
	return (num_rows(phi) > 0);
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
      if (flag == INIT_GRD_PHI)
350
	return (!grdPhi.empty());
351

Thomas Witkowski's avatar
Thomas Witkowski committed
352
      if (flag == INIT_D2_PHI)
353
354
355
356
	return (D2Phi != NULL);

      ERROR_EXIT("invalid flag\n");
      return false;
357
    }
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
    /// Returns \ref quadrature
    inline const Quadrature* getQuadrature() const 
    { 
362
      return quadrature; 
363
    }
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
    /// Returns \ref max_points
    inline int getMaxQuadPoints() 
    { 
368
      return max_points; 
369
    }
370

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
    /// Returns (*\ref grdPhi)[q][i][j]
    inline const double getGradient(int q, int i ,int j) const 
    {
380
      return (!grdPhi.empty()) ? grdPhi[q][i][j] : 0.0;
381
    }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    /// Returns (*\ref grdPhi)[q]
384
    inline const vector<mtl::dense_vector<double> >& getGradient(int q) const
Thomas Witkowski's avatar
Thomas Witkowski committed
385
    {
386
      return grdPhi[q];
387
    }
388

389
    inline const mtl::dense_vector<double>& getGradient(int q, int i) const
390
    {
391
      return grdPhi[q][i];
392
393
    }

394
    inline const mtl::dense2D<double>& getPhi() const
Thomas Witkowski's avatar
Thomas Witkowski committed
395
    {
396
      return phi;
397
    }
398

399
400
    /// Returns \ref phi[q][i]
    inline const double getPhi(int q, int i) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
401
    {
402
      return phi[q][i];
403
    }
404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
    /// Returns \ref quadrature ->integrateStdSimplex(f)
    inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) 
    {
408
      return quadrature->integrateStdSimplex(f); 
409
    }
410
  
Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
413
    /// Returns \ref quadrature ->getNumPoints()
    inline int getNumPoints() const 
    { 
414
      return quadrature->getNumPoints();
415
    }
416

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
    /// Returns \ref quadrature ->getWeight(p)
    inline double getWeight(int p) const 
    {
420
      return quadrature->getWeight(p);
421
    }
422

Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
425
    /// Returns \ref quadrature ->getDim()
    inline int getDim() const 
    { 
426
      return quadrature->getDim(); 
427
    }
428

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

Thomas Witkowski's avatar
Thomas Witkowski committed
435
    /// Returns \ref quadrature ->grdFAtQp(f, vec)
436
437
438
439
440
441
    inline const WorldVector<double> 
    *grdFAtQp(const AbstractFunction<WorldVector<double>, 
	      DimVec<double> >& f,
	      WorldVector<double>* vec) const 
    { 
      return quadrature->grdFAtQp(f, vec);
442
    }
443
  
Thomas Witkowski's avatar
Thomas Witkowski committed
444
    /// Returns \ref quadrature ->fAtQp(f, vec)
445
    inline const double *fAtQp(const AbstractFunction<double, 
446
			       DimVec<double> >& f, double *vec) const 
447
448
    {
      return quadrature->fAtQp(f, vec);    
449
    }
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453
    /// Returns \ref quadrature ->getLambda(a,b)
    inline double getLambda(int a,int b) const 
    {
454
      return quadrature->getLambda(a,b);
455
    }
456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
459
    /// Returns \ref quadrature ->getLambda(a)
    inline const DimVec<double>& getLambda(int a) const 
    { 
460
      return quadrature->getLambda(a); 
461
    }
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
465
    /// Returns \ref basisFunctions
    inline BasisFunction* getBasisFunctions() const 
    { 
466
      return basisFunctions; 
467
    }
468
469

  protected:
470
471
    /// Specifies which information should be pre-computed. Can be \ref INIT_PHI, 
    /// \ref INIT_GRD_PHI, or \ref INIT_D2_PHI
472
473
474
475
476
477
478
479
    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
     */
480
    mtl::dense2D<double> phi;
481
482
483
484
485
486
487
488

    /** \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
     */
489
    vector<vector<mtl::dense_vector<double> > > grdPhi; 
490
491
492
493
494
495
496
497
498
499

    /** \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;

500
    /// List of all used FastQuadratures
501
    static list<FastQuadrature*> fastQuadList;
502

503
504
505
    /// Maximal number of quadrature points for all yet initialised FastQuadrature
    /// objects. This value may change after a new initialisation of a
    /// FastQuadrature
506
507
    static int max_points;

Thomas Witkowski's avatar
Thomas Witkowski committed
508
    /// This FastQuadrature stores values for Quadrature quadrature
509
510
    Quadrature* quadrature;

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    /// Values stored for basis functions basisFunctions
512
513
514
515
516
517
518
    BasisFunction* basisFunctions;

  };

}

#endif  // AMDIS_QUADRATURE_H