Created a base class for deformation.

include/materials/LynxDamageDeformation.h

@@ -28,7 +28,7 @@ class LynxDamageDeformation;
 template <>
 InputParameters validParams<LynxDamageDeformation>();
 
-class LynxDamageDeformation : public LynxDeformation
+class LynxDamageDeformation : public LynxDeformationBase
 {
 public:
   LynxDamageDeformation(const InputParameters & parameters);
@@ -36,18 +36,23 @@ public:
 
 protected:
   virtual void initQpStatefulProperties() override;
-  virtual void computeQpDeformation() override;
   virtual void elasticModuli() override;
   virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) override;
-  virtual void convexPlasticCorrection(Real & pressure, RankTwoTensor & stress_dev);
-  virtual Real convexReferencePressure(const Real & p_tr, const Real & q_tr);
-  virtual Real convexPlasticYield2(const Real & rho);
+  virtual Real computePlasticityYield(const Real & pressure, const Real & eqv_stress);
+  virtual Real plasticIncrement(const Real & /*pressure*/, const Real & eqv_stress);
+  virtual Real computeConvexPlasticityYield2(const Real & pressure, const Real & eqv_stress);
+  virtual Real computeConvexPlasticityYield2(const Real & rho);
+  virtual Real convexPlasticIncrement(Real & vol_plastic_incr, Real & eqv_plastic_incr);
+  virtual void computeDamageProperties(const Real & pressure, const Real & eqv_stress);
+  virtual void updateDamageParameters();
+  virtual void updateDamageConvexParameters(const Real & pressure, const Real & eqv_stress);
+  virtual Real convexReferencePressure();
   virtual Real dConvexPlasticYield2(const Real & rho);
-  virtual Real convexPlasticYield2(const Real & pressure, const Real & eqv_stress);
   virtual Real dConvexPlasticYield2_dp(const Real & pressure, const Real & eqv_stress);
   virtual Real dConvexPlasticYield2_dq(const Real & pressure, const Real & eqv_stress);
   virtual Real getConvexProjection(const Real & x1, const Real & x2);
   virtual Real strainRatio(const RankTwoTensor & elastic_strain);
+
   // Coupled variables
   bool _coupled_dam;
   const VariableValue & _damage;
@@ -55,39 +60,22 @@ protected:
   bool _coupled_phi;
   const VariableValue & _porosity;
 
-  // Strain parameters
-
   // Elastic moduli parameters
   const std::vector<Real> _damage_modulus;
 
-  // Creep parameters
-
   // Plastic parameters
-  const std::vector<Real> _critical_pressure;
+  const std::vector<Real> _friction_angle;
+  const std::vector<Real> _porous_coeff;
+  const std::vector<Real> _porous_coeff_linear;
+  std::vector<Real> _one_on_plastic_eta;
+  std::vector<Real> _one_on_damage_eta;
 
-  // Convex yield parameters
-  Real _p_cr;
-  Real _xi0;
-  Real _p_r;
-  Real _q_r;
-  Real _dp_r_dp_tr;
-  Real _dp_r_dq_tr;
-  Real _p_tr;
-  Real _q_tr;
-  Real _rho_tr;
-  Real _rp;
-  Real _rq;
-  Real _p_k;
+  // Damage-Plasticity structure
+  damage_plasticity * _damage_plasticity;
 
   // Strain properties
   MaterialProperty<RankTwoTensor> & _elastic_strain;
   const MaterialProperty<RankTwoTensor> & _elastic_strain_old;
-
-  // Viscous properties
-
-  // Plastic properties
-
-  // Stress properties
 };
 
 #endif // LYNXDAMAGEDEFORMATION_H

include/materials/LynxDeformation.h

@@ -21,105 +21,28 @@
 #ifndef LYNXDEFORMATION_H
 #define LYNXDEFORMATION_H
 
