Major modification of almost all source files to use vector groups for various

variables. Not tested in all detail, but the standard features seem to work.
The version before all these changes was tagged with PRE_MULTI.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@126 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 90 insertions, 75 deletions

diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index 2bf3d7d..8fc9587 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -15,7 +15,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k, status
+  CCTK_INT :: i, j, k, l, status
   CCTK_REAL, dimension(3) :: xp, xpt
   CCTK_REAL, dimension(3,3) :: txyz
   CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc
@@ -35,13 +35,25 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
   cctk_time = last_time
 
+  ! Allocate the logical array containing the flag determining if the
+  ! corresponding level set should be re-initialized.
+
+  if ( allocated(reparam_this_level_set) ) then
+    deallocate ( reparam_this_level_set )
+  end if
+  allocate ( reparam_this_level_set(eh_number_level_sets) )
+  if ( allocated(reparam_undone) ) then
+    deallocate ( reparam_undone )
+  end if
+  allocate ( reparam_undone(eh_number_level_sets) )
+
   if ( evolve_generators .gt. 0 ) then
 
-    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" )
+    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::xg" )
     if ( status .lt. 0 ) then
       call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" )
     end if
-    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" )
     if ( status .lt. 0 ) then
       call CCTK_WARN ( 0, "cannot get local size for generator arrays" )
     end if
@@ -75,90 +87,93 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
     end if
   end if
 
-  ! If a sphere is requested...
-  if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then
-
-    ! Set up a sphere of radius initial_rad and translated 
-    ! (translate_x,translate_y,translate_z) away from the origin.
-    f = sqrt( ( x - translate_x )**2 + &
-              ( y - translate_y )**2 + &
-              ( z - translate_z )**2 ) - initial_rad
-    if ( evolve_generators .gt. 0 ) then
-      if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
-        do i = 1, lsh(1)
-           theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-           xg(i) = initial_rad * sin(theta) + translate_x
-           yg(i) = translate_y
-           zg(i) = initial_rad * cos(theta) + translate_z
-        end do
+  do l = 1, eh_number_level_sets
+
+    ! If a sphere is requested...
+    if ( CCTK_EQUALS( initial_f(l), 'sphere' ) ) then
+
+      ! Set up a sphere of radius initial_rad and translated 
+      ! (translate_x,translate_y,translate_z) away from the origin.
+      f(:,:,:,l) = sqrt( ( x - translate_x )**2 + &
+                         ( y - translate_y )**2 + &
+                         ( z - translate_z )**2 ) - initial_rad(l)
+      if ( evolve_generators .gt. 0 ) then
+        if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
+          do i = 1, lsh(1)
+             theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
+             xg(i,l) = initial_rad(l) * sin(theta) + translate_x
+             yg(i,l) = translate_y
+             zg(i,l) = initial_rad(l) * cos(theta) + translate_z
+          end do
+        end if
       end if
     end if
-  end if
 
