Coverage for /usr/local/lib/python3.10/site-packages/spam/mesh/unstructured.py: 93%

1# Library of SPAM functions for dealing with unstructured 3D meshes made of tetrahedra 

2# Copyright (C) 2020 SPAM Contributors 

3# 

4# This program is free software: you can redistribute it and/or modify it 

5# under the terms of the GNU General Public License as published by the Free 

6# Software Foundation, either version 3 of the License, or (at your option) 

7# any later version. 

8# 

9# This program is distributed in the hope that it will be useful, but WITHOUT 

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

12# more details. 

13# 

14# You should have received a copy of the GNU General Public License along with 

15# this program. If not, see <http://www.gnu.org/licenses/>. 

16""" 

17This module offers a set of tools for unstructured 3D meshes made of tetrahedra. 

18 

19>>> # import module 

20>>> import spam.mesh 

21>>> spam.mesh.createCuboid() 

22""" 

23 

24import multiprocessing 

25 

26try: 

27 multiprocessing.set_start_method("fork") 

28except RuntimeError: 

29 pass 

30 

31import gmsh 

32import numpy 

33import progressbar 

34import scipy 

35import spam.mesh 

36from spambind.mesh.meshToolkit import triangulateCGAL, countTetrahedraCGAL 

37 

38# Global number of processes 

39nProcessesDefault = multiprocessing.cpu_count() 

40 

41 

42def gmshToSpam(elementsByType, nodesByType): 

43 """Helper function 

44 Converts gmsh mesh data to SPAM format by: 

45 1. reordering by z y x 

46 2. keeping only tetrahedra 

47 

48 Parameters 

49 ---------- 

50 elementsByType: array 

51 Should be the output of `gmsh.model.mesh.getElementsByType(4)` 

52 nodesByType: array 

53 Should be the output of `gmsh.model.mesh.getNodesByType(4)` 

54 

55 Returns 

56 ------- 

57 points: 2D numpy array 

58 The coordinates of the mesh nodes (zyx) 

59 Each line is [zPos, yPos, xPos] 

60 

61 connectivity: 2D numpy array 

62 The connectivity matrix of the tetrahedra elements 

63 Each line is [node1, node2, node3, node4] 

64 """ 

65 

66 # Get connectivity and node coordinates 

67 element_tags, node_tags_element = elementsByType 

68 node_tags, coord, _ = nodesByType 

69 

70 # NOTE: gmsh returns coordinates in xyz. This means that: 

71 # 1. the coordinates array must be switched back to zyx 

72 # 2. the connectivity array must change so that the nodes of each tetrahedron are numbered 

73 # in such a way that the Jacobian is positive. Which means a permutation of the node numbering. 

74 # 3. in some cases the connectivity matrix has some holes in the node numerotations (some nodes are not used) 

75 # it needs to be rebuild for the projection 

76 

77 # get number of nodes and elements 

78 n_elements = element_tags.shape[0] 

79 nodes_set = set(node_tags) 

80 n_nodes = len(nodes_set) 

81 

82 # get new node numbering (without holes) 

83 new_node_numbering = {} 

84 for i, n in enumerate(nodes_set): 

85 new_node_numbering[n] = i 

86 

87 # reshape the connectivity matrix from flatten gmsh output 

88 connectivity = node_tags.reshape((n_elements, 4)) 

89 

90 # create nodes matrix 

91 nodes = numpy.zeros((n_nodes, 3)) 

92 for i, (nodes_4x1, coord_4x3) in enumerate(zip(connectivity, coord.reshape((n_elements, 4, 3)))): 

93 # change connectivity with new orderning 

94 connectivity[i] = [new_node_numbering[n] for n in nodes_4x1] 

95 

96 # fill node vector with new connectivity numbering and switch x&z 

97 for j, n in enumerate(connectivity[i]): 

98 nodes[n] = coord_4x3[j][::-1] 

99 

100 # rearange connectivity 

101 _ = connectivity.copy() 

102 connectivity[:, 1] = _[:, 3] 

103 connectivity[:, 3] = _[:, 1] 

104 

105 i = 0 

106 for e in list(set(connectivity.ravel())): 

107 if e != i: 

108 print("unused node {e}") 

109 i += 1 

110 

111 return nodes, connectivity 

112 

113 

114def getMeshCharacteristicLength(points, connectivity): 

115 """ 

116 Computes the average distance between two nodes of the edges of each elements. 

117 

118 Parameters 

119 ---------- 

120 points: Nx3 array 

121 List of coordinates of the mesh nodes. 

122 

123 connectivity: Mx4 array 

124 Connectivity matrix of the mesh. 

125 

126 Returns 

127 ------- 

128 float: the characteristic length 

129 """ 

130 

131 def _computeDist(n1, n2, points): 

132 d = [(x1 - x2) ** 2 for x1, x2 in zip(points[n1], points[n2])] 

133 return sum(d) ** 0.5 

134 

135 lc = 0.0 

136 for n1, n2, n3, n4 in connectivity: 

137 lc += sum( 

138 [_computeDist(n1, n2, points), _computeDist(n1, n3, points), _computeDist(n1, n4, points), _computeDist(n2, n3, points), _computeDist(n2, n4, points), _computeDist(n3, n4, points)] 

139 ) / (6.0 * len(connectivity)) 

140 

141 return lc 

142 

143 

