Plastic strain softening

docs/src/man/customrheology.md

@@ -1,4 +1,4 @@
-# `User-defined rheology`
+## `User-defined rheology`
 
 `CustomRheology` allows the user to interface with `GeoParams.jl` API for rheology calculations: 
 
@@ -13,11 +13,10 @@ end
 where `strain` and `stress` are functions that compute the strain rate and deviatoric stress tensors, respectively, and `args` is a `NamedTuple` containing any constant (i.e. **immutable**) physical parameters needed by our `CustomRheology`. 
 
 ## Example: Arrhenious-type rheology
-The deviatoric stress $\tau_{ij}$ and strain rate tensors $\dot{\varepsilon}_{ij}$ are simply
-    $$\tau_{ij} = 2 \eta \dot{\varepsilon}_{ij} \qquad \dot{\varepsilon}_{ij}  = {\tau_{ij}  \over 2 \eta}$$
+The deviatoric stress $\tau_{ij}$ and strain rate tensors $\dot{\varepsilon}_{ij}$ are simply $\tau_{ij} = 2 \eta \dot{\varepsilon}_{ij}$ and $\dot{\varepsilon}_{ij}  = {\tau_{ij}  \over 2 \eta}$
 
 and the viscosity $\eta$ is temperature-dependant
-    $$\eta = \eta_{0} * \exp\left(\frac{E}{T+T_{o}}-\frac{E}{T_{\eta}+T_{o}}\right)$$
+    $$\eta = \eta_{0}  \exp\left(\frac{E}{T+T_{o}}-\frac{E}{T_{\eta}+T_{o}}\right)$$
 
 where $\eta_0$ and $T_{\eta}$ are the respective reference viscosity and temperature, $T_o$ is the offset temperature, $T$ is the local temperature, and $E$ is activation energy. 
 

src/ConstitutiveRelationships.jl

@@ -22,10 +22,11 @@ const AxialCompression, SimpleShear, Invariant = 1, 2, 3
 @inline precision(::AbstractConstitutiveLaw{T}) where T = T
 
 include("Computations.jl")
-include("CreepLaw/CreepLaw.jl")              # viscous Creeplaws
-include("Elasticity/Elasticity.jl")          # elasticity
-include("Plasticity/Plasticity.jl")          # plasticity
-include("CompositeRheologies/CompositeRheologies.jl")            # composite constitutive relationships
+include("Softening/Softening.jl")                      # strain softening
+include("CreepLaw/CreepLaw.jl")                        # viscous Creeplaws
+include("Elasticity/Elasticity.jl")                    # elasticity
+include("Plasticity/Plasticity.jl")                    # plasticity
+include("CompositeRheologies/CompositeRheologies.jl")  # composite constitutive relationships
 
 export param_info,
     dεII_dτII,
@@ -58,7 +59,9 @@ export param_info,
     AxialCompression, SimpleShear, Invariant, 
     get_G, 
     get_Kb,
