From c4fd2855e6ff9252003d7b2607570ce6884938c4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 9 Oct 2015 09:20:15 +0700 Subject: [PATCH] Add support for general quadric surfaces --- docs/source/methods/geometry.rst | 4 + docs/source/usersguide/input.rst | 8 +- docs/source/usersguide/output/summary.rst | 3 +- src/input_xml.F90 | 14 +++ src/summary.F90 | 5 ++ src/surface_header.F90 | 102 +++++++++++++++++++++- 6 files changed, 131 insertions(+), 5 deletions(-) diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 9eece1e919..6592909298 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -117,6 +117,10 @@ to fully define the surface. | Cone parallel to the | z-cone | :math:`(x-x_0)^2 + (y-y_0)^2 | :math:`x_0 \; y_0 \; | | :math:`z`-axis | | = R^2(z-z_0)^2` | z_0 \; R^2` | +----------------------+------------+------------------------------+-------------------------+ + | General quadric | quadric | :math:`Ax^2 + By^2 + Cz^2 + | :math:`A \; B \; C \; D | + | surface | | Dxy + Eyz + Fxz + Gx + Hy + | \; E \; F \; G \; H \; | + | | | Jz + K` | J \; K` | + +----------------------+------------+------------------------------+-------------------------+ .. _universes: diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 44db252079..87b613a453 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -787,7 +787,7 @@ Each ```` element can have the following attributes or sub-elements: :type: The type of the surfaces. This can be "x-plane", "y-plane", "z-plane", "plane", "x-cylinder", "y-cylinder", "z-cylinder", "sphere", "x-cone", - "y-cone", or "z-cone". + "y-cone", "z-cone", or "quadric". *Default*: None @@ -855,6 +855,12 @@ The following quadratic surfaces can be modeled: R^2 (z - z_0)^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0 \: R^2`". + :quadric: + A general quadric surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + + Eyz + Fxz + Gx + Hy + Jz + K = 0` The coefficients specified are ":math:`A + \: B \: C \: D \: E \: F \: G \: H \: J \: K`". + + ```` Element ------------------ diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index d37b2172fa..f87f60c4a0 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -132,7 +132,8 @@ The current revision of the summary file format is 1. **/geometry/surfaces/surface /type** (*char[]*) Type of the surface. Can be 'x-plane', 'y-plane', 'z-plane', 'plane', - 'x-cylinder', 'y-cylinder', 'sphere', 'x-cone', 'y-cone', or 'z-cone'. + 'x-cylinder', 'y-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', or + 'quadric'. **/geometry/surfaces/surface /coefficients** (*double[]*) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6d15f0d204..40e0583daf 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1294,6 +1294,9 @@ contains case ('z-cone') coeffs_reqd = 4 allocate(SurfaceZCone :: surfaces(i)%obj) + case ('quadric') + coeffs_reqd = 10 + allocate(SurfaceQuadric :: surfaces(i)%obj) case default call fatal_error("Invalid surface type: " // trim(word)) end select @@ -1378,6 +1381,17 @@ contains s%y0 = coeffs(2) s%z0 = coeffs(3) s%r2 = coeffs(4) + type is (SurfaceQuadric) + s%A = coeffs(1) + s%B = coeffs(2) + s%C = coeffs(3) + s%D = coeffs(4) + s%E = coeffs(5) + s%F = coeffs(6) + s%G = coeffs(7) + s%H = coeffs(8) + s%J = coeffs(9) + s%K = coeffs(10) end select ! No longer need coefficients diff --git a/src/summary.F90 b/src/summary.F90 index 1bf57cc93f..d8ac1e3fa8 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -276,6 +276,11 @@ contains allocate(coeffs(4)) coeffs(:) = [s%x0, s%y0, s%z0, s%r2] + type is (SurfaceQuadric) + call write_dataset(surface_group, "type", "quadric") + allocate(coeffs(10)) + coeffs(:) = [s%A, s%B, s%C, s%D, s%E, s%F, s%G, s%H, s%J, s%K] + end select call write_dataset(surface_group, "coefficients", coeffs) deallocate(coeffs) diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 6fc1bbc237..a927ce4d40 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -1,6 +1,6 @@ module surface_header - use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT + use constants, only: ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT implicit none @@ -183,6 +183,15 @@ module surface_header procedure :: normal => z_cone_normal end type SurfaceZCone + type, extends(Surface) :: SurfaceQuadric + ! Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 + real(8) :: A, B, C, D, E, F, G, H, J, K + contains + procedure :: evaluate => quadric_evaluate + procedure :: distance => quadric_distance + procedure :: normal => quadric_normal + end type SurfaceQuadric + contains !=============================================================================== @@ -640,7 +649,7 @@ contains end function z_cylinder_normal !=============================================================================== -! SphereImplementation +! SurfaceSphere Implementation !=============================================================================== pure function sphere_evaluate(this, xyz) result(f) @@ -874,7 +883,7 @@ contains end function y_cone_normal !=============================================================================== -! SurfaceZConeImplementation +! SurfaceZCone Implementation !=============================================================================== pure function z_cone_evaluate(this, xyz) result(f) @@ -953,4 +962,91 @@ contains uvw(3) = -TWO*this%r2*(xyz(3) - this%z0) end function z_cone_normal +!=============================================================================== +! SurfaceQuadric Implementation +!=============================================================================== + + pure function quadric_evaluate(this, xyz) result(f) + class(SurfaceQuadric), intent(in) :: this + real(8), intent(in) :: xyz(3) + real(8) :: f + + associate (x => xyz(1), y => xyz(2), z => xyz(3)) + f = x*(this%A*x + this%D*y + this%G) + & + y*(this%B*y + this%E*z + this%H) + & + z*(this%C*z + this%F*x + this%J) + this%K + end associate + end function quadric_evaluate + + pure function quadric_distance(this, xyz, uvw, coincident) result(d) + class(SurfaceQuadric), intent(in) :: this + real(8), intent(in) :: xyz(3) + real(8), intent(in) :: uvw(3) + logical, intent(in) :: coincident + real(8) :: d + + real(8) :: a, k, c + real(8) :: quad, b + + associate (x => xyz(1), y => xyz(2), z => xyz(3), & + u => uvw(1), v => uvw(2), w => uvw(3)) + + a = this%A*u*u + this%B*v*v + this%C*w*w + this%D*u*v + this%E*v*w + & + this%F*u*w + k = (this%A*u*x + this%B*v*y + this%C*w*z + HALF*(this%D*(u*y + v*x) + & + this%E*(v*z + w*y) + this%F*(w*x + u*z) + this%G*u + this%H*v + & + this%J*w)) + c = this%A*x*x + this%B*y*y + this%C*z*z + this%D*x*y + this%E*y*z + & + this%F*y*z + this%G*x + this%H*y + this%J*z + K + quad = k*k - a*c + + if (quad < ZERO) then + ! no intersection with cone + + d = INFINITY + + elseif (coincident .or. abs(c) < FP_COINCIDENT) then + ! particle is on the cone, thus one distance is positive/negative and the + ! other is zero. The sign of k determines which distance is zero and which + ! is not. + + if (k >= ZERO) then + d = (-k - sqrt(quad))/a + else + d = (-k + sqrt(quad))/a + end if + + else + ! calculate both solutions to the quadratic + quad = sqrt(quad) + d = (-k - quad)/a + b = (-k + quad)/a + + ! determine the smallest positive solution + if (d < ZERO) then + if (b > ZERO) then + d = b + end if + else + if (b > ZERO) d = min(d, b) + end if + end if + + ! If the distance was negative, set boundary distance to infinity + if (d <= ZERO) d = INFINITY + end associate + end function quadric_distance + + pure function quadric_normal(this, xyz) result(uvw) + class(SurfaceQuadric), intent(in) :: this + real(8), intent(in) :: xyz(3) + real(8) :: uvw(3) + + associate (x => xyz(1), y => xyz(2), z => xyz(3)) + uvw(1) = TWO*this%A*x + this%D*y + this%F*z + this%G + uvw(2) = TWO*this%B*y + this%D*x + this%E*z + this%H + uvw(3) = TWO*this%C*z + this%E*y + this%F*x + this%J + end associate + end function quadric_normal + end module surface_header