Changes in methodology for filling in the interior of Cook-Pfeiffer data as

used in the "turducken" paper. This includes some spline blending from
Erik and the most recent method of filling in the interior using an elliptic
equation with the "good" points as boundary values. This method uses a
conjugate gradient solver (built in), that should be robust (i.e. it should
always work) but may not be the fastest. However, since this is only done on
initial data this shouldn't be an issue. The elliptic method currently support
a second order, fourth order and sixth order equation, that will give C0, C1
and C2 solutions, respectively, across the boundary of the fill in region.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@7 4ec1db94-0e4f-0410-ada3-8bed251432c9

1 files changed, 127 insertions, 79 deletions

diff --git a/src/overwrite.F90 b/src/overwrite.F90
index 6c6089f..3df9ffb 100644
--- a/src/overwrite.F90
+++ b/src/overwrite.F90
@@ -10,85 +10,114 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
   DECLARE_CCTK_PARAMETERS
-  
+
   CCTK_REAL, parameter :: zero=0
   CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
-  CCTK_REAL :: cx, cy, cz, rad, width
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
   integer   :: n
+
+  integer, parameter :: sm_linear = 1
+  integer, parameter :: sm_spline = 2
+  integer, parameter :: sm_cosine = 3
+  integer :: sm_type
+
+  if (CCTK_EQUALS(method,"old")) then
+    do n = 1, num_regions
+
+       cx = centre_x(n)
+       cy = centre_y(n)
+       cz = centre_z(n)
+       if (CCTK_EQUALS(region_shape(n), "sphere")) then
+          radx = radius(n)
+          rady = radius(n)
+          radz = radius(n)
+       else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+          radx = radius_x(n)
+          rady = radius_y(n)
+          radz = radius_z(n)
+       else
+          call CCTK_WARN (0, "internal error")
+       end if
+
+       if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+          sm_type = sm_linear
+       else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+          sm_type = sm_spline
+       else if (CCTK_EQUALS (smoothing_function(n), "cosine")) then
+          sm_type = sm_cosine
+       else
+          call CCTK_WARN (0, "internal error")
+       end if
+       width = smoothing_zone_width(n)
+
+       dist2 =   ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+            &  + ((z - cz) / radz)**2
+
+       if (overwrite_geometry(n) /= 0) then
+
+          if (conformal_state >= 1) then
+             where (dist2 <= 1)
+                psi = 1
+             end where
+          end if
+          if (conformal_state >= 2) then
+             where (dist2 <= 1)
+                psix = 0
+                psiy = 0
+                psiz = 0
+             end where
+          end if
+          if (conformal_state >= 3) then
+             where (dist2 <= 1)
+                psixx = 0
+                psixy = 0
+                psixz = 0
+                psiyy = 0
+                psiyz = 0
+                psizz = 0
+             end where
+          end if
   
