libdolfin-dev 2017.2.0.post0-2 / usr / include / dolfin / geometry / SimplexQuadrature.h

This file is indexed.

/usr/include/dolfin/geometry/SimplexQuadrature.h is in libdolfin-dev 2017.2.0.post0-2.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
// Copyright (C) 2014 Anders Logg
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// DOLFIN is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// First added:  2014-02-24
// Last changed: 2017-09-22

#ifndef __SIMPLEX_QUADRATURE_H
#define __SIMPLEX_QUADRATURE_H

#include <vector>
#include <Eigen/Dense>
#include "Point.h"

namespace dolfin
{

  // Forward declarations
  class Cell;

  /// This class defines quadrature rules for simplices.

  class SimplexQuadrature
  {
  public:

    /// Create SimplexQuadrature rules for reference simplex
    ///
    /// *Arguments*
    ///     tdim (std::size_t)
    ///         The topological dimension of the simplex.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    SimplexQuadrature(std::size_t tdim, std::size_t order);

    /// Compute quadrature rule for cell.
    ///
    /// *Arguments*
    ///     cell (Cell)
    ///         The cell.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    /// *Returns*
    ///     std::pair<std::vector<double>, std::vector<double>>
    ///         A flattened array of quadrature points and a
    ///         corresponding array of quadrature weights.
    std::pair<std::vector<double>, std::vector<double>>
    compute_quadrature_rule(const Cell& cell,
			    std::size_t order) const;

    /// Compute quadrature rule for simplex.
    ///
    /// *Arguments*
    ///     coordinates (std::vector<Point>)
    ///         Vertex coordinates for the simplex
    ///     tdim (std::size_t)
    ///         The topological dimension of the simplex.
    ///     gdim (std::size_t)
    ///         The geometric dimension.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    /// *Returns*
    ///     std::pair<std::vector<double>, std::vector<double>>
    ///         A flattened array of quadrature points and a
    ///         corresponding array of quadrature weights.
    std::pair<std::vector<double>, std::vector<double>>
    compute_quadrature_rule(const std::vector<Point>& coordinates,
			    std::size_t gdim,
                            std::size_t order) const;

    /// Compute quadrature rule for interval.
    ///
    /// *Arguments*
    ///     coordinates (std::vector<Point>)
    ///         Vertex coordinates for the simplex
    ///     gdim (std::size_t)
    ///         The geometric dimension.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    /// *Returns*
    ///     std::pair<std::vector<double>, std::vector<double>>
    ///         A flattened array of quadrature points and a
    ///         corresponding array of quadrature weights.
    std::pair<std::vector<double>, std::vector<double>>
    compute_quadrature_rule_interval(const std::vector<Point>& coordinates,
				     std::size_t gdim,
                                     std::size_t order) const;

    /// Compute quadrature rule for triangle.
    ///
    /// *Arguments*
    ///     coordinates (std::vector<Point>)
    ///         Vertex coordinates for the simplex
    ///     gdim (std::size_t)
    ///         The geometric dimension.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    /// *Returns*
    ///     std::pair<std::vector<double>, std::vector<double>>
    ///         A flattened array of quadrature points and a
    ///         corresponding array of quadrature weights.
    std::pair<std::vector<double>, std::vector<double>>
    compute_quadrature_rule_triangle(const std::vector<Point>& coordinates,
                                     std::size_t gdim,
				     std::size_t order) const;
    
    /// Compute quadrature rule for tetrahedron.
    ///
    /// *Arguments*
    ///     coordinates (std::vector<Point>)
    ///         Vertex coordinates for the simplex
    ///     gdim (std::size_t)
    ///         The geometric dimension.
    ///     order (std::size_t)
    ///         The order of convergence of the quadrature rule.
    ///
    /// *Returns*
    ///     std::pair<std::vector<double>, std::vector<double>>
    ///         A flattened array of quadrature points and a
    ///         corresponding array of quadrature weights.
    std::pair<std::vector<double>, std::vector<double>>
    compute_quadrature_rule_tetrahedron(const std::vector<Point>& coordinates,
					std::size_t gdim,
					std::size_t order) const;
    
