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

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 

26import gmsh 

27import numpy 

28import progressbar 

29import scipy 

30import spam.mesh 

31 

32# Global number of processes 

33nProcessesDefault = multiprocessing.cpu_count() 

34 

35 

36def gmshToSpam(elementsByType, nodesByType): 

37 """Helper function 

38 Converts gmsh mesh data to SPAM format by: 

39 1. reordering by z y x 

40 2. keeping only tetrahedra 

41 

42 Parameters 

43 ---------- 

44 elementsByType: array 

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

46 nodesByType: array 

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

48 

49 Returns 

50 ------- 

51 points: 2D numpy array 

52 The coordinates of the mesh nodes (zyx) 

53 Each line is [zPos, yPos, xPos] 

54 

55 connectivity: 2D numpy array 

56 The connectivity matrix of the tetrahedra elements 

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

58 """ 

59 

60 # Get connectivity and node coordinates 

61 element_tags, node_tags_element = elementsByType 

62 node_tags, coord, _ = nodesByType 

63 

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

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

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

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

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

69 # it needs to be rebuild for the projection 

70 

71 # get number of nodes and elements 

72 n_elements = element_tags.shape[0] 

73 nodes_set = set(node_tags) 

74 n_nodes = len(nodes_set) 

75 

76 # get new node numbering (without holes) 

77 new_node_numbering = {} 

78 for i, n in enumerate(nodes_set): 

79 new_node_numbering[n] = i 

80 

81 # reshape the connectivity matrix from flatten gmsh output 

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

83 

84 # create nodes matrix 

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

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

87 # change connectivity with new orderning 

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

89 

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

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

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

93 

94 # rearange connectivity 

95 _ = connectivity.copy() 

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

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

98 

99 i = 0 

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

101 if e != i: 

102 print("unused node {e}") 

103 i += 1 

104 

105 return nodes, connectivity 

106 

107 

108def createCuboid( 

109 lengths, 

110 lc, 

111 origin=[0.0, 0.0, 0.0], 

112 periodicity=False, 

113 verbosity=1, 

114 gmshFile=None, 

115 vtkFile=None, 

116 binary=False, 

117 skipOutput=False, 

118): 

119 """ 

120 Creates an unstructured mesh of tetrahedra inside a cuboid. 

121 

122 Parameters 

123 ---------- 

124 lengths: 1D array 

125 The lengths of the cuboid in the three directions 

126 The axis order is zyx 

127 

128 origin: 1D array 

129 The origin of the cuboid (zyx) 

130 

131 lc: float 

132 characteristic length of the elements of the mesh 

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

134 

135 periodicity: bool, optional 

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

137 Default = False 

138 

139 gmshFile: string, optional 

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

141 Default = None 

142 

143 vtkFile: string, optional 

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

145 Defaut = None 

146 

147 binary: bool, optional 

148 Save files in binary when possible 

149 Default = False 

150 

151 skipOutput: bool, optional 

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

153 Default = False 

154 

155 verbosity: int, optional 

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

157 0: silent except for fatal errors 

158 1: +errors 

159 2: +warnings 

160 3: +direct 

161 4: +information 

162 5: +status 

163 99: +debug 

164 Default = 1 

165 

166 Returns 

167 ------- 

168 points: 2D numpy array 

169 The coordinates of the mesh nodes (zyx) 

170 Each line is [zPos, yPos, xPos] 

171 

172 connectivity: 2D numpy array 

173 The connectivity matrix of the tetrahedra elements 

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

175 

176 Example 

177 ------- 

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

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

180 """ 

181 

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

183 lengths = lengths[::-1] 

184 origin = origin[::-1] 

185 

186 lx, ly, lz = lengths 

187 ox, oy, oz = origin 

188 

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

190 

191 # initialize mesh 

192 gmsh.initialize() 

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

194 

195 # set mesh length options 

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

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

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

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

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

201 

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

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

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

