Initializing improved differentiation interface

.github/workflows/CI.yml

@@ -19,7 +19,7 @@ jobs:
       matrix:
         version:
           - '1'
-          - '1.8'
+          - '1.11'
         os:
           - ubuntu-latest
           - macOS-latest

.typos.toml

@@ -1,2 +1,5 @@
 [files]
 extend-exclude = ["*.toml"]
+
+[default.extend-words]
+OT = "OT"

Project.toml

@@ -4,15 +4,21 @@ authors = ["Knut Andreas Meyer and contributors"]
 version = "0.2.4"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [compat]
-julia = "1.8"
+julia = "1.11"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2"
+
+[sources]
+MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"}
+#MaterialModelsTesting = {path = "/Users/knutan/.julia/dev/MaterialModelsTesting"}
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test", "FiniteDiff", "MaterialModelsTesting"]

docs/make.jl

@@ -19,6 +19,7 @@ makedocs(;
         "Conversion" => "conversion.md",
         "Differentiation" => "differentiation.md",
         "Implementation" => "implementing.md",
+        "Devdocs" => "devdocs.md"
     ],
 )
 

docs/src/conversion.md

@@ -2,25 +2,33 @@
 CurrentModule = MaterialModelsBase
 ```
 # Conversion
-`MaterialModelsBase` defines an interface for converting parameters and state variables
+`MaterialModelsBase` defines an interface for converting parameters, state variables, and other objects
 to and from `AbstractVector`s. This is useful when doing parameter identification, or 
 when interfacing with different languages that require information to be passed as arrays.
 
 These function can be divided into those providing information about the material and state variables, and those doing the actual conversion.
 
 ## Information functions
+The following should be defined for new materials,
 ```@docs
 get_tensorbase
+get_vector_length
+get_vector_eltype
+```
+
+Whereas the following already have default implementations that should work provided that the above are implemented.
+```@docs
 get_num_statevars
-get_statevar_eltype
-get_num_params
-get_params_eltype
 get_num_tensorcomponents
 ```
 
 ## Conversion functions
+The following should be defined for new materials, state variables, and alike
 ```@docs
 tovector!
 fromvector
