diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
index 608b56b9946c3c25f2e4f7c0cd9b4cf951a7d775..77e79851e0a8000807f6901373c8801c26130bb2 100644
--- a/src/algebra/IntegrationTools.hpp
+++ b/src/algebra/IntegrationTools.hpp
@@ -7,21 +7,13 @@
#include <utils/PugsMacros.hpp>
#include <utils/Types.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <analysis/QuadratureType.hpp>
#include <analysis/TensorialGaussLegendreQuadrature.hpp>
#include <analysis/TensorialGaussLobattoQuadrature.hpp>
#include <cmath>
-enum class QuadratureType : int8_t
-{
- QT__begin = 0,
- //
- gausslobatto = QT__begin,
- gausslegendre = 1,
- //
- QT__end
-};
-
template <size_t Order, typename T = double>
class IntegrationTools
{
@@ -165,14 +157,16 @@ class IntegrationTools
IntegrationTools(QuadratureType quadrature)
{
switch (quadrature) {
- case QuadratureType::gausslobatto: {
- TensorialGaussLobattoQuadrature<1> tensorial_gauss_lobatto(Order);
+ case QuadratureType::GaussLobatto: {
+ auto tensorial_gauss_lobatto = QuadratureManager::instance().getLineGaussLobattoFormula(Order);
+
m_quadrature_positions = tensorial_gauss_lobatto.pointList();
m_weights = tensorial_gauss_lobatto.weightList();
break;
}
- case QuadratureType::gausslegendre: {
- TensorialGaussLegendreQuadrature<1> tensorial_gauss_legendre(Order);
+ case QuadratureType::GaussLegendre: {
+ auto tensorial_gauss_legendre = QuadratureManager::instance().getLineGaussLegendreFormula(Order);
+
m_quadrature_positions = tensorial_gauss_legendre.pointList();
m_weights = tensorial_gauss_legendre.weightList();
break;
@@ -180,6 +174,8 @@ class IntegrationTools
case QuadratureType::QT__end: {
throw NormalError("QuadratureType is not defined!");
}
+ default: {
+ }
}
}
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index c818e70b3ee6e10c68b4a4c0d7000225e594ac0e..3cbded76d7fadead5565ba862289aec3dcdbaad8 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(
PugsAnalysis
CubeGaussQuadrature.cpp
+ QuadratureManager.cpp
PrismGaussQuadrature.cpp
PyramidGaussQuadrature.cpp
SquareGaussQuadrature.cpp
diff --git a/src/analysis/CubeGaussQuadrature.cpp b/src/analysis/CubeGaussQuadrature.cpp
index cabac8663385882f80d38c5e21ae9dd6c5881b1c..1d35bd875ee46c0faf5d492328c34d9ec7e6466c 100644
--- a/src/analysis/CubeGaussQuadrature.cpp
+++ b/src/analysis/CubeGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-CubeGaussQuadrature::_buildPointAndWeightLists()
+CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R3 = TinyVector<3>;
@@ -237,7 +237,7 @@ CubeGaussQuadrature::_buildPointAndWeightLists()
}
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -494,7 +494,8 @@ CubeGaussQuadrature::_buildPointAndWeightLists()
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 21 on cubes");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(CubeGaussQuadrature::max_degree) + " on cubes");
}
}
}
diff --git a/src/analysis/CubeGaussQuadrature.hpp b/src/analysis/CubeGaussQuadrature.hpp
index bda80470dc9fb635fb8378be64f04a92d9b8d7c6..bd4ed179dd3cf3745e0dc4512eaa543f1ba6e1b9 100644
--- a/src/analysis/CubeGaussQuadrature.hpp
+++ b/src/analysis/CubeGaussQuadrature.hpp
@@ -12,18 +12,24 @@
* "Addendum to the paper 'High-order symmetric cubature rules for
* tetrahedra and pyramids'" Jan Jaśkowiec & N. Sukumar (2020).
*/
-class CubeGaussQuadrature final : public QuadratureForumla<3>
+class CubeGaussQuadrature final : public QuadratureFormula<3>
{
+ public:
+ constexpr static size_t max_degree = 21;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
CubeGaussQuadrature(CubeGaussQuadrature&&) = default;
CubeGaussQuadrature(const CubeGaussQuadrature&) = default;
- explicit CubeGaussQuadrature(const size_t degree) : QuadratureForumla<3>(degree)
+ private:
+ friend class QuadratureManager;
+
+ explicit CubeGaussQuadrature(const size_t degree) : QuadratureFormula<3>{QuadratureType::Gauss}
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
CubeGaussQuadrature() = delete;
diff --git a/src/analysis/PrismGaussQuadrature.cpp b/src/analysis/PrismGaussQuadrature.cpp
index 2077f035368fb7f231e2e69c237f7f758ea105ca..c2cde582ba1c4e2b6b5077e976a43dea5deaca77 100644
--- a/src/analysis/PrismGaussQuadrature.cpp
+++ b/src/analysis/PrismGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-PrismGaussQuadrature::_buildPointAndWeightLists()
+PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R2 = TinyVector<2>;
using R3 = TinyVector<3>;
@@ -142,7 +142,7 @@ PrismGaussQuadrature::_buildPointAndWeightLists()
}
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -797,7 +797,8 @@ PrismGaussQuadrature::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 20 on prisms");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(PrismGaussQuadrature::max_degree) + "on prisms");
}
}
}
diff --git a/src/analysis/PrismGaussQuadrature.hpp b/src/analysis/PrismGaussQuadrature.hpp
index 96ae6fbee84fc22ac404153e44a8e1d49a846924..2750d76c29894fd960b89afbfd29ec908c54d03b 100644
--- a/src/analysis/PrismGaussQuadrature.hpp
+++ b/src/analysis/PrismGaussQuadrature.hpp
@@ -11,18 +11,24 @@
* 'High-order symmetric cubature rules for tetrahedra and pyramids'
* Jan Jaśkowiec & N. Sukumar (2020).
*/
-class PrismGaussQuadrature final : public QuadratureForumla<3>
+class PrismGaussQuadrature final : public QuadratureFormula<3>
{
+ public:
+ constexpr static size_t max_degree = 20;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
PrismGaussQuadrature(PrismGaussQuadrature&&) = default;
PrismGaussQuadrature(const PrismGaussQuadrature&) = default;
- explicit PrismGaussQuadrature(const size_t degree) : QuadratureForumla<3>(degree)
+ private:
+ friend class QuadratureManager;
+
+ explicit PrismGaussQuadrature(const size_t degree) : QuadratureFormula<3>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
PrismGaussQuadrature() = delete;
diff --git a/src/analysis/PyramidGaussQuadrature.cpp b/src/analysis/PyramidGaussQuadrature.cpp
index d0ed132aa8820da9839984454313b58b255f6092..869408a914c58c2402e2b9b8c7954549466a97f7 100644
--- a/src/analysis/PyramidGaussQuadrature.cpp
+++ b/src/analysis/PyramidGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-PyramidGaussQuadrature::_buildPointAndWeightLists()
+PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R3 = TinyVector<3>;
@@ -97,7 +97,7 @@ PyramidGaussQuadrature::_buildPointAndWeightLists()
}
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -969,7 +969,8 @@ PyramidGaussQuadrature::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 20 on pyramids");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(PyramidGaussQuadrature::max_degree) + " on pyramids");
}
}
}
diff --git a/src/analysis/PyramidGaussQuadrature.hpp b/src/analysis/PyramidGaussQuadrature.hpp
index c5bf15af625463e1a1daccf5805758ff8ba33264..b68c518040e39611a0f1670758f40563271abb11 100644
--- a/src/analysis/PyramidGaussQuadrature.hpp
+++ b/src/analysis/PyramidGaussQuadrature.hpp
@@ -11,18 +11,24 @@
* 'High-order symmetric cubature rules for tetrahedra and pyramids'
* Jan Jaśkowiec & N. Sukumar (2020).
*/
-class PyramidGaussQuadrature final : public QuadratureForumla<3>
+class PyramidGaussQuadrature final : public QuadratureFormula<3>
{
+ public:
+ constexpr static size_t max_degree = 20;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
PyramidGaussQuadrature(PyramidGaussQuadrature&&) = default;
PyramidGaussQuadrature(const PyramidGaussQuadrature&) = default;
- explicit PyramidGaussQuadrature(const size_t degree) : QuadratureForumla<3>(degree)
+ private:
+ friend class QuadratureManager;
+
+ explicit PyramidGaussQuadrature(const size_t degree) : QuadratureFormula<3>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
PyramidGaussQuadrature() = delete;
diff --git a/src/analysis/QuadratureFormula.hpp b/src/analysis/QuadratureFormula.hpp
index ece948fe790c29b7511428cfdc9b045b5d064c61..c472e3556a298d14df73d1cb45a8d32c9fe6dc28 100644
--- a/src/analysis/QuadratureFormula.hpp
+++ b/src/analysis/QuadratureFormula.hpp
@@ -2,20 +2,20 @@
#define QUADRATURE_FORMULA_HPP
#include <algebra/TinyVector.hpp>
+#include <analysis/QuadratureType.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <utils/SmallArray.hpp>
template <size_t Dimension>
-class QuadratureForumla
+class QuadratureFormula
{
+ public:
protected:
- size_t m_degree;
+ QuadratureType m_type;
SmallArray<const TinyVector<Dimension>> m_point_list;
SmallArray<const double> m_weight_list;
- virtual void _buildPointAndWeightLists() = 0;
-
public:
PUGS_INLINE size_t
numberOfPoints() const
@@ -52,15 +52,12 @@ class QuadratureForumla
return m_weight_list[i];
}
- QuadratureForumla(QuadratureForumla&&) = default;
- QuadratureForumla(const QuadratureForumla&) = default;
-
protected:
- explicit QuadratureForumla(const size_t degree) : m_degree{degree} {}
+ explicit QuadratureFormula(const QuadratureType type) : m_type{type} {}
public:
- QuadratureForumla() = delete;
- virtual ~QuadratureForumla() = default;
+ QuadratureFormula() = default;
+ virtual ~QuadratureFormula() = default;
};
#endif // QUADRATURE_FORMULA_HPP
diff --git a/src/analysis/QuadratureManager.cpp b/src/analysis/QuadratureManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a8063a2f6c758d16ec6b3ac1490e72e3b46fea8d
--- /dev/null
+++ b/src/analysis/QuadratureManager.cpp
@@ -0,0 +1,194 @@
+#include <analysis/QuadratureManager.hpp>
+
+#include <analysis/CubeGaussQuadrature.hpp>
+#include <analysis/PrismGaussQuadrature.hpp>
+#include <analysis/PyramidGaussQuadrature.hpp>
+#include <analysis/SquareGaussQuadrature.hpp>
+#include <analysis/TensorialGaussLegendreQuadrature.hpp>
+#include <analysis/TensorialGaussLobattoQuadrature.hpp>
+#include <analysis/TetrahedronGaussQuadrature.hpp>
+#include <analysis/TriangleGaussQuadrature.hpp>
+
+QuadratureManager* QuadratureManager::s_instance{nullptr};
+
+void
+QuadratureManager::create()
+{
+ Assert(s_instance == nullptr, "QuadratureManager is already created");
+ s_instance = new QuadratureManager;
+}
+
+void
+QuadratureManager::destroy()
+{
+ Assert(s_instance != nullptr, "QuadratureManager was not created!");
+
+ delete s_instance;
+ s_instance = nullptr;
+}
+
+Array<const QuadratureFormula<1>>
+QuadratureManager::_buildLineGaussLegendreFormulaList()
+{
+ Array<QuadratureFormula<1>> line_gauss_legendre_formula_list(TensorialGaussLegendreQuadrature<1>::max_degree / 2 + 1);
+ for (size_t i = 0; i < TensorialGaussLegendreQuadrature<1>::max_degree / 2 + 1; ++i) {
+ line_gauss_legendre_formula_list[i] = TensorialGaussLegendreQuadrature<1>(2 * i + 1);
+ }
+
+ return line_gauss_legendre_formula_list;
+}
+
+Array<const QuadratureFormula<1>>
+QuadratureManager::_buildLineGaussLobattoFormulaList()
+{
+ Array<QuadratureFormula<1>> line_gauss_lobatto_formula_list(TensorialGaussLobattoQuadrature<1>::max_degree / 2 + 1);
+ for (size_t i = 0; i < TensorialGaussLobattoQuadrature<1>::max_degree / 2 + 1; ++i) {
+ line_gauss_lobatto_formula_list[i] = TensorialGaussLobattoQuadrature<1>(2 * i + 1);
+ }
+
+ return line_gauss_lobatto_formula_list;
+}
+
+Array<const QuadratureFormula<2>>
+QuadratureManager::_buildSquareGaussFormulaList()
+{
+ Array<QuadratureFormula<2>> square_gauss_formula_list(SquareGaussQuadrature::max_degree / 2 + 1);
+ for (size_t i = 0; i < SquareGaussQuadrature::max_degree / 2 + 1; ++i) {
+ square_gauss_formula_list[i] = SquareGaussQuadrature(2 * i + 1);
+ }
+
+ return square_gauss_formula_list;
+}
+
+Array<const QuadratureFormula<2>>
+QuadratureManager::_buildSquareGaussLegendreFormulaList()
+{
+ Array<QuadratureFormula<2>> square_gauss_legendre_formula_list(TensorialGaussLegendreQuadrature<2>::max_degree / 2 +
+ 1);
+ for (size_t i = 0; i < TensorialGaussLegendreQuadrature<2>::max_degree / 2 + 1; ++i) {
+ square_gauss_legendre_formula_list[i] = TensorialGaussLegendreQuadrature<2>(2 * i + 1);
+ }
+
+ return square_gauss_legendre_formula_list;
+}
+
+Array<const QuadratureFormula<2>>
+QuadratureManager::_buildSquareGaussLobattoFormulaList()
+{
+ Array<QuadratureFormula<2>> square_gauss_lobatto_formula_list(TensorialGaussLobattoQuadrature<2>::max_degree / 2 + 1);
+ for (size_t i = 0; i < TensorialGaussLobattoQuadrature<2>::max_degree / 2 + 1; ++i) {
+ square_gauss_lobatto_formula_list[i] = TensorialGaussLobattoQuadrature<2>(2 * i + 1);
+ }
+
+ return square_gauss_lobatto_formula_list;
+}
+
+Array<const QuadratureFormula<2>>
+QuadratureManager::_buildTriangleGaussFormulaList()
+{
+ Array<QuadratureFormula<2>> triangle_gauss_formula_list(TriangleGaussQuadrature::max_degree);
+ for (size_t i = 0; i < TriangleGaussQuadrature::max_degree; ++i) {
+ triangle_gauss_formula_list[i] = TriangleGaussQuadrature(i + 1);
+ }
+
+ return triangle_gauss_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildCubeGaussFormulaList()
+{
+ Array<QuadratureFormula<3>> cube_gauss_formula_list(CubeGaussQuadrature::max_degree / 2 + 1);
+ for (size_t i = 0; i < CubeGaussQuadrature::max_degree / 2 + 1; ++i) {
+ cube_gauss_formula_list[i] = CubeGaussQuadrature(2 * i + 1);
+ }
+
+ return cube_gauss_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildCubeGaussLegendreFormulaList()
+{
+ Array<QuadratureFormula<3>> cube_gauss_legendre_formula_list(TensorialGaussLegendreQuadrature<3>::max_degree / 2 + 1);
+ for (size_t i = 0; i < TensorialGaussLegendreQuadrature<3>::max_degree / 2 + 1; ++i) {
+ cube_gauss_legendre_formula_list[i] = TensorialGaussLegendreQuadrature<3>(2 * i + 1);
+ }
+
+ return cube_gauss_legendre_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildCubeGaussLobattoFormulaList()
+{
+ Array<QuadratureFormula<3>> cube_gauss_lobatto_formula_list(TensorialGaussLobattoQuadrature<3>::max_degree / 2 + 1);
+ for (size_t i = 0; i < TensorialGaussLobattoQuadrature<3>::max_degree / 2 + 1; ++i) {
+ cube_gauss_lobatto_formula_list[i] = TensorialGaussLobattoQuadrature<3>(2 * i + 1);
+ }
+
+ return cube_gauss_lobatto_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildPrismGaussFormulaList()
+{
+ Array<QuadratureFormula<3>> prism_gauss_formula_list(PrismGaussQuadrature::max_degree);
+ for (size_t i = 0; i < PrismGaussQuadrature::max_degree; ++i) {
+ prism_gauss_formula_list[i] = PrismGaussQuadrature(i + 1);
+ }
+
+ return prism_gauss_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildPyramidGaussFormulaList()
+{
+ Array<QuadratureFormula<3>> pyramid_gauss_formula_list(PyramidGaussQuadrature::max_degree);
+ for (size_t i = 0; i < PyramidGaussQuadrature::max_degree; ++i) {
+ pyramid_gauss_formula_list[i] = PyramidGaussQuadrature(i + 1);
+ }
+
+ return pyramid_gauss_formula_list;
+}
+
+Array<const QuadratureFormula<3>>
+QuadratureManager::_buildTetrahedronGaussFormulaList()
+{
+ Array<QuadratureFormula<3>> tetrahedron_gauss_formula_list(TetrahedronGaussQuadrature::max_degree);
+ for (size_t i = 0; i < TetrahedronGaussQuadrature::max_degree; ++i) {
+ tetrahedron_gauss_formula_list[i] = TetrahedronGaussQuadrature(i + 1);
+ }
+
+ return tetrahedron_gauss_formula_list;
+}
+
+QuadratureManager::QuadratureManager()
+ : m_line_gauss_legendre_formula_list{this->_buildLineGaussLegendreFormulaList()},
+ m_line_gauss_lobatto_formula_list{this->_buildLineGaussLobattoFormulaList()},
+ //
+ m_square_gauss_formula_list{this->_buildSquareGaussFormulaList()},
+ m_square_gauss_legendre_formula_list{this->_buildSquareGaussLegendreFormulaList()},
+ m_square_gauss_lobatto_formula_list{this->_buildSquareGaussLobattoFormulaList()},
+ m_triangle_gauss_formula_list{this->_buildTriangleGaussFormulaList()},
+ //
+ m_cube_gauss_formula_list{this->_buildCubeGaussFormulaList()},
+ m_cube_gauss_legendre_formula_list{this->_buildCubeGaussLegendreFormulaList()},
+ m_cube_gauss_lobatto_formula_list{this->_buildCubeGaussLobattoFormulaList()},
+ m_prism_gauss_formula_list{this->_buildPrismGaussFormulaList()},
+ m_pyramid_gauss_formula_list{this->_buildPyramidGaussFormulaList()},
+ m_tetrahedron_gauss_formula_list{this->_buildTetrahedronGaussFormulaList()}
+{
+ Assert(this->maxGaussLegendreDegree() == TensorialGaussLegendreQuadrature<1>::max_degree);
+ Assert(this->maxGaussLegendreDegree() == TensorialGaussLegendreQuadrature<2>::max_degree);
+ Assert(this->maxGaussLegendreDegree() == TensorialGaussLegendreQuadrature<3>::max_degree);
+
+ Assert(this->maxGaussLobattoDegree() == TensorialGaussLobattoQuadrature<1>::max_degree);
+ Assert(this->maxGaussLobattoDegree() == TensorialGaussLobattoQuadrature<2>::max_degree);
+ Assert(this->maxGaussLobattoDegree() == TensorialGaussLobattoQuadrature<3>::max_degree);
+
+ Assert(this->maxSquareGaussDegree() == SquareGaussQuadrature::max_degree);
+ Assert(this->maxTriangleGaussDegree() == TriangleGaussQuadrature::max_degree);
+
+ Assert(this->maxCubeGaussDegree() == CubeGaussQuadrature::max_degree);
+ Assert(this->maxPrismGaussDegree() == PrismGaussQuadrature::max_degree);
+ Assert(this->maxPyramidGaussDegree() == PyramidGaussQuadrature::max_degree);
+ Assert(this->maxTetrahedronGaussDegree() == TetrahedronGaussQuadrature::max_degree);
+}
diff --git a/src/analysis/QuadratureManager.hpp b/src/analysis/QuadratureManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..875a2a3fdbc71a7b051d0b70cccaf1ade8822fd0
--- /dev/null
+++ b/src/analysis/QuadratureManager.hpp
@@ -0,0 +1,194 @@
+#ifndef QUADRATURE_MANAGER_HPP
+#define QUADRATURE_MANAGER_HPP
+
+#include <analysis/QuadratureFormula.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+class QuadratureManager
+{
+ private:
+ static QuadratureManager* s_instance;
+
+ Array<const QuadratureFormula<1>> m_line_gauss_legendre_formula_list;
+ Array<const QuadratureFormula<1>> m_line_gauss_lobatto_formula_list;
+
+ Array<const QuadratureFormula<2>> m_square_gauss_formula_list;
+ Array<const QuadratureFormula<2>> m_square_gauss_legendre_formula_list;
+ Array<const QuadratureFormula<2>> m_square_gauss_lobatto_formula_list;
+ Array<const QuadratureFormula<2>> m_triangle_gauss_formula_list;
+
+ Array<const QuadratureFormula<3>> m_cube_gauss_formula_list;
+ Array<const QuadratureFormula<3>> m_cube_gauss_legendre_formula_list;
+ Array<const QuadratureFormula<3>> m_cube_gauss_lobatto_formula_list;
+ Array<const QuadratureFormula<3>> m_prism_gauss_formula_list;
+ Array<const QuadratureFormula<3>> m_pyramid_gauss_formula_list;
+ Array<const QuadratureFormula<3>> m_tetrahedron_gauss_formula_list;
+
+ Array<const QuadratureFormula<1>> _buildLineGaussLobattoFormulaList();
+ Array<const QuadratureFormula<1>> _buildLineGaussLegendreFormulaList();
+
+ Array<const QuadratureFormula<2>> _buildSquareGaussFormulaList();
+ Array<const QuadratureFormula<2>> _buildSquareGaussLegendreFormulaList();
+ Array<const QuadratureFormula<2>> _buildSquareGaussLobattoFormulaList();
+ Array<const QuadratureFormula<2>> _buildTriangleGaussFormulaList();
+
+ Array<const QuadratureFormula<3>> _buildCubeGaussFormulaList();
+ Array<const QuadratureFormula<3>> _buildCubeGaussLegendreFormulaList();
+ Array<const QuadratureFormula<3>> _buildCubeGaussLobattoFormulaList();
+ Array<const QuadratureFormula<3>> _buildPrismGaussFormulaList();
+ Array<const QuadratureFormula<3>> _buildPyramidGaussFormulaList();
+ Array<const QuadratureFormula<3>> _buildTetrahedronGaussFormulaList();
+
+ public:
+ size_t
+ maxGaussLegendreDegree() const
+ {
+ Assert(m_square_gauss_legendre_formula_list.size() == m_line_gauss_legendre_formula_list.size());
+ Assert(m_cube_gauss_legendre_formula_list.size() == m_line_gauss_legendre_formula_list.size());
+
+ return m_line_gauss_legendre_formula_list.size() * 2 - 1;
+ }
+
+ size_t
+ maxGaussLobattoDegree() const
+ {
+ Assert(m_square_gauss_lobatto_formula_list.size() == m_line_gauss_lobatto_formula_list.size());
+ Assert(m_cube_gauss_lobatto_formula_list.size() == m_line_gauss_lobatto_formula_list.size());
+
+ return m_line_gauss_lobatto_formula_list.size() * 2 - 1;
+ }
+
+ const QuadratureFormula<1>&
+ getLineGaussLegendreFormula(const size_t degree) const
+ {
+ return m_line_gauss_legendre_formula_list[degree / 2];
+ }
+
+ const QuadratureFormula<1>&
+ getLineGaussLobattoFormula(const size_t degree) const
+ {
+ return m_line_gauss_lobatto_formula_list[degree / 2];
+ }
+
+ size_t
+ maxSquareGaussDegree() const
+ {
+ return m_square_gauss_formula_list.size() * 2 - 1;
+ }
+
+ const QuadratureFormula<2>&
+ getSquareGaussFormula(const size_t degree) const
+ {
+ return m_square_gauss_formula_list[degree / 2];
+ }
+
+ const QuadratureFormula<2>&
+ getSquareGaussLegendreFormula(const size_t degree) const
+ {
+ return m_square_gauss_legendre_formula_list[degree / 2];
+ }
+
+ const QuadratureFormula<2>&
+ getSquareGaussLobattoFormula(const size_t degree) const
+ {
+ return m_square_gauss_lobatto_formula_list[degree / 2];
+ }
+
+ size_t
+ maxTriangleGaussDegree() const
+ {
+ return m_triangle_gauss_formula_list.size();
+ }
+
+ const QuadratureFormula<2>&
+ getTriangleGaussFormula(const size_t degree) const
+ {
+ Assert(degree > 0);
+ return m_triangle_gauss_formula_list[degree - 1];
+ }
+
+ size_t
+ maxCubeGaussDegree() const
+ {
+ return m_cube_gauss_formula_list.size() * 2 - 1;
+ }
+
+ const QuadratureFormula<3>&
+ getCubeGaussFormula(const size_t degree) const
+ {
+ return m_cube_gauss_formula_list[degree / 2];
+ }
+
+ const QuadratureFormula<3>&
+ getCubeGaussLegendreFormula(const size_t degree) const
+ {
+ return m_cube_gauss_legendre_formula_list[degree / 2];
+ }
+
+ const QuadratureFormula<3>&
+ getCubeGaussLobattoFormula(const size_t degree) const
+ {
+ return m_cube_gauss_lobatto_formula_list[degree / 2];
+ }
+
+ size_t
+ maxPrismGaussDegree() const
+ {
+ return m_prism_gauss_formula_list.size();
+ }
+
+ const QuadratureFormula<3>&
+ getPrismGaussFormula(const size_t degree) const
+ {
+ Assert(degree > 0);
+ return m_prism_gauss_formula_list[degree - 1];
+ }
+
+ size_t
+ maxPyramidGaussDegree() const
+ {
+ return m_prism_gauss_formula_list.size();
+ }
+
+ const QuadratureFormula<3>&
+ getPyramidGaussFormula(const size_t degree) const
+ {
+ Assert(degree > 0);
+ return m_pyramid_gauss_formula_list[degree - 1];
+ }
+
+ size_t
+ maxTetrahedronGaussDegree() const
+ {
+ return m_tetrahedron_gauss_formula_list.size();
+ }
+
+ const QuadratureFormula<3>&
+ getTetrahedronGaussFormula(const size_t degree) const
+ {
+ Assert(degree > 0);
+ return m_tetrahedron_gauss_formula_list[degree - 1];
+ }
+
+ static void create();
+ static void destroy();
+
+ PUGS_INLINE
+ static const QuadratureManager&
+ instance()
+ {
+ Assert(s_instance != nullptr, "QuadratureManager was not created!");
+ return *s_instance;
+ }
+
+ private:
+ QuadratureManager(const QuadratureManager&) = delete;
+ QuadratureManager(QuadratureManager&&) = delete;
+
+ QuadratureManager();
+ ~QuadratureManager() = default;
+};
+
+#endif // QUADRATURE_MANAGER_HPP
diff --git a/src/analysis/QuadratureType.hpp b/src/analysis/QuadratureType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dca48806d3f00dc06e000bbeca34fd02453d50e1
--- /dev/null
+++ b/src/analysis/QuadratureType.hpp
@@ -0,0 +1,15 @@
+#ifndef QUADRATURE_TYPE_HPP
+#define QUADRATURE_TYPE_HPP
+
+enum class QuadratureType
+{
+ QT__begin = 0,
+ //
+ Gauss = QT__begin,
+ GaussLegendre = 1,
+ GaussLobatto = 2,
+ //
+ QT__end
+};
+
+#endif // QUADRATURE_TYPE_HPP
diff --git a/src/analysis/SquareGaussQuadrature.cpp b/src/analysis/SquareGaussQuadrature.cpp
index 8412638280a8fadc3e71714934c9fad1bbfb58ec..f0df6af9f2cbbc23318e861403dbaadaa844d749 100644
--- a/src/analysis/SquareGaussQuadrature.cpp
+++ b/src/analysis/SquareGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-SquareGaussQuadrature::_buildPointAndWeightLists()
+SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R2 = TinyVector<2>;
@@ -94,7 +94,7 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
Assert(k == point_list.size(), "invalid number of quadrature points");
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -320,7 +320,8 @@ SquareGaussQuadrature::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 21 on squares");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(SquareGaussQuadrature::max_degree) + " on squares");
}
}
}
diff --git a/src/analysis/SquareGaussQuadrature.hpp b/src/analysis/SquareGaussQuadrature.hpp
index 8dd1a00225658c0088c677702458e1353d903ad8..71ea25a7a7f25c65f862fbddc92f31f4d5057fa4 100644
--- a/src/analysis/SquareGaussQuadrature.hpp
+++ b/src/analysis/SquareGaussQuadrature.hpp
@@ -12,18 +12,24 @@
* ‘On the identification of symmetric quadrature rules for finite
* element methods‘ by F.D. Witherden & P.E. Vincent (2015).
*/
-class SquareGaussQuadrature final : public QuadratureForumla<2>
+class SquareGaussQuadrature final : public QuadratureFormula<2>
{
+ public:
+ constexpr static size_t max_degree = 21;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
SquareGaussQuadrature(SquareGaussQuadrature&&) = default;
SquareGaussQuadrature(const SquareGaussQuadrature&) = default;
- explicit SquareGaussQuadrature(const size_t degree) : QuadratureForumla<2>(degree)
+ private:
+ friend class QuadratureManager;
+
+ explicit SquareGaussQuadrature(const size_t degree) : QuadratureFormula<2>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
SquareGaussQuadrature() = delete;
diff --git a/src/analysis/TensorialGaussLegendreQuadrature.cpp b/src/analysis/TensorialGaussLegendreQuadrature.cpp
index ef3796f85bf45e9f369e4344abaa0130e69c28cd..9a4b524c0f70cbc3a68088bd5b14ad17d6c521ea 100644
--- a/src/analysis/TensorialGaussLegendreQuadrature.cpp
+++ b/src/analysis/TensorialGaussLegendreQuadrature.cpp
@@ -3,14 +3,14 @@
template <>
void
-TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists()
+TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points = m_degree / 2 + 1;
+ const size_t nb_points = degree / 2 + 1;
SmallArray<TinyVector<1>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
point_list[0] = 0;
@@ -228,7 +228,8 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss-Legendre quadrature formulae handle degrees up to 23");
+ throw NormalError("Gauss-Legendre quadrature formulae handle degrees up to " +
+ std::to_string(TensorialGaussLegendreQuadrature<1>::max_degree));
}
}
@@ -238,15 +239,15 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists()
template <>
void
-TensorialGaussLegendreQuadrature<2>::_buildPointAndWeightLists()
+TensorialGaussLegendreQuadrature<2>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points_1d = this->m_degree / 2 + 1;
+ const size_t nb_points_1d = degree / 2 + 1;
const size_t nb_points = nb_points_1d * nb_points_1d;
SmallArray<TinyVector<2>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- TensorialGaussLegendreQuadrature<1> gauss_legendre_1d(m_degree);
+ TensorialGaussLegendreQuadrature<1> gauss_legendre_1d(degree);
const auto& point_list_1d = gauss_legendre_1d.pointList();
const auto& weight_list_1d = gauss_legendre_1d.weightList();
@@ -264,15 +265,15 @@ TensorialGaussLegendreQuadrature<2>::_buildPointAndWeightLists()
template <>
void
-TensorialGaussLegendreQuadrature<3>::_buildPointAndWeightLists()
+TensorialGaussLegendreQuadrature<3>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points_1d = this->m_degree / 2 + 1;
+ const size_t nb_points_1d = degree / 2 + 1;
const size_t nb_points = nb_points_1d * nb_points_1d * nb_points_1d;
SmallArray<TinyVector<3>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- TensorialGaussLegendreQuadrature<1> gauss_legendre_1d(m_degree);
+ TensorialGaussLegendreQuadrature<1> gauss_legendre_1d(degree);
const auto& point_list_1d = gauss_legendre_1d.pointList();
const auto& weight_list_1d = gauss_legendre_1d.weightList();
diff --git a/src/analysis/TensorialGaussLegendreQuadrature.hpp b/src/analysis/TensorialGaussLegendreQuadrature.hpp
index 61ec916d7dba869ce47636b847b3ddfc0d48cfc7..51aadae6f64b13f2f8c3b68a59665174ad74150e 100644
--- a/src/analysis/TensorialGaussLegendreQuadrature.hpp
+++ b/src/analysis/TensorialGaussLegendreQuadrature.hpp
@@ -10,18 +10,27 @@
* \note formulae are extracted from High-order Finite Element Method [2004 - Chapman & Hall]
*/
template <size_t Dimension>
-class TensorialGaussLegendreQuadrature final : public QuadratureForumla<Dimension>
+class TensorialGaussLegendreQuadrature final : public QuadratureFormula<Dimension>
{
+ public:
+ constexpr static size_t max_degree = 23;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
TensorialGaussLegendreQuadrature(TensorialGaussLegendreQuadrature&&) = default;
TensorialGaussLegendreQuadrature(const TensorialGaussLegendreQuadrature&) = default;
- explicit TensorialGaussLegendreQuadrature(const size_t degree) : QuadratureForumla<Dimension>(degree)
+ private:
+ friend class QuadratureManager;
+
+ template <size_t D>
+ friend class TensorialGaussLegendreQuadrature;
+
+ explicit TensorialGaussLegendreQuadrature(const size_t degree) : QuadratureFormula<Dimension>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
TensorialGaussLegendreQuadrature() = delete;
diff --git a/src/analysis/TensorialGaussLobattoQuadrature.cpp b/src/analysis/TensorialGaussLobattoQuadrature.cpp
index fa1b3b5f352d757ad37ff15f9c801bad6458a74f..3ab97460b55dad3a2a48a24db56f5c575542361e 100644
--- a/src/analysis/TensorialGaussLobattoQuadrature.cpp
+++ b/src/analysis/TensorialGaussLobattoQuadrature.cpp
@@ -3,14 +3,14 @@
template <>
void
-TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists()
+TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points = m_degree / 2 + 2;
+ const size_t nb_points = degree / 2 + 2;
SmallArray<TinyVector<1>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
point_list[0] = -1;
@@ -117,7 +117,8 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss-Lobatto quadrature formulae handle degrees up to 13");
+ throw NormalError("Gauss-Lobatto quadrature formulae handle degrees up to " +
+ std::to_string(TensorialGaussLobattoQuadrature<1>::max_degree));
}
}
@@ -127,15 +128,15 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists()
template <>
void
-TensorialGaussLobattoQuadrature<2>::_buildPointAndWeightLists()
+TensorialGaussLobattoQuadrature<2>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points_1d = this->m_degree / 2 + 2;
+ const size_t nb_points_1d = degree / 2 + 2;
const size_t nb_points = nb_points_1d * nb_points_1d;
SmallArray<TinyVector<2>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- TensorialGaussLobattoQuadrature<1> gauss_lobatto_1d(m_degree);
+ TensorialGaussLobattoQuadrature<1> gauss_lobatto_1d(degree);
const auto& point_list_1d = gauss_lobatto_1d.pointList();
const auto& weight_list_1d = gauss_lobatto_1d.weightList();
@@ -153,15 +154,15 @@ TensorialGaussLobattoQuadrature<2>::_buildPointAndWeightLists()
template <>
void
-TensorialGaussLobattoQuadrature<3>::_buildPointAndWeightLists()
+TensorialGaussLobattoQuadrature<3>::_buildPointAndWeightLists(const size_t degree)
{
- const size_t nb_points_1d = this->m_degree / 2 + 2;
+ const size_t nb_points_1d = degree / 2 + 2;
const size_t nb_points = nb_points_1d * nb_points_1d * nb_points_1d;
SmallArray<TinyVector<3>> point_list(nb_points);
SmallArray<double> weight_list(nb_points);
- TensorialGaussLobattoQuadrature<1> gauss_lobatto_1d(m_degree);
+ TensorialGaussLobattoQuadrature<1> gauss_lobatto_1d(degree);
const auto& point_list_1d = gauss_lobatto_1d.pointList();
const auto& weight_list_1d = gauss_lobatto_1d.weightList();
diff --git a/src/analysis/TensorialGaussLobattoQuadrature.hpp b/src/analysis/TensorialGaussLobattoQuadrature.hpp
index fe63573f8b36532b6a03d6978cae0ad6083ae107..1a578bbbc1af1d7d4efdae4183cea1e3c6fb76ff 100644
--- a/src/analysis/TensorialGaussLobattoQuadrature.hpp
+++ b/src/analysis/TensorialGaussLobattoQuadrature.hpp
@@ -10,18 +10,27 @@
* \note formulae are extracted from High-order Finite Element Method [2004 - Chapman & Hall]
*/
template <size_t Dimension>
-class TensorialGaussLobattoQuadrature final : public QuadratureForumla<Dimension>
+class TensorialGaussLobattoQuadrature final : public QuadratureFormula<Dimension>
{
+ public:
+ constexpr static size_t max_degree = 13;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
TensorialGaussLobattoQuadrature(TensorialGaussLobattoQuadrature&&) = default;
TensorialGaussLobattoQuadrature(const TensorialGaussLobattoQuadrature&) = default;
- explicit TensorialGaussLobattoQuadrature(const size_t degree) : QuadratureForumla<Dimension>(degree)
+ private:
+ friend class QuadratureManager;
+
+ template <size_t D>
+ friend class TensorialGaussLobattoQuadrature;
+
+ explicit TensorialGaussLobattoQuadrature(const size_t degree) : QuadratureFormula<Dimension>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
TensorialGaussLobattoQuadrature() = delete;
diff --git a/src/analysis/TetrahedronGaussQuadrature.cpp b/src/analysis/TetrahedronGaussQuadrature.cpp
index 38c752e7f2b02bdfbbbe70211853613ec900e91f..644e2e29387821a274fa2ae8dd2861c60f2c4abc 100644
--- a/src/analysis/TetrahedronGaussQuadrature.cpp
+++ b/src/analysis/TetrahedronGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-TetrahedronGaussQuadrature::_buildPointAndWeightLists()
+TetrahedronGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R3 = TinyVector<3>;
@@ -145,7 +145,7 @@ TetrahedronGaussQuadrature::_buildPointAndWeightLists()
}
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -681,7 +681,8 @@ TetrahedronGaussQuadrature::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 20 on tetrahedra");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(TetrahedronGaussQuadrature::max_degree) + " on tetrahedra");
}
}
}
diff --git a/src/analysis/TetrahedronGaussQuadrature.hpp b/src/analysis/TetrahedronGaussQuadrature.hpp
index 0b1127cd986695946e4582a68c749d5bd866403c..352ff079a409952d99af6931bf75b8b84b06d622 100644
--- a/src/analysis/TetrahedronGaussQuadrature.hpp
+++ b/src/analysis/TetrahedronGaussQuadrature.hpp
@@ -11,18 +11,24 @@
* 'High-order symmetric cubature rules for tetrahedra and pyramids'
* Jan Jaśkowiec & N. Sukumar (2020).
*/
-class TetrahedronGaussQuadrature final : public QuadratureForumla<3>
+class TetrahedronGaussQuadrature final : public QuadratureFormula<3>
{
+ public:
+ constexpr static size_t max_degree = 20;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
TetrahedronGaussQuadrature(TetrahedronGaussQuadrature&&) = default;
TetrahedronGaussQuadrature(const TetrahedronGaussQuadrature&) = default;
- explicit TetrahedronGaussQuadrature(const size_t degree) : QuadratureForumla<3>(degree)
+ private:
+ friend class QuadratureManager;
+
+ explicit TetrahedronGaussQuadrature(const size_t degree) : QuadratureFormula<3>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
TetrahedronGaussQuadrature() = delete;
diff --git a/src/analysis/TriangleGaussQuadrature.cpp b/src/analysis/TriangleGaussQuadrature.cpp
index 6a8285ef5767b8f73e6d212979d6cbbe6a9b429d..df6a9a1ced3aa15d098aaa7696eb67af41c6090b 100644
--- a/src/analysis/TriangleGaussQuadrature.cpp
+++ b/src/analysis/TriangleGaussQuadrature.cpp
@@ -2,7 +2,7 @@
#include <utils/Exceptions.hpp>
void
-TriangleGaussQuadrature::_buildPointAndWeightLists()
+TriangleGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
{
using R2 = TinyVector<2>;
@@ -81,7 +81,7 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
Assert(k == point_list.size(), "invalid number of quadrature points");
};
- switch (m_degree) {
+ switch (degree) {
case 0:
case 1: {
constexpr size_t nb_points = 1;
@@ -483,7 +483,8 @@ TriangleGaussQuadrature::_buildPointAndWeightLists()
break;
}
default: {
- throw NormalError("Gauss quadrature formulae handle degrees up to 20 on triangles");
+ throw NormalError("Gauss quadrature formulae handle degrees up to " +
+ std::to_string(TriangleGaussQuadrature::max_degree) + " on triangles");
}
}
}
diff --git a/src/analysis/TriangleGaussQuadrature.hpp b/src/analysis/TriangleGaussQuadrature.hpp
index 58b0be03e2efb7b23e6247b2564e238ac1f01734..9e1edd3dedc01392f5cfae1ed8c0ab277b3ab2ba 100644
--- a/src/analysis/TriangleGaussQuadrature.hpp
+++ b/src/analysis/TriangleGaussQuadrature.hpp
@@ -11,18 +11,23 @@
* ‘On the identification of symmetric quadrature rules for finite
* element methods‘ by F.D. Witherden & P.E. Vincent (2015).
*/
-class TriangleGaussQuadrature final : public QuadratureForumla<2>
+class TriangleGaussQuadrature final : public QuadratureFormula<2>
{
+ public:
+ constexpr static size_t max_degree = 20;
+
private:
- void _buildPointAndWeightLists();
+ void _buildPointAndWeightLists(const size_t degree);
public:
TriangleGaussQuadrature(TriangleGaussQuadrature&&) = default;
TriangleGaussQuadrature(const TriangleGaussQuadrature&) = default;
- explicit TriangleGaussQuadrature(const size_t degree) : QuadratureForumla<2>(degree)
+ private:
+ friend class QuadratureManager;
+ explicit TriangleGaussQuadrature(const size_t degree) : QuadratureFormula<2>(QuadratureType::Gauss)
{
- this->_buildPointAndWeightLists();
+ this->_buildPointAndWeightLists(degree);
}
TriangleGaussQuadrature() = delete;
diff --git a/src/main.cpp b/src/main.cpp
index cf382d69b491e85968d977e9023ac094d81c3bc2..4b9dfffd985bc5f400ace5986355f55a591ff9f2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,3 +1,4 @@
+#include <analysis/QuadratureManager.hpp>
#include <language/PugsParser.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
@@ -13,6 +14,7 @@ main(int argc, char* argv[])
SynchronizerManager::create();
RandomEngine::create();
+ QuadratureManager::create();
MeshDataManager::create();
DiamondDualConnectivityManager::create();
DiamondDualMeshManager::create();
@@ -23,6 +25,7 @@ main(int argc, char* argv[])
DiamondDualConnectivityManager::destroy();
MeshDataManager::destroy();
RandomEngine::destroy();
+ QuadratureManager::destroy();
SynchronizerManager::destroy();
finalize();
diff --git a/tests/mpi_test_main.cpp b/tests/mpi_test_main.cpp
index 3f195f2ec81f455f46604d2e30cf5a3ab81bceaa..95b26e0a4964472335139a7afc0cdbf8205ff61f 100644
--- a/tests/mpi_test_main.cpp
+++ b/tests/mpi_test_main.cpp
@@ -2,6 +2,7 @@
#include <Kokkos_Core.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <language/utils/OperatorRepository.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
@@ -59,6 +60,7 @@ main(int argc, char* argv[])
SynchronizerManager::create();
RandomEngine::create();
+ QuadratureManager::create();
MeshDataManager::create();
DiamondDualConnectivityManager::create();
DiamondDualMeshManager::create();
@@ -91,6 +93,7 @@ main(int argc, char* argv[])
DiamondDualMeshManager::destroy();
DiamondDualConnectivityManager::destroy();
MeshDataManager::destroy();
+ QuadratureManager::destroy();
RandomEngine::destroy();
SynchronizerManager::destroy();
}
diff --git a/tests/test_CubeGaussQuadrature.cpp b/tests/test_CubeGaussQuadrature.cpp
index 299e0d6f91ebb59a9152b2a1897fb2a1043ac926..73d687a621c5975ec525d342067d5a14c3b18efb 100644
--- a/tests/test_CubeGaussQuadrature.cpp
+++ b/tests/test_CubeGaussQuadrature.cpp
@@ -5,35 +5,13 @@
#include <algebra/TinyMatrix.hpp>
#include <analysis/CubeGaussQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <utils/Exceptions.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("CubeGaussQuadrature", "[analysis]")
{
- auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
- auto point_list0 = quadrature_formula0.pointList();
- auto weight_list0 = quadrature_formula0.weightList();
-
- auto point_list1 = quadrature_formula1.pointList();
- auto weight_list1 = quadrature_formula1.weightList();
-
- bool is_same = true;
-
- for (size_t i = 0; i < point_list0.size(); ++i) {
- if (point_list0[i] != point_list1[i]) {
- return false;
- }
- }
- for (size_t i = 0; i < weight_list0.size(); ++i) {
- if (weight_list0[i] != weight_list1[i]) {
- return false;
- }
- }
-
- return is_same;
- };
-
auto integrate = [](auto f, auto quadrature_formula) {
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -224,7 +202,7 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 0 and 1")
{
- CubeGaussQuadrature l1(1);
+ const QuadratureFormula<3>& l1 = QuadratureManager::instance().getCubeGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -237,10 +215,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 2 and 3")
{
- CubeGaussQuadrature l2(2);
- CubeGaussQuadrature l3(2);
+ const QuadratureFormula<3>& l2 = QuadratureManager::instance().getCubeGaussFormula(2);
+ const QuadratureFormula<3>& l3 = QuadratureManager::instance().getCubeGaussFormula(3);
- REQUIRE(same_formulas(l2, l3));
+ REQUIRE(&l2 == &l3);
REQUIRE(l3.numberOfPoints() == 6);
@@ -255,10 +233,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 4 and 5")
{
- CubeGaussQuadrature l4(4);
- CubeGaussQuadrature l5(5);
+ const QuadratureFormula<3>& l4 = QuadratureManager::instance().getCubeGaussFormula(4);
+ const QuadratureFormula<3>& l5 = QuadratureManager::instance().getCubeGaussFormula(5);
- REQUIRE(same_formulas(l4, l5));
+ REQUIRE(&l4 == &l5);
REQUIRE(l5.numberOfPoints() == 14);
@@ -275,10 +253,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- CubeGaussQuadrature l6(6);
- CubeGaussQuadrature l7(7);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getCubeGaussFormula(6);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getCubeGaussFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 34);
@@ -297,10 +275,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 8 and 9")
{
- CubeGaussQuadrature l8(8);
- CubeGaussQuadrature l9(9);
+ const QuadratureFormula<3>& l8 = QuadratureManager::instance().getCubeGaussFormula(8);
+ const QuadratureFormula<3>& l9 = QuadratureManager::instance().getCubeGaussFormula(9);
- REQUIRE(same_formulas(l8, l9));
+ REQUIRE(&l8 == &l9);
REQUIRE(l9.numberOfPoints() == 58);
@@ -321,10 +299,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 10 and 11")
{
- CubeGaussQuadrature l10(10);
- CubeGaussQuadrature l11(11);
+ const QuadratureFormula<3>& l10 = QuadratureManager::instance().getCubeGaussFormula(10);
+ const QuadratureFormula<3>& l11 = QuadratureManager::instance().getCubeGaussFormula(11);
- REQUIRE(same_formulas(l10, l11));
+ REQUIRE(&l10 == &l11);
REQUIRE(l11.numberOfPoints() == 90);
@@ -347,10 +325,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 12 and 13")
{
- CubeGaussQuadrature l12(12);
- CubeGaussQuadrature l13(13);
+ const QuadratureFormula<3>& l12 = QuadratureManager::instance().getCubeGaussFormula(12);
+ const QuadratureFormula<3>& l13 = QuadratureManager::instance().getCubeGaussFormula(13);
- REQUIRE(same_formulas(l12, l13));
+ REQUIRE(&l12 == &l13);
REQUIRE(l13.numberOfPoints() == 154);
@@ -375,10 +353,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 14 and 15")
{
- CubeGaussQuadrature l14(14);
- CubeGaussQuadrature l15(15);
+ const QuadratureFormula<3>& l14 = QuadratureManager::instance().getCubeGaussFormula(14);
+ const QuadratureFormula<3>& l15 = QuadratureManager::instance().getCubeGaussFormula(15);
- REQUIRE(same_formulas(l14, l15));
+ REQUIRE(&l14 == &l15);
REQUIRE(l15.numberOfPoints() == 256);
@@ -405,10 +383,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 16 and 17")
{
- CubeGaussQuadrature l16(16);
- CubeGaussQuadrature l17(17);
+ const QuadratureFormula<3>& l16 = QuadratureManager::instance().getCubeGaussFormula(16);
+ const QuadratureFormula<3>& l17 = QuadratureManager::instance().getCubeGaussFormula(17);
- REQUIRE(same_formulas(l16, l17));
+ REQUIRE(&l16 == &l17);
REQUIRE(l17.numberOfPoints() == 346);
@@ -437,10 +415,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 18 and 19")
{
- CubeGaussQuadrature l18(18);
- CubeGaussQuadrature l19(19);
+ const QuadratureFormula<3>& l18 = QuadratureManager::instance().getCubeGaussFormula(18);
+ const QuadratureFormula<3>& l19 = QuadratureManager::instance().getCubeGaussFormula(19);
- REQUIRE(same_formulas(l18, l19));
+ REQUIRE(&l18 == &l19);
REQUIRE(l19.numberOfPoints() == 454);
@@ -471,10 +449,10 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
SECTION("degree 20 and 21")
{
- CubeGaussQuadrature l20(20);
- CubeGaussQuadrature l21(21);
+ const QuadratureFormula<3>& l20 = QuadratureManager::instance().getCubeGaussFormula(20);
+ const QuadratureFormula<3>& l21 = QuadratureManager::instance().getCubeGaussFormula(21);
- REQUIRE(same_formulas(l20, l21));
+ REQUIRE(&l20 == &l21);
REQUIRE(l21.numberOfPoints() == 580);
@@ -506,14 +484,15 @@ TEST_CASE("CubeGaussQuadrature", "[analysis]")
Catch::Approx(22));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(CubeGaussQuadrature(22), "error: Gauss quadrature formulae handle degrees up to 21 on cubes");
+ REQUIRE(QuadratureManager::instance().maxCubeGaussDegree() == CubeGaussQuadrature::max_degree);
}
SECTION("Access functions")
{
- CubeGaussQuadrature quadrature_formula(7);
+ const QuadratureFormula<3>& quadrature_formula = QuadratureManager::instance().getCubeGaussFormula(7);
+
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp
index 9a8839f5aada5171b4071d7cd787c1ddece8118c..5702f2a2a1da295b853c086bdb342d0106d62fce 100644
--- a/tests/test_FunctionEvaluation.cpp
+++ b/tests/test_FunctionEvaluation.cpp
@@ -12,7 +12,7 @@
TEST_CASE("IntegrationTools", "[integration]")
{
- IntegrationTools<2> integrationtools(QuadratureType::gausslobatto);
+ IntegrationTools<2> integrationtools(QuadratureType::GaussLobatto);
double a = 1;
double b = 2;
@@ -50,7 +50,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<5> integrationtools5(QuadratureType::gausslobatto);
+ IntegrationTools<5> integrationtools5(QuadratureType::GaussLobatto);
SECTION("TestOrder5")
{
double result = integrationtools5.testIntegrate(a, b);
@@ -60,7 +60,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<7> integrationtools7(QuadratureType::gausslobatto);
+ IntegrationTools<7> integrationtools7(QuadratureType::GaussLobatto);
SECTION("TestOrder7")
{
double result = integrationtools7.testIntegrate(a, b);
@@ -70,7 +70,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<9> integrationtools9(QuadratureType::gausslobatto);
+ IntegrationTools<9> integrationtools9(QuadratureType::GaussLobatto);
SECTION("TestOrder9")
{
double result = integrationtools9.testIntegrate(a, b);
@@ -80,7 +80,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<11> integrationtools11(QuadratureType::gausslobatto);
+ IntegrationTools<11> integrationtools11(QuadratureType::GaussLobatto);
SECTION("TestOrder11")
{
double result = integrationtools11.testIntegrate(a, b);
@@ -91,7 +91,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<2> integrationtoolslegendre(QuadratureType::gausslegendre);
+ IntegrationTools<2> integrationtoolslegendre(QuadratureType::GaussLegendre);
SECTION("Order-Legendre")
{
size_t order = integrationtools.order();
@@ -107,7 +107,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<5> integrationtoolslegendre5(QuadratureType::gausslegendre);
+ IntegrationTools<5> integrationtoolslegendre5(QuadratureType::GaussLegendre);
SECTION("TestOrder5-Legendre")
{
double result = integrationtoolslegendre5.testIntegrate(a, b);
@@ -117,7 +117,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<7> integrationtoolslegendre7(QuadratureType::gausslegendre);
+ IntegrationTools<7> integrationtoolslegendre7(QuadratureType::GaussLegendre);
SECTION("TestOrder7-Legendre")
{
double result = integrationtoolslegendre7.testIntegrate(a, b);
@@ -127,7 +127,7 @@ TEST_CASE("IntegrationTools", "[integration]")
REQUIRE(results[1] == 0);
REQUIRE(results[2] == Catch::Approx((-0.25 * (std::pow(2, 4) - std::pow(1, 4)))).epsilon(1E-14));
}
- IntegrationTools<9> integrationtoolslegendre9(QuadratureType::gausslegendre);
+ IntegrationTools<9> integrationtoolslegendre9(QuadratureType::GaussLegendre);
SECTION("TestOrder9-Legendre")
{
double result = integrationtoolslegendre9.testIntegrate(a, b);
diff --git a/tests/test_PrismGaussQuadrature.cpp b/tests/test_PrismGaussQuadrature.cpp
index c94f63c2d851b6299171fd651b111d8cdc397647..eb76c14e31349310616f7a9ad37e8063ec6d50e0 100644
--- a/tests/test_PrismGaussQuadrature.cpp
+++ b/tests/test_PrismGaussQuadrature.cpp
@@ -3,6 +3,7 @@
#include <catch2/matchers/catch_matchers_all.hpp>
#include <analysis/PrismGaussQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/AffineTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -211,7 +212,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 1")
{
- PrismGaussQuadrature l1(1);
+ const QuadratureFormula<3>& l1 = QuadratureManager::instance().getPrismGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -224,7 +225,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 2")
{
- PrismGaussQuadrature l2(2);
+ const QuadratureFormula<3>& l2 = QuadratureManager::instance().getPrismGaussFormula(2);
REQUIRE(l2.numberOfPoints() == 5);
@@ -239,7 +240,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 3")
{
- PrismGaussQuadrature l3(3);
+ const QuadratureFormula<3>& l3 = QuadratureManager::instance().getPrismGaussFormula(3);
REQUIRE(l3.numberOfPoints() == 8);
@@ -254,7 +255,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 4")
{
- PrismGaussQuadrature l4(4);
+ const QuadratureFormula<3>& l4 = QuadratureManager::instance().getPrismGaussFormula(4);
REQUIRE(l4.numberOfPoints() == 11);
@@ -270,7 +271,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 5")
{
- PrismGaussQuadrature l5(5);
+ const QuadratureFormula<3>& l5 = QuadratureManager::instance().getPrismGaussFormula(5);
REQUIRE(l5.numberOfPoints() == 16);
@@ -287,7 +288,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 6")
{
- PrismGaussQuadrature l6(6);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getPrismGaussFormula(6);
REQUIRE(l6.numberOfPoints() == 28);
@@ -306,7 +307,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 7")
{
- PrismGaussQuadrature l7(7);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getPrismGaussFormula(7);
REQUIRE(l7.numberOfPoints() == 35);
@@ -325,7 +326,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 8")
{
- PrismGaussQuadrature l8(8);
+ const QuadratureFormula<3>& l8 = QuadratureManager::instance().getPrismGaussFormula(8);
REQUIRE(l8.numberOfPoints() == 46);
@@ -346,7 +347,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 9")
{
- PrismGaussQuadrature l9(9);
+ const QuadratureFormula<3>& l9 = QuadratureManager::instance().getPrismGaussFormula(9);
REQUIRE(l9.numberOfPoints() == 59);
@@ -367,7 +368,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 10")
{
- PrismGaussQuadrature l10(10);
+ const QuadratureFormula<3>& l10 = QuadratureManager::instance().getPrismGaussFormula(10);
REQUIRE(l10.numberOfPoints() == 84);
@@ -390,7 +391,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 11")
{
- PrismGaussQuadrature l11(11);
+ const QuadratureFormula<3>& l11 = QuadratureManager::instance().getPrismGaussFormula(11);
REQUIRE(l11.numberOfPoints() == 99);
@@ -413,7 +414,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 12")
{
- PrismGaussQuadrature l12(12);
+ const QuadratureFormula<3>& l12 = QuadratureManager::instance().getPrismGaussFormula(12);
REQUIRE(l12.numberOfPoints() == 136);
@@ -438,7 +439,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 13")
{
- PrismGaussQuadrature l13(13);
+ const QuadratureFormula<3>& l13 = QuadratureManager::instance().getPrismGaussFormula(13);
REQUIRE(l13.numberOfPoints() == 162);
@@ -463,7 +464,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 14")
{
- PrismGaussQuadrature l14(14);
+ const QuadratureFormula<3>& l14 = QuadratureManager::instance().getPrismGaussFormula(14);
REQUIRE(l14.numberOfPoints() == 194);
@@ -491,7 +492,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 15")
{
- PrismGaussQuadrature l15(15);
+ const QuadratureFormula<3>& l15 = QuadratureManager::instance().getPrismGaussFormula(15);
REQUIRE(l15.numberOfPoints() == 238);
@@ -518,7 +519,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 16")
{
- PrismGaussQuadrature l16(16);
+ const QuadratureFormula<3>& l16 = QuadratureManager::instance().getPrismGaussFormula(16);
REQUIRE(l16.numberOfPoints() == 287);
@@ -548,7 +549,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 17")
{
- PrismGaussQuadrature l17(17);
+ const QuadratureFormula<3>& l17 = QuadratureManager::instance().getPrismGaussFormula(17);
REQUIRE(l17.numberOfPoints() == 338);
@@ -578,7 +579,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 18")
{
- PrismGaussQuadrature l18(18);
+ const QuadratureFormula<3>& l18 = QuadratureManager::instance().getPrismGaussFormula(18);
REQUIRE(l18.numberOfPoints() == 396);
@@ -611,7 +612,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 19")
{
- PrismGaussQuadrature l19(19);
+ const QuadratureFormula<3>& l19 = QuadratureManager::instance().getPrismGaussFormula(19);
REQUIRE(l19.numberOfPoints() == 420);
@@ -644,7 +645,7 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
SECTION("degree 20")
{
- PrismGaussQuadrature l20(20);
+ const QuadratureFormula<3>& l20 = QuadratureManager::instance().getPrismGaussFormula(20);
REQUIRE(l20.numberOfPoints() == 518);
@@ -677,8 +678,8 @@ TEST_CASE("PrismGaussQuadrature", "[analysis]")
Catch::Approx(22).margin(0.01));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(PrismGaussQuadrature(21), "error: Gauss quadrature formulae handle degrees up to 20 on prisms");
+ REQUIRE(QuadratureManager::instance().maxPrismGaussDegree() == PrismGaussQuadrature::max_degree);
}
}
diff --git a/tests/test_PrismTransformation.cpp b/tests/test_PrismTransformation.cpp
index c0e6618f0bf944a0a5a567e29a8ed7a82c3e3dd4..95f48a4bfc502d513fdbe055de8a35220c9b573f 100644
--- a/tests/test_PrismTransformation.cpp
+++ b/tests/test_PrismTransformation.cpp
@@ -1,7 +1,7 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-#include <analysis/PrismGaussQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/PrismTransformation.hpp>
// clazy:excludeall=non-pod-global-static
@@ -71,7 +71,8 @@ TEST_CASE("PrismTransformation", "[geometry]")
{
// due to the z component of the jacobian determinant, degree 3
// polynomials must be exactly integrated
- PrismGaussQuadrature gauss(3);
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPrismGaussFormula(3);
+
double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
@@ -92,7 +93,8 @@ TEST_CASE("PrismTransformation", "[geometry]")
};
// 5 is the minimum quadrature rule to integrate the polynomial on the prism
- PrismGaussQuadrature gauss(5);
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPrismGaussFormula(5);
+
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
diff --git a/tests/test_PyramidGaussQuadrature.cpp b/tests/test_PyramidGaussQuadrature.cpp
index 8cb76bcb776d85533fdfde62bb3a2bddd1b0d71e..38c04a3d6e48acc5fd8a604f4c19b619d838aed9 100644
--- a/tests/test_PyramidGaussQuadrature.cpp
+++ b/tests/test_PyramidGaussQuadrature.cpp
@@ -3,6 +3,7 @@
#include <catch2/matchers/catch_matchers_all.hpp>
#include <analysis/PyramidGaussQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/AffineTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -152,7 +153,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 1")
{
- PyramidGaussQuadrature l1(1);
+ const QuadratureFormula<3>& l1 = QuadratureManager::instance().getPyramidGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -163,7 +164,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 2")
{
- PyramidGaussQuadrature l2(2);
+ const QuadratureFormula<3>& l2 = QuadratureManager::instance().getPyramidGaussFormula(2);
REQUIRE(l2.numberOfPoints() == 5);
@@ -175,7 +176,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 3")
{
- PyramidGaussQuadrature l3(3);
+ const QuadratureFormula<3>& l3 = QuadratureManager::instance().getPyramidGaussFormula(3);
REQUIRE(l3.numberOfPoints() == 6);
@@ -188,7 +189,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 4")
{
- PyramidGaussQuadrature l4(4);
+ const QuadratureFormula<3>& l4 = QuadratureManager::instance().getPyramidGaussFormula(4);
REQUIRE(l4.numberOfPoints() == 10);
@@ -202,7 +203,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 5")
{
- PyramidGaussQuadrature l5(5);
+ const QuadratureFormula<3>& l5 = QuadratureManager::instance().getPyramidGaussFormula(5);
REQUIRE(l5.numberOfPoints() == 15);
@@ -217,7 +218,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 6")
{
- PyramidGaussQuadrature l6(6);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getPyramidGaussFormula(6);
REQUIRE(l6.numberOfPoints() == 23);
@@ -233,7 +234,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 7")
{
- PyramidGaussQuadrature l7(7);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getPyramidGaussFormula(7);
REQUIRE(l7.numberOfPoints() == 31);
@@ -250,7 +251,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 8")
{
- PyramidGaussQuadrature l8(8);
+ const QuadratureFormula<3>& l8 = QuadratureManager::instance().getPyramidGaussFormula(8);
REQUIRE(l8.numberOfPoints() == 47);
@@ -268,7 +269,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 9")
{
- PyramidGaussQuadrature l9(9);
+ const QuadratureFormula<3>& l9 = QuadratureManager::instance().getPyramidGaussFormula(9);
REQUIRE(l9.numberOfPoints() == 62);
@@ -287,7 +288,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 10")
{
- PyramidGaussQuadrature l10(10);
+ const QuadratureFormula<3>& l10 = QuadratureManager::instance().getPyramidGaussFormula(10);
REQUIRE(l10.numberOfPoints() == 80);
@@ -307,7 +308,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 11")
{
- PyramidGaussQuadrature l11(11);
+ const QuadratureFormula<3>& l11 = QuadratureManager::instance().getPyramidGaussFormula(11);
REQUIRE(l11.numberOfPoints() == 103);
@@ -328,7 +329,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 12")
{
- PyramidGaussQuadrature l12(12);
+ const QuadratureFormula<3>& l12 = QuadratureManager::instance().getPyramidGaussFormula(12);
REQUIRE(l12.numberOfPoints() == 127);
@@ -350,7 +351,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 13")
{
- PyramidGaussQuadrature l13(13);
+ const QuadratureFormula<3>& l13 = QuadratureManager::instance().getPyramidGaussFormula(13);
REQUIRE(l13.numberOfPoints() == 152);
@@ -373,7 +374,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 14")
{
- PyramidGaussQuadrature l14(14);
+ const QuadratureFormula<3>& l14 = QuadratureManager::instance().getPyramidGaussFormula(14);
REQUIRE(l14.numberOfPoints() == 184);
@@ -397,7 +398,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 15")
{
- PyramidGaussQuadrature l15(15);
+ const QuadratureFormula<3>& l15 = QuadratureManager::instance().getPyramidGaussFormula(15);
REQUIRE(l15.numberOfPoints() == 234);
@@ -422,7 +423,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 16")
{
- PyramidGaussQuadrature l16(16);
+ const QuadratureFormula<3>& l16 = QuadratureManager::instance().getPyramidGaussFormula(16);
REQUIRE(l16.numberOfPoints() == 285);
@@ -448,7 +449,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 17")
{
- PyramidGaussQuadrature l17(17);
+ const QuadratureFormula<3>& l17 = QuadratureManager::instance().getPyramidGaussFormula(17);
REQUIRE(l17.numberOfPoints() == 319);
@@ -475,7 +476,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 18")
{
- PyramidGaussQuadrature l18(18);
+ const QuadratureFormula<3>& l18 = QuadratureManager::instance().getPyramidGaussFormula(18);
REQUIRE(l18.numberOfPoints() == 357);
@@ -503,7 +504,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 19")
{
- PyramidGaussQuadrature l19(19);
+ const QuadratureFormula<3>& l19 = QuadratureManager::instance().getPyramidGaussFormula(19);
REQUIRE(l19.numberOfPoints() == 418);
@@ -532,7 +533,7 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
SECTION("degree 20")
{
- PyramidGaussQuadrature l20(20);
+ const QuadratureFormula<3>& l20 = QuadratureManager::instance().getPyramidGaussFormula(20);
REQUIRE(l20.numberOfPoints() == 489);
@@ -560,9 +561,8 @@ TEST_CASE("PyramidGaussQuadrature", "[analysis]")
REQUIRE(integrate(p21, l20) != Catch::Approx(796248143552124247176376796357. / 110883848060635650000000000000000.));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(PyramidGaussQuadrature(21),
- "error: Gauss quadrature formulae handle degrees up to 20 on pyramids");
+ REQUIRE(QuadratureManager::instance().maxPyramidGaussDegree() == PyramidGaussQuadrature::max_degree);
}
}
diff --git a/tests/test_PyramidTransformation.cpp b/tests/test_PyramidTransformation.cpp
index bfdcee769d97558a2ff8037752b8e060fb3a8c8b..38fb5fc1eea7918b157e7aede0324cc3f243890f 100644
--- a/tests/test_PyramidTransformation.cpp
+++ b/tests/test_PyramidTransformation.cpp
@@ -1,7 +1,7 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-#include <analysis/PyramidGaussQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/PyramidTransformation.hpp>
// clazy:excludeall=non-pod-global-static
@@ -62,7 +62,8 @@ TEST_CASE("PyramidTransformation", "[geometry]")
SECTION("volume calculation")
{
// The jacobian determinant is a degree 1 polynomial
- PyramidGaussQuadrature gauss(1);
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPyramidGaussFormula(1);
+
double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
@@ -83,8 +84,8 @@ TEST_CASE("PyramidTransformation", "[geometry]")
};
// 3 is the minimum quadrature rule to integrate the polynomial on the pyramid
- PyramidGaussQuadrature gauss(3);
- double integral = 0;
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getPyramidGaussFormula(3);
+ double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
}
diff --git a/tests/test_Q1Transformation.cpp b/tests/test_Q1Transformation.cpp
index 680b53e22c8f7ed83e1af660ad4dbe1161143c06..5649129330e3108bb18c77dbe48a9a32c0548591 100644
--- a/tests/test_Q1Transformation.cpp
+++ b/tests/test_Q1Transformation.cpp
@@ -1,7 +1,7 @@
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
-#include <analysis/TensorialGaussLegendreQuadrature.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/Q1Transformation.hpp>
// clazy:excludeall=non-pod-global-static
@@ -89,7 +89,8 @@ TEST_CASE("Q1Transformation", "[geometry]")
SECTION("Gauss order 1")
{
// One point is enough in 2d
- TensorialGaussLegendreQuadrature<2> gauss(1);
+ const QuadratureFormula<2>& gauss = QuadratureManager::instance().getSquareGaussLegendreFormula(1);
+
double surface = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
surface += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
@@ -110,7 +111,8 @@ TEST_CASE("Q1Transformation", "[geometry]")
// Jacbian determinant is a degree 1 polynomial, so the
// following formula is required to reach exactness
- TensorialGaussLegendreQuadrature<2> gauss(3);
+ const QuadratureFormula<2>& gauss = QuadratureManager::instance().getSquareGaussLegendreFormula(3);
+
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
@@ -194,7 +196,8 @@ TEST_CASE("Q1Transformation", "[geometry]")
SECTION("Gauss order 3")
{
// The jacobian determinant is of maximal degree 2 in each variable
- TensorialGaussLegendreQuadrature<3> gauss(2);
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeGaussLegendreFormula(2);
+
double volume = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
volume += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i));
@@ -216,7 +219,8 @@ TEST_CASE("Q1Transformation", "[geometry]")
// Jacbian determinant is a degree 2 polynomial, so the
// following formula is required to reach exactness
- TensorialGaussLegendreQuadrature<3> gauss(4);
+ const QuadratureFormula<3>& gauss = QuadratureManager::instance().getCubeGaussLegendreFormula(4);
+
double integral = 0;
for (size_t i = 0; i < gauss.numberOfPoints(); ++i) {
integral += gauss.weight(i) * t.jacobianDeterminant(gauss.point(i)) * p(gauss.point(i));
diff --git a/tests/test_SquareGaussQuadrature.cpp b/tests/test_SquareGaussQuadrature.cpp
index ee6d8d6d6346f647f7361ae6b7d4533677f812b0..1d2f799124d9296f51a7f093bc8e24c243efc5e4 100644
--- a/tests/test_SquareGaussQuadrature.cpp
+++ b/tests/test_SquareGaussQuadrature.cpp
@@ -4,6 +4,7 @@
#include <algebra/TinyMatrix.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <analysis/SquareGaussQuadrature.hpp>
#include <utils/Exceptions.hpp>
@@ -11,29 +12,6 @@
TEST_CASE("SquareGaussQuadrature", "[analysis]")
{
- auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
- auto point_list0 = quadrature_formula0.pointList();
- auto weight_list0 = quadrature_formula0.weightList();
-
- auto point_list1 = quadrature_formula1.pointList();
- auto weight_list1 = quadrature_formula1.weightList();
-
- bool is_same = true;
-
- for (size_t i = 0; i < point_list0.size(); ++i) {
- if (point_list0[i] != point_list1[i]) {
- return false;
- }
- }
- for (size_t i = 0; i < weight_list0.size(); ++i) {
- if (weight_list0[i] != weight_list1[i]) {
- return false;
- }
- }
-
- return is_same;
- };
-
auto integrate = [](auto f, auto quadrature_formula) {
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
@@ -195,7 +173,7 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 0 and 1")
{
- SquareGaussQuadrature l1(1);
+ const QuadratureFormula<2>& l1 = QuadratureManager::instance().getSquareGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -208,10 +186,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 2 and 3")
{
- SquareGaussQuadrature l2(2);
- SquareGaussQuadrature l3(3);
+ const QuadratureFormula<2>& l2 = QuadratureManager::instance().getSquareGaussFormula(2);
+ const QuadratureFormula<2>& l3 = QuadratureManager::instance().getSquareGaussFormula(3);
- REQUIRE(same_formulas(l2, l3));
+ REQUIRE(&l2 == &l3);
REQUIRE(l3.numberOfPoints() == 4);
@@ -226,10 +204,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 4 and 5")
{
- SquareGaussQuadrature l4(4);
- SquareGaussQuadrature l5(5);
+ const QuadratureFormula<2>& l4 = QuadratureManager::instance().getSquareGaussFormula(4);
+ const QuadratureFormula<2>& l5 = QuadratureManager::instance().getSquareGaussFormula(5);
- REQUIRE(same_formulas(l4, l5));
+ REQUIRE(&l4 == &l5);
REQUIRE(l5.numberOfPoints() == 8);
@@ -246,10 +224,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- SquareGaussQuadrature l6(6);
- SquareGaussQuadrature l7(7);
+ const QuadratureFormula<2>& l6 = QuadratureManager::instance().getSquareGaussFormula(6);
+ const QuadratureFormula<2>& l7 = QuadratureManager::instance().getSquareGaussFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 12);
@@ -268,10 +246,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 8 and 9")
{
- SquareGaussQuadrature l8(8);
- SquareGaussQuadrature l9(9);
+ const QuadratureFormula<2>& l8 = QuadratureManager::instance().getSquareGaussFormula(8);
+ const QuadratureFormula<2>& l9 = QuadratureManager::instance().getSquareGaussFormula(9);
- REQUIRE(same_formulas(l8, l9));
+ REQUIRE(&l8 == &l9);
REQUIRE(l9.numberOfPoints() == 20);
@@ -292,10 +270,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 10 and 11")
{
- SquareGaussQuadrature l10(10);
- SquareGaussQuadrature l11(11);
+ const QuadratureFormula<2>& l10 = QuadratureManager::instance().getSquareGaussFormula(10);
+ const QuadratureFormula<2>& l11 = QuadratureManager::instance().getSquareGaussFormula(11);
- REQUIRE(same_formulas(l10, l11));
+ REQUIRE(&l10 == &l11);
REQUIRE(l11.numberOfPoints() == 28);
@@ -318,10 +296,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 12 and 13")
{
- SquareGaussQuadrature l12(12);
- SquareGaussQuadrature l13(13);
+ const QuadratureFormula<2>& l12 = QuadratureManager::instance().getSquareGaussFormula(12);
+ const QuadratureFormula<2>& l13 = QuadratureManager::instance().getSquareGaussFormula(13);
- REQUIRE(same_formulas(l12, l13));
+ REQUIRE(&l12 == &l13);
REQUIRE(l13.numberOfPoints() == 37);
@@ -346,10 +324,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 14 and 15")
{
- SquareGaussQuadrature l14(14);
- SquareGaussQuadrature l15(15);
+ const QuadratureFormula<2>& l14 = QuadratureManager::instance().getSquareGaussFormula(14);
+ const QuadratureFormula<2>& l15 = QuadratureManager::instance().getSquareGaussFormula(15);
- REQUIRE(same_formulas(l14, l15));
+ REQUIRE(&l14 == &l15);
REQUIRE(l15.numberOfPoints() == 48);
@@ -376,10 +354,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 16 and 17")
{
- SquareGaussQuadrature l16(16);
- SquareGaussQuadrature l17(17);
+ const QuadratureFormula<2>& l16 = QuadratureManager::instance().getSquareGaussFormula(16);
+ const QuadratureFormula<2>& l17 = QuadratureManager::instance().getSquareGaussFormula(17);
- REQUIRE(same_formulas(l16, l17));
+ REQUIRE(&l16 == &l17);
REQUIRE(l17.numberOfPoints() == 60);
@@ -408,10 +386,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 18 and 19")
{
- SquareGaussQuadrature l18(18);
- SquareGaussQuadrature l19(19);
+ const QuadratureFormula<2>& l18 = QuadratureManager::instance().getSquareGaussFormula(18);
+ const QuadratureFormula<2>& l19 = QuadratureManager::instance().getSquareGaussFormula(19);
- REQUIRE(same_formulas(l18, l19));
+ REQUIRE(&l18 == &l19);
REQUIRE(l19.numberOfPoints() == 72);
@@ -442,10 +420,10 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
SECTION("degree 20 and 21")
{
- SquareGaussQuadrature l20(20);
- SquareGaussQuadrature l21(21);
+ const QuadratureFormula<2>& l20 = QuadratureManager::instance().getSquareGaussFormula(20);
+ const QuadratureFormula<2>& l21 = QuadratureManager::instance().getSquareGaussFormula(21);
- REQUIRE(same_formulas(l20, l21));
+ REQUIRE(&l20 == &l21);
REQUIRE(l21.numberOfPoints() == 85);
@@ -477,15 +455,15 @@ TEST_CASE("SquareGaussQuadrature", "[analysis]")
Catch::Approx(22).margin(0.01));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(SquareGaussQuadrature(22),
- "error: Gauss quadrature formulae handle degrees up to 21 on squares");
+ REQUIRE(QuadratureManager::instance().maxSquareGaussDegree() == SquareGaussQuadrature::max_degree);
}
SECTION("Access functions")
{
- SquareGaussQuadrature quadrature_formula(7);
+ const QuadratureFormula<2>& quadrature_formula = QuadratureManager::instance().getSquareGaussFormula(7);
+
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
diff --git a/tests/test_TensorialGaussLegendreQuadrature.cpp b/tests/test_TensorialGaussLegendreQuadrature.cpp
index 8fc97fdae479a13eb3625f6dc1629089bd667222..0f132c700f6402c48589993cadcf1d73174eb697 100644
--- a/tests/test_TensorialGaussLegendreQuadrature.cpp
+++ b/tests/test_TensorialGaussLegendreQuadrature.cpp
@@ -2,6 +2,7 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <analysis/TensorialGaussLegendreQuadrature.hpp>
#include <utils/Exceptions.hpp>
@@ -9,29 +10,6 @@
TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
{
- auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
- auto point_list0 = quadrature_formula0.pointList();
- auto weight_list0 = quadrature_formula0.weightList();
-
- auto point_list1 = quadrature_formula1.pointList();
- auto weight_list1 = quadrature_formula1.weightList();
-
- bool is_same = true;
-
- for (size_t i = 0; i < point_list0.size(); ++i) {
- if (point_list0[i] != point_list1[i]) {
- return false;
- }
- }
- for (size_t i = 0; i < weight_list0.size(); ++i) {
- if (weight_list0[i] != weight_list1[i]) {
- return false;
- }
- }
-
- return is_same;
- };
-
SECTION("1D")
{
auto is_symmetric_formula = [](auto quadrature_formula) {
@@ -179,12 +157,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
return 25 * std::pow(x, 24) + p23(X);
};
- SECTION("degree 0 and 1")
+ SECTION("degree 1")
{
- TensorialGaussLegendreQuadrature<1> l0(0);
- TensorialGaussLegendreQuadrature<1> l1(1);
+ const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineGaussLegendreFormula(1);
- REQUIRE(same_formulas(l0, l1));
REQUIRE(is_symmetric_formula(l1));
REQUIRE(l1.numberOfPoints() == 1);
@@ -198,10 +174,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 2 and 3")
{
- TensorialGaussLegendreQuadrature<1> l2(2);
- TensorialGaussLegendreQuadrature<1> l3(3);
+ const QuadratureFormula<1>& l2 = QuadratureManager::instance().getLineGaussLegendreFormula(2);
+ const QuadratureFormula<1>& l3 = QuadratureManager::instance().getLineGaussLegendreFormula(3);
- REQUIRE(same_formulas(l2, l3));
+ REQUIRE(&l2 == &l3);
REQUIRE(is_symmetric_formula(l3));
REQUIRE(l3.numberOfPoints() == 2);
@@ -217,10 +193,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 4 and 5")
{
- TensorialGaussLegendreQuadrature<1> l4(4);
- TensorialGaussLegendreQuadrature<1> l5(5);
+ const QuadratureFormula<1>& l4 = QuadratureManager::instance().getLineGaussLegendreFormula(4);
+ const QuadratureFormula<1>& l5 = QuadratureManager::instance().getLineGaussLegendreFormula(5);
- REQUIRE(same_formulas(l4, l5));
+ REQUIRE(&l4 == &l5);
REQUIRE(is_symmetric_formula(l5));
REQUIRE(l5.numberOfPoints() == 3);
@@ -238,10 +214,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLegendreQuadrature<1> l6(6);
- TensorialGaussLegendreQuadrature<1> l7(7);
+ const QuadratureFormula<1>& l6 = QuadratureManager::instance().getLineGaussLegendreFormula(6);
+ const QuadratureFormula<1>& l7 = QuadratureManager::instance().getLineGaussLegendreFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(is_symmetric_formula(l7));
REQUIRE(l7.numberOfPoints() == 4);
@@ -261,10 +237,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 8 and 9")
{
- TensorialGaussLegendreQuadrature<1> l8(8);
- TensorialGaussLegendreQuadrature<1> l9(9);
+ const QuadratureFormula<1>& l8 = QuadratureManager::instance().getLineGaussLegendreFormula(8);
+ const QuadratureFormula<1>& l9 = QuadratureManager::instance().getLineGaussLegendreFormula(9);
- REQUIRE(same_formulas(l8, l9));
+ REQUIRE(&l8 == &l9);
REQUIRE(is_symmetric_formula(l9));
REQUIRE(l9.numberOfPoints() == 5);
@@ -286,10 +262,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 10 and 11")
{
- TensorialGaussLegendreQuadrature<1> l10(10);
- TensorialGaussLegendreQuadrature<1> l11(11);
+ const QuadratureFormula<1>& l10 = QuadratureManager::instance().getLineGaussLegendreFormula(10);
+ const QuadratureFormula<1>& l11 = QuadratureManager::instance().getLineGaussLegendreFormula(11);
- REQUIRE(same_formulas(l10, l11));
+ REQUIRE(&l10 == &l11);
REQUIRE(is_symmetric_formula(l11));
REQUIRE(l11.numberOfPoints() == 6);
@@ -313,10 +289,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 12 and 13")
{
- TensorialGaussLegendreQuadrature<1> l12(12);
- TensorialGaussLegendreQuadrature<1> l13(13);
+ const QuadratureFormula<1>& l12 = QuadratureManager::instance().getLineGaussLegendreFormula(12);
+ const QuadratureFormula<1>& l13 = QuadratureManager::instance().getLineGaussLegendreFormula(13);
- REQUIRE(same_formulas(l12, l13));
+ REQUIRE(&l12 == &l13);
REQUIRE(is_symmetric_formula(l13));
REQUIRE(l13.numberOfPoints() == 7);
@@ -342,10 +318,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 14 and 15")
{
- TensorialGaussLegendreQuadrature<1> l14(14);
- TensorialGaussLegendreQuadrature<1> l15(15);
+ const QuadratureFormula<1>& l14 = QuadratureManager::instance().getLineGaussLegendreFormula(14);
+ const QuadratureFormula<1>& l15 = QuadratureManager::instance().getLineGaussLegendreFormula(15);
- REQUIRE(same_formulas(l14, l15));
+ REQUIRE(&l14 == &l15);
REQUIRE(is_symmetric_formula(l15));
REQUIRE(l15.numberOfPoints() == 8);
@@ -373,10 +349,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 16 and 17")
{
- TensorialGaussLegendreQuadrature<1> l16(16);
- TensorialGaussLegendreQuadrature<1> l17(17);
+ const QuadratureFormula<1>& l16 = QuadratureManager::instance().getLineGaussLegendreFormula(16);
+ const QuadratureFormula<1>& l17 = QuadratureManager::instance().getLineGaussLegendreFormula(17);
- REQUIRE(same_formulas(l16, l17));
+ REQUIRE(&l16 == &l17);
REQUIRE(is_symmetric_formula(l17));
REQUIRE(l17.numberOfPoints() == 9);
@@ -406,10 +382,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 18 and 19")
{
- TensorialGaussLegendreQuadrature<1> l18(18);
- TensorialGaussLegendreQuadrature<1> l19(19);
+ const QuadratureFormula<1>& l18 = QuadratureManager::instance().getLineGaussLegendreFormula(18);
+ const QuadratureFormula<1>& l19 = QuadratureManager::instance().getLineGaussLegendreFormula(19);
- REQUIRE(same_formulas(l18, l19));
+ REQUIRE(&l18 == &l19);
REQUIRE(is_symmetric_formula(l19));
REQUIRE(l19.numberOfPoints() == 10);
@@ -441,10 +417,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 20 and 21")
{
- TensorialGaussLegendreQuadrature<1> l20(20);
- TensorialGaussLegendreQuadrature<1> l21(21);
+ const QuadratureFormula<1>& l20 = QuadratureManager::instance().getLineGaussLegendreFormula(20);
+ const QuadratureFormula<1>& l21 = QuadratureManager::instance().getLineGaussLegendreFormula(21);
- REQUIRE(same_formulas(l20, l21));
+ REQUIRE(&l20 == &l21);
REQUIRE(is_symmetric_formula(l21));
REQUIRE(l21.numberOfPoints() == 11);
@@ -478,10 +454,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 22 and 23")
{
- TensorialGaussLegendreQuadrature<1> l22(22);
- TensorialGaussLegendreQuadrature<1> l23(23);
+ const QuadratureFormula<1>& l22 = QuadratureManager::instance().getLineGaussLegendreFormula(22);
+ const QuadratureFormula<1>& l23 = QuadratureManager::instance().getLineGaussLegendreFormula(23);
- REQUIRE(same_formulas(l22, l23));
+ REQUIRE(&l22 == &l23);
REQUIRE(is_symmetric_formula(l23));
REQUIRE(l23.numberOfPoints() == 12);
@@ -515,10 +491,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
REQUIRE(get_order(p24, l23, -1, 1, 26) == Catch::Approx(24).margin(0.05));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(TensorialGaussLegendreQuadrature<1>(24),
- "error: Gauss-Legendre quadrature formulae handle degrees up to 23");
+ REQUIRE(QuadratureManager::instance().maxGaussLegendreDegree() ==
+ TensorialGaussLegendreQuadrature<1>::max_degree);
}
}
@@ -558,10 +534,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLegendreQuadrature<2> l6(6);
- TensorialGaussLegendreQuadrature<2> l7(7);
+ const QuadratureFormula<2>& l6 = QuadratureManager::instance().getSquareGaussLegendreFormula(6);
+ const QuadratureFormula<2>& l7 = QuadratureManager::instance().getSquareGaussLegendreFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 4 * 4);
@@ -609,10 +585,10 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLegendreQuadrature<3> l6(6);
- TensorialGaussLegendreQuadrature<3> l7(7);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getCubeGaussLegendreFormula(6);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getCubeGaussLegendreFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 4 * 4 * 4);
@@ -621,19 +597,4 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
(0.8 - 0.2) * (std::pow(0.7, 8) - std::pow(-0.1, 8))));
}
}
-
- SECTION("Access functions")
- {
- TensorialGaussLegendreQuadrature<3> quadrature_formula(7);
- auto point_list = quadrature_formula.pointList();
- auto weight_list = quadrature_formula.weightList();
-
- REQUIRE(point_list.size() == quadrature_formula.numberOfPoints());
- REQUIRE(weight_list.size() == quadrature_formula.numberOfPoints());
-
- for (size_t i = 0; i < quadrature_formula.numberOfPoints(); ++i) {
- REQUIRE(&point_list[i] == &quadrature_formula.point(i));
- REQUIRE(&weight_list[i] == &quadrature_formula.weight(i));
- }
- }
}
diff --git a/tests/test_TensorialGaussLobattoQuadrature.cpp b/tests/test_TensorialGaussLobattoQuadrature.cpp
index 069073eea67a67c815b935188ed9142752033539..7744bb123f0e345b5d01211307fdc2be063b2b81 100644
--- a/tests/test_TensorialGaussLobattoQuadrature.cpp
+++ b/tests/test_TensorialGaussLobattoQuadrature.cpp
@@ -2,6 +2,7 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <analysis/TensorialGaussLobattoQuadrature.hpp>
#include <utils/Exceptions.hpp>
@@ -9,29 +10,6 @@
TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
{
- auto same_formulas = [](auto quadrature_formula0, auto quadrature_formula1) {
- auto point_list0 = quadrature_formula0.pointList();
- auto weight_list0 = quadrature_formula0.weightList();
-
- auto point_list1 = quadrature_formula1.pointList();
- auto weight_list1 = quadrature_formula1.weightList();
-
- bool is_same = true;
-
- for (size_t i = 0; i < point_list0.size(); ++i) {
- if (point_list0[i] != point_list1[i]) {
- return false;
- }
- }
- for (size_t i = 0; i < weight_list0.size(); ++i) {
- if (weight_list0[i] != weight_list1[i]) {
- return false;
- }
- }
-
- return is_same;
- };
-
SECTION("1D")
{
auto is_symmetric_formula = [](auto quadrature_formula) {
@@ -139,12 +117,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
return 15 * std::pow(x, 14) + p13(X);
};
- SECTION("degree 0 and 1")
+ SECTION("degree 1")
{
- TensorialGaussLobattoQuadrature<1> l0(0);
- TensorialGaussLobattoQuadrature<1> l1(1);
+ const QuadratureFormula<1>& l1 = QuadratureManager::instance().getLineGaussLobattoFormula(1);
- REQUIRE(same_formulas(l0, l1));
REQUIRE(is_symmetric_formula(l1));
REQUIRE(l1.numberOfPoints() == 2);
@@ -158,10 +134,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 2 and 3")
{
- TensorialGaussLobattoQuadrature<1> l2(2);
- TensorialGaussLobattoQuadrature<1> l3(3);
+ const QuadratureFormula<1>& l2 = QuadratureManager::instance().getLineGaussLobattoFormula(2);
+ const QuadratureFormula<1>& l3 = QuadratureManager::instance().getLineGaussLobattoFormula(3);
- REQUIRE(same_formulas(l2, l3));
+ REQUIRE(&l2 == &l3);
REQUIRE(is_symmetric_formula(l3));
REQUIRE(l3.numberOfPoints() == 3);
@@ -177,10 +153,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 4 and 5")
{
- TensorialGaussLobattoQuadrature<1> l4(4);
- TensorialGaussLobattoQuadrature<1> l5(5);
+ const QuadratureFormula<1>& l4 = QuadratureManager::instance().getLineGaussLobattoFormula(4);
+ const QuadratureFormula<1>& l5 = QuadratureManager::instance().getLineGaussLobattoFormula(5);
- REQUIRE(same_formulas(l4, l5));
+ REQUIRE(&l4 == &l5);
REQUIRE(is_symmetric_formula(l5));
REQUIRE(l5.numberOfPoints() == 4);
@@ -198,10 +174,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLobattoQuadrature<1> l6(6);
- TensorialGaussLobattoQuadrature<1> l7(7);
+ const QuadratureFormula<1>& l6 = QuadratureManager::instance().getLineGaussLobattoFormula(6);
+ const QuadratureFormula<1>& l7 = QuadratureManager::instance().getLineGaussLobattoFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(is_symmetric_formula(l7));
REQUIRE(l7.numberOfPoints() == 5);
@@ -221,10 +197,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 8 and 9")
{
- TensorialGaussLobattoQuadrature<1> l8(8);
- TensorialGaussLobattoQuadrature<1> l9(9);
+ const QuadratureFormula<1>& l8 = QuadratureManager::instance().getLineGaussLobattoFormula(8);
+ const QuadratureFormula<1>& l9 = QuadratureManager::instance().getLineGaussLobattoFormula(9);
- REQUIRE(same_formulas(l8, l9));
+ REQUIRE(&l8 == &l9);
REQUIRE(is_symmetric_formula(l9));
REQUIRE(l9.numberOfPoints() == 6);
@@ -246,10 +222,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 10 and 11")
{
- TensorialGaussLobattoQuadrature<1> l10(10);
- TensorialGaussLobattoQuadrature<1> l11(11);
+ const QuadratureFormula<1>& l10 = QuadratureManager::instance().getLineGaussLobattoFormula(10);
+ const QuadratureFormula<1>& l11 = QuadratureManager::instance().getLineGaussLobattoFormula(11);
- REQUIRE(same_formulas(l10, l11));
+ REQUIRE(&l10 == &l11);
REQUIRE(is_symmetric_formula(l11));
REQUIRE(l11.numberOfPoints() == 7);
@@ -273,10 +249,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 12 and 13")
{
- TensorialGaussLobattoQuadrature<1> l12(12);
- TensorialGaussLobattoQuadrature<1> l13(13);
+ const QuadratureFormula<1>& l12 = QuadratureManager::instance().getLineGaussLobattoFormula(12);
+ const QuadratureFormula<1>& l13 = QuadratureManager::instance().getLineGaussLobattoFormula(13);
- REQUIRE(same_formulas(l12, l13));
+ REQUIRE(&l12 == &l13);
REQUIRE(is_symmetric_formula(l13));
REQUIRE(l13.numberOfPoints() == 8);
@@ -300,10 +276,9 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
REQUIRE(get_order(p14, l13, -1, 1, 16) == Catch::Approx(14));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(TensorialGaussLobattoQuadrature<1>(14),
- "error: Gauss-Lobatto quadrature formulae handle degrees up to 13");
+ REQUIRE(QuadratureManager::instance().maxGaussLobattoDegree() == TensorialGaussLobattoQuadrature<1>::max_degree);
}
}
@@ -343,10 +318,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLobattoQuadrature<2> l6(6);
- TensorialGaussLobattoQuadrature<2> l7(7);
+ const QuadratureFormula<2>& l6 = QuadratureManager::instance().getSquareGaussLobattoFormula(6);
+ const QuadratureFormula<2>& l7 = QuadratureManager::instance().getSquareGaussLobattoFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 5 * 5);
@@ -394,10 +369,10 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("degree 6 and 7")
{
- TensorialGaussLobattoQuadrature<3> l6(6);
- TensorialGaussLobattoQuadrature<3> l7(7);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getCubeGaussLobattoFormula(6);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getCubeGaussLobattoFormula(7);
- REQUIRE(same_formulas(l6, l7));
+ REQUIRE(&l6 == &l7);
REQUIRE(l7.numberOfPoints() == 5 * 5 * 5);
@@ -409,7 +384,8 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
SECTION("Access functions")
{
- TensorialGaussLobattoQuadrature<3> quadrature_formula(7);
+ const QuadratureFormula<3>& quadrature_formula = QuadratureManager::instance().getCubeGaussLobattoFormula(7);
+
auto point_list = quadrature_formula.pointList();
auto weight_list = quadrature_formula.weightList();
diff --git a/tests/test_TetrahedronGaussQuadrature.cpp b/tests/test_TetrahedronGaussQuadrature.cpp
index 22250323132dbd5b788627be835232d462f32e75..de17a655cd21a19eab099f7f5cc06a1d6b37c1ab 100644
--- a/tests/test_TetrahedronGaussQuadrature.cpp
+++ b/tests/test_TetrahedronGaussQuadrature.cpp
@@ -2,6 +2,7 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <analysis/TetrahedronGaussQuadrature.hpp>
#include <geometry/AffineTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -204,7 +205,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 1")
{
- TetrahedronGaussQuadrature l1(1);
+ const QuadratureFormula<3>& l1 = QuadratureManager::instance().getTetrahedronGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -217,7 +218,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 2")
{
- TetrahedronGaussQuadrature l2(2);
+ const QuadratureFormula<3>& l2 = QuadratureManager::instance().getTetrahedronGaussFormula(2);
REQUIRE(l2.numberOfPoints() == 4);
@@ -231,7 +232,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 3")
{
- TetrahedronGaussQuadrature l3(3);
+ const QuadratureFormula<3>& l3 = QuadratureManager::instance().getTetrahedronGaussFormula(3);
REQUIRE(l3.numberOfPoints() == 8);
@@ -244,7 +245,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 4")
{
- TetrahedronGaussQuadrature l4(4);
+ const QuadratureFormula<3>& l4 = QuadratureManager::instance().getTetrahedronGaussFormula(4);
REQUIRE(l4.numberOfPoints() == 14);
@@ -260,7 +261,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 5")
{
- TetrahedronGaussQuadrature l5(5);
+ const QuadratureFormula<3>& l5 = QuadratureManager::instance().getTetrahedronGaussFormula(5);
REQUIRE(l5.numberOfPoints() == 14);
@@ -275,7 +276,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 6")
{
- TetrahedronGaussQuadrature l6(6);
+ const QuadratureFormula<3>& l6 = QuadratureManager::instance().getTetrahedronGaussFormula(6);
REQUIRE(l6.numberOfPoints() == 24);
@@ -293,7 +294,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 7")
{
- TetrahedronGaussQuadrature l7(7);
+ const QuadratureFormula<3>& l7 = QuadratureManager::instance().getTetrahedronGaussFormula(7);
REQUIRE(l7.numberOfPoints() == 35);
@@ -312,7 +313,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 8")
{
- TetrahedronGaussQuadrature l8(8);
+ const QuadratureFormula<3>& l8 = QuadratureManager::instance().getTetrahedronGaussFormula(8);
REQUIRE(l8.numberOfPoints() == 46);
@@ -333,7 +334,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 9")
{
- TetrahedronGaussQuadrature l9(9);
+ const QuadratureFormula<3>& l9 = QuadratureManager::instance().getTetrahedronGaussFormula(9);
REQUIRE(l9.numberOfPoints() == 59);
@@ -354,7 +355,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 10")
{
- TetrahedronGaussQuadrature l10(10);
+ const QuadratureFormula<3>& l10 = QuadratureManager::instance().getTetrahedronGaussFormula(10);
REQUIRE(l10.numberOfPoints() == 81);
@@ -377,7 +378,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 11")
{
- TetrahedronGaussQuadrature l11(11);
+ const QuadratureFormula<3>& l11 = QuadratureManager::instance().getTetrahedronGaussFormula(11);
REQUIRE(l11.numberOfPoints() == 110);
@@ -400,7 +401,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 12")
{
- TetrahedronGaussQuadrature l12(12);
+ const QuadratureFormula<3>& l12 = QuadratureManager::instance().getTetrahedronGaussFormula(12);
REQUIRE(l12.numberOfPoints() == 168);
@@ -425,7 +426,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 13")
{
- TetrahedronGaussQuadrature l13(13);
+ const QuadratureFormula<3>& l13 = QuadratureManager::instance().getTetrahedronGaussFormula(13);
REQUIRE(l13.numberOfPoints() == 172);
@@ -450,7 +451,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 14")
{
- TetrahedronGaussQuadrature l14(14);
+ const QuadratureFormula<3>& l14 = QuadratureManager::instance().getTetrahedronGaussFormula(14);
REQUIRE(l14.numberOfPoints() == 204);
@@ -477,7 +478,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 15")
{
- TetrahedronGaussQuadrature l15(15);
+ const QuadratureFormula<3>& l15 = QuadratureManager::instance().getTetrahedronGaussFormula(15);
REQUIRE(l15.numberOfPoints() == 264);
@@ -505,7 +506,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 16")
{
- TetrahedronGaussQuadrature l16(16);
+ const QuadratureFormula<3>& l16 = QuadratureManager::instance().getTetrahedronGaussFormula(16);
REQUIRE(l16.numberOfPoints() == 304);
@@ -534,7 +535,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 17")
{
- TetrahedronGaussQuadrature l17(17);
+ const QuadratureFormula<3>& l17 = QuadratureManager::instance().getTetrahedronGaussFormula(17);
REQUIRE(l17.numberOfPoints() == 364);
@@ -564,7 +565,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 18")
{
- TetrahedronGaussQuadrature l18(18);
+ const QuadratureFormula<3>& l18 = QuadratureManager::instance().getTetrahedronGaussFormula(18);
REQUIRE(l18.numberOfPoints() == 436);
@@ -597,7 +598,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 19")
{
- TetrahedronGaussQuadrature l19(19);
+ const QuadratureFormula<3>& l19 = QuadratureManager::instance().getTetrahedronGaussFormula(19);
REQUIRE(l19.numberOfPoints() == 487);
@@ -629,7 +630,7 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
SECTION("degree 20")
{
- TetrahedronGaussQuadrature l20(20);
+ const QuadratureFormula<3>& l20 = QuadratureManager::instance().getTetrahedronGaussFormula(20);
REQUIRE(l20.numberOfPoints() == 552);
@@ -661,9 +662,8 @@ TEST_CASE("TetrahedronGaussQuadrature", "[analysis]")
Catch::Approx(22).margin(0.01));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(TetrahedronGaussQuadrature(21),
- "error: Gauss quadrature formulae handle degrees up to 20 on tetrahedra");
+ REQUIRE(QuadratureManager::instance().maxTetrahedronGaussDegree() == TetrahedronGaussQuadrature::max_degree);
}
}
diff --git a/tests/test_TriangleGaussQuadrature.cpp b/tests/test_TriangleGaussQuadrature.cpp
index a02a8f7285f162c416c38e9dcb80fb54a730884d..4fa7604ce8b39b7f64097c2291a5d427780595da 100644
--- a/tests/test_TriangleGaussQuadrature.cpp
+++ b/tests/test_TriangleGaussQuadrature.cpp
@@ -2,6 +2,7 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <analysis/TriangleGaussQuadrature.hpp>
#include <geometry/AffineTransformation.hpp>
#include <utils/Exceptions.hpp>
@@ -165,7 +166,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 1")
{
- TriangleGaussQuadrature l1(1);
+ const QuadratureFormula<2>& l1 = QuadratureManager::instance().getTriangleGaussFormula(1);
REQUIRE(l1.numberOfPoints() == 1);
@@ -178,7 +179,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 2")
{
- TriangleGaussQuadrature l2(2);
+ const QuadratureFormula<2>& l2 = QuadratureManager::instance().getTriangleGaussFormula(2);
REQUIRE(l2.numberOfPoints() == 3);
@@ -193,8 +194,8 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 3 and 4")
{
- TriangleGaussQuadrature l3(3);
- TriangleGaussQuadrature l4(4);
+ const QuadratureFormula<2>& l3 = QuadratureManager::instance().getTriangleGaussFormula(3);
+ const QuadratureFormula<2>& l4 = QuadratureManager::instance().getTriangleGaussFormula(4);
bool is_same = true;
is_same &= (l3.numberOfPoints() == l4.numberOfPoints());
@@ -220,7 +221,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 5")
{
- TriangleGaussQuadrature l5(5);
+ const QuadratureFormula<2>& l5 = QuadratureManager::instance().getTriangleGaussFormula(5);
REQUIRE(l5.numberOfPoints() == 7);
@@ -237,7 +238,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 6")
{
- TriangleGaussQuadrature l6(6);
+ const QuadratureFormula<2>& l6 = QuadratureManager::instance().getTriangleGaussFormula(6);
REQUIRE(l6.numberOfPoints() == 12);
@@ -255,7 +256,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 7")
{
- TriangleGaussQuadrature l7(7);
+ const QuadratureFormula<2>& l7 = QuadratureManager::instance().getTriangleGaussFormula(7);
REQUIRE(l7.numberOfPoints() == 15);
@@ -274,7 +275,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 8")
{
- TriangleGaussQuadrature l8(8);
+ const QuadratureFormula<2>& l8 = QuadratureManager::instance().getTriangleGaussFormula(8);
REQUIRE(l8.numberOfPoints() == 16);
@@ -295,7 +296,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 9")
{
- TriangleGaussQuadrature l9(9);
+ const QuadratureFormula<2>& l9 = QuadratureManager::instance().getTriangleGaussFormula(9);
REQUIRE(l9.numberOfPoints() == 19);
@@ -316,7 +317,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 10")
{
- TriangleGaussQuadrature l10(10);
+ const QuadratureFormula<2>& l10 = QuadratureManager::instance().getTriangleGaussFormula(10);
REQUIRE(l10.numberOfPoints() == 25);
@@ -339,7 +340,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 11")
{
- TriangleGaussQuadrature l11(11);
+ const QuadratureFormula<2>& l11 = QuadratureManager::instance().getTriangleGaussFormula(11);
REQUIRE(l11.numberOfPoints() == 28);
@@ -362,7 +363,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 12")
{
- TriangleGaussQuadrature l12(12);
+ const QuadratureFormula<2>& l12 = QuadratureManager::instance().getTriangleGaussFormula(12);
REQUIRE(l12.numberOfPoints() == 33);
@@ -387,7 +388,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 13")
{
- TriangleGaussQuadrature l13(13);
+ const QuadratureFormula<2>& l13 = QuadratureManager::instance().getTriangleGaussFormula(13);
REQUIRE(l13.numberOfPoints() == 37);
@@ -412,7 +413,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 14")
{
- TriangleGaussQuadrature l14(14);
+ const QuadratureFormula<2>& l14 = QuadratureManager::instance().getTriangleGaussFormula(14);
REQUIRE(l14.numberOfPoints() == 42);
@@ -439,7 +440,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 15")
{
- TriangleGaussQuadrature l15(15);
+ const QuadratureFormula<2>& l15 = QuadratureManager::instance().getTriangleGaussFormula(15);
REQUIRE(l15.numberOfPoints() == 49);
@@ -466,7 +467,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 16")
{
- TriangleGaussQuadrature l16(16);
+ const QuadratureFormula<2>& l16 = QuadratureManager::instance().getTriangleGaussFormula(16);
REQUIRE(l16.numberOfPoints() == 55);
@@ -495,7 +496,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 17")
{
- TriangleGaussQuadrature l17(17);
+ const QuadratureFormula<2>& l17 = QuadratureManager::instance().getTriangleGaussFormula(17);
REQUIRE(l17.numberOfPoints() == 60);
@@ -524,7 +525,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 18")
{
- TriangleGaussQuadrature l18(18);
+ const QuadratureFormula<2>& l18 = QuadratureManager::instance().getTriangleGaussFormula(18);
REQUIRE(l18.numberOfPoints() == 67);
@@ -555,7 +556,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 19")
{
- TriangleGaussQuadrature l19(19);
+ const QuadratureFormula<2>& l19 = QuadratureManager::instance().getTriangleGaussFormula(19);
REQUIRE(l19.numberOfPoints() == 73);
@@ -586,7 +587,7 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
SECTION("degree 20")
{
- TriangleGaussQuadrature l20(20);
+ const QuadratureFormula<2>& l20 = QuadratureManager::instance().getTriangleGaussFormula(20);
REQUIRE(l20.numberOfPoints() == 79);
@@ -617,9 +618,8 @@ TEST_CASE("TriangleGaussQuadrature", "[analysis]")
REQUIRE(get_order(p21, l20, -520205676970121316953. / 215685489550781250000000.) == Catch::Approx(22).margin(0.5));
}
- SECTION("not implemented formulae")
+ SECTION("max implemented degree")
{
- REQUIRE_THROWS_WITH(TriangleGaussQuadrature(21),
- "error: Gauss quadrature formulae handle degrees up to 20 on triangles");
+ REQUIRE(QuadratureManager::instance().maxTriangleGaussDegree() == TriangleGaussQuadrature::max_degree);
}
}
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index 9c03ba9213ece8947688036ca7d5554c8e38a279..fa207c91b9b743ffa1986c20a3cdc6e1622c0a1a 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -2,6 +2,7 @@
#include <Kokkos_Core.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <language/utils/OperatorRepository.hpp>
#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/DiamondDualMeshManager.hpp>
@@ -36,6 +37,7 @@ main(int argc, char* argv[])
SynchronizerManager::create();
RandomEngine::create();
+ QuadratureManager::create();
MeshDataManager::create();
DiamondDualConnectivityManager::create();
DiamondDualMeshManager::create();
@@ -53,6 +55,7 @@ main(int argc, char* argv[])
DiamondDualMeshManager::destroy();
DiamondDualConnectivityManager::destroy();
MeshDataManager::destroy();
+ QuadratureManager::destroy();
RandomEngine::destroy();
SynchronizerManager::destroy();
}