From 8012962a9268aa06f1ffb67efeec14ebf13ccd36 Mon Sep 17 00:00:00 2001 From: Ivan Batistic Date: Thu, 28 Nov 2024 11:28:27 +0100 Subject: [PATCH] Added integration rules to triangleQuadrature --- .../triangleQuadrature/triangleQuadrature.C | 191 +++++++++++++++++- .../triangleQuadrature/triangleQuadrature.H | 48 +++-- 2 files changed, 220 insertions(+), 19 deletions(-) diff --git a/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C b/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C index 7a8a157cc..3d5fbb7d6 100644 --- a/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C +++ b/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C @@ -24,27 +24,204 @@ License namespace Foam { -const List>> TRI_QUADRATURE -( -); +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // +Map 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{vector(1.0/3.0, 1.0/3.0, 1.0/3.0)}, + List{1.0} + }); + + // 3 point quadrature, exact for polynomials up to 2 order + rules_.insert(3, quadratureRule{ + List + { + 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{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(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{-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(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 + { + 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(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 + { + 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(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 + { + 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 triangleQuadrature::mapGaussPoints(const List& pts) +{ + List gp; + gp.setSize(pts.size()); + + NotImplemented; + + return gp; +} + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +const Map& 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& triangleQuadrature::gaussPoints() const { + return gaussPoints_; +} + +const List& triangleQuadrature::weights() const +{ + return weights_; } + +const List triangleQuadrature::gaussPoints() +{ + return gaussPoints_; } + +const List triangleQuadrature::weights() +{ + return weights_; +} + + +} // End namespace Foam + // ************************************************************************* // diff --git a/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.H b/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.H index e03fbd15b..fdc57a758 100644 --- a/src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.H +++ b/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 gaussPoints; + List weights; + }; + + //- Static storage for quadrature rules + static Map rules_; + + //- Construct rules + static void constructRules(); + //- Number of Gauss points const label n_; - //- Mapped Gauss points and corresponding weights - List> gaussPointsWeights_; + //- Gaussian quadrature points in Cartesian coordinates + List gaussPoints_; + + //- Gaussian quadrature points weights + List weights_; // Private Member Functions //- Map integration points from local natural to cartesian coordinates - void mapGaussPoints(); + List mapGaussPoints(const List& pts); public: @@ -74,15 +90,23 @@ public: // Member Functions - // Properties + // Access + + //- Predefined integration rules + static const Map& rules(); + + //- Get reference to Gaussian quadrature points in global Cartesian + // coordinates + const List& gaussPoints() const; + + //- Get reference to Gaussian quadrature points weights + const List& weights() const; - //- Triangle quadrature points and weights in local natural coordinates - static const List - >> TRI_QUADRATURE; + //- Get Gaussian quadrature points in global Cartesian coordinates + const List gaussPoints(); + //- Get Gaussian quadrature points weights + const List weights(); }; } // End namespace Foam