-#include "LynxMaterialBase.h"
-#include "LynxRheologyStructures.h"
-#include "RankTwoTensor.h"
-#include "RankFourTensor.h"
+#include "LynxDeformationBase.h"
 
 class LynxDeformation;
 
 template <>
 InputParameters validParams<LynxDeformation>();
 
-class LynxDeformation : public LynxMaterialBase
+class LynxDeformation : public LynxDeformationBase
 {
 public:
   LynxDeformation(const InputParameters & parameters);
   virtual ~LynxDeformation() {}
-  static MooseEnum strainModel();
-  static MooseEnum viscousUpdate();
 
 protected:
   virtual void initQpStatefulProperties() override;
-  virtual void computeProperties() override;
-  virtual void computeQpProperties() override;
-  virtual void computeStrainIncrement();
-  virtual void calculateSmallStrain(const RankTwoTensor & grad_tensor,
-                                    const RankTwoTensor & grad_tensor_old);
-  virtual void calculateFiniteStrain(const RankTwoTensor & grad_tensor,
-                                     const RankTwoTensor & grad_tensor_old);
-  virtual void computeQpDeformation();
-  virtual void elasticModuli();
-  virtual void volumetricDeformation(Real & pressure);
-  virtual void deviatoricDeformation(const Real & pressure, RankTwoTensor & stress_dev);
-  virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev);
-  virtual void tangentOperator();
-  virtual void finiteTangentOperator();
-  virtual Real computeStokeEffectiveViscosity(const Real & pressure);
-  virtual RankTwoTensor computeStokeEffectiveViscosityDerivative();
-  virtual Real computeEffectiveViscosity(const Real & pressure, const Real & eqv_stress);
-  virtual Real computeCreepRate(const Real A, const Real n, const Real eqv_stress);
-  virtual Real viscousIncrement(const Real & pressure, const Real & eqv_stress);
+  virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) override;
   virtual Real computePlasticityYield(const Real & pressure, const Real & eqv_stress);
-  virtual Real plasticIncrement(const Real & pressure, const Real & eqv_stress);
-  virtual RankTwoTensor spinRotation(const RankTwoTensor & tensor);
-  virtual void updateSpinTangentModulus();
-  virtual void computeCreepProperties(const Real & pressure);
-  virtual void updateCreepParameters();
+  virtual Real plasticIncrement(const Real & /*pressure*/, const Real & eqv_stress);
   virtual void computePlasticityProperties(const Real & pressure);
   virtual void updatePlasticityParameters();
-  virtual Real rootBrent(iterative_viscous & viscous_model, const Real x1, const Real x2);
-  virtual Real rootNewtonSafe(iterative_viscous & viscous_model, const Real x1, const Real x2);
-  void computeQpThermalSources();
-
-  // Coupled variables
-  unsigned int _ndisp;
-  std::vector<const VariableGradient *> _grad_disp;
-  std::vector<const VariableGradient *> _grad_disp_old;
-  bool _coupled_temp;
-  const VariableValue & _temp;
-  bool _coupled_plith;
-  const VariableValue & _plith;
-  const VariableValue & _plith_old;
-  bool _coupled_pdyn;
-  const VariableValue & _pdyn;
-
-  // Strain parameters
-  const MooseEnum _strain_model;
-  bool _vol_locking_correction;
-  std::vector<RankTwoTensor> _deformation_gradient;
-  std::vector<RankTwoTensor> _deformation_gradient_old;
-  const Real & _current_elem_volume;
-
-  // Elastic moduli parameters
-  const bool _has_bulk_modulus;
-  const bool _has_shear_modulus;
-  const std::vector<Real> _bulk_modulus;
-  const std::vector<Real> _shear_modulus;
-
-  // Creep parameters
-  const bool _has_diffusion_creep;
-  const std::vector<Real> _A_diffusion;
-  const std::vector<Real> _E_diffusion;
-  const std::vector<Real> _V_diffusion;
-  const bool _has_dislocation_creep;
-  const std::vector<Real> _A_dislocation;
-  const std::vector<Real> _E_dislocation;
-  const std::vector<Real> _V_dislocation;
-  const std::vector<Real> _n_dislocation;
-  const Real _gas_constant;
-  const MooseEnum _viscous_update;
-  const Real _tol;
-  const unsigned int _itmax;
-  const bool _has_background_strain_rate;
-  bool _has_initial_viscosity;
-  const std::vector<Real> _initial_viscosity;
-  const Real _background_strain_rate;
-  const std::vector<Real> _eta_min;
-  const std::vector<Real> _eta_max;
 
   // Plastic parameters