    /// Compress a quadrature rule using algorithms from
    ///     Compression of multivariate discrete measures and applications
    ///     A. Sommariva, M. Vianello
    ///     Numerical Functional Analysis and Optimization
    ///     Volume 36, 2015 - Issue 9
    ///
    /// *Arguments*
    ///     qr (std::pair<std::vector<double>, std::vector<double>>)
    ///         The quadrature rule to be compressed
    ///     gdim (std::size_t)
    ///         The geometric dimension
    ///     quadrature_order (std::size_t)
    ///         The order of the quadrature rule
    ///
    /// *Returns*
    ///     std::vector<std::size_t>
    ///         The indices of the points that were kept (empty
    ///         if no compression was made)
    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);

  private:

    // Setup quadrature rule on a reference simplex
    void setup_qr_reference_interval(std::size_t order);
    void setup_qr_reference_triangle(std::size_t order);
    void setup_qr_reference_tetrahedron(std::size_t order);

    // Utility function for computing a Vandermonde type matrix in a
    // Chebyshev basis
    static Eigen::MatrixXd
      Chebyshev_Vandermonde_matrix
      (const std::pair<std::vector<double>, std::vector<double>>& qr,
       std::size_t gdim, std::size_t N);
    
    // Utility function for computing a Chebyshev basis
    static std::vector<Eigen::VectorXd>
      Chebyshev_polynomial(const Eigen::VectorXd& x, std::size_t N);
    
    // Utility function for creating a matrix with coefficients in
    // graded lexicographic order
    static std::vector<std::vector<std::size_t>>
      grlex(std::size_t gdim, std::size_t N);
    
    // Utility function for calculating all combinations (n over k)
    static std::size_t choose(std::size_t n, std::size_t k);
    
    // The following code has been copied from
    //
    // https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_dunavant_rule/triangle_dunavant_rule.cpp
    //
    // License: LGPL

    // Compute Duanvant quadrature rules for triangle

    static void dunavant_rule(std::size_t order,
			      std::vector<std::vector<double> >& p,
			      std::vector<double>& w);
    static std::size_t dunavant_order_num(std::size_t rule);
    static std::vector<std::size_t> dunavant_suborder(int rule, int suborder_num);
    static std::size_t dunavant_suborder_num(int rule);
    static void dunavant_subrule(std::size_t rule,
				 std::size_t suborder_num,
				 std::vector<double>& suborder_xyz,
				 std::vector<double>& w);
    static void dunavant_subrule_01(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_02(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_03(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_04(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_05(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_06(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_07(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_08(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_09(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_10(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_11(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_12(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_13(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_14(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_15(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_16(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_17(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_18(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_19(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static void dunavant_subrule_20(int suborder_num,
				    std::vector<double>& suborder_xyz,
				    std::vector<double>& suborder_w);
    static int i4_modp(int i, int j);
    static int i4_wrap(int ival, int ilo, int ihi);

    // The following code has been copied from
    //
    // https://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule_fast/legendre_rule_fast.cpp
    //
    // License: LGPL

    // Compute Gauss-Legendre quadrature rules for line

    static void legendre_compute_glr(std::size_t n,
				     std::vector<double>& x,
				     std::vector<double>& w);
    static void legendre_compute_glr0(std::size_t n,
				      double& p,
				      double& pp);
    static void legendre_compute_glr1(std::size_t n,
				      std::vector<double>& x,
				      std::vector<double>& w);
    static void legendre_compute_glr2(double pn0, int n, double& x1, double& d1);
    static double ts_mult(std::vector<double>& u, double h, int n);
    static double rk2_leg(double t1, double t2, double x, int n);

    // Quadrature rule on reference simplex (points and weights)
    std::vector<std::vector<double> > _p;
    std::vector<double> _w;

  };

}

#endif

APT Browse - Built by Thomas Orozco.