Coverage for pyxfem/solidmechanics/lefm.py: 0%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on Mon Mar 27 11:48:03 2023 

5 

6@author: nchevaug 

7Tools for Linear fracture mechanics. 

8""" 

9import functools 

10from numba.cuda import popc 

11import numpy as np 

12from scipy.integrate import quad 

13from pyxfem.core.mesh_simplex_one_level import SMesh 

14from pyxfem.solidmechanics.elasticity import planeStrainC2222CSR, to_tensor2d, to_voigt2d, \ 

15 PLANESTRAIN, LinearIsotropeElasticMaterial, ElasticField, P1FemElasticfield 

16from pyxfem.solidmechanics.mechanics_solver import ElasticitySolver 

17 

18 

19def cartesian2polar(x,y): 

20 ''' return r theta from x and y''' 

21 return np.hypot(x,y), np.arctan2(y,x) 

22 

23def polar2cartesian(r,theta): 

24 ''' return x,y from r, theta ''' 

25 return r*np.cos(theta), r*np.sin(theta) 

26 

27class crack_coordinate_system(): 

28 def __init__(self, crack_tip_xy, crack_angle): 

29 self.crack_tip_xy = crack_tip_xy 

30 self.crack_angle = crack_angle 

31 

32 def global_to_crack_xy(self, xy): 

33 c,s = (np.cos(self.crack_angle), np.sin(self.crack_angle)) 

34 mat = np.array([[c,-s],[s,c]]) 

35 return np.vecmat(xy-self.crack_tip_xy[np.newaxis], mat) 

36 

37 def global_to_crack_rtheta(self, xy): 

38 xyc = self.global_to_crack_xy(xy) 

39 r, theta = cartesian2polar(xyc[...,0], xyc[...,1]) 

40 return r, theta 

41 

42 def crack_xy_to_global(self, xyc): 

43 c,s = (np.cos(self.crack_angle), np.sin(self.crack_angle)) 

44 mat = np.array([[c,s],[-s,c]]) 

45 return np.vecmat(xyc, mat) + self.crack_tip_xy[np.newaxis] 

46 

47 def crack_rtheta_to_global(self, r, theta): 

48 x, y = polar2cartesian(r = r, theta = theta) 

49 return self.crack_xy_to_global(np.stack((x,y), axis =-1)) 

50 

51class WilliamsSerie(): 

52 ''' Tools for crack tip williams serie expansion of the elastic field 

53 Mode 1 williams serie with T and B stress 

54 A0, B0, B1 -> rigid body motion (x, y, theta respectivelly) 

55 ux = A0 * (kappa+1)/2/mu 

56 uy = B0 * (kappa+1)/2/mu 

57 uxr = (kappa+1)/2/mu B2(-rsin(theta)) 

58 uyr = (kappa+1)/2/mu B2( rcos(theta)) 

59 A1 = KI/sqrt(2*pi) 

60 B1 = -KII/sqrt(2*pi) 

61 A2 = 4T 

62 ''' 

63 def __init__(self, material_param : LinearIsotropeElasticMaterial): 

64 self.material_param = material_param 

65 @staticmethod 

66 def fI_tenseur2(n, theta, notation = "tensor2D"): # pylint: disable=C0103 

67 ''' return the part that depends only on theta and n (the term of the serie) 

68 of the stress tensor for mode 1 (openning) 

69 t2.shape = (3, theta.shape) 

70 using voight notation : 

71 t2[..., 0] : s11 

72 t2[;.., 1] : s22 

73 t2[..., 2] : s12 

74 ''' 

75 cn1 = np.cos((n/2-1)*theta) 

76 sn1 = np.sin((n/2-1)*theta) 

77 cn3 = np.cos((n/2-3)*theta) 

78 sn3 = np.sin((n/2-3)*theta) 

79 s11 = n/2*((2+(-1.)**n+n/2)*cn1 - (n/2-1)*cn3) 

80 s22 = n/2*((2-(-1.)**n-n/2)*cn1 + (n/2-1)*cn3) 

81 s12 = n/2*((n/2-1.)*sn3 - (n/2+ (-1.)**n)*sn1) 

82 if notation == "tensor2D" : 

83 return to_tensor2d(s11, s22, s12) 

84 if notation == "voigt2D" : 

85 return to_voigt2d(s11, s22, s12) 

86 raise ValueError 

87 

88 @staticmethod 

89 def fII_tenseur2(n, theta, notation = "tensor2D"): # pylint: disable=C0103 

90 ''' return the part that depends only on theta and n (the term of the serie) 

91 of the stress tensor for mode II (shear) 

92 t2.shape = (3, theta.shape) 

93 using voight notation : 

94 t2[0] : s11 

95 t2[1] : s22 

96 t2[2] : s12 

97 ''' 

98 cn1 = np.cos((n/2-1)*theta) 

99 sn1 = np.sin((n/2-1)*theta) 

100 cn3 = np.cos((n/2-3)*theta) 

101 sn3 = np.sin((n/2-3)*theta) 

102 s11 = n/2*((-2+(-1.)**n-n/2)*sn1 + (n/2-1)*sn3 ) 

103 s22 = n/2*((-2-(-1.)**n+n/2)*sn1 - (n/2-1)*sn3 ) 

104 s12 = n/2*(( n/2-1.)*cn3 - (n/2- (-1.)**n)*cn1 ) 

105 if notation == "voigt2D" : 

106 return to_voigt2d(s11, s22, s12) 

107 if notation == "tensor2D" : 

108 return to_tensor2d(s11, s22, s12) 

109 raise ValueError 

110 

111 @staticmethod 

112 def gI_vector2(kappa, n, theta): # pylint: disable=C0103 

113 ''' Return the part of the displacement field that depend on theta for symetric mode n ''' 

114 u1 = ( kappa + (-1.)**n + n/2)*np.cos(n/2*theta) - n/2*np.cos((n/2-2)*theta) 

115 u2 = ( kappa - (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta) 

116 return np.stack((u1, u2),-1) 

117 

118 @staticmethod 

119 def gII_vector2(kappa, n, theta): # pylint: disable=C0103 

120 ''' Return the part of the displacement field that depend on theta for 

121 anti symetric mode n ''' 

122 u1 = (-kappa + (-1.)**n - n/2)*np.sin(n/2*theta) + n/2*np.sin((n/2-2)*theta) 

123 u2 = ( kappa + (-1)**n - n/2)*np.cos(n/2*theta) + n/2*np.cos((n/2-2)*theta) 

124 return np.stack((u1, u2),-1) 

125 @staticmethod 

126 def u_rtheta_mode(mode, n, r, theta, kappa, mu): 

127 ''' return the displacement corresponding to mode 1 or 2, r^(n/2)''' 

128 r = np.atleast_1d(r) 

129 fr = r**(n/2)/2/mu 

130 if mode == 1: 

131 ftheta = WilliamsSerie.gI_vector2(kappa, n, theta) 

132 elif mode ==2 : 

133 ftheta = WilliamsSerie.gII_vector2(kappa, n, theta) 

134 else : 

135 raise ValueError 

136 return fr[..., np.newaxis]*ftheta 

137 @staticmethod 

138 def u_xy_mode(mode, n, x, y, kappa, mu): 

139 ''' return the displacement corresponding to mode 1 or 2, r^(n/2)''' 

140 r, theta = cartesian2polar(x, y) 

141 return WilliamsSerie.u_rtheta_mode(mode, n, r, theta, kappa, mu) 

142 

143 @staticmethod 

144 def stress_rtheta_mode(mode, n, r, theta): 

145 ''' return the stress corresponding to mode 1 or 2, r^(n/2)''' 

146 r = np.atleast_1d(r) 

147 fr = r**(n/2-1) 

148 if mode == 1: 

149 ftheta = WilliamsSerie.fI_tenseur2(n, theta) 

150 elif mode ==2 : 

151 ftheta = WilliamsSerie.fII_tenseur2(n, theta) 

152 else : 

153 raise ValueError 

154 return fr[..., np.newaxis, np.newaxis]*ftheta 

155 @staticmethod 

156 def stress_xy_mode(mode, n, x, y): 

157 ''' return the stress corresponding to mode 1 or 2, r^(n/2)''' 

158 r, theta = cartesian2polar(x, y) 

159 return WilliamsSerie.stress_rtheta_mode(mode, n, r, theta) 

160 

161 def stress2d_rtheta( self, r=1., theta=np.linspace(-np.pi, np.pi, 10), A=None, B=None): # pylint: disable=C0103 

162 ''' return the stress field for williams expansion defined by coefficient A and B 

163 sigm(r, theta) = An*fn(r)*fIn(theta) + Bn*fn(r)*fIIn(theta) 

164 (crack axi - coordinate) 

165 ''' 

166 if A is None : 

167 A = dict() 

168 if B is None : 

169 B = dict() 

170 r = np.atleast_1d(r) 

171 stress = np.zeros((*theta.shape, 2, 2)) 

172 for (n, a) in A.items() : 

173 stress += a*WilliamsSerie.stress_rtheta_mode(1, n, r, theta) 

174 for (n, b) in B.items() : 

175 stress += b*WilliamsSerie.stress_rtheta_mode(2, n, r, theta) 

176 return stress 

177 

178 def u_rtheta(self, r =1., theta=np.linspace(-np.pi, np.pi, 10), A=None, B=None): # pylint: disable=C0103 

179 ''' 

180 return the displacement field for williams expansion defined by coefficient A and B 

181 u(r, theta) = An*fn(r)*gIn(theta) + Bn*fn(r)*gIIn(theta) 

182 (crack axi - coordinate) 

183 ''' 

184 if A is None : 

185 A = dict() 

186 if B is None : 

187 B = dict() 

188 r = np.atleast_1d(r) 

189 theta = np.atleast_1d(theta) 

190 if np.isscalar(theta) : 

191 u = np.zeros(2) 

192 else : 

193 u = np.zeros(( *theta.shape, 2)) 

194 kappa = self.material_param.kappa 

195 mu = self.material_param.mu 

196 for (n, an) in A.items(): 

197 u+= an*WilliamsSerie.u_rtheta_mode(1, n, r, theta, kappa, mu) 

198 for (n, bn) in B.items() : 

199 u+= bn*WilliamsSerie.u_rtheta_mode(2, n, r, theta, kappa, mu) 

200 return u 

201 

202 def stress_xy(self, xy, A=None, B=None, compute_sigma33 = False): # pylint: disable=C0103 

203 ''' compute the stress from xy position and values A and B of the williams coefficients ''' 

204 if A is None : 

205 A = dict() 

206 if B is None : 

207 B = dict() 

208 r, theta = cartesian2polar(xy[...,0 ], xy[..., 1]) 

209 sigma2d = self.stress2d_rtheta(r, theta, A, B) 

210 if not compute_sigma33 : 

211 return sigma2d 

212 if self.material_param.hypothesis == PLANESTRAIN : 

213 lamb = self.material_param.lamb 

214 mu = self.material_param.mu 

215 sigma33 = lamb/(2*lamb + 2*mu) *(sigma2d[...,0,0] + sigma2d[...,1,1]) 

216 else : 

217 sigma33 = np.zeros(xy.shape[:-1]) 

218 return sigma2d, sigma33 

219 

220 def u_xy(self, xy, A=None, B=None): # pylint: disable=C0103 

221 ''' compute the displacement from xy position and values A and B 

222 of the williams coefficients ''' 

223 r, theta = cartesian2polar(xy[...,0 ], xy[..., 1]) 

224 return self.u_rtheta(r, theta, A, B) 

225 

226class WilliamsField(ElasticField): 

227 ''' Exact field from williams expansion (crack tip is at (0,0), x axis aligned with the crack) 

228 ^ y 

229 | 

230 | 

231 ______o---> x 

232 crack 

233 ''' 

234 def __init__(self, mat, coeffs_modeI=None, coeffs_modeII =None): # pylint: disable=C0103 

235 self.will = WilliamsSerie(mat) 

236 if coeffs_modeI is None : 

237 coeffs_modeI = dict() 

238 if coeffs_modeII is None : 

239 coeffs_modeII = dict() 

240 self.a = coeffs_modeI 

241 self.b = coeffs_modeII 

242 def displacement(self, xy): 

243 ''' Compute the displacement at xy ''' 

244 return self.will.u_xy(xy, self.a, self.b) 

245 def stress(self, xy): 

246 ''' Compute the stress at xy ''' 

247 return self.will.stress_xy(xy, self.a, self.b) 

248 def strain(self, xy): 

249 raise NotImplementedError 

250 

251def buckner_integrand(u1, t1, u2, t2): 

252 ''' return pointwise value of the brickner interaction integrand between two field 

253 on a line u1, t1 and u2, t2 are respectively the displacement and the tension in a direction 

254 normale to the line at the point. 

255 ''' 

256 return np.vecdot(t1,u2) - np.vecdot(t2,u1) 

257 

258def extract_williams_from_field_contour(field:ElasticField, cracktipxy, cracktheta, # pylint: disable=C0103,W0613 

259 mat, modeIterms, modeIIterms, r, 

260 method = 'trapeze', npt = 500, epsabs=1.49e-8, epsrel=1.49e-8): 

261 ''' Extract from the field the williams coefficients (Aterms for mode I, Bterms for mode II) 

262 by integrating the buckner integrand on a circle of radius r, 

263 using : 

264 - trapeze rule with npt points if method =='trapeze' 

265 - adaptive quadrature from scipy.integrate.quad if method =='adapt' 

266 (with parameters epsabs and epsrel are passed to quad) 

267 must pass an object field from which one can call stress(xy) and displacement(xy) 

268 ''' 

269 if modeIterms is None : 

270 modeIterms = [] 

271 if modeIIterms is None : 

272 modeIIterms = [] 

273 mu = mat.mu 

274 kappa = mat.kappa 

275 def normal(theta): 

276 theta = np.atleast_1d(theta) 

277 return np.column_stack ([ np.cos(theta), np.sin(theta)]) 

278 def field1(theta, normal): 

279 theta = np.atleast_1d(theta) 

280 xy = r*normal 

281 disp1 = field.displacement(xy) 

282 tension1 = np.matvec(field.stress(xy), normal) 

283 return disp1, tension1 

284 def field2(theta, normal, mode, n): 

285 disp2 = WilliamsSerie.u_rtheta_mode(mode, n, r, theta, kappa, mu) 

286 stress2 = WilliamsSerie.stress_rtheta_mode(mode, n, r, theta) 

287 tension2 = np.matvec(stress2, normal) 

288 return disp2, tension2 

289 if method == 'adapt' : 

290 def extract(mode, n): 

291 def buckner(theta): 

292 theta = np.atleast_1d(theta) 

293 norm = normal(theta) 

294 disp1, tension1 = field1(theta, norm) 

295 disp2, tension2 = field2(theta, norm, mode, -n) 

296 w = r/np.pi*mu/(kappa+1)/n/(-1)**n 

297 return buckner_integrand(disp1, tension1, disp2, tension2)*w 

298 return quad(buckner, -np.pi, np.pi, epsabs=epsabs, epsrel=epsrel) 

299 elif method == 'trapeze' : 

300 theta = np.linspace(-np.pi, np.pi, npt) 

301 theta_mid = (theta[:-1] + theta[1:])/2. 

302 norm = normal(theta_mid) 

303 disp1, tension1 = field1(theta_mid, norm) 

304 def extract_trap(mode, n): 

305 disp2, tension2 = field2(theta_mid, norm, mode, -n) 

306 w = 2*r*mu/(npt-1)/(kappa+1)/n/(-1)**n 

307 return float(np.sum(buckner_integrand(disp1, tension1, disp2, tension2)* w)) 

308 extract = extract_trap 

309 else : 

310 raise ValueError 

311 a = {n:extract(1,n) for n in modeIterms} 

312 b = {n:extract(2,n) for n in modeIIterms} 

313 return a, b #, disp1, tension1 

314 

315def extract_williams_from_mesh_contour(mesh:SMesh, vdisp, fstress, cracktipxy, cracktheta, mat, modeIterms=None, modeIIterms=None, r= 0.5, plot_bound = False): # pylint: disable=C0103 

316 import pyxfem.postpro.gmsh_post_pro_interface as gp # pylint: disable=C0415 

317 ___, boundary_edges, boundary_edges_flip = construct_domain_boundary(mesh, cracktipxy, r) 

318 edge_normal = mesh.enormal(eids=boundary_edges, flips=boundary_edges_flip) 

319 edge_cog = mesh.ecog(eids = boundary_edges) 

320 edge_lenght = mesh.elength(eids = boundary_edges) 

321 edge_normal = mesh.enormal(eids=boundary_edges, flips=boundary_edges_flip) 

322 edge_disp = np.sum(vdisp[mesh.e2v(eids=boundary_edges)], axis=1)/2. 

323 e2t = mesh.e2t() 

324 edge_face0 = e2t.data[e2t.index[boundary_edges]] 

325 edge_face1 = e2t.data[e2t.index[boundary_edges]+1] 

326 edge_stress = (fstress[edge_face0] + fstress[edge_face1])/2. 

327 modeIterms_res, modeIIterms_res = extract_williams_from_edge_list(edge_cog, edge_lenght, 

328 edge_normal, edge_disp, edge_stress, 

329 mat, modeIterms=modeIterms, modeIIterms = modeIIterms) # pylint: disable=C0103 

330 if (plot_bound): 

331 xy_edgeplot = mesh.e2v_coord(eids = boundary_edges) 

332 gp.listPlotFieldLine(xy_edgeplot, edge_disp, P0=True, 

333 viewname='Displacement on bound', VectorType = 'Arrow') 

334 gp.listPlotFieldLine(xy_edgeplot, edge_stress, P0=True, 

335 viewname='stress on bound', VectorType = 'Arrow') 

336 return modeIterms_res, modeIIterms_res 

337 

338def extract_williams_from_edge_list(edge_cog, edge_lenght, edge_normal, edge_disp, edge_stress, mat, 

339 modeIterms=None, modeIIterms=None): # pylint: disable=C0103 

340 ''' compute williams coefficient from values at the center of each edge of 

341 the contour integrale ''' 

342 if modeIterms is None : 

343 modeIterms = [] 

344 if modeIIterms is None : 

345 modeIIterms = [] 

346 ecog, elen, enormal, edisp1 = edge_cog, edge_lenght, edge_normal, edge_disp 

347 etens1 = np.matvec(edge_stress, enormal) 

348 er, etheta = cartesian2polar(ecog[:,0], ecog[:,1]) 

349 mu, kappa = mat.mu, mat.kappa 

350 def extract(mode, n): 

351 edisp2 = WilliamsSerie.u_rtheta_mode(mode, -n, er, etheta, kappa, mu) 

352 estress2 = WilliamsSerie.stress_rtheta_mode(mode, -n, er, etheta) 

353 etens2 = np.matvec(estress2, enormal) 

354 weight = elen*(mu/(kappa+1)/np.pi/n/(-1)**n) 

355 return float(np.dot(buckner_integrand(edisp1, etens1, edisp2, etens2),weight)) 

356 return {n:extract(1,n) for n in modeIterms}, {n:extract(2,n) for n in modeIIterms} 

357 

358def construct_face_domain(mesh:SMesh, cracktipxy, r0): 

359 ''' Return the list of edges on a contour of a domain of all the 2D triangle that has 

360 at least on node in the dis of radius r0 centered on cracktipxy. 

361 ''' 

362 r = np.linalg.norm( mesh.vcoord() - cracktipxy[np.newaxis, :], axis = 1) 

363 vertices = np.nonzero( r < r0)[0] 

364 v2t = mesh.v2t() 

365 f = np.array([], dtype = np.uint64) 

366 for v in vertices : 

367 fs = np.uint64(v2t.index[v]) 

368 fe = np.uint64(v2t.index[v+1]) 

369 nf =np.array(v2t.data[fs:fe], dtype = np.uint64) 

370 f = np.append(f, nf) 

371 f = np.unique(f) 

372 return f 

373 

374def get_boundary_edges(mesh:SMesh, faces): 

375 ''' return the list of boundary edge of the list of faces. Considered Boundary if inside 

376 the domain mesh (edge connected to 2 triangles) and on the boundary of the list of face 

377 (only one of the two faces of each edge is in the list faces) 

378 also, return orient that tell if each edge as the same orientation as the corresponding edge 

379 of it's connected face (usefull to compute the external normal) 

380 ''' 

381 istriin = np.zeros(mesh.ntris(), np.bool) 

382 istriin[faces] = True 

383 e2t = mesh.e2t() 

384 nbfaces = e2t.ndata_perentry() 

385 internaledges = np.nonzero(nbfaces==2)[0] 

386 e_face0 = e2t.data[ np.array(e2t.index[internaledges], dtype = np.uint64)] 

387 e_face1 = e2t.data[ np.array(e2t.index[internaledges]+1, dtype = np.uint64)] 

388 e_face0_in = istriin[e_face0] 

389 e_face1_in = istriin[e_face1] 

390 edges = np.nonzero( np.logical_xor(e_face0_in, e_face1_in))[0] 

391 orient = np.zeros(edges.shape, dtype = np.bool) 

392 boundary_edges_id = np.array(internaledges[edges], dtype =np.int64) 

393 nbe = boundary_edges_id.shape[0] 

394 orient = np.zeros(nbe, dtype = np.bool) 

395 e_f0 = e2t.data[np.array(e2t.index[boundary_edges_id], dtype = np.uint64)] 

396 e_f1 = e2t.data[np.array(e2t.index[boundary_edges_id]+1, dtype = np.uint64)] 

397 e_f_in = np.where( np.isin(e_f0, faces), e_f0, e_f1) 

398 t2e, t2eflip = mesh.t2e_flip() 

399 ief = np.nonzero(t2e[e_f_in] == boundary_edges_id[:, np.newaxis]) 

400 orient = t2eflip[e_f_in][ief] 

401 return boundary_edges_id, orient 

402 

403def construct_domain_boundary(mesh:SMesh, cracktipxy, r0): 

404 ''' return list of faces boundary edges and orientation of the boundary edge of a domain 

405 containing all the elements that have at least one node in the disk of radius r0 

406 ''' 

407 domain_faces = construct_face_domain(mesh, cracktipxy, r0) 

408 boundary_edges, boundary_edges_flip = get_boundary_edges(mesh, domain_faces) 

409 return domain_faces, boundary_edges, boundary_edges_flip 

410 

411 

412 

413 

414 

415 #edge_stress = field.stress(edge_cog) 

416 

417# def KMPotFourrier(imode, id, kappa): 

418# harm1 = np.array([id, (4-id), (-id)]) 

419# if imode==1 : 

420# amp1 = np.array([kappa, -id/2, id/2 + -1**id]) 

421# elif imode==2 : 

422# amp1 = np.array([-1j*kappa, -1j*(id/2), 1j*(id/2-(-1)**id)]) 

423# else : 

424# print("Error in KMPotFourrier. Mode should be 1 or 2") 

425# raise ValueError 

426# return harm1, amp1 

427 

428# def crack_strain(xs, ys, mode = 1, K = 1., T = 0., B= 0., 

429# elasticcoef = elastic_coeffs(), thetap = 0.): 

430# ''' thetap : angle de la fissure avec les axes ''' 

431# kappa = elasticcoef['kappa'] 

432# mu = elasticcoef['mu'] 

433# scal=2*mu*np.sqrt(2*np.pi) 

434# rs, ts = cartesian2polar(xs,ys) 

435# km= np.array([K,T/(4/np.sqrt(2*np.pi)),B])/scal 

436# epsxx=0 

437# epsyy=0 

438# epsxy=0 

439# for hi in [1,2,3] : 

440# harm1,amp1=KMPotFourrier(mode,hi,kappa) 

441# ft=np.zeros(rs.shape) 

442# fr=np.zeros(rs.shape) 

443# for amp, harm in zip(amp1, harm1): 

444# fr=fr+amp*np.exp(harm*ts*1j/2) 

445# ft=ft+harm*1j/2*amp*np.exp(harm*ts*1j/2) 

446# ft=ft*(rs**(hi/2-1)) 

447# fr=hi/2*fr*(rs**(hi/2-1)) 

448# dfndth = np.imag(ft) 

449# dfndr = np.imag(fr) 

450# dftdth = np.real(ft) 

451# dftdr = np.real(fr) 

452# drdt = np.cos(ts) 

453# drdn = np.sin(ts) 

454# dthdt = -np.sin(ts) 

455# dthdn = np.cos(ts) 

456# dfndt=dfndr*drdt+dfndth*dthdt 

457# dfndn=dfndr*drdn+dfndth*dthdn 

458# dftdt=dftdr*drdt+dftdth*dthdt 

459# dftdn=dftdr*drdn+dftdth*dthdn 

460# c0 = np.cos(thetap) 

461# s0 = np.sin(thetap) 

462# dftdx=dftdt*c0-dftdn*s0 

463# dftdy=dftdt*s0+dftdn*c0 

464# dfndx=dfndt*c0-dfndn*s0 

465# dfndy=dfndt*s0+dfndn*c0 

466# dfxdx=(dftdx*c0-dfndx*s0) 

467# dfxdy=(dftdy*c0-dfndy*s0) 

468# dfydx=(dftdx*s0+dfndx*c0) 

469# dfydy=(dftdy*s0+dfndy*c0) 

470# epsxx=epsxx+km[hi-1]*dfxdx 

471# epsyy=epsyy+km[hi-1]*dfydy 

472# epsxy=epsxy+km[hi-1]*(dfxdy+dfydx)/2 

473# return np.array([epsxx, epsyy, epsxy]) 

474 

475# def crack_VM(xs, ys, mode = 1, K = 1., T = 0., B= 0., 

476# elasticcoef = elastic_coeffs(), thetap = 0.): 

477# eps = crack_strain(xs, ys, mode, K, T, B, elasticcoef, thetap) 

478# sigma2d, sigmazz = stress_from_strain(eps, elasticcoef['lambda'], elasticcoef['mu']) 

479# return VonMises(sigma2d, sigmazz ) 

480 

481def plot_modes(mesh, mat, 

482 modeI:None|dict = None, # pylint: disable=C0103 

483 modeII:None|dict = None): # pylint: disable=C0103 

484 ''' plot each of the williams modes from n = 0 to 4, symetric then unsymetric mode ''' 

485 default_modes = {0:1., 1:1., 2:1., 3:1., 4:1.} 

486 # pylint: disable=C0415 

487 import pyxfem.core.gmsh_mesh_generator as mg 

488 import pyxfem.core.mesh_simplex_one_level as sm 

489 import pyxfem.core.function_spaces as fs 

490 import pyxfem.postpro.gmsh_post_pro_interface as gp 

491 from pyxfem.core.quadrature import T3_nodes, T3_gauss 

492 import pyxfem.core.sparse_tools as st 

493 fem_space = fs.FEMSpaceP1("DisplacementSpace", mesh, vecdim=2) 

494 will = WilliamsSerie(mat) 

495 xy = mesh.getTris2VerticesCoord() 

496 xycog = mesh.getTris2VerticesCog() 

497 NatVertices = fem_space.operator_dof2val_tri(T3_nodes().uv) # pylint: disable=C0103 

498 B = fem_space.operator_dof2grad_tri(T3_gauss(0).uv) # pylint: disable=C0103 

499 J = fem_space.get_tris_mapping().J # pylint: disable=C0103 

500 nf = len(J) 

501 C = planeStrainC2222CSR(mat = mat) # pylint: disable=C0103 

502 C = C.todense().flatten() # pylint: disable=C0103 

503 wD = np.einsum("i, ij -> ij", T3_gauss(0).w[0]*J, C[np.newaxis, :]).reshape((-1, 4, 4)) # pylint: disable=C0103 

504 wD = st.regular_block_diag_to_csr(wD) # pylint: disable=C0103 

505 D = np.broadcast_to(C, (nf, 16)).reshape((-1, 4, 4)) # pylint: disable=C0103 

506 D = st.regular_block_diag_to_csr(D) # pylint: disable=C0103 

507 stress_op = D@B 

508 kappa = mat.kappa 

509 mu = mat.mu 

510 def plots(mode, n, scal = 1.): 

511 def u_xy(mode, n): 

512 return lambda xy : scal*WilliamsSerie.u_xy_mode(mode, n, xy[...,0], xy[..., 1], kappa, mu) 

513 U = fem_space.interpolate(u_xy(mode, n)).flatten() # pylint: disable=C0103 

514 u = (NatVertices@U).reshape((-1, 3, 2)) 

515 stress_calc = (stress_op@U).reshape((-1, 1, 2, 2)) 

516 if mode ==1 : 

517 stress_exac = will.stress_xy(xycog, A = {n:scal}, compute_sigma33 = False) 

518 else : 

519 stress_exac = will.stress_xy(xycog, B = {n:scal}, compute_sigma33 = False) 

520 stress_exac = WilliamsSerie.stress_xy_mode(mode, n, xycog[..., 0], xycog[..., 1]) 

521 #print() 

522 gp.listPlotFieldTri(xy, u, P0=False, 

523 viewname= f"displacement mode {mode:d}, n = {n:d}", VectorType = "Displacement") 

524 gp.listPlotFieldTri(xy, stress_calc, P0=True, 

525 viewname=f'Stress calc mode {mode:d}, n = {n:d}', Range=[0.,10.]) 

526 gp.listPlotFieldTri(xy, stress_exac, P0=True, 

527 viewname=f'Stress exac mode {mode:d}, n = {n:d}', Range=[0.,10.]) 

528 if modeI is None: 

529 modeI = default_modes 

530 if modeII is None: 

531 modeII = default_modes 

532 for mode, coefs in zip([1,2], [modeI, modeII]): 

533 for n, scal in coefs.items(): 

534 print(n, scal) 

535 plots(mode, n, scal) 

536 

537def test_crack_coordinate_system(gmsh_popup = False): 

538 xy = np.array([[0.,0.], 

539 [1.,0.], 

540 [0.,1.]]) 

541 crack_axis = crack_coordinate_system(np.array([0.,0.]), np.pi/6.) 

542 import gmsh 

543 gmsh.initialize() 

544 from pyxfem.core.gmsh_mesh_generator import crack_in_square_contour, crack_in_square, square_with_inner_seg # pylint: disable=W0611 

545 from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh 

546 from pyxfem.postpro.gmsh_post_pro_interface import listPlotFieldTri 

547 

548 gmsh.option.setNumber("PostProcessing.AnimationCycle", 1) 

549 gmsh.option.setNumber('General.Terminal', 0) 

550 #model_name, ____ = crack_in_square_contour(relh = 0.1, r = 0.5) 

551 model_name, ____ = square_with_inner_seg(length=2., relh1=0.1, relh2=0.1, xc0=(-0.5, 0.), xc1=(0.5,0.)) 

552 if gmsh_popup: 

553 gmsh.fltk.run() 

554 

555 mesh = gmshmodel2smesh(model_name, edge_group_names= ['Crack']) 

556 

557 #mesh = gmshModel2sMesh(model_name, edgeGroupeNames= ['contour','bottom', 'top', 'left_top', 'left_bottom', 'right']) 

558 #dirichlet_vertices = mesh.get_vertices_on_edges_groups(['bottom', 'top', 'left_top', 'left_bottom', 'right']) 

559 xy = mesh.t2v_coord() 

560 data = np.ones(xy.shape[:-1]) 

561 listPlotFieldTri(xy, data) 

562 listPlotFieldTri(crack_axis.global_to_crack_xy(xy), data) 

563 

564 if gmsh_popup: 

565 gmsh.fltk.run() 

566 gmsh.finalize() 

567 # = crack_axis.global_to_crack_xy(xy) 

568 

569 # print(xyc) 

570 

571 # crack_axis = crack_coordinate_system(np.array([1.,1.]), np.pi/6.) 

572 # xyc = crack_axis.global_to_crack_xy(xy) 

573 # print(xyc) 

574 # xy = crack_axis.crack_xy_to_global(xyc) 

575 # print(xy) 

576 

577 

578 

579def check_extract_williams_circle(): 

580 ''' check williams coefficient extraction in case of exact elastic fields''' 

581 E = 1. # pylint: disable=C0103 

582 nu = 0.3 

583 se = 1. 

584 mat = LinearIsotropeElasticMaterial(E = E, nu = nu, elastic_limit = se) 

585 A = {1:1., 2:2., 3:3., 4:4., 5:5.} # pylint: disable=C0103 

586 B = {1:-1., 2:-2., 3:-3., 4:-4., 5:-5.} # pylint: disable=C0103 

587 field = WilliamsField(mat, A, B) 

588 print('Input values') 

589 print(A) 

590 print(B) 

591 print("Williams circle trapeze") 

592 extract = functools.partial(extract_williams_from_field_contour, field =field, 

593 cracktipxy=[0.,0.], cracktheta =0., 

594 mat=mat, modeIterms=[1,2,3,4,5,6], modeIIterms=[1,2,3], 

595 r=1., npt = 500, epsabs=1.49e-8, epsrel=1.49e-8) 

596 A, B = extract(field=field, method='trapeze') 

597 print(A) 

598 print(B) 

599 print("Williams circle adapt") 

600 A, B = extract(field=field, method='adapt') 

601 print(A) 

602 print(B) 

603 

604def check_extract_domain_boundary(showdom=False): 

605 ''' check domain construction for interactions integrales ''' 

606 # pylint: disable=C0415 

607 import pyxfem.core.gmsh_mesh_generator as mg 

608 from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh 

609 import pyxfem.postpro.gmsh_post_pro_interface as gp 

610 model_name, ____ = mg.crack_in_square_contour(relh= 0.025, r=0.45) 

611 mesh = gmshmodel2smesh(model_name, edge_group_names= ["contour"], 

612 tri_group_names = ["Omega_in", "Omega_ext"]) 

613 domain_faces, boundary_edges, boundary_edges_flip = \ 

614 construct_domain_boundary(mesh, np.array([0.,0.]), 0.601) 

615 xy_triplot = mesh.t2v_coord(fids = domain_faces) 

616 xy_edgeplot = mesh.e2v_coord(eids = boundary_edges) 

617 edge_normal = mesh.enormal(eids=boundary_edges, flips=boundary_edges_flip) 

618 if showdom: 

619 gp.listPlotFieldTri(xy_triplot, np.ones(domain_faces.shape[0]), P0=True, viewname='Domain') 

620 gp.listPlotFieldLine(xy_edgeplot, np.ones(boundary_edges.shape[0]), P0=True,