pydfnWorks
python wrapper for dfnWorks
map2continuum.py
Go to the documentation of this file.
1 """
2 .. module:: map2continuum.py
3  :synopsis: Produce octree-refined continuum mesh replacing dfn
4 .. moduleauthor:: Matthew Sweeney <msweeney2796@lanl.gov>
5 
6 """
7 
8 import os
9 import sys
10 import subprocess
11 import shutil
12 import numpy as np
13 from pydfnworks.dfnGen.meshing import mesh_dfn_helper as mh
14 import time
15 import multiprocessing as mp
16 #import Queue
17 
18 
19 def map_to_continuum(self, l, orl, path="./", dir_name="octree"):
20  """ This function generates an octree-refined continuum mesh using the
21  reduced_mesh.inp as input. To generate the reduced_mesh.inp, one must
22  turn visualization mode on in the DFN input card.
23 
24  Parameters
25  ----------
26  self : object
27  DFN Class
28  l : float
29  Size (m) of level-0 mesh element in the continuum mesh
30  orl : int
31  Number of total refinement levels in the octree
32  path : string
33  path to primary DFN directory
34  dir_name : string
35  name of directory where the octree mesh is created
36 
37  Returns
38  -------
39  None
40 
41  Notes
42  -----
43  octree_dfn.inp : Mesh file
44  Octree-refined continuum mesh
45  fracX.inp : Mesh files
46  Octree-refined continuum meshes, which contain intersection areas
47 
48  """
49  print('=' * 80)
50  print("Meshing Continuum Using LaGrit : Starting")
51  print('=' * 80)
52 
53  if type(orl) is not int or orl < 1:
54  error = "ERROR: orl must be positive integer. Exiting"
55  sys.stderr.write(error)
56  sys.exit(1)
57 
58  # Read in normal vectors and points from dfnWorks output
59  normal_vectors = np.genfromtxt(path + 'normal_vectors.dat', delimiter=' ')
60 
61  with open(path + 'translations.dat') as old, open('points.dat',
62  'w') as new:
63  old.readline()
64  for line in old:
65  if not 'R' in line:
66  new.write(line)
67  points = np.genfromtxt('points.dat', skip_header=0, delimiter=' ')
68 
69  try:
70  os.symlink(path + 'params.txt', 'params.txt')
71  except:
72  pass
73  num_poly, h, _, _, domain = mh.parse_params_file()
74 
75  # Extent of domain
76  x0 = 0 - (domain['x'] / 2.0)
77  x1 = 0 + (domain['x'] / 2.0)
78  y0 = 0 - (domain['y'] / 2.0)
79  y1 = 0 + (domain['y'] / 2.0)
80  z0 = 0 - (domain['z'] / 2.0)
81  z1 = 0 + (domain['z'] / 2.0)
82 
83  # Number of cell elements in each direction at coarse level
84  nx = domain['x'] / l + 1
85  ny = domain['y'] / l + 1
86  nz = domain['z'] / l + 1
87 
88  if nx * ny * nz > 1e8:
89  error = "ERROR: Number of elements > 1e8. Exiting"
90  sys.stderr.write(error)
91  sys.exit(1)
92 
93  print("\nCreating *.lgi files for octree mesh\n")
94  try:
95  os.mkdir(dir_name)
96  except OSError:
97  shutil.rmtree(dir_name)
98  os.mkdir(dir_name)
99 
100  lagrit_driver(dir_name, nx, ny, nz, num_poly, normal_vectors, points)
101  lagrit_parameters(dir_name, orl, x0, x1, y0, y1, z0, z1, nx, ny, nz, h)
102  lagrit_build(dir_name)
103  lagrit_intersect(dir_name)
104  lagrit_hex_to_tet(dir_name)
105  lagrit_remove(dir_name)
106  lagrit_run(path, dir_name)
107  lagrit_strip(num_poly)
108  driver_parallel(self, num_poly)
109 
110 
111 def lagrit_driver(dir_name, nx, ny, nz, num_poly, normal_vectors, points):
112  """ This function creates the main lagrit driver script, which calls all
113  lagrit scripts.
114 
115  Parameters
116  ----------
117  dir_name : string
118  Name of working directory
119  ni : int
120  Number of cells in each direction
121  num_poly : int
122  Number of fractures
123  normal_vectors : array
124  Array containing normal vectors of each fracture
125  points : array
126  Array containing a point on each fracture
127 
128  Returns
129  -------
130  None
131 
132  Notes
133  -----
134  None
135 
136  """
137  xn, yn, zn, xp, yp, zp = 0, 0, 0, 0, 0, 0
138  j = num_poly + 1
139  floop = ""
140  for i in range(1, int(num_poly + 2)):
141  if j != num_poly + 1:
142  xn = normal_vectors[j - 1][0]
143  yn = normal_vectors[j - 1][1]
144  zn = normal_vectors[j - 1][2]
145  xp = points[j - 1][0]
146  yp = points[j - 1][1]
147  zp = points[j - 1][2]
148  floop += """cmo / create / FRACTURE{0}
149  cmo / copy / FRACTURE{0} / MODFN
150  cmo / select / FRACTURE{0}
151  pset / pdel / attribute imt / 1 0 0 / ne / {0}
152  rmpoint pset, get, pdel
153  rmpoint compress
154  resetpts itp
155  pset / pdel / delete
156  eltset / edel / itetclr / ne / {0}
157  rmpoint element eltset, get, edel
158  rmpoint compress
159  resetpts itp
160  eltset / edel / delete
161 
162  intersect_elements / mohex2 / FRACTURE{0} / if_int
163  eltset / erefine / if_int / ne / 0
164  pset / pinter / eltset erefine
165 
166  cmo / create / temp{0}
167  cmo / copy / temp{0} / mohex2
168  cmo / select / temp{0}
169  pset / pdel / if_int / 1 0 0 / eq / 0
170  rmpoint, pset, get, pdel
171  rmpoint compress
172  resetpts itp
173  pset / pdel / delete
174  eltset / edel / if_int / eq / 0
175  rmpoint element eltset, get, edel
176  rmpoint compress
177  resetpts itp
178  eltset / edel / delete
179  cmo / setatt / temp{0} / itetclr / 1 0 0 / {0}
180  cmo / setatt / temp{0} / imt / 1 0 0 / {0}
181 
182  cmo / create / TETCOPY
183  cmo / copy / TETCOPY / MOTET
184  cmo / select / TETCOPY
185  interpolate / map / TETCOPY / itetclr / 1 0 0 / temp{0} / itetclr
186  compute / distance_field / TETCOPY / temp{0} /dfield
187  cmo / select / TETCOPY
188  pset / pfrac / attribute / dfield / 1 0 0 / le / 1.e-8
189  cmo / setatt / TETCOPY / imt / 1 0 0 / {1}
190  cmo / setatt / TETCOPY / imt / pset get pfrac / {0}
191  eltset / etemp / itetclr / ne / {0}
192  cmo / setatt / TETCOPY / itetclr / eltset get etemp / {1}
193  eltset / etemp / delete
194  pset / pfrac / delete
195 
196  cmo / set_id / MOTET_np1 / element / id_cell
197  cmo / set_id / MOTET_np1 / node / id_vertex
198 
199  extract / plane / ptnorm / &
200  {2} {3} {4} / &
201  {5} {6} {7} / &
202  1 0 0 / moext / MOTET_np1
203  cmo / status / brief
204 
205  createpts / median
206  cmo / DELATT / moext / id_cell
207  cmo / DELATT / moext / id_parent
208  dump / avs2 / ex_xyz{0}.table / moext 0 0 0 1
209 
210  cmo / addatt / moext / volume / area_tri
211  cmo / DELATT / moext / xmed
212  cmo / DELATT / moext / ymed
213  cmo / DELATT / moext / zmed
214  dump / avs2 / ex_area{0}.table / moext 0 0 0 1
215  cmo / delete / moext
216 
217  cmo / select / TETCOPY
218  cmo / DELATT / TETCOPY / icr1
219  cmo / DELATT / TETCOPY / itp1
220  cmo / DELATT / TETCOPY / isn1
221 
222  dump / avs2 / frac{0}.inp / TETCOPY
223  cmo / delete / temp{0}
224  cmo / delete / TETCOPY
225 
226  cmo / select / mohex2
227  cmo / setatt / mohex2 / itetclr / eltset get erefine / {0}
228  cmo / setatt / mohex2 / imt / pset get pinter / {0}
229  eltset / erefine / delete
230  pset / pinter / delete
231  cmo / DELATT / mohex2 / if_int
232  cmo / delete / FRACTURE{0}
233  """.format(j, num_poly + 1, xp, yp, zp, xn, yn, zn)
234  j = j - 1
235 
236  f_name = f'{dir_name}/driver_octree.lgi'
237  f = open(f_name, 'w')
238  fin = ("""#
239 # LaGriT control files to build an octree refined hex mesh with refinement
240 # based on intersection of hex mesh with a DFN triangulation mesh
241 #
242 # driver_octree.lgi
243 # parameters_octree_dfn.mlgi
244 # build_octree.mlgi
245 # intersect_refine.mlgi
246 # hex_to_tet.mlgi
247 # remove_cells.mlgi
248 #
249 # Define some parameters
250 #
251 infile parameters_octree_dfn.mlgi
252 #
253 # Read in DFN mesh
254 #
255 read / FTYPE / FNAME / MODFN
256 cmo / printatt / MODFN / -xyz- / minmax
257 #
258 # Octree refined orthogonal mesh based on intersection with DFN mesh
259 #
260 infile build_octree.mlgi
261 #
262 # Identify cells in hex mesh that are intersected by DFN mesh
263 #
264 # This is the last pass through intersect_elements in order to figure out
265 # which cells in the fully refined hex mesh are intersected by the dfn mesh
266 #
267 intersect_elements / MOHEX / MODFN / if_int
268 eltset / einter / if_int / ne / 0
269 pset / pinter / eltset einter
270 #
271 # Use the itetclr(cell) and imt(vertex) attribute to hold the information
272 #
273 cmo / setatt / MOHEX / itetclr / 1 0 0 / 1
274 cmo / setatt / MOHEX / itetclr / eltset get einter / 2
275 cmo / setatt / MOHEX / imt / 1 0 0 / 1
276 cmo / setatt / MOHEX / imt / pset get pinter / 2
277 #
278 # Output final hex mesh
279 #
280 #dump / avs2 / tmp_hex_refine.inp / MOHEX
281 #
282 # Same as above but for np1 hex mesh
283 #
284 intersect_elements / MOHEX_np1 / MODFN / if_int
285 eltset / einter / if_int / ne / 0
286 pset / pinter / eltset einter
287 #
288 # See above
289 #
290 cmo / setatt / MOHEX_np1 / itetclr / 1 0 0 / 1
291 cmo / setatt / MOHEX_np1 / itetclr / eltset get einter / 2
292 cmo / setatt / MOHEX_np1 / imt / 1 0 0 / 1
293 cmo / setatt / MOHEX_np1 / imt / pset get pinter / 2
294 #dump / avs2 / tmp_hex_np1_refine.inp / MOHEX_np1
295 #
296 # Convert the hex mesh to a tet mesh
297 #
298 infile hex_to_tet.mlgi
299 #
300 # Modify the hex data structure from a full octree data structure
301 # to one in which only the highest level of refined hex is maintained
302 # and all parent cells are stripped out of the data structure
303 #
304 grid2grid / tree_to_fe / mohex2 / MOHEX
305 #dump / avs / octree_hex_mesh.inp / MOHEX
306 #
307 cmo / delete / MOHEX
308 cmo / select / mohex2
309 #
310 # Remove all but the most refined hex cells
311 #
312 loop / do / NTIMES / 0 N_OCTREE_REFINE_M1 1 / loop_end &
313 infile remove_cells.mlgi
314 #
315 cmo / select / mohex2
316 cmo / DELATT / mohex2 / if_int
317 intersect_elements / mohex2 / MODFN / if_int
318 cmo / select / mohex2
319 eltset / edelete / if_int / eq / 0
320 rmpoint / element / eltset get edelete
321 eltset / edelete / release
322 rmpoint / compress
323 #
324 # NOTE: I commented out the following lines for unknown but seemingly
325 # important reason
326 #
327 #cmo / setatt / mohex2 / itetclr / 1 0 0 / 1
328 #cmo / setatt / mohex2 / imt / 1 0 0 / 1
329 """ + floop + """
330 #
331 #dump / avs / tmp_remove_cells.inp / mohex2
332 #dump / avs / tmp_tet.inp / MOTET
333 #
334 interpolate / map / MOTET / itetclr / 1 0 0 / mohex2 / itetclr
335 compute / distance_field / MOTET / mohex2 / dfield
336 cmo / select / MOTET
337 #
338 # This use of 1.e-8 is fine for now but might cause problems if the
339 # bounding box of the problem was very small (<1.e-7)
340 #
341 pset / pfrac / attribute / dfield / 1 0 0 / le / 1.e-8
342 cmo / setatt / MOTET / imt / 1 0 0 / 2
343 cmo / setatt / MOTET / imt / pset get pfrac / 1
344 #
345 cmo / modatt / MOTET / itp / ioflag / l
346 cmo / modatt / MOTET / isn / ioflag / l
347 cmo / modatt / MOTET / icr / ioflag / l
348 #
349 dump / pflotran / full_mesh / MOTET / nofilter_zero
350 dump / avs2 / octree_dfn.inp / MOTET
351 dump / coord / octree_dfn / MOTET
352 dump / stor / octree_dfn / MOTET
353 dump / zone_imt / octree_dfn / MOTET
354 dump / zone_outside / octree_dfn / MOTET
355 
356 define / ZONE / 1
357 define / FOUT / pboundary_top
358 pset / top / attribute / zic / 1 0 0 / gt / ZMAX
359 pset / top / zone / FOUT / ascii / ZONE
360 
361 define / ZONE / 2
362 define / FOUT / pboundary_bottom
363 pset / bottom / attribute / zic / 1 0 0 / lt / ZMIN
364 pset / bottom / zone / FOUT / ascii / ZONE
365 
366 define / ZONE / 3
367 define / FOUT / pboundary_left_w
368 pset / left_w / attribute / xic / 1 0 0 / lt / XMIN
369 pset / left_w / zone / FOUT / ascii / ZONE
370 
371 define / ZONE / 4
372 define / FOUT / pboundary_front_n
373 pset / front_n / attribute / yic / 1 0 0 / gt / YMAX
374 pset / front_n / zone / FOUT / ascii / ZONE
375 
376 define / ZONE / 5
377 define / FOUT / pboundary_right_e
378 pset / right_e / attribute / xic / 1 0 0 / gt / XMAX
379 pset / right_e / zone / FOUT / ascii / ZONE
380 
381 define / ZONE / 6
382 define / FOUT / pboundary_back_s
383 pset / back_s / attribute / yic / 1 0 0 / lt / YMIN
384 pset / back_s / zone / FOUT / ascii / ZONE
385 
386 #
387 # Work around for getting *.fehnm file
388 # Do we need this?
389 #
390 cmo / setatt / MOTET / itetclr / 1 0 0 / 1
391 cmo / setatt / MOTET / imt / 1 0 0 / 1
392 resetpts / itp
393 dump / fehm / tmp_tmp_ / MOTET
394 #
395 finish
396 """)
397  f.write(fin)
398  f.flush()
399  f.close()
400  print("Creating driver_octree.lgi file: Complete\n")
401 
402 
403 def lagrit_parameters(dir_name, orl, x0, x1, y0, y1, z0, z1, nx, ny, nz, h):
404  """ This function creates the parameters_octree_dfn.mlgi lagrit script.
405 
406  Parameters
407  ----------
408  dir_name : string
409  Name of working directory
410  orl : int
411  Number of total refinement levels in the octree
412  i0, i1 : float
413  Extent of domain in x, y, z directions
414  ni : int
415  Number of cells in each direction
416  Returns
417  -------
418  None
419 
420  Notes
421  -----
422  None
423 
424  """
425  f_name = f'{dir_name}/parameters_octree_dfn.mlgi'
426  f = open(f_name, 'w')
427  fin = """#
428 # Define some parameters
429 #
430 # Input DFN mesh
431 #
432 define / FTYPE / avs
433 define / FNAME / reduced_mesh.inp
434 #
435 define / MOHEX / mohex
436 define / MOTET / motet
437 #
438 # Set AMR refinement. 123 means refine in x,y,z directions
439 # See LaGriT refine command for more options
440 #
441 define / REFINE_AMR / 123
442 #
443  """
444  f.write(fin)
445  eps = h * 10**-3
446  f.write('define / N_OCTREE_REFINE / %d \n' % orl)
447  f.write('define / N_OCTREE_REFINE_M1 / %d \n' % (orl - 1))
448  f.write('define / N_OCTREE_np1 / %d \n' % (orl + 1))
449  f.write('define / N_OCTREE_np1_M1 / %d \n' % (orl))
450  f.write('define / X0 / %0.12f \n' % x0)
451  f.write('define / X1 / %0.12f \n' % x1)
452  f.write('define / Y0 / %0.12f \n' % y0)
453  f.write('define / Y1 / %0.12f \n' % y1)
454  f.write('define / Z0 / %0.12f \n' % z0)
455  f.write('define / Z1 / %0.12f \n' % z1)
456  f.write('define / XMAX / %0.12f \n' % (x1 - eps))
457  f.write('define / XMIN / %0.12f \n' % (x0 + eps))
458  f.write('define / YMAX / %0.12f \n' % (y1 - eps))
459  f.write('define / YMIN / %0.12f \n' % (y0 + eps))
460  f.write('define / ZMAX / %0.12f \n' % (z1 - eps))
461  f.write('define / ZMIN / %0.12f \n' % (z0 + eps))
462  f.write('define / NX / %d \n' % nx)
463  f.write('define / NY / %d \n' % ny)
464  f.write('define / NZ / %d \n' % nz)
465  f.write('finish\n')
466  f.flush()
467  f.close()
468  print("Creating parameters_octree_dfn.mlgi file: Complete\n")
469 
470 
471 def lagrit_build(dir_name):
472  """ This function creates the build_octree.mlgi lagrit script.
473 
474  Parameters
475  ----------
476  dir_name : string
477  name of working directory
478 
479  Returns
480  -------
481  None
482 
483  Notes
484  -----
485  None
486 
487  """
488  f_name = f'{dir_name}/build_octree.mlgi'
489  f = open(f_name, 'w')
490  fin = """cmo / create / MOHEX / / / hex
491 createpts / brick / xyz / NX NY NZ / X0 Y0 Z0 / X1 Y1 Z1 / 1 1 1
492 cmo / setatt / MOHEX / imt / 1 0 0 / 1
493 cmo / setatt / MOHEX / itetclr / 1 0 0 / 1
494 resetpts / itp
495 #
496 # Print to screen and logfiles the extents of the hex mesh and dfn mesh
497 #
498 cmo / printatt / MOHEX / -xyz- / minmax
499 cmo / select / MODFN
500 cmo / printatt / MODFN / -xyz- / minmax
501 #
502 # Generate copy of hex mesh for upscaling
503 #
504 cmo / copy / MOHEX_np1 / MOHEX
505 #
506 # Loop through steps to intersect dfn mesh (MODFN) with hex mesh (MOHEX)
507 # and refine hex mesh based on cell intersections. Loop through
508 # N_OCTREE_REFINE times
509 #
510 loop / do / NTIMES / 1 N_OCTREE_REFINE 1 / loop_end &
511 infile intersect_refine.mlgi
512 #
513 # See above - except do it once additional time ("np1")
514 #
515 loop / do / NTIMEs / 1 N_OCTREE_np1 1 / loop_end &
516 infile intersect_refine_np1.mlgi
517 #
518 finish
519  """
520  f.write(fin)
521  f.flush()
522  f.close()
523  print("Creating build_octree.mlgi file: Complete\n")
524 
525 
526 def lagrit_intersect(dir_name):
527  """ This function creates the intersect_refine.mlgi lagrit scripts.
528 
529  Parameters
530  ----------
531  dir_name : string
532  name of working directory
533 
534  Returns
535  -------
536  None
537 
538  Notes
539  -----
540  None
541 
542  """
543  f_name = f'{dir_name}/intersect_refine.mlgi'
544  f = open(f_name, 'w')
545  fin = """#
546 # Compute mesh to mesh intersection and refine hex mesh
547 #
548 intersect_elements / MOHEX / MODFN / if_int
549 eltset / erefine / if_int / ne / 0
550 pset / prefine / eltset erefine
551 #
552 # If one wants anisotropic mesh refinement, then the
553 # variable REFINE_AMR can be set in parameters_octree_dfn.mlgi
554 #
555 refine/constant/imt1/linear/element/pset,get,prefine/ &
556  -1.,0.,0./exclusive/amr REFINE_AMR
557 #
558 # Clean up eltset, pset, and if_int attribute
559 #
560 eltset / erefine / delete
561 pset / prefine / delete
562 cmo / DELATT / MOHEX / if_int
563 #
564 # Print out diagnostics
565 #
566 quality
567 cmo / status / brief
568 #
569 finish
570  """
571  f.write(fin)
572  f.flush()
573  f.close()
574  print("Creating intersect_refine.mlgi file: Complete\n")
575  f_name = f'{dir_name}/intersect_refine_np1.mlgi'
576  f = open(f_name, 'w')
577  fin = """#
578 # For comments see intersect_refine.mlgi
579 #
580 intersect_elements / MOHEX_np1 / MODFN / if_int
581 eltset / erefine / if_int / ne / 0
582 pset / prefine / eltset erefine
583 refine/constant/imt1/linear/element/pset,get,prefine/ &
584 -1.,0.,0./exclusive/amr REFINE_AMR
585 #
586 eltset / erefine / delete
587 pset / prefine / delete
588 cmo / DELATT / MOHEX_np1 / if_int
589 #
590 quality
591 cmo / status / brief
592 #
593 finish
594  """
595  f.write(fin)
596  f.flush()
597  f.close()
598  print("Creating intersect_refine_np1.mlgi file: Complete\n")
599 
600 
601 def lagrit_hex_to_tet(dir_name):
602  """ This function creates the hex_to_tet.mlgi lagrit script.
603 
604  Parameters
605  ----------
606  dir_name : string
607  name of working directory
608 
609  Returns
610  -------
611  None
612 
613  Notes
614  -----
615  None
616 
617  """
618  f_name = f'{dir_name}/hex_to_tet.mlgi'
619  f = open(f_name, 'w')
620  fin = """#
621 # Convert the octree hex mesh to tet mesh by connecting the
622 # octree vertex collection to a Delaunay mesh
623 #
624 cmo / create / MOTET
625 copypts / MOTET / MOHEX
626 cmo / select / MOTET
627 filter 1 0 0
628 rmpoint / compress
629 cmo / setatt / MOTET / imt / 1 0 0 / 1
630 cmo / setatt / MOTET / itp / 1 0 0 / 0
631 #
632 # Sort and reorder the nodes based on coordinates
633 #
634 sort / MOTET / index / ascending / ikey / xic yic zic
635 reorder / MOTET / ikey
636 cmo / DELATT / MOTET / ikey
637 #
638 # Connect
639 #
640 connect / noadd
641 cmo / setatt / MOTET / itetclr / 1 0 0 / 1
642 resetpts / itp
643 #
644 # Do the same for np1 mesh
645 #
646 cmo / create / MOTET_np1
647 copypts / MOTET_np1 / MOHEX_np1
648 cmo / select / MOTET_np1
649 filter 1 0 0
650 rmpoint / compress
651 cmo / setatt / MOTET_np1 / imt / 1 0 0 / 1
652 cmo / setatt / MOTET_np1 / itp / 1 0 0 / 0
653 #
654 # See above
655 #
656 sort / MOTET_np1 / index / ascending / ikey / xic yic zic
657 reorder / MOTET_np1 / ikey
658 cmo / DELATT / MOTET_np1 / ikey
659 #
660 connect / noadd
661 cmo / setatt / MOTET_np1 / itetclr / 1 0 0 / 1
662 resetpts / itp
663 #
664 finish
665  """
666  f.write(fin)
667  f.flush()
668  f.close()
669  print("Creating hex_to_tet.mlgi file: Complete\n")
670 
671 
672 def lagrit_remove(dir_name):
673  """ This function creates the remove_cells.mlgi lagrit script.
674 
675  Parameters
676  ----------
677  dir_name : string
678  name of working directory
679 
680  Returns
681  -------
682  None
683 
684  Notes
685  -----
686  None
687 
688  """
689  f_name = f'{dir_name}/remove_cells.mlgi'
690  f = open(f_name, 'w')
691  fin = """#
692 # Remove cells from hex mesh based on the level of refinement
693 # itetlev is the refinement level. Original mesh itetlev=0
694 #
695 eltset / edelete / itetlev / eq / NTIMES
696 rmpoint / element / eltset get edelete
697 eltset / edelete / release
698 rmpoint / compress
699 #
700 finish
701  """
702  f.write(fin)
703  f.flush()
704  f.close()
705  print("Creating remove_cells.mlgi file: Complete\n")
706 
707 
708 def lagrit_run(path, dir_name):
709  """ This function executes the lagrit scripts.
710 
711  Parameters
712  ----------
713  None
714 
715  Returns
716  -------
717  None
718 
719  Notes
720  -----
721  None
722 
723  """
724  # Run LaGriT
725  os.chdir(dir_name)
726 
727  if os.path.isfile(path + "reduced_mesh.inp"):
728  os.symlink(path + "reduced_mesh.inp", "reduced_mesh.inp")
729  elif os.path.isfile(path + "../" + "reduced_mesh.inp"):
730  os.symlink(path + "../" + "reduced_mesh.inp", "reduced_mesh.inp")
731  else:
732  error = "ERROR!!! Reduced Mesh not found. Please run mesh_dfn with visual_mode=True.\nExiting"
733  sys.stderr.write(error)
734  sys.exit(1)
735 
736  mh.run_lagrit_script("driver_octree.lgi")
737  # cmd = os.environ['LAGRIT_EXE'] + '< driver_octree.lgi'
738  # subprocess.call(cmd, shell=True)
739 
740 
741 def lagrit_strip(num_poly):
742  """ This function strips and replaces the headers of the files, which is
743  needed to assign the fracture areas to a mesh object.
744 
745  Parameters
746  ----------
747  num_poly : int
748  Number of fractures
749 
750  Returns
751  -------
752  None
753 
754  Notes
755  -----
756  None
757 
758  """
759  node_dict = {}
760  for i in range(1, num_poly + 1):
761  with open('ex_xyz{0}.table'.format(i), 'r') as infile:
762  for k, l in enumerate(infile):
763  pass
764  node_dict.setdefault(i, []).append(k - 4)
765  print(node_dict)
766 
767  for i in range(1, num_poly + 1):
768  with open('ex_xyz{0}.table'.format(i), 'r') as infile:
769  with open('ex_xyz{0}_2.inp'.format(i), 'w') as outfile:
770  for j, line in enumerate(infile):
771  if j == 0:
772  outfile.write("{0} 0 0 0 0\n".format(node_dict[i][0]))
773  elif j > 4:
774  outfile.write(line)
775  outfile.close()
776  infile.close()
777  with open('ex_area{0}.table'.format(i), 'r') as infile:
778  with open('ex_area{0}_2.table'.format(i), 'w') as outfile:
779  for j, line in enumerate(infile):
780  if j > 2:
781  outfile.write(line.split()[1])
782  outfile.write("\n")
783  outfile.close()
784  infile.close()
785 
786 
787 def driver_parallel(self, num_poly):
788  """ This function drives the parallelization of the area sums upscaling.
789 
790  Parameters
791  ----------
792  self : object
793  DFN Class
794  num_poly : int
795  Number of fractures
796 
797  Returns
798  -------
799  None
800 
801  Notes
802  -----
803  None
804 
805  """
806  frac_index = range(1, int(num_poly + 1))
807  number_of_task = len(frac_index)
808  number_of_processes = self.ncpu
809  tasks_to_accomplish = mp.Queue()
810  tasks_that_are_done = mp.Queue()
811  processes = []
812 
813  for i in range(number_of_task):
814  tasks_to_accomplish.put(i + 1)
815 
816  # Creating processes
817  for w in range(number_of_processes):
818  p = mp.Process(target=worker,
819  args=(tasks_to_accomplish, tasks_that_are_done))
820  processes.append(p)
821  p.start()
822  tasks_to_accomplish.put('STOP')
823 
824  for p in processes:
825  p.join()
826 
827  while not tasks_that_are_done.empty():
828  print(tasks_that_are_done.get())
829 
830  return True
831 
832 
833 def upscale_parallel(f_id):
834  """ Generates lagrit script that makes mesh files with area sums.
835 
836  Parameters
837  ----------
838  f_id : int
839  Fracture index
840 
841  Returns
842  -------
843  None
844 
845  Notes
846  -----
847  None
848 
849  """
850  fname = 'driver{0}.lgi'.format(f_id)
851  fin = ""
852  f = open(fname, 'w')
853  fin += """read / avs / ex_xyz{0}_2.inp / mo_vertex
854  cmo / addatt / mo_vertex / area_tri / vdouble / scalar / nnodes
855  cmo / readatt / mo_vertex / area_tri / 1 0 0 / ex_area{0}_2.table
856 
857  read / avs / frac{0}.inp / frac
858  cmo / addatt / frac / area_sum / vdouble / scalar / nnodes
859 
860  upscale / sum / frac, area_sum / 1 0 0 / mo_vertex, area_tri
861  dump / avs / area_sum{0}.inp / frac
862 
863  cmo / delete / mo_vertex
864  cmo / delete / frac
865  """.format(f_id)
866  fin += "finish"
867  print(fin)
868  f.write(fin)
869  f.flush()
870  f.close()
871 
872  # cmd = os.environ['LAGRIT_EXE'] + '< driver{0}.lgi'.format(f_id)
873  # subprocess.call(cmd, shell=True)
874  mh.run_lagrit_script("driver{0}.lgi".format(f_id))
875 
876 
877 def worker(tasks_to_accomplish, tasks_that_are_done):
878  """ Worker function for python parallel. See multiprocessing module
879  documentation for details.
880 
881  Parameters
882  ----------
883  tasks_to_accomplish : ?
884  Processes still in queue
885  tasks_that_are_done : ?
886  Processes complete
887 
888  Notes
889  -----
890  None
891 
892  """
893  try:
894  for f_id in iter(tasks_to_accomplish.get, 'STOP'):
895  upscale_parallel(f_id)
896  except:
897  pass
898  return True
pydfnworks.dfnGen.meshing.udfm.map2continuum.lagrit_parameters
def lagrit_parameters(dir_name, orl, x0, x1, y0, y1, z0, z1, nx, ny, nz, h)
Definition: map2continuum.py:403
pydfnworks.dfnGen.meshing.udfm.map2continuum.worker
def worker(tasks_to_accomplish, tasks_that_are_done)
Definition: map2continuum.py:877
pydfnworks.dfnGen.meshing.udfm.map2continuum.lagrit_driver
def lagrit_driver(dir_name, nx, ny, nz, num_poly, normal_vectors, points)
Definition: map2continuum.py:111
pydfnworks.dfnGen.meshing.udfm.map2continuum.map_to_continuum
def map_to_continuum(self, l, orl, path="./", dir_name="octree")
Definition: map2continuum.py:19