Added correct function with Gauss point stress and grad

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalLaws/linearGeometryLaws/linearElastic/linearElastic.C

@@ -398,5 +398,25 @@ void Foam::linearElastic::correct
     }
 }
 
+// TESTING
+void Foam::linearElastic::correct
+(
+    const List<List<tensor>>& gradDGPf,
+    List<List<symmTensor>>& sigmaGPf
+)
+{
+    forAll(sigmaGPf, faceI)
+    {
+        const List<symmTensor>& faceSigmaGP = sigmaGPf[faceI];
+
+        forAll(faceSigmaGP, gpI)
+        {
+            const symmTensor epsilon = symm(gradDGPf[faceI][gpI]);
+
+            sigmaGPf[faceI][gpI] =
+                (2.0*mu_*epsilon + lambda_*tr(epsilon)*symmTensor::I).value();
+        }
+    }
+}
 
 // ************************************************************************* //

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalLaws/linearGeometryLaws/linearElastic/linearElastic.H

@@ -153,6 +153,14 @@ public:
         (
             pointSymmTensorField& sigma, const pointTensorField& gradD
         );
+
+        // TESTING
+        //- Calculate the stress at Gauss points given the displacemen gradient
+        virtual void correct
+        (
+            const List<List<tensor>>& gradDGPf,
+            List<List<symmTensor>>& sigmaGPf
+        );
 };
 
 

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalLaws/mechanicalLaw/mechanicalLaw.C

@@ -1648,6 +1648,22 @@ void Foam::mechanicalLaw::correct
 }
 
 
+void Foam::mechanicalLaw::correct
+(
+    const List<List<tensor>>& gradDGPf,
+    List<List<symmTensor>>& sigmaGPf
+)
+{
+    notImplemented
+    (
+        type() + "::correct( const List<List<tensor>>& gradDGPf, "
+        "List<List<symmTensor>>& sigmaGPf)\n"
+        "The correct function is not implemented\n"
+        " for the " + type() + " mechanical law"
+    );
+}
+
+
 Foam::scalar Foam::mechanicalLaw::residual()
 {
     // Default to zero; this can be overwritten by any derived mechanical law

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalLaws/mechanicalLaw/mechanicalLaw.H

@@ -546,6 +546,14 @@ public:
             pointSymmTensorField& sigma, const pointTensorField& gradD
         );
 
+        // TESTING
+        // Added for high order stress calculation at Gauss points
+        virtual void correct
+        (
+            const List<List<tensor>>& gradDGPf,
+            List<List<symmTensor>>& sigmaGPf
+        );
+
         //- Return material residual i.e. a measured of how convergence of
         //  the material model
         virtual scalar residual();

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalModel.C

@@ -553,6 +553,29 @@ void Foam::mechanicalModel::correct
 }
 
 
+void Foam::mechanicalModel::correct
+(
+   const List<List<tensor>>& gradDGPf,
+   List<List<symmTensor>>& sigmaGPf
+)
+{
+    PtrList<mechanicalLaw>& laws = *this;
+
+    if (laws.size() == 1)
+    {
+        laws[0].correct(gradDGPf, sigmaGPf);
+    }
+    else
+    {
+        notImplemented
+        (
+            "mechanicalModel::correct(...): not implemented for more than "
+            "one material"
+        );
+    }
+}
+
+
 void Foam::mechanicalModel::mapGradToSubMeshes(const volTensorField& gradD)
 {
     const PtrList<mechanicalLaw>& laws = *this;

src/solids4FoamModels/materialModels/mechanicalModel/mechanicalModel.H

@@ -171,6 +171,12 @@ public:
                 pointSymmTensorField& sigma, const pointTensorField& gradD
             );
 
+            void correct
+            (
+                const List<List<tensor>>& gradDGPf,
+                List<List<symmTensor>>& sigmaGPf
+            );
+
             //- Map grad field to subMeshes
             void mapGradToSubMeshes(const volTensorField& gradD);
 

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

@@ -279,6 +279,9 @@ public:
 
     // Member Functions
 
+            //- Return number of Gauss points per triangle
+            label triQuadraturePtsNb() const {return triQuadraturePtsNb_;};
+
         // Edit
 
             //- Calculate the gradient of a volVectorField

src/solids4FoamModels/solidModels/linGeomTotalDispSolid/linGeomTotalDispSolid.C

@@ -476,7 +476,8 @@ linGeomTotalDispSolid::linGeomTotalDispSolid
         mesh(),
         dimensionedScalar("ds", (dimForce/dimVolume)/dimVelocity, 1.0)
     ),
-    hoGradPtr_()
+    hoGradPtr_(),
+    sigmaGPfPtr_()
 {
     DisRequired();
 
@@ -515,6 +516,23 @@ linGeomTotalDispSolid::linGeomTotalDispSolid
                 solidModelDict().subDict("higherOrderGradCoeffs")
             )
         );
+
+        sigmaGPfPtr_.set
+        (
+            new List<List<symmTensor>>(mesh().nFaces())
+        );
+
+        List<List<symmTensor>>& sigmaGPf = sigmaGPfPtr_.ref();
+
+        forAll(sigmaGPf, i)
+        {
+            List<symmTensor>& faceSigmaGP = sigmaGPf[i];
+
+            const label nbOfTriangles = mesh().faces()[i].size() ;
+            const label nbOfGaussPoints = nbOfTriangles*hoGradPtr_->triQuadraturePtsNb();
+
+            faceSigmaGP.setSize(nbOfGaussPoints);
+        }
     }
 
     if (predictor_)
@@ -729,6 +747,27 @@ tmp<vectorField> linGeomTotalDispSolid::residualMomentum
     // Traction vectors at the faces
     surfaceVectorField traction(n & fvc::interpolate(sigma()));
 
+    // TESTING
+    // Here we will calculate displacement gradient at Gauss points and call
+    // mechanical law to calculate sigma from it. To check implementation
+    // we can compare traction at Gauss points with traction at face centre
+    // from standard discretisation.
+    const List<List<tensor>> gradDGPf = hoGradPtr_->fGradGaussPoints(D);
+
+    List<List<symmTensor>>& sigmaGPf = sigmaGPfPtr_.ref();
+
+    mechanical().correct(gradDGPf, sigmaGPf);
+
+    // forAll(sigmaGPf, faceI)
+    // {
+    //     Info<<"Traction at gauss points: " << endl;
+    //     forAll(sigmaGPf[faceI], gpI)
+    //     {
+    //         Info<< (sigmaGPf[faceI][gpI] & n[faceI]) << endl;;
+    //     }
+    //     Info<<"Traction at face centre: " << traction[faceI] << nl << endl;
+    // }
+
     // Add stabilisation to the traction
     // We add this before enforcing the traction condition as the stabilisation
     // is set to zero on traction boundaries

src/solids4FoamModels/solidModels/linGeomTotalDispSolid/linGeomTotalDispSolid.H

@@ -92,6 +92,10 @@ class linGeomTotalDispSolid
         //- Higher order gradient calculator
         autoPtr<higherOrderGrad> hoGradPtr_;
 
+        //- Sigma at Gauss Points, this should live in solidModel. Added
+        //  here just for testing purpose.
+        mutable autoPtr<List<List<symmTensor>>> sigmaGPfPtr_;
+
 
     // Private Member Functions