DOLFIN
DOLFIN C++ interface
SimplexQuadrature.h
1 // Copyright (C) 2014 Anders Logg
2 //
3 // This file is part of DOLFIN.
4 //
5 // DOLFIN is free software: you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published by
7 // the Free Software Foundation, either version 3 of the License, or
8 // (at your option) any later version.
9 //
10 // DOLFIN is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17 //
18 // First added: 2014-02-24
19 // Last changed: 2017-12-06
20 
21 #ifndef __SIMPLEX_QUADRATURE_H
22 #define __SIMPLEX_QUADRATURE_H
23 
24 #include <vector>
25 #include <Eigen/Dense>
26 #include "Point.h"
27 
28 namespace dolfin
29 {
30 
31  // Forward declarations
32  class Cell;
33 
35 
37  {
38  public:
39 
48  SimplexQuadrature(std::size_t tdim, std::size_t order);
49 
60  std::pair<std::vector<double>, std::vector<double>>
61  compute_quadrature_rule(const Cell& cell) const;
62 
77  std::pair<std::vector<double>, std::vector<double>>
78  compute_quadrature_rule(const std::vector<Point>& coordinates,
79  std::size_t gdim) const;
80 
93  std::pair<std::vector<double>, std::vector<double>>
94  compute_quadrature_rule_interval(const std::vector<Point>& coordinates,
95  std::size_t gdim) const;
96 
109  std::pair<std::vector<double>, std::vector<double>>
110  compute_quadrature_rule_triangle(const std::vector<Point>& coordinates,
111  std::size_t gdim) const;
112 
125  std::pair<std::vector<double>, std::vector<double>>
126  compute_quadrature_rule_tetrahedron(const std::vector<Point>& coordinates,
127  std::size_t gdim) const;
128 
147  static std::vector<std::size_t>
148  compress(std::pair<std::vector<double>, std::vector<double>>& qr,
149  std::size_t gdim,
150  std::size_t quadrature_order);
151 
152  private:
153 
154  // Setup quadrature rule on a reference simplex
155  void setup_qr_reference_interval(std::size_t order);
156  void setup_qr_reference_triangle(std::size_t order);
157  void setup_qr_reference_tetrahedron(std::size_t order);
158 
159  // Utility function for computing a Vandermonde type matrix in a
160  // Chebyshev basis
161  static Eigen::MatrixXd
162  Chebyshev_Vandermonde_matrix
163  (const std::pair<std::vector<double>, std::vector<double>>& qr,
164  std::size_t gdim, std::size_t N);
165 
166  // Utility function for computing a Chebyshev basis
167  static std::vector<Eigen::VectorXd>
168  Chebyshev_polynomial(const Eigen::VectorXd& x, std::size_t N);
169 
170  // Utility function for creating a matrix with coefficients in
171  // graded lexicographic order
172  static std::vector<std::vector<std::size_t>>
173  grlex(std::size_t gdim, std::size_t N);
174 
175  // Utility function for calculating all combinations (n over k)
176  static std::size_t choose(std::size_t n, std::size_t k);
177 
178  // The following code has been copied from
179  //
180  // https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_dunavant_rule/triangle_dunavant_rule.cpp
181  //
182  // License: LGPL
183 
184  // Compute Duanvant quadrature rules for triangle
185 
186  static void dunavant_rule(std::size_t order,
187  std::vector<std::vector<double> >& p,
188  std::vector<double>& w);
189  static std::size_t dunavant_order_num(std::size_t rule);
190  static std::vector<std::size_t> dunavant_suborder(int rule, int suborder_num);
191  static std::size_t dunavant_suborder_num(int rule);
192  static void dunavant_subrule(std::size_t rule,
193  std::size_t suborder_num,
194  std::vector<double>& suborder_xyz,
195  std::vector<double>& w);
196  static void dunavant_subrule_01(int suborder_num,
197  std::vector<double>& suborder_xyz,
198  std::vector<double>& suborder_w);
199  static void dunavant_subrule_02(int suborder_num,
200  std::vector<double>& suborder_xyz,
201  std::vector<double>& suborder_w);
202  static void dunavant_subrule_03(int suborder_num,
203  std::vector<double>& suborder_xyz,
204  std::vector<double>& suborder_w);
205  static void dunavant_subrule_04(int suborder_num,
206  std::vector<double>& suborder_xyz,
207  std::vector<double>& suborder_w);
208  static void dunavant_subrule_05(int suborder_num,
209  std::vector<double>& suborder_xyz,
210  std::vector<double>& suborder_w);
211  static void dunavant_subrule_06(int suborder_num,
212  std::vector<double>& suborder_xyz,
213  std::vector<double>& suborder_w);
214  static void dunavant_subrule_07(int suborder_num,
215  std::vector<double>& suborder_xyz,
216  std::vector<double>& suborder_w);
217  static void dunavant_subrule_08(int suborder_num,
218  std::vector<double>& suborder_xyz,
219  std::vector<double>& suborder_w);
220  static void dunavant_subrule_09(int suborder_num,
221  std::vector<double>& suborder_xyz,
222  std::vector<double>& suborder_w);
223  static void dunavant_subrule_10(int suborder_num,
224  std::vector<double>& suborder_xyz,
225  std::vector<double>& suborder_w);
226  static void dunavant_subrule_11(int suborder_num,
227  std::vector<double>& suborder_xyz,
228  std::vector<double>& suborder_w);
229  static void dunavant_subrule_12(int suborder_num,
230  std::vector<double>& suborder_xyz,
231  std::vector<double>& suborder_w);
232  static void dunavant_subrule_13(int suborder_num,
233  std::vector<double>& suborder_xyz,
234  std::vector<double>& suborder_w);
235  static void dunavant_subrule_14(int suborder_num,
236  std::vector<double>& suborder_xyz,
237  std::vector<double>& suborder_w);
238  static void dunavant_subrule_15(int suborder_num,
239  std::vector<double>& suborder_xyz,
240  std::vector<double>& suborder_w);
241  static void dunavant_subrule_16(int suborder_num,
242  std::vector<double>& suborder_xyz,
243  std::vector<double>& suborder_w);
244  static void dunavant_subrule_17(int suborder_num,
245  std::vector<double>& suborder_xyz,
246  std::vector<double>& suborder_w);
247  static void dunavant_subrule_18(int suborder_num,
248  std::vector<double>& suborder_xyz,
249  std::vector<double>& suborder_w);
250  static void dunavant_subrule_19(int suborder_num,
251  std::vector<double>& suborder_xyz,
252  std::vector<double>& suborder_w);
253  static void dunavant_subrule_20(int suborder_num,
254  std::vector<double>& suborder_xyz,
255  std::vector<double>& suborder_w);
256  static int i4_modp(int i, int j);
257  static int i4_wrap(int ival, int ilo, int ihi);
258 
259  // The following code has been copied from
260  //
261  // https://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule_fast/legendre_rule_fast.cpp
262  //
263  // License: LGPL
264 
265  // Compute Gauss-Legendre quadrature rules for line
266 
267  static void legendre_compute_glr(std::size_t n,
268  std::vector<double>& x,
269  std::vector<double>& w);
270  static void legendre_compute_glr0(std::size_t n,
271  double& p,
272  double& pp);
273  static void legendre_compute_glr1(std::size_t n,
274  std::vector<double>& x,
275  std::vector<double>& w);
276  static void legendre_compute_glr2(double pn0, int n, double& x1, double& d1);
277  static double ts_mult(std::vector<double>& u, double h, int n);
278  static double rk2_leg(double t1, double t2, double x, int n);
279 
280  // Quadrature rule on reference simplex (points and weights)
281  std::vector<std::vector<double> > _p;
282  std::vector<double> _w;
283 
284  };
285 
286 }
287 
288 #endif
Definition: adapt.h:29
std::pair< std::vector< double >, std::vector< double > > compute_quadrature_rule(const Cell &cell) const
Definition: SimplexQuadrature.cpp:51
A Cell is a MeshEntity of topological codimension 0.
Definition: Cell.h:42
static std::vector< std::size_t > compress(std::pair< std::vector< double >, std::vector< double >> &qr, std::size_t gdim, std::size_t quadrature_order)
Definition: SimplexQuadrature.cpp:298
std::pair< std::vector< double >, std::vector< double > > compute_quadrature_rule_interval(const std::vector< Point > &coordinates, std::size_t gdim) const
Definition: SimplexQuadrature.cpp:102
std::pair< std::vector< double >, std::vector< double > > compute_quadrature_rule_tetrahedron(const std::vector< Point > &coordinates, std::size_t gdim) const
Definition: SimplexQuadrature.cpp:238
This class defines quadrature rules for simplices.
Definition: SimplexQuadrature.h:36
SimplexQuadrature(std::size_t tdim, std::size_t order)
Definition: SimplexQuadrature.cpp:28
std::pair< std::vector< double >, std::vector< double > > compute_quadrature_rule_triangle(const std::vector< Point > &coordinates, std::size_t gdim) const
Definition: SimplexQuadrature.cpp:170