Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/modelling/cracksol.py: 97%

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Analytical crack solutions for surface displacements and fault dislocations. 

8''' 

9 

10import numpy as num 

11import logging 

12 

13from pyrocko.guts import Float, Object 

14from pyrocko.guts_array import Array 

15 

16guts_prefix = 'modelling' 

17 

18logger = logging.getLogger(__name__) 

19 

20 

21class GriffithCrack(Object): 

22 ''' 

23 Analytical Griffith crack model. 

24 ''' 

25 

26 width = Float.T( 

27 help='Width equals to :math:`2 \\cdot a`.', 

28 default=1.) 

29 

30 poisson = Float.T( 

31 help='Poisson ratio :math:`\\nu`.', 

32 default=.25) 

33 

34 shearmod = Float.T( 

35 help='Shear modulus :math:`\\mu` [Pa].', 

36 default=1.e9) 

37 

38 stressdrop = Array.T( 

39 help='Stress drop array:' 

40 '[:math:`\\sigma_{3,r} - \\sigma_{3,c}`, ' 

41 ':math:`\\sigma_{2,r} - \\sigma_{2,c}`, ' 

42 ':math:`\\sigma_{1,r} - \\sigma_{1,c}`] :math:`=` ' 

43 '[:math:`\\Delta \\sigma_{strike}, ' 

44 '\\Delta \\sigma_{dip}, \\Delta \\sigma_{tensile}`].', 

45 default=num.array([0., 0., 0.])) 

46 

47 @property 

48 def a(self): 

49 ''' 

50 Half width of the crack in [m]. 

51 ''' 

52 

53 return self.width / 2. 

54 

55 def disloc_infinite2d(self, x_obs): 

56 ''' 

57 Calculation of dislocation at crack surface along x2 axis. 

58 

59 Follows equations by Pollard and Segall (1987) to calculate 

60 dislocations for an infinite 2D crack extended in x3 direction, 

61 opening in x1 direction and the crack extending in x2 direction. 

62 

63 :param x_obs: 

64 Observation point coordinates along x2-axis. 

65 If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output dislocations 

66 are zero. 

67 :type x_obs: 

68 :py:class:`~numpy.ndarray`: ``(N,)`` 

69 

70 :return: 

71 Dislocations at each observation point in strike, dip and 

72 tensile direction. 

73 :rtype: 

74 :py:class:`~numpy.ndarray`: ``(N, 3)`` 

75 ''' 

76 

77 if type(x_obs) is not num.ndarray: 

78 x_obs = num.array(x_obs) 

79 

80 factor = num.array([2. / self.shearmod]) 

81 factor = num.append( 

82 factor, num.tile( 

83 2. * (1. - self.poisson) / self.shearmod, (1, 2))) 

84 factor[1] *= -1. 

85 

86 crack_el = (x_obs > -self.a) | (x_obs < self.a) 

87 

88 disl = num.zeros((x_obs.shape[0], 3)) 

89 disl[crack_el, :] = \ 

90 self.stressdrop * num.sqrt( 

91 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \ 

92 factor 

93 

94 return disl 

95 

96 def disloc_circular(self, x_obs): 

97 ''' 

98 Calculation of dislocation at crack surface along x2 axis. 

99 

100 Follows equations by Pollard and Segall (1987) to calculate 

101 displacements for a circulat crack extended in x2 and x3 direction and 

102 opening in x1 direction. 

103 

104 :param x_obs: 

105 Observation point coordinates along axis through crack 

106 centre. If :math:`x_{obs} < -a` or :math:`x_{obs} > a`, output 

107 dislocations are zero. 

108 :type x_obs: 

109 :py:class:`~numpy.ndarray`: ``(N,)`` 

110 

111 :return: 

112 Dislocations at each observation point in strike, dip and 

113 tensile direction. 

114 :rtype: 

115 :py:class:`~numpy.ndarray`: ``(N, 3)`` 

116 ''' 

117 

118 if type(x_obs) is not num.ndarray: 

119 x_obs = num.array(x_obs) 

120 

121 factor = num.array([4. / (self.shearmod * num.pi)]) 

122 factor = num.append( 

123 factor, num.tile( 

124 4. * (1. - self.poisson) / (self.shearmod * num.pi), (1, 2))) 

125 factor[1] *= -1. 

126 

127 crack_el = (x_obs > -self.a) | (x_obs < self.a) 

128 

129 disl = num.zeros((x_obs.shape[0], 3)) 

130 disl[crack_el] = \ 

131 self.stressdrop * num.sqrt( 

132 self.a**2 - num.tile(x_obs[crack_el, num.newaxis], (1, 3))**2) * \ 

133 factor 

134 

135 return disl 

136 

137 def _displ_infinite2d_along_x1(self, x1_obs): 

138 ''' 

