diff --git a/MathLib/Integration/GaussLegendreTri.cpp b/MathLib/Integration/GaussLegendreTri.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b62888c78f92fab60f8de07617c2af6ba997daa6
--- /dev/null
+++ b/MathLib/Integration/GaussLegendreTri.cpp
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 "GaussLegendreTri.h"
+
+namespace MathLib
+{
+
+template <> const std::array<std::array<double, 2>, GaussLegendreTri<1>::NPoints> GaussLegendreTri<1>::X = {{{{1./3., 1./3.}}}};
+template <> double const GaussLegendreTri<1>::W[1] = {1.0};
+
+const std::array<std::array<double, 2>, GaussLegendreTri<2>::NPoints> GaussLegendreTri<2>::X = {{{{1./6., 1./6.}}, {{2./3., 1./6.}}, {{1./6., 2./3.}}}};
+double const GaussLegendreTri<2>::W[3] = {1./3., 1./3., 1./3.};
+
+const std::array<std::array<double, 2>, GaussLegendreTri<3>::NPoints> GaussLegendreTri<3>::X = {{{{1./3., 1./3.}}, {{1./5., 3./5.}}, {{1./5., 1./5.}}, {{3./5., 1./5.}}}};
+double const GaussLegendreTri<3>::W[4] = {-27./48., 25./48., 25./48., 25./48.};
+
+}
diff --git a/MathLib/Integration/GaussLegendreTri.h b/MathLib/Integration/GaussLegendreTri.h
new file mode 100644
index 0000000000000000000000000000000000000000..4dea7a9869ca497079b086a7ae99ee067ad3f556
--- /dev/null
+++ b/MathLib/Integration/GaussLegendreTri.h
@@ -0,0 +1,52 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 GAUSSLEGENDRETRI_H_
+#define GAUSSLEGENDRETRI_H_
+
+#include <array>
+
+namespace MathLib
+{
+
+/// Gauss-Legendre quadrature on triangles
+///
+/// \tparam ORDER integration order.
+template <unsigned ORDER>
+struct GaussLegendreTri {
+ static const unsigned Order = ORDER;
+ static const unsigned NPoints = ORDER;
+ static const std::array<std::array<double, 2>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendreTri<2> {
+ static const unsigned Order = 2;
+ static const unsigned NPoints = 3;
+ static const std::array<std::array<double, 2>, NPoints> X;
+ static const double W[NPoints];
+};
+
+template <>
+struct GaussLegendreTri<3> {
+ static const unsigned Order = 3;
+ static const unsigned NPoints = 4;
+ static const std::array<std::array<double, 2>, NPoints> X;
+ static const double W[NPoints];
+};
+
+}
+
+#endif //GAUSSLEGENDRETRI_H_
diff --git a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
index 68fb04d38e5a5ea927bc194c9200972f136540be..0b18ed706c237cc6ab34ae99d6c0daea14b02c18 100644
--- a/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
+++ b/NumLib/Fem/FiniteElement/C0IsoparametricElements.h
@@ -15,12 +15,15 @@
#define C0ISOPARAMETRICELEMENTS_H_
#include "MeshLib/Elements/Line.h"
+#include "MeshLib/Elements/Tri.h"
#include "MeshLib/Elements/Quad.h"
#include "MeshLib/Elements/Hex.h"
#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
#include "NumLib/Fem/Integration/IntegrationGaussRegular.h"
+#include "NumLib/Fem/Integration/IntegrationGaussTri.h"
#include "TemplateIsoparametric.h"
@@ -33,6 +36,12 @@ struct FeLINE2
typedef TemplateIsoparametric<MeshLib::Line, ShapeLine2, IntegrationGaussRegular<1>, T_SHAPE_VECTOR, T_DSHAPE_MATRIX, T_JACOBIAN_MATRIX> type;
};
+template <class T_SHAPE_VECTOR, class T_DSHAPE_MATRIX, class T_JACOBIAN_MATRIX>
+struct FeTRI3
+{
+ typedef TemplateIsoparametric<MeshLib::Tri, ShapeTri3, IntegrationGaussTri, T_SHAPE_VECTOR, T_DSHAPE_MATRIX, T_JACOBIAN_MATRIX> type;
+};
+
template <class T_SHAPE_VECTOR, class T_DSHAPE_MATRIX, class T_JACOBIAN_MATRIX>
struct FeQUAD4
{
diff --git a/NumLib/Fem/Integration/IntegrationGaussTri.h b/NumLib/Fem/Integration/IntegrationGaussTri.h
new file mode 100644
index 0000000000000000000000000000000000000000..52d3f9d6334926cecb927007a8c463a6a22dca8e
--- /dev/null
+++ b/NumLib/Fem/Integration/IntegrationGaussTri.h
@@ -0,0 +1,129 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 INTEGRATIONGAUSSTRI_H_
+#define INTEGRATIONGAUSSTRI_H_
+
+#include "MathLib/Integration/GaussLegendreTri.h"
+
+namespace NumLib
+{
+
+/**
+ * \brief Gauss quadrature rule for triangles
+ *
+ * Gauss quadrature rule for triangles is originally given as
+ * \f[
+ * \int F(x,y) dx dy = \int F(x(r, s), y(r, s)) j(r,s) dr ds \approx \frac{1}{2} \sum_i ( F(x(r_i, s_i), y(r_i, s_i)) w_i )
+ * \f]
+ *
+ * To make it consistent with other elements, we rewrite the above formula as
+ * \f[
+ * \int F(x,y) dx dy \approx \sum_i ( F(x(r_i, s_i), y(r_i, s_i)) w'_i )
+ * \f]
+ * by defining the new weight \f$ w'=\frac{1}{2} w \f$.
+ */
+class IntegrationGaussTri
+{
+ typedef typename MathLib::TemplateWeightedPoint<double, double, 2>
+ WeightedPoint;
+public:
+ /**
+ * Construct this object with the given integration order
+ *
+ * @param order integration order (default 2)
+ */
+ explicit IntegrationGaussTri(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::GaussLegendreTri<1> >(igp);
+ case 2: return getWeightedPoint<MathLib::GaussLegendreTri<2> >(igp);
+ case 3: return getWeightedPoint<MathLib::GaussLegendreTri<3> >(igp);
+ }
+ return WeightedPoint(std::array<double, 2>(), 0);
+ }
+
+ template <typename Method>
+ static WeightedPoint
+ getWeightedPoint(std::size_t igp)
+ {
+ return WeightedPoint(Method::X[igp], 0.5 * 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::GaussLegendreTri<1>::NPoints;
+ case 2: return MathLib::GaussLegendreTri<2>::NPoints;
+ case 3: return MathLib::GaussLegendreTri<3>::NPoints;
+ }
+ return 0;
+ }
+
+private:
+ std::size_t _order;
+ std::size_t _n_sampl_pt;
+};
+
+}
+
+#endif //INTEGRATIONGAUSSTRI_H_
diff --git a/NumLib/Fem/ShapeFunction/ShapeTri3-impl.h b/NumLib/Fem/ShapeFunction/ShapeTri3-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ddda9a86a78357911b0570477ee23736e080c4c
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeTri3-impl.h
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+namespace NumLib
+{
+
+template <class T_X, class T_N>
+void ShapeTri3::computeShapeFunction(const T_X &r, T_N &N)
+{
+ N[0] = 1. - r[0] - r[1];
+ N[1] = r[0];
+ N[2] = r[1];
+}
+
+template <class T_X, class T_N>
+void ShapeTri3::computeGradShapeFunction(const T_X &/*r*/, T_N &dN)
+{
+ //dN/dr
+ dN[0] = -1.0;
+ dN[1] = 1.0;
+ dN[2] = 0.0;
+ //dN/ds
+ dN[3] = -1.0;
+ dN[4] = 0.0;
+ dN[5] = 1.0;
+}
+
+}
+
diff --git a/NumLib/Fem/ShapeFunction/ShapeTri3.h b/NumLib/Fem/ShapeFunction/ShapeTri3.h
new file mode 100644
index 0000000000000000000000000000000000000000..523d039b63d3fab31e707bcbded89989b79e3994
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeTri3.h
@@ -0,0 +1,58 @@
+/**
+ * \file
+ * \author Norihiro Watanabe
+ * \date 2013-08-13
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2013, 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 SHAPETRI3_H_
+#define SHAPETRI3_H_
+
+namespace NumLib
+{
+
+/**
+ * Shape function for a triangle element of three nodes in natural coordinates
+ *
+ * 3 (0,1)
+ * *
+ * s | \
+ * | \
+ * | \
+ * *––––––* r
+ * 1 (0,0) 2 (1,0)
+ */
+class ShapeTri3
+{
+public:
+ /**
+ * Evaluate the shape function at the given point
+ *
+ * @param [in] r point coordinates
+ * @param [out] N a vector of calculated shape function.
+ */
+ template <class T_X, class T_N>
+ static void computeShapeFunction(const T_X &r, T_N &N3);
+
+ /**
+ * Evaluate derivatives of the shape function at the given point
+ *
+ * @param [in] r point coordinates
+ * @param [out] dN a matrix of the derivatives
+ */
+ template <class T_X, class T_N>
+ static void computeGradShapeFunction(const T_X &/*r*/, T_N &dN6);
+};
+
+}
+
+#include "ShapeTri3-impl.h"
+
+#endif //SHAPETRI3_H_
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestTri3.h b/Tests/NumLib/CoordinatesMappingTestData/TestTri3.h
new file mode 100644
index 0000000000000000000000000000000000000000..6223fceb42e3a2e0968e8bc718d200b9da4c5928
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestTri3.h
@@ -0,0 +1,99 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTTRI3_H_
+#define COORDINATESMAPPINGTESTDATA_TESTTRI3_H_
+
+#include "MeshLib/Elements/Tri.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestTri3
+{
+ public:
+ // Element information
+ typedef MeshLib::Tri ElementType;
+ typedef NumLib::ShapeTri3 ShapeFunctionType;
+ static const unsigned dim = 2; //ElementType::dimension;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+ // Coordinates where shape functions are evaluated
+ static const double r[dim];
+ // Expected results for natural shape
+ static const double nat_exp_N[e_nnodes];
+ static const double nat_exp_dNdr[e_nnodes*dim];
+ // Expected results for irregular shape
+ static const double ir_exp_J[dim*dim];
+ static const double ir_exp_invJ[dim*dim];
+ static const double ir_exp_detJ;
+ static const double ir_exp_dNdx[e_nnodes*dim];
+ // Expected results for clock-wise node ordering
+ static const double cl_exp_J[dim*dim];
+ // Expected results for zero volume
+ static const double cl_exp_detJ;
+ static const double ze_exp_J[dim*dim];
+
+ // element having shape identical to that in natural coordinates
+ MeshLib::Tri* createNaturalShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node( 0.0, 1.0, 0.0);
+ return new MeshLib::Tri(nodes);
+ }
+
+ // element having irregular or skew shape
+ MeshLib::Tri* createIrregularShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 0.1, 0.1, 0.0);
+ nodes[1] = new MeshLib::Node( 2.0, -0.6, 0.0);
+ nodes[2] = new MeshLib::Node( 1.5, 0.4, 0.0);
+ return new MeshLib::Tri(nodes);
+ }
+
+ // invalid case: clock wise node ordering
+ MeshLib::Tri* createClockWise()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 0.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 0.0, 1.0, 0.0);
+ return new MeshLib::Tri(nodes);
+ }
+
+ // invalid case: zero area
+ MeshLib::Tri* createZeroVolume()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ return new MeshLib::Tri(nodes);
+ }
+};
+
+const double TestTri3::r[dim] = {0.5, 0.5};
+const double TestTri3::nat_exp_N[e_nnodes] = {0.0, 0.5, 0.5};
+const double TestTri3::nat_exp_dNdr[e_nnodes*dim] = {-1, 1, 0, -1, 0, 1};
+const double TestTri3::ir_exp_J[dim*dim] = {1.9, -0.7, 1.4, 0.3};
+const double TestTri3::ir_exp_detJ = 1.55;
+const double TestTri3::ir_exp_invJ[dim*dim]
+ = {0.19354838709677424, 0.45161290322580649, -0.90322580645161299, 1.2258064516129032};
+const double TestTri3::ir_exp_dNdx[dim*e_nnodes]
+ = {-0.6451612903225807, 0.1935483870967742, 0.4516129032258065,
+ -0.3225806451612903, -0.903225806451613, 1.225806451612903};
+const double TestTri3::cl_exp_J[dim*dim] = {0.0, 1.0, 1.0, 0.0};
+const double TestTri3::cl_exp_detJ = -1.;
+const double TestTri3::ze_exp_J[dim*dim] = {1.0, 0.0, 1.0, 0.0};
+
+}
+
+#endif
diff --git a/Tests/NumLib/FeTestData/TestFeTRI3.h b/Tests/NumLib/FeTestData/TestFeTRI3.h
new file mode 100644
index 0000000000000000000000000000000000000000..077dd815adde7358e28379276b01190593c969cf
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeTRI3.h
@@ -0,0 +1,77 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, 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 TESTFETRI3_H_
+#define TESTFETRI3_H_
+
+#include "MeshLib/Elements/Tri.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for TRI3
+class TestFeTRI3
+{
+public:
+ // Fe type information
+ template <class T_MATRIX_TYPES>
+ struct FeType
+ {
+ typedef NumLib::FeTRI3<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+ };
+ typedef MeshLib::Tri MeshElementType;
+ static const unsigned dim = 2; //MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+ static const unsigned n_sample_pt_order2 = 3;
+ static const unsigned n_sample_pt_order3 = 4;
+
+ /// create a mesh element
+ MeshElementType* createMeshElement()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node( 0.0, 1.0, 0.0);
+ return new MeshLib::Tri(nodes);
+ }
+
+ /// set an expected mass matrix for 1m x 1m
+ template <class T_MATRIX>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./2;
+ m(1,1) = 1.0; m(1,2) = 1./2;
+ m(2,2) = 1.0;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= 1./6.*0.5;
+ }
+
+ /// set an expected laplace matrix for 1m x 1m
+ template <class T_MATRIX>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 2.0; m(0,1) = -1.0; m(0,2) = -1.0;
+ m(1,1) = 1.0; m(1,2) = 0.0;
+ m(2,2) = 1.0;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= k*0.5;
+ }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/TestCoordinatesMapping.cpp b/Tests/NumLib/TestCoordinatesMapping.cpp
index f4a5bf67f8985269b7a3478f84992a79418810e8..a3766ae222fb0a7d4216aa273207e913af12be5a 100644
--- a/Tests/NumLib/TestCoordinatesMapping.cpp
+++ b/Tests/NumLib/TestCoordinatesMapping.cpp
@@ -20,6 +20,7 @@
#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
#include "CoordinatesMappingTestData/TestLine2.h"
+#include "CoordinatesMappingTestData/TestTri3.h"
#include "CoordinatesMappingTestData/TestQuad4.h"
#include "CoordinatesMappingTestData/TestHex8.h"
@@ -35,6 +36,7 @@ namespace
// Element types to be tested
typedef ::testing::Types<
TestLine2,
+ TestTri3,
TestQuad4,
TestHex8
> TestTypes;
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 13b9fd48cf976f0ac91d7eb5d61b951184f86960..76c569bb06ead9a2f9dcc57b06d4696031219da4 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -19,6 +19,7 @@
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
#include "FeTestData/TestFeLINE2.h"
+#include "FeTestData/TestFeTRI3.h"
#include "FeTestData/TestFeQUAD4.h"
#include "FeTestData/TestFeHEX8.h"
@@ -60,10 +61,12 @@ struct TestCase
typedef ::testing::Types<
TestCase<TestFeLINE2>,
+ TestCase<TestFeTRI3>,
TestCase<TestFeQUAD4>,
TestCase<TestFeHEX8>
#ifdef OGS_USE_EIGEN
,TestCase<TestFeLINE2, EigenFixedMatrixTypes<TestFeLINE2> >
+ ,TestCase<TestFeTRI3, EigenFixedMatrixTypes<TestFeTRI3> >
,TestCase<TestFeQUAD4, EigenFixedMatrixTypes<TestFeQUAD4> >
,TestCase<TestFeHEX8, EigenFixedMatrixTypes<TestFeHEX8> >
#endif