phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Generalized floating point type. So far, ex1 works for a serial build without lapack.

Browse Source
This commit is contained in:
Dylan Copeland
2023-10-09 14:19:25 -07:00
parent c0867ca009
commit f106c03dd1
149 changed files with 6236 additions and 6227 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/fe/fe_base.cpp
+86 -86
View File
@@ -510,7 +510,7 @@ void ScalarFiniteElement::NodalLocalInterpolation(
ElementTransformation &Trans, DenseMatrix &I,
const ScalarFiniteElement &fine_fe) const
{
double v[Geometry::MaxDim];
fptype v[Geometry::MaxDim];
Vector vv(v, dim);
IntegrationPoint f_ip;
@@ -551,7 +551,7 @@ void ScalarFiniteElement::ScalarLocalInterpolation(
{
// General "interpolation", defined by L2 projection
double v[Geometry::MaxDim];
fptype v[Geometry::MaxDim];
Vector vv(v, dim);
IntegrationPoint f_ip;
@@ -591,7 +591,7 @@ void ScalarFiniteElement::ScalarLocalL2Restriction(
const ScalarFiniteElement &coarse_fe) const
{
// General "restriction", defined by L2 projection
double v[Geometry::MaxDim];
fptype v[Geometry::MaxDim];
Vector vv(v, dim);
const int cs = coarse_fe.GetDof(), fs = this->GetDof();
@@ -646,7 +646,7 @@ void NodalFiniteElement::ProjectCurl_2D(
{
fe.CalcCurlShape(Nodes.IntPoint(i), curl_shape);
double w = 1.0;
fptype w = 1.0;
if (GetMapType() == FiniteElement::VALUE)
{
Trans.SetIntPoint(&Nodes.IntPoint(i));
@@ -667,7 +667,7 @@ void InvertLinearTrans(ElementTransformation &trans,
p0.Set3(0, 0, 0);
trans.Transform(p0, x);
double store[3];
fptype store[3];
Vector v(store, x.Size());
pt.Get(store, x.Size());
v -= x;
@@ -860,7 +860,7 @@ void NodalFiniteElement::ProjectDiv(
const FiniteElement &fe, ElementTransformation &Trans,
DenseMatrix &div) const
{
double detJ;
fptype detJ;
Vector div_shape(fe.GetDof());
div.SetSize(dof, fe.GetDof());
@@ -978,10 +978,10 @@ void VectorFiniteElement::CalcVShape_ND(
}
void VectorFiniteElement::Project_RT(
const double *nk, const Array<int> &d2n,
const fptype *nk, const Array<int> &d2n,
VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
{
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
const int sdim = Trans.GetSpaceDim();
MFEM_ASSERT(vc.GetVDim() == sdim, "");
Vector xk(vk, sdim);
@@ -998,7 +998,7 @@ void VectorFiniteElement::Project_RT(
}
void VectorFiniteElement::Project_RT(
const double *nk, const Array<int> &d2n,
const fptype *nk, const Array<int> &d2n,
Vector &vc, ElementTransformation &Trans, Vector &dofs) const
{
const int sdim = Trans.GetSpaceDim();
@@ -1015,7 +1015,7 @@ void VectorFiniteElement::Project_RT(
}
void VectorFiniteElement::ProjectMatrixCoefficient_RT(
const double *nk, const Array<int> &d2n,
const fptype *nk, const Array<int> &d2n,
MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
{
// project the rows of the matrix coefficient in an RT space
@@ -1043,12 +1043,12 @@ void VectorFiniteElement::ProjectMatrixCoefficient_RT(
}
void VectorFiniteElement::Project_RT(
const double *nk, const Array<int> &d2n, const FiniteElement &fe,
const fptype *nk, const Array<int> &d2n, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &I) const
{
if (fe.GetRangeType() == SCALAR)
{
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
Vector shape(fe.GetDof());
int sdim = Trans.GetSpaceDim();
@@ -1064,7 +1064,7 @@ void VectorFiniteElement::Project_RT(
Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, vk);
if (fe.GetMapType() == INTEGRAL)
{
double w = 1.0/Trans.Weight();
fptype w = 1.0/Trans.Weight();
for (int d = 0; d < dim; d++)
{
vk[d] *= w;
@@ -1073,7 +1073,7 @@ void VectorFiniteElement::Project_RT(
for (int j = 0; j < shape.Size(); j++)
{
double s = shape(j);
fptype s = shape(j);
if (fabs(s) < 1e-12)
{
s = 0.0;
@@ -1090,7 +1090,7 @@ void VectorFiniteElement::Project_RT(
else
{
int sdim = Trans.GetSpaceDim();
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
DenseMatrix vshape(fe.GetDof(), sdim);
Vector vshapenk(fe.GetDof());
const bool square_J = (dim == sdim);
@@ -1118,7 +1118,7 @@ void VectorFiniteElement::Project_RT(
}
void VectorFiniteElement::ProjectGrad_RT(
const double *nk, const Array<int> &d2n, const FiniteElement &fe,
const fptype *nk, const Array<int> &d2n, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &grad) const
{
if (dim != 2)
@@ -1128,7 +1128,7 @@ void VectorFiniteElement::ProjectGrad_RT(
DenseMatrix dshape(fe.GetDof(), fe.GetDim());
Vector grad_k(fe.GetDof());
double tk[2];
fptype tk[2];
grad.SetSize(dof, fe.GetDof());
for (int k = 0; k < dof; k++)
@@ -1145,7 +1145,7 @@ void VectorFiniteElement::ProjectGrad_RT(
}
void VectorFiniteElement::ProjectCurl_ND(
const double *tk, const Array<int> &d2t, const FiniteElement &fe,
const fptype *tk, const Array<int> &d2t, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &curl) const
{
#ifdef MFEM_THREAD_SAFE
@@ -1183,7 +1183,7 @@ void VectorFiniteElement::ProjectCurl_ND(
}
void VectorFiniteElement::ProjectCurl_RT(
const double *nk, const Array<int> &d2n, const FiniteElement &fe,
const fptype *nk, const Array<int> &d2n, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &curl) const
{
DenseMatrix curl_shape(fe.GetDof(), dim);
@@ -1202,10 +1202,10 @@ void VectorFiniteElement::ProjectCurl_RT(
}
void VectorFiniteElement::Project_ND(
const double *tk, const Array<int> &d2t,
const fptype *tk, const Array<int> &d2t,
VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
{
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
Vector xk(vk, vc.GetVDim());
for (int k = 0; k < dof; k++)
@@ -1219,7 +1219,7 @@ void VectorFiniteElement::Project_ND(
}
void VectorFiniteElement::Project_ND(
const double *tk, const Array<int> &d2t,
const fptype *tk, const Array<int> &d2t,
Vector &vc, ElementTransformation &Trans, Vector &dofs) const
{
for (int k = 0; k < dof; k++)
@@ -1231,7 +1231,7 @@ void VectorFiniteElement::Project_ND(
}
void VectorFiniteElement::ProjectMatrixCoefficient_ND(
const double *tk, const Array<int> &d2t,
const fptype *tk, const Array<int> &d2t,
MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
{
// project the rows of the matrix coefficient in an ND space
@@ -1257,13 +1257,13 @@ void VectorFiniteElement::ProjectMatrixCoefficient_ND(
}
void VectorFiniteElement::Project_ND(
const double *tk, const Array<int> &d2t, const FiniteElement &fe,
const fptype *tk, const Array<int> &d2t, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &I) const
{
if (fe.GetRangeType() == SCALAR)
{
int sdim = Trans.GetSpaceDim();
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
Vector shape(fe.GetDof());
I.SetSize(dof, sdim*fe.GetDof());
@@ -1278,7 +1278,7 @@ void VectorFiniteElement::Project_ND(
Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
if (fe.GetMapType() == INTEGRAL)
{
double w = 1.0/Trans.Weight();
fptype w = 1.0/Trans.Weight();
for (int d = 0; d < sdim; d++)
{
vk[d] *= w;
@@ -1287,7 +1287,7 @@ void VectorFiniteElement::Project_ND(
for (int j = 0; j < shape.Size(); j++)
{
double s = shape(j);
fptype s = shape(j);
if (fabs(s) < 1e-12)
{
s = 0.0;
@@ -1304,7 +1304,7 @@ void VectorFiniteElement::Project_ND(
else
{
int sdim = Trans.GetSpaceDim();
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
DenseMatrix vshape(fe.GetDof(), sdim);
Vector vshapetk(fe.GetDof());
@@ -1330,7 +1330,7 @@ void VectorFiniteElement::Project_ND(
}
void VectorFiniteElement::ProjectGrad_ND(
const double *tk, const Array<int> &d2t, const FiniteElement &fe,
const fptype *tk, const Array<int> &d2t, const FiniteElement &fe,
ElementTransformation &Trans, DenseMatrix &grad) const
{
MFEM_ASSERT(fe.GetMapType() == VALUE, "");
@@ -1370,7 +1370,7 @@ void VectorFiniteElement::LocalL2Projection_RT(
for (int i = 0; i < ir.GetNPoints(); i++)
{
const IntegrationPoint &ip = ir.IntPoint(i);
double w = ip.weight;
fptype w = ip.weight;
this->CalcVShape(ip, fine_shape);
Trans.Transform(ip, v);
tr_ip.Set(v.GetData(), dim);
@@ -1381,7 +1381,7 @@ void VectorFiniteElement::LocalL2Projection_RT(
{
for (int j=0; j<cs; ++j)
{
double Mkj = 0.0;
fptype Mkj = 0.0;
for (int d1=0; d1<dim; ++d1)
{
for (int d2=0; d2<dim; ++d2)
@@ -1398,14 +1398,14 @@ void VectorFiniteElement::LocalL2Projection_RT(
}
void VectorFiniteElement::LocalInterpolation_RT(
const VectorFiniteElement &cfe, const double *nk, const Array<int> &d2n,
const VectorFiniteElement &cfe, const fptype *nk, const Array<int> &d2n,
ElementTransformation &Trans, DenseMatrix &I) const
{
MFEM_ASSERT(map_type == cfe.GetMapType(), "");
if (!is_nodal) { return LocalL2Projection_RT(cfe, Trans, I); }
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
Vector xk(vk, dim);
IntegrationPoint ip;
#ifdef MFEM_THREAD_SAFE
@@ -1428,7 +1428,7 @@ void VectorFiniteElement::LocalInterpolation_RT(
// I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
for (int j = 0; j < vshape.Height(); j++)
{
double Ikj = 0.;
fptype Ikj = 0.;
for (int i = 0; i < dim; i++)
{
Ikj += vshape(j, i) * vk[i];
@@ -1468,7 +1468,7 @@ void VectorFiniteElement::LocalL2Projection_ND(
{
for (int j=0; j<cs; ++j)
{
double Mkj = 0.0;
fptype Mkj = 0.0;
for (int d1=0; d1<dim; ++d1)
{
for (int d2=0; d2<dim; ++d2)
@@ -1485,12 +1485,12 @@ void VectorFiniteElement::LocalL2Projection_ND(
}
void VectorFiniteElement::LocalInterpolation_ND(
const VectorFiniteElement &cfe, const double *tk, const Array<int> &d2t,
const VectorFiniteElement &cfe, const fptype *tk, const Array<int> &d2t,
ElementTransformation &Trans, DenseMatrix &I) const
{
if (!is_nodal) { return LocalL2Projection_ND(cfe, Trans, I); }
double vk[Geometry::MaxDim];
fptype vk[Geometry::MaxDim];
Vector xk(vk, dim);
IntegrationPoint ip;
#ifdef MFEM_THREAD_SAFE
@@ -1513,7 +1513,7 @@ void VectorFiniteElement::LocalInterpolation_ND(
// I_k = vshape_k.J.t_k, k=1,...,Dof
for (int j = 0; j < vshape.Height(); j++)
{
double Ikj = 0.;
fptype Ikj = 0.;
for (int i = 0; i < dim; i++)
{
Ikj += vshape(j, i) * vk[i];
@@ -1524,10 +1524,10 @@ void VectorFiniteElement::LocalInterpolation_ND(
}
void VectorFiniteElement::LocalRestriction_RT(
const double *nk, const Array<int> &d2n, ElementTransformation &Trans,
const fptype *nk, const Array<int> &d2n, ElementTransformation &Trans,
DenseMatrix &R) const
{
double pt_data[Geometry::MaxDim];
fptype pt_data[Geometry::MaxDim];
IntegrationPoint ip;
Vector pt(pt_data, dim);
@@ -1537,7 +1537,7 @@ void VectorFiniteElement::LocalRestriction_RT(
Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
const DenseMatrix &J = Trans.Jacobian();
const double weight = Trans.Weight();
const fptype weight = Trans.Weight();
for (int j = 0; j < dof; j++)
{
InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
@@ -1549,7 +1549,7 @@ void VectorFiniteElement::LocalRestriction_RT(
pt /= weight;
for (int k = 0; k < dof; k++)
{
double R_jk = 0.0;
fptype R_jk = 0.0;
for (int d = 0; d < dim; d++)
{
R_jk += vshape(k,d)*pt_data[d];
@@ -1567,10 +1567,10 @@ void VectorFiniteElement::LocalRestriction_RT(
}
void VectorFiniteElement::LocalRestriction_ND(
const double *tk, const Array<int> &d2t, ElementTransformation &Trans,
const fptype *tk, const Array<int> &d2t, ElementTransformation &Trans,
DenseMatrix &R) const
{
double pt_data[Geometry::MaxDim];
fptype pt_data[Geometry::MaxDim];
IntegrationPoint ip;
Vector pt(pt_data, dim);
@@ -1590,7 +1590,7 @@ void VectorFiniteElement::LocalRestriction_ND(
Jinv.Mult(tk+dim*d2t[j], pt_data);
for (int k = 0; k < dof; k++)
{
double R_jk = 0.0;
fptype R_jk = 0.0;
for (int d = 0; d < dim; d++)
{
R_jk += vshape(k,d)*pt_data[d];
@@ -1608,7 +1608,7 @@ void VectorFiniteElement::LocalRestriction_ND(
}
Poly_1D::Basis::Basis(const int p, const double *nodes, EvalType etype)
Poly_1D::Basis::Basis(const int p, const fptype *nodes, EvalType etype)
: etype(etype), auxiliary_basis(NULL), scale_integrated(false)
{
switch (etype)
@@ -1636,7 +1636,7 @@ Poly_1D::Basis::Basis(const int p, const double *nodes, EvalType etype)
{
for (int j = 0; j < i; j++)
{
double xij = x(i) - x(j);
fptype xij = x(i) - x(j);
w(i) *= xij;
w(j) *= -xij;
}
@@ -1672,7 +1672,7 @@ Poly_1D::Basis::Basis(const int p, const double *nodes, EvalType etype)
}
}
void Poly_1D::Basis::Eval(const double y, Vector &u) const
void Poly_1D::Basis::Eval(const fptype y, Vector &u) const
{
switch (etype)
{
@@ -1685,7 +1685,7 @@ void Poly_1D::Basis::Eval(const double y, Vector &u) const
case Barycentric:
{
int i, k, p = x.Size() - 1;
double l, lk;
fptype l, lk;
if (p == 0)
{
@@ -1733,7 +1733,7 @@ void Poly_1D::Basis::Eval(const double y, Vector &u) const
}
}
void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d) const
void Poly_1D::Basis::Eval(const fptype y, Vector &u, Vector &d) const
{
switch (etype)
{
@@ -1747,7 +1747,7 @@ void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d) const
case Barycentric:
{
int i, k, p = x.Size() - 1;
double l, lp, lk, sk, si;
fptype l, lp, lk, sk, si;
if (p == 0)
{
@@ -1813,7 +1813,7 @@ void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d) const
}
}
void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d,
void Poly_1D::Basis::Eval(const fptype y, Vector &u, Vector &d,
Vector &d2) const
{
MFEM_VERIFY(etype == Barycentric,
@@ -1832,7 +1832,7 @@ void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d,
case Barycentric:
{
int i, k, p = x.Size() - 1;
double l, lp, lp2, lk, sk, si, sk2;
fptype l, lp, lp2, lk, sk, si, sk2;
if (p == 0)
{
@@ -1956,19 +1956,19 @@ const int *Poly_1D::Binom(const int p)
return binom[p];
}
void Poly_1D::ChebyshevPoints(const int p, double *x)
void Poly_1D::ChebyshevPoints(const int p, fptype *x)
{
for (int i = 0; i <= p; i++)
{
// x[i] = 0.5*(1. + cos(M_PI*(p - i + 0.5)/(p + 1)));
double s = sin(M_PI_2*(i + 0.5)/(p + 1));
fptype s = sin(M_PI_2*(i + 0.5)/(p + 1));
x[i] = s*s;
}
}
void Poly_1D::CalcMono(const int p, const double x, double *u)
void Poly_1D::CalcMono(const int p, const fptype x, fptype *u)
{
double xn;
fptype xn;
u[0] = xn = 1.;
for (int n = 1; n <= p; n++)
{
@@ -1976,9 +1976,9 @@ void Poly_1D::CalcMono(const int p, const double x, double *u)
}
}
void Poly_1D::CalcMono(const int p, const double x, double *u, double *d)
void Poly_1D::CalcMono(const int p, const fptype x, fptype *u, fptype *d)
{
double xn;
fptype xn;
u[0] = xn = 1.;
d[0] = 0.;
for (int n = 1; n <= p; n++)
@@ -1988,8 +1988,8 @@ void Poly_1D::CalcMono(const int p, const double x, double *u, double *d)
}
}
void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
double *u)
void Poly_1D::CalcBinomTerms(const int p, const fptype x, const fptype y,
fptype *u)
{
if (p == 0)
{
@@ -1999,7 +1999,7 @@ void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
{
int i;
const int *b = Binom(p);
double z = x;
fptype z = x;
for (i = 1; i < p; i++)
{
@@ -2017,8 +2017,8 @@ void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
}
}
void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
double *u, double *d)
void Poly_1D::CalcBinomTerms(const int p, const fptype x, const fptype y,
fptype *u, fptype *d)
{
if (p == 0)
{
@@ -2029,8 +2029,8 @@ void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
{
int i;
const int *b = Binom(p);
const double xpy = x + y, ptx = p*x;
double z = 1.;
const fptype xpy = x + y, ptx = p*x;
fptype z = 1.;
for (i = 1; i < p; i++)
{
@@ -2052,8 +2052,8 @@ void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
}
}
void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
double *d)
void Poly_1D::CalcDBinomTerms(const int p, const fptype x, const fptype y,
fptype *d)
{
if (p == 0)
{
@@ -2063,8 +2063,8 @@ void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
{
int i;
const int *b = Binom(p);
const double xpy = x + y, ptx = p*x;
double z = 1.;
const fptype xpy = x + y, ptx = p*x;
fptype z = 1.;
for (i = 1; i < p; i++)
{
@@ -2082,11 +2082,11 @@ void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
}
}
void Poly_1D::CalcLegendre(const int p, const double x, double *u)
void Poly_1D::CalcLegendre(const int p, const fptype x, fptype *u)
{
// use the recursive definition for [-1,1]:
// (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
double z;
fptype z;
u[0] = 1.;
if (p == 0) { return; }
u[1] = z = 2.*x - 1.;
@@ -2096,13 +2096,13 @@ void Poly_1D::CalcLegendre(const int p, const double x, double *u)
}
}
void Poly_1D::CalcLegendre(const int p, const double x, double *u, double *d)
void Poly_1D::CalcLegendre(const int p, const fptype x, fptype *u, fptype *d)
{
// use the recursive definition for [-1,1]:
// (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
// for the derivative use, z in [-1,1]:
// P'_{n+1}(z) = (2*n+1)*P_n(z)+P'_{n-1}(z)
double z;
fptype z;
u[0] = 1.;
d[0] = 0.;
if (p == 0) { return; }
@@ -2115,12 +2115,12 @@ void Poly_1D::CalcLegendre(const int p, const double x, double *u, double *d)
}
}
void Poly_1D::CalcChebyshev(const int p, const double x, double *u)
void Poly_1D::CalcChebyshev(const int p, const fptype x, fptype *u)
{
// recursive definition, z in [-1,1]
// T_0(z) = 1, T_1(z) = z
// T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
double z;
fptype z;
u[0] = 1.;
if (p == 0) { return; }
u[1] = z = 2.*x - 1.;
@@ -2130,7 +2130,7 @@ void Poly_1D::CalcChebyshev(const int p, const double x, double *u)
}
}
void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d)
void Poly_1D::CalcChebyshev(const int p, const fptype x, fptype *u, fptype *d)
{
// recursive definition, z in [-1,1]
// T_0(z) = 1, T_1(z) = z
@@ -2140,7 +2140,7 @@ void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d)
// U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
// U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
// T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
double z;
fptype z;
u[0] = 1.;
d[0] = 0.;
if (p == 0) { return; }
@@ -2153,8 +2153,8 @@ void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d)
}
}
void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d,
double *dd)
void Poly_1D::CalcChebyshev(const int p, const fptype x, fptype *u, fptype *d,
fptype *dd)
{
// recursive definition, z in [-1,1]
// T_0(z) = 1, T_1(z) = z
@@ -2165,7 +2165,7 @@ void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d,
// U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
// T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
// T''_{n+1}(z) = (n + 1)*(2*(n + 1)*T'_n(z) + z*T''_n(z)) / n
double z;
fptype z;
u[0] = 1.;
d[0] = 0.;
dd[0]= 0.;
@@ -2181,7 +2181,7 @@ void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d,
}
}
const double *Poly_1D::GetPoints(const int p, const int btype)
const fptype *Poly_1D::GetPoints(const int p, const int btype)
{
BasisType::Check(btype);
const int qtype = BasisType::GetQuadrature1D(btype);
@@ -2190,16 +2190,16 @@ const double *Poly_1D::GetPoints(const int p, const int btype)
if (points_container.find(btype) == points_container.end())
{
points_container[btype] = new Array<double*>(h_mt);
points_container[btype] = new Array<fptype*>(h_mt);
}
Array<double*> &pts = *points_container[btype];
Array<fptype*> &pts = *points_container[btype];
if (pts.Size() <= p)
{
pts.SetSize(p + 1, NULL);
}
if (pts[p] == NULL)
{
pts[p] = new double[p + 1];
pts[p] = new fptype[p + 1];
quad_func.GivePolyPoints(p+1, pts[p], qtype);
}
return pts[p];
@@ -2235,7 +2235,7 @@ Poly_1D::~Poly_1D()
for (PointsMap::iterator it = points_container.begin();
it != points_container.end() ; ++it)
{
Array<double*>& pts = *it->second;
Array<fptype*>& pts = *it->second;
for ( int i = 0 ; i < pts.Size() ; ++i )
{
delete [] pts[i];