Added integration rules to triangleQuadrature

src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C

@@ -24,27 +24,204 @@ License
 namespace Foam
 {
 
-const List<Tuple2<label, Tuple2<vectorList, scalarList>>> TRI_QUADRATURE
-(
-);
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+Map<triangleQuadrature::quadratureRule> triangleQuadrature::rules_;
+
+// Initialize and return quadrature rules
+void triangleQuadrature::constructRules()
+{
+    if (!rules_.empty())
+    {
+        FatalErrorInFunction
+            << "attempt to re-construct rules when they already exist"
+            << exit(FatalError);
+    }
+
+    // 1 point quadrature, exact for polynomials up to 1 order
+    rules_.insert(1, quadratureRule{
+        List<point>{vector(1.0/3.0, 1.0/3.0, 1.0/3.0)},
+        List<scalar>{1.0}
+    });
+
+    // 3 point quadrature, exact for polynomials up to 2 order
+    rules_.insert(3, quadratureRule{
+        List<point>
+        {
+            point(2.0/3.0, 1.0/6.0, 1.0/6.0),
+            point(1.0/6.0, 2.0/3.0, 1.0/6.0),
+            point(1.0/6.0, 1.0/6.0, 2.0/3.0)
+        },
+        List<scalar>{1.0/3.0, 1.0/3.0, 1.0/3.0}
+    });
+
+    // 4 point quadrature, exact for polynomials up to 3 order
+    rules_.insert(4, quadratureRule{
+        List<point>
+        {
+            point(1.0/3.0, 1.0/3.0, 1.0/3.0),
+            point(0.6, 0.2, 0.2),
+            point(0.2, 0.6, 0.2),
+            point(0.2, 0.2, 0.6)
+        },
+        List<scalar>{-9.0/16.0, 25.0/48.0, 25.0/48.0, 25.0/48.0}
+    });
+
+    // 6 point quadrature, exact for polynomials up to 4 order
+    rules_.insert(6, quadratureRule{
+        List<point>
+        {
+            point(0.108103018168070, 0.445948490915965, 0.445948490915965),
+            point(0.445948490915965, 0.108103018168070, 0.445948490915965),
+            point(0.445948490915965, 0.445948490915965, 0.108103018168070),
+            point(0.816847572980459, 0.091576213509771, 0.091576213509771),
+            point(0.091576213509771, 0.816847572980459, 0.091576213509771),
+            point(0.091576213509771, 0.091576213509771, 0.816847572980459)
+        },
+        List<scalar>
+        {
+            0.223381589678011,
+            0.223381589678011,
+            0.223381589678011,
+            0.109951743655322,
+            0.109951743655322,
+            0.109951743655322
+        }
+    });
+
+    // 7 point quadrature, exact for polynomials up to 5 order
+    rules_.insert(7, quadratureRule{
+        List<point>
+        {
+            point(1.0/3.0, 1.0/3.0, 1.0/3.0),
+            point(0.059715871789770, 0.470142064105115, 0.470142064105115),
+            point(0.470142064105115, 0.059715871789770, 0.470142064105115),
+            point(0.470142064105115, 0.470142064105115, 0.059715871789770),
+            point(0.797426985353087, 0.101286507323456, 0.101286507323456),
+            point(0.101286507323456, 0.797426985353087, 0.101286507323456),
+            point(0.101286507323456, 0.101286507323456, 0.797426985353087),
+        },
+        List<scalar>
+        {
+             0.225000000000000,
+             0.132394152788506,
+             0.132394152788506,
+             0.132394152788506,
+             0.125939180544827,
+             0.125939180544827,
+             0.125939180544827
+        }
+    });
+
+    // 12 point quadrature, exact for polynomials up to 6 order
+    rules_.insert(12, quadratureRule{
+        List<point>
+        {
+            point(0.501426509658179, 0.249286745170910, 0.249286745170910),
+            point(0.249286745170910, 0.501426509658179, 0.249286745170910),
+            point(0.249286745170910, 0.249286745170910, 0.501426509658179),
+            point(0.873821971016996, 0.063089014491502, 0.063089014491502),
+            point(0.063089014491502, 0.873821971016996, 0.063089014491502),
+            point(0.063089014491502, 0.063089014491502, 0.873821971016996),
+            point(0.053145049844817, 0.310352451033784, 0.636502499121399),
+            point(0.053145049844817, 0.636502499121399, 0.310352451033784),
+            point(0.310352451033784, 0.053145049844817, 0.636502499121399),
+            point(0.636502499121399, 0.053145049844817, 0.310352451033784),
+            point(0.310352451033784, 0.636502499121399, 0.053145049844817),
+            point(0.636502499121399, 0.310352451033784, 0.053145049844817)
+        },
+        List<scalar>
+        {
+            0.116786275726379,
+            0.116786275726379,
+            0.116786275726379,
+            0.050844906370207,
+            0.050844906370207,
+            0.050844906370207,
+            0.082851075618374,
+            0.082851075618374,
+            0.082851075618374,
+            0.082851075618374,
+            0.082851075618374,
+            0.082851075618374
+        }
+    });
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+List<point> triangleQuadrature::mapGaussPoints(const List<vector>& pts)
+{
+    List<point> gp;
+    gp.setSize(pts.size());
+
+    NotImplemented;
+
+    return gp;
+}
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+const Map<triangleQuadrature::quadratureRule>& triangleQuadrature::rules()
+{
+   if (rules_.empty())
+   {
+       constructRules();
+   }
+
+   return rules_;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 
 triangleQuadrature::triangleQuadrature(const triPoints& pts, const label& n)
 :
     triPoints(pts),
-    n_(n),
-    gaussPointsWeights_(n)
+    n_(n)
 {
-    mapGaussPoints();
+    constructRules();
+
+    // Check if the requested number of points exists in the rules
+    if (!rules().found(n_))
+    {
+         FatalErrorInFunction
+            << "Gaussian quadrature rule for " << n_ << " points not available!"
+            << abort(FatalError);
+    }
+
+    const quadratureRule& r = rules()[n_];
+    weights_ = r.weights;
+    gaussPoints_ = mapGaussPoints(r.gaussPoints);
 }
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void triangleQuadrature::mapGaussPoints()
+const List<point>& triangleQuadrature::gaussPoints() const
 {
+    return gaussPoints_;
+}
+
 
+const List<scalar>& triangleQuadrature::weights() const
+{
+    return weights_;
 }
 
+
+const List<point> triangleQuadrature::gaussPoints()
+{
+    return gaussPoints_;
 }
 
+
+const List<scalar> triangleQuadrature::weights()
+{
+     return weights_;
+}
+
+
+} // End namespace Foam
+
 // ************************************************************************* //

src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.H

@@ -19,7 +19,7 @@ Function
     triangleQuadrature
 
 Description
-    Gaussian quadrature integration for triangle
+    Gauss quadrature for triangle
 
     Source:
     D. A. Dunavant,
@@ -43,7 +43,7 @@ Author
 #include "Tuple2.H"
 #include "vectorList.H"
 #include "scalarList.H"
-
+#include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,16 +56,32 @@ class triangleQuadrature
 {
     // Private data
 
+        //- Quadrature rules storage
+        struct quadratureRule
+        {
+            List<vector> gaussPoints;
+            List<scalar> weights;
+        };
+
+        //- Static storage for quadrature rules
+        static Map<quadratureRule> rules_;
+
+        //- Construct rules
+        static void constructRules();
+
         //- Number of Gauss points
         const label n_;
 
-        //- Mapped Gauss points and corresponding weights
-        List<Tuple2<point, scalar>> gaussPointsWeights_;
+        //- Gaussian quadrature points in Cartesian coordinates
+        List<point> gaussPoints_;
+
+        //- Gaussian quadrature points weights
+        List<scalar> weights_;
 
     // Private Member Functions
 
         //- Map integration points from local natural to cartesian coordinates
-        void mapGaussPoints();
+        List<point> mapGaussPoints(const List<point>& pts);
 
 public:
 
@@ -74,15 +90,23 @@ public:
 
     // Member Functions
 
-    // Properties
+    // Access
+
+        //- Predefined integration rules
+        static const Map<quadratureRule>& rules();
+
+        //- Get reference to Gaussian quadrature points in global Cartesian
+        //  coordinates
+        const List<point>& gaussPoints() const;
+
+        //- Get reference to Gaussian quadrature points weights
+        const List<scalar>& weights() const;
 
-        //- Triangle quadrature points and weights in local natural coordinates
-        static const List<Tuple2
-            <
-                label,
-                Tuple2<vectorList, scalarList>
-            >> TRI_QUADRATURE;
+        //- Get Gaussian quadrature points in global Cartesian coordinates
+        const List<point> gaussPoints();
 
+        //- Get Gaussian quadrature points weights
+        const List<scalar> weights();
 };
 
 } // End namespace Foam