205 

206 # set general options 

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

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

209 

210 # Create cuboid geometry 

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

212 gmsh.model.occ.synchronize() 

213 

214 # set mesh density at all points of the box 

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

216 gmsh.model.occ.synchronize() 

217 

218 if periodicity: 

219 

220 def tr(t): 

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

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

223 return phi 

224 

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

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

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

228 gmsh.model.occ.synchronize() 

229 

230 # Generate mesh and optimize 

231 gmsh.model.mesh.generate(3) 

232 gmsh.model.occ.synchronize() 

233 

234 # write gmsh/vtk file 

235 if gmshFile is not None: 

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

237 

238 # can have additional nodes 

239 if vtkFile is not None: 

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

241 

242 # finish here if no output 

243 if skipOutput: 

244 gmsh.finalize() 

245 return None 

246 

247 # Get connectivity and node coordinates 

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

249 

250 # DEBUG: gmsh GUI 

251 # gmsh.fltk.run() 

252 

253 # done with gmsh 

254 gmsh.finalize() 

255 

256 # if vtkFile is not None: 

257 # meshio.write_points_cells( 

258 # f"{vtkFile}.vtk", 

259 # nodes, 

260 # {"tetra": connectivity}, 

261 # ) 

262 

263 # return coordinates and connectivity matrix 

264 return nodes, connectivity 

265 

266 

267def createCylinder( 

268 centre, 

269 radius, 

270 height, 

271 lc, 

272 zOrigin=0.0, 

273 verbosity=0, 

274 gmshFile=None, 

275 vtkFile=None, 

276 binary=False, 

277 skipOutput=False, 

278): 

279 """ 

280 Creates an unstructured mesh of tetrahedra inside a cylinder. 

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

282 

283 Parameters 

284 ---------- 

285 centre: 1D array 

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

287 

288 radius: float 

289 The radius of the base disk. 

290 

291 height: float 

292 The height of the cylinder. 

293 

294 lc: float 

295 characteristic length of the elements of the mesh 

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

297 

298 zOrigin: float, default = 0.0 

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

300 

301 gmshFile: string, optional 

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

303 Default = None 

304 

305 vtkFile: string, optional 

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

307 Defaut = None 

308 

309 binary: bool, optional 

310 Save files in binary when possible 

311 Default = False 

312 

313 skipOutput: bool, optional 

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

315 Default = False 

316 

317 verbosity: int, optional 

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

319 0: silent except for fatal errors 

320 1: +errors 

321 2: +warnings 

322 3: +direct 

323 4: +information 

324 5: +status 

325 99: +debug 

326 Default = 1 

327 

328 Returns 

329 ------- 

330 points: 2D numpy array 

331 The coordinates of the mesh nodes (zyx) 

332 Each line is [zPos, yPos, xPos] 

333 

334 connectivity: 2D numpy array 

335 The connectivity matrix of the tetrahedra elements 

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

337 

338 Example 

339 ------- 

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

341 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 

342 """ 

343 

344 # unpack 

345 cy, cx = centre 

346 r = radius 

347 

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

349 

350 # initialize mesh 

351 gmsh.initialize() 

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

353 

354 # set mesh length options 

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

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

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

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

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

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

361 

362 # set general options 

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

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

365 

366 # Create disk surface 

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

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

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

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

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

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

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

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

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

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

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

378 gmsh.model.geo.synchronize() 

379 

380 # Create volume surface 

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

382 gmsh.model.geo.synchronize() 

383 

384 # Generate mesh and optimize 

385 gmsh.model.mesh.generate() 

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

387 

388 # write gmsh/vtk file 

389 if gmshFile is not None: 

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

391 

392 # can have additional nodes 

393 if vtkFile is not None: 

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

395 

396 # finish here if no output 

397 if skipOutput: 

398 gmsh.finalize() 

399 return None 

400 

401 # Get connectivity and node coordinates 

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

403 

