# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
#
import os, sys
import numpy as np
try:
from genericROM.External.pyumat import py3umat as pyumat
except ImportError:# pragma: no cover
print("pyumat not available")
from Mordicus.Containers.ConstitutiveLaws.ConstitutiveLawBase import ConstitutiveLawBase
indices = [0,1,2,3,5,4]
[docs]
class ZmatConstitutiveLaw(ConstitutiveLawBase):
"""
Class containing a ZmatConstitutiveLaw
Attributes
----------
constitutiveLawVariables : dict
dictionary with variable names (str) as keys and variables as type
density : float
density of the material
rotation : np.ndarray((3,3))
rotation matrix to apply to the local coordinates (for order-2 tensors)
rotatedMaterial : bool
False if material is not rotated
behavior : str
Z-mat behavior keyword
"""
def __init__(self, set):
assert isinstance(set, str)
super(ZmatConstitutiveLaw, self).__init__(set, "mechanical")
self.constitutiveLawVariables = {}
self.density = None
self.rotation = None
self.rotatedMaterial = None
self.behavior = None
[docs]
def SetConstitutiveLawVariables(self, constitutiveLawVariables):
"""
Sets the constitutiveLawVariables dictionary
Parameters
----------
constitutiveLawVariables : dict
dictionary with variable names (str) as keys and variables as type
"""
self.constitutiveLawVariables = constitutiveLawVariables
[docs]
def SetOneConstitutiveLawVariable(self, var, value):
"""
Sets one variable of the constitutive law
Parameters
----------
var : str
name of the variable
value : custom_data_structure
variable of the constitutive law
"""
self.constitutiveLawVariables[var] = value
[docs]
def SetDensity(self, density):
"""
Sets the density of the constitutive law
Parameters
----------
density : float
density of the material
"""
self.density = density
[docs]
def SetRotation(self, Q):
"""
Sets the rotation matrix of the constitutive law
Parameters
----------
Q : np.ndarray((3,3))
Q of the material
"""
self.rotation = Q
if np.linalg.norm(self.rotation[:2,:2]-np.eye(2)) < 1.e-12:
self.rotatedMaterial = False
else:
self.rotatedMaterial = True
[docs]
def SetBehavior(self, behavior):
"""
Sets the name of the model behavior
Parameters
----------
behavior : str
name of the model behavior
"""
self.behavior = behavior
[docs]
def GetConstitutiveLawVariables(self):
"""
Returns
-------
dict
complete dictionary defining the constitutive law variables
"""
return self.constitutiveLawVariables
[docs]
def GetOneConstitutiveLawVariable(self, var):
"""
Returns one variable of the constitutive law
Parameters
----------
var : str
key of the dictionnary for storing the variable (e.g. name of the variable)
Returns
-------
custom_data_structure
variable of the constitutive law
"""
return self.constitutiveLawVariables[var]
[docs]
def GetDensity(self):
"""
Returns the density of the material
Returns
-------
float
density
"""
return self.density
[docs]
def ComputeConstitutiveLaw(self, temperature, dtemp, stran, dstran, statev):
"""
Main function of the class: computes a new material state using a constitutive law
solver from a previous material state and variations of temperature and strain
Parameters
----------
temperature : 1D np.ndarray or list
temperature at the previous state, at integration points (np.ndarray of size (nbIntPoints) or list of length nbIntPoints)
dtemp : 1D np.ndarray or list
variations of temperature between the previous state and the new state to compute,
at integration points (np.ndarray of size (nbIntPoints) or list of length nbIntPoints)
stran : np.ndarray
strain at the previous state, at integration points (np.ndarray of size (nbIntPoints,nbeOfDualComponents))
dstran : np.ndarray
variations of strain between the previous state and the new state to compute,
at integration points (np.ndarray of size (nbIntPoints,nbeOfDualComponents))
statev : np.ndarray
internal state variables at the previous state, at integration points (np.ndarray of size (nbIntPoints,nbeOfStateVariables))
Returns
-------
np.ndarray
of size (nbIntPoints, nbeOfDualComponents, nbeOfDualComponents) ddsdde: local tangent matrix at the new state
np.ndarray
of size (nbIntPoints, nbeOfDualComponents) stress: stress at the new state
np.ndarray
of size (nbIntPoints, nbeOfStateVariables) statev: internal state variables at the new state
"""
nbIntPoints = stran.shape[0]
stress = np.empty(stran.shape)
ddsdde = np.empty((nbIntPoints , stran.shape[1], stran.shape[1]))
stran = stran[:,indices]
dstran = dstran[:,indices]
if self.rotatedMaterial == True:
stran[:,3:6] *= 0.5
stranMat = VoigtToMatrix2ndOrder(stran)
stranMat = np.einsum('ik,mij,jl->mkl', self.rotation, stranMat, self.rotation, optimize = True)
stran = MatrixToVoigt2ndOrder(stranMat)
stran[:,3:6] *= 2
dstran[:,3:6] *= 0.5
dstranMat = VoigtToMatrix2ndOrder(dstran)
dstranMat = np.einsum('ik,mij,jl->mkl', self.rotation, dstranMat, self.rotation, optimize = True)
dstran = MatrixToVoigt2ndOrder(dstranMat)
dstran[:,3:6] *= 2
for k in range(nbIntPoints):
ddsdde[k,:,:] = self.PyumatCall(k, temperature, dtemp, stran, dstran, statev)
stress[k,:] = self.constitutiveLawVariables['stress']
if self.rotatedMaterial == True:
stressMat = VoigtToMatrix2ndOrder(stress)
stressMat = np.einsum('ki,mij,lj->mkl', self.rotation, stressMat, self.rotation, optimize = True)
stress = MatrixToVoigt2ndOrder(stressMat)
ddsddeMat = VoigtToMatrix4thOrder(ddsdde)
ddsddeMat = np.einsum('km,ln,omnpr,sp,tr->oklst', self.rotation, self.rotation, ddsddeMat, self.rotation, self.rotation, optimize = True)
ddsdde = MatrixToVoigt4thOrder(ddsddeMat)
stress = stress[:,indices]
ddsdde = ddsdde[:,:,indices]
ddsdde = ddsdde[:,indices,:]
return ddsdde, stress, statev
[docs]
def PyumatCall(self, k, temperature, dtemp, stran, dstran, statev):
"""
Computes a new material state using a constitutive law solver from a
previous material state and variations of temperature and strain at
the integration point ranked k
Parameters
----------
temperature : np.ndarray or list
temperature at the previous state, at integration points (np.ndarray of size (nbIntPoints) or list of length nbIntPoints)
dtemp : np.ndarray or list
variations of temperature between the previous state and the new state to compute,
at integration points (np.ndarray of size (nbIntPoints) or list of length nbIntPoints)
stran : np.ndarray
strain at the previous state, at integration points (np.ndarray of size (nbIntPoints,nbeOfDualComponents))
dstran : np.ndarray
variations of strain between the previous state and the new state to compute,
at integration points (np.ndarray of size (nbIntPoints,nbeOfDualComponents))
statev : np.ndarray
internal state variables at the previous state, at integration points (np.ndarray of size (nbIntPoints,nbeOfStateVariables))
Returns
-------
np.ndarray
of size (nbIntPoints, nbeOfDualComponents, nbeOfDualComponents) ddsdde: local tangent matrix at the new state
"""
return pyumat.umat(stress=self.constitutiveLawVariables['stress'],
statev=statev[k,:],
ddsdde=self.constitutiveLawVariables['ddsdde'],
sse=self.constitutiveLawVariables['sse'],
spd=self.constitutiveLawVariables['spd'],
scd=self.constitutiveLawVariables['scd'],
rpl=self.constitutiveLawVariables['rpl'],
ddsddt=self.constitutiveLawVariables['ddsddt'],
drplde=self.constitutiveLawVariables['drplde'],
drpldt=self.constitutiveLawVariables['drpldt'],
stran=stran[k,:],
dstran=dstran[k,:],
time=self.constitutiveLawVariables['timesim'],
dtime=self.constitutiveLawVariables['dtime'],
temp=temperature[k],
dtemp=dtemp[k],
predef=self.constitutiveLawVariables['predef'],
dpred=self.constitutiveLawVariables['dpred'],
cmname=self.constitutiveLawVariables['cmname'],
ndi=self.constitutiveLawVariables['ndi'],
nshr=self.constitutiveLawVariables['nshr'],
ntens=self.constitutiveLawVariables['ntens'],
nstatv=self.constitutiveLawVariables['nstatv'],
props=self.constitutiveLawVariables['props'],
nprops=self.constitutiveLawVariables['nprops'],
coords=self.constitutiveLawVariables['coords'],
drot=self.constitutiveLawVariables['drot'],
pnewdt=self.constitutiveLawVariables['pnewdt'],
celent=self.constitutiveLawVariables['celent'],
dfgrd0=self.constitutiveLawVariables['dfgrd0'],
dfgrd1=self.constitutiveLawVariables['dfgrd1'],
noel=self.constitutiveLawVariables['noel'],
npt=self.constitutiveLawVariables['npt'],
kslay=self.constitutiveLawVariables['kslay'],
kspt=self.constitutiveLawVariables['kspt'],
kstep=self.constitutiveLawVariables['kstep'],
kinc=self.constitutiveLawVariables['kinc'])
[docs]
def UpdateInternalState(self):
"""
Updates the state of the internal variables
Not applicable to the present implementation
"""
return
def __str__(self):
res = "Mechanical ZmatConstitutiveLaw on set "+self.set
return res
[docs]
def MatrixToVoigt2ndOrder(matrix):
#Only 3D tensors
size = matrix.shape[0]
voigt = np.zeros((size,6))
voigt[:,0] = matrix[:,0,0]
voigt[:,1] = matrix[:,1,1]
voigt[:,2] = matrix[:,2,2]
voigt[:,3] = matrix[:,0,1]
voigt[:,4] = matrix[:,0,2]
voigt[:,5] = matrix[:,1,2]
return voigt
[docs]
def VoigtToMatrix2ndOrder(voigt):
#Only 3D tensors
size = voigt.shape[0]
matrix = np.zeros((size, 3, 3))
for i in range(size):
matrix[i,:,:] = np.array([[voigt[i,0],voigt[i,3],voigt[i,4]],
[voigt[i,3],voigt[i,1],voigt[i,5]],
[voigt[i,4],voigt[i,5],voigt[i,2]]])
return matrix
[docs]
def VoigtToMatrix4thOrder(voigt):
#Only 3D tensors
size = voigt.shape[0]
matrix = np.zeros((size, 3, 3, 3, 3))
a = 1.
b = 1.#1./np.sqrt(2.)
c = 1.#0.5
matrix[:,0,0,0,0] = a*voigt[:,0,0]
matrix[:,0,0,1,1] = a*voigt[:,0,1]
matrix[:,0,0,2,2] = a*voigt[:,0,2]
matrix[:,0,0,0,1] = matrix[:,0,0,1,0] = b*voigt[:,0,3]
matrix[:,0,0,1,2] = matrix[:,0,0,2,1] = b*voigt[:,0,4]
matrix[:,0,0,0,2] = matrix[:,0,0,2,0] = b*voigt[:,0,5]
matrix[:,1,1,0,0] = a*voigt[:,1,0]
matrix[:,1,1,1,1] = a*voigt[:,1,1]
matrix[:,1,1,2,2] = a*voigt[:,1,2]
matrix[:,1,1,0,1] = matrix[:,1,1,1,0] = b*voigt[:,1,3]
matrix[:,1,1,1,2] = matrix[:,1,1,2,1] = b*voigt[:,1,4]
matrix[:,1,1,0,2] = matrix[:,1,1,2,0] = b*voigt[:,1,5]
matrix[:,2,2,0,0] = a*voigt[:,2,0]
matrix[:,2,2,1,1] = a*voigt[:,2,1]
matrix[:,2,2,2,2] = a*voigt[:,2,2]
matrix[:,2,2,0,1] = matrix[:,2,2,1,0] = b*voigt[:,2,3]
matrix[:,2,2,1,2] = matrix[:,2,2,2,1] = b*voigt[:,2,4]
matrix[:,2,2,0,2] = matrix[:,2,2,2,0] = b*voigt[:,2,5]
matrix[:,0,1,0,0] = matrix[:,1,0,0,0] = b*voigt[:,3,0]
matrix[:,0,1,1,1] = matrix[:,1,0,1,1] = b*voigt[:,3,1]
matrix[:,0,1,2,2] = matrix[:,1,0,2,2] = b*voigt[:,3,2]
matrix[:,0,1,0,1] = matrix[:,1,0,0,1] = matrix[:,0,1,1,0] = matrix[:,1,0,1,0] = c*voigt[:,3,3]
matrix[:,0,1,1,2] = matrix[:,0,1,2,1] = matrix[:,1,0,1,2] = matrix[:,1,0,2,1] = c*voigt[:,3,4]
matrix[:,0,1,0,2] = matrix[:,0,1,2,0] = matrix[:,1,0,0,2] = matrix[:,1,0,2,0] = c*voigt[:,3,5]
matrix[:,1,2,0,0] = matrix[:,2,1,0,0] = b*voigt[:,4,0]
matrix[:,1,2,1,1] = matrix[:,2,1,1,1] = b*voigt[:,4,1]
matrix[:,1,2,2,2] = matrix[:,2,1,2,2] = b*voigt[:,4,2]
matrix[:,1,2,0,1] = matrix[:,1,2,1,0] = matrix[:,2,1,0,1] = matrix[:,2,1,1,0] = c*voigt[:,4,3]
matrix[:,1,2,1,2] = matrix[:,1,2,2,1] = matrix[:,2,1,1,2] = matrix[:,2,1,2,1] = c*voigt[:,4,4]
matrix[:,1,2,0,2] = matrix[:,1,2,2,0] = matrix[:,2,1,0,2] = matrix[:,2,1,2,0] = c*voigt[:,4,5]
matrix[:,0,2,0,0] = matrix[:,2,0,0,0] = b*voigt[:,5,0]
matrix[:,0,2,1,1] = matrix[:,2,0,1,1] = b*voigt[:,5,1]
matrix[:,0,2,2,2] = matrix[:,2,0,2,2] = b*voigt[:,5,2]
matrix[:,0,2,0,1] = matrix[:,0,2,2,0] = matrix[:,2,0,0,1] = matrix[:,2,0,1,0] = c*voigt[:,5,3]
matrix[:,0,2,1,2] = matrix[:,0,2,2,1] = matrix[:,2,0,1,2] = matrix[:,2,0,2,1] = c*voigt[:,5,4]
matrix[:,0,2,0,2] = matrix[:,2,0,0,2] = matrix[:,0,2,2,0] = matrix[:,2,0,2,0] = c*voigt[:,5,5]
return matrix
[docs]
def MatrixToVoigt4thOrder(matrix):
#Only 3D tensors
size = matrix.shape[0]
voigt = np.zeros((size, 6, 6))
a = 1.
b = 1.#np.sqrt(2.)
c = 1.#2.
voigt[:,0,0] = a*matrix[:,0,0,0,0]
voigt[:,0,1] = a*matrix[:,0,0,1,1]
voigt[:,0,2] = a*matrix[:,0,0,2,2]
voigt[:,0,3] = b*matrix[:,0,0,0,1]
voigt[:,0,4] = b*matrix[:,0,0,1,2]
voigt[:,0,5] = b*matrix[:,0,0,0,2]
voigt[:,1,0] = a*matrix[:,1,1,0,0]
voigt[:,1,1] = a*matrix[:,1,1,1,1]
voigt[:,1,2] = a*matrix[:,1,1,2,2]
voigt[:,1,3] = b*matrix[:,1,1,0,1]
voigt[:,1,4] = b*matrix[:,1,1,1,2]
voigt[:,1,5] = b*matrix[:,1,1,0,2]
voigt[:,2,0] = a*matrix[:,2,2,0,0]
voigt[:,2,1] = a*matrix[:,2,2,1,1]
voigt[:,2,2] = a*matrix[:,2,2,2,2]
voigt[:,2,3] = b*matrix[:,2,2,0,1]
voigt[:,2,4] = b*matrix[:,2,2,1,2]
voigt[:,2,5] = b*matrix[:,2,2,0,2]
voigt[:,3,0] = b*matrix[:,0,1,0,0]
voigt[:,3,1] = b*matrix[:,0,1,1,1]
voigt[:,3,2] = b*matrix[:,0,1,2,2]
voigt[:,3,3] = c*matrix[:,0,1,0,1]
voigt[:,3,4] = c*matrix[:,0,1,1,2]
voigt[:,3,5] = c*matrix[:,0,1,0,2]
voigt[:,4,0] = b*matrix[:,1,2,0,0]
voigt[:,4,1] = b*matrix[:,1,2,1,1]
voigt[:,4,2] = b*matrix[:,1,2,2,2]
voigt[:,4,3] = c*matrix[:,1,2,0,1]
voigt[:,4,4] = c*matrix[:,1,2,1,2]
voigt[:,4,5] = c*matrix[:,1,2,0,2]
voigt[:,5,0] = b*matrix[:,0,2,0,0]
voigt[:,5,1] = b*matrix[:,0,2,1,1]
voigt[:,5,2] = b*matrix[:,0,2,2,2]
voigt[:,5,3] = c*matrix[:,0,2,0,1]
voigt[:,5,4] = c*matrix[:,0,2,1,2]
voigt[:,5,5] = c*matrix[:,0,2,0,2]
return voigt
if __name__ == "__main__":# pragma: no cover
from genericROM import RunTestFile
RunTestFile(__file__)