-  do n = 1, num_regions
-     
-     cx = centre_x(n)
-     cy = centre_y(n)
-     cz = centre_z(n)
-     rad = radius(n)
-     width = smoothing_zone_width(n)
-     
-     dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
-     
-     if (overwrite_geometry(n) /= 0) then
-        
-        if (conformal_state >= 1) then
-           where (dist2 <= rad**2)
-              psi = 1
-           end where
-        end if
-        if (conformal_state >= 2) then
-           where (dist2 <= rad**2)
-              psix = 0
-              psiy = 0
-              psiz = 0
-           end where
-        end if
-        if (conformal_state >= 3) then
-           where (dist2 <= rad**2)
-              psixx = 0
-              psixy = 0
-              psixz = 0
-              psiyy = 0
-              psiyz = 0
-              psizz = 0
-           end where
-        end if
-        
-        where (dist2 <= rad**2)
-           gxx = smooth (gxx,  Minkowski_scale(n), dist2)
-           gxy = smooth (gxy,  zero              , dist2)
-           gxz = smooth (gxz,  zero              , dist2)
-           gyy = smooth (gyy,  Minkowski_scale(n), dist2)
-           gyz = smooth (gyz,  zero              , dist2)
-           gzz = smooth (gzz,  Minkowski_scale(n), dist2)
-           kxx = smooth (kxx,  zero              , dist2)
-           kxy = smooth (kxy,  zero              , dist2)
-           kxz = smooth (kxz,  zero              , dist2)
-           kyy = smooth (kyy,  zero              , dist2)
-           kyz = smooth (kyz,  zero              , dist2)
-           kzz = smooth (kzz,  zero              , dist2)
-        end where
-        
-     end if
-     
-     if (overwrite_lapse(n) /= 0) then
-        
-        where (dist2 <= rad**2)
-           alp = smooth (alp, lapse_scale(n), dist2)
-        end where
-        
-     end if
-        
-     if (overwrite_shift(n) /= 0) then
-        
-        if (shift_state /= 0) then
-           where (dist2 <= rad**2)
-              betax = smooth (betax, zero, dist2)
-              betay = smooth (betay, zero, dist2)
-              betaz = smooth (betaz, zero, dist2)
-           end where
-        end if
-        
-     end if
-     
-  end do
+          where (dist2 <= 1)
+             gxx = smooth (gxx,  Minkowski_scale(n), dist2)
+             gxy = smooth (gxy,  zero              , dist2)
+             gxz = smooth (gxz,  zero              , dist2)
+             gyy = smooth (gyy,  Minkowski_scale(n), dist2)
+             gyz = smooth (gyz,  zero              , dist2)
+             gzz = smooth (gzz,  Minkowski_scale(n), dist2)
+             kxx = smooth (kxx,  zero              , dist2)
+             kxy = smooth (kxy,  zero              , dist2)
+             kxz = smooth (kxz,  zero              , dist2)
+             kyy = smooth (kyy,  zero              , dist2)
+             kyz = smooth (kyz,  zero              , dist2)
+             kzz = smooth (kzz,  zero              , dist2)
+          end where
+
+       end if
+
+       if (overwrite_lapse(n) /= 0) then
+
+          where (dist2 <= 1)
+             alp = smooth (alp, lapse_scale(n), dist2)
+          end where
+
+       end if
+
+       if (overwrite_shift(n) /= 0) then
+
+          if (shift_state /= 0) then
+             where (dist2 <= 1)
+                betax = smooth (betax, zero, dist2)
+                betay = smooth (betay, zero, dist2)
+                betaz = smooth (betaz, zero, dist2)
+             end where
+          end if
+
+       end if
+
+    end do
+    
+  end if
   
 contains
   
@@ -98,13 +127,32 @@ contains
     CCTK_REAL, intent(in) :: goodval
     CCTK_REAL, intent(in) :: dist2
     
-    if (dist2 <= (rad-width)**2) then
+    CCTK_REAL, parameter :: zero=0, two=2, half=1/two
+    integer,   parameter :: rk = kind(zero)
+    CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk
+    CCTK_REAL :: dist
+    CCTK_REAL :: x, f
+    
+    if (dist2 <= (1 - width)**2) then
        newval = goodval
-    else if (width == 0 .or. dist2 >= rad**2) then
+    else if (width == 0 .or. dist2 >= 1) then
        newval = oldval
     else
-       newval =   goodval &
-            &   + (oldval - goodval) * (sqrt(dist2) - (rad-width)) / width
+       dist = sqrt(dist2)
+       x = 1 - (1 - dist) / width
+       select case (sm_type)
+       case (sm_linear)
+          f = x
+       case (sm_spline)
+          if (x < half) then
+             f = 2 * x**2
+          else
+             f = 1 - 2 * (1-x)**2
+          end if
+       case (sm_cosine)
+          f = (1 - cos (pi * x)) / 2
+       end select
+       newval = oldval * f + goodval * (1 - f)
     end if
     
   end function smooth