use array dummy args for skew, more tests

autotest/TestGeomUtil.f90

@@ -22,7 +22,8 @@ subroutine collect_geomutil(testsuite)
                              test_point_in_polygon_tri), &
                 new_unittest("point_in_polygon_irr", &
                              test_point_in_polygon_irr), &
-                new_unittest("skew", test_skew) &
+                new_unittest("skew", test_skew), &
+                new_unittest("transform", test_transform) &
                 ]
   end subroutine collect_geomutil
 
@@ -292,6 +293,7 @@ subroutine test_point_in_polygon_irr(error)
     deallocate (face_pts)
   end subroutine test_point_in_polygon_irr
 
+  !> @brief Test 2D skew
   subroutine test_skew(error)
     type(error_type), allocatable, intent(out) :: error
     real(DP) :: v(2)
@@ -320,4 +322,9 @@ subroutine test_skew(error)
     call check(error, v(1) == -2.0_DP .and. v(2) == -2.0_DP, to_string(v(1)))
   end subroutine test_skew
 
+  subroutine test_transform(error)
+    type(error_type), allocatable, intent(out) :: error
+    ! todo
+  end subroutine test_transform
+
 end module TestGeomUtil

src/Solution/ParticleTracker/MethodSubcellTernary.f90

@@ -230,7 +230,7 @@ subroutine track_sub(this, subcell, particle, tmax)
       alp = 0d0
     end if
     loc = (/alp, bet/)
-    if (lbary) loc = skew(loc, sxx, sxy, syy, invert=.true.)
+    if (lbary) loc = skew(loc, (/sxx, sxy, syy/), invert=.true.)
     rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
     res = matmul(rot, loc) ! rotate vector
     x = res(1) + x0

src/Solution/ParticleTracker/Particle.f90

@@ -3,7 +3,7 @@ module ParticleModule
   use KindModule, only: DP, I4B, LGP
   use ConstantsModule, only: DZERO, DONE, LENMEMPATH, LENBOUNDNAME
   use GlobalDataModule
-  use GeomUtilModule, only: transform_coords, modify_transf
+  use GeomUtilModule, only: transform, compose
   implicit none
 
   private