-    iselastic
+    iselastic,
+    NoSoftening,
+    LinearSoftening
 
 # add methods programmatically 
 for myType in (:LinearViscous, :DiffusionCreep, :DislocationCreep, :ConstantElasticity, :DruckerPrager, :ArrheniusType, 

src/GeoParams.jl

@@ -196,6 +196,10 @@ export dεII_dτII,
     get_shearmodulus,
     get_bulkmodulus,
 
+    #       softening
+    NoSoftening,
+    LinearSoftening,
+
     #       Plasticity
     AbstractPlasticity,
     compute_yieldfunction,

src/Plasticity/DruckerPrager.jl

@@ -27,7 +27,9 @@ Plasticity is activated when ``F(\\tau_{II}^{trial})`` (the yield function compu
 where ``\\dot{\\lambda}`` is a (scalar) that is nonzero and chosen such that the resulting stress gives ``F(\\tau_{II}^{final})=0``, and ``\\sigma_{ij}=-P + \\tau_{ij}`` denotes the total stress tensor.   
 
 """
-@with_kw_noshow struct DruckerPrager{T, U, U1} <: AbstractPlasticity{T}
+@with_kw_noshow struct DruckerPrager{T, U, U1, S1<:AbstractSoftening, S2<:AbstractSoftening} <: AbstractPlasticity{T}
+    softening_ϕ::S1 = NoSoftening()
+    softening_C::S2 = NoSoftening()
     ϕ::GeoUnit{T,U} = 30NoUnits # Friction angle
     Ψ::GeoUnit{T,U} = 0NoUnits # Dilation angle
     sinϕ::GeoUnit{T,U} = sind(ϕ)NoUnits # Friction angle
@@ -36,7 +38,9 @@ where ``\\dot{\\lambda}`` is a (scalar) that is nonzero and chosen such that the
     cosΨ::GeoUnit{T,U} = cosd(Ψ)NoUnits # Dilation angle
     C::GeoUnit{T,U1} = 10e6Pa # Cohesion
 end
-DruckerPrager(args...) = DruckerPrager(convert.(GeoUnit, args)...)
+
+DruckerPrager(args...) = DruckerPrager(args[1:2]..., convert.(GeoUnit, args[3:end])...)
+DruckerPrager(softening_ϕ::AbstractSoftening, softening_C::AbstractSoftening, args...) = DruckerPrager(softening_ϕ, softening_C, convert.(GeoUnit, args)...)
 
 function isvolumetric(s::DruckerPrager)
     @unpack_val Ψ = s
@@ -50,10 +54,46 @@ function param_info(s::DruckerPrager) # info about the struct
 end
 
 # Calculation routines
-function (s::DruckerPrager{_T,U,U1})(;
+function (s::DruckerPrager{_T, U, U1, S, S})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T,U,U1,S<:AbstractSoftening}
+    @unpack_val sinϕ, cosϕ, ϕ, C = s
+    ϕ = s.softening_ϕ(ϕ, EII)
+    C = s.softening_C(C, EII)
+
+    cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))
+
+    F = τII - cosϕ * C - sinϕ * (P - Pf)   # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager{_T, U, U1, NoSoftening, S})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T,U,U1,S}
+    @unpack_val sinϕ, cosϕ, ϕ, C = s
+    C = s.softening_C(C, EII)
+
+    F = τII - cosϕ * C - sinϕ * (P - Pf)   # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager{_T, U, U1, S, NoSoftening})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T,U,U1,S}
+    @unpack_val sinϕ, cosϕ, ϕ, C = s
+    ϕ = s.softening_ϕ(ϕ, EII)
+
+    cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))
+
+    F = τII - cosϕ * C - sinϕ * (P - Pf)   # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager{_T, U, U1, NoSoftening, NoSoftening})(;
     P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), kwargs...
 ) where {_T,U,U1}
     @unpack_val sinϕ, cosϕ, ϕ, C = s
+
     F = τII - cosϕ * C - sinϕ * (P - Pf)   # with fluid pressure (set to zero by default)
     return F
 end
@@ -64,9 +104,9 @@ end
 Computes the plastic yield function `F` for a given second invariant of the deviatoric stress tensor `τII`,  `P` pressure, and `Pf` fluid pressure.
 """
 function compute_yieldfunction(
-    s::DruckerPrager{_T}; P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T)
+    s::DruckerPrager{_T}; P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), EII::_T=zero(_T)
 ) where {_T}
-    return s(; P=P, τII=τII, Pf=Pf)
+    return s(; P=P, τII=τII, Pf=Pf, EII=EII)
 end
 
 """
@@ -81,11 +121,12 @@ function compute_yieldfunction!(
     s::DruckerPrager{_T};
     P::AbstractArray{_T,N},
     τII::AbstractArray{_T,N},
-    Pf=zero(P)::AbstractArray{_T,N},
+    Pf::AbstractArray{_T,N} = zero(P),
+    EII::AbstractArray{_T,N} = zero(P),
     kwargs...,
 ) where {N,_T}
     @inbounds for i in eachindex(P)
-        F[i] = compute_yieldfunction(s; P=P[i], τII=τII[i], Pf=Pf[i])
+        F[i] = compute_yieldfunction(s; P=P[i], τII=τII[i], Pf=Pf[i], EII=EII[i])
     end
 
     return nothing
@@ -164,7 +205,8 @@ end
     `K` is the elastic bulk modulus, h is the harderning, and `τij`` is the stress tensor in Voigt notation.
     Equations from Duretz et al. 2019 G3
 """
