Support stress data

FSI/FSI.C

@@ -188,6 +188,15 @@ void preciceAdapter::FSI::FluidStructureInteraction::addWriters(std::string data
         );
         DEBUG(adapterInfo("Added writer: Displacement."));
     }
+    else if(dataName.find("Stress") == 0)
+    {
+      interface->addCouplingDataWriter
+          (
+              dataName,
+              new Stress(mesh_, runTime_.timeName(), solverType_) /* TODO: Add any other arguments here */
+              );
+      DEBUG(adapterInfo("Added writer: Stress."));
+    }
 
     // NOTE: If you want to couple another variable, you need
     // to add your new coupling data user as a coupling data
@@ -225,6 +234,16 @@ void preciceAdapter::FSI::FluidStructureInteraction::addReaders(std::string data
         );
         DEBUG(adapterInfo("Added reader: Displacement."));
     }
+    else if(dataName.find("Stress") == 0)
+    {
+      interface->addCouplingDataReader
+          (
+              dataName,
+              new Stress(mesh_, runTime_.timeName(), solverType_) /* TODO: Add any other arguments here */
+              );
+      DEBUG(adapterInfo("Added reader: Stress."));
+    }
+
     // NOTE: If you want to couple another variable, you need
     // to add your new coupling data user as a coupling data
     // writer here (and as a writer above).

FSI/FSI.H

@@ -6,6 +6,7 @@
 #include "FSI/Displacement.H"
 #include "FSI/DisplacementDelta.H"
 #include "FSI/Force.H"
+#include "FSI/Stress.H"
 
 #include "fvCFD.H"
 

FSI/Force.C

@@ -2,6 +2,16 @@
 
 using namespace Foam;
 
+preciceAdapter::FSI::Force::Force
+(
+        const Foam::fvMesh& mesh,
+        const std::string solverType
+)
+:
+mesh_(mesh),solverType_(solverType)
+{}
+
+
 preciceAdapter::FSI::Force::Force
 (
     const Foam::fvMesh& mesh,
@@ -13,8 +23,7 @@ preciceAdapter::FSI::Force::Force
     */
 )
 :
-mesh_(mesh),
-solverType_(solverType)
+Force(mesh, solverType)
 {
     //What about type "basic"?
     if (solverType_.compare("incompressible") != 0 && solverType_.compare("compressible") != 0) 

FSI/Force.H

@@ -5,7 +5,6 @@
 
 #include "fvCFD.H"
 
-// TEMPORARY
 #include "pointFields.H"
 #include "vectorField.H"
 #include "immiscibleIncompressibleTwoPhaseMixture.H"
@@ -31,7 +30,7 @@ private:
     //- Force field
     Foam::volVectorField * Force_;
 
-
+protected:
     //- Stress tensor (see the OpenFOAM "Forces" function object)
     Foam::tmp<Foam::volSymmTensorField> devRhoReff() const;
 
@@ -56,10 +55,17 @@ public:
         */
     );
 
-    //- Write the displacement values into the buffer
+    //- Constructor
+    Force
+    (
+        const Foam::fvMesh& mesh,
+        const std::string solverType
+    );
+
+    //- Write the forces values into the buffer
     void write(double * buffer, bool meshConnectivity, const unsigned int dim);
 
-    //- Read the displacement values from the buffer
+    //- Read the forces values from the buffer
     void read(double * buffer, const unsigned int dim);
 
     //- Destructor

FSI/Stress.C