404 # DEBUG: gmsh GUI 

405 # gmsh.fltk.run() 

406 

407 # done with gmsh 

408 gmsh.finalize() 

409 

410 # if vtkFile is not None: 

411 # meshio.write_points_cells( 

412 # f"{vtkFile}.vtk", 

413 # nodes, 

414 # {"tetra": connectivity}, 

415 # ) 

416 

417 # return coordinates and connectivity matrix 

418 return nodes, connectivity 

419 

420 

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

422 """ 

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

424 

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

426 

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

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

429 

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

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

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

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

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

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

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

437 

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

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

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

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

442 For more information visit: 

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

444 

445 Parameters 

446 ---------- 

447 points : Nx3 2D numpy array of floats 

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

449 

450 weights : numpy vector of N floats 

451 list of weights associated with points. 

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

453 

454 alpha : float 

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

456 Zero is no carving. 

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

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

459 Default = 0 (no carving). 

460 

461 Returns 

462 ------- 

463 connectivity : Mx4 numpy array of unsigned ints 

464 delaunay triangulation with each row containing point numbers 

465 

466 Note 

467 ---- 

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

469 

470 """ 

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

472 # in a smart way. 

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

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

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

476 

477 if weights is None: 

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

479 

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

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

482 

483 if alpha is None: 

484 alpha = numpy.array([0]) 

485 else: 

486 alpha = numpy.array([alpha]) 

487 

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

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

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

491 

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

493 nTet = spam.mesh.meshToolkit.countTetrahedraCGAL(points, weights, alpha) 

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

495 

496 # setup the return types and argument types 

497 spam.mesh.meshToolkit.triangulateCGAL(points, weights, connectivity, alpha) 

498 

499 return connectivity 

500 

501 

502def projectTetFieldToGrains(points, connectivity, tetField): 

503 """ 

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

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

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

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

508 

509 Parameters 

510 ---------- 

511 points: m x 3 numpy array 

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

513 

514 connectivity: n x 4 numpy array 

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

516 

517 tetField: n x 3 x 3 numpy array 

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

519 

520 Returns 

521 ------- 

522 grainField: m x 3 x 3 

523 array containing (3x3) arrays of strain 

524 

525 Example 

526 ------- 

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

528 Returns strains for each grain. 

529 

530 Notes 

531 ------ 

532 Function contributed by Ryan Hurley (Johns Hopkins University) 

533 

534 """ 

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

536 # RCH Oct 2018. 

537 

538 # from progressbar import ProgressBar 

539 # pbar = ProgressBar() 

540 

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

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

543 # print("done.") 

544 

545 # Initialize list of grain values 

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

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

548 

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

550 # Loop through grains... 

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

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

553 # print("label = ", label) 

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

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

556 

557 volTot = 0.0 

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

559 # Loop through list of tets, summing the strains 

560 for touchingTet in touchingTets: 

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

562 vol = tetVolumes2[touchingTet] 

563 # Track total volume 

564 volTot += vol 

565 

566 # Add volume-weighted field 

567 fieldTot += tetField[touchingTet] * vol 

568 

569 # Divide strainTot by volTot 

570 fieldTot = fieldTot / volTot 

571 

572 # Store in particles 

573 grainField[label] = fieldTot 

574 

575 return grainField 

576 

577 

578def BCFieldFromDVCField( 

579 points, 

580 dvcField, 

581 mask=None, 

582 pixelSize=1.0, 

583 meshType="cube", 

584 centre=[0, 0], 

585 radius=1.0, 

586 topBottom=False, 

587 tol=1e-6, 

588 neighbours=4, 

589): 

