''' Tools for Elasticity '''
from enum import Enum
from abc import ABC, abstractmethod
import warnings
import numpy as np
import scipy as sp
from pyxfem.postpro.interpolator import ConstantPerTriInterpolator, LinearPerTriInterpolator, get_trifinder
[docs]
class MaterialHypothesis(Enum):
""" Enum class to define the avalaible Material hypothesis"""
PLANESTRESS = 1
PLANESTRAIN = 2
DIM3 = 3
PLANESTRESS = MaterialHypothesis.PLANESTRESS
PLANESTRAIN = MaterialHypothesis.PLANESTRAIN
DIM3 = MaterialHypothesis.DIM3
HYP2DIM = {PLANESTRESS:2, PLANESTRAIN:2, DIM3:3}
[docs]
def YoungPoisson2LambdaMu(E, nu): # pylint: disable=C0103
""" Return Lame coefficients (Lambda, mu) from Young Modulus (E) and poisson ration nu """
lamb, mu = (E*nu/(1.+nu)/(1.-2*nu), E/2./(1+nu))
return lamb, mu
[docs]
class LinearIsotropeElasticMaterial():
''' Linear isotrope elastic material '''
def __init__(self, E:float = 1., nu:float= 0.3, # pylint: disable=C0103
hyp2d:MaterialHypothesis = PLANESTRAIN,
elastic_limit:float = 1.):
self.E = E # pylint: disable=C0103
self.nu = nu
self.lamb, self.mu = YoungPoisson2LambdaMu(E, nu)
self.hypothesis = hyp2d
self.elastic_limit = elastic_limit
if hyp2d == PLANESTRAIN :
kappa = 3-4*nu # kolossov constant en deformation plane (echantillon massif)
self.C2222CSR = planeStrainC2222CSR(lam = self.lamb, mu = self.mu) # pylint: disable=C0103
elif hyp2d == PLANESTRESS :
kappa = (3-nu)/(1+nu)
warnings.warn("warning, material tangent operator not definined yet for PLANESTRESS")
self.C2222CSR = None
else :
kappa = None
self.kappa = kappa
self.C2222CSR = planeStrainC2222CSR(lam = self.lamb, mu = self.mu)
[docs]
class Tensor2222CSR():
''' grouping CSR matrix repesentation of order 4 tensor with base space of dimension 2
For now we have, for A a second order tensor (repesented as a 2d matrix in rowmajor storage)
- ISym: .5 (AT + A) = ISym: A -> ISym.dot(A.flatten()).reshape((2, 2))
- IxI I*trace(A) = IxI: A -> IxI.dot(A.flatten()).reshape((2, 2))
'''
_data = np.array([1., .5, .5, .5, .5, 1.])
_row = np.array([0, 1, 2, 1, 2, 3])
_col = np.array([0, 1, 1, 2, 2, 3])
ISym = sp.sparse.csr_array((_data, (_row, _col)))
_data = np.array([1., 1., 1., 1.])
_col = np.array([0, 0, 3, 3])
_row = np.array([0, 3, 0, 3])
IxI = sp.sparse.csr_array((_data, (_row, _col)))
ISym = Tensor2222CSR.ISym
IxI = Tensor2222CSR.IxI
[docs]
def planeStrainC2222CSR(*, lam:float|None = None, mu:float|None = None, # pylint: disable=C0103
mat:LinearIsotropeElasticMaterial|None = None ):
''' Compute the C operator (tensor 4) for isotropic material
such as sigma = C:eps, for plane strain isotrope assumption
input : lambda and mu Lamé's coefficient or a linearIsotropeElasticMaterial instance.
'''
if lam is None or mu is None :
if mat is not None:
lam = mat.lamb
mu = mat.mu
else :
raise ValueError
return lam*IxI + 2*mu*ISym
[docs]
def stress_from_strain(strain, lamb, mu):
'''compute stress from plane deformation Strain. return sigma2d in voigt format and sigmazz'''
ltr = lamb*(strain[0] + strain[1])
sigmaxx = ltr + 2.*mu*strain[0]
sigmayy = ltr + 2.*mu*strain[1]
sigmaxy = 2.*mu*strain[2]
sigmazz = ltr
print(sigmaxx.shape, sigmayy.shape, sigmaxy.shape)
sigma2d = np.array([sigmaxx, sigmayy, sigmaxy])
return sigma2d, sigmazz
[docs]
def VonMises(stress2d, sigmazz): # pylint: disable=C0103
''' return the von mises stress '''
sxx = stress2d[0]
syy = stress2d[1]
sxy = stress2d[2]
szz =sigmazz
p = (sxx+syy+szz)/3
dxx = sxx-p
dyy = syy-p
dzz = szz-p
dxy = sxy
return np.sqrt(dyy**2+dxx**2+dzz**2+2*dxy**2)
[docs]
def to_voigt2d(s11, s22, s12):
''' stack the 3 stress component in voigt format '''
return np.stack([s11, s22, s12], axis =-1)
[docs]
def to_tensor2d(s11, s22, s12):
''' stack the 3 stress component in 2d tenso2r format '''
return np.stack([s11, s12, s12, s22], axis =-1).reshape( s11.shape + (2,2))
[docs]
class ElasticField(ABC):
''' base class for elastic field : thing thah can evaluate displacement and stress
(and strain ..) at each point xy
'''
[docs]
@abstractmethod
def displacement(self, xy):
''' return the displacement vector at points xy '''
raise NotImplementedError
[docs]
@abstractmethod
def stress(self, xy):
''' return the stress tensor at points xy '''
raise NotImplementedError
[docs]
@abstractmethod
def strain(self, xy):
''' return the strain tensor at points xy '''
raise NotImplementedError
[docs]
class P1FemElasticfield(ElasticField):
''' implementation of ElasticField that use the classical
P1 Fem interpolation (over triangles) '''
def __init__(self, xy, t2v, vdisp, estress, trifinder = None):
if trifinder is None:
trifinder = get_trifinder(xy=xy, t2v =t2v)
self.stress_interp = ConstantPerTriInterpolator(xy, t2v, estress, trifinder = trifinder)
self.disp_interp = LinearPerTriInterpolator(xy, t2v, values = vdisp,
trifinder = trifinder)
[docs]
def displacement(self, xy):
return self.disp_interp(xy)
[docs]
def stress(self, xy):
return self.stress_interp(xy)
[docs]
def strain(self, xy):
raise NotImplementedError