@@ -0,0 +1,141 @@
+#include "Stress.H"
+
+using namespace Foam;
+
+preciceAdapter::FSI::Stress::Stress
+(
+    const Foam::fvMesh& mesh,
+    const fileName& timeName,
+    const std::string solverType
+    /* TODO: We should add any required field names here.
+    /  They would need to be vector fields.
+    /  See FSI/Temperature.C for details.
+    */
+)
+:
+Force(mesh,solverType),
+mesh_(mesh),
+solverType_(solverType)
+{
+    if (solverType_.compare("incompressible") != 0 && solverType_.compare("compressible") != 0)
+    {
+        FatalErrorInFunction
+            << "Stresses calculation only supports "
+            << "compressible or incompressible solver type."
+            << exit(FatalError);
+    }
+
+    dataType_ = vector;
+
+    Stress_ = new volVectorField
+    (
+        IOobject
+        (
+            "Stress",
+            timeName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector
+        (
+            "pdim",
+            dimensionSet(1,-1,-2,0,0,0,0),
+            Foam::vector::zero
+        )
+    );
+}
+
+void preciceAdapter::FSI::Stress::write(double * buffer, bool meshConnectivity, const unsigned int dim)
+{
+    // Compute stress. See the Forces function object.
+
+    // Stress tensor boundary field
+    tmp<volSymmTensorField> tdevRhoReff = this->devRhoReff();
+    const volSymmTensorField::Boundary& devRhoReffb =
+        tdevRhoReff().boundaryField();
+
+    // Density boundary field
+    tmp<volScalarField> trho = this->rho();
+    const volScalarField::Boundary& rhob =
+        trho().boundaryField();
+
+    // Pressure boundary field
+    tmp<volScalarField> tp = mesh_.lookupObject<volScalarField>("p");
+    const volScalarField::Boundary& pb =
+        tp().boundaryField();        
+
+    int bufferIndex = 0;
+    // For every boundary patch of the interface
+    for (uint j = 0; j < patchIDs_.size(); j++)
+    {
+
+        int patchID = patchIDs_.at(j);
+
+        // Compute normal vectors on each patch
+        const vectorField nV = mesh_.boundary()[patchID].nf();
+
+        // Pressure constribution
+        if (solverType_.compare("incompressible") == 0)
+        {
+            Stress_->boundaryFieldRef()[patchID] =
+                nV * pb[patchID] * rhob[patchID];
+        }
+        else if (solverType_.compare("compressible") == 0)
+        {
+            Stress_->boundaryFieldRef()[patchID] =
+                nV * pb[patchID];
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Stress calculation does only support "
+                << "compressible or incompressible solver type."
+                << exit(FatalError);
+        }
+
+        // Viscous contribution
+        Stress_->boundaryFieldRef()[patchID] +=
+            nV & devRhoReffb[patchID];
+
+        // Write the stress to the preCICE buffer
+        // For every cell of the patch
+        forAll(Stress_->boundaryFieldRef()[patchID], i)
+        {
+            // Copy the stress into the buffer
+            // x-dimension
+            buffer[bufferIndex++]
+            = 
+            Stress_->boundaryFieldRef()[patchID][i].x();
+
+            // y-dimension
+            buffer[bufferIndex++]
+            =
+            Stress_->boundaryFieldRef()[patchID][i].y();
+
+            if(dim == 3)
+                // z-dimension
+                buffer[bufferIndex++]
+                        =
+                        Stress_->boundaryFieldRef()[patchID][i].z();
+        }
+    }
+}
+
+void preciceAdapter::FSI::Stress::read(double * buffer, const unsigned int dim)
+{
+    /* TODO: Implement
+    * We need two nested for-loops for each patch,
+    * the outer for the locations and the inner for the dimensions.
+    * See the preCICE readBlockVectorData() implementation.
+    */
+    FatalErrorInFunction
+        << "Reading stresses is not supported."
+        << exit(FatalError);
+}
+
+preciceAdapter::FSI::Stress::~Stress()
+{
+    delete Stress_;
+}

FSI/Stress.H

@@ -0,0 +1,63 @@
+#ifndef FSI_STRESS_H
+#define FSI_STRESS_H
+
+#include "CouplingDataUser.H"
+#include "Force.H"
+
+#include "fvCFD.H"
+
+#include "pointFields.H"
+#include "vectorField.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulentFluidThermoModel.H"
+#include "turbulentTransportModel.H"
+
+namespace preciceAdapter
+{
+namespace FSI
+{
+
+//- Class that writes and reads stress [N/m^2]:
+// This is essentially a force (in spatial coordinates) scaled by the deformed
+// cell face. Thus, a consistent quantity. Calculation concept has been copied
+// from the force module, but the scaled version here is commonly used in FEM
+// applications.
+class Stress : public Force {
+
+private:
+
+    //- OpenFOAM fvMesh object (we need to access the objects' registry multiple times)
+    const Foam::fvMesh& mesh_;
+
+    const std::string solverType_;
+
+    //- Stress field
+    Foam::volVectorField * Stress_;
+
+
+public:
+
+    //- Constructor
+    Stress
+    (
+        const Foam::fvMesh& mesh,
+        const fileName& timeName,
+        const std::string solverType
+        // We create an IOobject and we need the time directory
+        /* TODO: We should add any required field names here.
+        */
+    );
+
+    //- Write the stress values into the buffer
+    void write(double * buffer, bool meshConnectivity, const unsigned int dim);
+
+    //- Read the stress values from the buffer
+    void read(double * buffer, const unsigned int dim);
+
+    //- Destructor
+    ~Stress();
+};
+}
+}
+
+#endif

Make/files

@@ -14,6 +14,7 @@ CHT/CHT.C
 
 FSI/FSI.C
 FSI/Force.C
+FSI/Stress.C
 FSI/Displacement.C
 FSI/DisplacementDelta.C