diff --git a/InputFiles/c5g7_mox b/InputFiles/c5g7_mox index b8e86a4af..dbb696fc8 100644 --- a/InputFiles/c5g7_mox +++ b/InputFiles/c5g7_mox @@ -106,6 +106,7 @@ viz { centre (0.0 0.0 0.0); //width (25.0 25.0); axis z; + offset -17; res (500 500); } } @@ -119,6 +120,7 @@ nuclearData { water { temp 75675; + rgb (0 0 139); // Colour of water is dark blue moder { 1001.03 h-h2o.42; } composition { 1001.03 6.700E-002; diff --git a/NuclearData/materialMenu_mod.f90 b/NuclearData/materialMenu_mod.f90 index 3236083d2..729e209da 100644 --- a/NuclearData/materialMenu_mod.f90 +++ b/NuclearData/materialMenu_mod.f90 @@ -1,13 +1,14 @@ !! -!! Material Menu is a module (Singleton) that contains global definitions of diffrent materials +!! Material Menu is a module (Singleton) that contains global definitions of different materials !! !! It exists to make it easier for all databases to refer to the same materials by the same -!! name and index. This is necessary to avoid confusion resulting from diffrent materials with the -!! same name or index in diffrent databases. +!! name and index. This is necessary to avoid confusion resulting from different materials with the +!! same name or index in different databases. !! !! Public Members: !! materialDefs -> array of material definitions of type materialItem !! nameMap -> Map that maps material name to matIdx +!! colourMap -> Map that maps matIdx to 24bit colour (to use for visualisation) !! !! Interface: !! init -> Load material definitions from a dictionary @@ -21,8 +22,10 @@ module materialMenu_mod use numPrecision - use universalVariables, only : NOT_FOUND, VOID_MAT, OUTSIDE_MAT + use universalVariables, only : NOT_FOUND, VOID_MAT, OUTSIDE_MAT, UNDEF_MAT use genericProcedures, only : fatalError, charToInt, numToChar + use colours_func, only : rgb24bit + use intMap_class, only : intMap use charMap_class, only : charMap use dictionary_class, only : dictionary @@ -33,9 +36,9 @@ module materialMenu_mod !! Information about a single nuclide !! !! Based somewhat on MCNP conventions. - !! Atomic and Mass number identify cleary a nuclide species + !! Atomic and Mass number identify clearly a nuclide species !! Evaluation number T allows to refer to multiple states/evaluations of the same nuclide species - !! E.G. at a Diffrent temperature as in MCNP Library. + !! E.G. at a Different temperature as in MCNP Library. !! !! Public members: !! Z -> Atomic number @@ -83,6 +86,7 @@ module materialMenu_mod !! 5010.03 2.0E-005; !! } !! xsFile /home/uberMoffTarkin/XS/mat1.xs; + !! #rgb (255 0 0); # // RGB colour to be used in visualisation !! } !! !! NOTE: the moder dictionary is optional, necessary only if S(a,b) thermal scattering @@ -102,9 +106,16 @@ module materialMenu_mod procedure :: display => display_materialItem end type materialItem -!! MODULE COMPONENTS + !! Parameters + integer(shortInt), parameter :: COL_OUTSIDE = int(z'ffffff', shortInt) + integer(shortInt), parameter :: COL_VOID = int(z'000000', shortInt) + integer(shortInt), parameter :: COL_UNDEF = int(z'00ff00', shortInt) + + + !! MODULE COMPONENTS type(materialItem),dimension(:),allocatable,target,public :: materialDefs - type(charMap),target,public :: nameMap + type(charMap), target, public :: nameMap + type(intMap), public :: colourMap public :: init public :: kill @@ -131,7 +142,7 @@ subroutine init(dict) integer(shortInt) :: i character(nameLen) :: temp - ! Clean whatever may be alrady present + ! Clean whatever may be already present call kill() ! Load all material names @@ -142,17 +153,20 @@ subroutine init(dict) ! Load definitions do i=1,size(matNames) - call materialDefs(i) % init(matNames(i), dict % getDictPtr(matNames(i))) - materialDefs(i) % matIdx = i + call materialDefs(i) % init(matNames(i), i, dict % getDictPtr(matNames(i))) call nameMap % add(matNames(i), i) end do - ! Add special Material keywords to thedictionary + ! Add special Material keywords to the dictionary temp = 'void' call nameMap % add(temp, VOID_MAT) temp = 'outside' call nameMap % add(temp, OUTSIDE_MAT) + !! Load colours for the special materials + call colourMap % add(VOID_MAT, COL_VOID) + call colourMap % add(OUTSIDE_MAT, COL_OUTSIDE) + call colourMap % add(UNDEF_MAT, COL_UNDEF) end subroutine init @@ -204,7 +218,7 @@ end subroutine display !! Result: !! nameLen long character with material name !! - !! Erorrs: + !! Error: !! If idx is -ve or larger then number of defined materials !! Empty string '' is returned as its name !! @@ -250,25 +264,30 @@ end function matIdx !! !! Args: !! name [in] -> character with material name + !! idx [in] -> material index !! dict [in] -> dictionary with material definition !! !! Errors: !! FatalError if dictionary does not contain valid material definition. !! - subroutine init_materialItem(self, name, dict) - class(materialItem), intent(inout) :: self - character(nameLen),intent(in) :: name - class(dictionary), intent(in) :: dict - character(nameLen),dimension(:),allocatable :: keys, moderKeys - integer(shortInt) :: i - class(dictionary),pointer :: compDict, moderDict - logical(defBool) :: hasSab + subroutine init_materialItem(self, name, idx, dict) + class(materialItem), intent(inout) :: self + character(nameLen), intent(in) :: name + integer(shortInt), intent(in) :: idx + class(dictionary), intent(in) :: dict + character(nameLen), dimension(:), allocatable :: keys, moderKeys + integer(shortInt), dimension(:), allocatable :: temp + integer(shortInt) :: i + class(dictionary),pointer :: compDict, moderDict + logical(defBool) :: hasSab + character(100), parameter :: Here = 'init_materialItem (materialMenu_mod.f90)' ! Return to initial state call self % kill() ! Load easy components c self % name = name + self % matIdx = idx call dict % get(self % T,'temp') ! Get composition dictionary and load composition @@ -300,6 +319,17 @@ subroutine init_materialItem(self, name, dict) call self % nuclides(i) % init(keys(i)) end do + ! Add colour info if present + if(dict % isPresent('rgb')) then + call dict % get(temp, 'rgb') + + if (size(temp) /= 3) then + call fatalError(Here, "'rgb' keyword must have 3 values") + end if + + call colourMap % add(idx, rgb24bit(temp(1), temp(2), temp(3))) + end if + ! Save dictionary self % extraInfo = dict @@ -392,7 +422,7 @@ function isNucDefinition(key) result(isIt) za = verify(key(1:L),'0123456789') tt = verify(key(1:L),'0123456789', back = .true.) - ! Verify that the location of the dot is consistant + ! Verify that the location of the dot is consistent isIt = dot == za .and. dot == tt end function isNucDefinition @@ -437,7 +467,7 @@ end subroutine init_nuclideInfo !! None !! !! Result: - !! Character in format ZZAAA.TT that dscribes nuclide definition + !! Character in format ZZAAA.TT that describes nuclide definition !! !! Errors: !! None @@ -461,17 +491,17 @@ end function toChar_nuclideInfo !! Get pointer to a material definition under matIdx !! !! Args: - !! matIdx [in] -> Index of the material + !! idx [in] -> Index of the material !! !! Result: !! Pointer to a materialItem with the definition !! !! Errors: - !! FatalError if matIdx does not correspond to any defined material + !! FatalError if idx does not correspond to any defined material !! FatalError if material definitions were not loaded !! - function getMatPtr(matIdx) result(ptr) - integer(shortInt), intent(in) :: matIdx + function getMatPtr(idx) result(ptr) + integer(shortInt), intent(in) :: idx type(materialItem), pointer :: ptr character(100), parameter :: Here = 'getMatPtr (materialMenu_mod.f90)' @@ -481,13 +511,13 @@ function getMatPtr(matIdx) result(ptr) end if ! Verify matIdx - if( matIdx <= 0 .or. matIdx > nMat()) then - call fatalError(Here,"matIdx: "//numToChar(matIdx)// & + if( idx <= 0 .or. idx > nMat()) then + call fatalError(Here,"matIdx: "//numToChar(idx)// & " does not correspond to any defined material") end if ! Attach pointer - ptr => materialDefs(matIdx) + ptr => materialDefs(idx) end function getMatPtr diff --git a/SharedModules/CMakeLists.txt b/SharedModules/CMakeLists.txt index dcd1633ec..2a8baa4ce 100644 --- a/SharedModules/CMakeLists.txt +++ b/SharedModules/CMakeLists.txt @@ -12,7 +12,8 @@ add_sources( ./genericProcedures.f90 ./timer_mod.f90 ./charLib_func.f90 ./openmp_func.f90 - ./errors_mod.f90) + ./errors_mod.f90 + ./colours_func.f90) add_unit_tests( ./Tests/grid_test.f90 ./Tests/energyGrid_test.f90 @@ -21,4 +22,5 @@ add_unit_tests( ./Tests/grid_test.f90 ./Tests/hashFunctions_test.f90 ./Tests/timer_test.f90 ./Tests/conversions_test.f90 - ./Tests/charLib_test.f90) + ./Tests/charLib_test.f90 + ./Tests/colours_test.f90) diff --git a/SharedModules/Tests/colours_test.f90 b/SharedModules/Tests/colours_test.f90 new file mode 100644 index 000000000..3aba28b02 --- /dev/null +++ b/SharedModules/Tests/colours_test.f90 @@ -0,0 +1,26 @@ +module colours_test + use numPrecision + use colours_func, only : rgb24bit + use pFUnit_mod + + implicit none + +contains + + +@Test + subroutine testColourConversions() + + ! Test by comparison with some hex values + @assertEqual(int(z"123456"), rgb24bit(18, 52, 86)) + + @assertEqual(int(z"000000"), rgb24bit(0, 0, 0)) + @assertEqual(int(z"ffffff"), rgb24bit(255, 255, 255)) + + @assertEqual(int(z"ff0000"), rgb24bit(255, 0, 0)) + @assertEqual(int(z"00ff00"), rgb24bit(0, 255, 0)) + @assertEqual(int(z"0000ff"), rgb24bit(0, 0, 255)) + + end subroutine testColourConversions + +end module colours_test diff --git a/SharedModules/colours_func.f90 b/SharedModules/colours_func.f90 new file mode 100644 index 000000000..5e7fec61a --- /dev/null +++ b/SharedModules/colours_func.f90 @@ -0,0 +1,48 @@ +!! +!! Contains function to deal with colours +!! +module colours_func + use numPrecision + use genericProcedures, only: numToChar + use errors_mod, only: fatalError + implicit none + + public :: rgb24bit + +contains + + + !! + !! Return integer with 24bit colour value from r,g,b components + !! + !! Args: + !! r [in] -> red component in [0-255] range + !! g [in] -> green component in [0-255] range + !! b [in] -> blue component in [0-255] range + !! + function rgb24bit(r, g, b) result(col) + integer(shortInt), intent(in) :: r + integer(shortInt), intent(in) :: g + integer(shortInt), intent(in) :: b + integer(shortInt) :: col + character(100), parameter :: Here = "rgb24bit (colours_func.f90)" + + ! Check that the values are in the correct range + if (r < 0 .or. r > 255) then + call fatalError(Here, "Red value is out of [0-255] range: " // numToChar(r)) + end if + if (g < 0 .or. g > 255) then + call fatalError(Here, "Green value is out of [0-255] range: " // numToChar(g)) + end if + if (b < 0 .or. b > 255) then + call fatalError(Here, "Blue value is out of [0-255] range: " // numToChar(b)) + end if + + ! Compute the 24 bit colour + ! Note inversion, blue is the least significant bits + col = b + 256 * g + 256 * 256 * r + + end function rgb24bit + + +end module colours_func diff --git a/Visualisation/visualiser_class.f90 b/Visualisation/visualiser_class.f90 index 9641f7758..a0cdcd768 100644 --- a/Visualisation/visualiser_class.f90 +++ b/Visualisation/visualiser_class.f90 @@ -1,394 +1,398 @@ -module visualiser_class - - use numPrecision - use universalVariables - use genericProcedures, only : fatalError, numToChar - use hashFunctions_func, only : knuthHash - use imgBmp_func, only : imgBmp_toFile - use commandLineUI, only : getInputFile - use dictionary_class, only : dictionary - use geometry_inter, only : geometry - use outputVTK_class - - implicit none - private - - !! - !! Object responsible for controlling visualisation - !! - !! Object that creates images relating to SCONE geometries - !! Should be extensible for adding different visualisation methods - !! Recieves and generates data for visualisation - !! Requires a dictionary input which specifies the procedures to call - !! Presently supports: VTK voxel mesh creation - !! - !! Private members: - !! name -> name to be used for generating output files (corresponding to input) - !! geom -> pointer to geometry - !! vizDict -> dictionary containing visualisations to be generated - !! - !! Interface: - !! init -> initialises visualiser - !! makeViz -> constructs requested visualisations - !! kill -> cleans up visualiser - !! - !! Sample dictionary input: - !! viz{ - !! vizDict1{ } - !! #vizDict2{ }# - !! } - !! - !! NOTE: For details regarding contents of the vizDict dictionaries see the documentation - !! of 'makeVTK' and 'makeBmpImg' functions - !! - type, public :: visualiser - character(nameLen), private :: name - class(geometry), pointer, private :: geom => null() - type(dictionary), private :: vizDict - contains - procedure :: init - procedure :: makeViz - procedure :: kill - procedure, private :: makeVTK - procedure, private :: makeBmpImg - end type - -contains - - !! - !! Initialises visualiser - !! - !! Provides visualiser with filename for output, - !! geometry information, and the dictionary decribing - !! what is to be plotted - !! - !! Args: - !! geom [inout] -> pointer to the geometry - !! vizDict[in] -> dictionary containing what is to be visualised - !! - !! Result: - !! Initialised visualiser - !! - subroutine init(self, geom, vizDict) - class(visualiser), intent(inout) :: self - class(geometry), pointer, intent(inout) :: geom - class(dictionary), intent(in) :: vizDict - character(:), allocatable :: string - - ! Obtain file name - call getInputFile(string) - self % name = string - - ! Point to geometry - self % geom => geom - - ! Store visualisation dictionary - self % vizDict = vizDict - - end subroutine init - - !! - !! Generate all visualisations specified by vizDict - !! - !! Proceed through all dictionaries contained within vizDict - !! and perform all corresponding visualisations - !! - !! Result: - !! Visualisation outputs corresponding to dictionary contents - !! - !! Errors: - !! Returns an error if an unrecognised visualisation is requested - !! - subroutine makeViz(self) - class(visualiser), intent(inout) :: self - class(dictionary), pointer :: tempDict - character(nameLen),dimension(:), allocatable :: keysArr - integer(shortInt) :: i - character(nameLen) :: type - character(nameLen) :: here ='makeViz (visualiser_class.f90)' - - ! Loop through each sub-dictionary and generate visualisation - ! (if the visualisation method is available) - call self % vizDict % keys(keysArr,'dict') - - do i=1,size(keysArr) - tempDict => self % vizDict % getDictPtr(keysArr(i)) - call tempDict % get(type,'type') - select case(type) - case('vtk') - call self % makeVTK(tempDict) - - case('bmp') - call self % makeBmpImg(tempDict) - - case default - call fatalError(here, 'Unrecognised visualisation - presently only accept vtk') - - end select - - end do - - end subroutine makeViz - - !! - !! Generate a VTK output - !! - !! Creates the VTK file corresponding to the contents of dict - !! - !! Args: - !! dict [in] -> dictionary containing description of VTK file to be made - !! - !! Sample input dictionary: - !! VTK { - !! type vtk; - !! corner (-1.0 -1.0 -1.0); // lower corner of the plot volume - !! width (2.0 2.0 2.0); // width in each direction - !! vox (300 300 300); // Resolution in each direction - !! #what uniqueId;# // Plot target. 'material' or 'uniqueId'. Default: 'material' - !! } - !! - !! TODO: VTK output is placed in a input filename appended by '.vtk' extension. - !! This prevents multiple VTK visualistions (due to overriding). Might also become - !! weird for input files with extension e.g. 'input.dat'. - !! DEMAND USER TO GIVE OUTPUT NAME - !! - subroutine makeVTK(self, dict) - class(visualiser), intent(inout) :: self - class(dictionary), intent(in) :: dict - type(outputVTK) :: vtk - integer(shortInt), dimension(:,:,:), allocatable:: voxelMat - real(defReal), dimension(:), allocatable :: corner ! corner of the mesh - real(defReal), dimension(:), allocatable :: center ! center of the mesh - real(defReal), dimension(:), allocatable :: width ! corner of the mesh - integer(shortInt), dimension(:), allocatable :: nVox ! number of mesh voxels - character(nameLen) :: what - character(nameLen) :: here ='makeVTK (visualiser_class.f90)' - - call vtk % init(dict) - - ! Identify whether plotting 'material' or 'cellID' - call dict % getOrDefault(what, 'what', 'material') - - ! Obtain geometry data - call dict % get(corner, 'corner') - call dict % get(width, 'width') - center = corner + width/TWO - call dict % get(nVox, 'vox') - - if (size(corner) /= 3) then - call fatalError(here,'Voxel plot requires corner to have 3 values') - endif - if (size(width) /= 3) then - call fatalError(here,'Voxel plot requires width to have 3 values') - endif - if (size(nVox) /= 3) then - call fatalError(here,'Voxel plot requires vox to have 3 values') - endif - allocate(voxelMat(nVox(1), nVox(2), nVox(3))) - - ! Have geometry obtain data - call self % geom % voxelPlot(voxelMat, center, what, width) - - ! In principle, can add multiple data sets to VTK - not done here yet - ! VTK data set will use 'what' variable as a name - call vtk % addData(voxelMat, what) - call vtk % output(self % name) - call vtk % kill() - - end subroutine makeVTK - - !! - !! Generate a BMP slice image of the geometry - !! - !! Args: - !! dict [in] -> Dictionary with settings - !! - !! Sample dictionary input: - !! bmp_img { - !! type bmp; - !! #what uniqueID;# // Target of the plot. 'uniqueId' or 'material'. Default: 'material' - !! output img; // Name of output file without extension - !! centre (0.0 0.0 0.0); // Coordinates of the centre of the plot - !! axis x; // Must be 'x', 'y' or 'z' - !! res (300 300); // Resolution of the image - !! #width (1.0 2.0);# // Width of the plot from the centre - !! } - !! - !! NOTE: If 'width' is not given, the plot will extend to the bounds of the geometry. - !! This may result in the provided centre beeing moved to the center of the geoemtry in the - !! plot plane. However, the position on the plot axis will be unchanged. - !! - subroutine makeBmpImg(self, dict) - class(visualiser), intent(inout) :: self - class(dictionary), intent(in) :: dict - real(defReal), dimension(3) :: centre - real(defReal), dimension(2) :: width - character(1) :: dir - character(nameLen) :: tempChar - logical(defBool) :: useWidth - character(nameLen) :: what, outputFile - real(defReal), dimension(:), allocatable :: temp - integer(shortInt), dimension(:), allocatable :: tempInt - integer(shortInt), dimension(:,:), allocatable :: img - character(100), parameter :: Here = 'makeBmpImg (visualiser_class.f90)' - - ! Get plot parameters - - ! Identify whether plotting 'material' or 'cellID' - call dict % getOrDefault(what, 'what', 'material') - - ! Get name of the output file - call dict % get(outputFile, 'output') - outputFile = trim(outputFile) // '.bmp' - - ! Central point - call dict % get(temp, 'centre') - - if (size(temp) /= 3) then - call fatalError(Here, "'center' must have size 3. Has: "//numToChar(size(temp))) - end if - - centre = temp - - ! Axis - call dict % get(tempChar, 'axis') - - if (len_trim(tempChar) /= 1) then - call fatalError(Here, "'axis' must be x,y or z. Not: "//tempChar) - end if - - dir = tempChar(1:1) - - ! Resolution - call dict % get(tempInt, 'res') - - if (size(tempInt) /= 2) then - call fatalError(Here, "'res' must have size 2. Has: "//numToChar(size(tempInt))) - else if (any(tempInt <= 0)) then - call fatalError(Here, "Resolution must be +ve. There is 0 or -ve entry!") - end if - - allocate(img(tempInt(1), tempInt(2))) - - ! Optional width - useWidth = dict % isPresent('width') - if (useWidth) then - call dict % get(temp, 'width') - - ! Check for errors - if (size(temp) /= 2) then - call fatalError(Here, "'width' must have size 2. Has: "//numToChar((size(temp)))) - else if (any(temp <= ZERO)) then - call fatalError(Here, "'width' must be +ve. It isn't.") - end if - - width = temp - - end if - - ! Get plot - if (useWidth) then - call self % geom % slicePlot(img, centre, dir, what, width) - else - call self % geom % slicePlot(img, centre, dir, what) - end if - - ! Translate to an image - select case (what) - case ('material') - img = materialColor(img) - - case ('uniqueID') - img = uniqueIDColor(img) - - case default - call fatalError(Here, "Invalid request for plot target. Must be 'material' or 'uniqueID'& - & is: "//what) - end select - - ! Print image - call imgBmp_toFile(img, outputFile) - - end subroutine makeBmpImg - - !! - !! Terminates visualiser - !! - !! Cleans up remnants of visualiser once it is no longer needed - !! - !! Result: - !! An empty visualiser object - !! - subroutine kill(self) - class(visualiser), intent(inout) :: self - - self % name ='' - self % geom => null() - call self % vizDict % kill() - - end subroutine kill - - - !! - !! Convert matIdx to a 24bit color - !! - !! Special materials are associeted with special colors: - !! OUTSIDE_MAT -> white (#ffffff) - !! VOID_MAT -> black (#000000) - !! UNDEF_MAT -> green (#00ff00) - !! - !! Args: - !! matIdx [in] -> Value of the material index - !! - !! Result: - !! A 24-bit color specifing the material - !! - elemental function materialColor(matIdx) result(color) - integer(shortInt), intent(in) :: matIdx - integer(shortInt) :: color - integer(shortInt), parameter :: COL_OUTSIDE = int(z'ffffff', shortInt) - integer(shortInt), parameter :: COL_VOID = int(z'000000', shortInt) - integer(shortInt), parameter :: COL_UNDEF = int(z'00ff00', shortInt) - - select case (matIdx) - case (OUTSIDE_MAT) - color = COL_OUTSIDE - - case (VOID_MAT) - color = COL_VOID - - case (UNDEF_MAT) - color = COL_UNDEF - - case default - color = knuthHash(matIdx, 24) - - end select - - end function materialColor - - !! - !! Convert uniqueID to 24bit color - !! - !! An elemental wrapper over Knuth Hash - !! - !! We use a hash function to scatter colors accross all available. - !! Knuth multiplicative hash is very good at scattering integer - !! sequences e.g. {1, 2, 3...}. Thus, it is ideal for a colormap. - !! - !! Args: - !! uniqueID [in] -> Value of the uniqueID - !! - !! Result: - !! A 24-bit color specifing the uniqueID - !! - elemental function uniqueIDColor(uniqueID) result(color) - integer(shortInt), intent(in) :: uniqueID - integer(shortInt) :: color - - color = knuthHash(uniqueID, 24) - - end function uniqueIDColor - - -end module visualiser_class +module visualiser_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, numToChar + use hashFunctions_func, only : knuthHash, FNV_1 + use imgBmp_func, only : imgBmp_toFile + use commandLineUI, only : getInputFile + use dictionary_class, only : dictionary + use geometry_inter, only : geometry + use materialMenu_mod, only : mm_colourMap => colourMap + use outputVTK_class + + implicit none + private + + !! + !! Object responsible for controlling visualisation + !! + !! Object that creates images relating to SCONE geometries + !! Should be extensible for adding different visualisation methods + !! Receives and generates data for visualisation + !! Requires a dictionary input which specifies the procedures to call + !! Presently supports: VTK voxel mesh creation + !! + !! Private members: + !! name -> name to be used for generating output files (corresponding to input) + !! geom -> pointer to geometry + !! vizDict -> dictionary containing visualisations to be generated + !! + !! Interface: + !! init -> initialises visualiser + !! makeViz -> constructs requested visualisations + !! kill -> cleans up visualiser + !! + !! Sample dictionary input: + !! viz{ + !! vizDict1{ } + !! #vizDict2{ }# + !! } + !! + !! NOTE: For details regarding contents of the vizDict dictionaries see the documentation + !! of 'makeVTK' and 'makeBmpImg' functions + !! + type, public :: visualiser + character(nameLen), private :: name + class(geometry), pointer, private :: geom => null() + type(dictionary), private :: vizDict + contains + procedure :: init + procedure :: makeViz + procedure :: kill + procedure, private :: makeVTK + procedure, private :: makeBmpImg + end type + +contains + + !! + !! Initialises visualiser + !! + !! Provides visualiser with filename for output, + !! geometry information, and the dictionary describing + !! what is to be plotted + !! + !! Args: + !! geom [inout] -> pointer to the geometry + !! vizDict[in] -> dictionary containing what is to be visualised + !! + !! Result: + !! Initialised visualiser + !! + subroutine init(self, geom, vizDict) + class(visualiser), intent(inout) :: self + class(geometry), pointer, intent(inout) :: geom + class(dictionary), intent(in) :: vizDict + character(:), allocatable :: string + + ! Obtain file name + call getInputFile(string) + self % name = string + + ! Point to geometry + self % geom => geom + + ! Store visualisation dictionary + self % vizDict = vizDict + + end subroutine init + + !! + !! Generate all visualisations specified by vizDict + !! + !! Proceed through all dictionaries contained within vizDict + !! and perform all corresponding visualisations + !! + !! Result: + !! Visualisation outputs corresponding to dictionary contents + !! + !! Errors: + !! Returns an error if an unrecognised visualisation is requested + !! + subroutine makeViz(self) + class(visualiser), intent(inout) :: self + class(dictionary), pointer :: tempDict + character(nameLen),dimension(:), allocatable :: keysArr + integer(shortInt) :: i + character(nameLen) :: type + character(nameLen) :: here ='makeViz (visualiser_class.f90)' + + ! Loop through each sub-dictionary and generate visualisation + ! (if the visualisation method is available) + call self % vizDict % keys(keysArr,'dict') + + do i=1,size(keysArr) + tempDict => self % vizDict % getDictPtr(keysArr(i)) + call tempDict % get(type,'type') + select case(type) + case('vtk') + call self % makeVTK(tempDict) + + case('bmp') + call self % makeBmpImg(tempDict) + + case default + call fatalError(here, 'Unrecognised visualisation - presently only accept vtk') + + end select + + end do + + end subroutine makeViz + + !! + !! Generate a VTK output + !! + !! Creates the VTK file corresponding to the contents of dict + !! + !! Args: + !! dict [in] -> dictionary containing description of VTK file to be made + !! + !! Sample input dictionary: + !! VTK { + !! type vtk; + !! corner (-1.0 -1.0 -1.0); // lower corner of the plot volume + !! width (2.0 2.0 2.0); // width in each direction + !! vox (300 300 300); // Resolution in each direction + !! #what uniqueId;# // Plot target. 'material' or 'uniqueId'. Default: 'material' + !! } + !! + !! TODO: VTK output is placed in a input filename appended by '.vtk' extension. + !! This prevents multiple VTK visualisations (due to overriding). Might also become + !! weird for input files with extension e.g. 'input.dat'. + !! DEMAND USER TO GIVE OUTPUT NAME + !! + subroutine makeVTK(self, dict) + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + type(outputVTK) :: vtk + integer(shortInt), dimension(:,:,:), allocatable:: voxelMat + real(defReal), dimension(:), allocatable :: corner ! corner of the mesh + real(defReal), dimension(:), allocatable :: center ! center of the mesh + real(defReal), dimension(:), allocatable :: width ! corner of the mesh + integer(shortInt), dimension(:), allocatable :: nVox ! number of mesh voxels + character(nameLen) :: what + character(nameLen) :: here ='makeVTK (visualiser_class.f90)' + + call vtk % init(dict) + + ! Identify whether plotting 'material' or 'cellID' + call dict % getOrDefault(what, 'what', 'material') + + ! Obtain geometry data + call dict % get(corner, 'corner') + call dict % get(width, 'width') + center = corner + width/TWO + call dict % get(nVox, 'vox') + + if (size(corner) /= 3) then + call fatalError(here,'Voxel plot requires corner to have 3 values') + endif + if (size(width) /= 3) then + call fatalError(here,'Voxel plot requires width to have 3 values') + endif + if (size(nVox) /= 3) then + call fatalError(here,'Voxel plot requires vox to have 3 values') + endif + allocate(voxelMat(nVox(1), nVox(2), nVox(3))) + + ! Have geometry obtain data + call self % geom % voxelPlot(voxelMat, center, what, width) + + ! In principle, can add multiple data sets to VTK - not done here yet + ! VTK data set will use 'what' variable as a name + call vtk % addData(voxelMat, what) + call vtk % output(self % name) + call vtk % kill() + + end subroutine makeVTK + + !! + !! Generate a BMP slice image of the geometry + !! + !! Args: + !! dict [in] -> Dictionary with settings + !! + !! Sample dictionary input: + !! bmp_img { + !! type bmp; + !! #what uniqueID;# // Target of the plot. 'uniqueId' or 'material'. Default: 'material' + !! output img; // Name of output file without extension + !! centre (0.0 0.0 0.0); // Coordinates of the centre of the plot + !! axis x; // Must be 'x', 'y' or 'z' + !! res (300 300); // Resolution of the image + !! #offset 978; # // Parameter to 'randomize' the colour map + !! #width (1.0 2.0);# // Width of the plot from the centre + !! } + !! + !! NOTE: If 'width' is not given, the plot will extend to the bounds of the geometry. + !! This may result in the provided centre being moved to the center of the geometry in the + !! plot plane. However, the position on the plot axis will be unchanged. + !! + subroutine makeBmpImg(self, dict) + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal), dimension(3) :: centre + real(defReal), dimension(2) :: width + character(1) :: dir + character(nameLen) :: tempChar + logical(defBool) :: useWidth + character(nameLen) :: what, outputFile + real(defReal), dimension(:), allocatable :: temp + integer(shortInt), dimension(:), allocatable :: tempInt + integer(shortInt), dimension(:,:), allocatable :: img + integer(shortInt) :: offset + character(10) :: time + character(8) :: date + character(100), parameter :: Here = 'makeBmpImg (visualiser_class.f90)' + + ! Get plot parameters + + ! Identify whether plotting 'material' or 'cellID' + call dict % getOrDefault(what, 'what', 'material') + + ! Get name of the output file + call dict % get(outputFile, 'output') + outputFile = trim(outputFile) // '.bmp' + + ! Central point + call dict % get(temp, 'centre') + + if (size(temp) /= 3) then + call fatalError(Here, "'centre' must have size 3. Has: "//numToChar(size(temp))) + end if + + centre = temp + + ! Axis + call dict % get(tempChar, 'axis') + + if (len_trim(tempChar) /= 1) then + call fatalError(Here, "'axis' must be x,y or z. Not: "//tempChar) + end if + + dir = tempChar(1:1) + + ! Resolution + call dict % get(tempInt, 'res') + + if (size(tempInt) /= 2) then + call fatalError(Here, "'res' must have size 2. Has: "//numToChar(size(tempInt))) + else if (any(tempInt <= 0)) then + call fatalError(Here, "Resolution must be +ve. There is 0 or -ve entry!") + end if + + allocate(img(tempInt(1), tempInt(2))) + + ! Optional width + useWidth = dict % isPresent('width') + if (useWidth) then + call dict % get(temp, 'width') + + ! Check for errors + if (size(temp) /= 2) then + call fatalError(Here, "'width' must have size 2. Has: "//numToChar((size(temp)))) + else if (any(temp <= ZERO)) then + call fatalError(Here, "'width' must be +ve. It isn't.") + end if + + width = temp + + end if + + ! Colourmap offset + ! If not given select randomly + if (dict % isPresent('offset')) then + call dict % get(offset, 'offset') + + else + call date_and_time(date, time) + call FNV_1(date // time, offset) + + end if + + ! Get plot + if (useWidth) then + call self % geom % slicePlot(img, centre, dir, what, width) + else + call self % geom % slicePlot(img, centre, dir, what) + end if + + ! Translate to an image + select case (what) + case ('material') + img = materialColour(img, offset) + + case ('uniqueID') + img = uniqueIDColour(img) + + case default + call fatalError(Here, "Invalid request for plot target. Must be 'material' or 'uniqueID'& + & is: "//what) + end select + + ! Print image + call imgBmp_toFile(img, outputFile) + + end subroutine makeBmpImg + + !! + !! Terminates visualiser + !! + !! Cleans up remnants of visualiser once it is no longer needed + !! + !! Result: + !! An empty visualiser object + !! + subroutine kill(self) + class(visualiser), intent(inout) :: self + + self % name ='' + self % geom => null() + call self % vizDict % kill() + + end subroutine kill + + + !! + !! Convert matIdx to a 24bit colour + !! + !! Special materials are associated with special colours: + !! OUTSIDE_MAT -> white (#ffffff) + !! VOID_MAT -> black (#000000) + !! UNDEF_MAT -> green (#00ff00) + !! + !! Args: + !! matIdx [in] -> Value of the material index + !! offset [in] -> Offset to be used in the hash function + !! + !! Result: + !! A 24-bit colour specifying the material + !! + elemental function materialColour(matIdx, offset) result(colour) + integer(shortInt), intent(in) :: matIdx + integer(shortInt) :: colour + integer(shortInt), intent(in) :: offset + + ! Since Knuth hash is cheap we can compute it anyway even if colour from + ! the map will end up being used + colour = mm_colourMap % getOrDefault(matIdx, knuthHash(matIdx + offset, 24)) + + end function materialColour + + !! + !! Convert uniqueID to 24bit colour + !! + !! An elemental wrapper over Knuth Hash + !! + !! We use a hash function to scatter colours across all available. + !! Knuth multiplicative hash is very good at scattering integer + !! sequences e.g. {1, 2, 3...}. Thus, it is ideal for a colourmap. + !! + !! Args: + !! uniqueID [in] -> Value of the uniqueID + !! + !! Result: + !! A 24-bit colour specifying the uniqueID + !! + elemental function uniqueIDColour(uniqueID) result(colour) + integer(shortInt), intent(in) :: uniqueID + integer(shortInt) :: colour + + colour = knuthHash(uniqueID, 24) + + end function uniqueIDColour + + +end module visualiser_class diff --git a/docs/User Manual.rst b/docs/User Manual.rst index 1f941ba41..8b6e72922 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -275,9 +275,9 @@ neutronCEstd, to perform analog collision processing target nuclide movement. Target movement is sampled if target mass A < massThreshold. [Mn] * DBRCeMin (*optional*, default = 1.0e-08): minimum DBRC energy. [MeV] * DBRCeMax (*optional*, default = 2.0e-04): maximum DBRC energy. [MeV] - + Example: :: - + collisionOperator { neutronCE { type neutronCEstd; minEnergy 1.0e-12; maxEnergy 30.0; energyThreshold 200; massThreshold 2; DBRCeMin 1.0e-06; DBRCeMax 0.001; } } @@ -693,6 +693,9 @@ bmp * what (*optional*, default = material): defines what is highlighted in the plot; options are ``material`` and ``uniqueID``, where ``uniqueID`` highlights unique cell IDs +* offset (*optional*, default = random) An integer (positive or negative) that + shifts the sequence of colours assigned to materials. Allows to change colours + from the default sequence in a parametric way. Example: :: @@ -739,7 +742,7 @@ from ACE files. resonance probability tables treatment * DBRC (*optional*, default = no DBRC): list of ZAIDs of nuclides for which DBRC has to be applied. - + Example: :: ceData { type aceNuclearDatabase; aceLibrary ./myFolder/ACElib/JEF311.aceXS; @@ -748,7 +751,7 @@ Example: :: .. note:: If DBRC is applied, the 0K cross section ace files of the relevant nuclides must be included in the aceLibrary file. - + baseMgNeutronDatabase ##################### @@ -804,6 +807,9 @@ Other options are: * xsFile: needed for multi-group simulations. Must contain the path to the file where the multi-group cross sections are stored. +* rgb (*optional*): An array of three integers specifying the RGB colour e.g. ``(255 0 0)``. The + colour defined in this way will be used for visualisation of the material in the geometry plots. + Example 1: :: materials { @@ -814,6 +820,7 @@ Example 1: :: 8016.03 0.018535464; } } water { temp 273; + rgb (0 0 200); composition { 1001.03 0.0222222; 8016.03 0.00535; } @@ -1087,7 +1094,7 @@ Example: :: } } -.. note:: +.. note:: To calculate the average weight, one should divide weight moment 1 (weight1) by weight moment 0 (weight0). To calculate the variance of the weights, the tally results have to be post-processed as: var = weight2/weight0 - (weight1/weight0)^2