Flow123d  JS_before_hm-62-ge56f9d5
quadrature_lib.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file quadrature_lib.cc
15  * @brief Definitions of particular quadrature rules on simplices.
16  * @author Jan Stebel
17  */
18 
19 #include "system/global_defs.h"
20 #include "system/system.hh"
22 
23 
24 #define FLOAT double
25 #define SHORT int
26 #define _F(n) n
27 #include "quad.c"
28 #undef FLOAT
29 #undef SHORT
30 #undef _F
31 
32 
33 using namespace arma;
34 
35 
36 QGauss::QGauss(unsigned int dim, const unsigned int order)
37 : Quadrature(dim)
38 {
39  typedef QUAD* pQUAD;
40 
41  static const pQUAD
42  q1d[] = { QUAD_1D_P1,
50  },
51  q2d[] = { QUAD_2D_P1,
58  },
59  q3d[] = { QUAD_3D_P1,
64  };
65  static const double unit_cell_volume[] = { 1, 1, 0.5, 1./6 };
66  const pQUAD *q = nullptr;
67  unsigned int nquads = 0;
68 
69  switch (dim)
70  {
71  case 0:
72  this->quadrature_points.push_back({});
73  this->weights.push_back(1);
74  return;
75  case 1:
76  q = q1d;
77  nquads = sizeof(q1d) / sizeof(pQUAD);
78  break;
79  case 2:
80  q = q2d;
81  nquads = sizeof(q2d) / sizeof(pQUAD);
82  break;
83  case 3:
84  q = q3d;
85  nquads = sizeof(q3d) / sizeof(pQUAD);
86  break;
87  }
88 
89  ASSERT_DBG(order < nquads).error("Quadrature of given order is not implemented.");
90  ASSERT_PTR_DBG(q);
91 
92  vector<double> p(dim, 0);
93  for (int i=0; i<q[order]->npoints; i++)
94  {
95  for (unsigned int j=0; j<dim; j++)
96  p[j] = q[order]->points[i*(dim+1)+j];
97 
99  // The weights must be adjusted according to the volume of the unit cell:
100  // 1D cell: volume 1
101  // 2D cell: volume 1/2
102  // 3D cell: volume 1/6
103  this->weights.push_back(q[order]->weights[i]*unit_cell_volume[dim]);
104  }
105 }
106 
107 
unsigned int dim() const
Definition: quadrature.hh:71
#define QUAD_2D_P20
Definition: quad.h:115
#define QUAD_1D_P14
Definition: quad.h:63
#define QUAD_3D_P5
Definition: quad.h:129
#define QUAD_1D_P10
Definition: quad.h:57
#define QUAD_3D_P7
Definition: quad.h:133
#define QUAD_3D_P6
Definition: quad.h:131
#define QUAD_1D_P15
Definition: quad.h:64
#define QUAD_3D_P8
Definition: quad.h:135
#define QUAD_2D_P9
Definition: quad.h:93
#define QUAD_1D_P11
Definition: quad.h:58
#define QUAD_1D_P2
Definition: quad.h:45
#define QUAD_2D_P21
Definition: quad.h:117
#define QUAD_2D_P17
Definition: quad.h:109
#define QUAD_2D_P10
Definition: quad.h:95
Definition: quad.h:29
#define QUAD_1D_P18
Definition: quad.h:69
#define QUAD_2D_P19
Definition: quad.h:113
#define QUAD_2D_P2
Definition: quad.h:79
void push_back(const std::vector< Type > &p)
Definition: armor.hh:150
#define QUAD_1D_P5
Definition: quad.h:49
#define QUAD_1D_P21
Definition: quad.h:73
#define QUAD_3D_P13
Definition: quad.h:145
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
#define QUAD_2D_P14
Definition: quad.h:103
#define QUAD_2D_P6
Definition: quad.h:87
#define QUAD_2D_P12
Definition: quad.h:99
Global macros to enhance readability and debugging, general constants.
#define QUAD_3D_P2
Definition: quad.h:123
#define QUAD_2D_P8
Definition: quad.h:91
#define QUAD_2D_P16
Definition: quad.h:107
#define QUAD_2D_P1
Definition: quad.h:77
std::vector< double > weights
List of weights to the quadrature points.
Definition: quadrature.hh:141
#define QUAD_1D_P13
Definition: quad.h:61
#define QUAD_1D_P9
Definition: quad.h:55
#define QUAD_2D_P5
Definition: quad.h:85
#define QUAD_2D_P18
Definition: quad.h:111
#define QUAD_2D_P4
Definition: quad.h:83
#define QUAD_1D_P16
Definition: quad.h:66
#define QUAD_1D_P8
Definition: quad.h:54
#define QUAD_3D_P12
Definition: quad.h:143
#define QUAD_1D_P17
Definition: quad.h:67
#define QUAD_2D_P13
Definition: quad.h:101
#define QUAD_3D_P14
Definition: quad.h:147
#define QUAD_3D_P11
Definition: quad.h:141
#define QUAD_1D_P3
Definition: quad.h:46
#define QUAD_3D_P10
Definition: quad.h:139
#define ASSERT_DBG(expr)
QGauss(unsigned int dim, const unsigned int order)
Create a formula of given order.
Definitions of particular quadrature rules on simplices.
#define QUAD_1D_P19
Definition: quad.h:70
#define QUAD_3D_P1
Definition: quad.h:121
#define QUAD_2D_P7
Definition: quad.h:89
#define QUAD_1D_P4
Definition: quad.h:48
#define QUAD_1D_P1
Definition: quad.h:43
#define QUAD_1D_P12
Definition: quad.h:60
#define QUAD_2D_P11
Definition: quad.h:97
#define QUAD_2D_P15
Definition: quad.h:105
#define QUAD_1D_P20
Definition: quad.h:72
#define QUAD_2D_P3
Definition: quad.h:81
#define QUAD_1D_P7
Definition: quad.h:52
#define QUAD_3D_P4
Definition: quad.h:127
#define ASSERT_PTR_DBG(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:340
#define QUAD_3D_P9
Definition: quad.h:137
Armor::array quadrature_points
List of quadrature points.
Definition: quadrature.hh:134
#define QUAD_3D_P3
Definition: quad.h:125
#define QUAD_1D_P6
Definition: quad.h:51