144def createCuboid( 

145 lengths, 

146 lc, 

147 origin=[0.0, 0.0, 0.0], 

148 periodicity=False, 

149 verbosity=1, 

150 gmshFile=None, 

151 vtkFile=None, 

152 binary=False, 

153 skipOutput=False, 

154): 

155 """ 

156 Creates an unstructured mesh of tetrahedra inside a cuboid. 

157 

158 Parameters 

159 ---------- 

160 lengths: 1D array 

161 The lengths of the cuboid in the three directions 

162 The axis order is zyx 

163 

164 origin: 1D array 

165 The origin of the cuboid (zyx) 

166 

167 lc: float 

168 characteristic length of the elements of the mesh 

169 (`i.e.`, the average distance between two nodes) 

170 

171 periodicity: bool, optional 

172 if periodicity is True, the generated cube will have a periodicity of mesh on surfaces 

173 Default = False 

174 

175 gmshFile: string, optional 

176 If not None, save the gmsh file with name ``gmshFile`` and suffix ``.msh`` 

177 Default = None 

178 

179 vtkFile: string, optional 

180 If not None, save the vtk file with name ``vtkFile`` and suffix ``.vtk`` 

181 Defaut = None 

182 

183 binary: bool, optional 

184 Save files in binary when possible 

185 Default = False 

186 

187 skipOutput: bool, optional 

188 Returns None to save time (only write the files) 

189 Default = False 

190 

191 verbosity: int, optional 

192 Level of information printed on the terminal and the message console. 

193 0: silent except for fatal errors 

194 1: +errors 

195 2: +warnings 

196 3: +direct 

197 4: +information 

198 5: +status 

199 99: +debug 

200 Default = 1 

201 

202 Returns 

203 ------- 

204 points: 2D numpy array 

205 The coordinates of the mesh nodes (zyx) 

206 Each line is [zPos, yPos, xPos] 

207 

208 connectivity: 2D numpy array 

209 The connectivity matrix of the tetrahedra elements 

210 Each line is [node1, node2, node3, node4] 

211 

212 Example 

213 ------- 

214 >>> points, connectivity = spam.mesh.createCuboid((1.0,1.5,2.0), 0.5) 

215 create a mesh in a cuboid of size 1,1.5,2 with a characteristic length of 0.5 

216 """ 

217 

218 # We will switch the input arrays from zyx to xyz before passing them to gmsh 

219 lengths = lengths[::-1] 

220 origin = origin[::-1] 

221 

222 lx, ly, lz = lengths 

223 ox, oy, oz = origin 

224 

225 # https://gmsh.info/doc/texinfo/gmsh.pdf 

226 

227 # initialize mesh 

228 gmsh.initialize() 

229 gmsh.model.add("SPAM cuboid") 

230 

231 # set mesh length options 

232 # gmsh.option.setNumber("Mesh.MeshSizeMin", lc) 

233 # gmsh.option.setNumber("Mesh.MeshSizeMax", lc) 

234 gmsh.option.setNumber("Mesh.Algorithm", 1) # is the delaunay one (it's the gmsh default) 

235 gmsh.option.setNumber("Mesh.Optimize", True) 

236 gmsh.option.setNumber("Mesh.OptimizeNetgen", True) 

237 

238 # gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) 

239 # gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) 

240 # gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) 

241 

242 # set general options 

243 gmsh.option.setNumber("General.Verbosity", verbosity) 

244 gmsh.option.setNumber("Mesh.Binary", binary) 

245 

246 # Create cuboid geometry 

247 gmsh.model.occ.addBox(ox, oy, oz, lx, ly, lz) # create cube 

248 gmsh.model.occ.synchronize() 

249 

250 # set mesh density at all points of the box 

251 gmsh.model.mesh.setSize(gmsh.model.getEntities(0), lc) # set mesh density to the 8 vertices 

252 gmsh.model.occ.synchronize() 

253 

254 if periodicity: 

255 

256 def tr(t): 

257 phi = [1, 0, 0, t[0], 0, 1, 0, t[1], 0, 0, 1, t[2], 0, 0, 0, 1] 

258 phi = numpy.array(phi).astype("<f4") 

259 return phi 

260 

261 gmsh.model.mesh.setPeriodic(2, [2], [1], tr([lx, 0, 0])) # surface -> dim=2, surface 2 set as surface 1 based on translationX 

262 gmsh.model.mesh.setPeriodic(2, [4], [3], tr([0, ly, 0])) # surface -> dim=2, surface 4 set as surface 3 based on translationY 

263 gmsh.model.mesh.setPeriodic(2, [6], [5], tr([0, 0, lz])) # surface -> dim=2, surface 6 set as surface 5 based on translationZ 

264 gmsh.model.occ.synchronize() 

265 

266 # Generate mesh and optimize 

267 gmsh.model.mesh.generate(3) 

268 gmsh.model.occ.synchronize() 

269 

270 # write gmsh/vtk file 

271 if gmshFile is not None: 

272 gmsh.write(f"{gmshFile}.msh") 

273 

274 # can have additional nodes 

275 if vtkFile is not None: 

276 gmsh.write(f"{vtkFile}.vtk") 

277 

278 # finish here if no output 

279 if skipOutput: 

280 gmsh.finalize() 

281 return None 

282 

283 # Get connectivity and node coordinates 

284 nodes, connectivity = gmshToSpam(gmsh.model.mesh.getElementsByType(4), gmsh.model.mesh.getNodesByElementType(4)) 

285 