@@ -318,17 +318,17 @@ subroutine transf_coords(this, xorigin, yorigin, zorigin, &
     logical(LGP), intent(in), optional :: invert !< whether to invert
 
     ! -- Transform particle's coordinates
-    call transform_coords(this%x, this%y, this%z, &
-                          this%x, this%y, this%z, &
-                          xorigin, yorigin, zorigin, &
-                          sinrot, cosrot, invert)
+    call transform(this%x, this%y, this%z, &
+                   this%x, this%y, this%z, &
+                   xorigin, yorigin, zorigin, &
+                   sinrot, cosrot, invert)
 
     ! -- Modify transformation from model coordinates to particle's new
     ! -- local coordinates by incorporating this latest transformation
-    call modify_transf(this%xorigin, this%yorigin, this%zorigin, &
-                       this%sinrot, this%cosrot, &
-                       xorigin, yorigin, zorigin, &
-                       sinrot, cosrot, invert)
+    call compose(this%xorigin, this%yorigin, this%zorigin, &
+                 this%sinrot, this%cosrot, &
+                 xorigin, yorigin, zorigin, &
+                 sinrot, cosrot, invert)
 
     ! -- Set isTransformed flag to true. Note that there is no check
     ! -- to see whether the modification brings the coordinates back
@@ -361,9 +361,9 @@ subroutine get_model_coords(this, x, y, z)
 
     if (this%transformed) then
       ! -- Transform back from local to model coordinates
-      call transform_coords(this%x, this%y, this%z, x, y, z, &
-                            this%xorigin, this%yorigin, this%zorigin, &
-                            this%sinrot, this%cosrot, .true.)
+      call transform(this%x, this%y, this%z, x, y, z, &
+                     this%xorigin, this%yorigin, this%zorigin, &
+                     this%sinrot, this%cosrot, .true.)
     else
       ! -- Already in model coordinates
       x = this%x

src/Solution/ParticleTracker/TernarySolveTrack.f90

@@ -204,12 +204,12 @@ subroutine canonical(x0, y0, x1, y1, x2, y2, &
       alp1 = 1d0
       alp2 = 0d0
       bet2 = 1d0
-      cv0 = skew(cv0, sxx, sxy, syy)
-      cv1 = skew(cv1, sxx, sxy, syy)
-      cv2 = skew(cv2, sxx, sxy, syy)
+      cv0 = skew(cv0, (/sxx, sxy, syy/))
+      cv1 = skew(cv1, (/sxx, sxy, syy/))
+      cv2 = skew(cv2, (/sxx, sxy, syy/))
       ! todo turn alpi/beti into array
       res = (/alpi, beti/)
-      res = skew(res, sxx, sxy, syy)
+      res = skew(res, (/sxx, sxy, syy/))
       alpi = res(1)
       beti = res(2)
     end if

src/Utilities/GeomUtil.f90

@@ -143,6 +143,36 @@ subroutine get_jk(nodenumber, ncpl, nlay, icpl, ilay)
     end if
   end subroutine get_jk
 
+  !> @brief Skew a 2D vector along the x-axis.
+  pure function skew(v, s, invert) result(res)
+    ! -- dummy
+    real(DP), intent(in) :: v(2) !< vector
+    real(DP), intent(in) :: s(3) !< skew matrix entries (top left, top right, bottom left)
+    logical(LGP), intent(in), optional :: invert
+    real(DP) :: res(2)
+    ! -- local
+    logical(LGP) :: linvert
+    real(DP) :: sxx, sxy, syy
+
+    ! -- process optional arguments
+    if (present(invert)) then
+      linvert = invert
+    else
+      linvert = .false.
+    end if
+
+    sxx = s(1)
+    sxy = s(2)
+    syy = s(3)
+    if (.not. linvert) then
+      res(1) = sxx * v(1) + sxy * v(2)
+      res(2) = syy * v(2)
+    else
+      res(2) = v(2) / syy
+      res(1) = (v(1) - sxy * res(2)) / sxx
+    end if
+  end function skew
+
   !> @brief Apply a 3D translation and optional 2D rotation to coordinates.
   subroutine transform(xin, yin, zin, &
                        xout, yout, zout, &
@@ -334,228 +364,4 @@ subroutine defaults(xorigin, yorigin, zorigin, &
     if (present(invert_opt)) invert = invert_opt
   end subroutine defaults
 
-  !> @brief Skew a 2D vector (as an array, in-place) along the x-axis.
-  !!
-  !! Skew matrix entries sxx (upper left), sxy (upper right), and syy
-  !! (lower right) must be provided, and a 2-item array of the vector
-  !! components to skew. The transform can optionally be inverted.
-  !<
-  pure function skew(vec, sxx, sxy, syy, invert) result(res)
-    ! -- dummy
-    real(DP), intent(in) :: vec(2)
-    real(DP), intent(in) :: sxx, sxy, syy
-    logical(LGP), intent(in), optional :: invert
-    real(DP) :: res(2)
-    ! -- local
-    logical(LGP) :: linvert
-    real(DP) :: old(2)
-
-    ! -- process optional arguments
-    if (present(invert)) then
-      linvert = invert
-    else
-      linvert = .false.
-    end if
-
-    old = vec
-    if (.not. linvert) then
-      res(1) = sxx * old(1) + sxy * old(2)
-      res(2) = syy * old(2)
-    else
-      res(2) = old(2) / syy
-      res(1) = (old(1) - sxy * res(2)) / sxx
-    end if
-  end function skew
-
-  !> @brief Apply a 3D translation and optional 2D rotation to coordinates.
-  subroutine transform_coords(xin, yin, zin, &
-                              xout, yout, zout, &
-                              xorigin, yorigin, zorigin, &
-                              sinrot, cosrot, &
-                              invert)
-    ! -- dummy
-    real(DP) :: xin, yin, zin !< input coordinates
-    real(DP) :: xout, yout, zout !< output coordinates
-    real(DP), optional :: xorigin, yorigin, zorigin !< origin coordinates
-    real(DP), optional :: sinrot, cosrot !< sine and cosine of rotation
-    logical(LGP), optional :: invert !< whether to invert
-    ! -- local
-    logical(LGP) :: ltranslate, lrotate, linvert
-    real(DP) :: x, y
-    real(DP) :: lxorigin, lyorigin, lzorigin
-    real(DP) :: lsinrot, lcosrot
-
-    ! -- Process option arguments and set defaults and flags
-    call transf_opt_args_prep(lxorigin, lyorigin, lzorigin, &
-                              lsinrot, lcosrot, linvert, &
-                              ltranslate, lrotate, &
-                              xorigin, yorigin, zorigin, &
-                              sinrot, cosrot, invert)
-
-    ! -- Apply transformation or its inverse
-    if (.not. linvert) then
-      ! -- Apply transformation to coordinates
-      if (ltranslate) then
-        xout = xin - lxorigin
-        yout = yin - lyorigin
-        zout = zin - lzorigin
-      else
-        xout = lxorigin
-        yout = lyorigin
-        zout = lzorigin
-      end if
-      if (lrotate) then
-        x = xout
-        y = yout
-        xout = x * lcosrot + y * lsinrot
-        yout = -x * lsinrot + y * lcosrot
-      end if
-    else
-      ! -- Apply inverse of transformation to coordinates
-      if (lrotate) then
-        x = xin * lcosrot - yin * lsinrot
-        y = xin * lsinrot + yin * lcosrot
-      else
-        x = xin
-        y = yin
-      end if
-      if (ltranslate) then
-        xout = x + lxorigin
-        yout = y + lyorigin
-        zout = zin + lzorigin
-      end if
-    end if
-  end subroutine transform_coords
-
-  !> @brief Apply a 3D translation and 2D rotation to an existing transformation.
-  subroutine modify_transf(xorigin, yorigin, zorigin, &
-                           sinrot, cosrot, &
-                           xorigin_new, yorigin_new, zorigin_new, &
-                           sinrot_new, cosrot_new, &
-                           invert)
-    ! -- dummy
-    real(DP) :: xorigin, yorigin, zorigin !< origin coordinates (original)
-    real(DP) :: sinrot, cosrot !< sine and cosine of rotation (original)
-    real(DP), optional :: xorigin_new, yorigin_new, zorigin_new !< origin coordinates (new)
-    real(DP), optional :: sinrot_new, cosrot_new !< sine and cosine of rotation (new)
-    logical(LGP), optional :: invert !< whether to invert
-    ! -- local
-    logical(LGP) :: ltranslate, lrotate, linvert
-    real(DP) :: xorigin_add, yorigin_add, zorigin_add
-    real(DP) :: sinrot_add, cosrot_add
-    real(DP) :: x0, y0, z0, s0, c0
-
-    ! -- Process option arguments and set defaults and flags
-    call transf_opt_args_prep(xorigin_add, yorigin_add, zorigin_add, &
-                              sinrot_add, cosrot_add, linvert, &
-                              ltranslate, lrotate, &
-                              xorigin_new, yorigin_new, zorigin_new, &
-                              sinrot_new, cosrot_new, invert)
-
-    ! -- Copy existing transformation into working copy
-    x0 = xorigin
-    y0 = yorigin
-    z0 = zorigin
-    s0 = sinrot
-    c0 = cosrot
-
-    ! -- Modify transformation
-    if (.not. linvert) then
-      ! -- Apply additional transformation to existing transformation
-      if (ltranslate) then
-        ! -- Calculate modified origin, XOrigin + R^T XOrigin_add, where
-        ! -- XOrigin and XOrigin_add are the existing and additional origin
-        ! -- vectors, respectively, and R^T is the transpose of the existing
-        ! -- rotation matrix
-        call transform_coords(xorigin_add, yorigin_add, zorigin_add, &
-                              xorigin, yorigin, zorigin, &
-                              x0, y0, z0, s0, c0, .true.)
-      end if
-      if (lrotate) then
-        ! -- Calculate modified rotation matrix (represented by sinrot
-        ! -- and cosrot) as R_add R, where R and R_add are the existing
-        ! -- and additional rotation matrices, respectively
-        sinrot = cosrot_add * s0 + sinrot_add * c0
-        cosrot = cosrot_add * c0 - sinrot_add * s0
-      end if
-    else
-      ! -- Apply inverse of additional transformation to existing transformation
-      !
-      ! -- Calculate modified origin, R^T (XOrigin + R_add XOrigin_add), where
-      ! -- XOrigin and XOrigin_add are the existing and additional origin
-      ! -- vectors, respectively, R^T is the transpose of the existing rotation
-      ! -- matrix, and R_add is the additional rotation matrix
-      if (ltranslate) then
-        call transform_coords(-xorigin_add, -yorigin_add, zorigin_add, &
-                              x0, y0, z0, xorigin, yorigin, zorigin, &
-                              -sinrot_add, cosrot_add, .true.)
-      end if
-      xorigin = c0 * x0 - s0 * y0
-      yorigin = s0 * x0 + c0 * y0
-      zorigin = z0
-      if (lrotate) then
-        ! -- Calculate modified rotation matrix (represented by sinrot
-        ! -- and cosrot) as R_add^T R, where R and R_add^T are the existing
-        ! -- rotation matirx and the transpose of the additional rotation
-        ! -- matrix, respectively
-        sinrot = cosrot_add * s0 - sinrot_add * c0
-        cosrot = cosrot_add * c0 + sinrot_add * s0
-      end if
-    end if
-  end subroutine modify_transf
-
-  !> @brief Process optional arguments and set defaults/flags for transformations.
-  subroutine transf_opt_args_prep(xorigin, yorigin, zorigin, &
-                                  sinrot, cosrot, &
-                                  invert, translate, rotate, &
-                                  xorigin_opt, yorigin_opt, zorigin_opt, &
-                                  sinrot_opt, cosrot_opt, invert_opt)
-    ! -- dummy
-    real(DP) :: xorigin, yorigin, zorigin
-    real(DP) :: sinrot, cosrot
-    logical(LGP) :: invert, translate, rotate
-    real(DP), optional :: xorigin_opt, yorigin_opt, zorigin_opt
-    real(DP), optional :: sinrot_opt, cosrot_opt
-    logical(LGP), optional :: invert_opt
-
-    translate = .false.
-    xorigin = DZERO
-    if (present(xorigin_opt)) then
-      xorigin = xorigin_opt
-      translate = .true.
-    end if
-    yorigin = DZERO
-    if (present(yorigin_opt)) then
-      yorigin = yorigin_opt
-      translate = .true.
-    end if
-    zorigin = DZERO
-    if (present(zorigin_opt)) then
-      zorigin = zorigin_opt
-      translate = .true.
-    end if
-    rotate = .false.
-    sinrot = DZERO
-    cosrot = DONE
-    if (present(sinrot_opt)) then
-      sinrot = sinrot_opt
-      if (present(cosrot_opt)) then
-        cosrot = cosrot_opt
-      else
-        ! -- If sinrot_opt is specified but cosrot_opt is not,
-        ! -- default to corresponding non-negative cosrot_add
-        cosrot = dsqrt(DONE - sinrot * sinrot)
-      end if
-      rotate = .true.
-    else if (present(cosrot_opt)) then
-      cosrot = cosrot_opt
-      ! -- cosrot_opt is specified but sinrot_opt is not, so
-      ! -- default to corresponding non-negative sinrot_add
-      sinrot = dsqrt(DONE - cosrot * cosrot)
-      rotate = .true.
-    end if
-    invert = .false.
-    if (present(invert_opt)) invert = invert_opt
-  end subroutine transf_opt_args_prep
-
 end module GeomUtilModule