Add support for general quadric surfaces

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:
 

docs/source/usersguide/input.rst

@@ -787,7 +787,7 @@ Each ``<surface>`` 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`".
+
+
 ``<cell>`` Element
 ------------------
 

docs/source/usersguide/output/summary.rst

@@ -132,7 +132,8 @@ The current revision of the summary file format is 1.
 **/geometry/surfaces/surface <uid>/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 <uid>/coefficients** (*double[]*)
 

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

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)

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