286 # DEBUG: gmsh GUI 

287 # gmsh.fltk.run() 

288 

289 # done with gmsh 

290 gmsh.finalize() 

291 

292 # if vtkFile is not None: 

293 # meshio.write_points_cells( 

294 # f"{vtkFile}.vtk", 

295 # nodes, 

296 # {"tetra": connectivity}, 

297 # ) 

298 

299 # return coordinates and connectivity matrix 

300 return nodes, connectivity 

301 

302 

303def createCylinder( 

304 centre, 

305 radius, 

306 height, 

307 lc, 

308 zOrigin=0.0, 

309 membrane=0.0, 

310 verbosity=0, 

311 gmshFile=None, 

312 vtkFile=None, 

313 binary=False, 

314 skipOutput=False, 

315): 

316 """ 

317 Creates an unstructured mesh of tetrahedra inside a cylinder. 

318 The height of the cylinder is along the z axis. 

319 

320 Parameters 

321 ---------- 

322 centre: 1D array 

323 The two coordinates of the centre of the base disk (yx). 

324 

325 radius: float 

326 The radius of the base disk. 

327 

328 height: float 

329 The height of the cylinder. 

330 

331 lc: float 

332 characteristic length of the elements of the mesh 

333 (`i.e.`, the average distance between two nodes) 

334 

335 zOrigin: float, default = 0.0 

336 Translate the points coordinates by zOrigin in the z direction. 

337 

338 membrane: float, default = 0.0 

339 Radius of the membrane (pipe added outside of the cylinder). 

340 If membrane < lc the membrane is ignored. 

341 

342 gmshFile: string, optional 

343 If not None, save the gmsh file with name ``gmshFile`` and suffix ``.msh`` 

344 Default = None 

345 

346 vtkFile: string, optional 

347 If not None, save the vtk file with name ``vtkFile`` and suffix ``.vtk`` 

348 Defaut = None 

349 

350 binary: bool, optional 

351 Save files in binary when possible 

352 Default = False 

353 

354 skipOutput: bool, optional 

355 Returns None to save time (only write the files) 

356 Default = False 

357 

358 verbosity: int, optional 

359 GMSH level of information printed on the terminal and the message console. 

360 0: silent except for fatal errors 

361 1: +errors 

362 2: +warnings 

363 3: +direct 

364 4: +information 

365 5: +status 

366 99: +debug 

367 Default = 1 

368 

369 Returns 

370 ------- 

371 points: 2D numpy array 

372 The coordinates of the mesh nodes (zyx) 

373 Each line is [zPos, yPos, xPos] 

374 

375 connectivity: 2D numpy array 

376 The connectivity matrix of the tetrahedra elements 

377 Each line is [node1, node2, node3, node4] 

378 

379 Example 

380 ------- 

381 >>> points, connectivity = spam.mesh.createCylinder((0.0,0.0), 0.5, 2.0, 0.5) 

382 create a mesh in a cylinder of centre 0,0,0 radius, 0.5 and height 2.0 with a characteristic length of 0.5 

383 """ 

384 

385 # unpack 

386 cy, cx = centre 

387 r = radius 

388 

389 # https://gmsh.info/doc/texinfo/gmsh.pdf 

390 

391 # initialize mesh 

392 gmsh.initialize() 

393 gmsh.model.add("SPAM cylinder") 

394 

395 # set mesh length options 

396 gmsh.option.setNumber("Mesh.MeshSizeMin", lc) 

397 gmsh.option.setNumber("Mesh.MeshSizeMax", lc) 

398 gmsh.option.setNumber("Mesh.Algorithm", 5) # is the delaunay one 

399 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) 

400 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) 

401 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) 

402 

403 # set general options 

404 gmsh.option.setNumber("General.Verbosity", verbosity) 

405 gmsh.option.setNumber("Mesh.Binary", binary) 

406 

407 # Create disk surface 

408 p0 = gmsh.model.geo.addPoint(cx, cy, zOrigin, lc, 1) 

409 p1 = gmsh.model.geo.addPoint(cx + r, cy, zOrigin, lc, 2) 

410 p2 = gmsh.model.geo.addPoint(cx, cy + r, zOrigin, lc, 3) 

411 p3 = gmsh.model.geo.addPoint(cx - r, cy, zOrigin, lc, 4) 

412 p4 = gmsh.model.geo.addPoint(cx, cy - r, zOrigin, lc, 5) 

413 c1 = gmsh.model.geo.addCircleArc(p1, p0, p2) 

414 c2 = gmsh.model.geo.addCircleArc(p2, p0, p3) 

415 c3 = gmsh.model.geo.addCircleArc(p3, p0, p4) 

416 c4 = gmsh.model.geo.addCircleArc(p4, p0, p1) 

417 l1 = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4]) 

418 s1 = gmsh.model.geo.addPlaneSurface([l1]) 

419 

420 # Create membrane surface 

421 if membrane > lc: 

422 r += membrane 

423 p5 = gmsh.model.geo.addPoint(cx + r, cy, zOrigin, lc, 6) 

424 p6 = gmsh.model.geo.addPoint(cx, cy + r, zOrigin, lc, 7) 

425 p7 = gmsh.model.geo.addPoint(cx - r, cy, zOrigin, lc, 8) 

426 p8 = gmsh.model.geo.addPoint(cx, cy - r, zOrigin, lc, 9) 

427 c5 = gmsh.model.geo.addCircleArc(p5, p0, p6) 

