diff --git a/Apps/scone.f90 b/Apps/scone.f90 index 6e45f8ec0..8502b226b 100644 --- a/Apps/scone.f90 +++ b/Apps/scone.f90 @@ -10,6 +10,7 @@ program scone use physicsPackage_inter, only : physicsPackage use physicsPackageFactory_func, only : new_physicsPackage use vizPhysicsPackage_class, only : vizPhysicsPackage + use renderPhysicsPackage_class, only : renderPhysicsPackage use timer_mod , only : registerTimer, timerStart, timerStop, timerTime, secToChar implicit none @@ -27,6 +28,8 @@ program scone call addClOption('--omp', 1, ['int'], & 'Number of OpenMP threads in a parallel calculation') #endif + call addClOption('--render', 0, ['int'],& + 'Executes command-line ray plotting specified by a viz dict in the input file') ! Get path to input file call getInputFile(inputPath) @@ -54,6 +57,9 @@ program scone if (clOptionIsPresent('--plot')) then allocate(vizPhysicsPackage :: core) call core % init(input) + else if (clOptionIsPresent('--render')) then + allocate(renderPhysicsPackage :: core) + call core % init(input) else allocate( core, source = new_physicsPackage(input)) endif diff --git a/Geometry/geometry_inter.f90 b/Geometry/geometry_inter.f90 index c473f5194..d14c0d8d2 100644 --- a/Geometry/geometry_inter.f90 +++ b/Geometry/geometry_inter.f90 @@ -1,7 +1,7 @@ module geometry_inter use numPrecision - use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INFINITY, COLL_EV + use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INF, COLL_EV use genericProcedures, only : fatalError use dictionary_class, only : dictionary use charMap_class, only : charMap @@ -476,8 +476,9 @@ end subroutine voxelPlot !! mats [in] -> List of materials which are transparent to rays !! fov [in] -> Field-of-view in the horizontal axis in radians !! ambient [in] -> Fraction of illumination which is ambient, rather than from the light + !! bounds [in] -> [-x, -y, -z, +x, +y, +z] bounds outside of which a ray hit can't occur !! - subroutine rayPlot(self, brightness, matIDs, camera, light, M, mats, fov, ambient) + subroutine rayPlot(self, brightness, matIDs, camera, light, M, mats, fov, ambient, bounds) class(geometry), intent(in) :: self real(defReal), dimension(:,:), intent(inout) :: brightness integer(shortInt), dimension(:,:), intent(inout) :: matIDs @@ -487,6 +488,7 @@ subroutine rayPlot(self, brightness, matIDs, camera, light, M, mats, fov, ambien integer(shortInt), dimension(:), intent(in) :: mats real(defReal), intent(in) :: fov real(defReal), intent(in) :: ambient + real(defReal), dimension(6), intent(in) :: bounds integer(shortInt) :: iVert, nHoriz, nVert integer(shortInt), save :: iHoriz, matIdx real(defReal), save :: bright @@ -516,7 +518,7 @@ subroutine rayPlot(self, brightness, matIDs, camera, light, M, mats, fov, ambien call ray % init(camera, matmul(M,vec)) ! Get local material and its illumination - call self % phongTrace(ray, matIdx, bright, mats, ambient, light) + call self % phongTrace(ray, matIdx, bright, mats, ambient, light, bounds) matIDs(iHoriz, iVert) = matIdx brightness(iHoriz, iVert) = bright @@ -530,7 +532,7 @@ end subroutine rayPlot !! Procedure for tracing from a camera until an opaque object is hit. !! Then calculates the luminous contribution from a light source. !! - subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) + subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light, bounds) class(geometry), intent(in) :: self type(coordlist), intent(inout) :: ray integer(shortInt), intent(out) :: matIdx @@ -538,32 +540,53 @@ subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) integer(shortInt), dimension(:), intent(in) :: mats real(defReal), intent(in) :: ambient real(defReal), dimension(3), intent(in) :: light - logical(defBool) :: lightTrace - real(defReal) :: dist, dotP, dNudge + real(defReal), dimension(6), intent(in) :: bounds + logical(defBool) :: lightTrace, intersects + real(defReal) :: dist, dotP, dNudge, & + dMin, dMax integer(shortInt) :: event, matIdx0 - real(defReal), dimension(3) :: dir, normal0, normal + real(defReal), dimension(3) :: dir, normal0, normal, & + nMin, nMax lightTrace = .true. - call self % placeCoord(ray) + ! Compute allowable distances between which the ray is inside the box + call boxIntersection(ray % lvl(1) % r, ray % lvl(1) % dir, bounds, & + dMin, dMax, nMin, nMax, intersects) + ! Move ray straight to box intersection point + ray % lvl(1) % r = ray % lvl(1) % r + dMin * ray % lvl(1) % dir + call self % placeCoord(ray) matIdx = ray % matIdx - ! Trace until an opaque material is hit + normal0 = nMin + + ! Early escape if the ray never intersects the opaque box + if (.not. intersects) then + bright = ZERO + return + end if + + ! Trace until an opaque material or the tracking boundary is hit + dMax = dMax - dMin do while (any(mats == matIdx)) - dist = INFINITY + dist = dMax call self % moveNoBC(ray, dist, event, normal0) matIdx = ray % matIdx - ! If distance is infinite or particle collided, not pointing towards anything opaque. + ! If particle collided outside or distance is infinite, + ! it is not pointing towards anything opaque. ! Hence, short circuit the trace - if ((dist == INFINITY) .or. (event == COLL_EV)) then + if (event == COLL_EV .or. dist == INF) then lightTrace = .false. exit end if + ! Decrement the distance to the boundary + dMax = dMax - dist + end do bright = ZERO @@ -588,25 +611,37 @@ subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) ! Compute product betwen normal and light direction dotP = max(dot_product(normal0, dir), ZERO) + + ! Compute the short-circuit distance until ray escapes the box + call boxIntersection(ray % lvl(1) % r, ray % lvl(1) % dir, bounds, & + dMin, dMax, nMin, nMax, intersects) ! Does the ray fly directly to the light? - matIdx0 = ray % matIdx - do while (any(mats == matIdx0)) + ! Goes through the box towards light + if (intersects) then + matIdx0 = ray % matIdx + do while (any(mats == matIdx0)) - ! Find maximum flight distance - dist = norm2(ray % lvl(1) % r - light) - call self % moveNoBC(ray, dist, event, normal) + ! Find maximum flight distance + dist = min(norm2(ray % lvl(1) % r - light), dMax) + call self % moveNoBC(ray, dist, event, normal) - matIdx0 = ray % matIdx + matIdx0 = ray % matIdx - ! The ray flew straight to the light - if (event == COLL_EV) then - bright = (ONE - ambient) * dotP + ! The ray flew straight to the light + if (event == COLL_EV) then + bright = (ONE - ambient) * dotP + exit + end if - exit - end if + ! Decrement the distance until the box + dMax = dMax - dist - end do + end do + ! Needn't cross the box - goes straight to the light + else + bright = (ONE - ambient) * dotP + end if ! Add ambient contribution bright = bright + ambient @@ -615,4 +650,82 @@ subroutine phongTrace(self, ray, matIdx, bright, mats, ambient, light) end subroutine phongTrace + !! + !! Returns the distance and normals for the two intersections + !! of a ray on a box, as well as whether it intersects at all. + !! Helper function for phongTrace. + !! + subroutine boxIntersection(r, omega, bounds, dMin, dMax, nMin, nMax, intersects) + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: omega + real(defReal), dimension(6), intent(in) :: bounds + real(defReal), intent(out) :: dMin + real(defReal), intent(out) :: dMax + real(defReal), dimension(3), intent(out) :: nMin + real(defReal), dimension(3), intent(out) :: nMax + logical(defBool), intent(out) :: intersects + integer(shortInt) :: i + real(defReal) :: d1, d2, invDir + real(defReal), dimension(3) :: n1, n2 + + dMin = -INF + dMax = INF + nMin = ZERO + nMax = ZERO + intersects = .true. + + do i = 1, 3 + + if (abs(omega(i)) < floatTol) then + ! Ray is parallel to slab + if (r(i) < bounds(i) .or. r(i) > bounds(i + 3)) then + dMax = -INF + exit + end if + else + invDir = ONE / omega(i) + + d1 = (bounds(i) - r(i)) * invDir + d2 = (bounds(i + 3) - r(i)) * invDir + + n1 = ZERO + n2 = ZERO + + if (omega(i) > ZERO) then + n1(i) = -ONE + n2(i) = ONE + else + n1(i) = ONE + n2(i) = -ONE + invDir = d1 + d1 = d2 + d2 = invDir + end if + + if (d1 > dMin) then + dMin = d1 + nMin = n1 + end if + + if (d2 < dMax) then + dMax = d2 + nMax = n2 + end if + + ! No intersection - escape + if (dMin > dMax) then + dMin = 0 + dMax = INF + intersects = .false. + return + end if + end if + + end do + + if (dMin < ZERO) dMin = ZERO + if (dMax < dMin) intersects = .false. + + end subroutine boxIntersection + end module geometry_inter diff --git a/PhysicsPackages/CMakeLists.txt b/PhysicsPackages/CMakeLists.txt index 24723bd94..5f18802ac 100644 --- a/PhysicsPackages/CMakeLists.txt +++ b/PhysicsPackages/CMakeLists.txt @@ -8,5 +8,6 @@ add_sources( ./physicsPackage_inter.f90 ./kineticPhysicsPackage_class.f90 ./alphaPhysicsPackage_class.f90 ./rayVolPhysicsPackage_class.f90 + ./renderPhysicsPackage_class.f90 ) #./dynamPhysicsPackage_class.f90) diff --git a/PhysicsPackages/physicsPackageFactory_func.f90 b/PhysicsPackages/physicsPackageFactory_func.f90 index cb0040205..25ba26c8c 100644 --- a/PhysicsPackages/physicsPackageFactory_func.f90 +++ b/PhysicsPackages/physicsPackageFactory_func.f90 @@ -17,6 +17,7 @@ module physicsPackageFactory_func use rayVolPhysicsPackage_class, only : rayVolPhysicsPackage use kineticPhysicsPackage_class, only : kineticPhysicsPackage use alphaPhysicsPackage_class, only : alphaPhysicsPackage + use renderPhysicsPackage_class, only : renderPhysicsPackage ! use dynamPhysicsPackage_class, only : dynamPhysicsPackage implicit none @@ -31,6 +32,7 @@ module physicsPackageFactory_func 'vizPhysicsPackage ',& 'kineticPhysicsPackage ',& 'alphaPhysicsPackage ',& + 'renderPhysicsPackage ',& 'rayVolPhysicsPackage '] !! @@ -73,6 +75,9 @@ function new_physicsPackage(dict) result(new) case('rayVolPhysicsPackage') allocate( rayVolPhysicsPackage :: new) + case('renderPhysicsPackage') + allocate( renderPhysicsPackage :: new) + case default print *, AVAILABLE_physicsPackages call fatalError(Here, 'Unrecognised type of Physics Package : ' // trim(type)) diff --git a/PhysicsPackages/renderPhysicsPackage_class.f90 b/PhysicsPackages/renderPhysicsPackage_class.f90 new file mode 100644 index 000000000..a08b1ab28 --- /dev/null +++ b/PhysicsPackages/renderPhysicsPackage_class.f90 @@ -0,0 +1,377 @@ +module renderPhysicsPackage_class + + use numPrecision + use universalVariables + use genericProcedures, only : fatalError, rotateAroundAxis, crossProduct, numToChar + use dictionary_class, only : dictionary + + ! Physics package interface + use physicsPackage_inter, only : physicsPackage + + ! Geometry + use geometry_inter, only : geometry + use geometryReg_mod, only : gr_geomPtr => geomPtr, gr_geomIdx => geomIdx + use geometryFactory_func, only : new_geometry + + ! Nuclear Data + use materialMenu_mod, only : mm_nMat => nMat, mm_matIdx => matIdx, mm_matName => matName + use nuclearDataReg_mod, only : ndReg_init => init ,& + ndReg_getMatNames => getMatNames + use nuclearDatabase_inter, only : nuclearDatabase + + ! Visualisation + use imgBmp_func, only : imgBmp_toFile + use visualiser_class, only : visualiser, materialColour, brightnessScale, & + getRayPlotInfo, getRayPlotTransformation + + implicit none + private + + !! + !! Physics Package for live rendering of ray plots + !! As well as a geometry, requires only a visualiser dictionary + !! containing a ray plot. + !! + !! The rendering allows modifying the ray plot live. Requires the presence + !! of an app for viewing/updating images (eog works well). + !! + type, public,extends(physicsPackage) :: renderPhysicsPackage + private + ! Building blocks + class(geometry), pointer :: geom => null() + integer(shortInt) :: geomIdx = 0 + type(visualiser) :: viz + type(dictionary) :: dict + character(nameLen) :: displayApp = 'eog' + + contains + procedure :: init + procedure :: run + procedure :: kill + + end type renderPhysicsPackage + +contains + + !! + !! Calls visualiser to generate visualisation + !! + subroutine run(self) + class(renderPhysicsPackage), intent(inout) :: self + real(defReal) :: fov, ambient + integer(shortInt), dimension(2) :: res + integer(shortInt), dimension(:), allocatable :: mats, temp + integer(shortInt), dimension(:,:), allocatable :: img, matIDs + real(defReal), dimension(:,:), allocatable :: lum + real(defReal), dimension(3,3) :: M + integer(shortInt) :: offset, ios, matIdx, pos, n, i, j + character(nameLen) :: outputFile, outputLive + real(defReal) :: a, b, c, d, e, f + real(defReal), dimension(3) :: dir, dirV, dirH, centre, camera,& + light, up, right + real(defReal), dimension(6) :: bounds + character(100) :: line + character(nameLen) :: cmd, matName, argStr + character(100), parameter :: Here = 'run (renderPhysicsPackage_class.f90)' + + ! Extract ray info from the viz dictionary + call getRayPlotInfo(self % dict, outputFile, centre, camera, light, up, M, fov, ambient, & + res, mats, offset, bounds) + + allocate(matIDs(res(1), res(2))) + allocate(lum(res(1), res(2))) + + dirV = M(:,3) + dirH = M(:,2) + + ! Overwrite output + outputFile = 'output_temp.bmp' + outputLive = 'output.bmp' + + print *, 'COMMANDS:' + print *, '1D Pan camera: w/a/s/d ' + print *, '2D Pan camera: pan ' + print *, 'Rotate camera: rot ' + print *, 'Zoom camera: +/- ' + print *, 'Set "up": up <3D direction>' + print *, 'Set light: light <3D position>' + print *, 'Set camera: camera <3D position>' + print *, 'Set centre: centre <3D position>' + print *, 'Set ambient light: amb ' + print *, 'Make a material opaque: opaq ' + print *, 'Make a material transparent: transp ' + print *, 'Set rendered volume: bounds <-x -y -z +x +y +z>' + print *, 'Set resolution in x and y directions: res ' + print *, 'Provide list of material names: mats' + print *, 'Provide the positions of the camera, centre, and light: pos' + print *, 'Quit: q' + + ! lum contains luminosity values, matIDs identifies which materials were hit + call self % geom % rayPlot(lum, matIDs, camera, light, M, mats, fov, ambient, bounds) + + ! Translate to an image. + ! Obtain material colours and scale by luminosity + matIDs = materialColour(matIDs, offset) + img = matIDs + + ! Scale image by brightness + img = brightnessScale(img, lum) + + ! Print image + call imgBmp_toFile(img, outputFile) + call execute_command_line("mv "//outputFile//" "//outputLive) + call execute_command_line(self % displayApp//" "//outputLive//" &") + + print *,'Input command:' + do + + a = ZERO + b = ZERO + c = ZERO + cmd = '' + read(*,'(A)') line + read(line, *) cmd + + pos = index(line, ' ') + if (pos > 0) then + argStr = adjustl(line(pos+1:)) + else + argStr = '' + end if + + select case (trim(cmd)) + case ("+") + read(argStr, *, iostat=ios) a + dir = centre - camera + dir = dir / norm2(dir) + camera = camera + a * dir + case ("-") + read(argStr, *, iostat=ios) a + dir = centre - camera + dir = dir / norm2(dir) + camera = camera - a * dir + case ("w") + read(argStr, *, iostat=ios) a + centre = centre + a * dirV + camera = camera + a * dirV + case ("s") + read(argStr, *, iostat=ios) a + centre = centre - a * dirV + camera = camera - a * dirV + case ("a") + read(argStr, *, iostat=ios) a + centre = centre - a * dirH + camera = camera - a * dirH + case ("d") + read(argStr, *, iostat=ios) a + centre = centre + a * dirH + camera = camera + a * dirH + case ("pan") + read(argStr, *, iostat=ios) a, b + centre = centre + a * dirH + b * dirV + camera = camera + a * dirH + b * dirV + case ("rot") + read(argStr, *, iostat=ios) a, b + + dir = camera - centre + + ! Yaw + dir = rotateAroundAxis(dir, up, a * PI / 180) + + ! Pitch + right = crossProduct(dir, up) + right = right / norm2(right) + dir = rotateAroundAxis(dir, right, b * PI / 180) + + camera = centre + dir + + M = getRayPlotTransformation(camera, centre, up) + dirV = M(:,3) + dirH = M(:,2) + + case ("up") + read(argStr, *, iostat=ios) a, b, c + up = [a, b, c] + if (all(up == [ZERO, ZERO, ZERO])) up = [ZERO, ZERO, ONE] + up = up /norm2(up) + M = getRayPlotTransformation(camera, centre, up) + dirV = M(:,3) + dirH = M(:,2) + case ("camera") + read(argStr, *, iostat=ios) a, b, c + camera = [a, b, c] + M = getRayPlotTransformation(camera, centre, up) + dirV = M(:,3) + dirH = M(:,2) + case ("centre") + read(argStr, *, iostat=ios) a, b, c + centre = [a, b, c] + M = getRayPlotTransformation(camera, centre, up) + dirV = M(:,3) + dirH = M(:,2) + case ("light") + read(argStr, *, iostat=ios) a, b, c + light = [a, b, c] + case ("bounds") + read(argStr, *, iostat=ios) a, b, c, d, e, f + bounds = [a, b, c, d, e, f] + case ("amb") + read(argStr, *, iostat=ios) a + a = max(a, ZERO) + a = min(a, ONE) + ambient = a + + case ("transp") + matName = trim(argStr) + matIdx = mm_matIdx(matName) + + if (matIdx == NOT_FOUND) then + print *,'Material '//trim(matName)//' is not present in the simulation.' + cycle + end if + + ! Check mats to see if matIdx is present + if (any(mats == matIdx)) then + print *,'Material '//trim(matName)//' is already transparent.' + cycle + end if + + n = size(mats) + + allocate(temp(n+1)) + temp(1:n) = mats + temp(n+1) = matIdx + + call move_alloc(temp, mats) + + case ("opaq") + matName = trim(argStr) + matIdx = mm_matIdx(matName) + + if (matIdx == NOT_FOUND) then + print *,'Material '//trim(matName)//' is not present in the simulation.' + cycle + end if + + ! Check mats to see if matIdx is present + if (all(mats /= matIdx)) then + print *,'Material '//trim(matName)//' is already opaque.' + cycle + end if + + mats = pack(mats, mats /= matIdx) + + case("mats") + print *,'MATERIAL NAMES:' + do i = 1, mm_nMat() + print *,trim(mm_matName(i)) + end do + cycle + + case("res") + read(argStr, *, iostat=ios) i, j + + if (i < 1 .or. j < 1) then + print *,'Invalid pixels values provided' + cycle + end if + + res(1) = i + res(2) = j + deallocate(matIDs, lum) + allocate(matIDs(res(1), res(2))) + allocate(lum(res(1), res(2))) + + case ("pos") + print *,'Camera position: '//numToChar(camera) + print *,'Centre position: '//numToChar(centre) + print *,'Light position: '//numToChar(light) + cycle + case ("q") + exit + case default + print *,'Unrecognised command' + cycle + end select + + ! lum contains luminosity values, matIDs identifies which materials were hit + call self % geom % rayPlot(lum, matIDs, camera, light, M, mats, fov, ambient, bounds) + + ! Translate to an image. + ! Obtain material colours and scale by luminosity + matIDs = materialColour(matIDs, offset) + img = matIDs + + ! Scale image by brightness + img = brightnessScale(img, lum) + + ! Print image + call imgBmp_toFile(img, outputFile) + + ! Display image + call execute_command_line("mv "//outputFile//" "//outputLive) + + end do + + end subroutine run + + !! + !! Initialise from individual components and dictionaries + !! + subroutine init(self, dict) + class(renderPhysicsPackage), intent(inout) :: self + class(dictionary), intent(inout) :: dict + class(dictionary),pointer :: tempDict, vizDict + character(nameLen) :: geomName + logical(defBool) :: found + integer(shortInt) :: i + character(nameLen), dimension(:), allocatable :: keys + character(nameLen) :: vizType + character(100), parameter :: Here ='init (renderPhysicsPackage_class.f90)' + + ! Read the app used to display the live-updated image + call dict % getOrDefault(self % displayApp, "app", "eog") + + ! Build Nuclear Data + call ndReg_init(dict % getDictPtr("nuclearData")) + + ! Build geometry + tempDict => dict % getDictPtr('geometry') + geomName = 'visualGeom' + call new_geometry(tempDict, geomName) + self % geomIdx = gr_geomIdx(geomName) + self % geom => gr_geomPtr(self % geomIdx) + + ! Find ray dictionary + if (.not. dict % isPresent('viz')) call fatalError(here,'Must provide viz dict for plotting') + + ! Go through dictionary until a subdictionary has a rayPlot + vizDict => dict % getDictPtr('viz') + call vizDict % keys(keys) + + found = .false. + do i = 1, size(keys) + tempDict => vizDict % getDictPtr(keys(i)) + + call tempDict % get(vizType, 'type') + if (vizType == 'ray') then + found = .true. + self % dict = tempDict + end if + end do + + if (.not. found) call fatalError(Here,'Must provide a rayPlot dictionary') + + end subroutine init + + !! + !! Deallocate memory + !! + subroutine kill(self) + class(renderPhysicsPackage), intent(inout) :: self + + ! TODO: This subroutine + + end subroutine kill + +end module renderPhysicsPackage_class diff --git a/SharedModules/genericProcedures.f90 b/SharedModules/genericProcedures.f90 index e256253c4..556659d6c 100644 --- a/SharedModules/genericProcedures.f90 +++ b/SharedModules/genericProcedures.f90 @@ -1139,6 +1139,25 @@ subroutine rotationMatrix(matrix, phi, theta, psi) end subroutine rotationMatrix + !! + !! Perform a rotation of vector dir about a given axis over angle (in RADIANS) + !! to produce newDir. Useful for camera rotations + !! + function rotateAroundAxis(dir, axis, angle) result(newDir) + real(defReal), dimension(3), intent(in) :: dir + real(defReal), dimension(3), intent(in) :: axis + real(defReal), intent(in) :: angle + real(defReal), dimension(3) :: newDir + real(defReal) :: c, s + real(defReal), dimension(3) :: k + + k = axis / norm2(axis) + c = cos(angle) + s = sin(angle) + + newDir = dir * c + crossProduct(k, dir) * s + k * dot_product(k, dir)*(ONE - c) + + end function rotateAroundAxis !! !! Dot product for 3D vector !! diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index 69d92f068..74fb59995 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -18,9 +18,7 @@ module universalVariables integer(shortInt), parameter, public :: MAX_COL = 70 ! Maximum number of columns in console display ! Define variables which are important for tracking neutrons in the geometry - real(defReal), parameter, public :: INFINITY = 2.0_defReal**63, & - surface_tol = 1.0e-12_defReal, & ! Tol. on closeness to surface - SURF_TOL = 1.0E-12_defReal, & + real(defReal), parameter, public :: SURF_TOL = 1.0E-12_defReal, & INF = 2.0_defReal**63, & NUDGE = 1.0e-8_defReal ! Distance to poke neutrons across boundaries for surface tracking diff --git a/TransportOperator/transportOperatorHT_class.f90 b/TransportOperator/transportOperatorHT_class.f90 index c1461b4a0..1bc9c0015 100644 --- a/TransportOperator/transportOperatorHT_class.f90 +++ b/TransportOperator/transportOperatorHT_class.f90 @@ -195,8 +195,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else diff --git a/TransportOperator/transportOperatorST_class.f90 b/TransportOperator/transportOperatorST_class.f90 index f73d27563..80a35b7ef 100644 --- a/TransportOperator/transportOperatorST_class.f90 +++ b/TransportOperator/transportOperatorST_class.f90 @@ -70,8 +70,8 @@ subroutine surfaceTracking(self, p, tally, thisCycle, nextCycle) ! This branch is called in the case of voids with no imposed XS if (sigmaTrack < tol) then - dist = INFINITY - invSigmaTrack = INFINITY + dist = INF + invSigmaTrack = INF sigmaT = ZERO else diff --git a/Visualisation/visualiser_class.f90 b/Visualisation/visualiser_class.f90 index b7aaae9fa..70b1691ea 100644 --- a/Visualisation/visualiser_class.f90 +++ b/Visualisation/visualiser_class.f90 @@ -14,6 +14,11 @@ module visualiser_class implicit none private + public :: materialColour + public :: brightnessScale + public :: getRayPlotInfo + public :: getRayPlotTransformation + !! !! Object responsible for controlling visualisation !! @@ -361,147 +366,31 @@ end subroutine makeBmpImg !! #fov 70;# // Field-of-view in the horizontal axis in degrees !! #offset 978; # // Parameter to 'randomize' the colour map !! #transparent (mat1, mat2, ... );# // Names of transparent materials + !! #bounds (-2 -5 10 4 -3 12);# // Bounds outside of which everything is transparent !! } !! !! subroutine makeRayPlot(self, dict) - class(visualiser), intent(inout) :: self - class(dictionary), intent(in) :: dict - real(defReal) :: fov, ambient - real(defReal), dimension(3) :: centre, up, camera, light, cv, ch, d - real(defReal), dimension(:), allocatable :: temp - integer(shortInt), dimension(:), allocatable :: tempInt - character(nameLen), dimension(:), allocatable :: matNames + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal) :: fov, ambient + real(defReal), dimension(3) :: centre, camera, light, up + integer(shortInt), dimension(2) :: res integer(shortInt), dimension(:), allocatable :: mats integer(shortInt), dimension(:,:), allocatable :: img, matIDs real(defReal), dimension(:,:), allocatable :: lum real(defReal), dimension(3,3) :: M - integer(shortInt) :: i, offset - character(10) :: time - character(8) :: date + integer(shortInt) :: offset + real(defReal), dimension(6) :: bounds character(nameLen) :: outputFile - character(100), parameter :: Here = 'makeRayPlot (visualiser_class.f90)' - ! Get name of the output file - call dict % get(outputFile, 'output') - outputFile = trim(outputFile) // '.bmp' - - ! Central/target 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 - deallocate(temp) - - ! Camera origin - call dict % get(temp, 'camera') - - if (size(temp) /= 3) then - call fatalError(Here, "'camera' must have size 3. Has: "//numToChar(size(temp))) - end if - camera = temp - deallocate(temp) + call getRayPlotInfo(dict, outputFile, centre, camera, light, up, M, fov, ambient, & + res, mats, offset, bounds) - ! Get light location or default to camera location - if (dict % isPresent('light')) then - call dict % get(temp, 'light') - - if (size(temp) /= 3) then - call fatalError(Here, "'light' must have size 3. Has: "//numToChar(size(temp))) - end if - light = temp - deallocate(temp) + allocate(matIDs(res(1), res(2))) + allocate(lum(res(1), res(2))) - else - light = camera - end if - - ! The up direction, which sets the camera rotation - if (dict % isPresent('up')) then - call dict % get(temp, 'up') - - if (size(temp) /= 3) then - call fatalError(Here, "'up' must have size 3. Has: "//numToChar(size(temp))) - end if - up = temp - up = up / norm2(up) - deallocate(temp) - - else - up = [0, 0, 1] - end if - - ! Optional field of view in degrees - ! Convert to radians - call dict % getOrDefault(fov, 'fov', 70.0_defReal) - fov = fov * PI / 180 - - ! Optional fraction of diffuse light - call dict % getOrDefault(ambient, 'ambient', 0.3_defReal) - if ((ambient > 1) .or. (ambient < 0)) call fatalError(Here,'The fraction of diffuse light must be between 0 and 1') - - ! Create governing vectors for the plot - d = (centre - camera) - d = d/norm2(d) - - ! Ensure that up is not colinear with the view direction - if (all(abs(crossProduct(up, d)) < 1E-6)) call fatalError(Here,"View direction is co-linear with 'up'.") - - cv = crossProduct(d, up) - cv = cv / norm2(cv) - - ch = crossProduct(cv, d) - ch = ch / norm2(ch) - - ! Create coordinate matrix - M(:,1) = d - M(:,2) = cv - M(:,3) = ch - - ! 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(matIDs(tempInt(1), tempInt(2))) - allocate(lum(tempInt(1), tempInt(2))) - - ! Transparent materials - ! Defaults to an array containing only VOID and OUTSIDE - if (dict % isPresent('transparent')) then - call dict % get(matNames, 'transparent') - allocate(mats(size(matNames) + 2)) - do i = 1, size(matNames) - mats(i) = mm_nameMap % get(matNames(i)) - end do - mats(i) = OUTSIDE_MAT - mats(i+1) = VOID_MAT - else - allocate(mats(2)) - mats(1) = OUTSIDE_MAT - mats(2) = VOID_MAT - 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 - - ! Img contains luminosity values, matIDs identifies which materials were hit - call self % geom % rayPlot(lum, matIDs, camera, light, M, mats, fov, ambient) + call self % geom % rayPlot(lum, matIDs, camera, light, M, mats, fov, ambient, bounds) ! Translate to an image. ! Obtain material colours and scale by luminosity @@ -608,4 +497,176 @@ elemental function brightnessScale(img0, brightness) result(img) end function brightnessScale + !! + !! Helper function to extract input for the ray plot for use by + !! the live renderer + !! + subroutine getRayPlotInfo(dict, outputFile, centre, camera, light, up, M, fov, ambient, & + res, mats, offset, bounds) + class(dictionary), intent(in) :: dict + character(nameLen), intent(out) :: outputFile + real(defReal), dimension(3), intent(out) :: centre + real(defReal), dimension(3), intent(out) :: camera + real(defReal), dimension(3), intent(out) :: light + real(defReal), dimension(3) :: up + real(defReal), dimension(3, 3), intent(out) :: M + real(defReal), intent(out) :: fov + real(defReal), intent(out) :: ambient + integer(shortInt), dimension(2), intent(out) :: res + integer(shortInt), dimension(:), allocatable, intent(out) :: mats + integer(shortInt), intent(out) :: offset + real(defReal), dimension(6), intent(out) :: bounds + real(defReal), dimension(:), allocatable :: temp + integer(shortInt), dimension(:), allocatable :: tempInt + character(nameLen), dimension(:), allocatable :: matNames + character(10) :: time + character(8) :: date + integer(shortInt) :: i + character(100), parameter :: Here = 'getRayPlotInfo (visualiser_class.f90)' + + ! Get name of the output file + call dict % get(outputFile, 'output') + outputFile = trim(outputFile) // '.bmp' + + ! Central/target 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 + deallocate(temp) + + ! Camera origin + call dict % get(temp, 'camera') + + if (size(temp) /= 3) then + call fatalError(Here, "'camera' must have size 3. Has: "//numToChar(size(temp))) + end if + camera = temp + deallocate(temp) + + ! Get light location or default to camera location + if (dict % isPresent('light')) then + call dict % get(temp, 'light') + + if (size(temp) /= 3) then + call fatalError(Here, "'light' must have size 3. Has: "//numToChar(size(temp))) + end if + light = temp + deallocate(temp) + + else + light = camera + end if + + ! The up direction, which sets the camera rotation + if (dict % isPresent('up')) then + call dict % get(temp, 'up') + + if (size(temp) /= 3) then + call fatalError(Here, "'up' must have size 3. Has: "//numToChar(size(temp))) + end if + up = temp + up = up / norm2(up) + deallocate(temp) + + else + up = [0, 0, 1] + end if + + ! Optional field of view in degrees + ! Convert to radians + call dict % getOrDefault(fov, 'fov', 70.0_defReal) + fov = fov * PI / 180 + + ! Optional fraction of diffuse light + call dict % getOrDefault(ambient, 'ambient', 0.3_defReal) + if ((ambient > 1) .or. (ambient < 0)) call fatalError(Here,'The fraction of diffuse light must be between 0 and 1') + + M = getRayPlotTransformation(camera, centre, up) + + ! 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 + + res = tempInt + + ! Transparent materials + ! Defaults to an array containing only VOID and OUTSIDE + if (dict % isPresent('transparent')) then + call dict % get(matNames, 'transparent') + allocate(mats(size(matNames) + 2)) + do i = 1, size(matNames) + mats(i) = mm_nameMap % get(matNames(i)) + end do + mats(i) = OUTSIDE_MAT + mats(i+1) = VOID_MAT + else + allocate(mats(2)) + mats(1) = OUTSIDE_MAT + mats(2) = VOID_MAT + 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 + + ! Bounds of the transparent region + if (dict % isPresent('bounds')) then + call dict % get(temp, 'bounds') + if (size(temp) /= 6) call fatalError(Here,'Size of "bounds" must be 6 (-x -y -z +x +y +z)') + bounds = temp + else + bounds = [-INF, -INF, -INF, INF, INF, INF] + end if + + end subroutine getRayPlotInfo + + !! + !! Returns the transformation matrix for a ray plotted view given + !! the three defining vectors of camera position, view target/centre, and + !! the upwards direction of the view. + !! + function getRayPlotTransformation(camera, centre, up) result(M) + real(defReal), dimension(3), intent(in) :: camera + real(defReal), dimension(3), intent(in) :: centre + real(defReal), dimension(3), intent(in) :: up + real(defReal), dimension(3,3) :: M + real(defReal), dimension(3) :: d, cv, ch + character(100), parameter :: Here = 'getRayPlotTransformation (visualiser_class.f90)' + + ! Create governing vectors for the plot + d = (centre - camera) + d = d/norm2(d) + + ! Ensure that up is not colinear with the view direction + if (all(abs(crossProduct(up, d)) < 1E-6)) call fatalError(Here,"View direction is co-linear with 'up'.") + + cv = crossProduct(d, up) + cv = cv / norm2(cv) + + ch = crossProduct(cv, d) + ch = ch / norm2(ch) + + ! Create coordinate matrix + M(:,1) = d + M(:,2) = cv + M(:,3) = ch + + end function getRayPlotTransformation + end module visualiser_class diff --git a/docs/Input Manual.rst b/docs/Input Manual.rst index 530aa0821..549c7eace 100644 --- a/docs/Input Manual.rst +++ b/docs/Input Manual.rst @@ -296,6 +296,29 @@ Example: :: geometry { } viz { } +.. note:: + SCONE can be run to visualise geometry without requiring an explicit vizPhysicsPackage by + including ``--plot`` when running the application. + +renderPhysicsPackage +#################### + +renderPhysicsPackage, used for performing ray plots which can be modified +directly from the command line. Looks for a visualiser dictionary which contains +a ray plot. Optionally allows the system name for an external visualiser app to +be given. eog is recommended and is the default. + +Example: :: + + type renderPhysicsPackage; + app eog; + geometry { } + viz { } + +.. note:: + SCONE can be run to begin a live render without requiring an explicit renderPhysicsPackage by + including ``--render`` when running the application. + Source ------