Added support for finding surface area. However only for one surface and

in full mode. It is parallelized though. Need to put in a routine to
detect how many surface are present, their symmetries and handle various
symmetries.


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

1 files changed, 30 insertions, 13 deletions

diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90
index c46a43a..04a2a46 100644
--- a/src/EHFinder_Integrate.F90
+++ b/src/EHFinder_Integrate.F90
@@ -28,7 +28,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
   CCTK_REAL :: dfsquare, dfsquare_i
   CCTK_REAL, dimension(3) :: vec1, vec2
   CCTK_REAL :: tmp1, tmp2, tmp3
-  CCTK_REAL :: idetg, sqrtdg2
+  CCTK_REAL :: detg, idetg, sqrtdg2
   CCTK_REAL :: a, b, c, r1, r2, r3, q
 
   interface
@@ -63,8 +63,8 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
             tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
   
-            idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + &
-                                              gxz(i,j,k)*tmp3 )
+            detg = gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3
+            idetg = one / detg
 
 ! Invert the spatial 3-metric
 
@@ -84,8 +84,10 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
                               ginv(2,2) * dfy(i,j,k)**2 + &
                               two * ginv(2,3) * dfy(i,j,k)*dfz(i,j,k) + &
                               ginv(3,3) *  dfz(i,j,k)**2 )
+!            print*,dfsquare
 
             dfsquare_i = one / dfsquare
+!            print*,dfsquare_i
 !            dfsquare_i = one / sqrt ( dfx(i,j,k)**2 + &
 !                                      dfy(i,j,k)**2 + &
 !                                      dfz(i,j,k)**2 )
@@ -93,6 +95,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             norm(1) = dfx(i,j,k) * dfsquare_i
             norm(2) = dfy(i,j,k) * dfsquare_i
             norm(3) = dfz(i,j,k) * dfsquare_i
+!            print*,norm
 
 !            norm(1) = dfx(i,j,k)
 !            norm(2) = dfy(i,j,k)
@@ -131,6 +134,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
             g2yy = gyy(i,j,k) - norm(2)**2
             g2yz = gyz(i,j,k) - norm(2) * norm(3)
             g2zz = gzz(i,j,k) - norm(3)**2
+!            print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
 
 !            a = - ( g2xx + g2yy + g2zz )
             b = g2xx * g2yy + g2xx * g2zz + g2yy * g2zz - &
@@ -139,11 +143,14 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
 !                  g2xx * g2yy * g2zz - two * g2xy * g2xz * g2yz
 
 !            q = - half * ( a + sign ( one, a ) * sqrt ( a**2 - 4 * b ) )
-!            print*,q,b/q,sqrt(b)
-            sqrtdg2 = sqrt ( sqrt (b) )
-!            pause
+!            print*,'factors : ',a,b,c,q,b/q,sqrt(b)
+!            print*
+!            sqrtdg2 = sqrt ( sqrt (b) )
+            sqrtdg2 = sqrt (b)
 !            call cubic ( a, b, c, r1, r2, r3 )
-!           print*
+!            print*,'roots : ',r1,r2,r3
+!            print*
+!            pause
 !           print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
 !           print*
 
@@ -189,6 +196,7 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
           tmp1 = g2yy*g2zz - g2yz**2
           tmp2 = g2xz*g2yz - g2xy*g2zz
           tmp3 = g2xy*g2yz - g2xz*g2yy        
+!          print*,tmp1,tmp2,tmp3
 
 !          print*,tmp1,tmp2,tmp3
 !
@@ -198,23 +206,25 @@ subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
 !            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
 !            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
 
-!           print*,sqrtdg2 !,g22xx,g22xy,g22yy
+!           print*,sqrtdg2,g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3
 !           pause
 
             local_area = local_area + sqrtdg2 * dirac(f(i,j,k),eps,eps_inv) * &
-                         sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
+                                      dfsquare * sqrt(detg)
+!                         sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
           end if
         end do
       end do
     end do
   end if
   local_area = dv * local_area
+!  print*,'local_area = ',local_area,sqrtdg2
 
   call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
                                     local_area, area, CCTK_VARIABLE_REAL )
 
   eh_area = area
-!  print*,'Area = ',eh_area
+  print*,'Area = ',eh_area
 
   return
 end subroutine EHFinder_Integrate
@@ -312,12 +322,19 @@ subroutine cubic ( a, b, c, x1, x2, x3 )
 
   CCTK_REAL, intent(in) :: a, b, c
   CCTK_REAL, intent(out) :: x1, x2, x3
-  CCTK_REAL :: q, r
+  CCTK_REAL :: q, r, theta
+  CCTK_REAL, parameter :: third = 1.0d0 / 3.0d0, &
+                          pi2 = 6.2831853071795864769d0
 
   q = ( a**2 - 3.0d0 * b ) / 9.0d0
   r = ( 2.0d0 * a**3 - 9.0d0 * a * b + 27.0d0 * c ) / 54.0d0
+  theta = acos ( r / sqrt( q**3 ) )
+
+  x1 = - 2.0d0 * sqrt( q ) * cos( third * theta ) - third * a
+  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta + pi2 ) ) - third * a 
+  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta - pi2 ) ) - third * a 
 
   print*,a,b,c
-  print*,q,r
-  pause
+  print*,q,r,q**3-r**2
+!  pause
 end subroutine cubic