2974 lines
118 KiB
C++
2974 lines
118 KiB
C++
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
|
||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||
//
|
||
// This file is part of the MFEM library. For more information and source code
|
||
// availability visit https://mfem.org.
|
||
//
|
||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||
// CONTRIBUTING.md for details.
|
||
|
||
// Implementation of IntegrationRule(s) classes
|
||
|
||
// Acknowledgment: Some of the high-precision triangular and tetrahedral
|
||
// quadrature rules below were obtained from the Encyclopaedia of Cubature
|
||
// Formulas at http://nines.cs.kuleuven.be/research/ecf/ecf.html
|
||
//
|
||
// Positive-weight simplex quadrature rules from:
|
||
//
|
||
// [1] F.D. Witherden, P.E. Vincent, "On the identification of symmetric
|
||
// quadrature rules for finite element methods", Computers & Mathematics
|
||
// with Applications, 69(10):1232-1241, 2015.
|
||
// Used for: triangles (d=0-20), tetrahedra (d=0-13).
|
||
//
|
||
// [2] G. Chuluunbaatar, O. Chuluunbaatar, A.A. Gusev, S.I. Vinitsky,
|
||
// "PI-type fully symmetric quadrature rules on the 3-,...,6-simplexes",
|
||
// Computers & Mathematics with Applications, 124:89-97, 2022.
|
||
// Used for: tetrahedra (d=14-20).
|
||
|
||
#include "fem.hpp"
|
||
#include "../mesh/nurbs.hpp"
|
||
#include <cmath>
|
||
|
||
#ifdef MFEM_USE_MPFR
|
||
#include <mpfr.h>
|
||
#endif
|
||
|
||
using namespace std;
|
||
|
||
namespace mfem
|
||
{
|
||
|
||
IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry)
|
||
{
|
||
int i, j, nx, ny;
|
||
|
||
nx = irx.GetNPoints();
|
||
ny = iry.GetNPoints();
|
||
SetSize(nx * ny);
|
||
SetPointIndices();
|
||
Order = std::min(irx.GetOrder(), iry.GetOrder());
|
||
|
||
for (j = 0; j < ny; j++)
|
||
{
|
||
IntegrationPoint &ipy = iry.IntPoint(j);
|
||
for (i = 0; i < nx; i++)
|
||
{
|
||
IntegrationPoint &ipx = irx.IntPoint(i);
|
||
IntegrationPoint &ip = IntPoint(j*nx+i);
|
||
|
||
ip.x = ipx.x;
|
||
ip.y = ipy.x;
|
||
ip.weight = ipx.weight * ipy.weight;
|
||
}
|
||
}
|
||
}
|
||
|
||
IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
|
||
IntegrationRule &irz)
|
||
{
|
||
const int nx = irx.GetNPoints();
|
||
const int ny = iry.GetNPoints();
|
||
const int nz = irz.GetNPoints();
|
||
SetSize(nx*ny*nz);
|
||
SetPointIndices();
|
||
Order = std::min({irx.GetOrder(), iry.GetOrder(), irz.GetOrder()});
|
||
|
||
for (int iz = 0; iz < nz; ++iz)
|
||
{
|
||
IntegrationPoint &ipz = irz.IntPoint(iz);
|
||
for (int iy = 0; iy < ny; ++iy)
|
||
{
|
||
IntegrationPoint &ipy = iry.IntPoint(iy);
|
||
for (int ix = 0; ix < nx; ++ix)
|
||
{
|
||
IntegrationPoint &ipx = irx.IntPoint(ix);
|
||
IntegrationPoint &ip = IntPoint(iz*nx*ny + iy*nx + ix);
|
||
|
||
ip.x = ipx.x;
|
||
ip.y = ipy.x;
|
||
ip.z = ipz.x;
|
||
ip.weight = ipx.weight*ipy.weight*ipz.weight;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
const Array<real_t> &IntegrationRule::GetWeights() const
|
||
{
|
||
if (weights.Size() != GetNPoints())
|
||
{
|
||
weights.SetSize(GetNPoints());
|
||
for (int i = 0; i < GetNPoints(); i++)
|
||
{
|
||
weights[i] = IntPoint(i).weight;
|
||
}
|
||
}
|
||
return weights;
|
||
}
|
||
|
||
void IntegrationRule::SetPointIndices()
|
||
{
|
||
for (int i = 0; i < Size(); i++)
|
||
{
|
||
IntPoint(i).index = i;
|
||
}
|
||
}
|
||
|
||
void IntegrationRule::GrundmannMollerSimplexRule(int s, int n)
|
||
{
|
||
// for pow on older compilers
|
||
using std::pow;
|
||
const int d = 2*s + 1;
|
||
Vector fact(d + n + 1);
|
||
Array<int> beta(n), sums(n);
|
||
|
||
fact(0) = 1.;
|
||
for (int i = 1; i < fact.Size(); i++)
|
||
{
|
||
fact(i) = fact(i - 1)*i;
|
||
}
|
||
|
||
// number of points is \binom{n + s + 1}{n + 1}
|
||
int np = 1, f = 1;
|
||
for (int i = 0; i <= n; i++)
|
||
{
|
||
np *= (s + i + 1), f *= (i + 1);
|
||
}
|
||
np /= f;
|
||
SetSize(np);
|
||
SetPointIndices();
|
||
Order = 2*s + 1;
|
||
|
||
int pt = 0;
|
||
for (int i = 0; i <= s; i++)
|
||
{
|
||
real_t weight;
|
||
|
||
weight = pow(2., -2*s)*pow(static_cast<real_t>(d + n - 2*i),
|
||
d)/fact(i)/fact(d + n - i);
|
||
if (i%2)
|
||
{
|
||
weight = -weight;
|
||
}
|
||
|
||
// loop over all beta : beta_0 + ... + beta_{n-1} <= s - i
|
||
int k = s - i;
|
||
beta = 0;
|
||
sums = 0;
|
||
while (true)
|
||
{
|
||
IntegrationPoint &ip = IntPoint(pt++);
|
||
ip.weight = weight;
|
||
ip.x = real_t(2*beta[0] + 1)/(d + n - 2*i);
|
||
ip.y = real_t(2*beta[1] + 1)/(d + n - 2*i);
|
||
if (n == 3)
|
||
{
|
||
ip.z = real_t(2*beta[2] + 1)/(d + n - 2*i);
|
||
}
|
||
|
||
int j = 0;
|
||
while (sums[j] == k)
|
||
{
|
||
beta[j++] = 0;
|
||
if (j == n)
|
||
{
|
||
goto done_beta;
|
||
}
|
||
}
|
||
beta[j]++;
|
||
sums[j]++;
|
||
for (j--; j >= 0; j--)
|
||
{
|
||
sums[j] = sums[j+1];
|
||
}
|
||
}
|
||
done_beta:
|
||
;
|
||
}
|
||
}
|
||
|
||
IntegrationRule*
|
||
IntegrationRule::ApplyToKnotIntervals(KnotVector const& kv) const
|
||
{
|
||
const int np = this->GetNPoints();
|
||
const int ne = kv.GetNE();
|
||
|
||
IntegrationRule *kvir = new IntegrationRule(ne * np);
|
||
kvir->SetOrder(GetOrder());
|
||
|
||
real_t x0 = kv[0];
|
||
real_t x1 = x0;
|
||
|
||
int id = 0;
|
||
for (int e=0; e<ne; ++e)
|
||
{
|
||
x0 = x1;
|
||
|
||
if (e == ne-1)
|
||
{
|
||
x1 = kv[kv.Size() - 1];
|
||
}
|
||
else
|
||
{
|
||
// Find the next unique knot
|
||
while (id < kv.Size() - 1)
|
||
{
|
||
id++;
|
||
if (kv[id] != x0)
|
||
{
|
||
x1 = kv[id];
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
|
||
const real_t s = x1 - x0;
|
||
|
||
for (int j=0; j<this->GetNPoints(); ++j)
|
||
{
|
||
const real_t x = x0 + (s * (*this)[j].x);
|
||
(*kvir)[(e * np) + j].Set1w(x, (*this)[j].weight);
|
||
}
|
||
}
|
||
|
||
return kvir;
|
||
}
|
||
|
||
IntegrationRule IntegrationRule::Reorder(const Array<int> &ordering) const
|
||
{
|
||
const int np = GetNPoints();
|
||
MFEM_VERIFY(np == ordering.Size(), "Invalid permutation size");
|
||
IntegrationRule ir(np);
|
||
ir.SetOrder(GetOrder());
|
||
|
||
for (int i = 0; i < np; i++)
|
||
{
|
||
IntegrationPoint &ip_new = ir.IntPoint(i);
|
||
const IntegrationPoint &ip_old = IntPoint(ordering[i]);
|
||
ip_new.Set(ip_old.x, ip_old.y, ip_old.z, ip_old.weight);
|
||
}
|
||
|
||
return ir;
|
||
}
|
||
|
||
IntegrationRule DuffyTrans(const IntegrationRule &ir, int dim)
|
||
{
|
||
IntegrationRule ir_mapped(ir.GetNPoints());
|
||
ir_mapped.SetOrder(ir.GetOrder());
|
||
|
||
if (dim == 2)
|
||
{
|
||
for (int i = 0; i < ir.GetNPoints(); i++)
|
||
{
|
||
IntegrationPoint &ip_mapped = ir_mapped.IntPoint(i);
|
||
ip_mapped.y = ir.IntPoint(i).y * (1 - ir.IntPoint(i).x);
|
||
ip_mapped.x = ir.IntPoint(i).x;
|
||
ip_mapped.weight = ir.IntPoint(i).weight;
|
||
}
|
||
return ir_mapped;
|
||
}
|
||
else if (dim == 3)
|
||
{
|
||
for (int i = 0; i < ir.GetNPoints(); i++)
|
||
{
|
||
IntegrationPoint &ip_mapped = ir_mapped.IntPoint(i);
|
||
ip_mapped.z = ir.IntPoint(i).z * (1 - ir.IntPoint(i).x) * (1 - ir.IntPoint(
|
||
i).y);
|
||
ip_mapped.y = ir.IntPoint(i).y * (1 - ir.IntPoint(i).x);
|
||
ip_mapped.x = ir.IntPoint(i).x;
|
||
ip_mapped.weight = ir.IntPoint(i).weight;
|
||
}
|
||
return ir_mapped;
|
||
}
|
||
else
|
||
{
|
||
MFEM_ABORT("Duffy transformation not implemented for this dimension!");
|
||
}
|
||
}
|
||
|
||
#ifdef MFEM_USE_MPFR
|
||
|
||
// Class for computing hi-precision (HP) quadrature in 1D
|
||
class HP_Quadrature1D
|
||
{
|
||
protected:
|
||
mpfr_t pi, z, pp, p1, p2, p3, dz, w, rtol;
|
||
|
||
public:
|
||
static const mpfr_rnd_t rnd = GMP_RNDN;
|
||
static const int default_prec = 128;
|
||
|
||
// prec = MPFR precision in bits
|
||
HP_Quadrature1D(const int prec = default_prec)
|
||
{
|
||
mpfr_inits2(prec, pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
|
||
mpfr_const_pi(pi, rnd);
|
||
mpfr_set_si_2exp(rtol, 1, -32, rnd); // 2^(-32) < 2.33e-10
|
||
}
|
||
|
||
// set rtol = 2^exponent
|
||
// this is a tolerance for the last correction of x_i in Newton's algorithm;
|
||
// this gives roughly rtol^2 accuracy for the final x_i.
|
||
void SetRelTol(const int exponent = -32)
|
||
{
|
||
mpfr_set_si_2exp(rtol, 1, exponent, rnd);
|
||
}
|
||
|
||
// n - number of quadrature points
|
||
// k - index of the point to compute, 0 <= k < n
|
||
// see also: QuadratureFunctions1D::GaussLegendre
|
||
void ComputeGaussLegendrePoint(const int n, const int k)
|
||
{
|
||
MFEM_ASSERT(n > 0 && 0 <= k && k < n, "invalid n = " << n
|
||
<< " and/or k = " << k);
|
||
|
||
int i = (k < (n+1)/2) ? k+1 : n-k;
|
||
|
||
// Initial guess for the x-coordinate:
|
||
// set z = cos(pi * (i - 0.25) / (n + 0.5)) =
|
||
// = sin(pi * ((n+1-2*i) / (2*n+1)))
|
||
mpfr_set_si(z, n+1-2*i, rnd);
|
||
mpfr_div_si(z, z, 2*n+1, rnd);
|
||
mpfr_mul(z, z, pi, rnd);
|
||
mpfr_sin(z, z, rnd);
|
||
|
||
bool done = false;
|
||
while (1)
|
||
{
|
||
mpfr_set_si(p2, 1, rnd);
|
||
mpfr_set(p1, z, rnd);
|
||
for (int j = 2; j <= n; j++)
|
||
{
|
||
mpfr_set(p3, p2, rnd);
|
||
mpfr_set(p2, p1, rnd);
|
||
// p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
|
||
mpfr_mul_si(p1, z, 2*j-1, rnd);
|
||
mpfr_mul_si(p3, p3, j-1, rnd);
|
||
mpfr_fms(p1, p1, p2, p3, rnd);
|
||
mpfr_div_si(p1, p1, j, rnd);
|
||
}
|
||
// p1 is Legendre polynomial
|
||
|
||
// derivative of the Legendre polynomial:
|
||
// pp = n * (z*p1-p2) / (z*z - 1);
|
||
mpfr_fms(pp, z, p1, p2, rnd);
|
||
mpfr_mul_si(pp, pp, n, rnd);
|
||
mpfr_sqr(p2, z, rnd);
|
||
mpfr_sub_si(p2, p2, 1, rnd);
|
||
mpfr_div(pp, pp, p2, rnd);
|
||
|
||
if (done) { break; }
|
||
|
||
// set delta_z: dz = p1/pp;
|
||
mpfr_div(dz, p1, pp, rnd);
|
||
// compute absolute tolerance: atol = rtol*(1-z)
|
||
mpfr_t &atol = w;
|
||
mpfr_si_sub(atol, 1, z, rnd);
|
||
mpfr_mul(atol, atol, rtol, rnd);
|
||
if (mpfr_cmpabs(dz, atol) <= 0)
|
||
{
|
||
done = true;
|
||
// continue the computation: get pp at the new point, then exit
|
||
}
|
||
// update z = z - dz
|
||
mpfr_sub(z, z, dz, rnd);
|
||
}
|
||
|
||
// map z to (0,1): z = (1 - z)/2
|
||
mpfr_si_sub(z, 1, z, rnd);
|
||
mpfr_div_2si(z, z, 1, rnd);
|
||
|
||
// weight: w = 1/(4*z*(1 - z)*pp*pp)
|
||
mpfr_sqr(w, pp, rnd);
|
||
mpfr_mul_2si(w, w, 2, rnd);
|
||
mpfr_mul(w, w, z, rnd);
|
||
mpfr_si_sub(p1, 1, z, rnd); // p1 = 1-z
|
||
mpfr_mul(w, w, p1, rnd);
|
||
mpfr_si_div(w, 1, w, rnd);
|
||
|
||
if (k >= (n+1)/2) { mpfr_swap(z, p1); }
|
||
}
|
||
|
||
// n - number of quadrature points
|
||
// k - index of the point to compute, 0 <= k < n
|
||
// see also: QuadratureFunctions1D::GaussLobatto
|
||
void ComputeGaussLobattoPoint(const int n, const int k)
|
||
{
|
||
MFEM_ASSERT(n > 1 && 0 <= k && k < n, "invalid n = " << n
|
||
<< " and/or k = " << k);
|
||
|
||
int i = (k < (n+1)/2) ? k : n-1-k;
|
||
|
||
if (i == 0)
|
||
{
|
||
mpfr_set_si(z, 0, rnd);
|
||
mpfr_set_si(p1, 1, rnd);
|
||
mpfr_set_si(w, n*(n-1), rnd);
|
||
mpfr_si_div(w, 1, w, rnd); // weight = 1/(n*(n-1))
|
||
return;
|
||
}
|
||
// initial guess is the corresponding Chebyshev point, z:
|
||
// z = -cos(pi * i/(n-1)) = sin(pi * (2*i-n+1)/(2*n-2))
|
||
mpfr_set_si(z, 2*i-n+1, rnd);
|
||
mpfr_div_si(z, z, 2*(n-1), rnd);
|
||
mpfr_mul(z, pi, z, rnd);
|
||
mpfr_sin(z, z, rnd);
|
||
bool done = false;
|
||
for (int iter = 0 ; true ; ++iter)
|
||
{
|
||
// build Legendre polynomials, up to P_{n}(z)
|
||
mpfr_set_si(p1, 1, rnd);
|
||
mpfr_set(p2, z, rnd);
|
||
|
||
for (int l = 1 ; l < (n-1) ; ++l)
|
||
{
|
||
// P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1}(x) ]/(l+1)
|
||
mpfr_mul_si(p1, p1, l, rnd);
|
||
mpfr_mul_si(p3, z, 2*l+1, rnd);
|
||
mpfr_fms(p3, p3, p2, p1, rnd);
|
||
mpfr_div_si(p3, p3, l+1, rnd);
|
||
|
||
mpfr_set(p1, p2, rnd);
|
||
mpfr_set(p2, p3, rnd);
|
||
}
|
||
if (done) { break; }
|
||
// compute dz = resid/deriv = (z*p2 - p1) / (n*p2);
|
||
mpfr_fms(dz, z, p2, p1, rnd);
|
||
mpfr_mul_si(p3, p2, n, rnd);
|
||
mpfr_div(dz, dz, p3, rnd);
|
||
// update: z = z - dz
|
||
mpfr_sub(z, z, dz, rnd);
|
||
// compute absolute tolerance: atol = rtol*(1 + z)
|
||
mpfr_t &atol = w;
|
||
mpfr_add_si(atol, z, 1, rnd);
|
||
mpfr_mul(atol, atol, rtol, rnd);
|
||
// check for convergence
|
||
if (mpfr_cmpabs(dz, atol) <= 0)
|
||
{
|
||
done = true;
|
||
// continue the computation: get p2 at the new point, then exit
|
||
}
|
||
// If the iteration does not converge fast, something is wrong.
|
||
MFEM_VERIFY(iter < 8, "n = " << n << ", i = " << i
|
||
<< ", dz = " << mpfr_get_d(dz, rnd));
|
||
}
|
||
// Map to the interval [0,1] and scale the weights
|
||
mpfr_add_si(z, z, 1, rnd);
|
||
mpfr_div_2si(z, z, 1, rnd);
|
||
// set the symmetric point
|
||
mpfr_si_sub(p1, 1, z, rnd);
|
||
// w = 1/[ n*(n-1)*[P_{n-1}(z)]^2 ]
|
||
mpfr_sqr(w, p2, rnd);
|
||
mpfr_mul_si(w, w, n*(n-1), rnd);
|
||
mpfr_si_div(w, 1, w, rnd);
|
||
|
||
if (k >= (n+1)/2) { mpfr_swap(z, p1); }
|
||
}
|
||
|
||
real_t GetPoint() const { return mpfr_get_d(z, rnd); }
|
||
real_t GetSymmPoint() const { return mpfr_get_d(p1, rnd); }
|
||
real_t GetWeight() const { return mpfr_get_d(w, rnd); }
|
||
|
||
const mpfr_t &GetHPPoint() const { return z; }
|
||
const mpfr_t &GetHPSymmPoint() const { return p1; }
|
||
const mpfr_t &GetHPWeight() const { return w; }
|
||
|
||
~HP_Quadrature1D()
|
||
{
|
||
mpfr_clears(pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
|
||
mpfr_free_cache();
|
||
}
|
||
};
|
||
|
||
#endif // MFEM_USE_MPFR
|
||
|
||
|
||
void QuadratureFunctions1D::GaussJacobi(const int np, const real_t alpha,
|
||
const real_t beta, IntegrationRule* ir)
|
||
{
|
||
/* The np-point Gauss-Jacobi quadrature rule is exact for polynomials of
|
||
degree 2np - 1 with weight function w(x) = (1-x)^alpha * x^beta. The
|
||
nodes are the zeros of the Jacobi polynomial P_{np}^{alpha,beta} and
|
||
the weights are
|
||
|
||
w_i = C / [(1 - x_i^2) * P'_{np}^{alpha,beta}(x_i)^2]
|
||
C = 2^{alpha + beta + 1} * Gamma(np + alpha + 1) * Gamma(np + beta + 1)
|
||
/ [Gamma(np + alpha + beta + 1) * Gamma(np + 1)].
|
||
|
||
The nodes are computed via nonlinear solve (Newton's method) with an
|
||
initial guess corresponding to Gatteschi's asymptotic expansions of the
|
||
Jacobi polynomial roots [1].
|
||
|
||
The current initial guess has been tested and performs well for
|
||
np <= 200 and -1 <= alpha, beta <= 4. For larger np, it may be necessary
|
||
utilize different initial guesses in the vicinity of x = -1,+1 [2].
|
||
|
||
[1] Gautschi, W., & Giordano, C. (2008). Luigi Gatteschi’s work on
|
||
asymptotics of special functions and their zeros. Numerical Algorithms,
|
||
49, 11-31.
|
||
[2] Hale, N., & Townsend, A. (2013). Fast and accurate computation of
|
||
Gauss--Legendre and Gauss--Jacobi quadrature nodes and weights.
|
||
SIAM Journal on Scientific Computing, 35(2), A652-A674.
|
||
*/
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->SetOrder(2*np - 1);
|
||
|
||
if (alpha <= -1.0 || beta <= -1.0)
|
||
{
|
||
MFEM_ABORT("Gauss-Jacobi quadrature only defined for alpha > -1 and beta > -1");
|
||
}
|
||
// Jacobi weight function is undefined whenever alpha <= -1 or beta <= -1
|
||
|
||
if (alpha > 4.0 || beta > 4.0)
|
||
{
|
||
MFEM_ABORT("Current Gauss-Jacobi quadrature implementation only tested for alpha <= 4 and beta <= 4");
|
||
}
|
||
// current asymptotic expansions for initial guess may perform poorly for large alpha, beta
|
||
|
||
switch (np)
|
||
{
|
||
case 1:
|
||
real_t x = (beta - alpha) / (alpha + beta + 2);
|
||
real_t w = pow(2, alpha + beta + 1) * tgamma(alpha + 2) * tgamma(
|
||
beta + 2) / (tgamma(alpha + beta + 2));
|
||
w = 0.5 * w / pow(2, alpha + beta);
|
||
// map weight to to [0,1], with additional 1/(2^(alpha + beta)) factor coming from mapping
|
||
// the weight (1-x)^alpha * (1+x)^beta to [0,1] as well.
|
||
ir->IntPoint(0).Set1w(0.5 * x + 0.5,
|
||
4.0 * w / ((1.0 - x*x) * (alpha + beta + 2) * (alpha + beta + 2)));
|
||
return;
|
||
}
|
||
|
||
#ifdef MFEM_USE_MPFR
|
||
MFEM_WARNING("MPFR implementation of Gauss-Jacobi quadrature not implemented yet. Falling "
|
||
"back to double precision implementation...");
|
||
#endif
|
||
|
||
const int n = np;
|
||
// common constants for Jacobi polynomials
|
||
real_t ab = alpha + beta;
|
||
real_t a2_minus_b2 = (alpha - beta) * (alpha + beta);
|
||
|
||
// roots of P^(alpha,beta)_n in the interval [-1,1]
|
||
for (int i = 1; i <= n; i++)
|
||
{
|
||
// rather than using Chebyshev points for initial guess, use Gatteschi's asymptotic expansion for roots of Jacobi
|
||
// polynomials
|
||
real_t n_ab_plus_1 = 2 * n + alpha + beta + 1;
|
||
real_t v = (2 * i + alpha - 0.5) * M_PI / n_ab_plus_1;
|
||
real_t theta = v + 1.0 / (n_ab_plus_1*n_ab_plus_1) * ((0.25 - alpha*alpha) *
|
||
1.0/tan(0.5*v) - (0.25 - beta*beta) * tan(0.5*v));
|
||
real_t z = cos(theta);
|
||
|
||
real_t pp, p1, dz, xi = 0.;
|
||
bool done = false;
|
||
while (1)
|
||
{
|
||
real_t p2 = 1;
|
||
p1 = ((alpha-beta) + (alpha + beta + 2) * z) / 2;
|
||
for (int j = 1; j <= n-1; j++)
|
||
{
|
||
real_t p3 = p2;
|
||
p2 = p1;
|
||
|
||
real_t jx2_ab = 2 * j + ab;
|
||
real_t an = (jx2_ab) * (jx2_ab + 2);
|
||
real_t bn = a2_minus_b2;
|
||
real_t cn = 2 * (j + alpha) * (j + beta) * (jx2_ab + 2) / (jx2_ab + 1);
|
||
|
||
real_t D = (jx2_ab + 1) / (2 * (j + 1) * (j + ab + 1) * (jx2_ab));
|
||
p1 = ((an * z + bn) * p2 - cn * p3) * D;
|
||
}
|
||
// p1 is Jacobi polynomial
|
||
pp = n * (alpha - beta - (2 * n + ab) * z) * p1 + 2 * (n + alpha) *
|
||
(n + beta) * p2;
|
||
pp = pp / ((2 * n + ab) * (1 - z*z));
|
||
// derivative of the Jacobi polynomial
|
||
if (done) { break; }
|
||
|
||
dz = p1/pp;
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (std::abs(dz) < 1e-7)
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (std::abs(dz) < std::numeric_limits<real_t>::epsilon())
|
||
// this seems to cause trouble if we try std::abs(dz) < 1e-16
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
// if (std::abs(dz) < 1e-16)
|
||
#endif
|
||
{
|
||
done = true;
|
||
xi = z - dz;
|
||
}
|
||
z -= dz;
|
||
}
|
||
real_t c0 = exp(lgamma(n + alpha + 1) - lgamma(n + ab + 1)) * exp(lgamma(
|
||
n + beta + 1) - lgamma(n + 1));
|
||
// ratio of gamma functions prone to overflow for large n, so compute logarithms
|
||
// of Gamma function instead, i.e. Gamma(a)/Gamma(b) = exp(lgamma(a) - lgamma(b))
|
||
ir->IntPoint(n-i).x = 0.5 * xi + 0.5;
|
||
ir->IntPoint(n-i).weight = 0.5 * c0 * pow(2.0,
|
||
ab + 1) / ((1.0 - xi*xi)*pp*pp) / pow(2, ab);
|
||
// map nodes and weights to the interval [0,1]
|
||
}
|
||
}
|
||
|
||
|
||
void QuadratureFunctions1D::GaussLegendre(const int np, IntegrationRule* ir)
|
||
{
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->SetOrder(2*np - 1);
|
||
|
||
switch (np)
|
||
{
|
||
case 1:
|
||
ir->IntPoint(0).Set1w(0.5, 1.0);
|
||
return;
|
||
case 2:
|
||
ir->IntPoint(0).Set1w(0.21132486540518711775, 0.5);
|
||
ir->IntPoint(1).Set1w(0.78867513459481288225, 0.5);
|
||
return;
|
||
case 3:
|
||
ir->IntPoint(0).Set1w(0.11270166537925831148, 5./18.);
|
||
ir->IntPoint(1).Set1w(0.5, 4./9.);
|
||
ir->IntPoint(2).Set1w(0.88729833462074168852, 5./18.);
|
||
return;
|
||
}
|
||
|
||
const int n = np;
|
||
const int m = (n+1)/2;
|
||
|
||
#ifndef MFEM_USE_MPFR
|
||
|
||
for (int i = 1; i <= m; i++)
|
||
{
|
||
real_t z = cos(M_PI * (i - 0.25) / (n + 0.5));
|
||
real_t pp, p1, dz, xi = 0.;
|
||
bool done = false;
|
||
while (1)
|
||
{
|
||
real_t p2 = 1;
|
||
p1 = z;
|
||
for (int j = 2; j <= n; j++)
|
||
{
|
||
real_t p3 = p2;
|
||
p2 = p1;
|
||
p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
|
||
}
|
||
// p1 is Legendre polynomial
|
||
|
||
pp = n * (z*p1-p2) / (z*z - 1);
|
||
if (done) { break; }
|
||
|
||
dz = p1/pp;
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (std::abs(dz) < 1e-7)
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (std::abs(dz) < 1e-16)
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
if (std::abs(dz) < 1e-16)
|
||
#endif
|
||
{
|
||
done = true;
|
||
// map the new point (z-dz) to (0,1):
|
||
xi = ((1 - z) + dz)/2; // (1 - (z - dz))/2 has bad round-off
|
||
// continue the computation: get pp at the new point, then exit
|
||
}
|
||
// update: z = z - dz
|
||
z -= dz;
|
||
}
|
||
|
||
ir->IntPoint(i-1).x = xi;
|
||
ir->IntPoint(n-i).x = 1 - xi;
|
||
ir->IntPoint(i-1).weight =
|
||
ir->IntPoint(n-i).weight = 1./(4*xi*(1 - xi)*pp*pp);
|
||
}
|
||
|
||
#else // MFEM_USE_MPFR is defined
|
||
|
||
HP_Quadrature1D hp_quad;
|
||
for (int i = 1; i <= m; i++)
|
||
{
|
||
hp_quad.ComputeGaussLegendrePoint(n, i-1);
|
||
|
||
ir->IntPoint(i-1).x = hp_quad.GetPoint();
|
||
ir->IntPoint(n-i).x = hp_quad.GetSymmPoint();
|
||
ir->IntPoint(i-1).weight = ir->IntPoint(n-i).weight = hp_quad.GetWeight();
|
||
}
|
||
|
||
#endif // MFEM_USE_MPFR
|
||
|
||
}
|
||
|
||
void QuadratureFunctions1D::GaussLobatto(const int np, IntegrationRule* ir)
|
||
{
|
||
/* An np point Gauss-Lobatto quadrature has (np - 2) free abscissa the other
|
||
(2) abscissa are the interval endpoints.
|
||
|
||
The interior x_i are the zeros of P'_{np-1}(x). The weights of the
|
||
interior points on the interval [-1,1] are:
|
||
|
||
w_i = 2/(np*(np-1)*[P_{np-1}(x_i)]^2)
|
||
|
||
The end point weights (on [-1,1]) are: w_{end} = 2/(np*(np-1)).
|
||
|
||
The interior abscissa are found via a nonlinear solve, the initial guess
|
||
for each point is the corresponding Chebyshev point.
|
||
|
||
After we find all points on the interval [-1,1], we will map and scale the
|
||
points and weights to the MFEM natural interval [0,1].
|
||
|
||
References:
|
||
[1] E. E. Lewis and W. F. Millier, "Computational Methods of Neutron
|
||
Transport", Appendix A
|
||
[2] the QUADRULE software by John Burkardt,
|
||
https://people.sc.fsu.edu/~jburkardt/cpp_src/quadrule/quadrule.cpp
|
||
*/
|
||
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
if ( np == 1 )
|
||
{
|
||
ir->IntPoint(0).Set1w(0.5, 1.0);
|
||
ir->SetOrder(1);
|
||
}
|
||
else
|
||
{
|
||
ir->SetOrder(2*np - 3);
|
||
|
||
#ifndef MFEM_USE_MPFR
|
||
|
||
// endpoints and respective weights
|
||
ir->IntPoint(0).x = 0.0;
|
||
ir->IntPoint(np-1).x = 1.0;
|
||
ir->IntPoint(0).weight = ir->IntPoint(np-1).weight = 1.0/(np*(np-1));
|
||
|
||
// interior points and weights
|
||
// use symmetry and compute just half of the points
|
||
for (int i = 1 ; i <= (np-1)/2 ; ++i)
|
||
{
|
||
// initial guess is the corresponding Chebyshev point, x_i:
|
||
// x_i = -cos(\pi * (i / (np-1)))
|
||
real_t x_i = std::sin(M_PI * ((real_t)(i)/(np-1) - 0.5));
|
||
real_t z_i = 0., p_l;
|
||
bool done = false;
|
||
for (int iter = 0 ; true ; ++iter)
|
||
{
|
||
// build Legendre polynomials, up to P_{np}(x_i)
|
||
real_t p_lm1 = 1.0;
|
||
p_l = x_i;
|
||
|
||
for (int l = 1 ; l < (np-1) ; ++l)
|
||
{
|
||
// The Legendre polynomials can be built by recursion:
|
||
// x * P_l(x) = 1/(2*l+1)*[ (l+1)*P_{l+1}(x) + l*P_{l-1} ], i.e.
|
||
// P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1} ]/(l+1)
|
||
real_t p_lp1 = ( (2*l + 1)*x_i*p_l - l*p_lm1)/(l + 1);
|
||
|
||
p_lm1 = p_l;
|
||
p_l = p_lp1;
|
||
}
|
||
if (done) { break; }
|
||
// after this loop, p_l holds P_{np-1}(x_i)
|
||
// resid = (x^2-1)*P'_{np-1}(x_i)
|
||
// but use the recurrence relationship
|
||
// (x^2 -1)P'_l(x) = l*[ x*P_l(x) - P_{l-1}(x) ]
|
||
// thus, resid = (np-1) * (x_i*p_l - p_lm1)
|
||
|
||
// The derivative of the residual is:
|
||
// \frac{d}{d x} \left[ (x^2 -1)P'_l(x) ] \right] =
|
||
// l * (l+1) * P_l(x), with l = np-1,
|
||
// therefore, deriv = np * (np-1) * p_l;
|
||
|
||
// compute dx = resid/deriv
|
||
real_t dx = (x_i*p_l - p_lm1) / (np*p_l);
|
||
#ifdef MFEM_USE_SINGLE
|
||
if (std::abs(dx) < 1e-7)
|
||
#elif defined MFEM_USE_DOUBLE
|
||
if (std::abs(dx) < 1e-16)
|
||
#else
|
||
MFEM_ABORT("Floating point type undefined");
|
||
if (std::abs(dx) < 1e-16)
|
||
#endif
|
||
{
|
||
done = true;
|
||
// Map the point to the interval [0,1]
|
||
z_i = ((1.0 + x_i) - dx)/2;
|
||
// continue the computation: get p_l at the new point, then exit
|
||
}
|
||
// If the iteration does not converge fast, something is wrong.
|
||
MFEM_VERIFY(iter < 8, "np = " << np << ", i = " << i
|
||
<< ", dx = " << dx);
|
||
// update x_i:
|
||
x_i -= dx;
|
||
}
|
||
// Map to the interval [0,1] and scale the weights
|
||
IntegrationPoint &ip = ir->IntPoint(i);
|
||
ip.x = z_i;
|
||
// w_i = (2/[ n*(n-1)*[P_{n-1}(x_i)]^2 ]) / 2
|
||
ip.weight = (real_t)(1.0 / (np*(np-1)*p_l*p_l));
|
||
|
||
// set the symmetric point
|
||
IntegrationPoint &symm_ip = ir->IntPoint(np-1-i);
|
||
symm_ip.x = 1.0 - z_i;
|
||
symm_ip.weight = ip.weight;
|
||
}
|
||
|
||
#else // MFEM_USE_MPFR is defined
|
||
|
||
HP_Quadrature1D hp_quad;
|
||
// use symmetry and compute just half of the points
|
||
for (int i = 0 ; i <= (np-1)/2 ; ++i)
|
||
{
|
||
hp_quad.ComputeGaussLobattoPoint(np, i);
|
||
ir->IntPoint(i).x = hp_quad.GetPoint();
|
||
ir->IntPoint(np-1-i).x = hp_quad.GetSymmPoint();
|
||
ir->IntPoint(i).weight =
|
||
ir->IntPoint(np-1-i).weight = hp_quad.GetWeight();
|
||
}
|
||
|
||
#endif // MFEM_USE_MPFR
|
||
|
||
}
|
||
}
|
||
|
||
void QuadratureFunctions1D::OpenUniform(const int np, IntegrationRule* ir)
|
||
{
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->SetOrder(np - 1 + np%2);
|
||
|
||
// The Newton-Cotes quadrature is based on weights that integrate exactly the
|
||
// interpolatory polynomial through the equally spaced quadrature points.
|
||
for (int i = 0; i < np ; ++i)
|
||
{
|
||
ir->IntPoint(i).x = real_t(i+1) / real_t(np + 1);
|
||
}
|
||
|
||
CalculateUniformWeights(ir, Quadrature1D::OpenUniform);
|
||
}
|
||
|
||
void QuadratureFunctions1D::ClosedUniform(const int np,
|
||
IntegrationRule* ir)
|
||
{
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->SetOrder(np - 1 + np%2);
|
||
if ( np == 1 ) // allow this case as "closed"
|
||
{
|
||
ir->IntPoint(0).Set1w(0.5, 1.0);
|
||
return;
|
||
}
|
||
|
||
for (int i = 0; i < np ; ++i)
|
||
{
|
||
ir->IntPoint(i).x = real_t(i) / (np-1);
|
||
}
|
||
|
||
CalculateUniformWeights(ir, Quadrature1D::ClosedUniform);
|
||
}
|
||
|
||
void QuadratureFunctions1D::OpenHalfUniform(const int np, IntegrationRule* ir)
|
||
{
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->SetOrder(np - 1 + np%2);
|
||
|
||
// Open half points: the centers of np uniform intervals
|
||
for (int i = 0; i < np ; ++i)
|
||
{
|
||
ir->IntPoint(i).x = real_t(2*i+1) / (2*np);
|
||
}
|
||
|
||
CalculateUniformWeights(ir, Quadrature1D::OpenHalfUniform);
|
||
}
|
||
|
||
void QuadratureFunctions1D::ClosedGL(const int np, IntegrationRule* ir)
|
||
{
|
||
ir->SetSize(np);
|
||
ir->SetPointIndices();
|
||
ir->IntPoint(0).x = 0.0;
|
||
ir->IntPoint(np-1).x = 1.0;
|
||
ir->SetOrder(np - 1 + np%2); // Is this the correct order?
|
||
|
||
if ( np > 2 )
|
||
{
|
||
IntegrationRule gl_ir;
|
||
GaussLegendre(np-1, &gl_ir);
|
||
|
||
for (int i = 1; i < np-1; ++i)
|
||
{
|
||
ir->IntPoint(i).x = (gl_ir.IntPoint(i-1).x + gl_ir.IntPoint(i).x)/2;
|
||
}
|
||
}
|
||
|
||
CalculateUniformWeights(ir, Quadrature1D::ClosedGL);
|
||
}
|
||
|
||
void QuadratureFunctions1D::GivePolyPoints(const int np, real_t *pts,
|
||
const int type)
|
||
{
|
||
IntegrationRule ir(np);
|
||
|
||
switch (type)
|
||
{
|
||
case Quadrature1D::GaussLegendre:
|
||
{
|
||
GaussLegendre(np,&ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::GaussLobatto:
|
||
{
|
||
GaussLobatto(np, &ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::OpenUniform:
|
||
{
|
||
OpenUniform(np,&ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::ClosedUniform:
|
||
{
|
||
ClosedUniform(np,&ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::OpenHalfUniform:
|
||
{
|
||
OpenHalfUniform(np, &ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::ClosedGL:
|
||
{
|
||
ClosedGL(np, &ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::Invalid:
|
||
{
|
||
MFEM_ABORT("Asking for an unknown type of 1D Quadrature points, "
|
||
"type = " << type);
|
||
}
|
||
}
|
||
|
||
for (int i = 0 ; i < np ; ++i)
|
||
{
|
||
pts[i] = ir.IntPoint(i).x;
|
||
}
|
||
}
|
||
|
||
void QuadratureFunctions1D::CalculateUniformWeights(IntegrationRule *ir,
|
||
const int type)
|
||
{
|
||
/* The Lagrange polynomials are:
|
||
p_i = \prod_{j \neq i} {\frac{x - x_j }{x_i - x_j}}
|
||
|
||
The weight associated with each abscissa is the integral of p_i over
|
||
[0,1]. To calculate the integral of p_i, we use a Gauss-Legendre
|
||
quadrature rule. This approach does not suffer from bad round-off/
|
||
cancellation errors for large number of points.
|
||
*/
|
||
const int n = ir->Size();
|
||
switch (n)
|
||
{
|
||
case 1:
|
||
ir->IntPoint(0).weight = 1.;
|
||
return;
|
||
case 2:
|
||
ir->IntPoint(0).weight = .5;
|
||
ir->IntPoint(1).weight = .5;
|
||
return;
|
||
}
|
||
|
||
#ifndef MFEM_USE_MPFR
|
||
|
||
// This algorithm should work for any set of points, not just uniform
|
||
const IntegrationRule &glob_ir = IntRules.Get(Geometry::SEGMENT, n-1);
|
||
const int m = glob_ir.GetNPoints();
|
||
Vector xv(n);
|
||
for (int j = 0; j < n; j++)
|
||
{
|
||
xv(j) = ir->IntPoint(j).x;
|
||
}
|
||
Poly_1D::Basis basis(n-1, xv.GetData()); // nodal basis, with nodes at 'xv'
|
||
Vector w(n);
|
||
// Integrate all nodal basis functions using 'glob_ir':
|
||
w = 0.0;
|
||
for (int i = 0; i < m; i++)
|
||
{
|
||
const IntegrationPoint &ip = glob_ir.IntPoint(i);
|
||
basis.Eval(ip.x, xv);
|
||
w.Add(ip.weight, xv); // w += ip.weight * xv
|
||
}
|
||
for (int j = 0; j < n; j++)
|
||
{
|
||
ir->IntPoint(j).weight = w(j);
|
||
}
|
||
|
||
#else // MFEM_USE_MPFR is defined
|
||
|
||
static const mpfr_rnd_t rnd = HP_Quadrature1D::rnd;
|
||
HP_Quadrature1D hp_quad;
|
||
mpfr_t l, lk, w0, wi, tmp, *weights;
|
||
mpfr_inits2(hp_quad.default_prec, l, lk, w0, wi, tmp, (mpfr_ptr) 0);
|
||
weights = new mpfr_t[n];
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
mpfr_init2(weights[i], hp_quad.default_prec);
|
||
mpfr_set_si(weights[i], 0, rnd);
|
||
}
|
||
hp_quad.SetRelTol(-48); // rtol = 2^(-48) ~ 3.5e-15
|
||
const int p = n-1;
|
||
const int m = p/2+1; // number of points for Gauss-Legendre quadrature
|
||
int hinv = 0, ihoffset = 0; // x_i = (i+ihoffset/2)/hinv
|
||
switch (type)
|
||
{
|
||
case Quadrature1D::ClosedUniform:
|
||
// x_i = i/p, i=0,...,p
|
||
hinv = p;
|
||
ihoffset = 0;
|
||
break;
|
||
case Quadrature1D::OpenUniform:
|
||
// x_i = (i+1)/(p+2), i=0,...,p
|
||
hinv = p+2;
|
||
ihoffset = 2;
|
||
break;
|
||
case Quadrature1D::OpenHalfUniform:
|
||
// x_i = (i+1/2)/(p+1), i=0,...,p
|
||
hinv = p+1;
|
||
ihoffset = 1;
|
||
break;
|
||
case Quadrature1D::GaussLegendre:
|
||
case Quadrature1D::GaussLobatto:
|
||
case Quadrature1D::ClosedGL:
|
||
case Quadrature1D::Invalid:
|
||
MFEM_ABORT("invalid Quadrature1D type: " << type);
|
||
}
|
||
// set w0 = (-1)^p*(p!)/(hinv^p)
|
||
mpfr_fac_ui(w0, p, rnd);
|
||
mpfr_ui_pow_ui(tmp, hinv, p, rnd);
|
||
mpfr_div(w0, w0, tmp, rnd);
|
||
if (p%2) { mpfr_neg(w0, w0, rnd); }
|
||
|
||
for (int j = 0; j < m; j++)
|
||
{
|
||
hp_quad.ComputeGaussLegendrePoint(m, j);
|
||
|
||
// Compute l = \prod_{i=0}^p (x-x_i) and lk = l/(x-x_k), where
|
||
// x = hp_quad.GetHPPoint(), x_i = (i+ihoffset/2)/hinv, and x_k is the
|
||
// node closest to x, i.e. k = min(max(round(x*hinv-ihoffset/2),0),p)
|
||
mpfr_mul_si(tmp, hp_quad.GetHPPoint(), hinv, rnd);
|
||
mpfr_sub_d(tmp, tmp, 0.5*ihoffset, rnd);
|
||
mpfr_round(tmp, tmp);
|
||
int k = min(max((int)mpfr_get_si(tmp, rnd), 0), p);
|
||
mpfr_set_si(lk, 1, rnd);
|
||
for (int i = 0; i <= p; i++)
|
||
{
|
||
mpfr_set_si(tmp, 2*i+ihoffset, rnd);
|
||
mpfr_div_si(tmp, tmp, 2*hinv, rnd);
|
||
mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
|
||
if (i != k)
|
||
{
|
||
mpfr_mul(lk, lk, tmp, rnd);
|
||
}
|
||
else
|
||
{
|
||
mpfr_set(l, tmp, rnd);
|
||
}
|
||
}
|
||
mpfr_mul(l, l, lk, rnd);
|
||
mpfr_set(wi, w0, rnd);
|
||
for (int i = 0; true; i++)
|
||
{
|
||
if (i != k)
|
||
{
|
||
// tmp = l/(wi*(x - x_i))
|
||
mpfr_set_si(tmp, 2*i+ihoffset, rnd);
|
||
mpfr_div_si(tmp, tmp, 2*hinv, rnd);
|
||
mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
|
||
mpfr_mul(tmp, tmp, wi, rnd);
|
||
mpfr_div(tmp, l, tmp, rnd);
|
||
}
|
||
else
|
||
{
|
||
// tmp = lk/wi
|
||
mpfr_div(tmp, lk, wi, rnd);
|
||
}
|
||
// weights[i] += hp_quad.weight*tmp
|
||
mpfr_mul(tmp, tmp, hp_quad.GetHPWeight(), rnd);
|
||
mpfr_add(weights[i], weights[i], tmp, rnd);
|
||
|
||
if (i == p) { break; }
|
||
|
||
// update wi *= (i+1)/(i-p)
|
||
mpfr_mul_si(wi, wi, i+1, rnd);
|
||
mpfr_div_si(wi, wi, i-p, rnd);
|
||
}
|
||
}
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
ir->IntPoint(i).weight = mpfr_get_d(weights[i], rnd);
|
||
mpfr_clear(weights[i]);
|
||
}
|
||
delete [] weights;
|
||
mpfr_clears(l, lk, w0, wi, tmp, (mpfr_ptr) 0);
|
||
|
||
#endif // MFEM_USE_MPFR
|
||
|
||
}
|
||
|
||
|
||
int Quadrature1D::CheckClosed(int type)
|
||
{
|
||
switch (type)
|
||
{
|
||
case GaussLobatto:
|
||
case ClosedUniform:
|
||
case ClosedGL:
|
||
return type;
|
||
default:
|
||
return Invalid;
|
||
}
|
||
}
|
||
|
||
int Quadrature1D::CheckOpen(int type)
|
||
{
|
||
switch (type)
|
||
{
|
||
case GaussLegendre:
|
||
case GaussLobatto:
|
||
case OpenUniform:
|
||
case ClosedUniform:
|
||
case OpenHalfUniform:
|
||
case ClosedGL:
|
||
return type; // all types can work as open
|
||
default:
|
||
return Invalid;
|
||
}
|
||
}
|
||
|
||
|
||
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre);
|
||
|
||
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre);
|
||
|
||
IntegrationRules::IntegrationRules(int ref, int type)
|
||
: quad_type(type)
|
||
{
|
||
refined = ref;
|
||
|
||
if (refined < 0) { own_rules = 0; return; }
|
||
|
||
own_rules = 1;
|
||
|
||
const MemoryType h_mt = MemoryType::HOST;
|
||
PointIntRules.SetSize(2, h_mt);
|
||
PointIntRules = NULL;
|
||
|
||
SegmentIntRules.SetSize(32, h_mt);
|
||
SegmentIntRules = NULL;
|
||
|
||
// TriangleIntegrationRule() assumes that this size is >= 26
|
||
TriangleIntRules.SetSize(32, h_mt);
|
||
TriangleIntRules = NULL;
|
||
|
||
SquareIntRules.SetSize(32, h_mt);
|
||
SquareIntRules = NULL;
|
||
|
||
// TetrahedronIntegrationRule() assumes that this size is >= 22
|
||
TetrahedronIntRules.SetSize(32, h_mt);
|
||
TetrahedronIntRules = NULL;
|
||
|
||
PyramidIntRules.SetSize(32, h_mt);
|
||
PyramidIntRules = NULL;
|
||
|
||
PrismIntRules.SetSize(32, h_mt);
|
||
PrismIntRules = NULL;
|
||
|
||
CubeIntRules.SetSize(32, h_mt);
|
||
CubeIntRules = NULL;
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
IntRuleLocks.SetSize(Geometry::NUM_GEOMETRIES, h_mt);
|
||
for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
|
||
{
|
||
omp_init_lock(&IntRuleLocks[i]);
|
||
}
|
||
#endif
|
||
}
|
||
|
||
const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
|
||
{
|
||
Array<IntegrationRule *> *ir_array = NULL;
|
||
|
||
switch (GeomType)
|
||
{
|
||
case Geometry::POINT: ir_array = &PointIntRules; Order = 0; break;
|
||
case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
|
||
case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
|
||
case Geometry::SQUARE: ir_array = &SquareIntRules; break;
|
||
case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
|
||
case Geometry::CUBE: ir_array = &CubeIntRules; break;
|
||
case Geometry::PRISM: ir_array = &PrismIntRules; break;
|
||
case Geometry::PYRAMID: ir_array = &PyramidIntRules; break;
|
||
case Geometry::INVALID:
|
||
case Geometry::NUM_GEOMETRIES:
|
||
MFEM_ABORT("Unknown type of reference element!");
|
||
}
|
||
|
||
if (Order < 0)
|
||
{
|
||
Order = 0;
|
||
}
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_set_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
|
||
if (!HaveIntRule(*ir_array, Order))
|
||
{
|
||
IntegrationRule *ir = GenerateIntegrationRule(GeomType, Order);
|
||
#ifdef MFEM_DEBUG
|
||
int RealOrder = Order;
|
||
while (RealOrder+1 < ir_array->Size() && (*ir_array)[RealOrder+1] == ir)
|
||
{
|
||
RealOrder++;
|
||
}
|
||
MFEM_VERIFY(RealOrder == ir->GetOrder(), "internal error");
|
||
#else
|
||
MFEM_CONTRACT_VAR(ir);
|
||
#endif
|
||
}
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_unset_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
|
||
return *(*ir_array)[Order];
|
||
}
|
||
|
||
void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
|
||
{
|
||
Array<IntegrationRule *> *ir_array = NULL;
|
||
|
||
switch (GeomType)
|
||
{
|
||
case Geometry::POINT: ir_array = &PointIntRules; break;
|
||
case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
|
||
case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
|
||
case Geometry::SQUARE: ir_array = &SquareIntRules; break;
|
||
case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
|
||
case Geometry::CUBE: ir_array = &CubeIntRules; break;
|
||
case Geometry::PRISM: ir_array = &PrismIntRules; break;
|
||
case Geometry::PYRAMID: ir_array = &PyramidIntRules; break;
|
||
case Geometry::INVALID:
|
||
case Geometry::NUM_GEOMETRIES:
|
||
MFEM_ABORT("Unknown type of reference element!");
|
||
}
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_set_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
|
||
if (HaveIntRule(*ir_array, Order))
|
||
{
|
||
MFEM_ABORT("Overwriting set rules is not supported!");
|
||
}
|
||
|
||
AllocIntRule(*ir_array, Order);
|
||
|
||
(*ir_array)[Order] = &IntRule;
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_unset_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
}
|
||
|
||
void IntegrationRules::DeleteIntRuleArray(
|
||
Array<IntegrationRule *> &ir_array) const
|
||
{
|
||
// Many of the intrules have multiple contiguous copies in the ir_array
|
||
// so we have to be careful to not delete them twice.
|
||
IntegrationRule *ir = NULL;
|
||
for (int i = 0; i < ir_array.Size(); i++)
|
||
{
|
||
if (ir_array[i] != NULL && ir_array[i] != ir)
|
||
{
|
||
ir = ir_array[i];
|
||
delete ir;
|
||
}
|
||
}
|
||
}
|
||
|
||
IntegrationRules::~IntegrationRules()
|
||
{
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
|
||
{
|
||
omp_destroy_lock(&IntRuleLocks[i]);
|
||
}
|
||
#endif
|
||
|
||
if (!own_rules) { return; }
|
||
|
||
DeleteIntRuleArray(PointIntRules);
|
||
DeleteIntRuleArray(SegmentIntRules);
|
||
DeleteIntRuleArray(TriangleIntRules);
|
||
DeleteIntRuleArray(SquareIntRules);
|
||
DeleteIntRuleArray(TetrahedronIntRules);
|
||
DeleteIntRuleArray(CubeIntRules);
|
||
DeleteIntRuleArray(PrismIntRules);
|
||
DeleteIntRuleArray(PyramidIntRules);
|
||
}
|
||
|
||
|
||
IntegrationRule *IntegrationRules::GenerateIntegrationRule(int GeomType,
|
||
int Order)
|
||
{
|
||
switch (GeomType)
|
||
{
|
||
case Geometry::POINT:
|
||
return PointIntegrationRule(Order);
|
||
case Geometry::SEGMENT:
|
||
return SegmentIntegrationRule(Order);
|
||
case Geometry::TRIANGLE:
|
||
return TriangleIntegrationRule(Order);
|
||
case Geometry::SQUARE:
|
||
return SquareIntegrationRule(Order);
|
||
case Geometry::TETRAHEDRON:
|
||
return TetrahedronIntegrationRule(Order);
|
||
case Geometry::CUBE:
|
||
return CubeIntegrationRule(Order);
|
||
case Geometry::PRISM:
|
||
return PrismIntegrationRule(Order);
|
||
case Geometry::PYRAMID:
|
||
return PyramidIntegrationRule(Order);
|
||
case Geometry::INVALID:
|
||
case Geometry::NUM_GEOMETRIES:
|
||
MFEM_ABORT("Unknown type of reference element!");
|
||
}
|
||
return NULL;
|
||
}
|
||
|
||
|
||
// Integration rules for a point
|
||
IntegrationRule *IntegrationRules::PointIntegrationRule(int Order)
|
||
{
|
||
if (Order > 1)
|
||
{
|
||
MFEM_ABORT("Point Integration Rule of Order > 1 not defined");
|
||
return NULL;
|
||
}
|
||
|
||
IntegrationRule *ir = new IntegrationRule(1);
|
||
ir->IntPoint(0).x = .0;
|
||
ir->IntPoint(0).weight = 1.;
|
||
ir->SetOrder(1);
|
||
|
||
PointIntRules[1] = PointIntRules[0] = ir;
|
||
|
||
return ir;
|
||
}
|
||
|
||
// Integration rules for line segment [0,1]
|
||
IntegrationRule *IntegrationRules::SegmentIntegrationRule(int Order)
|
||
{
|
||
int RealOrder = GetSegmentRealOrder(Order); // RealOrder >= Order
|
||
// Order is one of {RealOrder-1,RealOrder}
|
||
AllocIntRule(SegmentIntRules, RealOrder);
|
||
|
||
IntegrationRule *ir = new IntegrationRule;
|
||
|
||
int n = 0;
|
||
// n is the number of points to achieve the exact integral of a
|
||
// degree Order polynomial
|
||
switch (quad_type)
|
||
{
|
||
case Quadrature1D::GaussLegendre:
|
||
{
|
||
// Gauss-Legendre is exact for 2*n-1
|
||
n = Order/2 + 1;
|
||
QuadratureFunctions1D::GaussLegendre(n, ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::GaussLobatto:
|
||
{
|
||
// Gauss-Lobatto is exact for 2*n-3
|
||
n = Order/2 + 2;
|
||
QuadratureFunctions1D::GaussLobatto(n, ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::OpenUniform:
|
||
{
|
||
// Open Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
|
||
n = Order | 1; // n is always odd
|
||
QuadratureFunctions1D::OpenUniform(n, ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::ClosedUniform:
|
||
{
|
||
// Closed Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
|
||
n = Order | 1; // n is always odd
|
||
QuadratureFunctions1D::ClosedUniform(n, ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::OpenHalfUniform:
|
||
{
|
||
// Open half Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
|
||
n = Order | 1; // n is always odd
|
||
QuadratureFunctions1D::OpenHalfUniform(n, ir);
|
||
break;
|
||
}
|
||
case Quadrature1D::Invalid:
|
||
{
|
||
MFEM_ABORT("unknown Quadrature1D type: " << quad_type);
|
||
}
|
||
}
|
||
if (refined)
|
||
{
|
||
// Effectively passing memory management to SegmentIntegrationRules
|
||
IntegrationRule *refined_ir = new IntegrationRule(2*n);
|
||
refined_ir->SetOrder(ir->GetOrder());
|
||
for (int j = 0; j < n; j++)
|
||
{
|
||
refined_ir->IntPoint(j).x = ir->IntPoint(j).x/2.0;
|
||
refined_ir->IntPoint(j).weight = ir->IntPoint(j).weight/2.0;
|
||
refined_ir->IntPoint(j+n).x = 0.5 + ir->IntPoint(j).x/2.0;
|
||
refined_ir->IntPoint(j+n).weight = ir->IntPoint(j).weight/2.0;
|
||
}
|
||
delete ir;
|
||
ir = refined_ir;
|
||
}
|
||
SegmentIntRules[RealOrder-1] = SegmentIntRules[RealOrder] = ir;
|
||
return ir;
|
||
}
|
||
|
||
// Triangle rules from Witherden & Vincent [1].
|
||
// Orbit data from PyFR (https://pyfr.org), licensed under CC-BY 4.0.
|
||
IntegrationRule *IntegrationRules::TriangleIntegrationRule(int Order)
|
||
{
|
||
IntegrationRule *ir = NULL;
|
||
|
||
switch (Order)
|
||
{
|
||
case 0:
|
||
case 1:
|
||
ir = new IntegrationRule(1);
|
||
ir->AddTriMidPoint(0, 0.5);
|
||
ir->SetOrder(1);
|
||
TriangleIntRules[0] =
|
||
TriangleIntRules[1] = ir;
|
||
return ir;
|
||
|
||
case 2:
|
||
ir = new IntegrationRule(3);
|
||
ir->AddTriPoints3(0, 1./6., 1./6.);
|
||
ir->SetOrder(2);
|
||
TriangleIntRules[2] = ir;
|
||
return ir;
|
||
|
||
case 3:
|
||
case 4:
|
||
ir = new IntegrationRule(6);
|
||
ir->AddTriPoints3(0, 4.45948490915964890213e-01, 1.11690794839005735906e-01);
|
||
ir->AddTriPoints3(3, 9.15762135097707430376e-02, 5.49758718276609353870e-02);
|
||
ir->SetOrder(4);
|
||
TriangleIntRules[3] =
|
||
TriangleIntRules[4] = ir;
|
||
return ir;
|
||
|
||
case 5:
|
||
ir = new IntegrationRule(7);
|
||
ir->AddTriMidPoint(0, 0.1125);
|
||
ir->AddTriPoints3(1, 1.01286507323456342888e-01, 6.29695902724135697648e-02);
|
||
ir->AddTriPoints3(4, 4.70142064105115109474e-01, 6.61970763942530959767e-02);
|
||
ir->SetOrder(5);
|
||
TriangleIntRules[5] = ir;
|
||
return ir;
|
||
|
||
case 6:
|
||
ir = new IntegrationRule(12);
|
||
ir->AddTriPoints3(0, 6.30890144915022266225e-02, 2.54224531851034094010e-02);
|
||
ir->AddTriPoints3(3, 2.49286745170910428726e-01, 5.83931378631896841336e-02);
|
||
ir->AddTriPoints6(6, 6.36502499121398668258e-01, 3.10352451033784393353e-01,
|
||
4.14255378091867854096e-02);
|
||
ir->SetOrder(6);
|
||
TriangleIntRules[6] = ir;
|
||
return ir;
|
||
|
||
case 7:
|
||
ir = new IntegrationRule(15);
|
||
ir->AddTriPoints3(0, 3.37306485545878498300e-02, 8.27252505539606552976e-03);
|
||
ir->AddTriPoints3(3, 2.41577382595403566956e-01, 6.39720856150777922311e-02);
|
||
ir->AddTriPoints3(6, 4.74309692504718327655e-01, 3.85433230929930342734e-02);
|
||
ir->AddTriPoints6(9, 7.54280040550053154647e-01, 1.98683314797351684433e-01,
|
||
2.79393664515998896292e-02);
|
||
ir->SetOrder(7);
|
||
TriangleIntRules[7] = ir;
|
||
return ir;
|
||
|
||
case 8:
|
||
ir = new IntegrationRule(16);
|
||
ir->AddTriMidPoint(0, 7.21578038388935860681e-02);
|
||
ir->AddTriPoints3(1, 4.59292588292723236165e-01, 4.75458171336423096598e-02);
|
||
ir->AddTriPoints3(4, 1.70569307751760268488e-01, 5.16086852673591223173e-02);
|
||
ir->AddTriPoints3(7, 5.05472283170309566458e-02, 1.62292488115990396480e-02);
|
||
ir->AddTriPoints6(10, 7.28492392955404244326e-01, 2.63112829634638112353e-01,
|
||
1.36151570872174963733e-02);
|
||
ir->SetOrder(8);
|
||
TriangleIntRules[8] = ir;
|
||
return ir;
|
||
|
||
case 9:
|
||
ir = new IntegrationRule(19);
|
||
ir->AddTriMidPoint(0, 4.85678981413994181882e-02);
|
||
ir->AddTriPoints3(1, 4.37089591492936690997e-01, 3.89137705023871391385e-02);
|
||
ir->AddTriPoints3(4, 1.88203535619032802373e-01, 3.98238694636051243636e-02);
|
||
ir->AddTriPoints3(7, 4.89682519198737620236e-01, 1.56673501135695357467e-02);
|
||
ir->AddTriPoints3(10, 4.47295133944527467662e-02, 1.27888378293490156262e-02);
|
||
ir->AddTriPoints6(13, 7.41198598784498008385e-01, 2.21962989160765733487e-01,
|
||
2.16417696886446880855e-02);
|
||
ir->SetOrder(9);
|
||
TriangleIntRules[9] = ir;
|
||
return ir;
|
||
|
||
case 10:
|
||
ir = new IntegrationRule(25);
|
||
ir->AddTriMidPoint(0, 4.08716645731429864541e-02);
|
||
ir->AddTriPoints3(1, 3.20553732169435168231e-02, 6.67648440657478327992e-03);
|
||
ir->AddTriPoints3(4, 1.42161101056564431744e-01, 2.29789818023723654838e-02);
|
||
ir->AddTriPoints6(7, 5.30054118927343997925e-01, 3.21812995288835446139e-01,
|
||
3.19524531982120219009e-02);
|
||
ir->AddTriPoints6(13, 6.01233328683459244957e-01, 3.69146781827810910315e-01,
|
||
1.70923240814797143539e-02);
|
||
ir->AddTriPoints6(19, 8.07930600922879049719e-01, 1.63701733737182442141e-01,
|
||
1.26488788536441923438e-02);
|
||
ir->SetOrder(10);
|
||
TriangleIntRules[10] = ir;
|
||
return ir;
|
||
|
||
case 11:
|
||
ir = new IntegrationRule(28);
|
||
ir->AddTriMidPoint(0, 4.28805898661121093207e-02);
|
||
ir->AddTriPoints3(1, 2.84854176143718995640e-02, 5.21593525644734826857e-03);
|
||
ir->AddTriPoints3(4, 2.10219956703178278978e-01, 3.52578420558582877886e-02);
|
||
ir->AddTriPoints3(7, 1.02635482712246428605e-01, 1.93153796185096607307e-02);
|
||
ir->AddTriPoints3(10, 4.95891900965890919384e-01, 8.30313652729268436570e-03);
|
||
ir->AddTriPoints3(13, 4.38465926764352253997e-01, 3.36580770397341480504e-02);
|
||
ir->AddTriPoints6(16, 8.43349783661853091843e-01, 1.49324788652082374174e-01,
|
||
5.14514478647663895533e-03);
|
||
ir->AddTriPoints6(22, 6.64408374196864159877e-01, 2.89581125637705882880e-01,
|
||
2.01662383202502772106e-02);
|
||
ir->SetOrder(11);
|
||
TriangleIntRules[11] = ir;
|
||
return ir;
|
||
|
||
case 12:
|
||
ir = new IntegrationRule(33);
|
||
ir->AddTriPoints3(0, 4.88203750945541581352e-01, 1.21334190407260157640e-02);
|
||
ir->AddTriPoints3(3, 1.09257827659354322947e-01, 1.42430260344387719235e-02);
|
||
ir->AddTriPoints3(6, 2.71462507014926135440e-01, 3.12706065979513822550e-02);
|
||
ir->AddTriPoints3(9, 2.46463634363356387524e-02, 3.96582125498681943576e-03);
|
||
ir->AddTriPoints3(12, 4.40111648658593201944e-01, 2.49591674640304711508e-02);
|
||
ir->AddTriPoints6(15, 6.85310163906391878186e-01, 2.91655679738340944951e-01,
|
||
1.08917925193037796322e-02);
|
||
ir->AddTriPoints6(21, 6.28249751683556123538e-01, 2.55454228638517299999e-01,
|
||
2.16136818297071042760e-02);
|
||
ir->AddTriPoints6(27, 8.51337792510240110033e-01, 1.27279717233589384495e-01,
|
||
7.54183878825571887144e-03);
|
||
ir->SetOrder(12);
|
||
TriangleIntRules[12] = ir;
|
||
return ir;
|
||
|
||
case 13:
|
||
ir = new IntegrationRule(37);
|
||
ir->AddTriMidPoint(0, 3.39800182934158201409e-02);
|
||
ir->AddTriPoints3(1, 4.89076946452539351728e-01, 1.19972009644473652512e-02);
|
||
ir->AddTriPoints3(4, 2.21372286291832920391e-01, 2.91392425595999905730e-02);
|
||
ir->AddTriPoints3(7, 4.26941414259800422482e-01, 2.78009837652266646180e-02);
|
||
ir->AddTriPoints3(10, 2.15096811088433259584e-02, 3.02616855176958583773e-03);
|
||
ir->AddTriPoints6(13, 7.48507115899952224503e-01, 1.63597401067850478640e-01,
|
||
1.20895199057969096601e-02);
|
||
ir->AddTriPoints6(19, 8.64707770295442768038e-01, 1.10922042803463405392e-01,
|
||
7.48270055258283377231e-03);
|
||
ir->AddTriPoints6(25, 6.23545995553675513889e-01, 3.08441760892117777804e-01,
|
||
1.73206380704241866275e-02);
|
||
ir->AddTriPoints6(31, 7.22357793124188019007e-01, 2.72515817773429591675e-01,
|
||
4.79534050177163155559e-03);
|
||
ir->SetOrder(13);
|
||
TriangleIntRules[13] = ir;
|
||
return ir;
|
||
|
||
case 14:
|
||
ir = new IntegrationRule(42);
|
||
ir->AddTriPoints3(0, 1.77205532412543442788e-01, 2.10812943684965080349e-02);
|
||
ir->AddTriPoints3(3, 4.17644719340453940415e-01, 1.63941767720626740967e-02);
|
||
ir->AddTriPoints3(6, 6.17998830908725871325e-02, 7.21684983488833382143e-03);
|
||
ir->AddTriPoints3(9, 4.88963910362178677538e-01, 1.09417906847144447147e-02);
|
||
ir->AddTriPoints3(12, 2.73477528308838646609e-01, 2.58870522536457925433e-02);
|
||
ir->AddTriPoints3(15, 1.93909612487010996063e-02, 2.46170180120004094063e-03);
|
||
ir->AddTriPoints6(18, 6.86980167808087793802e-01, 2.98372882136257788765e-01,
|
||
7.21815405676692022074e-03);
|
||
ir->AddTriPoints6(24, 7.70608554774996457049e-01, 1.72266687821355679588e-01,
|
||
1.23328766062818367955e-02);
|
||
ir->AddTriPoints6(30, 5.70222290846683188548e-01, 3.36861459796344964168e-01,
|
||
1.92857553935303419057e-02);
|
||
ir->AddTriPoints6(36, 8.79757171370171064950e-01, 1.18974497696956893478e-01,
|
||
2.50511441925033596229e-03);
|
||
ir->SetOrder(14);
|
||
TriangleIntRules[14] = ir;
|
||
return ir;
|
||
|
||
case 15:
|
||
ir = new IntegrationRule(49);
|
||
ir->AddTriMidPoint(0, 2.21676936910920364954e-02);
|
||
ir->AddTriPoints3(1, 4.05362214133975495844e-01, 2.13568907857302828224e-02);
|
||
ir->AddTriPoints3(4, 7.01735528999860580512e-02, 8.22236878131258133728e-03);
|
||
ir->AddTriPoints3(7, 4.74170681438019769871e-01, 8.69807400038170690226e-03);
|
||
ir->AddTriPoints3(10, 2.26378713420349653163e-01, 2.33916808643548149171e-02);
|
||
ir->AddTriPoints3(13, 4.94996956769126195130e-01, 4.78692309123004283711e-03);
|
||
ir->AddTriPoints3(16, 1.58117262509887002153e-02, 1.48038731895268772104e-03);
|
||
ir->AddTriPoints6(19, 6.66975644801868106093e-01, 3.14648242812450851247e-01,
|
||
7.80128641528798211224e-03);
|
||
ir->AddTriPoints6(25, 9.19912157726236134891e-01, 7.09486052364554087291e-02,
|
||
2.01492668600904969653e-03);
|
||
ir->AddTriPoints6(31, 7.15222356931450642392e-01, 1.90535589476393929509e-01,
|
||
1.43602934626006709801e-02);
|
||
ir->AddTriPoints6(37, 8.13292641049419229304e-01, 1.68068645222414381202e-01,
|
||
5.83631059078792285844e-03);
|
||
ir->AddTriPoints6(43, 5.65252664877114230357e-01, 3.38950611475277163720e-01,
|
||
1.56577381424846430458e-02);
|
||
ir->SetOrder(15);
|
||
TriangleIntRules[15] = ir;
|
||
return ir;
|
||
|
||
case 16:
|
||
ir = new IntegrationRule(55);
|
||
ir->AddTriMidPoint(0, 2.26322830369093952463e-02);
|
||
ir->AddTriPoints3(1, 2.45990070467141719313e-01, 2.05464615718494759966e-02);
|
||
ir->AddTriPoints3(4, 4.15584896885420551627e-01, 2.03559166562126796218e-02);
|
||
ir->AddTriPoints3(7, 8.53555665867003487968e-02, 7.39081734511220188322e-03);
|
||
ir->AddTriPoints3(10, 1.61918644191271221544e-01, 1.47092048494940497855e-02);
|
||
ir->AddTriPoints3(13, 5.00000000000000000000e-01, 2.20927315607528452004e-03);
|
||
ir->AddTriPoints3(16, 4.75280727545942083268e-01, 1.29871666491385793357e-02);
|
||
ir->AddTriPoints6(19, 7.54170061444767725334e-01, 1.91074763640529221576e-01,
|
||
9.46913623220784969603e-03);
|
||
ir->AddTriPoints6(25, 9.68244368030958701965e-01, 2.32034277688137335893e-02,
|
||
8.27233357417524097638e-04);
|
||
ir->AddTriPoints6(31, 6.49303698245446425652e-01, 3.31764523474147643434e-01,
|
||
7.50430089214290316213e-03);
|
||
ir->AddTriPoints6(37, 9.00273703270429548340e-01, 8.06961669858730079596e-02,
|
||
3.97379696669624901673e-03);
|
||
ir->AddTriPoints6(43, 5.89148840564247877616e-01, 3.08244969196354023921e-01,
|
||
1.59918050396850343342e-02);
|
||
ir->AddTriPoints6(49, 8.06621867499395683865e-01, 1.87441782483782071189e-01,
|
||
2.69559355842440570919e-03);
|
||
ir->SetOrder(16);
|
||
TriangleIntRules[16] = ir;
|
||
return ir;
|
||
|
||
case 17:
|
||
ir = new IntegrationRule(60);
|
||
ir->AddTriPoints3(0, 4.17103444361599295931e-01, 1.36554632640510532210e-02);
|
||
ir->AddTriPoints3(3, 1.47554916607539610141e-02, 1.38694378881882109979e-03);
|
||
ir->AddTriPoints3(6, 4.65597871618890324363e-01, 1.25097254752486782697e-02);
|
||
ir->AddTriPoints3(9, 1.80358116266370605008e-01, 1.31563152940089925225e-02);
|
||
ir->AddTriPoints3(12, 6.66540634795969033632e-02, 6.22950040115272107855e-03);
|
||
ir->AddTriPoints3(15, 2.85706502436586629035e-01, 1.88581185763976415248e-02);
|
||
ir->AddTriPoints6(18, 8.24790070165088096132e-01, 1.59192287472792681768e-01,
|
||
3.98915010296479674579e-03);
|
||
ir->AddTriPoints6(24, 6.26369030386452196879e-01, 3.06281591746186521164e-01,
|
||
1.12438862733455335191e-02);
|
||
ir->AddTriPoints6(30, 5.71294867944684092720e-01, 4.15475459295228999324e-01,
|
||
5.19921997791976831654e-03);
|
||
ir->AddTriPoints6(36, 7.53235145936458128091e-01, 1.68722513495259462957e-01,
|
||
1.02789491602272593102e-02);
|
||
ir->AddTriPoints6(42, 7.15072259110642427515e-01, 2.71791870055354878311e-01,
|
||
4.34610725050059605590e-03);
|
||
ir->AddTriPoints6(48, 9.15919353297816929427e-01, 7.25054707990024915887e-02,
|
||
2.29217420086793351869e-03);
|
||
ir->AddTriPoints6(54, 5.43275579596159796658e-01, 2.99218942476970228839e-01,
|
||
1.30858129676684944304e-02);
|
||
ir->SetOrder(17);
|
||
TriangleIntRules[17] = ir;
|
||
return ir;
|
||
|
||
case 18:
|
||
ir = new IntegrationRule(67);
|
||
ir->AddTriMidPoint(0, 1.81778676507133342410e-02);
|
||
ir->AddTriPoints3(1, 3.99955628067576229867e-01, 1.66522350166950668104e-02);
|
||
ir->AddTriPoints3(4, 4.87580301574869645620e-01, 6.02332381699985548729e-03);
|
||
ir->AddTriPoints3(7, 4.61809506406449243876e-01, 9.47458575338943308208e-03);
|
||
ir->AddTriPoints3(10, 2.42264702514271956790e-01, 1.82375447044718190515e-02);
|
||
ir->AddTriPoints3(13, 3.88302560886856218403e-02, 3.56466300985948522304e-03);
|
||
ir->AddTriPoints3(16, 9.19477421216432500017e-02, 8.27957997600162372287e-03);
|
||
ir->AddTriPoints6(19, 7.70372376214675247397e-01, 1.83822707925463957324e-01,
|
||
6.87980811747110256732e-03);
|
||
ir->AddTriPoints6(25, 6.70953985194234547862e-01, 2.06349257433837918185e-01,
|
||
1.18909554500764153007e-02);
|
||
ir->AddTriPoints6(31, 6.00418954634256873959e-01, 3.95683434332269712286e-01,
|
||
2.26526725112853252742e-03);
|
||
ir->AddTriPoints6(37, 8.78342189467521738955e-01, 1.08195793791033278985e-01,
|
||
3.42005505980359086893e-03);
|
||
ir->AddTriPoints6(43, 6.39988092004714625993e-01, 3.19751624525377309283e-01,
|
||
8.87374455101020212511e-03);
|
||
ir->AddTriPoints6(49, 7.58929479855198430016e-01, 2.35772184958191743931e-01,
|
||
2.50533043728986106261e-03);
|
||
ir->AddTriPoints6(55, 9.72360728962795684005e-01, 2.70909109951620319379e-02,
|
||
6.11474063480544911126e-04);
|
||
ir->AddTriPoints6(61, 5.45918775386194599086e-01, 3.33493529449880754534e-01,
|
||
1.27410876559122202695e-02);
|
||
ir->SetOrder(18);
|
||
TriangleIntRules[18] = ir;
|
||
return ir;
|
||
|
||
case 19:
|
||
ir = new IntegrationRule(73);
|
||
ir->AddTriMidPoint(0, 1.72346988520061666916e-02);
|
||
ir->AddTriPoints3(1, 5.25238903512089683190e-02, 3.55462829889906543543e-03);
|
||
ir->AddTriPoints3(4, 4.92512675041336889237e-01, 5.16087757147214078873e-03);
|
||
ir->AddTriPoints3(7, 1.11448873323021391268e-01, 7.61717554650914990128e-03);
|
||
ir->AddTriPoints3(10, 4.59194201039543670184e-01, 1.14917950133708035576e-02);
|
||
ir->AddTriPoints3(13, 4.03969722551901222474e-01, 1.57687674465774863020e-02);
|
||
ir->AddTriPoints3(16, 1.78170104781764315760e-01, 1.23259574240954274116e-02);
|
||
ir->AddTriPoints3(19, 1.16394611837894457196e-02, 8.82661388221423837477e-04);
|
||
ir->AddTriPoints3(22, 2.55161632913607716588e-01, 1.58765096830015377261e-02);
|
||
ir->AddTriPoints6(25, 8.30156464400275351245e-01, 1.30697676268032414448e-01,
|
||
4.84774224342752330097e-03);
|
||
ir->AddTriPoints6(31, 5.59369805720300927732e-01, 3.11317629809541251973e-01,
|
||
1.31731609886953666272e-02);
|
||
ir->AddTriPoints6(37, 6.33313293128784149388e-01, 3.64617780974611060962e-01,
|
||
1.64103827591790965915e-03);
|
||
ir->AddTriPoints6(43, 7.04004819966042139079e-01, 2.21434885432331141075e-01,
|
||
9.05397246560622585843e-03);
|
||
ir->AddTriPoints6(49, 8.52566954376889230005e-01, 1.42425757365756355810e-01,
|
||
1.46315755173510018451e-03);
|
||
ir->AddTriPoints6(55, 6.05083979068707922266e-01, 3.54028009735275261960e-01,
|
||
8.05108138201205379703e-03);
|
||
ir->AddTriPoints6(61, 7.43181368957436361278e-01, 2.41894578960579587079e-01,
|
||
4.22794374976824798712e-03);
|
||
ir->AddTriPoints6(67, 9.30137698876805085746e-01, 6.00862753223067036501e-02,
|
||
1.66360068142969402642e-03);
|
||
ir->SetOrder(19);
|
||
TriangleIntRules[19] = ir;
|
||
return ir;
|
||
|
||
case 20:
|
||
ir = new IntegrationRule(79);
|
||
ir->AddTriMidPoint(0, 1.39101107014531159140e-02);
|
||
ir->AddTriPoints3(1, 2.54579267673339160183e-01, 1.40832013075202471669e-02);
|
||
ir->AddTriPoints3(4, 1.09761410283977789426e-02, 7.98840791066619858654e-04);
|
||
ir->AddTriPoints3(7, 1.09383596711714603522e-01, 7.83023077607453328597e-03);
|
||
ir->AddTriPoints3(10, 1.86294997744540946627e-01, 9.17346297425291473671e-03);
|
||
ir->AddTriPoints3(13, 4.45551056955924895675e-01, 9.45239993323244813428e-03);
|
||
ir->AddTriPoints3(16, 3.73108805988847103130e-02, 2.16127541066557732688e-03);
|
||
ir->AddTriPoints3(19, 3.93425347817099924086e-01, 1.37880506290704585304e-02);
|
||
ir->AddTriPoints3(22, 4.76245611540499047543e-01, 7.10182530340844071076e-03);
|
||
ir->AddTriPoints6(25, 8.33295511838236246938e-01, 1.59133707657067247077e-01,
|
||
2.20289741855849742491e-03);
|
||
ir->AddTriPoints6(31, 7.54921502863547422280e-01, 1.98518132228788335425e-01,
|
||
5.98639857895469015836e-03);
|
||
ir->AddTriPoints6(37, 9.31054476783942153162e-01, 6.40905856084340586065e-02,
|
||
1.12986960212586558597e-03);
|
||
ir->AddTriPoints6(43, 6.11877703547425655373e-01, 3.33134817309587605294e-01,
|
||
8.66722556721933289070e-03);
|
||
ir->AddTriPoints6(49, 8.61684018936486717521e-01, 9.99522962881386756173e-02,
|
||
4.14571152761385782609e-03);
|
||
ir->AddTriPoints6(55, 6.78165737889635522606e-01, 2.15607057390094447591e-01,
|
||
7.72260782209923009323e-03);
|
||
ir->AddTriPoints6(61, 5.70144692890973359134e-01, 4.20023758816224113133e-01,
|
||
3.69568150025529782929e-03);
|
||
ir->AddTriPoints6(67, 5.42331804172428100230e-01, 3.17860123835772001577e-01,
|
||
1.16917457318277372147e-02);
|
||
ir->AddTriPoints6(73, 7.08681375720323636358e-01, 2.80581411423665327831e-01,
|
||
3.57820023845768515197e-03);
|
||
ir->SetOrder(20);
|
||
TriangleIntRules[20] = ir;
|
||
return ir;
|
||
|
||
case 21:
|
||
case 22:
|
||
case 23:
|
||
case 24:
|
||
case 25:
|
||
ir = new IntegrationRule(126);
|
||
ir->AddTriPoints3b(0, 0.0279464830731742, 0.0040027909400102085);
|
||
ir->AddTriPoints3b(3, 0.131178601327651467, 0.00797353841619525);
|
||
ir->AddTriPoints3b(6, 0.220221729512072267, 0.006554570615397765);
|
||
ir->AddTriPoints3 (9, 0.298443234019804467, 0.00979150048281781);
|
||
ir->AddTriPoints3(12, 0.2340441723373718, 0.008235442720768635);
|
||
ir->AddTriPoints3(15, 0.151468334609017567, 0.00427363953704605);
|
||
ir->AddTriPoints3(18, 0.112733893545993667, 0.004080942928613246);
|
||
ir->AddTriPoints3(21, 0.0777156920915263, 0.0030605732699918895);
|
||
ir->AddTriPoints3(24, 0.034893093614297, 0.0014542491324683325);
|
||
ir->AddTriPoints3(27, 0.00725818462093236667, 0.00034613762283099815);
|
||
ir->AddTriPoints6(30, 0.0012923527044422, 0.227214452153364077,
|
||
0.0006241445996386985);
|
||
ir->AddTriPoints6(36, 0.0053997012721162, 0.435010554853571706,
|
||
0.001702376454401511);
|
||
ir->AddTriPoints6(42, 0.006384003033975, 0.320309599272204437,
|
||
0.0016798271630320255);
|
||
ir->AddTriPoints6(48, 0.00502821150199306667, 0.0917503222800051889,
|
||
0.000858078269748377);
|
||
ir->AddTriPoints6(54, 0.00682675862178186667, 0.0380108358587243835,
|
||
0.000740428158357803);
|
||
ir->AddTriPoints6(60, 0.0100161996399295333, 0.157425218485311668,
|
||
0.0017556563053643425);
|
||
ir->AddTriPoints6(66, 0.02575781317339, 0.239889659778533193,
|
||
0.003696775074853242);
|
||
ir->AddTriPoints6(72, 0.0302278981199158, 0.361943118126060531,
|
||
0.003991543738688279);
|
||
ir->AddTriPoints6(78, 0.0305049901071620667, 0.0835519609548285602,
|
||
0.0021779813065790205);
|
||
ir->AddTriPoints6(84, 0.0459565473625693333, 0.148443220732418205,
|
||
0.003682528350708916);
|
||
ir->AddTriPoints6(90, 0.0674428005402775333, 0.283739708727534955,
|
||
0.005481786423209775);
|
||
ir->AddTriPoints6(96, 0.0700450914159106, 0.406899375118787573,
|
||
0.00587498087177056);
|
||
ir->AddTriPoints6(102, 0.0839115246401166, 0.194113987024892542,
|
||
0.005007800356899285);
|
||
ir->AddTriPoints6(108, 0.120375535677152667, 0.32413434700070316,
|
||
0.00665482039381434);
|
||
ir->AddTriPoints6(114, 0.148066899157366667, 0.229277483555980969,
|
||
0.00707722325261307);
|
||
ir->AddTriPoints6(120, 0.191771865867325067, 0.325618122595983752,
|
||
0.007440689780584005);
|
||
ir->SetOrder(25);
|
||
TriangleIntRules[21] =
|
||
TriangleIntRules[22] =
|
||
TriangleIntRules[23] =
|
||
TriangleIntRules[24] =
|
||
TriangleIntRules[25] = ir;
|
||
return ir;
|
||
|
||
default:
|
||
// Grundmann-Moller fallback for orders beyond tabulated rules
|
||
int i = (Order / 2) * 2 + 1; // closest odd >= Order
|
||
AllocIntRule(TriangleIntRules, i);
|
||
ir = new IntegrationRule;
|
||
ir->GrundmannMollerSimplexRule(i/2, 2);
|
||
if (!TriangleIntRules[i-1]) { TriangleIntRules[i-1] = ir; }
|
||
TriangleIntRules[i] = ir;
|
||
return ir;
|
||
}
|
||
}
|
||
|
||
// Integration rules for unit square
|
||
IntegrationRule *IntegrationRules::SquareIntegrationRule(int Order)
|
||
{
|
||
int RealOrder = GetSegmentRealOrder(Order);
|
||
// Order is one of {RealOrder-1,RealOrder}
|
||
if (!HaveIntRule(SegmentIntRules, RealOrder))
|
||
{
|
||
SegmentIntegrationRule(RealOrder);
|
||
}
|
||
AllocIntRule(SquareIntRules, RealOrder); // RealOrder >= Order
|
||
SquareIntRules[RealOrder-1] =
|
||
SquareIntRules[RealOrder] =
|
||
new IntegrationRule(*SegmentIntRules[RealOrder],
|
||
*SegmentIntRules[RealOrder]);
|
||
return SquareIntRules[Order];
|
||
}
|
||
|
||
// Tet rules d=0-13 from Witherden & Vincent [1], orbit data from PyFR, CC-BY 4.0.
|
||
// Tet rules d=14-20 from Chuluunbaatar et al. [2], supplementary data.
|
||
IntegrationRule *IntegrationRules::TetrahedronIntegrationRule(int Order)
|
||
{
|
||
IntegrationRule *ir = NULL;
|
||
|
||
switch (Order)
|
||
{
|
||
case 0:
|
||
case 1:
|
||
ir = new IntegrationRule(1);
|
||
ir->AddTetMidPoint(0, 1./6.);
|
||
ir->SetOrder(1);
|
||
TetrahedronIntRules[0] =
|
||
TetrahedronIntRules[1] = ir;
|
||
return ir;
|
||
|
||
case 2:
|
||
ir = new IntegrationRule(4);
|
||
ir->AddTetPoints4(0, 1.38196601125010531952e-01, 1./24.);
|
||
ir->SetOrder(2);
|
||
TetrahedronIntRules[2] = ir;
|
||
return ir;
|
||
|
||
case 3:
|
||
ir = new IntegrationRule(8);
|
||
ir->AddTetPoints4(0, 3.28163302516381705232e-01, 2.27029737561812265667e-02);
|
||
ir->AddTetPoints4(4, 1.08047249898428621151e-01, 1.89636929104854412564e-02);
|
||
ir->SetOrder(3);
|
||
TetrahedronIntRules[3] = ir;
|
||
return ir;
|
||
|
||
case 4:
|
||
case 5:
|
||
ir = new IntegrationRule(14);
|
||
ir->AddTetPoints4(0, 3.10885919263300669613e-01, 1.87813209530026427319e-02);
|
||
ir->AddTetPoints4(4, 9.27352503108912484819e-02, 1.22488405193936587129e-02);
|
||
ir->AddTetPoints6(8, 4.54496295874350364485e-01, 7.09100346284691120807e-03);
|
||
ir->SetOrder(5);
|
||
TetrahedronIntRules[4] =
|
||
TetrahedronIntRules[5] = ir;
|
||
return ir;
|
||
|
||
case 6:
|
||
ir = new IntegrationRule(24);
|
||
ir->AddTetPoints4(0, 4.06739585346113652342e-02, 1.67953517588677390775e-03);
|
||
ir->AddTetPoints4(4, 3.22337890142275540484e-01, 9.22619692394245453915e-03);
|
||
ir->AddTetPoints4(8, 2.14602871259152117034e-01, 6.65379170969458179352e-03);
|
||
ir->AddTetPoints12(12, 6.36610018750174977420e-02, 6.03005664791649187428e-01,
|
||
8.03571428571428492127e-03);
|
||
ir->SetOrder(6);
|
||
TetrahedronIntRules[6] = ir;
|
||
return ir;
|
||
|
||
case 7:
|
||
ir = new IntegrationRule(35);
|
||
ir->AddTetMidPoint(0, 1.59142149106884754628e-02);
|
||
ir->AddTetPoints4(1, 3.15701149778202794227e-01, 7.05493020166117132397e-03);
|
||
ir->AddTetPoints6(5, 4.49510177401603649994e-01, 5.31615463880959638471e-03);
|
||
ir->AddTetPoints12(11, 1.88833831026001153219e-01, 5.75171637586999962011e-01,
|
||
6.20118845472243662709e-03);
|
||
ir->AddTetPoints12(23, 2.12654725414832546093e-02, 8.10830241098548620826e-01,
|
||
1.35179513831722359664e-03);
|
||
ir->SetOrder(7);
|
||
TetrahedronIntRules[7] = ir;
|
||
return ir;
|
||
|
||
case 8:
|
||
ir = new IntegrationRule(46);
|
||
ir->AddTetPoints4(0, 1.07952724962210866444e-01, 4.40444181806813866292e-03);
|
||
ir->AddTetPoints4(4, 1.85109487782586568105e-01, 8.67195792728975463348e-03);
|
||
ir->AddTetPoints4(8, 4.23165436847673381848e-02, 1.25420935892336655841e-03);
|
||
ir->AddTetPoints4(12, 3.14181709124039088010e-01, 6.96063047615581593358e-03);
|
||
ir->AddTetPoints6(16, 4.35591328583830206256e-01, 6.04682171021813687217e-03);
|
||
ir->AddTetPoints12(22, 2.14339301271305737728e-02, 7.17464063426308307214e-01,
|
||
1.19281714847407210867e-03);
|
||
ir->AddTetPoints12(34, 2.04139333876029116510e-01, 5.83797378302144398532e-01,
|
||
2.57558102516005586052e-03);
|
||
ir->SetOrder(8);
|
||
TetrahedronIntRules[8] = ir;
|
||
return ir;
|
||
|
||
case 9:
|
||
ir = new IntegrationRule(59);
|
||
ir->AddTetMidPoint(0, 9.66842481874670943431e-03);
|
||
ir->AddTetPoints4(1, 6.19817086544571793638e-10, 1.07198802932093984424e-05);
|
||
ir->AddTetPoints4(5, 1.60774535395261597426e-01, 3.86222307707090968185e-03);
|
||
ir->AddTetPoints4(9, 3.22276521821420969260e-01, 4.92715205590488116577e-03);
|
||
ir->AddTetPoints4(13, 4.51089183454135844720e-02, 1.34399666326936377374e-03);
|
||
ir->AddTetPoints6(17, 3.87703453995623947836e-01, 6.35568001728374458448e-03);
|
||
ir->AddTetPoints12(23, 4.58871448752459276665e-01, 7.97025232620401369310e-02,
|
||
1.39740369971642539558e-03);
|
||
ir->AddTetPoints12(35, 3.37758706853386048152e-02, 7.18350326442074527122e-01,
|
||
1.70575989212422133613e-03);
|
||
ir->AddTetPoints12(47, 1.83641369809927956780e-01, 5.98301349801968918030e-01,
|
||
3.42081932799802312939e-03);
|
||
ir->SetOrder(9);
|
||
TetrahedronIntRules[9] = ir;
|
||
return ir;
|
||
|
||
case 10:
|
||
ir = new IntegrationRule(81);
|
||
ir->AddTetMidPoint(0, 7.89996225933678984654e-03);
|
||
ir->AddTetPoints4(1, 3.12250068695188676138e-01, 4.48950999871145037950e-03);
|
||
ir->AddTetPoints4(5, 1.14309653857346149586e-01, 1.64485995279889710662e-03);
|
||
ir->AddTetPoints12(9, 4.10430739218965501269e-01, 1.65486025619611065718e-01,
|
||
1.89898020336587186781e-03);
|
||
ir->AddTetPoints12(21, 6.13800882479076381770e-03, 9.42988767345204870196e-01,
|
||
6.03240573898756009806e-05);
|
||
ir->AddTetPoints12(33, 1.21050181145589408338e-01, 4.77190379904280370660e-01,
|
||
4.28995533007601147213e-03);
|
||
ir->AddTetPoints12(45, 3.27794682164426753879e-02, 5.94256269480006982242e-01,
|
||
1.68931194662596552945e-03);
|
||
ir->AddTetPoints12(57, 3.24852815648231096901e-02, 8.01177284658344368573e-01,
|
||
1.09602454617265063913e-03);
|
||
ir->AddTetPoints12(69, 1.74979342183939068356e-01, 6.28071845475365986289e-01,
|
||
2.15117263314366490706e-03);
|
||
ir->SetOrder(10);
|
||
TetrahedronIntRules[10] = ir;
|
||
return ir;
|
||
|
||
case 11:
|
||
ir = new IntegrationRule(96);
|
||
ir->AddTetPoints4(0, 2.71527207067321363354e-02, 3.30755017786941475644e-04);
|
||
ir->AddTetPoints4(4, 7.29513610462571016058e-02, 1.27724462275054500421e-03);
|
||
ir->AddTetPoints4(8, 1.16306248902001030388e-01, 2.22195840281977797376e-03);
|
||
ir->AddTetPoints4(12, 1.79873804986097840519e-01, 3.55549424791121128006e-03);
|
||
ir->AddTetPoints4(16, 2.90224794862315171873e-01, 4.27767411104971236741e-03);
|
||
ir->AddTetPoints4(20, 3.25420936748619160639e-01, 2.29560465583227143668e-03);
|
||
ir->AddTetPoints6(24, 4.99998725049884129579e-01, 1.95152894059845476377e-04);
|
||
ir->AddTetPoints6(30, 3.94300142842090972639e-01, 4.13762713030314983886e-03);
|
||
ir->AddTetPoints12(36, 1.53994139264412854828e-02, 8.20202176629804657892e-01,
|
||
3.44113207868302869216e-04);
|
||
ir->AddTetPoints12(48, 4.36843254717693696421e-02, 6.27516751622257062948e-01,
|
||
2.00540889524405963051e-03);
|
||
ir->AddTetPoints12(60, 1.32316796082697751835e-01, 7.35366407834604496330e-01,
|
||
5.92160675031106853265e-04);
|
||
ir->AddTetPoints12(72, 2.14430354900043917965e-01, 5.31595425719235903372e-01,
|
||
3.02373698028425954079e-03);
|
||
ir->AddTetPoints12(84, 4.39586615093850330283e-01, 1.15789732843376125260e-01,
|
||
1.10416876556284249307e-03);
|
||
ir->SetOrder(11);
|
||
TetrahedronIntRules[11] = ir;
|
||
return ir;
|
||
|
||
case 12:
|
||
ir = new IntegrationRule(123);
|
||
ir->AddTetMidPoint(0, 3.73841522662751247000e-03);
|
||
ir->AddTetPoints4(1, 1.87550512633127830497e-02, 1.42695998696545141987e-04);
|
||
ir->AddTetPoints4(5, 1.08129536920462676619e-01, 2.21101382522646480733e-03);
|
||
ir->AddTetPoints4(9, 2.00131676822545012673e-01, 1.02716311841773611131e-03);
|
||
ir->AddTetPoints4(13, 3.00854293538076578152e-01, 3.84627572131096594557e-03);
|
||
ir->AddTetPoints4(17, 3.33333333333333259318e-01, 3.93263480259983290444e-04);
|
||
ir->AddTetPoints6(21, 4.61659950214442116323e-01, 3.64372815936188636666e-04);
|
||
ir->AddTetPoints12(27, 1.44707549187619299857e-02, 8.12825119403836504617e-01,
|
||
2.78228183082975693598e-04);
|
||
ir->AddTetPoints12(39, 1.93543398987769954545e-02, 6.01428742996530152354e-01,
|
||
5.44744358748572190913e-04);
|
||
ir->AddTetPoints12(51, 7.79277628532308863640e-02, 8.27642519452021385717e-01,
|
||
4.96603607085240776955e-04);
|
||
ir->AddTetPoints12(63, 1.22055870746741623734e-01, 4.77806520042316273944e-01,
|
||
3.50252598711278933380e-03);
|
||
ir->AddTetPoints12(75, 2.47870739372197945727e-01, 4.77761799116294072487e-01,
|
||
1.92951947508030146987e-03);
|
||
ir->AddTetPoints12(87, 4.29731509588804683197e-01, 1.17747435901101149547e-01,
|
||
1.59458000732484511328e-03);
|
||
ir->AddTetPoints24(99, 6.53037808968305988344e-01, 2.26776739658831050228e-01,
|
||
9.77349816032284102185e-02, 1.25441443948160597475e-03);
|
||
ir->SetOrder(12);
|
||
TetrahedronIntRules[12] = ir;
|
||
return ir;
|
||
|
||
case 13:
|
||
ir = new IntegrationRule(145);
|
||
ir->AddTetMidPoint(0, 4.65163625751287973520e-03);
|
||
ir->AddTetPoints4(1, 1.83047574861928130652e-02, 1.11328544338301012079e-04);
|
||
ir->AddTetPoints4(5, 1.79015082630022803745e-01, 3.21788310144650391634e-03);
|
||
ir->AddTetPoints4(9, 3.29615853754448240309e-01, 1.11613482625956810662e-03);
|
||
ir->AddTetPoints6(13, 4.84258919196047465938e-01, 4.22777636660235864308e-04);
|
||
ir->AddTetPoints6(19, 4.37693799377281589358e-01, 1.42741910394635984974e-03);
|
||
ir->AddTetPoints12(25, 1.68580250819503341120e-02, 7.16622722155387803511e-01,
|
||
3.27226375772531215529e-04);
|
||
ir->AddTetPoints12(37, 2.35257182448596058322e-02, 8.50600927605402956644e-01,
|
||
4.21131904215649394228e-04);
|
||
ir->AddTetPoints12(49, 6.94418950495464537553e-02, 7.08574314820604622689e-01,
|
||
9.98940835257684659268e-04);
|
||
ir->AddTetPoints12(61, 9.37005272821476720146e-02, 5.63456727822489344959e-01,
|
||
1.88623398086488818989e-03);
|
||
ir->AddTetPoints12(73, 1.25885360164042447995e-01, 7.40105985667665167149e-01,
|
||
5.36847192210583730974e-04);
|
||
ir->AddTetPoints12(85, 2.17955270547737667286e-01, 5.35290625276012344003e-01,
|
||
1.79931440889047528954e-03);
|
||
ir->AddTetPoints12(97, 3.54663455472783883948e-01, 2.00014910210617791186e-01,
|
||
2.91862732938839245650e-03);
|
||
ir->AddTetPoints12(109, 4.16689287657038387458e-01, 1.52054854976777675812e-01,
|
||
9.62750257012783515476e-04);
|
||
ir->AddTetPoints24(121, 6.06652560730350010054e-01, 3.04158433676372519372e-01,
|
||
7.79465622318078477093e-02, 6.21649861415869242447e-04);
|
||
ir->SetOrder(13);
|
||
TetrahedronIntRules[13] = ir;
|
||
return ir;
|
||
|
||
case 14:
|
||
// Chuluunbaatar et al. 2022: 175 pts, 1xCent + 6xS31 + 1xS22 + 10xS211 + 1xS1111
|
||
ir = new IntegrationRule(175);
|
||
ir->AddTetMidPoint(0, 2.79630622899013732072e-03);
|
||
ir->AddTetPoints4(1, 3.33328696010048830534e-01, 1.46917540892973303920e-04);
|
||
ir->AddTetPoints4(5, 2.03700979179134489261e-01, 1.62781576883158380503e-03);
|
||
ir->AddTetPoints4(9, 4.23119120487503441730e-02, 3.41825853298758448786e-04);
|
||
ir->AddTetPoints4(13, 1.66911321524259963212e-02, 6.86936532495300726303e-05);
|
||
ir->AddTetPoints4(17, 1.64429779556425403886e-01, 2.04589809259575743788e-03);
|
||
ir->AddTetPoints4(21, 3.05243130480787605574e-01, 2.91353212326472864671e-03);
|
||
ir->AddTetPoints6(25, 3.64287147870284933049e-01, 2.86344364876423311192e-03);
|
||
ir->AddTetPoints12(31, 2.15024351638664623643e-01, 5.09377822427890203372e-01,
|
||
1.86604702493200981690e-03);
|
||
ir->AddTetPoints12(43, 4.08453557824531576781e-01, 2.80291884809145824820e-02,
|
||
1.37079721642533879263e-03);
|
||
ir->AddTetPoints12(55, 2.59921479331125596102e-02, 7.58920443676681433232e-01,
|
||
5.79340182017268356431e-04);
|
||
ir->AddTetPoints12(67, 1.49228115767079897586e-02, 6.12483481308898292106e-01,
|
||
3.15838010189281473503e-04);
|
||
ir->AddTetPoints12(79, 8.37923462693035414617e-02, 8.22678981568125355928e-01,
|
||
3.37570057080462679680e-04);
|
||
ir->AddTetPoints12(91, 2.67966909251860618824e-01, 1.21450733299050331326e-02,
|
||
7.33968440516890967273e-04);
|
||
ir->AddTetPoints12(103, 7.22614743275835913483e-02, 2.93829968686904419162e-01,
|
||
1.73940500456261186446e-03);
|
||
ir->AddTetPoints12(115, 4.61690356122462508548e-01, 6.24312988244191680032e-02,
|
||
5.50156605537168688809e-04);
|
||
ir->AddTetPoints12(127, 1.02556247843651599492e-05, 9.09652220362116237240e-01,
|
||
3.74818592914694638193e-05);
|
||
ir->AddTetPoints12(139, 1.30667193397036723868e-01, 6.88760930866849863108e-01,
|
||
1.25403778742792012396e-03);
|
||
ir->AddTetPoints24(151, 4.96082264783182565887e-03, 1.13153535288820022986e-01,
|
||
2.55337379104889128367e-01, 5.28969173366363918341e-04);
|
||
ir->SetOrder(14);
|
||
TetrahedronIntRules[14] = ir;
|
||
return ir;
|
||
|
||
case 15:
|
||
// Chuluunbaatar et al. 2022: 209 pts, 1xCent + 4xS31 + 2xS22 + 11xS211 + 2xS1111
|
||
ir = new IntegrationRule(209);
|
||
ir->AddTetMidPoint(0, 1.41781886024826123995e-03);
|
||
ir->AddTetPoints4(1, 3.28314281102506377863e-01, 8.78039594754075579386e-04);
|
||
ir->AddTetPoints4(5, 5.95315181460130682378e-02, 6.49096346642090296988e-04);
|
||
ir->AddTetPoints4(9, 1.79953296856689010097e-01, 2.51315014419990637520e-03);
|
||
ir->AddTetPoints4(13, 2.87467267398706316506e-01, 1.49510722611222407301e-03);
|
||
ir->AddTetPoints6(17, 1.63820429539269674102e-01, 1.59662635597863640340e-03);
|
||
ir->AddTetPoints6(23, 4.49691037286174599696e-01, 5.56101982388976408267e-04);
|
||
ir->AddTetPoints12(29, 4.00736540413628217205e-01, 1.44084693588290586180e-02,
|
||
7.41098096825775301023e-04);
|
||
ir->AddTetPoints12(41, 9.30917130105696349895e-02, 5.00389099764625755462e-01,
|
||
1.48911776902807633377e-03);
|
||
ir->AddTetPoints12(53, 8.58969713610300000806e-02, 6.54216648304909331735e-01,
|
||
1.11933741359922559432e-03);
|
||
ir->AddTetPoints12(65, 2.10555592437809635520e-01, 6.11819979966393701076e-02,
|
||
1.62698522585349981094e-03);
|
||
ir->AddTetPoints12(77, 8.52418251507089524965e-02, 9.51322614055337932581e-03,
|
||
1.79269162929391381043e-04);
|
||
ir->AddTetPoints12(89, 3.39062766530538045595e-02, 1.00939010271196565223e-03,
|
||
6.94170071688395131987e-05);
|
||
ir->AddTetPoints12(101, 3.61923170322900333851e-01, 6.94744848139675630350e-02,
|
||
1.59447426536825733780e-03);
|
||
ir->AddTetPoints12(113, 4.82786943073602314858e-01, 2.98248878093377127463e-02,
|
||
1.85989656618890347554e-04);
|
||
ir->AddTetPoints12(125, 1.67482762532157707092e-02, 8.28143120169573809797e-01,
|
||
2.46974837616148544094e-04);
|
||
ir->AddTetPoints12(137, 1.77918303634979659000e-02, 2.95058306317452390122e-01,
|
||
3.64978786118624204133e-04);
|
||
ir->AddTetPoints12(149, 2.22052218944333024098e-01, 5.47335431979886655185e-01,
|
||
5.97708840203148013097e-04);
|
||
ir->AddTetPoints24(161, 7.01737933129022994905e-01, 1.90468590405707266511e-01,
|
||
9.00527571062145620884e-02, 6.09486799192013577881e-04);
|
||
ir->AddTetPoints24(185, 9.30456155647334665071e-02, 3.39714197260826189506e-01,
|
||
1.75628396157984228987e-02, 7.07458692200529210524e-04);
|
||
ir->SetOrder(15);
|
||
TetrahedronIntRules[15] = ir;
|
||
return ir;
|
||
|
||
case 16:
|
||
// Chuluunbaatar et al. 2022: 248 pts, 8xS31 + 2xS22 + 11xS211 + 3xS1111
|
||
ir = new IntegrationRule(248);
|
||
ir->AddTetPoints4(0, 3.27237393634992601577e-01, 1.02720766161859349518e-03);
|
||
ir->AddTetPoints4(4, 1.70006239733430930539e-01, 1.65526299995553852033e-03);
|
||
ir->AddTetPoints4(8, 1.15524594427552973475e-01, 9.38857005487496389280e-04);
|
||
ir->AddTetPoints4(12, 2.91444830780401391290e-02, 6.51307095609853886316e-05);
|
||
ir->AddTetPoints4(16, 2.99333264802760234957e-01, 1.67595204112568215392e-03);
|
||
ir->AddTetPoints4(20, 3.08156348381425804206e-01, 5.31861976759345615219e-04);
|
||
ir->AddTetPoints4(24, 1.50096498602994791322e-02, 5.70730781338430579624e-05);
|
||
ir->AddTetPoints4(28, 2.15377318942399170743e-01, 1.98341225672248695419e-03);
|
||
ir->AddTetPoints6(32, 4.31794349434656055120e-01, 1.33750640425309239717e-03);
|
||
ir->AddTetPoints6(38, 3.51744151127164061954e-01, 1.76194827177588689456e-03);
|
||
ir->AddTetPoints12(44, 8.29503118270854405969e-03, 1.02752005688054329213e-01,
|
||
7.68034210054447698257e-05);
|
||
ir->AddTetPoints12(56, 4.80211290069074772657e-02, 1.31255041819827861227e-01,
|
||
6.11232548742599826781e-04);
|
||
ir->AddTetPoints12(68, 1.58330800578366723275e-02, 7.34396180497698725098e-01,
|
||
2.28071918651308634196e-04);
|
||
ir->AddTetPoints12(80, 2.33860521976982954628e-01, 5.24417475091189966285e-01,
|
||
5.40957583749601829924e-04);
|
||
ir->AddTetPoints12(92, 4.03140399383019043533e-01, 1.61507733387133249614e-02,
|
||
7.92724342585451077596e-04);
|
||
ir->AddTetPoints12(104, 4.62664186044103586948e-01, 8.10608626417101511830e-03,
|
||
2.79154547643383486085e-04);
|
||
ir->AddTetPoints12(116, 9.54740058566225929804e-02, 2.29727696324949964835e-01,
|
||
9.08138512826854073234e-04);
|
||
ir->AddTetPoints12(128, 5.22644355194657739272e-02, 2.65023714835432855352e-01,
|
||
5.22825831653034723938e-04);
|
||
ir->AddTetPoints12(140, 1.54579380484822833525e-02, 3.89592797271225033118e-01,
|
||
2.70158998322436651619e-04);
|
||
ir->AddTetPoints12(152, 6.51049846146104782552e-02, 1.06703843590993385781e-02,
|
||
1.88285026959176838299e-04);
|
||
ir->AddTetPoints12(164, 1.54103896531334966236e-01, 6.61715834360067312048e-01,
|
||
9.04088475565851504463e-04);
|
||
ir->AddTetPoints24(176, 3.06401667521507548031e-01, 9.60517630472854377910e-02,
|
||
5.82737849734082380415e-01, 5.64378019946093099565e-04);
|
||
ir->AddTetPoints24(200, 7.52465510383990981991e-02, 7.45338308409307703783e-01,
|
||
1.68928693209220324653e-04, 1.50966565378719797755e-04);
|
||
ir->AddTetPoints24(224, 6.67856954025341370551e-02, 1.77429946613237937703e-01,
|
||
4.69596074566060506239e-01, 1.47055596469915315222e-03);
|
||
ir->SetOrder(16);
|
||
TetrahedronIntRules[16] = ir;
|
||
return ir;
|
||
|
||
case 17:
|
||
// Chuluunbaatar et al. 2022: 284 pts, 8xS31 + 2xS22 + 14xS211 + 3xS1111
|
||
ir = new IntegrationRule(284);
|
||
ir->AddTetPoints4(0, 7.70317217555786387662e-02, 5.61182432136912863994e-04);
|
||
ir->AddTetPoints4(4, 3.33178098937441047322e-01, 1.00689770519544758830e-04);
|
||
ir->AddTetPoints4(8, 4.70791056455278841830e-02, 3.67647177915098600219e-05);
|
||
ir->AddTetPoints4(12, 3.04818016813530989761e-01, 1.65889979279151768450e-03);
|
||
ir->AddTetPoints4(16, 1.30944391509640850613e-01, 1.44088986921434188127e-03);
|
||
ir->AddTetPoints4(20, 1.92535395691919936079e-01, 9.15395452245713805300e-04);
|
||
ir->AddTetPoints4(24, 2.76657577444746005657e-01, 1.54901732083429399985e-03);
|
||
ir->AddTetPoints4(28, 1.25812395975189866837e-02, 3.37712309004839736381e-05);
|
||
ir->AddTetPoints6(32, 3.79480026881957605706e-03, 6.09736907157372445766e-05);
|
||
ir->AddTetPoints6(38, 6.63574479091031538269e-02, 9.84476156008183339238e-04);
|
||
ir->AddTetPoints12(44, 1.55274078991054307469e-02, 1.62632006416904367763e-01,
|
||
1.72804764846827806904e-04);
|
||
ir->AddTetPoints12(56, 1.50826606630984655366e-01, 2.85203367985520928052e-01,
|
||
1.47139558752681400000e-03);
|
||
ir->AddTetPoints12(68, 2.15452134252546806392e-01, 4.98848586166224794436e-01,
|
||
1.39306363547645069810e-03);
|
||
ir->AddTetPoints12(80, 8.76098177043343750992e-02, 5.57666868254345748923e-01,
|
||
1.10287787958478406859e-03);
|
||
ir->AddTetPoints12(92, 2.76970146665180327883e-01, 4.27355049241211759625e-01,
|
||
7.68445912031412610432e-04);
|
||
ir->AddTetPoints12(104, 7.84854004483451911378e-02, 8.28808438721931994841e-01,
|
||
2.88969312852391376975e-04);
|
||
ir->AddTetPoints12(116, 4.16031236751370603333e-01, 5.86692256093802096822e-03,
|
||
3.64077087845785364682e-04);
|
||
ir->AddTetPoints12(128, 3.11418578536613735799e-03, 2.87446864291519776913e-01,
|
||
3.72979694270385682886e-05);
|
||
ir->AddTetPoints12(140, 1.48006973486492082737e-01, 3.23952556344107162056e-02,
|
||
8.00305827281323160088e-04);
|
||
ir->AddTetPoints12(152, 1.43033359605475689225e-02, 6.71799737186767331742e-02,
|
||
1.04156882366980408738e-04);
|
||
ir->AddTetPoints12(164, 4.67432600299733047589e-01, 1.06651363448138844503e-02,
|
||
3.28661937749209721524e-04);
|
||
ir->AddTetPoints12(176, 3.84385292133539946402e-01, 1.72438767641856061097e-01,
|
||
1.30629277435961670649e-03);
|
||
ir->AddTetPoints12(188, 4.91933575124020999736e-02, 1.88887052205218147760e-01,
|
||
6.32377286180225181393e-04);
|
||
ir->AddTetPoints12(200, 2.12339226453523544080e-01, 4.66228370030305223209e-03,
|
||
2.79290741478308349437e-04);
|
||
ir->AddTetPoints24(212, 5.56903978599597615506e-01, 3.03486206973905936479e-01,
|
||
1.18763501187465259079e-01, 6.89373805061415279201e-04);
|
||
ir->AddTetPoints24(236, 7.72626757757003540528e-02, 7.24630018034681633310e-01,
|
||
1.98007306189710574618e-01, 1.50128885728721429014e-04);
|
||
ir->AddTetPoints24(260, 3.94151894733046209707e-02, 1.37686205384732439361e-02,
|
||
3.27821577026260191356e-01, 2.69135394730690453036e-04);
|
||
ir->SetOrder(17);
|
||
TetrahedronIntRules[17] = ir;
|
||
return ir;
|
||
|
||
case 18:
|
||
// Chuluunbaatar et al. 2022: 343 pts, 1xCent + 6xS31 + 1xS22 + 18xS211 + 4xS1111
|
||
ir = new IntegrationRule(343);
|
||
ir->AddTetMidPoint(0, 1.50320520665968271855e-03);
|
||
ir->AddTetPoints4(1, 1.48031283019549930735e-01, 1.29183701010426432026e-03);
|
||
ir->AddTetPoints4(5, 9.18424577295562372115e-02, 5.65165600131921114398e-04);
|
||
ir->AddTetPoints4(9, 1.21731006846268821620e-02, 2.96638661291863846498e-05);
|
||
ir->AddTetPoints4(13, 2.96287086243479214076e-01, 1.68021599975901742008e-03);
|
||
ir->AddTetPoints4(17, 3.26360945420223702573e-01, 5.68549421501223146459e-04);
|
||
ir->AddTetPoints4(21, 2.16789137320780644913e-01, 5.55087845531530756776e-04);
|
||
ir->AddTetPoints6(25, 4.02614199568341046831e-01, 9.23257682834909085973e-04);
|
||
ir->AddTetPoints12(31, 4.40000294606430919497e-01, 2.67115252330815747261e-02,
|
||
5.89236226881148018354e-04);
|
||
ir->AddTetPoints12(43, 4.12210360146149590310e-01, 1.74700018752273367184e-01,
|
||
2.20201645262970173086e-04);
|
||
ir->AddTetPoints12(55, 3.72351275734651154803e-01, 2.01036053559202482210e-01,
|
||
1.19114067751767893286e-03);
|
||
ir->AddTetPoints12(67, 2.69242251358920825499e-01, 4.49869807112305730712e-01,
|
||
4.91465288515259167076e-04);
|
||
ir->AddTetPoints12(79, 9.76627002277863226487e-02, 5.25608170470769464622e-01,
|
||
9.64486286268097844226e-04);
|
||
ir->AddTetPoints12(91, 7.95369840699704060138e-03, 9.21057037109250575924e-01,
|
||
3.50539800669945179434e-05);
|
||
ir->AddTetPoints12(103, 1.87267495112264620305e-01, 6.22654876406694146596e-01,
|
||
2.16974411183503681708e-04);
|
||
ir->AddTetPoints12(115, 1.11072676172302167719e-01, 7.68063140393386190041e-01,
|
||
2.34282910136374296871e-04);
|
||
ir->AddTetPoints12(127, 7.28479245819699250397e-02, 6.49722949416110195919e-01,
|
||
4.87753175028903700802e-04);
|
||
ir->AddTetPoints12(139, 4.35242030819264283381e-02, 1.17645467635705727738e-01,
|
||
3.42108745586317604826e-04);
|
||
ir->AddTetPoints12(151, 4.75773867838972852606e-01, 4.07808063375027368691e-02,
|
||
1.54370387724311262630e-04);
|
||
ir->AddTetPoints12(163, 1.34493684207502642303e-02, 7.10473230757141527292e-01,
|
||
1.53084245505138950502e-04);
|
||
ir->AddTetPoints12(175, 1.55188126953707594691e-01, 4.80062353992693063853e-02,
|
||
7.41762550201790354931e-04);
|
||
ir->AddTetPoints12(187, 9.79636346189126545891e-03, 3.98937047199835026490e-01,
|
||
9.70189072253606300464e-05);
|
||
ir->AddTetPoints12(199, 2.22972555180978554423e-01, 4.76498622927509052349e-01,
|
||
1.23606250864963739498e-03);
|
||
ir->AddTetPoints12(211, 4.21260289724883496554e-02, 5.74400719147259319897e-01,
|
||
5.42935360416797459064e-04);
|
||
ir->AddTetPoints12(223, 5.33405089760383491204e-02, 8.80640211270697026436e-01,
|
||
1.37609067910292992955e-04);
|
||
ir->AddTetPoints12(235, 1.61516953402295604381e-01, 3.97825734805866804145e-01,
|
||
1.38360407936932868454e-03);
|
||
ir->AddTetPoints24(247, 1.45407574632878761056e-01, 2.91164813093958863011e-01,
|
||
2.48113931665566757323e-02, 7.08660104525119741853e-04);
|
||
ir->AddTetPoints24(271, 3.05950727633923398596e-03, 8.24555339410965038027e-01,
|
||
2.68565028804892240444e-02, 7.44949831353088352858e-05);
|
||
ir->AddTetPoints24(295, 6.97926614756394281258e-01, 2.13661027507205059095e-01,
|
||
1.66604565123493510159e-02, 3.48774332060316084436e-04);
|
||
ir->AddTetPoints24(319, 7.44110286479428006956e-02, 3.30529031579615995007e-01,
|
||
5.94802303330229431566e-01, 1.27738536486342427233e-04);
|
||
ir->SetOrder(18);
|
||
TetrahedronIntRules[18] = ir;
|
||
return ir;
|
||
|
||
case 19:
|
||
// Chuluunbaatar et al. 2022: 383 pts, 1xCent + 7xS31 + 3xS22 + 18xS211 + 5xS1111
|
||
ir = new IntegrationRule(383);
|
||
ir->AddTetMidPoint(0, 1.63415516118113374362e-03);
|
||
ir->AddTetPoints4(1, 1.99329047511150186933e-01, 1.40358969707281050661e-03);
|
||
ir->AddTetPoints4(5, 3.19645815662209620278e-01, 7.59091701642911130359e-04);
|
||
ir->AddTetPoints4(9, 1.28431938718745111000e-02, 3.39837309564992975973e-05);
|
||
ir->AddTetPoints4(13, 4.48982322308715264825e-02, 2.64426237344559213801e-04);
|
||
ir->AddTetPoints4(17, 1.36956642483832574664e-01, 7.02856031175181001670e-04);
|
||
ir->AddTetPoints4(21, 2.90023140379987831583e-01, 1.36565197140369625796e-03);
|
||
ir->AddTetPoints4(25, 9.74727471850612287030e-02, 5.49967641751253539552e-04);
|
||
ir->AddTetPoints6(29, 9.71505706534206842084e-02, 8.76433367036966496331e-04);
|
||
ir->AddTetPoints6(35, 3.43201766155373844125e-01, 1.42791869569969145752e-03);
|
||
ir->AddTetPoints6(41, 4.81048682529186311108e-01, 2.14104659346860782621e-04);
|
||
ir->AddTetPoints12(47, 1.65183633185142647593e-01, 6.20343271948313068620e-01,
|
||
7.66624121398570487762e-04);
|
||
ir->AddTetPoints12(59, 6.34076230480007801971e-02, 1.00712056526788096511e-06,
|
||
5.83669782740338187726e-05);
|
||
ir->AddTetPoints12(71, 1.14358224224068777061e-03, 4.00863727471189756901e-01,
|
||
2.15179213304474853136e-05);
|
||
ir->AddTetPoints12(83, 1.72841898844304682481e-02, 2.80974353964125900252e-01,
|
||
2.11660243173845957895e-04);
|
||
ir->AddTetPoints12(95, 3.72064110603054443160e-01, 1.99705629942815821032e-01,
|
||
8.75485553265252145448e-04);
|
||
ir->AddTetPoints12(107, 5.20242736384851911513e-02, 1.37531998164526964024e-01,
|
||
2.41875978548309640689e-04);
|
||
ir->AddTetPoints12(119, 1.60391938988630439189e-01, 6.93856268828658799552e-03,
|
||
2.41421983084743403854e-04);
|
||
ir->AddTetPoints12(131, 2.04009311788170112981e-03, 1.85612033720361502276e-01,
|
||
2.06414634761428247161e-05);
|
||
ir->AddTetPoints12(143, 5.27148535043930888122e-02, 3.66275984765305206992e-01,
|
||
5.92721175866269185152e-04);
|
||
ir->AddTetPoints12(155, 2.43716736043213727525e-01, 4.92302726125094514131e-01,
|
||
5.87805221231469288319e-04);
|
||
ir->AddTetPoints12(167, 2.35953750535251915998e-01, 4.36147999140300668408e-01,
|
||
1.13870552683716466137e-03);
|
||
ir->AddTetPoints12(179, 3.76396525093455058819e-01, 2.39487568074092604942e-01,
|
||
3.83776606890474675932e-04);
|
||
ir->AddTetPoints12(191, 1.01345934656195998946e-01, 7.75724008450901503231e-01,
|
||
3.13399374751491786861e-04);
|
||
ir->AddTetPoints12(203, 1.03978576485678226443e-02, 9.15876371248755760668e-01,
|
||
5.00269118618190970684e-05);
|
||
ir->AddTetPoints12(215, 4.34825905512194077485e-01, 1.54748703679767566493e-02,
|
||
5.11025086732051357814e-04);
|
||
ir->AddTetPoints12(227, 1.21014927057761692564e-01, 5.16611295841134299245e-01,
|
||
1.01828426898903845119e-03);
|
||
ir->AddTetPoints12(239, 2.02221519510974792611e-02, 1.24848861415461573343e-01,
|
||
1.43839758373707761420e-04);
|
||
ir->AddTetPoints12(251, 6.45349753497565792326e-02, 6.58797428082353198064e-01,
|
||
5.65395165437599170333e-04);
|
||
ir->AddTetPoints24(263, 5.41751356566992220420e-02, 1.84132753560833639650e-01,
|
||
7.55314661874476711567e-01, 1.46482060260151955187e-04);
|
||
ir->AddTetPoints24(287, 8.81662316742434920558e-02, 6.34269181940578685719e-01,
|
||
1.71630199679625762565e-02, 3.35574461050504140036e-04);
|
||
ir->AddTetPoints24(311, 2.99876324799488391815e-01, 5.12011805461436986242e-01,
|
||
1.39945756017624101109e-01, 6.87105878169750922298e-04);
|
||
ir->AddTetPoints24(335, 5.43033048543508201078e-01, 3.98941499659789464149e-03,
|
||
2.96197429831241032527e-01, 1.95466436126630728823e-04);
|
||
ir->AddTetPoints24(359, 5.64748804926985426000e-01, 5.04233375578317308263e-02,
|
||
5.38328049657370907161e-03, 1.64230458281612438400e-04);
|
||
ir->SetOrder(19);
|
||
TetrahedronIntRules[19] = ir;
|
||
return ir;
|
||
|
||
case 20:
|
||
// Chuluunbaatar et al. 2022: 441 pts, 1xCent + 8xS31 + 4xS22 + 20xS211 + 6xS1111
|
||
ir = new IntegrationRule(441);
|
||
ir->AddTetMidPoint(0, 1.18189152746071531389e-03);
|
||
ir->AddTetPoints4(1, 1.44398440418483348102e-01, 9.66793543666131616025e-04);
|
||
ir->AddTetPoints4(5, 7.58537944731913719304e-03, 7.60496594911488061041e-06);
|
||
ir->AddTetPoints4(9, 2.92814880923072839991e-01, 1.15590090717039508002e-03);
|
||
ir->AddTetPoints4(13, 3.21283351882928780441e-01, 6.88636083704610908220e-04);
|
||
ir->AddTetPoints4(17, 1.99126561548209762842e-01, 1.13432796616464034133e-03);
|
||
ir->AddTetPoints4(21, 9.94395843777090698845e-02, 4.86016284904959409725e-04);
|
||
ir->AddTetPoints4(25, 5.64411054542164752901e-02, 3.22610846546138882538e-04);
|
||
ir->AddTetPoints4(29, 2.29655958568571322287e-02, 5.81298999656466752642e-05);
|
||
ir->AddTetPoints6(33, 1.79141319969986889671e-01, 1.09529243625682828553e-03);
|
||
ir->AddTetPoints6(39, 1.28302591207222288494e-01, 5.73999808941768881014e-04);
|
||
ir->AddTetPoints6(45, 1.45390267831490144212e-02, 1.87025723866918939901e-04);
|
||
ir->AddTetPoints6(51, 4.21353125034043096697e-01, 7.00971843927932197066e-04);
|
||
ir->AddTetPoints12(57, 2.31009838674045342444e-01, 4.40089079468850674637e-01,
|
||
9.76079610667369892280e-04);
|
||
ir->AddTetPoints12(69, 5.34968850305569260106e-03, 1.93845423004036843118e-01,
|
||
3.05362840831512903964e-05);
|
||
ir->AddTetPoints12(81, 1.24450062776632008887e-01, 2.68040213994591769442e-01,
|
||
8.86819612544235340128e-04);
|
||
ir->AddTetPoints12(93, 2.04325810970697151203e-02, 2.58308114426872181824e-01,
|
||
1.92969809694827680807e-04);
|
||
ir->AddTetPoints12(105, 4.82828231821251577238e-02, 5.80396926577266859815e-03,
|
||
6.19875387085781813859e-05);
|
||
ir->AddTetPoints12(117, 2.89149787325270579696e-01, 4.18198665244226719384e-01,
|
||
1.83280536557154684341e-04);
|
||
ir->AddTetPoints12(129, 4.46052961749180063022e-02, 1.63432451082933805075e-01,
|
||
3.33781952353497124181e-04);
|
||
ir->AddTetPoints12(141, 2.70746696244779198881e-03, 9.47473449934223443947e-01,
|
||
9.76902870004440518315e-06);
|
||
ir->AddTetPoints12(153, 1.81599434536722420530e-01, 5.60679256816125279328e-02,
|
||
7.16163033436471277611e-04);
|
||
ir->AddTetPoints12(165, 7.98241847716316815786e-02, 2.24797220241435974364e-01,
|
||
6.67198343081130593700e-04);
|
||
ir->AddTetPoints12(177, 1.46751424308700709198e-02, 1.05996716787954137207e-01,
|
||
1.02112432519835052438e-04);
|
||
ir->AddTetPoints12(189, 1.47988480867858707146e-01, 6.98383249617804735543e-01,
|
||
1.65331156535892409174e-04);
|
||
ir->AddTetPoints12(201, 4.52421884404454466289e-01, 8.04532263339534647884e-02,
|
||
3.75572061072812404120e-04);
|
||
ir->AddTetPoints12(213, 4.82988925439242506449e-03, 3.56408995372303527560e-01,
|
||
3.70372995907359520520e-05);
|
||
ir->AddTetPoints12(225, 3.96075105354866952023e-01, 1.93643489260524576112e-01,
|
||
4.20355142354041121707e-04);
|
||
ir->AddTetPoints12(237, 2.11473197416018027228e-01, 5.66485319771568907044e-01,
|
||
3.20737452814357003415e-04);
|
||
ir->AddTetPoints12(249, 2.56953469781508958558e-01, 4.54914546979866607490e-01,
|
||
5.44699807144207316482e-04);
|
||
ir->AddTetPoints12(261, 3.63762446007509787638e-01, 7.07704644682126682298e-02,
|
||
8.71724711789707298014e-04);
|
||
ir->AddTetPoints12(273, 4.34084693566413395982e-02, 5.61826818484091217165e-01,
|
||
4.66716966955154613419e-04);
|
||
ir->AddTetPoints12(285, 1.26096840063810999855e-01, 4.15826803180920218095e-02,
|
||
4.05374493922667372432e-04);
|
||
ir->AddTetPoints24(297, 3.19830665436792060952e-01, 4.36032333155288651105e-02,
|
||
5.02224495466759290885e-01, 7.04309544520175120040e-04);
|
||
ir->AddTetPoints24(321, 6.93104765295092647634e-04, 7.61588104432530443866e-01,
|
||
4.88127934219945436300e-02, 6.22004488693067759562e-05);
|
||
ir->AddTetPoints24(345, 3.01059509765821443905e-03, 4.68001522754562040984e-02,
|
||
5.92923309754523120141e-01, 1.11516752988944994277e-04);
|
||
ir->AddTetPoints24(369, 8.05046770737637640281e-01, 1.15008082676269607347e-01,
|
||
6.71446262920421810261e-02, 1.37946350993749402127e-04);
|
||
ir->AddTetPoints24(393, 3.27837344763098725853e-01, 1.41845459805815643506e-01,
|
||
3.73982332962941683638e-03, 1.75610557117139728066e-04);
|
||
ir->AddTetPoints24(417, 1.69431918115453827856e-02, 9.30651836894259287813e-02,
|
||
6.46141331822991715761e-01, 3.76832469454361519961e-04);
|
||
ir->SetOrder(20);
|
||
TetrahedronIntRules[20] = ir;
|
||
return ir;
|
||
|
||
default:
|
||
// Grundmann-Moller fallback for orders beyond tabulated rules
|
||
int i = (Order / 2) * 2 + 1; // closest odd >= Order
|
||
AllocIntRule(TetrahedronIntRules, i);
|
||
ir = new IntegrationRule;
|
||
ir->GrundmannMollerSimplexRule(i/2, 3);
|
||
if (!TetrahedronIntRules[i-1]) { TetrahedronIntRules[i-1] = ir; }
|
||
TetrahedronIntRules[i] = ir;
|
||
return ir;
|
||
}
|
||
}
|
||
|
||
// Integration rules for reference pyramid
|
||
IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
|
||
{
|
||
// This is a simple integration rule adapted from an integration
|
||
// rule for a cube which seems to be adequate for now. We should continue
|
||
// to search for a more appropriate integration rule designed specifically
|
||
// for pyramid elements.
|
||
const IntegrationRule &irc = Get(Geometry::CUBE, Order);
|
||
int npts = irc.GetNPoints();
|
||
AllocIntRule(PyramidIntRules, Order);
|
||
PyramidIntRules[Order] = new IntegrationRule(npts);
|
||
PyramidIntRules[Order]->SetOrder(Order);
|
||
|
||
if (npts == 1)
|
||
{
|
||
// We handle this as a special case because with only one integration
|
||
// point we cannot accurately integrate the quadratic factor
|
||
// pow(1.0 - ipc.z, 2) and the resulting weight does not match the volume
|
||
// of the reference element.
|
||
IntegrationPoint &ipp = PyramidIntRules[Order]->IntPoint(0);
|
||
ipp.x = 0.375;
|
||
ipp.y = 0.375;
|
||
ipp.z = 0.25;
|
||
ipp.weight = 1.0 / 3.0;
|
||
}
|
||
else
|
||
{
|
||
for (int k=0; k<npts; k++)
|
||
{
|
||
const IntegrationPoint &ipc = irc.IntPoint(k);
|
||
IntegrationPoint &ipp = PyramidIntRules[Order]->IntPoint(k);
|
||
ipp.x = ipc.x * (1.0 - ipc.z);
|
||
ipp.y = ipc.y * (1.0 - ipc.z);
|
||
ipp.z = ipc.z;
|
||
ipp.weight = ipc.weight * pow(1.0 - ipc.z, 2);
|
||
}
|
||
}
|
||
return PyramidIntRules[Order];
|
||
}
|
||
|
||
// Integration rules for reference prism
|
||
IntegrationRule *IntegrationRules::PrismIntegrationRule(int Order)
|
||
{
|
||
const IntegrationRule &irt = Get(Geometry::TRIANGLE, Order);
|
||
const IntegrationRule &irs = Get(Geometry::SEGMENT, Order);
|
||
int nt = irt.GetNPoints();
|
||
int ns = irs.GetNPoints();
|
||
AllocIntRule(PrismIntRules, Order);
|
||
PrismIntRules[Order] = new IntegrationRule(nt * ns);
|
||
PrismIntRules[Order]->SetOrder(std::min(irt.GetOrder(), irs.GetOrder()));
|
||
while (Order < std::min(irt.GetOrder(), irs.GetOrder()))
|
||
{
|
||
AllocIntRule(PrismIntRules, ++Order);
|
||
PrismIntRules[Order] = PrismIntRules[Order-1];
|
||
}
|
||
|
||
for (int ks=0; ks<ns; ks++)
|
||
{
|
||
const IntegrationPoint &ips = irs.IntPoint(ks);
|
||
for (int kt=0; kt<nt; kt++)
|
||
{
|
||
int kp = ks * nt + kt;
|
||
const IntegrationPoint &ipt = irt.IntPoint(kt);
|
||
IntegrationPoint &ipp = PrismIntRules[Order]->IntPoint(kp);
|
||
ipp.x = ipt.x;
|
||
ipp.y = ipt.y;
|
||
ipp.z = ips.x;
|
||
ipp.weight = ipt.weight * ips.weight;
|
||
}
|
||
}
|
||
return PrismIntRules[Order];
|
||
}
|
||
|
||
// Integration rules for reference cube
|
||
IntegrationRule *IntegrationRules::CubeIntegrationRule(int Order)
|
||
{
|
||
int RealOrder = GetSegmentRealOrder(Order);
|
||
if (!HaveIntRule(SegmentIntRules, RealOrder))
|
||
{
|
||
SegmentIntegrationRule(RealOrder);
|
||
}
|
||
AllocIntRule(CubeIntRules, RealOrder);
|
||
CubeIntRules[RealOrder-1] =
|
||
CubeIntRules[RealOrder] =
|
||
new IntegrationRule(*SegmentIntRules[RealOrder],
|
||
*SegmentIntRules[RealOrder],
|
||
*SegmentIntRules[RealOrder]);
|
||
return CubeIntRules[Order];
|
||
}
|
||
|
||
StroudIntegrationRules StroudIntRules;
|
||
|
||
StroudIntegrationRules::StroudIntegrationRules()
|
||
{
|
||
const MemoryType h_mt = MemoryType::HOST;
|
||
SquareStroudIntRules.SetSize(32, h_mt);
|
||
SquareStroudIntRules = NULL;
|
||
|
||
TriangleStroudIntRules.SetSize(32, h_mt);
|
||
TriangleStroudIntRules = NULL;
|
||
|
||
CubeStroudIntRules.SetSize(32, h_mt);
|
||
CubeStroudIntRules = NULL;
|
||
|
||
TetrahedronStroudIntRules.SetSize(32, h_mt);
|
||
TetrahedronStroudIntRules = NULL;
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
IntRuleLocks.SetSize(Geometry::NUM_GEOMETRIES, h_mt);
|
||
for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
|
||
{
|
||
omp_init_lock(&IntRuleLocks[i]);
|
||
}
|
||
#endif
|
||
}
|
||
|
||
const IntegrationRule &StroudIntegrationRules::Get(int GeomType, int Order)
|
||
{
|
||
Array<IntegrationRule *> *ir_array = NULL;
|
||
|
||
switch (GeomType)
|
||
{
|
||
case Geometry::TRIANGLE: ir_array = &TriangleStroudIntRules; break;
|
||
case Geometry::TETRAHEDRON: ir_array = &TetrahedronStroudIntRules; break;
|
||
case Geometry::INVALID:
|
||
case Geometry::NUM_GEOMETRIES:
|
||
MFEM_ABORT("Unknown type of reference element!");
|
||
default:
|
||
MFEM_ABORT("Stroud rules only valid for triangular and tetrahedral elements!");
|
||
}
|
||
|
||
if (Order < 0)
|
||
{
|
||
Order = 0;
|
||
}
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_set_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
|
||
if (!HaveIntRule(*ir_array, Order))
|
||
{
|
||
IntegrationRule *ir = GenerateIntegrationRule(GeomType, Order);
|
||
#ifdef MFEM_DEBUG
|
||
int RealOrder = Order;
|
||
while (RealOrder+1 < ir_array->Size() && (*ir_array)[RealOrder+1] == ir)
|
||
{
|
||
RealOrder++;
|
||
}
|
||
MFEM_VERIFY(RealOrder == ir->GetOrder(), "internal error");
|
||
#else
|
||
MFEM_CONTRACT_VAR(ir);
|
||
#endif
|
||
}
|
||
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
omp_unset_lock(&IntRuleLocks[GeomType]);
|
||
#endif
|
||
|
||
return *(*ir_array)[Order];
|
||
}
|
||
|
||
void StroudIntegrationRules::DeleteIntRuleArray(
|
||
Array<IntegrationRule *> &ir_array) const
|
||
{
|
||
// Many of the intrules have multiple contiguous copies in the ir_array
|
||
// so we have to be careful to not delete them twice.
|
||
IntegrationRule *ir = NULL;
|
||
for (int i = 0; i < ir_array.Size(); i++)
|
||
{
|
||
if (ir_array[i] != NULL && ir_array[i] != ir)
|
||
{
|
||
ir = ir_array[i];
|
||
delete ir;
|
||
}
|
||
}
|
||
}
|
||
|
||
StroudIntegrationRules::~StroudIntegrationRules()
|
||
{
|
||
#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
|
||
for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
|
||
{
|
||
omp_destroy_lock(&IntRuleLocks[i]);
|
||
}
|
||
#endif
|
||
DeleteIntRuleArray(SquareStroudIntRules);
|
||
DeleteIntRuleArray(TriangleStroudIntRules);
|
||
DeleteIntRuleArray(CubeStroudIntRules);
|
||
DeleteIntRuleArray(TetrahedronStroudIntRules);
|
||
}
|
||
|
||
|
||
IntegrationRule *StroudIntegrationRules::GenerateIntegrationRule(int GeomType,
|
||
int Order)
|
||
{
|
||
switch (GeomType)
|
||
{
|
||
case Geometry::TRIANGLE:
|
||
return TriangleStroudIntegrationRule(Order);
|
||
case Geometry::TETRAHEDRON:
|
||
return TetrahedronStroudIntegrationRule(Order);
|
||
case Geometry::INVALID:
|
||
case Geometry::NUM_GEOMETRIES:
|
||
MFEM_ABORT("Unknown type of reference element!");
|
||
default:
|
||
MFEM_ABORT("Stroud rules only valid for triangular and tetrahedral elements!");
|
||
}
|
||
return NULL;
|
||
}
|
||
|
||
/* Integration rule in reference triangle according to tensor product Gauss-Jacobi rule.
|
||
The nodes and weights are used in the original form defined on the reference
|
||
square to evaluate the component 1D basis functions. Mapping to the reference
|
||
triangle via IntegrationRule::DuffyTrans() occurs only in evaluation of coefficient
|
||
vectors, see e.g. MassIntegrator::AssemblePASimplex. */
|
||
IntegrationRule *StroudIntegrationRules::TriangleStroudIntegrationRule(
|
||
int Order)
|
||
{
|
||
int RealOrder = GetSegmentRealOrder(Order);
|
||
// Order is one of {RealOrder-1,RealOrder}
|
||
// if (!HaveIntRule(SegmentIntRules, RealOrder))
|
||
// {
|
||
// SegmentIntegrationRule(RealOrder);
|
||
// }
|
||
IntegrationRule ir_0_0;
|
||
// Gauss-Jacobi is exact for 2*n-1
|
||
int n = RealOrder/2 + 1;
|
||
QuadratureFunctions1D::GaussJacobi(n, 0.0, 0.0, &ir_0_0);
|
||
|
||
IntegrationRule ir_1_0;
|
||
QuadratureFunctions1D::GaussJacobi(n, 1.0, 0.0, &ir_1_0);
|
||
|
||
AllocIntRule(TriangleStroudIntRules, RealOrder); // RealOrder >= Order
|
||
// create rule in unit square
|
||
TriangleStroudIntRules[RealOrder-1] =
|
||
TriangleStroudIntRules[RealOrder] =
|
||
new IntegrationRule(ir_1_0, ir_0_0);
|
||
// map rule to reference triangle
|
||
// TriangleStroudIntRules[RealOrder-1]->DuffyTrans(2);
|
||
*TriangleStroudIntRules[RealOrder-1] =
|
||
DuffyTrans(*TriangleStroudIntRules[RealOrder-1], 2);
|
||
return TriangleStroudIntRules[Order];
|
||
}
|
||
|
||
/* Integration rule in reference tetrahedron according to tensor product Gauss-Jacobi rule.
|
||
The nodes and weights are used in the original form defined on the reference
|
||
square to evaluate the component 1D basis functions. Mapping to the reference
|
||
triangle via IntegrationRule::DuffyTrans() occurs only in evaluation of coefficient
|
||
vectors, see e.g. MassIntegrator::AssemblePASimplex. */
|
||
IntegrationRule *StroudIntegrationRules::TetrahedronStroudIntegrationRule(
|
||
int Order)
|
||
{
|
||
int RealOrder = GetSegmentRealOrder(Order);
|
||
// Order is one of {RealOrder-1,RealOrder}
|
||
|
||
IntegrationRule ir_0_0;
|
||
int n = RealOrder/2 + 1;
|
||
QuadratureFunctions1D::GaussJacobi(n, 0.0, 0.0, &ir_0_0);
|
||
|
||
IntegrationRule ir_1_0;
|
||
QuadratureFunctions1D::GaussJacobi(n, 1.0, 0.0, &ir_1_0);
|
||
|
||
IntegrationRule ir_2_0;
|
||
QuadratureFunctions1D::GaussJacobi(n, 2.0, 0.0, &ir_2_0);
|
||
|
||
AllocIntRule(TetrahedronStroudIntRules, RealOrder); // RealOrder >= Order
|
||
// create rule in unit cube
|
||
TetrahedronStroudIntRules[RealOrder-1] =
|
||
TetrahedronStroudIntRules[RealOrder] =
|
||
new IntegrationRule(ir_2_0, ir_1_0, ir_0_0);
|
||
// map rule to reference tetrahedron
|
||
// TetrahedronStroudIntRules[RealOrder-1]->DuffyTrans(3);
|
||
*TetrahedronStroudIntRules[RealOrder-1] =
|
||
DuffyTrans(*TetrahedronStroudIntRules[RealOrder-1], 3);
|
||
return TetrahedronStroudIntRules[Order];
|
||
}
|
||
|
||
IntegrationRule& NURBSMeshRules::GetElementRule(const int elem,
|
||
const int patch, const int *ijk,
|
||
Array<const KnotVector*> const& kv) const
|
||
{
|
||
// First check whether a rule has been assigned to element index elem.
|
||
auto search = elementToRule.find(elem);
|
||
if (search != elementToRule.end())
|
||
{
|
||
return *elementRule[search->second];
|
||
}
|
||
|
||
#ifndef MFEM_THREAD_SAFE
|
||
// If no prescribed rule is given for the current element, a temporary one is
|
||
// formed by restricting a tensor-product of 1D rules to the element. The
|
||
// ownership model for this temporary rule is not thread-safe.
|
||
|
||
MFEM_VERIFY(patchRules1D.NumRows(),
|
||
"Undefined rule in NURBSMeshRules::GetElementRule");
|
||
|
||
// Use a tensor product of rules on the patch.
|
||
MFEM_VERIFY(kv.Size() == dim, "");
|
||
|
||
int np = 1;
|
||
std::vector<std::vector<real_t>> el(dim);
|
||
|
||
std::vector<int> npd;
|
||
npd.assign(3, 0);
|
||
|
||
for (int d=0; d<dim; ++d)
|
||
{
|
||
const int order = kv[d]->GetOrder();
|
||
|
||
const real_t kv0 = (*kv[d])[order + ijk[d]];
|
||
const real_t kv1 = (*kv[d])[order + ijk[d] + 1];
|
||
|
||
const bool rightEnd = (order + ijk[d] + 1) == (kv[d]->Size() - 1);
|
||
|
||
for (int i=0; i<patchRules1D(patch,d)->Size(); ++i)
|
||
{
|
||
const IntegrationPoint& ip = (*patchRules1D(patch,d))[i];
|
||
if (kv0 <= ip.x && (ip.x < kv1 || rightEnd))
|
||
{
|
||
const real_t x = (ip.x - kv0) / (kv1 - kv0);
|
||
el[d].push_back(x);
|
||
el[d].push_back(ip.weight);
|
||
}
|
||
}
|
||
|
||
npd[d] = static_cast<int>(el[d].size() / 2);
|
||
np *= npd[d];
|
||
}
|
||
|
||
temporaryElementRule.SetSize(np);
|
||
|
||
// Set temporaryElementRule[i + j*npd[0] + k*npd[0]*npd[1]] =
|
||
// (el[0][2*i], el[1][2*j], el[2][2*k])
|
||
|
||
MFEM_VERIFY(npd[0] > 0 && npd[1] > 0, "Assuming 2D or 3D");
|
||
|
||
for (int i = 0; i < npd[0]; ++i)
|
||
{
|
||
for (int j = 0; j < npd[1]; ++j)
|
||
{
|
||
for (int k = 0; k < std::max(npd[2], 1); ++k)
|
||
{
|
||
const int id = i + j*npd[0] + k*npd[0]*npd[1];
|
||
temporaryElementRule[id].x = el[0][2*i];
|
||
temporaryElementRule[id].y = el[1][2*j];
|
||
|
||
temporaryElementRule[id].weight = el[0][(2*i)+1];
|
||
temporaryElementRule[id].weight *= el[1][(2*j)+1];
|
||
|
||
if (npd[2] > 0)
|
||
{
|
||
temporaryElementRule[id].z = el[2][2*k];
|
||
temporaryElementRule[id].weight *= el[2][(2*k)+1];
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
return temporaryElementRule;
|
||
#else
|
||
MFEM_ABORT("Temporary integration rules on NURBS elements "
|
||
"are not thread-safe.");
|
||
#endif
|
||
}
|
||
|
||
void NURBSMeshRules::GetIntegrationPointFrom1D(const int patch, int i, int j,
|
||
int k, IntegrationPoint & ip)
|
||
{
|
||
MFEM_VERIFY(patchRules1D.NumRows() > 0,
|
||
"Assuming patchRules1D is set.");
|
||
|
||
ip.weight = (*patchRules1D(patch,0))[i].weight;
|
||
ip.x = (*patchRules1D(patch,0))[i].x;
|
||
|
||
if (dim > 1)
|
||
{
|
||
ip.weight *= (*patchRules1D(patch,1))[j].weight;
|
||
ip.y = (*patchRules1D(patch,1))[j].x; // 1D rule only has x
|
||
}
|
||
|
||
if (dim > 2)
|
||
{
|
||
ip.weight *= (*patchRules1D(patch,2))[k].weight;
|
||
ip.z = (*patchRules1D(patch,2))[k].x; // 1D rule only has x
|
||
}
|
||
}
|
||
|
||
void NURBSMeshRules::Finalize(Mesh const& mesh)
|
||
{
|
||
if ((int) pointToElem.size() == npatches) { return; } // Already set
|
||
|
||
MFEM_VERIFY(elementToRule.empty() && patchRules1D.NumRows() > 0
|
||
&& npatches > 0, "Assuming patchRules1D is set.");
|
||
MFEM_VERIFY(mesh.NURBSext, "");
|
||
MFEM_VERIFY(mesh.Dimension() == dim, "");
|
||
|
||
pointToElem.resize(npatches);
|
||
patchRules1D_KnotSpan.resize(npatches);
|
||
|
||
// First, find all the elements in each patch.
|
||
std::vector<std::vector<int>> patchElements(npatches);
|
||
|
||
for (int e=0; e<mesh.GetNE(); ++e)
|
||
{
|
||
patchElements[mesh.NURBSext->GetElementPatch(e)].push_back(e);
|
||
}
|
||
|
||
Array<int> ijk(3);
|
||
Array<int> maxijk(3);
|
||
Array<int> np(3); // Number of points in each dimension
|
||
ijk = 0;
|
||
|
||
Array<const KnotVector*> pkv;
|
||
|
||
for (int p=0; p<npatches; ++p)
|
||
{
|
||
patchRules1D_KnotSpan[p].resize(dim);
|
||
|
||
// For each patch, get the range of ijk.
|
||
mesh.NURBSext->GetPatchKnotVectors(p, pkv);
|
||
MFEM_VERIFY((int) pkv.Size() == dim, "");
|
||
|
||
maxijk = 1;
|
||
np = 1;
|
||
for (int d=0; d<dim; ++d)
|
||
{
|
||
maxijk[d] = pkv[d]->GetNKS();
|
||
np[d] = patchRules1D(p,d)->Size();
|
||
}
|
||
|
||
// For each patch, set a map from ijk to element index.
|
||
Array3D<int> ijk2elem(maxijk[0], maxijk[1], maxijk[2]);
|
||
ijk2elem = -1;
|
||
|
||
for (auto elem : patchElements[p])
|
||
{
|
||
mesh.NURBSext->GetElementIJK(elem, ijk);
|
||
MFEM_VERIFY(ijk2elem(ijk[0], ijk[1], ijk[2]) == -1, "");
|
||
ijk2elem(ijk[0], ijk[1], ijk[2]) = elem;
|
||
}
|
||
|
||
// For each point, find its ijk and from that its element index.
|
||
// It is assumed here that the NURBSFiniteElement kv the same as the
|
||
// patch kv.
|
||
|
||
for (int d=0; d<dim; ++d)
|
||
{
|
||
patchRules1D_KnotSpan[p][d].SetSize(patchRules1D(p,d)->Size());
|
||
|
||
for (int r=0; r<patchRules1D(p,d)->Size(); ++r)
|
||
{
|
||
const IntegrationPoint& ip = (*patchRules1D(p,d))[r];
|
||
|
||
const int order = pkv[d]->GetOrder();
|
||
|
||
// Find ijk_d such that ip.x is in the corresponding knot-span.
|
||
int ijk_d = 0;
|
||
bool found = false;
|
||
while (!found)
|
||
{
|
||
const real_t kv0 = (*pkv[d])[order + ijk_d];
|
||
const real_t kv1 = (*pkv[d])[order + ijk_d + 1];
|
||
|
||
const bool rightEnd = (order + ijk_d + 1) == (pkv[d]->Size() - 1);
|
||
|
||
if (kv0 <= ip.x && (ip.x < kv1 || rightEnd))
|
||
{
|
||
found = true;
|
||
}
|
||
else
|
||
{
|
||
ijk_d++;
|
||
}
|
||
}
|
||
|
||
patchRules1D_KnotSpan[p][d][r] = ijk_d;
|
||
}
|
||
}
|
||
|
||
pointToElem[p].SetSize(np[0], np[1], np[2]);
|
||
for (int i=0; i<np[0]; ++i)
|
||
for (int j=0; j<np[1]; ++j)
|
||
for (int k=0; k<np[2]; ++k)
|
||
{
|
||
const int elem = ijk2elem(patchRules1D_KnotSpan[p][0][i],
|
||
patchRules1D_KnotSpan[p][1][j],
|
||
patchRules1D_KnotSpan[p][2][k]);
|
||
MFEM_VERIFY(elem >= 0, "");
|
||
pointToElem[p](i,j,k) = elem;
|
||
}
|
||
} // Loop (p) over patches
|
||
}
|
||
|
||
void NURBSMeshRules::SetPatchRules1D(const int patch,
|
||
std::vector<const IntegrationRule*> & ir1D)
|
||
{
|
||
MFEM_VERIFY((int) ir1D.size() == dim, "Wrong dimension");
|
||
|
||
for (int i=0; i<dim; ++i)
|
||
{
|
||
patchRules1D(patch,i) = ir1D[i];
|
||
}
|
||
}
|
||
|
||
NURBSMeshRules::~NURBSMeshRules()
|
||
{
|
||
for (int i=0; i<patchRules1D.NumRows(); ++i)
|
||
for (int j=0; j<patchRules1D.NumCols(); ++j)
|
||
{
|
||
delete patchRules1D(i, j);
|
||
}
|
||
}
|
||
|
||
}
|