-  const bool _has_plasticity;
   const std::vector<Real> _friction_angle_0;
   const std::vector<Real> _cohesion_0;
   const std::vector<Real> _friction_angle_res;
@@ -130,52 +53,12 @@ protected:
   std::vector<Real> _one_on_plastic_eta;
   const bool _has_hardening;
 
-  // Rheology boolean
-  const bool _has_elastic;
-  const bool _has_viscous;
-
-  // Creep structures
-  diffusion_creep * _diffusion_creep;
-  dislocation_creep * _dislocation_creep;
-  iterative_viscous * _iterative_viscous;
-
   // Plasticity structure
   plasticity * _plasticity;
 
-  // Viscous scalars for tangent modulus
-  Real _stress_corr_v;
-  Real _dq_dq_tr_v;
-
-  // Plastic scalars for tangent modulus
-  Real _stress_corr_p;
-  Real _dp_dp_tr_p;
-  Real _dp_dq_tr_p;
-  Real _dq_dp_tr_p;
-  Real _dq_dq_tr_p;
-
-  // Strain properties
-  MaterialProperty<RankTwoTensor> & _strain_increment;
-  MaterialProperty<RankTwoTensor> & _spin_tensor;
-  const MaterialProperty<RankTwoTensor> & _thermal_strain_incr;
-
-  // Viscous properties
-  MaterialProperty<RankTwoTensor> & _viscous_strain_incr;
-  MaterialProperty<Real> & _eta_eff;
-
   // Plastic properties
-  MaterialProperty<RankTwoTensor> & _plastic_strain_incr;
-  MaterialProperty<Real> & _plastic_yield_function;
   MaterialProperty<Real> * _intnl;
   const MaterialProperty<Real> * _intnl_old;
-
-  // Stress properties
-  MaterialProperty<RankTwoTensor> & _stress;
-  const MaterialProperty<RankTwoTensor> & _stress_old;
-  MaterialProperty<Real> & _K;
-  MaterialProperty<Real> & _G;
-  MaterialProperty<RankFourTensor> & _tangent_modulus;
-  MaterialProperty<Real> & _inelastic_heat;
-  MaterialProperty<RankTwoTensor> & _dinelastic_heat_dstrain;
 };
 
 #endif // LYNXDEFORMATION_H

include/materials/LynxDeformationBase.h

