fGradGaussPoints func development

src/solids4FoamModels/Make/files.openfoam

@@ -236,5 +236,6 @@ numerics/pointPointLeastSquaresVectors/pointPointLeastSquaresVectors.C
 numerics/sparseMatrix/sparseMatrix.C
 numerics/sparseMatrix/SparseMatrixTemplates.C
 numerics/sparseMatrix/sparseMatrixTools.C
+numerics/triangleQuadrature/triangleQuadrature.C
 
 LIB = $(FOAM_MODULE_LIBBIN)/libsolids4FoamModels

src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.C

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+License
+    This file is part of solids4foam.
+
+    solids4foam is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    solids4foam is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with solids4foam.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triangleQuadrature.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+const List<Tuple2<label, Tuple2<vectorList, scalarList>>> TRI_QUADRATURE
+(
+);
+
+
+
+triangleQuadrature::triangleQuadrature(const triPoints& pts, const label& n)
+:
+    triPoints(pts),
+    n_(n),
+    gaussPointsWeights_(n)
+{
+    mapGaussPoints();
+}
+
+
+void triangleQuadrature::mapGaussPoints()
+{
+
+}
+
+}
+
+// ************************************************************************* //

src/solids4FoamModels/numerics/triangleQuadrature/triangleQuadrature.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+License
+    This file is part of solids4foam.
+
+    solids4foam is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    solids4foam is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with solids4foam.  If not, see <http://www.gnu.org/licenses/>.
+
+Function
+    triangleQuadrature
+
+Description
+    Gaussian quadrature integration for triangle
+
+    Source:
+    D. A. Dunavant,
+    High degree efficient symmetrical Gaussian quadrature rule for the triangle.
+    International Journal for Numerical Methods in Engineering,
+    21(6):1129–1148, jun 1985.
+
+SourceFile
+    triangleQuadrature.C
+
+Author
+    Ivan Batistic, FAMENA. All rights reserved.
+    Philip Cardiff, UCD. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triangleQuadrature_H
+#define triangleQuadrature_H
+
+#include "triPoints.H"
+#include "Tuple2.H"
+#include "vectorList.H"
+#include "scalarList.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class triangleQuadrature
+:
+    public triPoints
+{
+    // Private data
+
+        //- Number of Gauss points
+        const label n_;
+
+        //- Mapped Gauss points and corresponding weights
+        List<Tuple2<point, scalar>> gaussPointsWeights_;
+
+    // Private Member Functions
+
+        //- Map integration points from local natural to cartesian coordinates
+        void mapGaussPoints();
+
+public:
+
+    // Constructors
+    triangleQuadrature(const triPoints& pts, const label& n);
+
+    // Member Functions
+
+    // Properties
+
+        //- Triangle quadrature points and weights in local natural coordinates
+        static const List<Tuple2
+            <
+                label,
+                Tuple2<vectorList, scalarList>
+            >> TRI_QUADRATURE;
+
+};
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/solids4FoamModels/solidModels/linGeomTotalDispSolid/higherOrderGrad/higherOrderGrad.C