428 c6 = gmsh.model.geo.addCircleArc(p6, p0, p7) 

429 c7 = gmsh.model.geo.addCircleArc(p7, p0, p8) 

430 c8 = gmsh.model.geo.addCircleArc(p8, p0, p5) 

431 l2 = gmsh.model.geo.addCurveLoop([c5, c6, c7, c8]) 

432 s2 = gmsh.model.geo.addPlaneSurface([l1, l2]) 

433 # elif membrane > : 

434 # warnings.warn(f"membrane thickness is smaller than the characteristic length of the mesh ({membrane} < {lc})") 

435 

436 gmsh.model.geo.synchronize() 

437 

438 # Create volume 

439 gmsh.model.geo.extrude([(2, s1)], 0, 0, height) 

440 if membrane > lc: 

441 gmsh.model.geo.extrude([(2, s2)], 0, 0, height) 

442 

443 # Generate mesh and optimize 

444 gmsh.model.geo.synchronize() 

445 gmsh.model.mesh.generate() 

446 gmsh.model.mesh.optimize("Netgen", True) 

447 

448 # DEBUG: gmsh GUI 

449 # gmsh.fltk.run() 

450 

451 # write gmsh/vtk file 

452 if gmshFile is not None: 

453 gmsh.write(f"{gmshFile}.msh") 

454 

455 # can have additional nodes 

456 if vtkFile is not None: 

457 gmsh.write(f"{vtkFile}.vtk") 

458 

459 # finish here if no output 

460 if skipOutput: 

461 gmsh.finalize() 

462 return None 

463 

464 # Get connectivity and node coordinates 

465 nodes, connectivity = gmshToSpam(gmsh.model.mesh.getElementsByType(4), gmsh.model.mesh.getNodesByElementType(4)) 

466 

467 # done with gmsh 

468 gmsh.finalize() 

469 

470 # if vtkFile is not None: 

471 # meshio.write_points_cells( 

472 # f"{vtkFile}.vtk", 

473 # nodes, 

474 # {"tetra": connectivity}, 

475 # ) 

476 

477 # return coordinates and connectivity matrix 

478 return nodes, connectivity 

479 

480 

481def triangulate(points, alpha=None, weights=None): 

482 """ 

483 Takes a series of points and optionally their weights and returns a tetrahedral connectivity matrix. 

484 

485 If a completely regular grid is passed, there will be trouble, add some tiny noise. 

486 

487 This function uses CGAL's Regular Triangulation in the background (with all weights=1 if they are not passed). 

488 Weights are passed to CGAL's Regular Triangulation directly and so should be squared if they are radii of particles. 

489 

490 Users can optionally pass an alpha value to the function with the goal of carving flat boundary tetrahedra away from 

491 the final mesh. Typical use of the alpha shape could be the removal of flat (almost 0 volume) boundary tetrahedra 

492 from concave/convex boundaries. Another example might be the removal of tetrahedra from an internal void that exists 

493 within the domain. In all cases, the user can imagine the alpha tool as a spherical "scoop" that will remove any 

494 tetrahedra that it is capable of entering. It follows that flat tetrahedra have a very large open face which are 

495 easily entered by large alpha "scoops". Thus, the user should imagine using a very large alpha value (try starting 

496 with 5*domain size) to avoid letting the alpha "scoop" enter well formed tetrahedra. 

497 

498 Consider the following CGAL analogy: The full mesh is an ice cream sundae and the vertices are the 

499 chocolate chips. The value of alpha is the squared radius of the icecream scoop (following the mesh coordinate units) 

500 that will go in and remove any tetrahedra that it can pass into. Positive alpha is user defined, negative alpha 

501 allows CGAL to automatically select a continuous solid (not recommended). 

502 For more information visit: 

503 https://doc.cgal.org/latest/Alpha_shapes_3/index.html 

504 

505 Parameters 

506 ---------- 

507 points : Nx3 2D numpy array of floats 

508 N points to triangulate (Z, Y, X) 

509 

510 weights : numpy vector of N floats 

511 list of weights associated with points. 

512 Default = None (which means all weights = 1). 

513 

514 alpha : float 

515 size of the alpha shape used for carving the mesh. 

516 Zero is no carving. 

517 Negative is CGAL-autocarve which gives an overcut mesh. 

518 positive is a user-selected size, try 5*domain size. 

519 Default = 0 (no carving). 

520 

521 Returns 

522 ------- 

523 connectivity : Mx4 numpy array of unsigned ints 

524 delaunay triangulation with each row containing point numbers 

525 

526 Note 

527 ---- 

528 Function contributed by Robert Caulk (Laboratoire 3SR, Grenoble) 

529 

530 """ 

531 # There are two passes here -- CGAL is all based in templates, and just stores the triangulation 

532 # in a smart way. 

533 # We want the connectivity matrix, which CGAL doesn't know about, so we use a first pass to get the 

534 # number of tetrahedra, so that we can allocate the connectivity matrix in python, and pass it to 

535 # the next pass through the code, which fills the connectivity matrix 

536 

537 if weights is None: 

538 weights = numpy.ones(len(points)) 

539 

540 elif weights.shape[0] != points.shape[0]: 

541 raise Exception("weights array dim1 != points array dim1") 

542 

543 if alpha is None: 

544 alpha = numpy.array([0]) 

545 else: 

546 alpha = numpy.array([alpha]) 

547 

548 points = points.astype("<f4") 

