Making the constraint calculation aware of LegoExcision.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@53 b7a48df3-cbbf-4440-997f-b4b717c9f7fc

1 files changed, 39 insertions, 0 deletions

diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F
index 2157094..529bdc4 100644
--- a/src/ADMConstraints.F
+++ b/src/ADMConstraints.F
@@ -39,6 +39,10 @@
       integer :: i,j,k
       integer :: nx,ny,nz
 
+#ifdef EXCISION_LEGOEXCISION
+      CCTK_REAL, allocatable, dimension (:,:,:) :: dirx,diry,dirz,aux
+#endif
+
 c     Stencil width used for calculating constraints
 c     (for outer boundary condition)
 
@@ -218,6 +222,41 @@ c                  = - (T_{i0} - beta^j T_{ij})/alpha
 #include "CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h"
 #include "CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h"
 
+c     LegoExcision (must be done before symmetries are applied).
+
+      if (excise==1) then
+
+#ifdef EXCISION_LEGOEXCISION
+
+         allocate(dirx(nx,ny,nz),diry(nx,ny,nz),dirz(nx,ny,nz))
+         allocate(aux(nx,ny,nz))
+
+         aux = 0.0D0
+
+         call excision_findboundary(ierr,emask,nx,ny,nz)
+
+         call CCTK_SyncGroup(cctkGH,"einstein::mask")
+         call CartSymGN(ierr,cctkGH,"einstein::mask")
+
+         call excision_findnormals(ierr,emask,dirx,diry,dirz,nx,ny,nz)
+
+         call excision_extrapolate(ierr,ham ,aux,emask,
+     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
+         call excision_extrapolate(ierr,momx,aux,emask,
+     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
+         call excision_extrapolate(ierr,momy,aux,emask,
+     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
+         call excision_extrapolate(ierr,momz,aux,emask,
+     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
+
+         deallocate(dirx,diry,dirz)
+
+#else
+         call CCTK_WARN(0,"You have not compiled with LegoExcision")
+#endif
+
+      end if
+
 c     Apply symmetry boundary conditions.
 
       call CartSymGN(ierr,cctkGH,"admconstraints::hamiltonian")