+/******************************************************************************/
+/*                       LYNX, a MOOSE-based application                      */
+/*                                                                            */
+/*          Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace         */
+/*             GFZ Potsdam, German Research Centre for Geosciences            */
+/*                                                                            */
+/*    This program 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.                    */
+/*                                                                            */
+/*       This program 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 this program. If not, see <http://www.gnu.org/licenses/>     */
+/******************************************************************************/
+
+#ifndef LYNXDEFORMATIONBASE_H
+#define LYNXDEFORMATIONBASE_H
+
+#include "LynxMaterialBase.h"
+#include "LynxRheologyStructures.h"
+#include "RankTwoTensor.h"
+#include "RankFourTensor.h"
+
+class LynxDeformationBase;
+
+template <>
+InputParameters validParams<LynxDeformationBase>();
+
+class LynxDeformationBase : public LynxMaterialBase
+{
+public:
+  LynxDeformationBase(const InputParameters & parameters);
+  virtual ~LynxDeformationBase() {}
+  static MooseEnum strainModel();
+  static MooseEnum viscousUpdate();
+
+protected:
+  virtual void initQpStatefulProperties() override;
+  virtual void computeProperties() override;
+  virtual void computeQpProperties() override;
+  virtual void computeStrainIncrement();
+  virtual void calculateSmallStrain(const RankTwoTensor & grad_tensor,
+                                    const RankTwoTensor & grad_tensor_old);
+  virtual void calculateFiniteStrain(const RankTwoTensor & grad_tensor,
+                                     const RankTwoTensor & grad_tensor_old);
+  virtual void computeQpDeformation();
+  virtual void initializeQpDeformation();
+  virtual void elasticModuli();
+  virtual void volumetricDeformation(Real & pressure);
+  virtual void deviatoricDeformation(const Real & pressure, RankTwoTensor & stress_dev);
+  virtual void reformStressTensor(const Real & pressure, const RankTwoTensor & stress_dev);
+  virtual void plasticCorrection(Real & pressure, RankTwoTensor & stress_dev) = 0;
+  virtual void tangentOperator();
+  virtual void plasticTangentOperator(const RankTwoTensor & flow_direction,
+                                      const RankFourTensor & flow_direction_dyad);
+  virtual void finiteTangentOperator();
+  virtual Real computeStokeEffectiveViscosity(const Real & pressure);
+  virtual RankTwoTensor computeStokeEffectiveViscosityDerivative();
+  virtual Real computeEffectiveViscosity(const Real & pressure, const Real & eqv_stress);
+  virtual Real computeCreepRate(const Real A, const Real n, const Real eqv_stress);
+  virtual Real viscousIncrement(const Real & pressure, const Real & eqv_stress);
+  virtual RankTwoTensor spinRotation(const RankTwoTensor & tensor);
+  virtual void updateSpinTangentModulus();
+  virtual void computeCreepProperties(const Real & pressure);
+  virtual void updateCreepParameters();
+  virtual Real rootBrent(iterative_viscous & viscous_model, const Real x1, const Real x2);
+  virtual Real rootNewtonSafe(iterative_viscous & viscous_model, const Real x1, const Real x2);
+  void computeQpThermalSources();
+
+  // Coupled variables
+  unsigned int _ndisp;
+  std::vector<const VariableGradient *> _grad_disp;
+  std::vector<const VariableGradient *> _grad_disp_old;
+  bool _coupled_temp;
+  const VariableValue & _temp;
+  bool _coupled_plith;
+  const VariableValue & _plith;
+  const VariableValue & _plith_old;
+  bool _coupled_pdyn;
+  const VariableValue & _pdyn;
+
+  // Strain parameters
+  const MooseEnum _strain_model;
+  bool _vol_locking_correction;
+  std::vector<RankTwoTensor> _deformation_gradient;
+  std::vector<RankTwoTensor> _deformation_gradient_old;
+  const Real & _current_elem_volume;
+
+  // Elastic moduli parameters
+  const bool _has_bulk_modulus;
+  const bool _has_shear_modulus;
+  const std::vector<Real> _bulk_modulus;
+  const std::vector<Real> _shear_modulus;
+
+  // Creep parameters
+  const bool _has_diffusion_creep;
+  const std::vector<Real> _A_diffusion;
+  const std::vector<Real> _E_diffusion;
+  const std::vector<Real> _V_diffusion;
+  const bool _has_dislocation_creep;
+  const std::vector<Real> _A_dislocation;
+  const std::vector<Real> _E_dislocation;
+  const std::vector<Real> _V_dislocation;
+  const std::vector<Real> _n_dislocation;
+  const Real _gas_constant;
+  const MooseEnum _viscous_update;
+  const Real _tol;
+  const unsigned int _itmax;
+  const bool _has_background_strain_rate;
+  bool _has_initial_viscosity;
+  const std::vector<Real> _initial_viscosity;
+  const Real _background_strain_rate;
+  const std::vector<Real> _eta_min;
+  const std::vector<Real> _eta_max;
+
+  // Rheology boolean
+  const bool _has_elastic;
+  const bool _has_viscous;
+  bool _has_plasticity;
+
+  // Creep structures
+  diffusion_creep * _diffusion_creep;
+  dislocation_creep * _dislocation_creep;
+  iterative_viscous * _iterative_viscous;
+
+  // RankFourTensor utilities
+  RankTwoTensor _identity_two;
+  RankFourTensor _volumetric_four;
+  RankFourTensor _identity_four;
+  RankFourTensor _deviatoric_four;
+
+  // Viscous scalars for tangent modulus
+  Real _stress_corr_v;
+  Real _dq_dq_tr_v;
+
+  // Plastic scalars for tangent modulus
+  Real _stress_corr_p;
+  Real _dp_dp_tr_p;
+  Real _dp_dq_tr_p;
+  Real _dq_dp_tr_p;
+  Real _dq_dq_tr_p;
+
+  // Strain properties
+  MaterialProperty<RankTwoTensor> & _strain_increment;
+  MaterialProperty<RankTwoTensor> & _spin_tensor;
+  const MaterialProperty<RankTwoTensor> & _thermal_strain_incr;
+
+  // Viscous properties
+  MaterialProperty<RankTwoTensor> & _viscous_strain_incr;
+  MaterialProperty<Real> & _eta_eff;
+
+  // Plastic properties
+  MaterialProperty<RankTwoTensor> & _plastic_strain_incr;
+  MaterialProperty<Real> & _plastic_yield_function;
+
+  // Stress properties
+  MaterialProperty<RankTwoTensor> & _stress;
+  const MaterialProperty<RankTwoTensor> & _stress_old;
+  MaterialProperty<Real> & _K;
+  MaterialProperty<Real> & _G;
+  MaterialProperty<RankFourTensor> & _tangent_modulus;
+  MaterialProperty<Real> & _inelastic_heat;
+  MaterialProperty<RankTwoTensor> & _dinelastic_heat_dstrain;
+};
+
+#endif // LYNXDEFORMATIONBASE_H