549 weights = weights.astype("<f4") 

550 alpha = alpha.astype("<f4") 

551 

552 # get the number of tetrahedra so we can properly size our connectivityMatrix 

553 nTet = countTetrahedraCGAL(points, weights, alpha) 

554 connectivity = numpy.zeros([nTet, 4], dtype="<u4") 

555 

556 # setup the return types and argument types 

557 triangulateCGAL(points, weights, connectivity, alpha) 

558 

559 return connectivity 

560 

561 

562def projectTetFieldToGrains(points, connectivity, tetField): 

563 """ 

564 Projects/coarse-grains any field defined on tetrahedra onto grains by volume-averaging over 

565 all tetrahedra for which a given grain is a node. 

566 This can be useful for smoothing out a noisy strain field and will not affect the overall agreement between 

567 the average of strains and the macroscopically observed strains (R.C. Hurley has verified this in a 2017 paper). 

568 

569 Parameters 

570 ---------- 

571 points: m x 3 numpy array 

572 M Particles' coordinates (in deformed configuration for strain field) 

573 

574 connectivity: n x 4 numpy array 

575 Delaunay triangulation connectivity generated by spam.mesh.triangulate for example 

576 

577 tetField: n x 3 x 3 numpy array 

578 Any field defined on tetrahedra (e.g., Bagi strains from bagiStrain). 

579 

580 Returns 

581 ------- 

582 grainField: m x 3 x 3 

583 array containing (3x3) arrays of strain 

584 

585 Example 

586 ------- 

587 grainStrain = spam.mesh.projectBagiStrainToGrains(connectivity,bagiStrain[0],points) 

588 Returns strains for each grain. 

589 

590 Notes 

591 ------ 

592 Function contributed by Ryan Hurley (Johns Hopkins University) 

593 

594 """ 

595 # grainStrainVoigt: Ng array containing (6x1) arrays of strain in Voigt notation. 

596 # RCH Oct 2018. 

597 

598 # from progressbar import ProgressBar 

599 # pbar = ProgressBar() 

600 

601 # print("spam.mesh.projectTetFieldToGrains(): Pre-computing volumes...", end='') 

602 tetVolumes2 = spam.mesh.tetVolumes(points, connectivity) 

603 # print("done.") 

604 

605 # Initialize list of grain values 

606 grainField = numpy.zeros(([points.shape[0]] + list(tetField.shape[1:])), dtype="<f4") 

607 # Get all the unique grain IDs in the Deluanay triangulation 

608 

609 # print("spam.mesh.projectTetFieldToGrains(): Projecting tetrahedal to grains...") 

610 # Loop through grains... 

611 # for label in pbar(range(points.shape[0])): 

612 for label in range(points.shape[0]): 

613 # print("label = ", label) 

614 # Get the list of tets for which this label is a node 

615 touchingTets = numpy.where(connectivity == label)[0] 

616 

617 volTot = 0.0 

618 fieldTot = numpy.zeros(list(tetField.shape[1:]), dtype="<f4") 

619 # Loop through list of tets, summing the strains 

620 for touchingTet in touchingTets: 

621 # print("\ttet = ", touchingTet) 

622 vol = tetVolumes2[touchingTet] 

623 # Track total volume 

624 volTot += vol 

625 

626 # Add volume-weighted field 

627 fieldTot += tetField[touchingTet] * vol 

628 

629 # Divide strainTot by volTot 

630 fieldTot = fieldTot / volTot 

631 

632 # Store in particles 

633 grainField[label] = fieldTot 

634 

635 return grainField 

636 

637 

638def BCFieldFromDVCField( 

639 points, 

640 dvcField, 

641 mask=None, 

642 pixelSize=1.0, 

643 meshType="cube", 

644 centre=[0, 0], 

645 radius=1.0, 

646 topBottom=False, 

647 tol=1e-6, 

648 neighbours=4, 

649): 

650 """ 

651 This function imposes boundary conditions coming from a DVC result to the nodes of an unstructured FE mesh. 

652 

653 Parameters 

654 ---------- 

655 points: 2D numpy array 

656 Array of ``n`` node positions of the unstructured mesh 

657 Each line is [nodeNumber, z, y, x] 

658 

659 dvcField: 2D numpy array 

660 Array of ``m`` points of the dvc field 

661 Each line is [zPos, yPos, xPos, zDisp, yDisp, xDisp] 

662 

663 mask: 2D numpy array, optional 

664 Boolean array of ``m`` points of the dvc field 

665 Points with 0 will be ignored for the field interpolation 

666 Default = None (`i.e.` interpolate based on all of the dvc points) 

667 

668 pixelSize: float, optional 

669 physical size of a pixel (`i.e.` 1mm/px) 

670 Default = 1.0 

671 

672 meshType: string, optional 

673 For the moment, valid inputs are ``cube`` and ``cylinder`` 

674 The axis of a cylinder is considered to be ``z`` 

675 Note that if a cylindrical mesh is passed, ``centre`` and ``radius`` are highly recommended 

676 Default = `cube` 

677 

678 centre: float, optional 

679 The centre of the cylinder [y, x] in physical units (`i.e.` mm) 

680 Default = [0, 0] 

681 

682 radius: float, optional 

683 The radius of the cylinder in physical units (`i.e.` mm) 

684 Default = 1.0 

685 

686 topBottom: bool, optional 

687 If boundary conditions are passed only for the top (`i.e.` z=zmax) and bottom (`i.e.` z=zmin) surfaces of the mesh 

688 Default = False 

689 

690 tol: float, optional 

691 Numerical tolerance for floats equality 

692 Default = 1e-6 

693 

694 neighbours: int, , optional 

695 Neighbours for field interpolation 

696 Default = 4 

697 Returns 

698 ------- 

699 bc: 2D numpy array 

700 Boundary node displacements 

701 Each line is [nodeNumber, zPos, yPos, xPos, zDisp, yDisp, xDisp] 

702 

703 WARNING 

704 ------- 

705 1. All coordinates and displacement arrays are ``z``, ``y``, ``x`` 

706 2. The axis of a cylinder is considered to be ``z`` 

707 """ 

