diff --git a/MathLib/Integration/GaussLegendrePyramid.cpp b/MathLib/Integration/GaussLegendrePyramid.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f769a73abcb0998d88625c8e146bcd9ef898fd7
--- /dev/null
+++ b/MathLib/Integration/GaussLegendrePyramid.cpp
@@ -0,0 +1,60 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "GaussLegendrePyramid.h"
+
+namespace MathLib
+{
+
+template <> const std::array<std::array<double, 3>, GaussLegendrePyramid<1>::NPoints>
+GaussLegendrePyramid<1>::X = {{{{1./3., 1./3., 1./3.}}}};
+template <> double const GaussLegendrePyramid<1>::W[1] = {1.0};
+
+const std::array<std::array<double, 3>, GaussLegendrePyramid<2>::NPoints>
+GaussLegendrePyramid<2>::X = {{ {{-0.584237394672177188, -0.584237394672177188, -2./3.}},
+ {{ 0.584237394672177188, -0.584237394672177188, -2./3.}},
+ {{ 0.584237394672177188, 0.584237394672177188, -2./3.}},
+ {{-0.584237394672177188, 0.584237394672177188, -2./3.}},
+ {{0., 0., 2./5.}}
+ }};
+double const GaussLegendrePyramid<2>::W[5] = {81./100., 81./100., 81./100., 81./100., 125./27.};
+
+const std::array<std::array<double, 3>, GaussLegendrePyramid<3>::NPoints>
+GaussLegendrePyramid<3>::X = {{
+ {{-0.673931986207731726, -0.673931986207731726, -0.142857142857142857}},
+ {{ 0.673931986207731726, -0.673931986207731726, -0.142857142857142857}},
+ {{ 0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
+ {{-0.673931986207731726, 0.673931986207731726, -0.142857142857142857}},
+ {{-0.610639618865075532, 0.0, -0.321428571428571429}},
+ {{ 0.610639618865075532, 0.0, -0.321428571428571429}},
+ {{ 0.0, -0.610639618865075532, -0.321428571428571429}},
+ {{ 0.0, 0.610639618865075532, -0.321428571428571429}},
+ {{ 0.0, 0.0, 0.524394036075370072}},
+ {{-0.580939660561084423, -0.580939660561084423, -0.830065359477124183}},
+ {{ 0.580939660561084423, -0.580939660561084423, -0.830065359477124183}},
+ {{ 0.580939660561084423, 0.580939660561084423, -0.830065359477124183}},
+ {{-0.580939660561084423, 0.580939660561084423, -0.830065359477124183}}
+ }};
+double const GaussLegendrePyramid<3>::W[13] = {
+ 0.515003019323671498,
+ 0.515003019323671498,
+ 0.515003019323671498,
+ 0.515003019323671498,
+ 0.257183745242064659,
+ 0.257183745242064659,
+ 0.257183745242064659,
+ 0.257183745242064659,
+ 2.474004977113405936,
+ 0.419515737191525950,
+ 0.419515737191525950,
+ 0.419515737191525950,
+ 0.419515737191525950
+ };
+
+}
diff --git a/MathLib/Integration/GaussLegendrePyramid.h b/MathLib/Integration/GaussLegendrePyramid.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f826b7b3c48c91172247d44300d3db316b02458
--- /dev/null
+++ b/MathLib/Integration/GaussLegendrePyramid.h
@@ -0,0 +1,47 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef GAUSSLEGENDREPYRAMID_H_
+#define GAUSSLEGENDREPYRAMID_H_
+
+#include <array>
+
+namespace MathLib
+{
+
+/// Gauss-Legendre quadrature on pyramid
+///
+/// \tparam ORDER integration order.
+template <unsigned ORDER>
+struct GaussLegendrePyramid {
+ static const unsigned Order = ORDER;
+ static const unsigned NPoints = ORDER;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendrePyramid<2> {
+ static const unsigned Order = 2;
+ static const unsigned NPoints = 5;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendrePyramid<3> {
+ static const unsigned Order = 3;
+ static const unsigned NPoints = 13;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+}
+
+#endif //GAUSSLEGENDREPYRAMID_H_
diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ad827751cf406b340520f272d52a9d5dc367df7
--- /dev/null
+++ b/MathLib/Integration/GaussLegendreTet.cpp
@@ -0,0 +1,60 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "GaussLegendreTet.h"
+
+namespace MathLib
+{
+
+template <> const std::array<std::array<double, 3>, GaussLegendreTet<1>::NPoints>
+GaussLegendreTet<1>::X = {{{{1./3., 1./3., 1./3.}}}};
+template <> double const GaussLegendreTet<1>::W[1] = {1.0};
+
+const std::array<std::array<double, 3>, GaussLegendreTet<2>::NPoints>
+GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
+ {{1./6., 1./6., 1./6.}},
+ {{1./2., 1./6., 1./6.}},
+ {{1./6., 1./2., 1./6.}},
+ {{1./6., 1./6., 1./2.}} }};
+double const GaussLegendreTet<2>::W[5] = {-0.133333333333333, 0.075, 0.075, 0.075, 0.075};
+
+const std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints>
+GaussLegendreTet<3>::X = {{ {{1./2., 1./2., 1./2.}},
+ {{0.09197107805272303, 0.09197107805272303, 0.09197107805272303}},
+ {{0.72408676584183096, 0.09197107805272303, 0.09197107805272303}},
+ {{0.09197107805272303, 0.72408676584183096, 0.09197107805272303}},
+ {{0.09197107805272303, 0.09197107805272303, 0.72408676584183096}},
+ {{0.44364916731037080, 0.05635083268962915, 0.05635083268962915}},
+ {{0.05635083268962915, 0.44364916731037080, 0.05635083268962915}},
+ {{0.05635083268962915, 0.05635083268962915, 0.44364916731037080}},
+ {{0.05635083268962915, 0.44364916731037080, 0.44364916731037080}},
+ {{0.44364916731037080, 0.05635083268962915, 0.44364916731037080}},
+ {{0.44364916731037080, 0.44364916731037080, 0.05635083268962915}},
+ {{0.31979362782962989, 0.31979362782962989, 0.31979362782962989}},
+ {{0.04061911651111023, 0.31979362782962989, 0.31979362782962989}},
+ {{0.31979362782962989, 0.04061911651111023, 0.31979362782962989}},
+ {{0.31979362782962989, 0.31979362782962989, 0.04061911651111023}} }};
+double const GaussLegendreTet<3>::W[15] = {0.019753086419753086,
+ 0.011989513963169772,
+ 0.011989513963169772,
+ 0.011989513963169772,
+ 0.011989513963169772,
+ 0.008818342151675485,
+ 0.008818342151675485,
+ 0.008818342151675485,
+ 0.008818342151675485,
+ 0.008818342151675485,
+ 0.008818342151675485,
+ 0.011511367871045397,
+ 0.011511367871045397,
+ 0.011511367871045397,
+ 0.011511367871045397
+ };
+
+}
diff --git a/MathLib/Integration/GaussLegendreTet.h b/MathLib/Integration/GaussLegendreTet.h
new file mode 100644
index 0000000000000000000000000000000000000000..11b2e670730e7788bd07e4147f72710d9e08d48c
--- /dev/null
+++ b/MathLib/Integration/GaussLegendreTet.h
@@ -0,0 +1,47 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef GAUSSLEGENDRETET_H_
+#define GAUSSLEGENDRETET_H_
+
+#include <array>
+
+namespace MathLib
+{
+
+/// Gauss-Legendre quadrature on tetrahedrals
+///
+/// \tparam ORDER integration order.
+template <unsigned ORDER>
+struct GaussLegendreTet {
+ static const unsigned Order = ORDER;
+ static const unsigned NPoints = ORDER;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendreTet<2> {
+ static const unsigned Order = 2;
+ static const unsigned NPoints = 5;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendreTet<3> {
+ static const unsigned Order = 3;
+ static const unsigned NPoints = 15;
+ static const std::array<std::array<double, 3>, NPoints> X;
+ static const double W[NPoints];
+};
+
+}
+
+#endif //GAUSSLEGENDRETET_H_
diff --git a/NumLib/Fem/Integration/GaussIntegrationPolicy.h b/NumLib/Fem/Integration/GaussIntegrationPolicy.h
index 5d424d97e219fff15b38aa9acc5315dbefc558d0..496c91c5c59f59c59812a7cf746a8c4dae7efda0 100644
--- a/NumLib/Fem/Integration/GaussIntegrationPolicy.h
+++ b/NumLib/Fem/Integration/GaussIntegrationPolicy.h
@@ -12,8 +12,14 @@
#include "NumLib/Fem/Integration/IntegrationGaussRegular.h"
#include "NumLib/Fem/Integration/IntegrationGaussTri.h"
+#include "NumLib/Fem/Integration/IntegrationGaussTet.h"
+#include "NumLib/Fem/Integration/IntegrationGaussPrism.h"
+#include "NumLib/Fem/Integration/IntegrationGaussPyramid.h"
#include "MeshLib/Elements/Tri.h"
+#include "MeshLib/Elements/Tet.h"
+#include "MeshLib/Elements/Prism.h"
+#include "MeshLib/Elements/Pyramid.h"
namespace NumLib
{
@@ -41,6 +47,55 @@ struct GaussIntegrationPolicy<MeshLib::Tri>
using IntegrationMethod = NumLib::IntegrationGaussTri;
};
+template<>
+struct GaussIntegrationPolicy<MeshLib::Tri6>
+{
+ using MeshElement = MeshLib::Tri6;
+ using IntegrationMethod = NumLib::IntegrationGaussTri;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Tet>
+{
+ using MeshElement = MeshLib::Tri;
+ using IntegrationMethod = NumLib::IntegrationGaussTet;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Tet10>
+{
+ using MeshElement = MeshLib::Tet10;
+ using IntegrationMethod = NumLib::IntegrationGaussTet;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Prism>
+{
+ using MeshElement = MeshLib::Prism;
+ using IntegrationMethod = NumLib::IntegrationGaussPrism;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Prism15>
+{
+ using MeshElement = MeshLib::Prism15;
+ using IntegrationMethod = NumLib::IntegrationGaussPrism;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Pyramid>
+{
+ using MeshElement = MeshLib::Pyramid;
+ using IntegrationMethod = NumLib::IntegrationGaussPyramid;
+};
+
+template<>
+struct GaussIntegrationPolicy<MeshLib::Pyramid13>
+{
+ using MeshElement = MeshLib::Pyramid13;
+ using IntegrationMethod = NumLib::IntegrationGaussPyramid;
+};
+
} // namespace NumLib
#endif // NUMLIB_GAUSSINTEGRATIONPOLICY_H_
diff --git a/NumLib/Fem/Integration/IntegrationGaussPrism.h b/NumLib/Fem/Integration/IntegrationGaussPrism.h
new file mode 100644
index 0000000000000000000000000000000000000000..97f01b65b4bae81771eeced3995c7e7904f55f31
--- /dev/null
+++ b/NumLib/Fem/Integration/IntegrationGaussPrism.h
@@ -0,0 +1,110 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef INTEGRATIONGAUSSPRISM_H_
+#define INTEGRATIONGAUSSPRISM_H_
+
+#include "MathLib/Integration/GaussLegendre.h"
+#include "MathLib/Integration/GaussLegendreTri.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Gauss quadrature rule for prisms
+ */
+class IntegrationGaussPrism
+{
+ typedef MathLib::TemplateWeightedPoint<double, double, 3>
+ WeightedPoint;
+public:
+ /**
+ * Construct this object with the given integration order
+ *
+ * @param order integration order (default 2)
+ */
+ explicit IntegrationGaussPrism(std::size_t order = 2)
+ : _order(2), _n_sampl_pt(0)
+ {
+ this->setIntegrationOrder(order);
+ }
+
+ /// Change the integration order.
+ void setIntegrationOrder(std::size_t /*order*/)
+ {
+ _order = 2; // fixed
+ _n_sampl_pt = getNPoints(_order);
+ }
+
+ /// return current integration order.
+ std::size_t getIntegrationOrder() const {return _order;}
+
+ /// return the number of sampling points
+ std::size_t getNPoints() const {return _n_sampl_pt;}
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param igp The integration point index
+ * @return a weighted point
+ */
+ WeightedPoint getWeightedPoint(std::size_t igp)
+ {
+ return getWeightedPoint(getIntegrationOrder(), igp);
+ }
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param order the number of integration points
+ * @param pt_id the sampling point id
+ * @return weight
+ */
+ static WeightedPoint
+ getWeightedPoint(std::size_t /*order*/, std::size_t igp)
+ {
+ const unsigned gp_r = igp % 3;
+ const unsigned gp_t = (unsigned)(igp/3);
+ std::array<double, 3> rst;
+ rst[0] = MathLib::GaussLegendreTri<2>::X[gp_r][0];
+ rst[1] = MathLib::GaussLegendreTri<2>::X[gp_r][1];
+ rst[2] = MathLib::GaussLegendre<2>::X[gp_t];
+ double w = MathLib::GaussLegendreTri<2>::W[gp_r] * 0.5 * MathLib::GaussLegendre<2>::W[gp_t];
+ return WeightedPoint(rst, w);
+ }
+
+ template <typename Method>
+ static WeightedPoint
+ getWeightedPoint(std::size_t igp)
+ {
+ return WeightedPoint(Method::X[igp], Method::W[igp]);
+ }
+
+
+ /**
+ * get the number of integration points
+ *
+ * @param order the number of integration points
+ * @return the number of points
+ */
+ static std::size_t
+ getNPoints(std::size_t order)
+ {
+ if (order==2) return 6;
+ return 0;
+ }
+
+private:
+ std::size_t _order;
+ std::size_t _n_sampl_pt;
+};
+
+}
+
+#endif //INTEGRATIONGAUSSPRISM_H_
diff --git a/NumLib/Fem/Integration/IntegrationGaussPyramid.h b/NumLib/Fem/Integration/IntegrationGaussPyramid.h
new file mode 100644
index 0000000000000000000000000000000000000000..45f3b1482ea4e5d036751186dd7b1f1bded7e264
--- /dev/null
+++ b/NumLib/Fem/Integration/IntegrationGaussPyramid.h
@@ -0,0 +1,113 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef INTEGRATIONGAUSSPYRAMID_H_
+#define INTEGRATIONGAUSSPYRAMID_H_
+
+#include "MathLib/Integration/GaussLegendrePyramid.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Gauss quadrature rule for pyramid
+ */
+class IntegrationGaussPyramid
+{
+ typedef MathLib::TemplateWeightedPoint<double, double, 3>
+ WeightedPoint;
+public:
+ /**
+ * Construct this object with the given integration order
+ *
+ * @param order integration order (default 2)
+ */
+ explicit IntegrationGaussPyramid(std::size_t order = 2)
+ : _order(order), _n_sampl_pt(0)
+ {
+ this->setIntegrationOrder(order);
+ }
+
+ /// Change the integration order.
+ void setIntegrationOrder(std::size_t order)
+ {
+ _order = order;
+ _n_sampl_pt = getNPoints(_order);
+ }
+
+ /// return current integration order.
+ std::size_t getIntegrationOrder() const {return _order;}
+
+ /// return the number of sampling points
+ std::size_t getNPoints() const {return _n_sampl_pt;}
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param igp The integration point index
+ * @return a weighted point
+ */
+ WeightedPoint getWeightedPoint(std::size_t igp)
+ {
+ return getWeightedPoint(getIntegrationOrder(), igp);
+ }
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param order the number of integration points
+ * @param pt_id the sampling point id
+ * @return weight
+ */
+ static WeightedPoint
+ getWeightedPoint(std::size_t order, std::size_t igp)
+ {
+ switch (order)
+ {
+ case 1: return getWeightedPoint<MathLib::GaussLegendrePyramid<1> >(igp);
+ case 2: return getWeightedPoint<MathLib::GaussLegendrePyramid<2> >(igp);
+ case 3: return getWeightedPoint<MathLib::GaussLegendrePyramid<3> >(igp);
+ }
+ return WeightedPoint(std::array<double, 3>(), 0);
+ }
+
+ template <typename Method>
+ static WeightedPoint
+ getWeightedPoint(std::size_t igp)
+ {
+ return WeightedPoint(Method::X[igp], Method::W[igp]);
+ }
+
+
+ /**
+ * get the number of integration points
+ *
+ * @param order the number of integration points
+ * @return the number of points
+ */
+ static std::size_t
+ getNPoints(std::size_t order)
+ {
+ switch (order)
+ {
+ case 1: return MathLib::GaussLegendrePyramid<1>::NPoints;
+ case 2: return MathLib::GaussLegendrePyramid<2>::NPoints;
+ case 3: return MathLib::GaussLegendrePyramid<3>::NPoints;
+ }
+ return 0;
+ }
+
+private:
+ std::size_t _order;
+ std::size_t _n_sampl_pt;
+};
+
+}
+
+#endif //INTEGRATIONGAUSSPYRAMID_H_
diff --git a/NumLib/Fem/Integration/IntegrationGaussTet.h b/NumLib/Fem/Integration/IntegrationGaussTet.h
new file mode 100644
index 0000000000000000000000000000000000000000..70372b6ccc3b745825b882ab90d69a15fe431ce0
--- /dev/null
+++ b/NumLib/Fem/Integration/IntegrationGaussTet.h
@@ -0,0 +1,113 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef INTEGRATIONGAUSSTET_H_
+#define INTEGRATIONGAUSSTET_H_
+
+#include "MathLib/Integration/GaussLegendreTet.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Gauss quadrature rule for tetrahedrals
+ */
+class IntegrationGaussTet
+{
+ typedef MathLib::TemplateWeightedPoint<double, double, 3>
+ WeightedPoint;
+public:
+ /**
+ * Construct this object with the given integration order
+ *
+ * @param order integration order (default 2)
+ */
+ explicit IntegrationGaussTet(std::size_t order = 2)
+ : _order(order), _n_sampl_pt(0)
+ {
+ this->setIntegrationOrder(order);
+ }
+
+ /// Change the integration order.
+ void setIntegrationOrder(std::size_t order)
+ {
+ _order = order;
+ _n_sampl_pt = getNPoints(order);
+ }
+
+ /// return current integration order.
+ std::size_t getIntegrationOrder() const {return _order;}
+
+ /// return the number of sampling points
+ std::size_t getNPoints() const {return _n_sampl_pt;}
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param igp The integration point index
+ * @return a weighted point
+ */
+ WeightedPoint getWeightedPoint(std::size_t igp)
+ {
+ return getWeightedPoint(getIntegrationOrder(), igp);
+ }
+
+ /**
+ * get coordinates of a integration point
+ *
+ * @param order the number of integration points
+ * @param pt_id the sampling point id
+ * @return weight
+ */
+ static WeightedPoint
+ getWeightedPoint(std::size_t order, std::size_t igp)
+ {
+ switch (order)
+ {
+ case 1: return getWeightedPoint<MathLib::GaussLegendreTet<1> >(igp);
+ case 2: return getWeightedPoint<MathLib::GaussLegendreTet<2> >(igp);
+ case 3: return getWeightedPoint<MathLib::GaussLegendreTet<3> >(igp);
+ }
+ return WeightedPoint(std::array<double, 3>(), 0);
+ }
+
+ template <typename Method>
+ static WeightedPoint
+ getWeightedPoint(std::size_t igp)
+ {
+ return WeightedPoint(Method::X[igp], Method::W[igp]);
+ }
+
+
+ /**
+ * get the number of integration points
+ *
+ * @param order the number of integration points
+ * @return the number of points
+ */
+ static std::size_t
+ getNPoints(std::size_t order)
+ {
+ switch (order)
+ {
+ case 1: return MathLib::GaussLegendreTet<1>::NPoints;
+ case 2: return MathLib::GaussLegendreTet<2>::NPoints;
+ case 3: return MathLib::GaussLegendreTet<3>::NPoints;
+ }
+ return 0;
+ }
+
+private:
+ std::size_t _order;
+ std::size_t _n_sampl_pt;
+};
+
+}
+
+#endif //INTEGRATIONGAUSSTET_H_