11struct MaterialDerivatives{T}
22 dσdϵ:: Matrix{T}
3- dσdⁿs:: Matrix{T}
3+ # dσdⁿs::Matrix{T}
44 dσdp:: Matrix{T}
55 dsdϵ:: Matrix{T}
6- dsdⁿs:: Matrix{T}
6+ # dsdⁿs::Matrix{T}
77 dsdp:: Matrix{T}
88end
99
10+ function Base. getproperty (d:: MaterialDerivatives , key:: Symbol )
11+ if key === :dσdⁿs || key === :dsdⁿs
12+ error (" You are probably assuming MaterialModelsBase v0.2 behavior for differentiation"
13+ else
14+ @inline getfield (d, key)
15+ end
16+ end
17+
1018"""
1119 MaterialDerivatives(m::AbstractMaterial)
1220
1321A struct that saves all derivative information using a `Matrix{T}` for each derivative,
14- where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`,
15- `get_num_statevars`, and `get_num_params`. The values should be updated in `differentiate_material!`
16- by direct access of the fields, where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current
22+ where `T=get_vector_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`,
23+ `get_num_statevars`, and `get_vector_length`. `m` must support `tovector` and `fromvector`, while
24+ the output of `initial_material_state` must support `tovector`, and in addition the element type
25+ of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`.
26+
27+ The values should be updated in `differentiate_material!` by direct access of the fields,
28+ where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current
1729and old state variables, and `p` the material parameter vector.
1830
1931* `dσdϵ`
@@ -25,17 +37,16 @@ and old state variables, and `p` the material parameter vector.
2537
2638"""
2739function MaterialDerivatives (m:: AbstractMaterial )
28- T = get_params_eltype (m)
40+ T = get_vector_eltype (m)
2941 n_tensor = get_num_tensorcomponents (m)
3042 n_state = get_num_statevars (m)
31- n_params = get_num_params (m)
43+ n_params = get_vector_length (m)
44+ dsdp = ForwardDiff. jacobian (p -> tovector (initial_material_state (fromvector (p, m))), tovector (m))
3245 return MaterialDerivatives (
3346 zeros (T, n_tensor, n_tensor), # dσdϵ
34- zeros (T, n_tensor, n_state), # dσdⁿs
3547 zeros (T, n_tensor, n_params), # dσdp
3648 zeros (T, n_state, n_tensor), # dsdϵ
37- zeros (T, n_state, n_state), # dsdⁿs
38- zeros (T, n_state, n_params) # dsdp
49+ dsdp
3950 )
4051end
4152
@@ -65,5 +76,108 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput()
6576
6677Calculate the derivatives and save them in `diff`, see
6778[`MaterialDerivatives`](@ref) for a description of the fields in `diff`.
79+
80+ differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...)
81+
82+ For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`,
83+ this method will work automatically by
84+ 1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted)
85+ 2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})`
86+ 3) Updating `ssd` according to the constraints imposed by the `stress_state`.
87+
88+ For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly
89+ to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the
90+ (total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated.
91+ The exact interface of the `StressStateDerivatives` datastructure is still not fixed, and may change in the future.
6892"""
69- function differentiate_material! end
93+ function differentiate_material! end
94+
95+ struct StressStateDerivatives{T}
96+ mderiv:: MaterialDerivatives{T}
97+ dϵdp:: Matrix{T}
98+ dσdp:: Matrix{T}
99+ ϵindex:: SMatrix{3, 3, Int} # To allow indexing by (i, j) into only
100+ σindex:: SMatrix{3, 3, Int} # saved values to avoid storing unused rows.
101+ # TODO : Reduce the dimensions, for now all entries (even those that are zero) are stored.
102+ end
103+
104+ function StressStateDerivatives (:: AbstractStressState , m:: AbstractMaterial )
105+ mderiv = MaterialDerivatives (m)
106+ np = get_vector_length (m)
107+ nt = get_num_tensorcomponents (m) # Should be changed to only save non-controlled entries
108+ dϵdp = zeros (nt, np)
109+ dσdp = zeros (nt, np)
110+ # Should be changed to only save non-controlled entries
111+ vo = Tensors. DEFAULT_VOIGT_ORDER[3 ]
112+ if nt == 6
113+ index = SMatrix {3, 3, Int} (min (vo[i, j], vo[j, i]) for i in 1 : 3 , j in 1 : 3 )
114+ else
115+ index = SMatrix {3, 3, Int} (vo)
116+ end
117+ return StressStateDerivatives (mderiv, dϵdp, dσdp, index, index)
118+ end
119+
120+ function differentiate_material! (ssd:: StressStateDerivatives , stress_state:: AbstractStressState , m:: AbstractMaterial , ϵ:: AbstractTensor , args:: Vararg{Any,N} ) where {N}
121+ σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response (stress_state, m, ϵ, args... )
122+ differentiate_material! (ssd. mderiv, m, ϵ_full, args... , dσdϵ_full)
123+
124+ if isa (stress_state, NoIterationState)
125+ copy! (ssd. dσdp, ssd. mderiv. dσdp)
126+ fill! (ssd. dϵdp, 0 )
127+ else
128+ sc = stress_controlled_indices (stress_state, ϵ)
129+ ec = strain_controlled_indices (stress_state, ϵ)
130+ dσᶠdϵᶠ_inv = inv (get_unknowns (stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations
131+ ssd. dϵdp[sc, :] .= .- dσᶠdϵᶠ_inv * ssd. mderiv. dσdp[sc, :]
132+ ssd. dσdp[ec, :] .= ssd. mderiv. dσdp[ec, :] .+ ssd. mderiv. dσdϵ[ec, sc] * ssd. dϵdp[sc, :]
133+ ssd. mderiv. dsdp .+ = ssd. mderiv. dsdϵ[:, sc] * ssd. dϵdp[sc, :]
134+ end
135+ return reduce_tensordim (stress_state, σ_full), reduce_stiffness (stress_state, dσdϵ_full), state, ϵ_full
136+ end
137+
138+ """
139+ stress_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
140+
141+ Get the `N` indices that are stress-controlled in `stress_state`. The tensor input is used to
142+ determine if a symmetric or full tensor is used.
143+ """
144+ function stress_controlled_indices end
145+
146+ """
147+ strain_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
148+
149+ Get the `N` indices that are strain-controlled in `stress_state`. The tensor input is used to
150+ determine if a symmetric or full tensor is used.
151+ """
152+ function strain_controlled_indices end
153+
154+ # NoIterationState
155+ stress_controlled_indices (:: NoIterationState , :: AbstractTensor ) = SVector {0,Int} ()
156+ strain_controlled_indices (:: NoIterationState , :: SymmetricTensor ) = @SVector ([1 , 2 , 3 , 4 , 5 , 6 ])
157+ strain_controlled_indices (:: NoIterationState , :: Tensor ) = @SVector ([1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 ])
158+
159+ # UniaxialStress
160+ stress_controlled_indices (:: UniaxialStress , :: SymmetricTensor ) = @SVector ([2 , 3 , 4 , 5 , 6 ])
161+ stress_controlled_indices (:: UniaxialStress , :: Tensor ) = @SVector ([2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 ])
162+ strain_controlled_indices (:: UniaxialStress , :: AbstractTensor ) = @SVector ([1 ])
163+
164+ # UniaxialNormalStress
165+ stress_controlled_indices (:: UniaxialNormalStress , :: AbstractTensor ) = @SVector ([2 ,3 ])
166+ strain_controlled_indices (:: UniaxialNormalStress , :: SymmetricTensor ) = @SVector ([1 , 4 , 5 , 6 ])
167+ strain_controlled_indices (:: UniaxialNormalStress , :: Tensor ) = @SVector ([1 , 4 , 5 , 6 , 7 , 8 , 9 ])
168+
169+ # PlaneStress 12 -> 6, 21 -> 9
170+ stress_controlled_indices (:: PlaneStress , :: SymmetricTensor ) = @SVector ([3 , 4 , 5 ])
171+ stress_controlled_indices (:: PlaneStress , :: Tensor ) = @SVector ([3 , 4 , 5 , 7 , 8 ])
172+ strain_controlled_indices (:: PlaneStress , :: SymmetricTensor ) = @SVector ([1 , 2 , 6 ])
173+ strain_controlled_indices (:: PlaneStress , :: Tensor ) = @SVector ([1 , 2 , 6 , 9 ])
174+
175+ # GeneralStressState
176+ stress_controlled_indices (ss:: GeneralStressState{Nσ} , :: AbstractTensor ) where Nσ = controlled_indices_from_tensor (ss. σ_ctrl, true , Val (Nσ))
177+ function strain_controlled_indices (ss:: GeneralStressState{Nσ,TT} , :: AbstractTensor ) where {Nσ,TT}
178+ N = Tensors. n_components (Tensors. get_base (TT)) - Nσ
179+ return controlled_indices_from_tensor (ss. σ_ctrl, false , Val (N))
180+ end
181+ function controlled_indices_from_tensor (ctrl:: AbstractTensor , return_if, :: Val{N} ) where N
182+ return SVector {N} (i for (i, v) in pairs (tovoigt (SVector, ctrl)) if v == return_if)
183+ end
0 commit comments