-@inline function lambda(F::T, p::DruckerPrager, ηve::T, ηvp::T; K=zero(T), dt=zero(T), h=zero(T), τij=(one(T), one(T), one(T))) where T
-    F * inv(ηve + ηvp + K * dt * p.sinΨ * p.sinϕ + h * p.cosϕ * plastic_strain(p, τij, zero(T)))
+@inline function lambda(F::T, p::DruckerPrager, ηve::T, ηvp::T; K=zero(T), dt=zero(T), EII=zero(T), τij=(one(T), one(T), one(T))) where T
+    ϕ = s.softening_ϕ(p.cosϕ, EII)
+    F * inv(ηve + ηvp + K * dt * p.sinΨ * p.sinϕ + h * cosdϕ * plastic_strain(p, τij, zero(T)))
 end
 #-------------------------------------------------------------------------

src/Plasticity/DruckerPrager_regularised.jl

@@ -27,17 +27,21 @@ Plasticity is activated when ``F(\\tau_{II}^{trial})`` (the yield function compu
 where ``\\dot{\\lambda}`` is a (scalar) that is nonzero and chosen such that the resulting stress gives ``F(\\tau_{II}^{final})=0``, and ``\\sigma_{ij}=-P + \\tau_{ij}`` denotes the total stress tensor.   
 
 """
-@with_kw_noshow struct DruckerPrager_regularised{T, U, U1, U2} <: AbstractPlasticity{T}
-    ϕ::GeoUnit{T,U} = 30NoUnits # Friction angle
-    Ψ::GeoUnit{T,U} = 0NoUnits # Dilation angle
+@with_kw_noshow struct DruckerPrager_regularised{T, U, U1, U2, S1<:AbstractSoftening, S2<:AbstractSoftening} <: AbstractPlasticity{T}
+    softening_ϕ::S1 = NoSoftening()
+    softening_C::S2 = NoSoftening()
+    ϕ::GeoUnit{T,U} = 30NoUnits         # Friction angle
+    Ψ::GeoUnit{T,U} = 0NoUnits          # Dilation angle
     sinϕ::GeoUnit{T,U} = sind(ϕ)NoUnits # Friction angle
     cosϕ::GeoUnit{T,U} = cosd(ϕ)NoUnits # Friction angle
     sinΨ::GeoUnit{T,U} = sind(Ψ)NoUnits # Dilation angle
     cosΨ::GeoUnit{T,U} = cosd(Ψ)NoUnits # Dilation angle
-    C::GeoUnit{T,U1} = 10e6Pa # Cohesion
-    η_vp::GeoUnit{T,U2} = 1e20Pa*s # regularisation viscosity
+    C::GeoUnit{T,U1} = 10e6Pa           # Cohesion
+    η_vp::GeoUnit{T,U2} = 1e20Pa*s      # regularisation viscosity
 end
-DruckerPrager_regularised(args...) = DruckerPrager_regularised(convert.(GeoUnit, args)...)
+
+DruckerPrager_regularised(args...) = DruckerPrager_regularised(args[1:2]..., convert.(GeoUnit, args[3:end])...)
+DruckerPrager_regularised(softening_ϕ::AbstractSoftening, softening_C::AbstractSoftening, args::Vararg{GeoUnit, N}) where N = DruckerPrager_regularised(softening_ϕ, softening_C, convert.(GeoUnit, args)...)
 
 function isvolumetric(s::DruckerPrager_regularised)
     @unpack_val Ψ = s
@@ -51,24 +55,58 @@ function param_info(s::DruckerPrager_regularised) # info about the struct
 end
 
 # Calculation routines
-function (s::DruckerPrager_regularised{_T,U,U1})(;
-    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T= zero(_T),kwargs...
+function (s::DruckerPrager_regularised{_T})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T= zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T}
+    @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
+    ϕ = s.softening_ϕ(ϕ, EII)
+    C = s.softening_C(C, EII)
+    cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))
+    ε̇II_pl = λ*∂Q∂τII(s, τII)  # plastic strainrate
+    F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2 * η_vp * ε̇II_pl # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening,S})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T= zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T,U,U1,S}
+    @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
+    C = s.softening_C(C, EII)
+    ε̇II_pl = λ*∂Q∂τII(s, τII)  # plastic strainrate
+    F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2 * η_vp * ε̇II_pl # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager_regularised{_T,U,U1,S,NoSoftening})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T= zero(_T), EII::_T=zero(_T), kwargs...
+) where {_T,U,U1,S}
+    @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
+    ϕ = s.softening_ϕ(ϕ, EII)
+    cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))
+    ε̇II_pl = λ*∂Q∂τII(s, τII)  # plastic strainrate
+    F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2 * η_vp * ε̇II_pl # with fluid pressure (set to zero by default)
+    return F
+end
+
+function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening,NoSoftening})(;
+    P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T= zero(_T), kwargs...
 ) where {_T,U,U1}
     @unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
     ε̇II_pl = λ*∂Q∂τII(s, τII)  # plastic strainrate
-    F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2*η_vp*ε̇II_pl # with fluid pressure (set to zero by default)
+    F = τII - cosϕ * C - sinϕ * (P - Pf)  - 2 * η_vp * ε̇II_pl # with fluid pressure (set to zero by default)
     return F
 end
 
+
 """
     compute_yieldfunction(s::DruckerPrager_regularised; P, τII, Pf, λ, kwargs...) 
 
 Computes the plastic yield function `F` for a given second invariant of the deviatoric stress tensor `τII`,  `P` pressure, and `Pf` fluid pressure.
 """
 function compute_yieldfunction(
-    s::DruckerPrager_regularised{_T}; P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T=zero(_T)
+    s::DruckerPrager_regularised{_T}; P::_T=zero(_T), τII::_T=zero(_T), Pf::_T=zero(_T), λ::_T=zero(_T), EII::_T=zero(_T)
 ) where {_T}
-    return s(; P=P, τII=τII, Pf=Pf, λ=λ)
+    return s(; P=P, τII=τII, Pf=Pf, λ=λ, EII=EII)
 end
 
 """
@@ -85,10 +123,11 @@ function compute_yieldfunction!(
     τII::AbstractArray{_T,N},
     Pf=zero(P)::AbstractArray{_T,N},
     λ=zero(P)::AbstractArray{_T,N},
+    EII::AbstractArray{_T,N} = zero(P),
     kwargs...,
 ) where {N,_T}
     @inbounds for i in eachindex(P)
-        F[i] = compute_yieldfunction(s; P=P[i], τII=τII[i], Pf=Pf[i], λ=λ)
+        F[i] = compute_yieldfunction(s; P=P[i], τII=τII[i], Pf=Pf[i], λ=λ[i], EII=EII[i])
     end
 
     return nothing
@@ -100,9 +139,9 @@ end
 ∂Q∂P(p::DruckerPrager_regularised, args; kwargs...) = -NumValue(p.sinΨ)
 
 # Derivatives of yield function
-∂F∂τII(p::DruckerPrager_regularised, τII::_T; P=zero(_T), kwargs...) where _T  = _T(1)
-∂F∂P(p::DruckerPrager_regularised, P::_T; τII=zero(_T), kwargs...) where _T    = -NumValue(p.sinϕ)
-∂F∂λ(p::DruckerPrager_regularised, τII::_T; P=zero(_T), kwargs...) where _T    = -2*NumValue(p.η_vp)*∂Q∂τII(p, τII, P=P) 
+∂F∂τII(p::DruckerPrager_regularised, τII::_T; kwargs...) where _T             = _T(1)
+∂F∂P(p::DruckerPrager_regularised,     P::_T; kwargs...) where _T             = -NumValue(p.sinϕ)
+∂F∂λ(p::DruckerPrager_regularised,   τII::_T; P=zero(_T), kwargs...) where _T = -2 * NumValue(p.η_vp)*∂Q∂τII(p, τII, P=P) 
 
 
 # Derivatives w.r.t stress tensor
@@ -111,16 +150,16 @@ end
 for t in (:NTuple,:SVector)
     @eval begin
         ## 3D derivatives 
-        ∂Q∂τxx(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[1] / second_invariant(τij)
-        ∂Q∂τyy(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[2] / second_invariant(τij)
-        ∂Q∂τzz(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[3] / second_invariant(τij)
-        ∂Q∂τyz(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[4] / second_invariant(τij)
-        ∂Q∂τxz(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[5] / second_invariant(τij)
-        ∂Q∂τxy(p::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[6] / second_invariant(τij) 
+        ∂Q∂τxx(::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[1] / second_invariant(τij)
+        ∂Q∂τyy(::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[2] / second_invariant(τij)
+        ∂Q∂τzz(::DruckerPrager_regularised, τij::$(t){6, T}) where T = 0.5 * τij[3] / second_invariant(τij)
+        ∂Q∂τyz(::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[4] / second_invariant(τij)
+        ∂Q∂τxz(::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[5] / second_invariant(τij)
+        ∂Q∂τxy(::DruckerPrager_regularised, τij::$(t){6, T}) where T = τij[6] / second_invariant(τij) 
         ## 2D derivatives 
-        ∂Q∂τxx(p::DruckerPrager_regularised, τij::$(t){3, T}) where T = 0.5 * τij[1] / second_invariant(τij)
-        ∂Q∂τyy(p::DruckerPrager_regularised, τij::$(t){3, T}) where T = 0.5 * τij[2] / second_invariant(τij)
-        ∂Q∂τxy(p::DruckerPrager_regularised, τij::$(t){3, T}) where T = τij[3] / second_invariant(τij) 
+        ∂Q∂τxx(::DruckerPrager_regularised, τij::$(t){3, T}) where T = 0.5 * τij[1] / second_invariant(τij)
+        ∂Q∂τyy(::DruckerPrager_regularised, τij::$(t){3, T}) where T = 0.5 * τij[2] / second_invariant(τij)
+        ∂Q∂τxy(::DruckerPrager_regularised, τij::$(t){3, T}) where T = τij[3] / second_invariant(τij) 
     end
 end
 

src/Softening/Softening.jl

@@ -0,0 +1,29 @@
+
+abstract type AbstractSoftening end
+
+struct NoSoftening <: AbstractSoftening end
+
+@inline (softening::NoSoftening)(max_value, ::Vararg{Any, N}) where N = max_value
+
+struct LinearSoftening{T} <: AbstractSoftening
+    hi::T
+    lo::T
+    max_value::T
+    min_value::T
+    damage::T
+
+    function LinearSoftening(min_value::T, max_value::T, lo::T, hi::T) where T
+        damage = max_value + (max_value - min_value) / (hi - lo)
+        new{T}(hi, lo, max_value, min_value, damage)
+    end
+end
+
+LinearSoftening(min_max_values::NTuple{2, T}, lo_hi::NTuple{2, T}) where T = LinearSoftening(min_max_values..., lo_hi...)
+
+function (softening::LinearSoftening)(max_value, softening_var::T) where T
+
+    softening_var ≥ softening.hi && return softening.min_value
+    softening_var ≤ softening.lo && return T(max_value)
+
+    return softening_var * softening.damage
+end

test/test_Softening.jl

@@ -0,0 +1,55 @@
+using GeoParams, Test
+
+@testset "LinearSoftening" begin
+
+    # Test NoSoftening
+    soft = NoSoftening()
+    x = rand()
+    @test x === soft(x, rand())
+
+    # Test LinearSoftening
+    min_v, max_v = 15e0, 30e0
+    lo, hi = 0.0, 1.0
+
+    @test LinearSoftening(min_v, max_v, lo, hi) === LinearSoftening((min_v, max_v), (lo, hi))
+
+    soft_ϕ = LinearSoftening(min_v, max_v, lo, hi)
+
+    @test soft_ϕ(30, 1) == min_v
+    @test soft_ϕ(30, 0) == max_v
+    @test soft_ϕ(30, 0.5) == 0.5 * (min_v + max_v)
+
+    # test Drucker-Prager with softening
+    τII = 20e6
+    P = 1e6
+    args = (P=P, τII=τII)
+
+    p = DruckerPrager()
+    @test compute_yieldfunction(p, args) ≈ 1.0839745962155614e7
+
+    args = (P=P, τII=τII, EII=1e0)
+
+    p1 = DruckerPrager(; softening_ϕ = soft_ϕ)
+    p2 = DruckerPrager(; softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+    p3 = DruckerPrager(; softening_ϕ = soft_ϕ, softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+
+    @test compute_yieldfunction(p1, args) ≈ 1.0081922692006797e7
+    @test compute_yieldfunction(p2, args) ≈ 1.95e7
+    @test compute_yieldfunction(p3, args) ≈ 1.974118095489748e7
+
+    # test regularized Drucker-Prager with softening
+    p = DruckerPrager_regularised()
+    @test compute_yieldfunction(p, args) ≈ 1.0839745962155614e7
+
+    args = (P=P, τII=τII, EII=1e0)
+
+    p1 = DruckerPrager(; softening_ϕ = soft_ϕ)
+    p2 = DruckerPrager(; softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+    p3 = DruckerPrager(; softening_ϕ = soft_ϕ, softening_C = LinearSoftening((0e0, 10e6), (lo, hi)))
+
+    @test compute_yieldfunction(p1, args) ≈ 1.0081922692006797e7
+    @test compute_yieldfunction(p2, args) ≈ 1.95e7
+    @test compute_yieldfunction(p3, args) ≈ 1.974118095489748e7
+
+end
+