@@ -22,7 +22,8 @@ License
 #include "surfaceFields.H"
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
-
+#include "triangle.H"
+#include "triFace.H"
 
 namespace Foam
 {
@@ -2228,6 +2229,8 @@ higherOrderGrad::~higherOrderGrad()
 
 tmp<volTensorField> higherOrderGrad::grad(const volVectorField& D) const
 {
+    auto test =  fGradGaussPoints(D);
+
     if (useQRDecomposition_)
     {
         if (useGlobalStencils_)
@@ -2358,6 +2361,111 @@ tmp<surfaceTensorField> higherOrderGrad::fGrad(const volVectorField& D) const
 }
 
 
+refPtr<List<List<tensor>>> higherOrderGrad::fGradGaussPoints(const volVectorField& D) const
+{
+    // Number of quadrature points per triangle.
+    const label triQuadraturePtsNb = 3;
+
+    const fvMesh& mesh = mesh_;
+
+    const pointField& pts = mesh.points();
+
+    // Prepare the return field
+    // Philip: I'm using refPtr, tmp is more appropriate but I had a compilation
+    // problems, when it is used for List<List<tensor>
+    refPtr<List<List<tensor>>> tgradDGP(new List<List<tensor>>(mesh.nFaces()));
+    List<List<tensor>>& gradDGP = tgradDGP.ref();
+
+    forAll(gradDGP, i)
+    {
+        List<tensor>& faceGradGP = gradDGP[i];
+
+        const label nbOfTriangles = mesh.faces()[i].size();
+        const label nbOfGaussPoints = nbOfTriangles * triQuadraturePtsNb;
+
+        faceGradGP.setSize(nbOfGaussPoints);
+    }
+
+    // 1. Stage - decompose faces into triangles. Store triangle points
+
+    // Triangulate each face and store points of each triangle
+    List<List<triPoints>> faceTri(mesh.nFaces());
+    forAll(faceTri, i)
+    {
+        List<triPoints>& fT = faceTri[i];
+
+        fT.setSize(mesh.faces()[i].size());
+    }
+
+    // Loop over faces and decompose each face, store triangles of each face
+    for (label faceI = 0; faceI < mesh.nFaces(); ++faceI)
+    {
+        const face& f = mesh.faces()[faceI];
+
+        const point fc = f.centre(pts);
+
+        const label nPoints = f.size();
+
+        label nextpI;
+        for (label pI = 0; pI<nPoints; ++pI)
+        {
+            if (pI < f.size() - 1)
+            {
+                nextpI = pI + 1;
+            }
+            else
+            {
+                nextpI = 0;
+            }
+
+            const triPoints tri
+            (
+                pts[f[pI]],
+                pts[f[nextpI]],
+                fc
+            );
+
+            faceTri[faceI][pI] = tri;
+        }
+    }
+
+    // 2. Stage - for each triangle calculate Gauss point locations and store
+    //            corresponding weights
+
+    // Gauss point locations on each face
+    List<List<point>> faceGP(mesh.nFaces());
+    forAll(faceGP, i)
+    {
+        List<point>& fGP = faceGP[i];
+        fGP.setSize(mesh.faces()[i].size() * triQuadraturePtsNb);
+    }
+
+    // Gauss point weights
+    List<List<scalar>> faceGPW(mesh.nFaces());
+    forAll(faceGPW, i)
+    {
+        List<scalar>& fGPW = faceGPW[i];
+        fGPW.setSize(mesh.faces()[i].size() * triQuadraturePtsNb);
+    }
+
+    // Loop over faces; loop over triangles. Store Gauss point locations and
+    // weight
+    forAll(faceTri, faceI)
+    {
+        List<triPoints>& fT = faceTri[faceI];
+
+        forAll(fT, tI)
+        {
+            const triPoints& tp = fT[tI];
+
+            // Sad bi tu isla funkcija koja mi vraca Gaussove tocke i tezine
+        }
+    }
+
+    return tgradDGP;
+}
+
+
 tmp<volTensorField> higherOrderGrad::gradQR(const volVectorField& D) const
 {
     if (debug)

src/solids4FoamModels/solidModels/linGeomTotalDispSolid/higherOrderGrad/higherOrderGrad.H

@@ -41,6 +41,7 @@ SourceFiles
 #include "fvMesh.H"
 #include "globalIndex.H"
 #include <Eigen/Dense>
+#include "tensorList.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -261,6 +262,9 @@ public:
             //- Calculate the face gradient of a volVectorField
             tmp<surfaceTensorField> fGrad(const volVectorField& D) const;
 
+            //- Calculate the face Gauss points gradient of a volVectorField
+            refPtr<List<List<tensor>>> fGradGaussPoints(const volVectorField& D) const;
+
             //- Calculate the gradient of a volVectorField using QR
             //  decomposition
             tmp<volTensorField> gradQR(const volVectorField& D) const;