"""Drucker-Prager elasto-plastic constitutive model with isotropic hardening.
This module includes the implementation of the Drucker-Prager constitutive
model with non-associative plastic flor rule and associative isotropic strain
hardening.
Classes
-------
DruckerPrager
Drucker-Prager constitutive model with isotropic strain hardening.
"""
#
# Modules
# =============================================================================
# Standard
import math
import copy
# Third-party
import torch
# Local
from simulators.fetorch.material.models.interface import ConstitutiveModel
from simulators.fetorch.material.models.standard.elastic import Elastic
from simulators.fetorch.math.matrixops import get_problem_type_parameters, \
vget_tensor_mf, vget_tensor_from_mf, vget_state_3Dmf_from_2Dmf, \
vget_state_2Dmf_from_3Dmf
from simulators.fetorch.math.tensorops import get_id_operators, dyad22_1
from utilities.type_conversion import convert_dict_to_tensor, \
convert_tensor_to_float64, convert_dict_to_float64, \
convert_dict_to_float32, convert_tensor_to_float32
#
# Authorship & Credits
# =============================================================================
__author__ = 'Bernardo Ferreira (bernardo_ferreira@brown.edu)'
__credits__ = ['Bernardo Ferreira', ]
__status__ = 'Stable'
# =============================================================================
#
# =============================================================================
[docs]class DruckerPrager(ConstitutiveModel):
"""Drucker-Prager constitutive model with isotropic strain hardening.
Attributes
----------
_name : str
Constitutive model name.
_strain_type : {'infinitesimal', 'finite', 'finite-kinext'}
Material constitutive model strain formulation: infinitesimal strain
formulation ('infinitesimal'), finite strain formulation ('finite') or
finite strain formulation through kinematic extension
('finite-kinext').
_model_parameters : dict
Material constitutive model parameters.
_n_dim : int
Problem number of spatial dimensions.
_comp_order_sym : list
Strain/Stress components symmetric order.
_comp_order_nsym : list
Strain/Stress components nonsymmetric order.
_is_su_float64 : bool
If True, then state update is locally computed in floating-point
double precision. If False, then default floating-point precision
is assumed.
_device_type : {'cpu', 'cuda'}
Type of device on which torch.Tensor is allocated.
_device : torch.device
Device on which torch.Tensor is allocated.
Methods
-------
get_required_model_parameters()
Get required material constitutive model parameters.
state_init(self)
Get initialized material constitutive model state variables.
state_update(self, inc_strain, state_variables_old)
Perform material constitutive model state update.
"""
[docs] def __init__(self, strain_formulation, problem_type, model_parameters,
is_su_float64=True, device_type='cpu'):
"""Constitutive model constructor.
Parameters
----------
strain_formulation: {'infinitesimal', 'finite'}
Problem strain formulation.
problem_type : int
Problem type: 2D plane strain (1), 2D plane stress (2),
2D axisymmetric (3) and 3D (4).
model_parameters : dict
Material constitutive model parameters.
device_type : {'cpu', 'cuda'}, default='cpu'
Type of device on which torch.Tensor is allocated.
"""
# Set material constitutive model name
self._name = 'drucker_prager'
# Set constitutive model strain formulation
self._strain_type = 'finite-kinext'
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set initialization parameters
self._strain_formulation = strain_formulation
self._problem_type = problem_type
self._model_parameters = convert_dict_to_tensor(model_parameters,
is_inplace=True)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set state update floating-point precision
self._is_su_float64 = is_su_float64
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set device
self.set_device(device_type)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Get problem type parameters
self._n_dim, self._comp_order_sym, self._comp_order_nsym = \
get_problem_type_parameters(problem_type)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Get elastic symmetry
elastic_symmetry = model_parameters['elastic_symmetry']
# Check finite strains formulation
if self._strain_formulation == 'finite' and \
elastic_symmetry != 'isotropic':
raise RuntimeError('The Drucker-Prager constitutive model is only '
'available under finite strains for the '
'elastic isotropic case.')
# Compute technical constants of elasticity
if elastic_symmetry == 'isotropic':
# Compute technical constants of elasticity
technical_constants = Elastic.get_technical_from_elastic_moduli(
elastic_symmetry, model_parameters)
# Assemble technical constants of elasticity
self._model_parameters.update(technical_constants)
else:
raise RuntimeError('The Drucker-Prager constitutive model is '
'currently only available for the elastic '
'isotropic case.')
# -------------------------------------------------------------------------
[docs] @staticmethod
def get_required_model_parameters():
"""Get required material constitutive model parameters.
Model parameters:
- 'elastic_symmetry' : Elastic symmetry (str, {'isotropic',
'transverse_isotropic', 'orthotropic', 'monoclinic', 'triclinic'});
- 'elastic_moduli' : Elastic moduli (dict, {'Eijkl': float});
- 'euler_angles' : Euler angles (degrees) sorted according with Bunge
convention (tuple[float]).
- 'yield_cohesion_parameter' : Yield surface cohesion parameter
- 'yield_pressure_parameter' : Yield surface pressure parameter
- 'flow_pressure_parameter' : Plastic flow rule pressure parameter
- 'hardening_law' : Isotropic hardening law (function)
- 'hardening_parameters' : Isotropic hardening law parameters (dict)
Returns
-------
model_parameters_names : tuple[str]
Material constitutive model parameters names (str).
"""
# Set material properties names
model_parameters_names = ('elastic_symmetry', 'elastic_moduli',
'euler_angles',
'yield_cohesion_parameter',
'yield_pressure_parameter',
'flow_pressure_parameter',
'hardening_law', 'hardening_parameters')
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
return model_parameters_names
# -------------------------------------------------------------------------
[docs] def state_init(self):
"""Get initialized material constitutive model state variables.
Constitutive model state variables:
* ``e_strain_mf``
* *Infinitesimal strains*: Elastic infinitesimal strain tensor
(matricial form).
* *Finite strains*: Elastic spatial logarithmic strain tensor
(matricial form).
* *Symbol*: :math:`\\boldsymbol{\\varepsilon^{e}}` /
:math:`\\boldsymbol{\\varepsilon^{e}}`
* ``acc_p_strain``
* Accumulated plastic strain.
* *Symbol*: :math:`\\bar{\\varepsilon}^{p}`
* ``strain_mf``
* *Infinitesimal strains*: Infinitesimal strain tensor
(matricial form).
* *Finite strains*: Spatial logarithmic strain tensor
(matricial form).
* *Symbol*: :math:`\\boldsymbol{\\varepsilon}` /
:math:`\\boldsymbol{\\varepsilon}`
* ``stress_mf``
* *Infinitesimal strains*: Cauchy stress tensor (matricial form).
* *Finite strains*: Kirchhoff stress tensor (matricial form) within
:py:meth:`state_update`, first Piola-Kirchhoff stress tensor
(matricial form) otherwise.
* *Symbol*: :math:`\\boldsymbol{\\sigma}` /
(:math:`\\boldsymbol{\\tau}`, :math:`\\boldsymbol{P}`)
* ``is_plastic``
* Plastic step flag.
* ``is_su_fail``
* State update failure flag.
* ``is_apex_return``
* Return-mapping to apex flag.
----
Returns
-------
state_variables_init : dict
Initialized material constitutive model state variables.
"""
# Initialize constitutive model state variables
state_variables_init = dict()
# Initialize strain tensors
state_variables_init['e_strain_mf'] = vget_tensor_mf(
torch.zeros((self._n_dim, self._n_dim), device=self._device),
self._n_dim, self._comp_order_sym)
state_variables_init['strain_mf'] = \
state_variables_init['e_strain_mf'].clone()
# Initialize stress tensors
if self._strain_formulation == 'infinitesimal':
# Cauchy stress tensor (symmetric)
state_variables_init['stress_mf'] = vget_tensor_mf(
torch.zeros((self._n_dim, self._n_dim), device=self._device),
self._n_dim, self._comp_order_sym)
else:
# First Piola-Kirchhoff stress tensor (nonsymmetric)
state_variables_init['stress_mf'] = vget_tensor_mf(
torch.zeros((self._n_dim, self._n_dim), device=self._device),
self._n_dim, self._comp_order_nsym)
# Initialize internal variables
state_variables_init['acc_p_strain'] = \
torch.tensor(0.0, device=self._device)
# Initialize state flags
state_variables_init['is_plast'] = \
torch.tensor(False, device=self._device)
state_variables_init['is_su_fail'] = \
torch.tensor(False, device=self._device)
state_variables_init['is_apex_return'] = \
torch.tensor(False, device=self._device)
# Set additional out-of-plane strain and stress components
if self._problem_type == 1:
state_variables_init['e_strain_33'] = \
torch.tensor(0.0, device=self._device)
state_variables_init['stress_33'] = \
torch.tensor(0.0, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Return
return state_variables_init
# -------------------------------------------------------------------------
[docs] def state_update(self, inc_strain, state_variables_old):
"""Perform material constitutive model state update.
Parameters
----------
inc_strain : torch.Tensor(2d)
Incremental strain second-order tensor.
state_variables_old : dict
Last converged constitutive model material state variables.
Returns
-------
state_variables : dict
Material constitutive model state variables.
consistent_tangent_mf : torch.Tensor(2d)
Material constitutive model consistent tangent modulus stored in
matricial form.
"""
# Get model parameters
model_parameters = self._model_parameters
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize floating-point precision conversion flag
is_precision_conversion = False
# Handle state update floating-point precision
if torch.get_default_dtype() == torch.float32 and self._is_su_float64:
# Set floating-point precision conversion flag
is_precision_conversion = True
# Set default floating-point precision
torch.set_default_dtype(torch.float64)
# Perform floating-point precision conversion
model_parameters = convert_dict_to_float64(model_parameters,
is_inplace=False)
inc_strain = convert_tensor_to_float64(inc_strain)
state_variables_old = convert_dict_to_float64(state_variables_old,
is_inplace=False)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set state update convergence tolerance
su_conv_tol = 1e-6
# Set state update maximum number of iterations
su_max_n_iterations = 20
# Set minimum threshold to handle values close or equal to zero
small = 1e-8
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Build incremental strain tensor matricial form
inc_strain_mf = vget_tensor_mf(inc_strain, self._n_dim,
self._comp_order_sym,
device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Get material properties
E = model_parameters['E']
v = model_parameters['v']
xi = model_parameters['yield_cohesion_parameter']
etay = model_parameters['yield_pressure_parameter']
etaf = model_parameters['flow_pressure_parameter']
# Get material isotropic strain hardening law
hardening_law = model_parameters['hardening_law']
hardening_parameters = model_parameters['hardening_parameters']
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute bulk and shear modulus
K = E/(3.0*(1.0 - 2.0*v))
G = E/(2.0*(1.0 + v))
# Compute Lamé parameters
lam = (E*v)/((1.0 + v)*(1.0 - 2.0*v))
miu = E/(2.0*(1.0 + v))
# Compute material parameters
alpha = xi/etaf
beta = xi/etay
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Get last increment converged state variables
e_strain_old_mf = state_variables_old['e_strain_mf']
p_strain_old_mf = state_variables_old['strain_mf'] - e_strain_old_mf
acc_p_strain_old = state_variables_old['acc_p_strain']
if self._problem_type == 1:
e_strain_33_old = state_variables_old['e_strain_33']
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize state update failure flag
is_su_fail = torch.tensor(False, device=self._device)
# Initialize plastic step flag
is_plast = torch.tensor(False, device=self._device)
# Initialize return-mapping to apex flag
is_apex_return = torch.tensor(False, device=self._device)
#
# 2D > 3D conversion
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# When the problem type corresponds to a 2D analysis, perform the state
# update and consistent tangent computation as in the 3D case,
# considering the appropriate out-of-plain strain and stress components
if self._problem_type == 4:
n_dim = self._n_dim
comp_order_sym = self._comp_order_sym
else:
# Set 3D problem parameters
n_dim, comp_order_sym, _ = get_problem_type_parameters(4)
# Build strain tensors (matricial form) by including the
# appropriate out-of-plain components
inc_strain_mf = vget_state_3Dmf_from_2Dmf(
inc_strain_mf, comp_33=0.0, device=self._device)
e_strain_old_mf = vget_state_3Dmf_from_2Dmf(
e_strain_old_mf, e_strain_33_old, device=self._device)
#
# State update
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set required second-order and fourth-order tensors
soid, _, _, fosym, fodiagtrace, _, fodevprojsym = \
get_id_operators(n_dim, device=self._device)
soid_mf = vget_tensor_mf(soid, n_dim, comp_order_sym,
device=self._device)
fodevprojsym_mf = vget_tensor_mf(fodevprojsym, n_dim, comp_order_sym,
device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute elastic trial strain
e_trial_strain_mf = e_strain_old_mf + inc_strain_mf
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute elastic consistent tangent modulus according to problem type
# and store it in matricial form
if self._problem_type in [1, 4]:
e_consistent_tangent = lam*fodiagtrace + 2.0*miu*fosym
e_consistent_tangent_mf = vget_tensor_mf(e_consistent_tangent,
n_dim, comp_order_sym,
device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute trial stress
trial_stress_mf = torch.matmul(e_consistent_tangent_mf,
e_trial_strain_mf)
# Compute trial pressure
trial_pressure = \
(1.0/3.0)*torch.trace(vget_tensor_from_mf(trial_stress_mf,
n_dim, comp_order_sym,
device=self._device))
# Compute deviatoric trial stress
dev_trial_stress_mf = torch.matmul(fodevprojsym_mf, trial_stress_mf)
# Compute second invariant of deviatoric trial stress
j2_dev_trial_stress = 0.5*torch.norm(dev_trial_stress_mf)**2
# Compute trial accumulated plastic strain
acc_p_trial_strain = acc_p_strain_old
# Compute trial cohesion
trial_cohesion, _ = \
hardening_law(hardening_parameters, acc_p_trial_strain)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute yield function
yield_function = (torch.sqrt(j2_dev_trial_stress) + etay*trial_pressure
- xi*trial_cohesion)
# Check plastic consistency
if torch.abs(trial_cohesion) < small:
plastic_consistency = yield_function
else:
plastic_consistency = yield_function/torch.abs(trial_cohesion)
# If the trial stress state lies inside the Druger-Prager yield
# function, then the state update is purely elastic and coincident with
# the elastic trial state. Otherwise, the state update is elastoplastic
# and the return-mapping system of nonlinear equations must be solved
# in order to update the state variables
if plastic_consistency <= su_conv_tol:
# Update elastic strain
e_strain_mf = e_trial_strain_mf
# Update stress
stress_mf = trial_stress_mf
# Update accumulated plastic strain
acc_p_strain = acc_p_strain_old
else:
# Set plastic step flag
is_plast = torch.tensor(True, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set incremental plastic multiplier initial iterative guess
inc_p_mult = torch.tensor(0.0, device=self._device)
# Compute initial hardening modulus
cohesion, H = hardening_law(hardening_parameters,
acc_p_strain_old + xi*inc_p_mult)
# Initialize Newton-Raphson iteration counter
nr_iter = 0
# Compute return-mapping residual (cone surface)
residual = (torch.sqrt(j2_dev_trial_stress) + etay*trial_pressure
- xi*cohesion)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Start Newton-Raphson iterative loop
while True:
# Compute return-mapping Jacobian (scalar)
jacobian = -G - K*etaf*etay - (xi**2)*H
# Solve return-mapping linearized equation
d_iter = -residual/jacobian
# Update incremental plastic multiplier
inc_p_mult = inc_p_mult + d_iter
# Compute current cohesion and hardening modulus
cohesion, H = hardening_law(hardening_parameters,
acc_p_strain_old + xi*inc_p_mult)
# Compute return-mapping residual (cone surface)
residual = (torch.sqrt(j2_dev_trial_stress) - G*inc_p_mult
+ etay*(trial_pressure - K*etaf*inc_p_mult)
- xi*cohesion)
# Compute residual convergence norm
if torch.abs(cohesion) < small:
conv_norm_res = torch.abs(residual)
else:
conv_norm_res = torch.abs(residual/cohesion)
# Check Newton-Raphson iterative procedure convergence
is_converged = (conv_norm_res < su_conv_tol
and nr_iter > 0)
# Control Newton-Raphson iteration loop flow
if is_converged:
# Leave Newton-Raphson iterative loop (converged solution)
break
elif nr_iter == su_max_n_iterations:
# Update state update failure flag
is_su_fail = torch.tensor(True, device=self._device)
# Leave Newton-Raphson iterative loop (failed solution)
break
else:
# Increment iteration counter
nr_iter = nr_iter + 1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute pressure
pressure = trial_pressure - K*etaf*inc_p_mult
# Compute deviatoric trial stress second invariant factor
if torch.sqrt(j2_dev_trial_stress) > small:
dev_stress_factor = \
(1.0 - ((G*inc_p_mult)/torch.sqrt(j2_dev_trial_stress)))
else:
dev_stress_factor = torch.zeros(1, device=self._device)
# Compute deviatoric stress
dev_stress_mf = dev_stress_factor*dev_trial_stress_mf
# Update stress
stress_mf = pressure*soid_mf + dev_stress_mf
# Update accumulated plastic strain
acc_p_strain = acc_p_strain_old + xi*inc_p_mult
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Check return-mapping validity
if not is_su_fail:
is_cone_surface = \
torch.sqrt(j2_dev_trial_stress) - G*inc_p_mult >= 0
else:
# If convergence of return-mapping to cone surface failed, then
# avoid return-mapping to cone apex
is_cone_surface = torch.tensor(True, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute return-mapping to cone apex
if not is_su_fail and not is_cone_surface:
# Set return-mapping to apex flag
is_apex_return = torch.tensor(True, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set incremental plastic volumetric strain initial iterative
# guess
inc_vol_p_strain = torch.tensor(0.0, device=self._device)
# Compute initial hardening modulus
cohesion, H = hardening_law(
hardening_parameters,
acc_p_strain_old + alpha*inc_vol_p_strain)
# Initialize Newton-Raphson iteration counter
nr_iter = 0
# Compute return-mapping residual (cone apex)
residual = cohesion*beta - trial_pressure
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Start Newton-Raphson iterative loop
while True:
# Compute return-mapping Jacobian (scalar)
jacobian = alpha*beta*H + K
# Solve return-mapping linearized equation
d_iter = -residual/jacobian
# Update incremental plastic volumetric strain
inc_vol_p_strain = inc_vol_p_strain + d_iter
# Compute current cohesion and hardening modulus
cohesion, H = hardening_law(
hardening_parameters,
acc_p_strain_old + alpha*inc_vol_p_strain)
# Compute return-mapping residual (cone apex)
residual = \
cohesion*beta - (trial_pressure - K*inc_vol_p_strain)
# Compute residual convergence norm
if torch.abs(cohesion) < small:
conv_norm_res = torch.abs(residual)
else:
conv_norm_res = torch.abs(residual/cohesion)
# Check Newton-Raphson iterative procedure convergence
is_converged = (conv_norm_res < su_conv_tol
and nr_iter > 0)
# Control Newton-Raphson iteration loop flow
if is_converged:
# Leave Newton-Raphson iterative loop (converged
# solution)
break
elif nr_iter == su_max_n_iterations:
# Update state update failure flag
is_su_fail = torch.tensor(True, device=self._device)
# Leave Newton-Raphson iterative loop (failed solution)
break
else:
# Increment iteration counter
nr_iter = nr_iter + 1
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute pressure
pressure = trial_pressure - K*inc_vol_p_strain
# Compute deviatoric stress
dev_stress_mf = torch.zeros_like(dev_trial_stress_mf,
device=self._device)
# Update stress
stress_mf = pressure*soid_mf
# Update accumulated plastic strain
acc_p_strain = acc_p_strain_old + alpha*inc_vol_p_strain
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Update elastic strain
e_strain_mf = \
(1.0/(3.0*K))*pressure*soid_mf + (1.0/(2.0*G))*dev_stress_mf
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set state variables to NaN if state update fails
if is_su_fail:
# Set elastic strain to NaN if state update fails
e_strain_mf = torch.full(e_strain_mf.shape, torch.nan,
device=self._device)
# Set stress to NaN if state update fails
stress_mf = torch.full(stress_mf.shape, torch.nan,
device=self._device)
# Set accumulated plastic strain to NaN if state update fails
acc_p_strain = torch.tensor(torch.nan, device=self._device)
# Set incremental plastic multiplier to NaN if state update
# fails
inc_p_mult = torch.tensor(torch.nan, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Get the out-of-plane strain and stress components
if self._problem_type == 1:
e_strain_33 = e_strain_mf[comp_order_sym.index('33')]
stress_33 = stress_mf[comp_order_sym.index('33')]
#
# Consistent tangent modulus
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute consistent tangent modulus if state update converged,
# otherwise set to NaN
if is_su_fail:
# Set consistent tangent modulus to NaN
consistent_tangent_mf = torch.full(e_consistent_tangent_mf.shape,
torch.nan, device=self._device)
else:
# If the state update was purely elastic, then the consistent
# tangent modulus is the elastic consistent tangent modulus.
# Otherwise, compute the elastoplastic consistent tangent modulus
if is_plast:
if is_apex_return:
# Compute elastoplastic consistent tangent modulus
# (cone surface)
consistent_tangent = \
K*(1.0 - (K/(K + alpha*beta*H)))*dyad22_1(soid, soid)
else:
# Compute deviatoric elastic trial strain
dev_trial_e_strain = vget_tensor_from_mf(
torch.matmul(fodevprojsym_mf, e_trial_strain_mf),
n_dim, comp_order_sym, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute deviatoric elastic trial strain norm divison
# factor
if torch.norm(e_trial_strain_mf) > small:
norm_div_factor = 1.0/torch.norm(e_trial_strain_mf)
else:
norm_div_factor = torch.zeros(1, device=self._device)
# Compute deviatoric elastic strain unit vector
trial_unit = norm_div_factor*dev_trial_e_strain
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute common scalar terms
s1 = norm_div_factor*(inc_p_mult/math.sqrt(2))
s2 = 1.0/(G + K*etay*etaf + H*xi**2)
# Compute elastoplastic consistent tangent modulus
# (cone apex)
consistent_tangent = 2.0*G*(1.0 - s1)*fodevprojsym \
+ 2.0*G*(s1 - G*s2)*dyad22_1(trial_unit, trial_unit) \
- math.sqrt(2)*G*s2*K*(etay*dyad22_1(trial_unit, soid)
+ etaf*dyad22_1(soid, trial_unit)) \
+ K*(1.0 - K*etay*etaf*s2)*dyad22_1(soid, soid)
else:
consistent_tangent = e_consistent_tangent
# Build consistent tangent modulus matricial form
consistent_tangent_mf = vget_tensor_mf(consistent_tangent, n_dim,
comp_order_sym,
device=self._device)
#
# 3D > 2D Conversion
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# When the problem type corresponds to a 2D analysis, build the 2D
# strain and stress tensors (matricial form) once the state update has
# been performed
if self._problem_type == 1:
# Builds 2D strain and stress tensors (matricial form) from the
# associated 3D counterparts
e_trial_strain_mf = vget_state_2Dmf_from_3Dmf(
e_trial_strain_mf, device=self._device)
e_strain_mf = vget_state_2Dmf_from_3Dmf(
e_strain_mf, device=self._device)
stress_mf = vget_state_2Dmf_from_3Dmf(
stress_mf, device=self._device)
#
# Update state variables
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize state variables dictionary
state_variables = self.state_init()
# Store updated state variables
state_variables['e_strain_mf'] = e_strain_mf
state_variables['acc_p_strain'] = acc_p_strain
state_variables['strain_mf'] = e_trial_strain_mf + p_strain_old_mf
state_variables['stress_mf'] = stress_mf
state_variables['is_su_fail'] = is_su_fail
state_variables['is_plast'] = is_plast
state_variables['is_apex_return'] = is_apex_return
if self._problem_type == 1:
state_variables['e_strain_33'] = e_strain_33
state_variables['stress_33'] = stress_33
#
# 3D > 2D Conversion
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# When the problem type corresponds to a 2D analysis, build the 2D
# consistent tangent modulus (matricial form) from the 3D counterpart
if self._problem_type == 1:
consistent_tangent_mf = vget_state_2Dmf_from_3Dmf(
consistent_tangent_mf, device=self._device)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Restore floating-point precision
if is_precision_conversion:
# Reset default floating-point precision
torch.set_default_dtype(torch.float32)
# Perform floating-point precision conversion
state_variables = convert_dict_to_float32(state_variables,
is_inplace=True)
consistent_tangent_mf = \
convert_tensor_to_float32(consistent_tangent_mf)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
return state_variables, consistent_tangent_mf