Coverage for pyxfem/solidmechanics/elasticity.py: 46%

1''' Tools for Elasticity ''' 

2from enum import Enum 

3from abc import ABC, abstractmethod 

4import warnings 

5import numpy as np 

6import scipy as sp 

7from pyxfem.postpro.interpolator import ConstantPerTriInterpolator, LinearPerTriInterpolator, get_trifinder 

8 

9 

10class MaterialHypothesis(Enum): 

11 """ Enum class to define the avalaible Material hypothesis""" 

12 PLANESTRESS = 1 

13 PLANESTRAIN = 2 

14 DIM3 = 3 

15 

16PLANESTRESS = MaterialHypothesis.PLANESTRESS 

17PLANESTRAIN = MaterialHypothesis.PLANESTRAIN 

18DIM3 = MaterialHypothesis.DIM3 

19HYP2DIM = {PLANESTRESS:2, PLANESTRAIN:2, DIM3:3} 

20 

21def YoungPoisson2LambdaMu(E, nu): # pylint: disable=C0103 

22 """ Return Lame coefficients (Lambda, mu) from Young Modulus (E) and poisson ration nu """ 

23 lamb, mu = (E*nu/(1.+nu)/(1.-2*nu), E/2./(1+nu)) 

24 return lamb, mu 

25 

26class LinearIsotropeElasticMaterial(): 

27 ''' Linear isotrope elastic material ''' 

28 def __init__(self, E:float = 1., nu:float= 0.3, # pylint: disable=C0103 

29 hyp2d:MaterialHypothesis = PLANESTRAIN, 

30 elastic_limit:float = 1.): 

31 self.E = E # pylint: disable=C0103 

32 self.nu = nu 

33 self.lamb, self.mu = YoungPoisson2LambdaMu(E, nu) 

34 self.hypothesis = hyp2d 

35 self.elastic_limit = elastic_limit 

36 if hyp2d == PLANESTRAIN : 

37 kappa = 3-4*nu # kolossov constant en deformation plane (echantillon massif) 

38 self.C2222CSR = planeStrainC2222CSR(lam = self.lamb, mu = self.mu) # pylint: disable=C0103 

39 elif hyp2d == PLANESTRESS : 

40 kappa = (3-nu)/(1+nu) 

41 warnings.warn("warning, material tangent operator not definined yet for PLANESTRESS") 

42 self.C2222CSR = None 

43 else : 

44 kappa = None 

45 self.kappa = kappa 

46 self.C2222CSR = planeStrainC2222CSR(lam = self.lamb, mu = self.mu) 

47 

48class Tensor2222CSR(): 

49 ''' grouping CSR matrix repesentation of order 4 tensor with base space of dimension 2 

50 

51 For now we have, for A a second order tensor (repesented as a 2d matrix in rowmajor storage) 

52 - ISym: .5 (AT + A) = ISym: A -> ISym.dot(A.flatten()).reshape((2, 2)) 

53 - IxI I*trace(A) = IxI: A -> IxI.dot(A.flatten()).reshape((2, 2)) 

54 ''' 

55 _data = np.array([1., .5, .5, .5, .5, 1.]) 

56 _row = np.array([0, 1, 2, 1, 2, 3]) 

57 _col = np.array([0, 1, 1, 2, 2, 3]) 

58 ISym = sp.sparse.csr_array((_data, (_row, _col))) 

59 _data = np.array([1., 1., 1., 1.]) 

60 _col = np.array([0, 0, 3, 3]) 

61 _row = np.array([0, 3, 0, 3]) 

62 IxI = sp.sparse.csr_array((_data, (_row, _col))) 

63 

64ISym = Tensor2222CSR.ISym 

65IxI = Tensor2222CSR.IxI 

66 

67def planeStrainC2222CSR(*, lam:float|None = None, mu:float|None = None, # pylint: disable=C0103 

68 mat:LinearIsotropeElasticMaterial|None = None ): 

69 ''' Compute the C operator (tensor 4) for isotropic material 

70 such as sigma = C:eps, for plane strain isotrope assumption 

71 input : lambda and mu Lamé's coefficient or a linearIsotropeElasticMaterial instance. 

72 ''' 

73 if lam is None or mu is None : 

74 if mat is not None: 

75 lam = mat.lamb 

76 mu = mat.mu 

77 else : 

78 raise ValueError 

79 return lam*IxI + 2*mu*ISym 

80 

81def stress_from_strain(strain, lamb, mu): 

82 '''compute stress from plane deformation Strain. return sigma2d in voigt format and sigmazz''' 

83 ltr = lamb*(strain[0] + strain[1]) 

84 sigmaxx = ltr + 2.*mu*strain[0] 

85 sigmayy = ltr + 2.*mu*strain[1] 

86 sigmaxy = 2.*mu*strain[2] 

87 sigmazz = ltr 

88 print(sigmaxx.shape, sigmayy.shape, sigmaxy.shape) 

89 sigma2d = np.array([sigmaxx, sigmayy, sigmaxy]) 

90 return sigma2d, sigmazz 

91 

92def VonMises(stress2d, sigmazz): # pylint: disable=C0103 

93 ''' return the von mises stress ''' 

94 sxx = stress2d[0] 

95 syy = stress2d[1] 

96 sxy = stress2d[2] 

97 szz =sigmazz 

98 p = (sxx+syy+szz)/3 

99 dxx = sxx-p 

100 dyy = syy-p 

101 dzz = szz-p 

102 dxy = sxy 

103 return np.sqrt(dyy**2+dxx**2+dzz**2+2*dxy**2) 

104 

105def to_voigt2d(s11, s22, s12): 

106 ''' stack the 3 stress component in voigt format ''' 

107 return np.stack([s11, s22, s12], axis =-1) 

108 

109def to_tensor2d(s11, s22, s12): 

110 ''' stack the 3 stress component in 2d tenso2r format ''' 

111 return np.stack([s11, s12, s12, s22], axis =-1).reshape( s11.shape + (2,2)) 

112 

113class ElasticField(ABC): 

114 ''' base class for elastic field : thing thah can evaluate displacement and stress 

115 (and strain ..) at each point xy 

116 ''' 

117 @abstractmethod 

118 def displacement(self, xy): 

119 ''' return the displacement vector at points xy ''' 

120 raise NotImplementedError 

121 @abstractmethod 

122 def stress(self, xy): 

123 ''' return the stress tensor at points xy ''' 

124 raise NotImplementedError 

125 @abstractmethod 

126 def strain(self, xy): 

127 ''' return the strain tensor at points xy ''' 

128 raise NotImplementedError 

129 

130class P1FemElasticfield(ElasticField): 

131 ''' implementation of ElasticField that use the classical 

132 P1 Fem interpolation (over triangles) ''' 

133 def __init__(self, xy, t2v, vdisp, estress, trifinder = None): 

134 if trifinder is None: 

135 trifinder = get_trifinder(xy=xy, t2v =t2v) 

136 self.stress_interp = ConstantPerTriInterpolator(xy, t2v, estress, trifinder = trifinder) 

137 self.disp_interp = LinearPerTriInterpolator(xy, t2v, values = vdisp, 

138 trifinder = trifinder) 

139 def displacement(self, xy): 

140 return self.disp_interp(xy) 

141 def stress(self, xy): 

142 return self.stress_interp(xy) 

143 def strain(self, xy): 

144 raise NotImplementedError