139 Calculation of displacement at crack surface along x1-axis. 

140 

141 Follows equations by Pollard and Segall (1987) to calculate 

142 displacements for an infinite 2D crack extended in x3 direction, 

143 opening in x1 direction and the crack tip in x2 direction. 

144 

145 :param x1_obs: 

146 Observation point coordinates along x1-axis. 

147 :type x1_obs: 

148 :py:class:`~numpy.ndarray`: ``(N,)`` 

149 

150 :return: 

151 Displacements at each observation point in strike, dip and 

152 tensile direction. 

153 :rtype: 

154 :py:class:`~numpy.ndarray`: ``(M, 3)`` 

155 ''' 

156 displ = num.zeros((x1_obs.shape[0], 3)) 

157 

158 if self.stressdrop[0] != 0.: 

159 sign = num.sign(x1_obs) 

160 x1_ratio = x1_obs / self.a 

161 

162 displ[:, 0] = ( 

163 num.sqrt((x1_ratio)**2 + 1.) - num.abs( 

164 x1_ratio) 

165 ) * self.stressdrop[0] * self.a / self.shearmod * sign 

166 

167 return displ 

168 

169 def _displ_infinite2d_along_x2(self, x2_obs): 

170 ''' 

171 Calculation of displacement at crack surface along x2-axis. 

172 

173 Follows equations by Pollard and Segall (1987) to calculate 

174 displacements for an infinite 2D crack extended in x3 direction, 

175 opening in x1 direction and the crack tip in x2 direction. 

176 

177 :param x2_obs: 

178 Observation point coordinates along x2-axis. 

179 :type x2_obs: 

180 :py:class:`~numpy.ndarray`: ``(N,)`` 

181 

182 :return: 

183 Displacements at each observation point in strike, dip and 

184 tensile direction. 

185 :rtype: 

186 :py:class:`~numpy.ndarray`: ``(N, 3)`` 

187 ''' 

188 

189 crack_el = (x2_obs >= -self.a) & (x2_obs <= self.a) 

190 

191 displ = num.zeros((x2_obs.shape[0], 3)) 

192 

193 if self.stressdrop[1] != 0.: 

194 factor = (1. - 2. * self.poisson) / (2. * self.shearmod) 

195 

196 displ[crack_el, 2] = \ 

197 self.stressdrop[1] * factor * x2_obs[crack_el] 

198 

199 sign = num.sign(x2_obs) 

200 

201 displ[~crack_el, 2] = \ 

202 self.stressdrop[1] * factor * self.a * sign[~crack_el] * ( 

203 num.abs(x2_obs[~crack_el] / self.a) - 

204 num.sqrt(x2_obs[~crack_el]**2 / self.a**2 - 1.)) 

205 

206 return displ 

207 

208 def displ_infinite2d(self, x1_obs, x2_obs): 

209 ''' 

210 Calculation of displacement at crack surface along different axis. 

211 

212 Follows equations by Pollard and Segall (1987) to calculate 

213 displacements for an infinite 2D crack extended in x3 direction, 

214 opening in x1 direction and the crack tip in x2 direction. 

215 

216 :param x1_obs: 

217 Observation point coordinates along x1-axis. 

218 If :math:`x1_obs = 0.`, displacment is calculated along x2-axis. 

219 :type x1_obs: 

220 :py:class:`~numpy.ndarray`: ``(M,)`` 

221 

222 :param x2_obs: 

223 Observation point coordinates along x2-axis. 

224 If :math:`x2_obs = 0.`, displacment is calculated along x1-axis. 

225 :type x2_obs: 

226 :py:class:`~numpy.ndarray`: ``(N,)`` 

227 

228 :return: 

229 Displacements at each observation point in strike, dip and 

230 tensile direction. 

231 :rtype: 

232 :py:class:`~numpy.ndarray`: ``(M, 3)`` or ``(N, 3)`` 

233 ''' 

234 

235 if type(x1_obs) is not num.ndarray: 

236 x1_obs = num.array(x1_obs) 

237 if type(x2_obs) is not num.ndarray: 

238 x2_obs = num.array(x2_obs) 

239 

240 if (x1_obs == 0.).all(): 

241 return self._displ_infinite2d_along_x2(x2_obs) 

242 elif (x2_obs == 0.).all(): 

243 return self._displ_infinite2d_along_x1(x1_obs) 

244 

245 

246__all__ = [ 

247 'GriffithCrack']