feat: Added option to postprocess Deviatoric stress

Author: giavancini

Material/Elasticity/TPZElasticityTH.cpp

@@ -304,6 +304,8 @@ int TPZElasticityTH::VariableIndex(const std::string &name) const
         return EStrain;
     if (!strcmp("VonMises", name.c_str()))
         return EVonMises;
+    if (!strcmp("DeviatoricStress", name.c_str()))
+        return EDeviatoricStress;
 
     std::cout << "\n\nVar index not implemented\n\n";
     DebugStop();
@@ -327,6 +329,7 @@ int TPZElasticityTH::NSolutionVariables(int var) const
         break;
     case EStress: // stress tensor
     case EStrain: // strain tensor
+    case EDeviatoricStress: // deviatoric stress tensor
         aux = 9;
         break;
     default:
@@ -475,6 +478,35 @@ void TPZElasticityTH::Solution(const TPZVec<TPZMaterialDataT<STATE>> &datavec, i
         Solout[0] = sqrt(0.5 * Solout[0]);
         break;
     }
+    case EDeviatoricStress:
+    {
+        TPZFNMatrix<6, STATE> sigmavoight(n, 1, 0.0);
+        DeviatoricStressTensor(gradU, sigmavoight);
+
+        TPZFNMatrix<9, STATE> devsigma(3, 3, 0.0);
+        int cont = fdimension - 1;
+        for (int i = 0; i < fdimension; i++)
+        {
+            devsigma(i, i) = sigmavoight(i, 0);
+            for (int j = i + 1; j < fdimension; j++)
+            {
+                devsigma(i, j) = sigmavoight(++cont, 0);
+                devsigma(j, i) = sigmavoight(cont, 0);
+            }
+        }
+
+        if (fAnalysisType == AnalysisType::EPlaneStrain)
+            devsigma(2, 2) = fthickness * fpoisson * (devsigma(0, 0) + devsigma(1, 1));
+
+        for (int i = 0; i < 3; i++)
+        {
+            for (int j = 0; j < 3; j++)
+            {
+                Solout[i * 3 + j] = devsigma(i, j);
+            }
+        }
+        break;
+    }
 
     default:
     {

Material/Elasticity/TPZElasticityTH.h

@@ -130,7 +130,9 @@ class TPZElasticityTH : public TPZMatBase<STATE, TPZMatCombinedSpacesT<STATE>, T
 
     virtual void StressTensor(const TPZFNMatrix<10, STATE>& gradU, TPZFNMatrix<6,REAL>& sigma, REAL pressure); //This function computes sigma through its deviatoric and hidrostatic counterparts
 
-    enum SolutionVars {ENone = -1, EPressure = 0, EDisplacement = 1, EForce = 2, EStress = 3, EStrain = 4, EVonMises = 5};
+    enum SolutionVars {ENone = -1, EPressure = 0, EDisplacement = 1, EForce = 2, EStress = 3, EStrain = 4, EVonMises = 5, EDeviatoricStress = 6};
+
+    
 };
 
 #endif

Material/Elasticity/TPZHybridMixedElasticityUP.cpp

@@ -388,6 +388,7 @@ int TPZHybridMixedElasticityUP::VariableIndex(const std::string& name) const {
     if(!strcmp("Stress", name.c_str())) return EStress;
     if(!strcmp("Strain", name.c_str())) return EStrain;
     if(!strcmp("VonMises", name.c_str())) return EVonMises;
+    if(!strcmp("DeviatoricStress", name.c_str())) return EDeviatoricStress;
 
     std::cout << "\n\nVar index not implemented\n\n";
     DebugStop();
@@ -409,6 +410,7 @@ int TPZHybridMixedElasticityUP::NSolutionVariables(int var) const{
             break;
         case EStress: // stress tensor
         case EStrain: // strain tensor
+        case EDeviatoricStress: // deviatoric stress tensor
             aux = 9;
             break;
         default:
@@ -555,6 +557,35 @@ void TPZHybridMixedElasticityUP::Solution(const TPZVec<TPZMaterialDataT<STATE>>&
             std::cout << "\n\nVar index not implemented\n\n";
             DebugStop();
         }
+        case EDeviatoricStress:
+        {
+            TPZFNMatrix<6, STATE> sigmavoight(n, 1, 0.0);
+            DeviatoricStressTensor(gradU, sigmavoight);
+
+            TPZFNMatrix<9, STATE> devsigma(3, 3, 0.0);
+            int cont = fdimension-1;
+            for (int i = 0; i < fdimension; i++)
+            {
+                devsigma(i,i) = sigmavoight(i,0);
+                for (int j = i+1; j < fdimension; j++)
+                {
+                    devsigma(i,j) = sigmavoight(++cont,0);
+                    devsigma(j,i) = sigmavoight(cont,0);
+                }
+            }
+            
+            if (fAnalysisType == AnalysisType::EPlaneStrain)
+                devsigma(2,2) = fthickness * fpoisson * (devsigma(0,0) + devsigma(1,1));
+
+            for (int i = 0; i < 3; i++)
+            {
+                for (int j = 0; j < 3; j++)
+                {
+                    Solout[i * 3 + j] = devsigma(i, j);
+                }
+            }
+            break;
+        }
     }
 }
 

Material/Elasticity/TPZHybridMixedElasticityUP.h

@@ -132,7 +132,7 @@ class TPZHybridMixedElasticityUP : public TPZMatBase<STATE, TPZMatCombinedSpaces
 
     [[nodiscard]] virtual int IntegrationRuleOrder(const TPZVec<int> &elPMaxOrder) const override;
 
-    enum SolutionVars {ENone = -1, EPressure = 0, EDisplacement = 1, EForce = 2, EStress = 3, EStrain = 4, EVonMises = 5};
+    enum SolutionVars {ENone = -1, EPressure = 0, EDisplacement = 1, EForce = 2, EStress = 3, EStrain = 4, EVonMises = 5, EDeviatoricStress = 6};
 };
 
 #endif