Loading...
Searching...
No Matches
SimplexQuadrature.h
// 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-12-06
#ifndef __SIMPLEX_QUADRATURE_H
#define __SIMPLEX_QUADRATURE_H
#include <vector>
#include <Eigen/Dense>
#include "Point.h"
namespace dolfin
{
// Forward declarations
class Cell;
class SimplexQuadrature
{
public:
SimplexQuadrature(std::size_t tdim, std::size_t order);
std::pair<std::vector<double>, std::vector<double>>
compute_quadrature_rule(const Cell& cell) const;
std::pair<std::vector<double>, std::vector<double>>
compute_quadrature_rule(const std::vector<Point>& coordinates,
std::size_t gdim) const;
std::pair<std::vector<double>, std::vector<double>>
compute_quadrature_rule_interval(const std::vector<Point>& coordinates,
std::size_t gdim) const;
std::pair<std::vector<double>, std::vector<double>>
compute_quadrature_rule_triangle(const std::vector<Point>& coordinates,
std::size_t gdim) const;
std::pair<std::vector<double>, std::vector<double>>
compute_quadrature_rule_tetrahedron(const std::vector<Point>& coordinates,
std::size_t gdim) const;
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
dolfin::Cell
A Cell is a MeshEntity of topological codimension 0.
Definition Cell.h:43
dolfin::SimplexQuadrature::SimplexQuadrature
SimplexQuadrature(std::size_t tdim, std::size_t order)
Definition SimplexQuadrature.cpp:28
dolfin::SimplexQuadrature::compute_quadrature_rule_triangle
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
dolfin::SimplexQuadrature::compress
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
dolfin::SimplexQuadrature::compute_quadrature_rule
std::pair< std::vector< double >, std::vector< double > > compute_quadrature_rule(const Cell &cell) const
Definition SimplexQuadrature.cpp:51
dolfin::SimplexQuadrature::compute_quadrature_rule_tetrahedron
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
dolfin::SimplexQuadrature::compute_quadrature_rule_interval
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
dolfin
Definition adapt.h:30