708 

709 # STEP 1: find the edge nodes 

710 posMax = [points[:, i].max() for i in range(1, 4)] 

711 posMin = [points[:, i].min() for i in range(1, 4)] 

712 bcNodes = [] 

713 

714 if meshType == "cube": 

715 # extract edge nodes from the mesh 

716 for point in points: 

717 if topBottom: 

718 testMin = abs(point[1] - posMin[0]) < tol 

719 testMax = abs(point[1] - posMax[0]) < tol 

720 else: 

721 testMin = abs(point[1] - posMin[0]) < tol or abs(point[2] - posMin[1]) < tol or abs(point[3] - posMin[2]) < tol 

722 testMax = abs(point[1] - posMax[0]) < tol or abs(point[2] - posMax[1]) < tol or abs(point[3] - posMax[2]) < tol 

723 if testMin or testMax: 

724 bcNodes.append(point) 

725 elif meshType == "cylinder": 

726 # extract edge nodes from the mesh 

727 for point in points: 

728 testZ = abs(point[1] - posMin[0]) < tol or abs(point[1] - posMax[0]) < tol 

729 testXY = False 

730 if not topBottom: 

731 testXY = (point[2] - centre[0]) ** 2 + (point[3] - centre[1]) ** 2 >= (radius - tol) ** 2 

732 if testZ or testXY: 

733 bcNodes.append(point) 

734 

735 bcNodes = numpy.array(bcNodes) 

736 m = bcNodes.shape[0] 

737 

738 # STEP 2: convert dvc field to physical unit 

739 dvcField *= pixelSize 

740 

741 # STEP 3: Interpolate the disp values of FE using a weighted influence of the k nearest neighbours from DVC coord 

742 bcDisp = numpy.zeros((m, 3)) 

743 # create the k-d tree of the coordinates of DVC good points 

744 if mask is None: 

745 mask = numpy.ones(dvcField.shape[0]) 

746 

747 goodPoints = numpy.where(mask == 1) 

748 treeCoord = scipy.spatial.KDTree(dvcField[:, :3][goodPoints]) 

749 

750 for point in range(m): 

751 distance, ind = treeCoord.query(bcNodes[point, 1:], k=neighbours) 

752 

753 # Check if we've hit the same point 

754 if numpy.any(distance == 0): 

755 bcDisp[point, 0] = dvcField[goodPoints][ind][numpy.where(distance == 0)][0, 3] 

756 bcDisp[point, 1] = dvcField[goodPoints][ind][numpy.where(distance == 0)][0, 4] 

757 bcDisp[point, 2] = dvcField[goodPoints][ind][numpy.where(distance == 0)][0, 5] 

758 

759 else: 

760 weightSumInv = sum(1 / distance) 

761 # loop over neighbours 

762 for neighbour in range(neighbours): 

763 # calculate its weight 

764 weightInv = (1 / distance[neighbour]) / float(weightSumInv) 

765 # replace the disp values the weighted disp components of the ith nearest neighbour 

766 bcDisp[point, 0] += dvcField[goodPoints][ind[neighbour], 3] * weightInv 

767 bcDisp[point, 1] += dvcField[goodPoints][ind[neighbour], 4] * weightInv 

768 bcDisp[point, 2] += dvcField[goodPoints][ind[neighbour], 5] * weightInv 

769 

770 # return node number and displacements 

771 return numpy.hstack((bcNodes, bcDisp)) 

772 

773 

774def tetVolumes(points, connectivity): 

775 """ 

776 This function computes volumes of the tetrahedra passed with a connectivity matrix. 

777 Using algorithm in https://en.wikipedia.org/wiki/Tetrahedron#Volume 

778 

779 Parameters 

780 ---------- 

781 points : Nx3 array 

782 Array of ``N`` coordinates of the points 

783 

784 connectivity : Mx4 array 

785 Array of ``M`` none numbers that are connected as tetrahedra (e.g., the output from triangulate) 

786 

787 Returns 

788 ------- 

789 volumes : vector of length M 

790 Volumes of tetrahedra 

791 

792 Note 

793 ----- 

794 Pure python function. 

795 """ 

796 

797 # Sanity checks on lengths: 

798 # if connectivity.shape[1] != 4 or points.shape[1] != 3 or connectivity.max() > points.shape[0]: 

799 # print("spam.mesh.tetVolumes(): Dimensionality problem, not running") 

800 # print("connectivity.max()", connectivity.max(), "points.shape[0]", points.shape[0]) 

801 # return 

802 

803 volumes = numpy.zeros(connectivity.shape[0], dtype="<f4") 

804 connectivity = connectivity.astype(numpy.uint) 

805 

806 # loop through tetrahedra 