590 """ 

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

592 

593 Parameters 

594 ---------- 

595 points: 2D numpy array 

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

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

598 

599 dvcField: 2D numpy array 

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

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

602 

603 mask: 2D numpy array, optional 

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

605 Points with 0 will be ignored for the field interpolation 

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

607 

608 pixelSize: float, optional 

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

610 Default = 1.0 

611 

612 meshType: string, optional 

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

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

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

616 Default = `cube` 

617 

618 centre: float, optional 

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

620 Default = [0, 0] 

621 

622 radius: float, optional 

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

624 Default = 1.0 

625 

626 topBottom: bool, optional 

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

628 Default = False 

629 

630 tol: float, optional 

631 Numerical tolerance for floats equality 

632 Default = 1e-6 

633 

634 neighbours: int, , optional 

635 Neighbours for field interpolation 

636 Default = 4 

637 Returns 

638 ------- 

639 bc: 2D numpy array 

640 Boundary node displacements 

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

642 

643 WARNING 

644 ------- 

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

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

647 """ 

648 

649 # STEP 1: find the edge nodes 

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

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

652 bcNodes = [] 

653 

654 if meshType == "cube": 

655 # extract edge nodes from the mesh 

656 for point in points: 

657 if topBottom: 

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

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

660 else: 

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

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

663 if testMin or testMax: 

664 bcNodes.append(point) 

665 elif meshType == "cylinder": 

666 # extract edge nodes from the mesh 

667 for point in points: 

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

669 testXY = False 

670 if not topBottom: 

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

672 if testZ or testXY: 

673 bcNodes.append(point) 

674 

675 bcNodes = numpy.array(bcNodes) 

676 m = bcNodes.shape[0] 

677 

678 # STEP 2: convert dvc field to physical unit 

679 dvcField *= pixelSize 

680 

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

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

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

684 if mask is None: 

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

686 

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

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

689 

690 for point in range(m): 

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

692 

693 # Check if we've hit the same point 

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

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

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

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

698 

699 else: 

700 weightSumInv = sum(1 / distance) 

701 # loop over neighbours 

702 for neighbour in range(neighbours): 

703 # calculate its weight 

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

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

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

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

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

709 

710 # return node number and displacements 

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

712 

713 

714def tetVolumes(points, connectivity): 

715 """ 

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

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

718 

719 Parameters 

720 ---------- 

721 points : Nx3 array 

722 Array of ``N`` coordinates of the points 

723 

724 connectivity : Mx4 array 

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

726 

727 Returns 

728 ------- 

729 volumes : vector of length M 

730 Volumes of tetrahedra 

731 

732 Note 

733 ----- 

734 Pure python function. 

735 """ 

736 

737 # Sanity checks on lengths: 

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

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

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

741 # return 

742 

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

744 connectivity = connectivity.astype(numpy.uint) 

745 

746 # loop through tetrahedra 

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

748 fourNodes = connectivity[tetNumber] 

749 

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

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

752 else: 

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

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

755 else: 

756 a = points[fourNodes[0]] 

757 b = points[fourNodes[1]] 

758 c = points[fourNodes[2]] 

759 d = points[fourNodes[3]] 

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

761 

762 return volumes 

763 

764 

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

766# import pygmsh 

767# import meshio 

768# 

769# lx, ly, lz = lengths 

770# ox, oy, oz = origin 

771# 

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

773# lengths = lengths[::-1] 

774# origin = origin[::-1] 

775# 

776# # raw code 

777# code = [] 

778# 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

816# 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)) 

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

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

819# 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)) 

820# 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)) 

821# 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)) 

822# 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)) 

823# 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)) 

824# 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)) 

825# 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)) 

826# 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)) 

827# 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)) 

828# 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)) 

829# 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)) 

830# 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)) 

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

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

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

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

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

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

837# 

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

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

840# geom = pygmsh.opencascade.Geometry() 

841# geom2 = pygmsh.opencascade.Geometry() 

842# 

843# # add raw code to geometry 

844# geom.add_raw_code(code) 

845# geom2.add_raw_code(code2) 

846# 

847# # mesh 

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

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

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

851# points = mesh.points 

852# cells = mesh.cells 