include/utils/LynxRheologyStructures.h

@@ -206,4 +206,55 @@ struct plasticity
   }
 };
 
-#endif // LYNXRHEOLOGYSTRUCTURES_H
\ No newline at end of file
+struct damage_plasticity
+{
+  Real _xi0;
+  Real _gamma;
+  Real _p_cr;
+  Real _eta_p; // actually one on eta
+  Real _eta_d; // actually one on eta
+  Real _alpha0;
+  Real _p_tr;        // trial pressure
+  Real _q_tr;        // trial eqv_stress
+  Real _p_r;         // reference pressure for convex yield
+  Real _q_r;         // reference eqv_stress for convex yield
+  Real _rho_tr;      // trial distance to reference point
+  Real _rp;          // direction for projection
+  Real _rq;          // direction for projection
+  Real _xi_cr;       // critical strain ratio for capped yield
+  Real _alpha2;      // square of alpha parameter for capped yield
+  Real _dxi_cr_dp;   // derivative wrt pressure of the critical strain ratio
+  Real _dxi_cr_dq;   // derivative wrt eqv_stress of the critical strain ratio
+  Real _dmu2_dxi_cr; // derivative wrt the critical strain ratio of the square of alpha parameter
+  damage_plasticity()
+    : _xi0(-std::sqrt(3.0)),
+      _gamma(0.0),
+      _p_cr(0.0),
+      _eta_p(0.0),
+      _eta_d(0.0),
+      _alpha0(0.0),
+      _p_tr(0.0),
+      _q_tr(0.0),
+      _p_r(0.0),
+      _q_r(0.0),
+      _rho_tr(0.0),
+      _rp(0.0),
+      _rq(0.0),
+      _xi_cr(-std::sqrt(3.0)),
+      _alpha2(0.0),
+      _dxi_cr_dp(0.0),
+      _dxi_cr_dq(0.0),
+      _dmu2_dxi_cr(0.0)
+  {
+  }
+  void fill(const Real xi0, const Real gamma, const Real p_cr, const Real eta_p, const Real eta_d)
+  {
+    _xi0 = xi0;
+    _gamma = gamma;
+    _p_cr = p_cr;
+    _eta_p = eta_p;
+    _eta_d = eta_d;
+  }
+};
+
+#endif // LYNXRHEOLOGYSTRUCTURES_H