807 for tetNumber in range(connectivity.shape[0]): 

808 fourNodes = connectivity[tetNumber] 

809 

810 if max(fourNodes) >= points.shape[0]: 

811 print("spam.mesh.unstructured.tetVolumes(): this tet has node > points.shape[0], skipping.") 

812 else: 

813 if numpy.isfinite(points[fourNodes]).sum() != 3 * 4: 

814 print("spam.mesh.unstructured.bagiStrain(): nans in position, skipping") 

815 else: 

816 a = points[fourNodes[0]] 

817 b = points[fourNodes[1]] 

818 c = points[fourNodes[2]] 

819 d = points[fourNodes[3]] 

820 volumes[tetNumber] = numpy.abs(numpy.dot((a - d), numpy.cross((b - d), (c - d)))) / 6.0 

821 

822 return volumes 

823 

824 

825# def create2Patche(lengths, mesh_size_1, mesh_size_2, origin=[0., 0., 0.], verbose=False, gmshFile=None, vtkFile=None): 

826# import pygmsh 

827# import meshio 

828# 

829# lx, ly, lz = lengths 

830# ox, oy, oz = origin 

831# 

832# # We will switch the input arrays from zyx to xyz before passing them to pygmsh 

833# lengths = lengths[::-1] 

834# origin = origin[::-1] 

835# 

836# # raw code 

837# code = [] 

838# 

839# code.append("SetFactory(\"OpenCASCADE\");") 

840# code.append("Mesh.RandomSeed = 0.5;") 

841# code.append("Point(1) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox, y=oy, z=oz, mesh_size=mesh_size_1)) 

842# code.append("Point(2) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox+lx, y=oy, z=oz, mesh_size=mesh_size_1)) 

843# code.append("Point(3) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox+lx, y=oy+ly, z=oz, mesh_size=mesh_size_1)) 

844# code.append("Point(4) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox, y=oy+ly, z=oz, mesh_size=mesh_size_1)) 

845# code.append("Point(5) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox, y=oy, z=oz+lz, mesh_size=mesh_size_1)) 

846# code.append("Point(6) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox+lx, y=oy, z=oz+lz, mesh_size=mesh_size_1)) 

847# code.append("Point(7) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox+lx, y=oy+ly, z=oz+lz, mesh_size=mesh_size_1)) 

848# code.append("Point(8) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox, y=oy+ly, z=oz+lz, mesh_size=mesh_size_1)) 

849# code.append("Line(1) = {1,2};") 

850# code.append("Line(2) = {2,3};") 

851# code.append("Line(3) = {3,4};") 

852# code.append("Line(4) = {4,1};") 

853# code.append("Line(5) = {5,6};") 

854# code.append("Line(6) = {6,7};") 

855# code.append("Line(7) = {7,8};") 

856# code.append("Line(8) = {8,5};") 

857# code.append("Line(9) = {5,1};") 

858# code.append("Line(10) = {2,6};") 

859# code.append("Line(11) = {7,3};") 

860# code.append("Line(12) = {8,4};") 

861# code.append("Line Loop(13) = { 7, 8, 5, 6 };") 

862# code.append("Plane Surface(14) = {13};") 

863# code.append("Line Loop(15) = {1, 2, 3, 4};") 

864# code.append("Plane Surface(16) = {15};") 

865# code.append("Line Loop(17) = {8, 9, -4, -12};") 

866# code.append("Plane Surface(18) = {17};") 

867# code.append("Line Loop(19) = {1, 10, -5, 9};") 

868# code.append("Plane Surface(20) = {19};") 

869# code.append("Line Loop(21) = {10, 6, 11, -2};") 

870# code.append("Plane Surface(22) = {21};") 

871# code.append("Line Loop(23) = {11, 3, -12, -7};") 

872# code.append("Plane Surface(24) = {23};") 

873# code.append("Surface Loop(25) = {14, 24, 22, 20, 16, 18};") 

874# code.append("Volume(26) = {25};") 

875# code2 = code[:] # this one is for auxiliary (coarse) patch 

876# code2.append("Point(28) = {{ {x}, {y}, {z}, {mesh_size} }};".format(x=ox+lx/2, y=oy+ly/2, z=oz+lz/2, mesh_size=mesh_size_2)) 

877# code2.append("Point{28} In Volume{26};") 

878# e = 1e-6 #### this part is to enforce periodicity conditions for the patch mes 

879# code.append("back_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz+lz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

880# code2.append("back_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz+lz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

881# code.append("front_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+e)) 

882# code2.append("front_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+e)) 

883# code.append("left_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

884# code2.append("left_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

885# code.append("right_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox+lx-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

886# code2.append("right_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox+lx-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

887# code.append("down_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+e,zmax=oz+lz+e)) 

888# code2.append("down_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+e,zmax=oz+lz+e)) 

889# code.append("up_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy+ly-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

890# code2.append("up_surface() = Surface In BoundingBox{{{xmin},{ymin},{zmin},{xmax},{ymax},{zmax}}};".format(xmin=ox-e,ymin=oy+ly-e,zmin=oz-e,xmax=ox+lx+e,ymax=oy+ly+e,zmax=oz+lz+e)) 

891# code.append("Periodic Surface{ right_surface() } = { left_surface() } Translate{ "+str(lx)+",0,0 };") 

892# code2.append("Periodic Surface{ right_surface() } = { left_surface() } Translate{ "+str(lx)+",0,0 };") 

