Apply boundary condition after setting mask to a fixed region.

Allow two fixed regions.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@43 f75ba9e5-694f-0410-ac2c-87ea7ce7132b

5 files changed, 99 insertions, 15 deletions

diff --git a/interface.ccl b/interface.ccl
index 7010b37..792f43e 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -39,3 +39,11 @@ SUBROUTINE ExcisionFindNormals(CCTK_INT OUT ierr, \
 			       CCTK_INT IN nk)
 PROVIDES FUNCTION ExcisionFindNormals WITH excision_findnormals \
                   LANGUAGE Fortran
+
+
+
+CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER IN GH, \
+  CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \
+  CCTK_STRING IN group_name, CCTK_STRING IN bc_name)
+
+USES FUNCTION Boundary_SelectGroupForBC
diff --git a/param.ccl b/param.ccl
index 20176da..10f914c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -8,6 +8,17 @@ KEYWORD fixed_excision "Excise a fixed region?"
   "none"   :: "No fixed excision"
 } "none"
 
+
+
+INT num_fixed_regions "Number of fixed excision regions"
+{
+  0:* :: ""
+} 1
+
+
+
+# First fixed excision region
+
 REAL fixed_size "The size (radius or corner location) of the fixed region"
 {
   0:* :: "Positive"
@@ -27,3 +38,27 @@ REAL fixed_origin_z "The z coordinate of the centre of the fixed region"
 {
   *:* :: ""
 } 0.0
+
+
+
+# Second fixed excision region
+
+REAL fixed2_size "The size (radius or corner location) of the second fixed region"
+{
+  0:* :: "Positive"
+} 1.0
+
+REAL fixed2_origin_x "The x coordinate of the centre of the second fixed region"
+{
+  *:* :: ""
+} 0.0
+
+REAL fixed2_origin_y "The y coordinate of the centre of the second fixed region"
+{
+  *:* :: ""
+} 0.0
+
+REAL fixed2_origin_z "The z coordinate of the centre of the second fixed region"
+{
+  *:* :: ""
+} 0.0
diff --git a/schedule.ccl b/schedule.ccl
index 169e015..961b7e0 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,5 +4,9 @@
 schedule Lego_FixedSphere in CCTK_POSTINITIAL
 {
   LANG:Fortran
+  SYNC:SpaceMask::mask
 } "Excise a fixed lego sphere"
 
+schedule GROUP ApplyBCs as Lego_ApplyBCs in CCTK_POSTINITIAL after Lego_FixedSphere
+{
+} "Apply boundary conditions controlled by thorn Boundary"
diff --git a/src/fixedsphere.F90 b/src/fixedsphere.F90
index dafca7d..fd38a6a 100644
--- a/src/fixedsphere.F90
+++ b/src/fixedsphere.F90
@@ -1,3 +1,4 @@
+! $Header$
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -12,29 +13,60 @@ subroutine Lego_FixedSphere(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
   
+  CCTK_INT, parameter :: ione = 1
+  CCTK_INT, parameter :: faces=CCTK_ALL_FACES
+  
   integer ierr
 
   if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then
      call CCTK_WARN (0, "The routine Lego_FixedSphere was called, but thorn " // CCTK_THORNSTRING // " is not active")
   end if
 
-  ! Only excise if there would be more than a single grid point excised
-  if (fixed_size < 2*minval(CCTK_DELTA_SPACE(:))) return
-
-  if (CCTK_EQUALS(fixed_excision,"sphere")) then
-    where ((x - fixed_origin_x)**2 + (y - fixed_origin_y)**2 + (z - fixed_origin_z)**2 < fixed_size**2)
-      emask = 0.d0
-    elsewhere
-      emask = 1.d0
-    end where
-  else if (CCTK_EQUALS(fixed_excision,"cube")) then
-    where (max(abs(x - fixed_origin_x), abs(y - fixed_origin_y), abs(z - fixed_origin_z)) < fixed_size)
-      emask = 0.d0
-    elsewhere
-      emask = 1.d0
-    end where
+  ! Default: no excision
+  emask = 1
+
+  ! First excision region
+  if (num_fixed_regions >= 1) then
+
+     ! Only excise if there would be more than a single grid point excised
+     if (fixed_size >= 2*minval(CCTK_DELTA_SPACE(:))) then
+
+        if (CCTK_EQUALS(fixed_excision,"sphere")) then
+           where ((x - fixed_origin_x)**2 + (y - fixed_origin_y)**2 + (z - fixed_origin_z)**2 < fixed_size**2)
+              emask = 0
+           end where
+        else if (CCTK_EQUALS(fixed_excision,"cube")) then
+           where (max(abs(x - fixed_origin_x), abs(y - fixed_origin_y), abs(z - fixed_origin_z)) < fixed_size)
+              emask = 0
+           end where
+        end if
+
+     end if
+
+  end if
+
+  ! Second excision region
+  if (num_fixed_regions >= 2) then
+
+     ! Only excise if there would be more than a single grid point excised
+     if (fixed2_size >= 2*minval(CCTK_DELTA_SPACE(:))) then
+
+        if (CCTK_EQUALS(fixed_excision,"sphere")) then
+           where ((x - fixed2_origin_x)**2 + (y - fixed2_origin_y)**2 + (z - fixed2_origin_z)**2 < fixed2_size**2)
+              emask = 0
+           end where
+        else if (CCTK_EQUALS(fixed_excision,"cube")) then
+           where (max(abs(x - fixed2_origin_x), abs(y - fixed2_origin_y), abs(z - fixed2_origin_z)) < fixed2_size)
+              emask = 0
+           end where
+        end if
+
+     end if
+
   end if
   
   call excision_findboundary (ierr, emask, cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  ierr = Boundary_SelectGroupForBC (cctkGH, faces, ione, -ione, "spacemask::mask", "None")
+  if (ierr/=0) call CCTK_WARN (0, "internal error")
   
 end subroutine Lego_FixedSphere
diff --git a/src/maskvalues.h b/src/maskvalues.h
index bbaa95b..af1097d 100644
--- a/src/maskvalues.h
+++ b/src/maskvalues.h
@@ -2,6 +2,9 @@
 
 /* Values of the mask */
 
+#ifndef MASKVALUES_H
+#define MASKVALUES_H
+
 #ifdef CCODE
 
 #define MASK_EXCISED  0.0
@@ -17,3 +20,5 @@
 #define MASK_ACTIVE   1.0D0
 
 #endif
+
+#endif /* #ifdef MASKVALUES_H */