+```
+whereas the following already have default implementations that should work provided that the above functions are implemented,
+```@docs
 tovector
 ```

docs/src/devdocs.md

@@ -0,0 +1,8 @@
+# Internal documentation
+The following documentation returns to internals of the package and may change
+at any time.
+
+```@docs
+MaterialModelsBase.stress_controlled_indices
+MaterialModelsBase.strain_controlled_indices
+```

src/MaterialModelsBase.jl

@@ -1,5 +1,6 @@
 module MaterialModelsBase
 using Tensors, StaticArrays
+using ForwardDiff: ForwardDiff
 import Base: @kwdef
 
 # General interface for <:AbstractMaterial
@@ -24,8 +25,9 @@ export update_stress_state!                                 # For nonzero stress
 # For parameter identification and differentiation of materials
 export tovector, tovector!, fromvector                      # Convert to/from `AbstractVector`s
 export get_num_tensorcomponents, get_num_statevars          # Information about the specific material
-export get_num_params, get_params_eltype                    # 
-export MaterialDerivatives, differentiate_material!         # Differentiation routines
+export get_vector_length, get_vector_eltype                 # Get the length and type when converting objects to vectors
+export MaterialDerivatives, StressStateDerivatives          # Derivative collections
+export differentiate_material!                              # Differentiation routines
 export allocate_differentiation_output                      # 
 
 abstract type AbstractMaterial end
@@ -112,14 +114,16 @@ end
 Return the (default) initial `state::AbstractMaterialState` 
 of the material `m`. 
 
-Defaults to the empty `NoMaterialState()`
+Defaults to the empty `NoMaterialState{T}()`, where 
+`T = get_vector_eltype(m)` (which defaults to `Float64`).
 """
-function initial_material_state(::AbstractMaterial) 
-    return NoMaterialState()
+function initial_material_state(m::AbstractMaterial)
+    T = get_vector_eltype(m)
+    return NoMaterialState{T}()
 end
 
 abstract type AbstractMaterialState end
-struct NoMaterialState <: AbstractMaterialState end
+struct NoMaterialState{T} <: AbstractMaterialState end
 
 # Material cache 
 """
@@ -153,8 +157,8 @@ abstract type AbstractExtraOutput end
 struct NoExtraOutput <: AbstractExtraOutput end
 
 include("vector_conversion.jl")
-include("differentiation.jl")
 include("stressiterations.jl")
+include("differentiation.jl")
 include("ErrorExceptions.jl")
 
 end

src/differentiation.jl

@@ -1,19 +1,31 @@
 struct MaterialDerivatives{T}
     dσdϵ::Matrix{T}
-    dσdⁿs::Matrix{T}
+    #dσdⁿs::Matrix{T}
     dσdp::Matrix{T}
     dsdϵ::Matrix{T}
-    dsdⁿs::Matrix{T}
+    #dsdⁿs::Matrix{T}
     dsdp::Matrix{T}
 end
 
+function Base.getproperty(d::MaterialDerivatives, key::Symbol)
+    if key === :dσdⁿs || key === :dsdⁿs
+        error("You are probably assuming MaterialModelsBase v0.2 behavior for differentiation")
+    else
+        @inline getfield(d, key)
+    end
+end
+
 """
     MaterialDerivatives(m::AbstractMaterial)
 
 A struct that saves all derivative information using a `Matrix{T}` for each derivative,
-where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, 
-`get_num_statevars`, and `get_num_params`. The values should be updated in `differentiate_material!`
-by direct access of the fields, where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
+where `T=get_vector_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, 
+`get_num_statevars`, and `get_vector_length`. `m` must support `tovector` and `fromvector`, while 
+the output of `initial_material_state` must support `tovector`, and in addition the element type
+of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`.
+
+The values should be updated in `differentiate_material!` by direct access of the fields, 
+where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
 and old state variables, and `p` the material parameter vector.
 
 * `dσdϵ`
@@ -25,17 +37,16 @@ and old state variables, and `p` the material parameter vector.
 
 """
 function MaterialDerivatives(m::AbstractMaterial)
-    T = get_params_eltype(m)
+    T = get_vector_eltype(m)
     n_tensor = get_num_tensorcomponents(m)
     n_state = get_num_statevars(m)
-    n_params = get_num_params(m)
+    n_params = get_vector_length(m)
+    dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m))
     return MaterialDerivatives(
         zeros(T, n_tensor, n_tensor),  # dσdϵ
-        zeros(T, n_tensor, n_state),   # dσdⁿs
         zeros(T, n_tensor, n_params),  # dσdp
         zeros(T, n_state, n_tensor),   # dsdϵ
-        zeros(T, n_state, n_state),    # dsdⁿs
-        zeros(T, n_state, n_params)    # dsdp
+        dsdp
         )
 end
 
@@ -65,5 +76,108 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput()
 
 Calculate the derivatives and save them in `diff`, see
 [`MaterialDerivatives`](@ref) for a description of the fields in `diff`.
+
+    differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...)
+
+For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`,
+this method will work automatically by
+1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted)
+2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})`
+3) Updating `ssd` according to the constraints imposed by the `stress_state`.
+
+For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly
+to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the 
+(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. 
+The exact interface of the `StressStateDerivatives` datastructure is still not fixed, and may change in the future.
 """
-function differentiate_material! end
+function differentiate_material! end
+
+struct StressStateDerivatives{T}
+    mderiv::MaterialDerivatives{T}
+    dϵdp::Matrix{T}
+    dσdp::Matrix{T}
+    ϵindex::SMatrix{3, 3, Int} # To allow indexing by (i, j) into only
+    σindex::SMatrix{3, 3, Int} # saved values to avoid storing unused rows.
+    # TODO: Reduce the dimensions, for now all entries (even those that are zero) are stored.
+end
+
+function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial)
+    mderiv = MaterialDerivatives(m)
+    np = get_vector_length(m)
+    nt = get_num_tensorcomponents(m) # Should be changed to only save non-controlled entries
+    dϵdp = zeros(nt, np)
+    dσdp = zeros(nt, np)
+    # Should be changed to only save non-controlled entries
+    vo = Tensors.DEFAULT_VOIGT_ORDER[3]
+    if nt == 6
+        index = SMatrix{3, 3, Int}(min(vo[i, j], vo[j, i]) for i in 1:3, j in 1:3)
+    else
+        index = SMatrix{3, 3, Int}(vo)
+    end
+    return StressStateDerivatives(mderiv, dϵdp, dσdp, index, index)
+end
+
+function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where {N}
+    σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, ϵ, args...)
+    differentiate_material!(ssd.mderiv, m, ϵ_full, args..., dσdϵ_full)
+
+    if isa(stress_state, NoIterationState)
+        copy!(ssd.dσdp, ssd.mderiv.dσdp)
+        fill!(ssd.dϵdp, 0)
+    else
+        sc = stress_controlled_indices(stress_state, ϵ)
+        ec = strain_controlled_indices(stress_state, ϵ)
+        dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations
+        ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :]
+        ssd.dσdp[ec, :] .= ssd.mderiv.dσdp[ec, :] .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :]
+        ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :]
+    end
+    return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full
+end
+
+"""
+    stress_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are stress-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function stress_controlled_indices end
+
+"""
+    strain_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are strain-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function strain_controlled_indices end
+
+# NoIterationState
+stress_controlled_indices(::NoIterationState, ::AbstractTensor) = SVector{0,Int}()
+strain_controlled_indices(::NoIterationState, ::SymmetricTensor) = @SVector([1, 2, 3, 4, 5, 6])
+strain_controlled_indices(::NoIterationState, ::Tensor) = @SVector([1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+# UniaxialStress
+stress_controlled_indices(::UniaxialStress, ::SymmetricTensor) = @SVector([2, 3, 4, 5, 6])
+stress_controlled_indices(::UniaxialStress, ::Tensor) = @SVector([2, 3, 4, 5, 6, 7, 8, 9])
+strain_controlled_indices(::UniaxialStress, ::AbstractTensor) = @SVector([1])
+
+# UniaxialNormalStress
+stress_controlled_indices(::UniaxialNormalStress, ::AbstractTensor) = @SVector([2,3])
+strain_controlled_indices(::UniaxialNormalStress, ::SymmetricTensor) = @SVector([1, 4, 5, 6])
+strain_controlled_indices(::UniaxialNormalStress, ::Tensor) = @SVector([1, 4, 5, 6, 7, 8, 9])
+
+# PlaneStress 12 -> 6, 21 -> 9
+stress_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([3, 4, 5])
+stress_controlled_indices(::PlaneStress, ::Tensor) = @SVector([3, 4, 5, 7, 8])
+strain_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([1, 2, 6])
+strain_controlled_indices(::PlaneStress, ::Tensor) = @SVector([1, 2, 6, 9])
+
+# GeneralStressState
+stress_controlled_indices(ss::GeneralStressState{Nσ}, ::AbstractTensor) where Nσ = controlled_indices_from_tensor(ss.σ_ctrl, true, Val(Nσ))
+function strain_controlled_indices(ss::GeneralStressState{Nσ,TT}, ::AbstractTensor) where {Nσ,TT}
+    N = Tensors.n_components(Tensors.get_base(TT)) - Nσ
+    return controlled_indices_from_tensor(ss.σ_ctrl, false, Val(N))
+end
+function controlled_indices_from_tensor(ctrl::AbstractTensor, return_if, ::Val{N}) where N
+    return SVector{N}(i for (i, v) in pairs(tovoigt(SVector, ctrl)) if v == return_if)
+end