893# code.append("Periodic Surface{ up_surface() } = { down_surface() } Translate{ 0,"+str(ly)+",0 };") 

894# code2.append("Periodic Surface{ up_surface() } = { down_surface() } Translate{ 0,"+str(ly)+",0 };") 

895# code.append("Periodic Surface{ back_surface() } = { front_surface() } Translate{ 0,0,"+str(lz)+" };") 

896# code2.append("Periodic Surface{ back_surface() } = { front_surface() } Translate{ 0,0,"+str(lz)+" };") 

897# 

898# # geom = pygmsh.opencascade.Geometry(characteristic_length_min=mesh_size, characteristic_length_max=mesh_size,) 

899# # geom = pygmsh.geo.Geometry() 

900# geom = pygmsh.opencascade.Geometry() 

901# geom2 = pygmsh.opencascade.Geometry() 

902# 

903# # add raw code to geometry 

904# geom.add_raw_code(code) 

905# geom2.add_raw_code(code2) 

906# 

907# # mesh 

908# # points, cells, _, _, _ = geom.generate_mesh(verbose=verbose, extra_gmsh_arguments=["-optimize_netgen"]) 

909# # points, cells, _, _, _ = geom.generate_mesh(verbose=verbose, extra_gmsh_arguments=["-optimize","-algo","del3d","-clmin",str(mesh_size),"-clmax",str(mesh_size)]) 

910# mesh = geom.generate_mesh(verbose=verbose, extra_gmsh_arguments=["-3", "-optimize", "-algo", "del3d"]) 

911# points = mesh.points 

912# cells = mesh.cells 

913# connectivity = cells['tetra'] 

914# meshaux = geom2.generate_mesh(verbose=verbose, extra_gmsh_arguments=["-3", "-optimize", "-algo", "del3d"]) 

915# pointsaux = meshaux.points 

916# cellsaux = meshaux.cells 

917# connectivityaux = cellsaux['tetra'] 

918# 

919# # write gmsh/vtk file 

920# if gmshFile is not None: 

921# meshio.write_points_cells("{}_fin.msh".format(gmshFile), points, cells, file_format='gmsh22') 

922# meshio.write_points_cells("{}_aux.msh".format(gmshFile), pointsaux, {'tetra': connectivityaux}, file_format='gmsh22') 

923# if vtkFile is not None: 

924# meshio.write_points_cells("{}_fin.vtk".format(vtkFile), points, {'tetra': connectivity}, file_format='vtk') 

925# meshio.write_points_cells("{}_aux.vtk".format(vtkFile), pointsaux, {'tetra': connectivityaux}, file_format='vtk') 

926# 

927# ### TEST pour 1 translation de lx (ox = ox + lx) dans la direction x 

928# # on veut produire : 

929# ### 2 fichiers msh "fin" domain_fin_i.msh (et aussi domain_fin_i.vtk pour la visu) 

930# ### 2 fichiers msh "aux" domain_aux_i.msh (et aussi vtk pour la visu) 

931# ### 1 fichier msh "global" qui reunit les deux fichiers aux (et aussi vtk pour la visu) 

932# points[:, 0] += lx # pour les fichiers "fin" 

933# 

934# temppoints = copy.deepcopy(pointsaux) 

935# pointsaux[:, 0] += lx #translate aux mesh 

936# glob_points = numpy.concatenate((temppoints,pointsaux), axis = 0) #create an array for global mesh (union of aux meshes) 

937# tempconnec = copy.deepcopy(connectivityaux) 

938# connectivityaux[:, :] += pointsaux.shape[0] #translate aux connectivity 

939# glob_connectivity = numpy.concatenate((tempconnec,connectivityaux), axis = 0) #create an array for global mesh (union of aux meshes) 

940# 

941# # write gmsh/vtk file for second fine and aux mesh 

942# if gmshFile is not None: 

943# meshio.write_points_cells("{}_fin_2.msh".format(gmshFile), points, cells, file_format='gmsh22') 

944# meshio.write_points_cells("{}_aux_2.msh".format(gmshFile), pointsaux, {'tetra': tempconnec}, file_format='gmsh22') 

945# if vtkFile is not None: 

946# meshio.write_points_cells("{}_fin_2.vtk".format(vtkFile), points, {'tetra': connectivity}, file_format='vtk') 

947# meshio.write_points_cells("{}_aux_2.vtk".format(vtkFile), pointsaux, {'tetra': tempconnec}, file_format='vtk') # attention ici on ne decale pas la connectivité 

948# 

949# # now try to generate global mesh 

950# # its a three step process 

951# # first, make a list of a list of all nodes that appear more than once in glob_points. 

952# # second, replace each one of those "double" node by the smaller number in the connectivity (glob_connectivity). At this stage all those double nodes are now unlinked to any element. 

953# # third, remove the double nodes from the node list and modify the connectivity accordingly (most difficult step!) 

954# 

955# 

956# print('Start to build GLOBAL mesh - step 1') 

957# double = [] 

958# for i,node1 in enumerate(glob_points): # look at every point in the mesh 

959# for j,node2 in enumerate(glob_points[i+1:]): #and check for existing node at the same place 

960# if ((node1[0]==node2[0]) and (node1[1]==node2[1]) and (node1[2]==node2[2])): 

961# print('Finding double node of coordinates: ', i+1, '=',glob_points[i], j+i+2, '=',glob_points[j+i+1])