diff --git a/fem/fe/fe_h1.cpp b/fem/fe/fe_h1.cpp index bf1936de16..227551fcad 100644 --- a/fem/fe/fe_h1.cpp +++ b/fem/fe/fe_h1.cpp @@ -1243,12 +1243,12 @@ void H1_FuentesPyramidElement::calcBasis(const int p, int o = 0; // Vertices - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { u[0] = u[1] = u[2] = u[3] = 0.25 * (1.0 - z); u[4] = z; - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y; @@ -1439,7 +1439,7 @@ void H1_FuentesPyramidElement::calcGradBasis(const int p, int o = 0; // Vertices - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { du(0,0) = -0.5; du(0,1) = -0.5; du(0,2) = -0.75; du(1,0) = 0.5; du(1,1) = -0.5; du(1,2) = -0.25; @@ -1447,7 +1447,7 @@ void H1_FuentesPyramidElement::calcGradBasis(const int p, du(3,0) = -0.5; du(3,1) = 0.5; du(3,2) = -0.25; du(4,0) = 0.0; du(4,1) = 0.0; du(4,2) = 1.0; - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y; diff --git a/fem/fe/fe_nd.cpp b/fem/fe/fe_nd.cpp index b1945ef652..116fdc3af9 100644 --- a/fem/fe/fe_nd.cpp +++ b/fem/fe/fe_nd.cpp @@ -1904,9 +1904,9 @@ void ND_FuentesPyramidElement::calcBasis(const int p, Vector xy({x,y}), dmu(3); real_t mu, mu2; - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y; @@ -2170,9 +2170,9 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p, real_t mu, mu2; - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y; diff --git a/fem/fe/fe_pyramid.hpp b/fem/fe/fe_pyramid.hpp index 74156769b6..a4ed471804 100644 --- a/fem/fe/fe_pyramid.hpp +++ b/fem/fe/fe_pyramid.hpp @@ -58,6 +58,9 @@ private: mutable DenseMatrix V_R_mtmp; #endif +protected: + const real_t apex_tol = 1e-8; + public: FuentesPyramid() = default; diff --git a/fem/fe/fe_rt.cpp b/fem/fe/fe_rt.cpp index 4b7a2ffe5e..887bf97d2d 100644 --- a/fem/fe/fe_rt.cpp +++ b/fem/fe/fe_rt.cpp @@ -1521,9 +1521,9 @@ void RT_FuentesPyramidElement::calcBasis(const int p, Vector xy({x,y}); real_t mu; - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y; @@ -1782,10 +1782,10 @@ void RT_FuentesPyramidElement::calcDivBasis(const int p, real_t mu; bool limz1 = false; - if (std::fabs(1.0 - z) < 1e-8) + if (std::fabs(1.0 - z) < apex_tol) { limz1 = true; - z = 1.0 - 1e-8; + z = 1.0 - apex_tol; y = 0.5 * (1.0 - z); x = 0.5 * (1.0 - z); xy(0) = x; xy(1) = y;