Fix stress used in principal_stress.cc

geodynamics · Nov 7, 2023 · ce7921a · ce7921a
1 parent 2cf3bc7
commit ce7921a
Showing 1 changed file with 21 additions and 20 deletions.
diff --git a/source/postprocess/visualization/principal_stress.cc b/source/postprocess/visualization/principal_stress.cc
@@ -115,36 +115,37 @@ namespace aspect
 
         for (unsigned int q=0; q<n_quadrature_points; ++q)
           {
-            const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
-            const SymmetricTensor<2,dim> deviatoric_strain_rate
-              = (this->get_material_model().is_compressible()
-                 ?
-                 strain_rate - 1./3. * trace(strain_rate) * unit_symmetric_tensor<dim>()
-                 :
-                 strain_rate);
-
-            const double eta = out.viscosities[q];
-
-            // Compressive stress is positive in geoscience applications
-            SymmetricTensor<2,dim> stress = -2. * eta * deviatoric_strain_rate + in.pressure[q] * unit_symmetric_tensor<dim>();
+            SymmetricTensor<2,dim> stress;
 
             // Add elastic stresses if existent
             if (this->get_parameters().enable_elasticity == true)
               {
-                stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
-                stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
-                stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
+                stress[0][0] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
+                stress[1][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
+                stress[0][1] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
 
                 if (dim == 3)
                   {
-                    stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
-                    stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
-                    stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
+                    stress[2][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
+                    stress[0][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
+                    stress[1][2] = in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
                   }
               }
+            else
+              {
+                const SymmetricTensor<2,dim> deviatoric_strain_rate
+                  = (this->get_material_model().is_compressible()
+                     ?
+                     in.strain_rate[q] - 1./3. * trace(in.strain_rate[q]) * unit_symmetric_tensor<dim>()
+                     :
+                     in.strain_rate[q]);
+
+                // Compressive stress is positive in geoscience applications
+                stress = -2. * out.viscosities[q] * deviatoric_strain_rate;
+              }
 
-            if (use_deviatoric_stress == true)
-              stress -= 1./dim * trace(stress) * unit_symmetric_tensor<dim>();
+            if (use_deviatoric_stress == false)
+              stress += in.pressure[q] * unit_symmetric_tensor<dim>();
 
             const std::array<std::pair<double, Tensor<1,dim>>, dim> principal_stresses_and_directions = eigenvectors(stress);