853# connectivity = cells['tetra'] 

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

855# pointsaux = meshaux.points 

856# cellsaux = meshaux.cells 

857# connectivityaux = cellsaux['tetra'] 

858# 

859# # write gmsh/vtk file 

860# if gmshFile is not None: 

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

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

863# if vtkFile is not None: 

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

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

866# 

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

868# # on veut produire : 

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

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

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

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

873# 

874# temppoints = copy.deepcopy(pointsaux) 

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

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

877# tempconnec = copy.deepcopy(connectivityaux) 

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

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

880# 

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

882# if gmshFile is not None: 

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

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

885# if vtkFile is not None: 

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

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

888# 

889# # now try to generate global mesh 

890# # its a three step process 

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

892# # 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. 

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

894# 

895# 

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

897# double = [] 

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

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

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

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

902# double.append(i+1) 

903# double.append(j+i+2) 

904# 

905# 

906# print('Start to build GLOBAL mesh - step 2') 

907# #here we should replace the nodes written in the double list in the glob_connectivity 

908# for node1,node2 in zip(double[0::2], double[1::2]): 

909# for k,elem in enumerate(glob_connectivity): 

910# if elem[0] == node2-1: 

911# print('Replacing double node ', node2, ' by ', node1, 'in element ', k+1, 'component ', 1) 

912# glob_connectivity[k][0] = node1-1 

913# if elem[1] == node2-1: 

914# print('Replacing double node ', node2, ' by ', node1, 'in element ', k+1, 'component ', 2) 

915# glob_connectivity[k][1] = node1-1 

916# if elem[2] == node2-1: 

917# print('Replacing double node ', node2, ' by ', node1, 'in element ', k+1, 'component ', 3) 

918# glob_connectivity[k][2] = node1-1 

919# if elem[3] == node2-1: 

920# print('Replacing double node ', node2, ' by ', node1, 'in element ', k+1, 'component ', 4) 

921# glob_connectivity[k][3] = node1-1 

922# 

923# print('Start to build GLOBAL mesh - step 3') 

924# #here we should erase the double nodes written in the node list and shift -1 the glob_connectivity 

925# toberemoved = [] 

926# for node1,node2 in zip(double[0::2], double[1::2]): 

927# if len(toberemoved) == 0: 

928# toberemoved.append(node2) 

929# elif toberemoved[-1] != node2: 

930# toberemoved.append(node2) 

931# toberemoved.sort(reverse=True) 

932# # print('toberemoved : ', toberemoved) 

933# for node in toberemoved: 

934# glob_points = numpy.delete(glob_points, node-1,0) # Point removing 

935# for k,elem in enumerate(glob_connectivity): 

936# if elem[0] > node-1: 

937# glob_connectivity[k][0] -= 1 

938# if elem[1] > node-1: 

939# glob_connectivity[k][1] -= 1 

940# if elem[2] > node-1: 

941# glob_connectivity[k][2] -= 1 

942# if elem[3] > node-1: 

943# glob_connectivity[k][3] -= 1 

944# 

945# meshio.write_points_cells("global.msh", glob_points, {'tetra': glob_connectivity}, file_format='gmsh22') 

946# meshio.write_points_cells("global.vtk", glob_points, {'tetra': glob_connectivity}, file_format='vtk') 

947# 

948# 

949# return 

950# 

951# 

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

953# # NOT USED 

954# import pygmsh 

955# import meshio 

956# 

957# lx, ly, lz = lengths 

958# ox, oy, oz = origin 

959# 

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

961# lengths = lengths[::-1] 

962# origin = origin[::-1] 

963# 

964# # raw code 

965# code = [] 

966# 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1004# 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)) 

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

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

1007# 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)) 

1008# 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)) 

1009# 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)) 

1010# 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)) 

1011# 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)) 

1012# 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)) 

1013# 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)) 

1014# 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)) 

1015# 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)) 

1016# 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))