summaryrefslogtreecommitdiff
path: root/PyGan/data/rDragView_proc
diff options
context:
space:
mode:
Diffstat (limited to 'PyGan/data/rDragView_proc')
-rw-r--r--PyGan/data/rDragView_proc/3DHEX.c2m222
-rw-r--r--PyGan/data/rDragView_proc/README.md30
-rw-r--r--PyGan/data/rDragView_proc/get_lcm.py105
-rw-r--r--PyGan/data/rDragView_proc/utils.py56
-rw-r--r--PyGan/data/rDragView_proc/vtu_tools.py463
5 files changed, 876 insertions, 0 deletions
diff --git a/PyGan/data/rDragView_proc/3DHEX.c2m b/PyGan/data/rDragView_proc/3DHEX.c2m
new file mode 100644
index 0000000..52cb024
--- /dev/null
+++ b/PyGan/data/rDragView_proc/3DHEX.c2m
@@ -0,0 +1,222 @@
+*-----------------------------------------------
+* Benchmark 3D TAKEDA MODEL 4 CR FULLY INSERTED
+*-----------------------------------------------
+PARAMETER Geom Macr Trck Syst Flux ::
+ ::: LINKED_LIST Geom Macr Trck Syst Flux ; ;
+*
+MODULE GEO: SNT: MAC: ASM: FLU: OUT: GREP: END: ;
+*
+SEQ_ASCII OutGeom :: FILE './_GEOM' ;
+SEQ_ASCII OutMacr :: FILE './_MACR' ;
+SEQ_ASCII OutTrck :: FILE './_TRCK' ;
+SEQ_ASCII OutSyst :: FILE './_ASMR' ;
+SEQ_ASCII OutFlux :: FILE './_FLUX' ;
+*
+*
+* spl = region split
+* sch = scheme parameter
+* spa = spatial order
+* ang = sn order
+* sca = scattering
+*
+*
+INTEGER spl := 2 ;
+INTEGER sch := 1 ;
+INTEGER spa := 0 ;
+INTEGER ang := 4 ;
+INTEGER sca := 1 ;
+*
+*
+Geom := GEO: :: HEXZ 169 5
+ EDIT 0
+ HBC COMPLETE VOID
+ SIDE 7.50
+ Z- VOID Z+ VOID
+ MESHZ 0.0 45.0 65.0 125.0 145.0 190.0
+ SPLITZ 2 1 3 1 2
+ SPLITL <<spl>>
+ MIX
+
+ (* LAYER 1: STEEL *)
+ 1
+ 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
+ 7 7 7 7 7 7
+ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
+ 8 8 8 8 8 8 8 8 8 8 8 8
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10
+ (* LAYER 2: AXIAL BLANKET AND REFLECTOR *)
+ 2
+ 2 2 2 2 2 2
+ 3 1 3 1 3 1 3 1 3 1 3 1
+ 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
+ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
+ 7 7 7 7 7 7
+ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
+ 8 8 8 8 8 8 8 8 8 8 8 8
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10
+ (* LAYER 3: TEST ZONE *)
+ 4
+ 4 4 4 4 4 4
+ 5 12 5 12 5 12 5 12 5 12 5 12
+ 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
+ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
+ 7 7 7 7 7 7
+ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
+ 8 8 8 8 8 8 8 8 8 8 8 8
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10
+ (* LAYER 4: AXIAL BLANKET AND REFLECTOR *)
+ 2
+ 2 2 2 2 2 2
+ 3 1 3 1 3 1 3 1 3 1 3 1
+ 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
+ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
+ 7 7 7 7 7 7
+ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
+ 8 8 8 8 8 8 8 8 8 8 8 8
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10
+ (* LAYER 5: STEEL *)
+ 1
+ 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+ 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
+ 7 7 7 7 7 7
+ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
+ 8 8 8 8 8 8 8 8 8 8 8 8
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
+ 10 10 10 10 10 10
+ ;
+OutGeom := Geom ;
+*
+*
+Macr := MAC: ::
+ EDIT 0 NGRO 4 NMIX 12 NIFI 1
+ READ INPUT
+ MIX 1 (* STEEL *)
+ TOTAL 9.83638E-2 1.35140E-1 2.24749E-1 2.83117E-1
+ SCAT 1 1 9.06050E-2
+ 2 2 1.30581E-1 7.42377E-3
+ 3 3 2.19547E-1 4.35250E-3 1.18163E-4
+ 4 4 2.80707E-1 4.64594E-3 3.41675E-7 8.25890E-7
+ MIX 2 (* AXIAL BLANKET *)
+ TOTAL 1.40462E-1 2.25534E-1 3.27065E-1 3.41224E-1
+ NUSIGF 2.96101E-3 6.56171E-5 1.14630E-4 4.93483E-4
+ CHI 0.908564 0.087307 0.004129 0.0
+ SCAT 1 1 1.23805E-1
+ 2 2 2.17260E-1 1.45483E-2
+ 3 3 3.17948E-1 6.78885E-3 1.70276E-4
+ 4 4 3.31281E-1 4.38782E-3 6.04793E-6 9.37083E-7
+ MIX 3 (* AXIAL REFLECTOR *)
+ TOTAL 1.32933E-1 1.78531E-1 2.83151E-1 4.62167E-1
+ SCAT 1 1 1.22995E-1
+ 2 2 1.73095E-1 9.41231E-3
+ 3 3 2.77194E-1 5.09881E-3 1.93791E-4
+ 4 4 4.58598E-1 5.09601E-3 7.05075E-7 1.39307E-6
+ MIX 4 (* TEST ZONE *)
+ TOTAL 1.24526E-1 2.01025E-1 2.86599E-1 3.68772E-1
+ NUSIGF 1.79043E-2 1.59961E-2 2.40856E-2 7.33104E-2
+ CHI 0.908564 0.087307 0.004129 0.0
+ SCAT 1 1 1.05964E-1
+ 2 2 1.89370E-1 1.12738E-2
+ 3 3 2.70207E-1 3.64847E-3 1.46192E-4
+ 4 4 3.18960E-1 1.80479E-3 1.06888E-6 9.62178E-7
+ MIX 5 (* DRIVER WITHOUT MODERATOR *)
+ TOTAL 1.40226E-1 2.28245E-1 3.25806E-1 4.18327E-1
+ NUSIGF 1.59878E-2 1.64446E-2 2.71451E-2 8.45807E-2
+ CHI 0.908564 0.087307 0.004129 0.0
+ SCAT 1 1 1.19887E-1
+ 2 2 2.15213E-1 1.30790E-2
+ 3 3 3.06885E-1 4.00117E-3 1.59938E-4
+ 4 4 3.60906E-1 1.67341E-3 1.82716E-6 1.07166E-6
+ MIX 6 (* DRIVER WITH MODERATOR *)
+ TOTAL 1.41428E-1 2.45394E-1 3.98255E-1 4.35990E-1
+ NUSIGF 1.01663E-2 9.46359E-3 1.87325E-2 8.25335E-2
+ CHI 0.908564 0.087307 0.004129 0.0
+ SCAT 1 1 1.14337E-1
+ 2 2 2.12006E-1 2.09664E-2
+ 3 3 3.52093E-1 2.67269E-2 1.39132E-3
+ 4 4 3.70872E-1 3.29030E-2 1.08186E-3 6.10281E-5
+ MIX 7 (* REFLECTOR WITHOUT MODERATOR *)
+ TOTAL 1.59346E-1 2.16355E-1 3.48692E-1 6.24249E-1
+ SCAT 1 1 1.47969E-1
+ 2 2 2.10410E-1 1.06607E-2
+ 3 3 3.42085E-1 5.46711E-3 2.49956E-4
+ 4 4 6.19306E-1 5.36879E-3 1.00157E-6 1.82565E-6
+ MIX 8 (* REFLECTOR WITH MODERATOR *)
+ TOTAL 1.39164E-1 2.46993E-1 4.52425E-1 5.36256E-1
+ SCAT 1 1 1.05911E-1
+ 2 2 1.84820E-1 2.96485E-2
+ 3 3 3.73072E-1 5.91780E-2 3.06502E-3
+ 4 4 5.12103E-1 7.81326E-2 2.69229E-3 1.41697E-4
+ MIX 9 (* KNK1 REFLECTOR*)
+ TOTAL 1.51644E-1 1.42382E-1 1.65132E-1 8.04845E-1
+ SCAT 1 1 1.38427E-1
+ 2 2 1.37502E-1 1.23901E-2
+ 3 3 1.60722E-1 4.41927E-3 3.66930E-4
+ 4 4 7.98932E-1 3.33075E-3 1.63280E-6 1.69036E-6
+ MIX 10 (* SODIUM STEEL ZONE *)
+ TOTAL 9.65097E-2 9.87095E-2 1.34200E-1 4.12670E-1
+ SCAT 1 1 8.83550E-2
+ 2 2 9.52493E-2 7.73409E-3
+ 3 3 1.30756E-1 3.22568E-3 1.94719E-4
+ 4 4 4.09632E-1 2.90481E-3 7.98494E-7 8.89615E-7
+ MIX 11 (* SODIUM FILLED CRP *)
+ TOTAL 7.27587E-2 1.00218E-1 1.60703E-1 1.51576E-1
+ SCAT 1 1 6.63634E-2
+ 2 2 9.61236E-2 6.23393E-3
+ 3 3 1.56016E-1 4.01375E-3 7.02121E-5
+ 4 4 1.50368E-1 4.49111E-3 1.26939E-7 4.16388E-7
+ MIX 12 (* CONTROL ROD *)
+ TOTAL 1.39085E-1 2.28152E-1 3.18806E-1 6.27366E-1
+ SCAT 1 1 1.17722E-1
+ 2 2 1.94699E-1 1.26066E-2
+ 3 3 2.44352E-1 4.32219E-3 1.33314E-4
+ 4 4 3.14816E-1 3.68781E-4 1.85491E-7 1.08839E-6
+ ;
+OutMacr := Macr ;
+*
+*
+Trck := SNT: Geom ::
+ TITLE '3D TAKEDA MODEL 4 CR FULLY INSERTED'
+ EDIT 0 MAXR 20000000
+ SCHM <<sch>>
+ DIAM <<spa>>
+ SN <<ang>>
+ SCAT <<sca>>
+ NLIVO NDSA ONEG
+ MAXI 25 EPSI 1.0E-5
+ QUAD 2 (*Level-Symmetric*) ;
+OutTrck := Trck ;
+*
+*
+Syst := ASM: Macr Trck :: EDIT 0 ARM ;
+*
+*
+Flux := FLU: Macr Trck Syst ::
+ EDIT 1 TYPE K EXTE 100 1.0E-5 ;
+OutFlux := Flux ;
+*
+*
+END: ;
diff --git a/PyGan/data/rDragView_proc/README.md b/PyGan/data/rDragView_proc/README.md
new file mode 100644
index 0000000..792f99e
--- /dev/null
+++ b/PyGan/data/rDragView_proc/README.md
@@ -0,0 +1,30 @@
+# DragView
+
+`DragView` aims to help with the visualisation of DRAGON/DONJON output using the ParaView software. While it currently only works for the SN type calculation for 2D and 3D Cartesian or hexagonal geometries, it is hoped that this is flexible enough that other types of geometries and calculations can be added.
+
+### Visualisation with ParaView
+
+ParaView can be downloaded at the following [link](https://www.paraview.org/download/). Once installed, click on
+```
+File -> Open
+```
+and navigate to the folder containing the vtu files. Select the desired file and _*click on 'Apply'*_ (usually found in the left panel in the 'Properties' tab).
+
+There are countless visualisation (e.g. clip, slice, extract) and data manipulation options available natively in ParaView, and the [wiki](https://www.paraview.org/Wiki/ParaView) is a good place to start.
+
+## Contributing
+
+If you make improvements to this code or have suggestions, please do not hesitate to fork the repository or submit bug reports on github. The repository's URL is:
+```
+https://github.com/atyab191/DragView
+```
+
+## License
+
+Copyright (C) 2025 Polytechnique Montreal
+
+This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version.
+
+## Acknowledgement
+
+This work was completed during a PhD and postdoc fellowship funded by NSERC.
diff --git a/PyGan/data/rDragView_proc/get_lcm.py b/PyGan/data/rDragView_proc/get_lcm.py
new file mode 100644
index 0000000..d83512e
--- /dev/null
+++ b/PyGan/data/rDragView_proc/get_lcm.py
@@ -0,0 +1,105 @@
+#!/usr/bin/env python
+
+# Purpose:
+# Read or create lcm objects.
+#
+# Copyright:
+# Copyright (C) 2025 Polytechnique Montreal
+# This library is free software; you can redistribute it and/or modify it under the terms
+# of the GNU Lesser General Public License as published by the Free Software Foundation;
+# either version 2.1 of the License, or (at your option) any later version
+#
+# Author(s): Atyab A. Calloo
+
+from pathlib import Path
+
+import os
+import lifo
+import lcm
+import cle2000
+
+
+def read_lcm_obj(verbose):
+ """
+ Reads Geom, Macr, Trck, Flux LCM objects from current working directory.
+ Expects file name to begin with an underscore and be capitalised.
+
+ Args:
+ verbose : controls write output level
+ Returns:
+ Geom : geometry LCM object
+ Macr : macrolib LCM object
+ Trck : tracking LCM object
+ Flux : flux LCM object
+ """
+ ### FETCH VERBOSE LEVEL
+ lcm_verbose = verbose-2 if verbose > 2 else 0
+
+ ### PICK UP LCM OBJECTS
+ Geom=lcm.new(pytype='LCM_INP',name='GEOM',iact=0,impx=lcm_verbose)
+ Macr=lcm.new(pytype='LCM_INP',name='MACR',iact=0,impx=lcm_verbose)
+ Trck=lcm.new(pytype='LCM_INP',name='TRCK',iact=0,impx=lcm_verbose)
+ Flux=lcm.new(pytype='LCM_INP',name='FLUX',iact=0,impx=lcm_verbose)
+
+ if verbose > 0:
+ print('>>>>>>> SUCCESSFULLY LOADED LCM OBJECTS.')
+ print('-'*74 + '\n')
+
+ return Geom, Macr, Trck, Flux
+
+
+def run_pydragon(test_name,verbose):
+ """
+ Executes c2m procedure using DRAGON/DONJON through the cle200 class in PyGan
+ and returns Geom, Macr, Trck, Flux LCM objects.
+
+ Args:
+ test_name : file name of the c2m procedure to be computed
+ verbose : controls write output level
+ Returns:
+ Geom : geometry LCM object
+ Macr : macrolib LCM object
+ Trck : tracking LCM object
+ Flux : flux LCM object
+ Raises:
+ FileNotFound : if provided file name (test_name) does not match any of the
+ files in current working directory.
+ """
+
+ ### CHECK IF FILE EXISTS
+ extnsn = '.c2m'
+ template_file = test_name+extnsn
+ path_file = Path(template_file)
+ if not path_file.exists():
+ for entry in os.scandir('.'):
+ if entry.is_file():
+ print(entry.name)
+ raise FileNotFoundError('FILE '+template_file+' NOT FOUND IN CURRENT DIRECTORY.')
+
+ ### CONSTRUCT LIFO STACK
+ ipLifo1=lifo.new()
+ ipLifo1.pushEmpty("Geom", "LCM")
+ ipLifo1.pushEmpty("Macr", "LCM")
+ ipLifo1.pushEmpty("Trck", "LCM")
+ ipLifo1.pushEmpty("Syst", "LCM")
+ ipLifo1.pushEmpty("Flux", "LCM")
+
+ ### RUN CLE-2000 PROCEDURE
+ CLE2KFILE = cle2000.new(test_name,ipLifo1,1)
+ CLE2KFILE.exec()
+
+ ### RECOVER OUTPUT LCM OBJECTS
+ Geom = ipLifo1.node("Geom")
+ Macr = ipLifo1.node("Macr")
+ Trck = ipLifo1.node("Trck")
+ Flux = ipLifo1.node("Flux")
+
+ ### EMPTY LIFO STACK
+ while ipLifo1.getMax() > 0:
+ ipLifo1.pop()
+
+ if verbose > 0:
+ print('>>>>>>> SUCCESSFULLY RAN ',template_file,' DRAGON/DONJON CALCULATION.')
+ print('-'*74 + '\n')
+
+ return Geom, Macr, Trck, Flux \ No newline at end of file
diff --git a/PyGan/data/rDragView_proc/utils.py b/PyGan/data/rDragView_proc/utils.py
new file mode 100644
index 0000000..1b96c67
--- /dev/null
+++ b/PyGan/data/rDragView_proc/utils.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+
+# Purpose:
+# Utility functions
+#
+# Copyright:
+# Copyright (C) 2025 Polytechnique Montreal
+# This library is free software; you can redistribute it and/or modify it under the terms
+# of the GNU Lesser General Public License as published by the Free Software Foundation;
+# either version 2.1 of the License, or (at your option) any later version
+#
+# Author(s): Atyab A. Calloo
+
+import argparse
+
+
+def parserfunc():
+ """
+ Parse the input on the command line to pick up relevant parameters
+
+ Args:
+ None
+ Returns:
+ dict: of arguments
+ Raises:
+ None
+ Notes:
+ None
+ """
+ parser = argparse.ArgumentParser(prog='rDragView', description='Create vtu objects for '
+ 'visualisation of DRAGON5 flux with Paraview. Currently works for SN 2D/D in Cartesian and '
+ 'hexagonal geometries. Default run assumes that the required LCM objects are in the '
+ 'directory from where this is run. ')
+ parser.add_argument('-c', '--calc_type', required=False,
+ type=str, action="store", nargs=1,
+ default=["view"], choices={"view", "dragon_view"},
+ help='Specify type of calculation. "view": (default) use already available '
+ 'LCM objects to compute vtu files. "dragon_view": run dragon followed by '
+ 'vtu creation.')
+
+ parser.add_argument('-n', '--test_name', required=False, metavar='CALC_NAME',
+ type=str, action="store", nargs=1,
+ help='Specify a test name if desired. Used for naming vtu files. ')
+
+ parser.add_argument('-v', '--verbose', required=False, metavar='VERB',
+ type=int, action="store", nargs=1,
+ default=[1],
+ help='Verbose level, ie, level of print output.')
+
+ args = parser.parse_args()
+
+ ### ENSURE THAT TEST NAME GIVEN FOR dragon_view CALCULATION
+ if args.calc_type[0]=='dragon_view' and (args.test_name is None):
+ parser.error("dragon_view CALCULATION REQUIRES A TEST NAME.")
+
+ return args
diff --git a/PyGan/data/rDragView_proc/vtu_tools.py b/PyGan/data/rDragView_proc/vtu_tools.py
new file mode 100644
index 0000000..392c89a
--- /dev/null
+++ b/PyGan/data/rDragView_proc/vtu_tools.py
@@ -0,0 +1,463 @@
+#!/usr/bin/env python
+
+# Purpose:
+# Functions related to creating vtu files.
+#
+# Copyright:
+# Copyright (C) 2025 Polytechnique Montreal
+# This library is free software; you can redistribute it and/or modify it under the terms
+# of the GNU Lesser General Public License as published by the Free Software Foundation;
+# either version 2.1 of the License, or (at your option) any later version
+#
+# Author(s): Atyab A. Calloo
+
+import numpy as np
+import vtk
+from math import sqrt
+
+
+def create_vtu(Geom, Mcro, Trck, Flux, test_name, verbose):
+ """
+ Uses Geom, Macr, Trck, Flux LCM objects to create the material and flux vtu
+ objects that can then be viewed in Paraview.
+
+ Args:
+ Geom : geometry LCM object
+ Macr : macrolib LCM object
+ Trck : tracking LCM object
+ Flux : flux LCM object
+ test_name : file name of the c2m procedure to be computed
+ verbose : controls write output level
+ """
+
+ ### RECOVER STATEVECTORS FROM LCM OBJECTS
+ SV_geom = Geom["STATE-VECTOR"]
+ SV_mcro = Mcro["STATE-VECTOR"]
+ SV_trck = Trck["STATE-VECTOR"]
+ SV_flux = Flux["STATE-VECTOR"]
+ GV_mcro = Mcro["GROUP"]
+
+ ### RECOVER RELEVANT MATERIAL, TRACKING AND FLUX DATA
+ n_mat = SV_mcro[1]
+ n_reg = SV_trck[0]
+ n_unk = SV_trck[1]
+ n_dim = SV_trck[8]
+ spltl = SV_trck[25]
+ n_grp = SV_flux[0]
+ geom_type = SV_trck[5] # 5 - Car2D; 7 - Car3D; 8 - Hex2D; 9 - Hex3D
+ ### RECOVER OTHER DATA FROM LCM OBJECTS
+ mat_array = Trck["MATCOD"]
+
+ ### INITIALIZE FLUX AND NUSIGF ARRAYS
+ flx_array = np.zeros((n_grp,n_unk), dtype='f')
+ nsf_mat_array = np.zeros((n_grp,n_mat), dtype='f')
+ nsf_reg_array = np.zeros((n_grp,n_reg), dtype='f')
+ ### PICK UP FLUX AND NUSIGF ARRAYS
+ for i_g in range(n_grp):
+ flx_array[i_g,:] = Flux["FLUX"][i_g]
+ nsf_mat_array[i_g,:] = GV_mcro[i_g]["NUSIGF"]
+ ### CREATE NUSIGF ARRAY BY REGION
+ for i_reg in range(n_reg):
+ ind_mat = mat_array[i_reg]
+ if ind_mat == 0 :
+ nsf_reg_array[i_g,i_reg] = 0
+ elif ind_mat > 0 :
+ ind_mat=ind_mat-1
+ nsf_reg_array[i_g,i_reg] = nsf_mat_array[i_g,ind_mat]
+
+ ### REBUILD MESHES
+ lx, meshx, ly, meshy, lz, meshz = rebuild_meshes(Geom)
+
+ ### FIND HEXAGON CENTRES IF NEEDED
+ if geom_type==8 or geom_type==9:
+ side = Geom["SIDE"][0]
+ n_hex=int(n_reg/(lz*3*spltl**2))
+ n_ring = int((sqrt( float((4*n_hex-1)/3) )+1.)/2.)
+ hex_centres = compute_hex_centres(n_hex,n_ring,side)
+
+ ### CALL TO MAKE_VTU GEOM-DEPENDENT
+ nscat = SV_trck[6]
+ ielem = SV_trck[7]
+ if (geom_type == 5) or (geom_type == 7):
+ if verbose > 0:
+ print('\n>>>>>>> GENERATING VTUs FOR ',test_name,'.')
+ print('-'*74 + '\n')
+ make_vtu_cartesian_geom(lx,meshx,ly,meshy,lz,meshz,mat_array,"material",test_name )
+ for i_g in range(n_grp):
+ make_vtu_cartesian_geom(lx,meshx,ly,meshy,lz,meshz,flx_array[i_g,:],"flux",test_name ,nscat=nscat,ielem=ielem,i_g=i_g)
+ elif (geom_type == 8) or (geom_type == 9):
+ if verbose > 0:
+ print('\n>>>>>>> GENERATING VTUs FOR ',test_name,'.')
+ print('-'*74 + '\n')
+ make_vtu_hexagonal_geom(n_hex,side,hex_centres,mat_array,"material",test_name ,spltl=spltl,meshz=meshz,lz=lz)
+ for i_g in range(n_grp):
+ make_vtu_hexagonal_geom(n_hex,side,hex_centres,flx_array[i_g,:],"flux",test_name ,meshz=meshz,lz=lz,nscat=nscat,ielem=ielem,spltl=spltl,i_g=i_g)
+
+ print("\n>>>>>>> VTU CREATION COMPLETED!")
+ print('-'*74 + '\n')
+
+
+########################################################################
+def rebuild_meshes(Geom):
+ """
+ Uses Geom LCM objects to return number of mesh elements and mesh arrays along each
+ axis, taking split into account.
+
+ Args:
+ Geom : geometry LCM object
+ Returns:
+ lx : number of mesh elements along x
+ meshx : mesh array along x
+ ly : number of mesh elements along y
+ meshy : mesh array along y
+ lz : number of mesh elements along z
+ meshz : mesh array along z
+ """
+
+ SV_geom = Geom["STATE-VECTOR"]
+ geom_type = SV_geom[0] # 5 - Car2D; 7 - Car3D; 8 - Hex2D; 9 - Hex3D
+
+ mesh_dict = {}
+
+ ### FOR CARTESIAN GEOMETRIES
+ if geom_type==5 or geom_type==7:
+ lx = SV_geom[2]
+ meshx = Geom["MESHX"]
+ spltx = Geom["SPLITX"]
+ lx, meshx = get_mesh_oneaxis(lx,meshx,spltx)
+ ly = SV_geom[3]
+ meshy = Geom["MESHY"]
+ splty = Geom["SPLITY"]
+ ly, meshy = get_mesh_oneaxis(ly,meshy,splty)
+ else:
+ lx =1
+ meshx=1
+ ly =1
+ meshy=1
+
+ ### FOR 3D GEOMETRIES
+ if geom_type==7 or geom_type==9:
+ lz = SV_geom[4]
+ meshz = Geom["MESHZ"]
+ spltz = Geom["SPLITZ"]
+ lz, meshz = get_mesh_oneaxis(lz,meshz,spltz)
+ else:
+ lz =1
+ meshz=1
+ spltz=1
+
+ return lx, meshx, ly, meshy, lz, meshz
+
+
+########################################################################
+def get_mesh_oneaxis(lx,mesh,splt):
+ """
+ Compute number of mesh elements and mesh arrays taking split into account.
+
+ Args:
+ lx : number of mesh elements along axis w/o splits
+ meshx : mesh array along axis w/o splits
+ splt : split value along axis
+ Returns:
+ lx : number of mesh elements along axis w/ splits
+ meshx : mesh array along axis w/ splits
+ """
+
+ tmp_mesh=np.empty(1, dtype='f')
+ val = 0.
+ tmp_mesh[0] = val
+ for i in range(len(splt)):
+ step=(mesh[i+1]-mesh[i])/splt[i]
+ for j in range(splt[i]):
+ val += step
+ tmp_mesh = np.append(tmp_mesh,val)
+
+ lx=len(tmp_mesh)-1
+
+ return lx, tmp_mesh
+
+
+########################################################################
+def compute_hex_centres(n_hex,n_ring,side):
+ """
+ Compute centres of hexagons within hexagonal domain. Centre hexagon is zeroth
+ ring.
+ Args:
+ n_hex : number of hexagons in domain
+ n_ring : number of rings in domain; centre hexagon is zeroth ring
+ side : length of one side of hexagon
+ Returns:
+ hex_centres : 2D array with x and y coordinates of hexagon centres
+ """
+
+ hex_centres = np.zeros((2,n_hex), dtype='f')
+ ind_hex =0
+ hex_centres[0, 0] = 0.0
+ hex_centres[1, 0] = 0.0
+
+ ### LOOP OVER EACH RING OF HEXAGONS
+ for ind_ring in range(1, n_ring):
+ numhex_ring = ind_ring * 6
+ xcoord = ind_ring * 1.5 * side
+ ycoord = ind_ring * (sqrt(3) / 2) * side
+
+ hex_centres[0, ind_hex + 1] = xcoord
+ hex_centres[1, ind_hex + 1] = ycoord
+ ind_hex += 1
+
+ ind_dir = 1
+ ind_step = 0
+ numhex_perdir = ind_ring
+
+ ### LOOP OVER HEXAGONS IN CURRENT RING
+ ### For each of the six major axis directions, move in the direction
+ ### of that axis for numhex_perdir hexagons (which happens to be the
+ ### same as the ring index).
+ for i_hex in range(1, numhex_ring):
+ ind_step += 1
+ if ind_dir == 1:
+ xcoord -= 1.5 * side
+ ycoord += (sqrt(3) / 2) * side
+ elif ind_dir == 2:
+ xcoord -= 1.5 * side
+ ycoord -= (sqrt(3) / 2) * side
+ elif ind_dir == 3:
+ ycoord -= sqrt(3) * side
+ elif ind_dir == 4:
+ xcoord += 1.5 * side
+ ycoord -= (sqrt(3) / 2) * side
+ elif ind_dir == 5:
+ xcoord += 1.5 * side
+ ycoord += (sqrt(3) / 2) * side
+ elif ind_dir == 6:
+ ycoord += sqrt(3) * side
+
+ hex_centres[0, ind_hex + 1] = xcoord
+ hex_centres[1, ind_hex + 1] = ycoord
+ ind_hex += 1
+
+ if ind_step == numhex_perdir:
+ ind_dir += 1
+ ind_step = 0
+
+ return hex_centres
+
+
+########################################################################
+def writeXML(ugrid, file_name):
+ """
+ Save unstructured grid vtk object (ugrid) to file.
+ Args:
+ ugrid : vtkUnstructuredGrid object
+ file_name : file name under which to save ugrid
+ Returns:
+ None
+ """
+ writer = vtk.vtkXMLUnstructuredGridWriter()
+ writer.SetFileName(file_name)
+ writer.SetInputData(ugrid)
+ writer.Write()
+
+
+########################################################################
+def make_vtu_cartesian_geom(
+ lx,meshx,ly,meshy,lz,meshz,
+ data_np,
+ dataset_name,test_name,
+ **kwargs
+ ):
+ """
+ Build vtk unstructured grid object for 2D or 3D Cartesian mesh and populate with
+ given data, material number or flux.
+ Args:
+ lx : number of mesh elements along x
+ meshx : mesh array along x
+ ly : number of mesh elements along y
+ meshy : mesh array along y
+ lz : number of mesh elements along z
+ meshz : mesh array along z
+ data_np : given data, material number or flux
+ dataset_name : data name, usually 'material' or 'flux'
+ test_name : file name of the computed c2m procedure in case of a prior
+ dragon run or given name by user or default name
+ Returns:
+ None
+ """
+
+ ### ASSIGN OPTIONAL ARGS
+ nscat = kwargs.get("nscat", 1)
+ ielem = kwargs.get("ielem", 1)
+ i_g = kwargs.get("i_g", 0)
+ if dataset_name=="rate":
+ nsf_reg=kwargs.get("nusigf", 1)
+
+ ### 2D CONSIDERATIONS
+ if lz==1: meshz = [0,0]
+ ndim = 2 if lz==1 else 3
+
+ ### INITIALISE VTK OBJECTS
+ points = vtk.vtkPoints()
+ element = vtk.vtkVoxel()
+ cells = vtk.vtkCellArray()
+ ugrid = vtk.vtkUnstructuredGrid()
+ ### INITIALISE SCALAR CELL DATA
+ DRAGONdataset = vtk.vtkFloatArray()
+ DRAGONdataset.SetName(dataset_name)
+ ugrid.GetCellData().SetScalars(DRAGONdataset)
+
+ ### COMPUTE CONSTANTS
+ flux_offset = 1 if dataset_name=="material" else nscat*(ielem)**ndim
+
+ i_cell=-1
+ ### LOOP OVER Z-AXIS MESH
+ for kk in range(lz):
+
+ ### LOOP OVER Y-AXIS MESH
+ for jj in range(ly):
+ ycoord = [meshy[jj], meshy[jj], meshy[jj+1], meshy[jj+1]]
+
+ ### LOOP OVER X-AXIS MESH
+ for ii in range(lx):
+ xcoord = [meshx[ii], meshx[ii+1], meshx[ii+1], meshx[ii]]
+
+ ### INSERT NEW POINTS AND CELLS (DO NOT MERGE THE TWO LOOPS BELOW!!!)
+ for i_crd in range(0,4):
+ points.InsertNextPoint( xcoord[i_crd],ycoord[i_crd], meshz[kk])
+ for i_crd in range(0,4):
+ points.InsertNextPoint( xcoord[i_crd],ycoord[i_crd], meshz[kk+1])
+ for i in range(0, 8):
+ i_cell += 1
+ element.GetPointIds().SetId(i, i_cell)
+ cells.InsertNextCell(element)
+
+ ugrid.SetPoints(points)
+ ugrid.SetCells(vtk.VTK_HEXAHEDRON, cells)
+
+ ### SET SCALAR CELL DATA
+ ind1 = (kk*lx*ly) + (jj*lx) + ii
+ ind2 = (kk*lx*ly)*flux_offset + (jj*lx)*flux_offset + ii*flux_offset
+ if dataset_name=="rate":
+ DRAGONdataset.InsertTuple1(ind1, nsf_reg[ind1]*data_np[ind2])
+ else:
+ DRAGONdataset.InsertTuple1(ind1, data_np[ind2])
+
+ print(test_name, dataset_name,' PROCESSED. SAVING ...')
+ filename=test_name+'_'+dataset_name+str(i_g)+".vtu"
+ ### SAVE VTU FILE
+ writeXML(ugrid, filename)
+ return ugrid
+
+
+########################################################################
+def make_vtu_hexagonal_geom(
+ n_hex,side,hex_centres,
+ data_np,
+ dataset_name,test_name,
+ **kwargs
+ ):
+ """
+ Build vtk unstructured grid object for 2D or 3D hexagon mesh and populate with
+ given data, material number or flux.
+ Args:
+ n_hex : number of hexagons in domain
+ side : length of one side of hexagon
+ hex_centres : 2D array with x and y coordinates of hexagon centres
+ data_np : given data, material number or flux
+ dataset_name : data name, usually 'material' or 'flux'
+ test_name : file name of the computed c2m procedure in case of a prior
+ dragon run or given name by user or default name
+ Returns:
+ None
+ """
+ ### ASSIGN OPTIONAL ARGS
+ spltl = kwargs.get("spltl", 1)
+ lz = kwargs.get("lz", 1)
+ meshz = kwargs.get("meshz", 1)
+ nscat = kwargs.get("nscat", 1)
+ ielem = kwargs.get("ielem", 1)
+ i_g = kwargs.get("i_g", 0)
+
+ ### 2D CONSIDERATIONS
+ if lz==1: meshz = [0,0]
+ ndim = 2 if lz==1 else 3
+
+ ### INITIALISE VTK OBJECTS
+ points = vtk.vtkPoints()
+ lozenge = vtk.vtkVoxel()
+ cells = vtk.vtkCellArray()
+ ugrid = vtk.vtkUnstructuredGrid()
+ ### INITIALISE SCALAR CELL DATA
+ DRAGONdataset = vtk.vtkFloatArray()
+ DRAGONdataset.SetName(dataset_name)
+ ugrid.GetCellData().SetScalars(DRAGONdataset)
+
+ ### COMPUTE CONSTANTS
+ x_step=(1.0)*(side/spltl)
+ y_step=(sqrt(3)/2)*(side/spltl)
+ n_loz = 3*spltl**2
+ flux_offset = 1 if dataset_name=="material" else nscat*(ielem)**ndim
+
+ i_cell = -1
+ ### LOOP OVER Z-AXIS LAYERS
+ for i_z in range(0,lz):
+ ### LOOP OVER HEXAGONS
+ for i_hex in range(n_hex):
+ ### FETCH CENTRE OF EACH HEX
+ cx, cy = hex_centres[0, i_hex], hex_centres[1, i_hex]
+
+ ### LOOP OVER EACH HEX SECTION (EACH OF THE THREE MAIN LOZENGES)
+ for hex_section in range(3):
+ ### SET INITIAL VALUES OF XCOORD AND YCOORD
+ if hex_section == 0:
+ x_coords = [0.0, x_step*0.5, x_step, x_step*0.5]
+ y_coords = [0.0, -y_step, 0.0, y_step]
+ elif hex_section == 1:
+ x_coords = [-side, -side+x_step, -side+(x_step*1.5), -side+(x_step * 0.5)]
+ y_coords = [0.0, 0.0, y_step, y_step]
+ else:
+ x_coords = [-side/2, (-side/2)+x_step, (-side/2)+(x_step*0.5), (-side/2)-(x_step*0.5)]
+ y_coords = [-side*sqrt(3)/2, -side*sqrt(3)/2, (-side*sqrt(3)/2)+y_step, (-side*sqrt(3)/2)+y_step]
+
+ ### LOOP OVER Y-MESH INSIDE LOZENGE, INCREMENT XCOORD, YCOORD
+ for i_mshy in range(spltl):
+ if i_mshy != 0:
+ if hex_section == 0:
+ x_coords = x_coords - ((x_step/2)*(spltl-2))
+ y_coords = y_coords + y_step*(spltl)
+ elif hex_section == 1:
+ x_coords = x_coords - ((x_step)*(spltl-1.5))
+ y_coords = y_coords + y_step
+ else:
+ x_coords = x_coords - ((x_step)*(spltl-0.5))
+ y_coords = y_coords + y_step
+
+ ### LOOP OVER X-MESH INSIDE LOZENGE, INCREMENT XCOORD, YCOORD
+ for i_mshx in range(spltl):
+ if i_mshx != 0:
+ x_coords = x_coords + (x_step/2) if hex_section == 0 else x_coords + (x_step)
+ y_coords = y_coords - (y_step) if hex_section == 0 else y_coords
+
+ ### INSERT NEW POINTS AND CELLS (DO NOT MERGE THE TWO LOOPS BELOW!!!)
+ for i_crd in range(0,4):
+ points.InsertNextPoint(x_coords[i_crd]+cx, y_coords[i_crd]+cy, meshz[i_z])
+ for i_crd in range(0,4):
+ points.InsertNextPoint(x_coords[i_crd]+cx, y_coords[i_crd]+cy, meshz[i_z+1])
+ for i in range(0, 8):
+ i_cell += 1
+ lozenge.GetPointIds().SetId(i, i_cell)
+ cells.InsertNextCell(lozenge)
+
+ ugrid.SetPoints(points)
+ ugrid.SetCells(vtk.VTK_HEXAHEDRON, cells)
+
+ ### SET SCALAR CELL DATA
+ for i_loz in range(0,n_loz):
+ ind1 = (i_z*n_hex*n_loz) + (i_hex*n_loz) + i_loz
+ ind2 = (i_z*n_hex*n_loz*flux_offset) + (i_hex*n_loz*flux_offset) + (i_loz*flux_offset)
+ DRAGONdataset.InsertTuple1(ind1, data_np[ind2])
+
+ print(test_name, dataset_name,' PROCESSED. SAVING ...')
+ filename=test_name+'_'+dataset_name+str(i_g)+".vtu"
+ ### SAVE VTU FILE
+ writeXML(ugrid, filename)
+ return ugrid
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 18:08:43 +0000