test/tests/plastic/convex_compression.i

@@ -138,7 +138,8 @@
     bulk_modulus = 2.0e+09
     shear_modulus = 2.0e+09
     friction_angle = 30
-    critical_pressure = 8.0e+07
+    porous_coeff = 2.1246 # critical pressure 8.0e+07
+    # critical_pressure = 8.0e+07
   [../]
 []
 
@@ -185,4 +186,3 @@
   exodus = true
   csv = true
 []
-

test/tests/plastic/convex_compression_vp.i

@@ -138,7 +138,8 @@
     bulk_modulus = 2.0e+09
     shear_modulus = 2.0e+09
     friction_angle = 30
-    critical_pressure = 8.0e+07
+    porous_coeff = 2.1246 # critical pressure 8.0e+07
+    # critical_pressure = 8.0e+07
     plastic_viscosity = 5.0e+10
   [../]
 []
@@ -186,4 +187,3 @@
   exodus = true
   csv = true
 []
-

test/tests/plastic/convex_oedometer.i

-[Mesh]
-  type = GeneratedMesh
-  dim = 3
-  nx = 4
-  ny = 4
-  nz = 4
-  xmin = 0
-  xmax = 1
-  ymin = 0
-  ymax = 1
-  zmin = 0
-  zmax = 1
-[]
-
-[GlobalParams]
-  displacements = 'disp_x disp_y disp_z'
-[]
-
-[Variables]
-  [./disp_x]
-    order = FIRST
-    family = LAGRANGE
-  [../]
-  [./disp_y]
-    order = FIRST
-    family = LAGRANGE
-  [../]
-  [./disp_z]
-    order = FIRST
-    family = LAGRANGE
-  [../]
-[]
-
-[Kernels]
-  [./mech_x]
-    type = LynxSolidMomentum
-    variable = disp_x
-    component = 0
-  [../]
-  [./mech_y]
-    type = LynxSolidMomentum
-    variable = disp_y
-    component = 1
-  [../]
-  [./mech_z]
-    type = LynxSolidMomentum
-    variable = disp_z
-    component = 2
-  [../]
-[]
-
-[AuxVariables]
-  [./strain_xx]
-    order = CONSTANT
-    family = MONOMIAL
-  [../]
-  [./plastic_strain_xx]
-    order = CONSTANT
-    family = MONOMIAL
-  [../]
-  [./stress_xx]
-    order = CONSTANT
-    family = MONOMIAL
-  [../]
-  [./yield]
-    order = CONSTANT
-    family = MONOMIAL
-  [../]
-[]
-
-[AuxKernels]
-  [./strain_xx_aux]
-    type = LynxStrainAux
-    variable = strain_xx
-    index_i = 0
-    index_j = 0
-  [../]
-  [./plastic_strain_xx_aux]
-    type = LynxStrainAux
-    variable = plastic_strain_xx
-    strain_type = plastic
-    index_i = 0
-    index_j = 0
-  [../]
-  [./stress_xx_aux]
-    type = LynxStressAux
-    variable = stress_xx
-    index_i = 0
-    index_j = 0
-  [../]
-  [./yield_aux]
-    type = MaterialRealAux
-    variable = yield
-    property = plastic_yield_function
-  [../]
-[]
-
-[BCs]
-  [./no_x_left]
-    type = PresetBC
-    variable = disp_x
-    boundary = left
-    value = 0.0
-  [../]
-  [./load_x_right]
-    type = LynxVelocityBC
-    variable = disp_x
-    boundary = right
-    value = -1.0e-05
-  [../]
-  [./no_y]
-    type = PresetBC
-    variable = disp_y
-    boundary = 'bottom top'
-    value = 0.0
-  [../]
-  [./no_z_back]