DOLFIN
DOLFIN C++ interface
Loading...
Searching...
No Matches
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
28namespace 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
A Cell is a MeshEntity of topological codimension 0.
Definition Cell.h:43
This class defines quadrature rules for simplices.
Definition SimplexQuadrature.h:37
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
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(const Cell &cell) const
Definition SimplexQuadrature.cpp:51
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
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
Definition adapt.h:30