New routine to do ''integrals'' over the surface without actually finding

it. Still in testing phase.


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

3 files changed, 259 insertions, 0 deletions

diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90
new file mode 100644
index 0000000..e8ebbf9
--- /dev/null
+++ b/src/EHFinder_Integrate.F90
@@ -0,0 +1,257 @@
+! Integration over the surface
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k
+  CCTK_REAL :: idx, idy, idz
+  CCTK_REAL :: al, ar, bl, br, cl, cr
+  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
+  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
+  CCTK_REAL :: dv, local_area, area, eps, eps_inv
+  CCTK_REAL, dimension(3,3) :: ginv
+  CCTK_REAL :: g2xx, g2xy, g2xz, g2yy, g2yz, g2zz
+  CCTK_REAL :: g22xx, g22xy, g22yy
+  CCTK_REAL, dimension(3) :: norm
+  CCTK_REAL :: dfsquare, dfsquare_i
+  CCTK_REAL, dimension(3) :: vec1, vec2
+  CCTK_REAL :: tmp1, tmp2, tmp3
+  CCTK_REAL :: idetg, sqrtdg2
+
+  interface
+    CCTK_REAL function dirac ( x, eps, eps_inv )
+      CCTK_REAL, intent(in) :: x, eps, eps_inv
+    end function dirac
+  end interface
+
+# include "include/physical_part.h"
+
+  idx = half / cctk_delta_space(1)
+  idy = half / cctk_delta_space(2)
+  idz = half / cctk_delta_space(3)
+
+  dv = cctk_delta_space(1)*cctk_delta_space(2)*cctk_delta_space(3)
+  eps = minval(cctk_delta_space)
+  eps_inv = one / eps
+  
+  local_area = zero
+
+  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+          if ( abs(f(i,j,k)) .lt. eps ) then
+
+#include "include/upwind_second2.h"
+
+            tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
+            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 )
+
+            ginv(1,1) = tmp1 * idetg
+            ginv(1,2) = tmp2 * idetg
+            ginv(1,3) = tmp3 * idetg
+            ginv(2,1) = ginv(1,2)
+            ginv(2,2) = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
+            ginv(2,3) = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) )*idetg
+            ginv(3,1) = ginv(1,3)
+            ginv(3,2) = ginv(2,3)
+            ginv(3,3) = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
+
+            dfsquare = sqrt ( ginv(1,1) * dfx(i,j,k)**2 + &
+                              two * ginv(1,2) * dfx(i,j,k) * dfy(i,j,k) + &
+                              two * ginv(1,3) * dfx(i,j,k) * dfz(i,j,k) + &
+                              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 )
+
+            dfsquare_i = one / dfsquare
+
+            norm(1) = dfx(i,j,k) * dfsquare_i
+            norm(2) = dfy(i,j,k) * dfsquare_i
+            norm(3) = dfz(i,j,k) * dfsquare_i
+
+            call find_orthogonal ( norm, vec1, vec2, ginv, idetg )
+!            print*,norm(1),norm(2),norm(3)
+!            print*,vec1
+!            print*,vec2
+!            print*,norm(1)*vec1(1)+norm(2)*vec1(2)+norm(3)*vec1(3)
+!            print*,norm(1)*vec2(1)+norm(2)*vec2(2)+norm(3)*vec2(3)
+!            print*,gxx(i,j,k)*vec1(1)*vec2(1)+gxy(i,j,k)*vec1(1)*vec2(2)+&
+!                   gxz(i,j,k)*vec1(1)*vec2(3)+gxy(i,j,k)*vec1(2)*vec2(1)+&
+!                   gyy(i,j,k)*vec1(2)*vec2(2)+gyz(i,j,k)*vec1(2)*vec2(3)+&
+!                   gxz(i,j,k)*vec1(3)*vec2(1)+gyz(i,j,k)*vec1(3)*vec2(2)+&
+!                   gzz(i,j,k)*vec1(3)*vec2(3)
+!            print*,gxx(i,j,k)*vec1(1)*vec1(1)+gxy(i,j,k)*vec1(1)*vec1(2)+&
+!                   gxz(i,j,k)*vec1(1)*vec1(3)+gxy(i,j,k)*vec1(2)*vec1(1)+&
+!                   gyy(i,j,k)*vec1(2)*vec1(2)+gyz(i,j,k)*vec1(2)*vec1(3)+&
+!                   gxz(i,j,k)*vec1(3)*vec1(1)+gyz(i,j,k)*vec1(3)*vec1(2)+&
+!                   gzz(i,j,k)*vec1(3)*vec1(3)
+!            print*,gxx(i,j,k)*vec2(1)*vec2(1)+gxy(i,j,k)*vec2(1)*vec2(2)+&
+!                   gxz(i,j,k)*vec2(1)*vec2(3)+gxy(i,j,k)*vec2(2)*vec2(1)+&
+!                   gyy(i,j,k)*vec2(2)*vec2(2)+gyz(i,j,k)*vec2(2)*vec2(3)+&
+!                   gxz(i,j,k)*vec2(3)*vec2(1)+gyz(i,j,k)*vec2(3)*vec2(2)+&
+!                   gzz(i,j,k)*vec2(3)*vec2(3)
+ 
+!            print*,gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
+!                   gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)
+            g2xx = gxx(i,j,k) - norm(1)**2
+            g2xy = gxy(i,j,k) - norm(1) * norm(2)
+            g2xz = gxz(i,j,k) - norm(1) * norm(3)
+            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
+
+            g22xx = vec1(1)**2*ginv(1,1) + two*vec1(1)*vec1(2)*ginv(1,2) +&
+                    two*vec1(1)*vec1(3)*ginv(1,3) + vec1(2)**2*ginv(2,2) +&
+                    two*vec1(2)*vec1(3)*ginv(2,3) + vec1(3)**2 *ginv(3,3)
+            g22xy = vec1(1)*vec2(1)*ginv(1,1) + vec1(1)*vec2(2)*ginv(1,2) +&
+                    vec1(1)*vec2(3)*ginv(1,3) + vec1(2)*vec2(1)*ginv(2,1) +&
+                    vec1(2)*vec2(2)*ginv(2,2) + vec1(2)*vec2(3)*ginv(2,3) +&
+                    vec1(3)*vec2(1)*ginv(3,1) + vec1(3)*vec2(2)*ginv(3,2) +&
+                    vec1(3)*vec2(3)*ginv(3,3)
+            g22yy = vec2(1)**2*ginv(1,1) + two*vec2(1)*vec2(2)*ginv(1,2) +&
+                    two*vec2(1)*vec2(3)*ginv(1,3) + vec2(2)**2*ginv(2,2) +&
+                    two*vec2(2)*vec2(3)*ginv(2,3) + vec2(3)**2*ginv(3,3)
+
+!            print*
+!            print*,g22xx,g22xy,g22yy
+!            print*,g22xx*g22yy-g22xy**2
+
+!          tmp1 = g2yy*g2zz - g2yz**2
+!          tmp2 = g2xz*g2yz - g2xy*g2zz
+!          tmp3 = g2xy*g2yz - g2xz*g2yy        
+!
+!          print*,tmp1,tmp2,tmp3
+!
+!          print*, g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3
+!          sqrtdg2 = sqrt ( g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3 )
+            sqrtdg2 = one / sqrt ( g22xx*g22yy-g22xy**2 )
+!            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
+
+!            print*,sqrtdg2
+!            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 )
+          end if
+        end do
+      end do
+    end do
+  end if
+  local_area = dv * local_area
+
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                                    local_area, area, CCTK_VARIABLE_REAL )
+
+  eh_area = area
+!  print*,'Area = ',eh_area
+
+  return
+end subroutine EHFinder_Integrate
+
+
+subroutine find_orthogonal ( n, v1, v2, ginv, idetg )
+
+  use EHFinder_Constants
+
+  CCTK_REAL, dimension(3) :: n, v1, v2
+  CCTK_REAL, dimension(3) :: v1l
+  CCTK_REAL, dimension(3,3) :: ginv
+  CCTK_REAL :: idetg
+  CCTK_REAL, dimension(1) :: maxn
+  CCTK_REAL :: nv1, nv2
+  CCTK_INT, dimension(1) :: maxp
+  CCTK_INT :: first, second, third
+
+  maxn = maxval (n)
+  maxp = maxloc (n)
+
+!  print*,'Orth'
+!  print*,n
+!  print*,maxn
+!  print*,maxp
+  third = maxp(1)
+  if ( third .eq. 1 ) then
+    first = 2
+    second =3
+  else if ( third .eq. 2 ) then
+    first = 3
+    second = 1
+  else
+    first = 1
+    second = 2
+  end if
+
+!  print*,first,second,third
+  v1l(first) = -n(third)
+  v1l(second) = n(first)
+  v1l(third) = - ( ginv(first,first) * n(first) * v1l(first) + &
+                  ginv(first,second) * n(first) * v1l(second) + &
+                  ginv(second,first) * n(second) * v1l(first) + &
+                  ginv(second,second) * n(second) * v1l(second) + &
+                  ginv(third,first) * n(third) * v1l(first) + &
+                  ginv(third,second) * n(third) * v1l(second) ) / &
+                ( ginv(first,third) * n(first) + &
+                  ginv(second,third) * n(second) + &
+                  ginv(third,third) * n(third) )
+
+  nv1 = sqrt ( ginv(1,1) * v1l(1)**2 + &
+               two * ginv(1,2) * v1l(1) * v1l(2) + &
+               two * ginv(1,3) * v1l(1) * v1l(3) + &
+               ginv(2,2) * v1l(2)**2 + &
+               two * ginv(2,3) * v1l(2) * v1l(3) + &
+               ginv(3,3) * v1l(3)**2 )
+
+  v1l = v1l / nv1
+  v1(1) = ginv(1,1) * v1l(1) + ginv(1,2) * v1l(2) + ginv(1,3) * v1l(3)
+  v1(2) = ginv(2,1) * v1l(1) + ginv(2,2) * v1l(2) + ginv(2,3) * v1l(3)
+  v1(3) = ginv(3,1) * v1l(1) + ginv(3,2) * v1l(2) + ginv(3,3) * v1l(3)
+
+!  print*, v1(1)*v1l(1)+v1(2)*v1l(2)+v1(3)*v1l(3)
+!  print*
+
+  v2(1) = n(2) * v1(3) - n(3) * v1(2)
+  v2(2) = n(3) * v1(1) - n(1) * v1(3)
+  v2(3) = n(1) * v1(2) - n(2) * v1(1)
+  v2 = v2 * sqrt(idetg)
+
+!  print*, v2
+!  print*
+!  print*,v1(1)*n(1)+v1(2)*n(2)+v1(3)*n(3)
+!  print*,v2(1)*n(1)+v2(2)*n(2)+v2(3)*n(3)
+!  print*,v2(1)*v1l(1)+v2(2)*v1l(2)+v2(3)*v1l(3)
+!  print*,'Orth end'
+!  pause 
+  return
+end subroutine find_orthogonal
+  
+CCTK_REAL function dirac ( x, eps, eps_inv )
+
+  use EHFinder_Constants
+
+  CCTK_REAL, intent(in) :: x, eps, eps_inv
+
+  if ( abs(x) .lt. eps ) then
+    dirac = half * ( one + cos(pi*x*eps_inv) ) * eps_inv
+  else
+    dirac = zero
+  end if
+end function dirac
diff --git a/src/make.code.defn b/src/make.code.defn
index 33d3132..8f9c8c8 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -14,6 +14,7 @@ SRCS = EHFinder_mod.F90 \
        EHFinder_SetMask.F90 \
        EHFinder_SetSym.F90 \
        EHFinder_Check.F90 \
+       EHFinder_Integrate.F90 \
        EHFinder_ReadData.F90
 
 # Subdirectories containing source files
diff --git a/src/make.code.deps b/src/make.code.deps
index 907b00e..07978ca 100644
--- a/src/make.code.deps
+++ b/src/make.code.deps
@@ -12,3 +12,4 @@ EHFinder_SetMask.F90.o:		EHFinder_mod.F90.o
 EHFinder_SetSym.F90.o:		EHFinder_mod.F90.o
 EHFinder_ReadData.F90.o:	EHFinder_mod.F90.o
 EHFinder_Check.F90.o:		EHFinder_mod.F90.o
+EHFinder_Integrate.F90.o:	EHFinder_mod.F90.o