-  ! If an ellipsoid is requested...
-  if ( CCTK_EQUALS( initial_f, 'ellipsoid' ) ) then
-
-    ! Calculate sines and cosines of the rotation parameters.
-    cosa = cos(rotation_alpha)
-    sina = sin(rotation_alpha)
-    cosb = cos(rotation_beta)
-    sinb = sin(rotation_beta)
-    cosc = cos(rotation_gamma)
-    sinc = sin(rotation_gamma)
-
-    ! Set up the rotation matrix. The order is alpha around the z-axis,
-    ! beta around the y-axis and finally gamma around the x-axis.
-    txyz(1,1) = cosa * cosb
-    txyz(1,2) = sina * cosb
-    txyz(1,3) = -sinb
-    txyz(2,1) = cosa * sinb * sinc - sina * cosc
-    txyz(2,2) = sina * sinb * sinc + cosa * cosc
-    txyz(2,3) = cosb * sinc
-    txyz(3,1) = cosa * sinb * cosc + sina * sinc
-    txyz(3,2) = sina * sinb * cosc - cosa * sinc
-    txyz(3,3) = cosb * cosc
-!    print*,txyz
+    ! If an ellipsoid is requested...
+    if ( CCTK_EQUALS( initial_f(l), 'ellipsoid' ) ) then
+
+      ! Calculate sines and cosines of the rotation parameters.
+      cosa = cos(rotation_alpha)
+      sina = sin(rotation_alpha)
+      cosb = cos(rotation_beta)
+      sinb = sin(rotation_beta)
+      cosc = cos(rotation_gamma)
+      sinc = sin(rotation_gamma)
+
+      ! Set up the rotation matrix. The order is alpha around the z-axis,
+      ! beta around the y-axis and finally gamma around the x-axis.
+      txyz(1,1) = cosa * cosb
+      txyz(1,2) = sina * cosb
+      txyz(1,3) = -sinb
+      txyz(2,1) = cosa * sinb * sinc - sina * cosc
+      txyz(2,2) = sina * sinb * sinc + cosa * cosc
+      txyz(2,3) = cosb * sinc
+      txyz(3,1) = cosa * sinb * cosc + sina * sinc
+      txyz(3,2) = sina * sinb * cosc - cosa * sinc
+      txyz(3,3) = cosb * cosc
+!      print*,txyz
 
     ! Apply the rotations and translation for all points on the grid.
     ! Even though at first glance it looks like the translation is done
     ! first, the opposite is actually true.
-    do k = 1, nz
-      do j = 1, ny
-        do i = 1, nx
-          xp(1) = x(i,j,k) - translate_x
-          xp(2) = y(i,j,k) - translate_y
-          xp(3) = z(i,j,k) - translate_z
-          xpt = matmul ( txyz, xp )
-          f(i,j,k) = sqrt( xpt(1)**2 / initial_a**2 + &
-                           xpt(2)**2 / initial_b**2 + &
-                           xpt(3)**2 / initial_c**2) - 1.0
+      do k = 1, nz
+        do j = 1, ny
+          do i = 1, nx
+            xp(1) = x(i,j,k) - translate_x
+            xp(2) = y(i,j,k) - translate_y
+            xp(3) = z(i,j,k) - translate_z
+            xpt = matmul ( txyz, xp )
+            f(i,j,k,l) = sqrt( xpt(1)**2 / initial_a**2 + &
+                             xpt(2)**2 / initial_b**2 + &
+                             xpt(3)**2 / initial_c**2) - 1.0
+          end do
         end do
       end do
-    end do
-
-    if ( evolve_generators .gt. 0 ) then
-      if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
-        do i = 1, lsh(1)
-          theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
-          r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + &
-                         cos(theta)**2 / initial_c**2 ) )
-          xp(1) = r_el * sin(theta) + translate_x
-          xp(2) = translate_y
-          xp(3) = r_el * cos(theta) + translate_z
-          xpt = matmul ( txyz, xp )
-          xg(i) = xpt(1)
-          yg(i) = xpt(2)
-          zg(i) = xpt(3)
-        end do
+
+      if ( evolve_generators .gt. 0 ) then
+        if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
+          do i = 1, lsh(1)
+            theta = thetamin + dtheta * ( i + lbnd(1) - 1 )
+            r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + &
+                           cos(theta)**2 / initial_c**2 ) )
+            xp(1) = r_el * sin(theta) + translate_x
+            xp(2) = translate_y
+            xp(3) = r_el * cos(theta) + translate_z
+            xpt = matmul ( txyz, xp )
+            xg(i,l) = xpt(1)
+            yg(i,l) = xpt(2)
+            zg(i,l) = xpt(3)
+          end do
+        end if
       end if
     end if
-  end if
 
-  ! if an ovaloid of Cassini is requested...
-  if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
-    f = (x**2+y**2+z**2)**2 + cas_a**4 - &
-        2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4
-  end if
+    ! if an ovaloid of Cassini is requested...
+    if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
+      f(:,:,:,l) = (x**2+y**2+z**2)**2 + cas_a**4 - &
+            2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4
+    end if
+  end do
 
   ! Initialise the internal mask.
   eh_mask = 0