Added 4-3 Full restricted norm operators. Added parameter to choose the type

of differences (diagonal, full restricted). Added routine to return the
coefficient of the norm at the boundary. Added routine to return the
mask for the norm calculation.


git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@10 f69c4107-0314-4c4f-9ad4-17e986b73f4a

7 files changed, 538 insertions, 76 deletions

diff --git a/src/Derivatives_4_3.F90 b/src/Derivatives_4_3.F90
new file mode 100644
index 0000000..011a074
--- /dev/null
+++ b/src/Derivatives_4_3.F90
@@ -0,0 +1,242 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_4_3 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_INT, intent(IN) :: ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
+  CCTK_INT, intent(IN) :: dir
+  CCTK_INT, intent(IN) :: bb(2)
+  CCTK_INT, intent(IN) :: gsize
+  CCTK_REAL, intent(IN) :: delta
+  CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar
+
+  CCTK_REAL, dimension(2), save :: a
+  CCTK_REAL, dimension(7,5), save :: q 
+  CCTK_REAL :: idel
+  CCTK_REAL :: f0, f1, f2, f3
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+  logical, save :: first = .true.
+
+  if ( first ) then
+    a(1) = 2.0_wp/3.0_wp; a(2) = -1.0_wp/12.0_wp
+
+    f0 = sqrt(26116897.0_wp)
+    f1 = -56764003702447356523.0_wp + 8154993476273221.0_wp * f0
+    f2 = -55804550303.0_wp + 9650225.0_wp * f0
+    f3 = 3262210757.0_wp + 271861.0_wp * f0
+
+    q(1,1) = -11._wp/6.0_wp; q(2,1) = 3.0_wp
+    q(3,1) = -3.0_wp/2.0_wp; q(4,1) = 1.0_wp/3.0_wp
+    q(5,1) = zero; q(6,1) = zero; q(7,1) = zero
+
+    q(1,2) =  -24.0_wp * ( -779042810827742869.0_wp + &
+                            104535124033147.0_wp * f0 ) / f1
+    q(2,2) = -( -176530817412806109689.0_wp + &
+                 29768274816875927.0_wp * f0 ) / ( 6.0_wp * f1 )
+    q(3,2) = 343.0_wp * ( -171079116122226871.0_wp + &
+                           27975630462649.0_wp * f0 ) / f1
+    q(4,2) = -3.0_wp * ( -7475554291248533227.0_wp + &
+                          1648464218793925.0_wp * f0 ) / ( 2.0_wp * f1 )
+    q(5,2) = ( -2383792768180030915.0_wp + &
+                1179620587812973.0_wp * f0 ) / ( 3.0_wp * f1 )
+    q(6,2) = -1232.0_wp * ( -115724529581315.0_wp + 37280576429.0_wp * f0 ) / f1
+    q(7,2) = zero
+
+    q(1,3) = -12.0_wp * ( -380966843.0_wp + 86315.0_wp * f0 ) / f2
+    q(2,3) = ( 5024933015.0_wp + 2010631.0_wp * f0 ) / ( 3.0_wp * f2 )
+    q(3,3) = -231.0_wp * ( -431968921.0_wp + 86711.0_wp * f0 ) / ( 2.0_wp * f2 )
+    q(4,3) = ( -65931742559.0_wp + 12256337.0_wp * f0 ) / f2
+    q(5,3) = -( -50597298167.0_wp + 9716873.0_wp * f0 ) / ( 6.0_wp * f2 )
+    q(6,3) = -88.0_wp * ( -15453061.0_wp + 2911.0_wp * f0 ) / f2
+    q(7,3) = zero
+
+    q(1,4) = 48.0_wp * ( -56020909845192541.0_wp + &
+                          9790180507043.0_wp * f0 ) / f1
+    q(2,4) = ( -9918249049237586011.0_wp + &
+                1463702013196501.0_wp * f0 ) / ( 6.0_wp * f1)
+    q(3,4) = -13.0_wp * ( -4130451756851441723.0_wp + &
+                           664278707201077.0_wp * f0 ) / f1
+    q(4,4) = 3.0_wp * ( -26937108467782666617.0_wp + &
+                         5169063172799767.0_wp * f0 ) / ( 2.0_wp * f1 )
+    q(5,4) = -( 6548308508012371315.0_wp + &
+                3968886380989379.0_wp * f0 ) / ( 3.0_wp * f1 )
+    q(6,4) = 88.0_wp * ( -91337851897923397.0_wp + &
+                          19696768305507.0_wp * f0 ) / f1
+    q(7,4) = 242.0_wp * ( -120683.0_wp + 15.0_wp * f0 ) / f3
+
+    q(1,5) = 264.0_wp * ( -120683.0_wp + 15.0_wp * f0 ) / f3
+    q(2,5) = ( -43118111.0_wp + 23357.0_wp * f0 ) / ( 3.0_wp * f3 )
+    q(3,5) = -47.0_wp * ( -28770085.0_wp + 2259.0_wp * f0 ) / ( 2.0_wp * f3 ) 
+    q(4,5) = -3.0_wp * ( 1003619433.0_wp + 11777.0_wp * f0 ) / f3
+    q(5,5) = -11.0_wp * ( -384168269.0_wp + 65747.0_wp * f0 ) / ( 6.0_wp * f3 )
+    q(6,5) = 22.0_wp * ( 87290207.0_wp + 10221.0_wp * f0 ) / f3
+    q(7,5) = -66.0_wp * ( 3692405.0_wp + 419.0_wp * f0 ) / f3
+    first = .false.
+  end if
+
+  idel = 1.0_wp / delta
+
+  direction: select case (dir)
+  case (0) direction
+    if ( bb(1) == 0 ) then
+      il = 1 + gsize
+    else
+      dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
+                      q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel
+      dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
+                      q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + &
+                      q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) ) * idel
+      dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+                      q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + &
+                      q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) ) * idel
+      dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + &
+                      q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + &
+                      q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + &
+                      q(7,4) * var(7,:,:) ) * idel
+      dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + &
+                      q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + &
+                      q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + &
+                      q(7,5) * var(7,:,:) ) * idel
+      il = 6
+    end if
+    if ( bb(2) == 0 ) then
+      ir = ni - gsize
+    else
+      dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + &
+                           q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + &
+                           q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + &
+                           q(7,5) * var(ni-6,:,:) ) * idel
+      dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + &
+                           q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + &
+                           q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + &
+                           q(7,4) * var(ni-6,:,:) ) * idel
+      dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + &
+                           q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + &
+                           q(5,3) * var(ni-4,:,:) + &
+                           q(6,3) * var(ni-5,:,:) ) * idel
+      dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + &
+                           q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + &
+                           q(5,2) * var(ni-4,:,:) + &
+                           q(6,2) * var(ni-5,:,:) ) * idel
+      dvar(ni,:,:)   = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + &
+                           q(3,1) * var(ni-2,:,:) + &
+                           q(4,1) * var(ni-3,:,:) ) * idel
+      ir = ni - 5
+    end if
+    dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - &
+                                 var(il-1:ir-1,:,:) ) + &
+                        a(2) * ( var(il+2:ir+2,:,:) - &
+                                 var(il-2:ir-2,:,:) ) ) * idel
+  case (1) direction
+    if ( bb(1) == 0 ) then
+      jl = 1 + gsize
+    else
+      dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+                      q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel
+      dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
+                      q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + &
+                      q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) ) * idel
+      dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+                      q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + &
+                      q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) ) * idel
+      dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + &
+                      q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + &
+                      q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + &
+                      q(7,4) * var(:,7,:) ) * idel
+      dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + &
+                      q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + &
+                      q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + &
+                      q(7,5) * var(:,7,:) ) * idel
+      jl = 6
+    end if
+    if ( bb(2) == 0 ) then
+      jr = nj - gsize
+    else
+      dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + &
+                           q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + &
+                           q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + &
+                           q(7,5) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + &
+                           q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + &
+                           q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + &
+                           q(7,4) * var(:,nj-6,:) ) * idel
+      dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + &
+                           q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + &
+                           q(5,3) * var(:,nj-4,:) + &
+                           q(6,3) * var(:,nj-5,:) ) * idel
+      dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + &
+                           q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + &
+                           q(5,2) * var(:,nj-4,:) + &
+                           q(6,2) * var(:,nj-5,:) ) * idel
+      dvar(:,nj,:)   = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + &
+                           q(3,1) * var(:,nj-2,:) + &
+                           q(4,1) * var(:,nj-3,:) ) * idel
+      jr = nj - 5
+    end if
+    dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - &
+                                 var(:,jl-1:jr-1,:) ) + &
+                        a(2) * ( var(:,jl+2:jr+2,:) - &
+                                 var(:,jl-2:jr-2,:) ) ) * idel
+  case (2) direction
+    if ( bb(1) == 0 ) then
+      kl = 1 + gsize
+    else
+      dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+                      q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel
+      dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
+                      q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + &
+                      q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) ) * idel
+      dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+                      q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + &
+                      q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) ) * idel
+      dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + &
+                      q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + &
+                      q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + &
+                      q(7,4) * var(:,:,7) ) * idel
+      dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + &
+                      q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + &
+                      q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + &
+                      q(7,5) * var(:,:,7) ) * idel
+      kl = 6
+    end if
+    if ( bb(2) == 0 ) then
+      kr = nk - gsize
+    else 
+      dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + &
+                           q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + &
+                           q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + &
+                           q(7,5) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + &
+                           q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + &
+                           q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + &
+                           q(7,4) * var(:,:,nk-6) ) * idel
+      dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + &
+                           q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + &
+                           q(5,3) * var(:,:,nk-4) + &
+                           q(6,3) * var(:,:,nk-5) ) * idel
+      dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + &
+                           q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + &
+                           q(5,2) * var(:,:,nk-4) + &
+                           q(6,2) * var(:,:,nk-5) ) * idel
+      dvar(:,:,nk)   = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + &
+                           q(3,1) * var(:,:,nk-2) + &
+                           q(4,1) * var(:,:,nk-3) ) * idel
+      kr = nk - 5
+    end if
+    dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - &
+                                 var(:,:,kl-1:kr-1) ) + &
+                        a(2) * ( var(:,:,kl+2:kr+2) - &
+                                 var(:,:,kl-2:kr-2) ) ) * idel
+  end select direction
+end subroutine deriv_gf_4_3
diff --git a/src/Get_Coeff.F90 b/src/Get_Coeff.F90
new file mode 100644
index 0000000..7b3fcf5
--- /dev/null
+++ b/src/Get_Coeff.F90
@@ -0,0 +1,32 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+CCTK_REAL function GetCoeff ()
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+
+  if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
+    select case (order)
+    case (2)
+      GetCoeff = 1.0_wp / 2.0_wp
+    case (4)
+      GetCoeff = 17.0_wp/48.0_wp
+    case (6)
+      GetCoeff = 13649.0_wp/43200.0_wp
+    case (8)
+      GetCoeff =  1498139.0_wp/5080320.0_wp
+    end select
+  else
+    select case (order)
+    case(4)
+      GetCoeff = 3.0_wp/11.0_wp
+    end select
+  end if
+end function GetCoeff
diff --git a/src/call_derivs.c b/src/call_derivs.c
index e8c051a..f8807c2 100644
--- a/src/call_derivs.c
+++ b/src/call_derivs.c
@@ -33,6 +33,15 @@ void DiffGv ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
                                        const CCTK_INT *gsize,
                                        const CCTK_REAL *delta,
                                        CCTK_REAL *dvar);
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_4_3)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
+                                       const CCTK_INT *dir,
+                                       const CCTK_INT *bb,
+                                       const CCTK_INT *gsize,
+                                       const CCTK_REAL *delta,
+                                       CCTK_REAL *dvar);
   void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(const CCTK_REAL *var,
                                        const CCTK_INT *ni,
                                        const CCTK_INT *nj,
@@ -78,24 +87,35 @@ void DiffGv ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
       assert(0);
   }
 
-  switch(order) {
-  case 2: {
-    CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  case 4: {
-    CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  case 6: {
-    CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  case 8: {
-    CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  default:
-    assert(0);
+  if ( CCTK_Equals(norm_type,"Diagonal") ) {
+    switch(order) {
+    case 2: {
+      CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 4: {
+      CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 6: {
+      CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 8: {
+      CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    default:
+      assert(0);
+    }
+  } else {
+    switch(order) {
+    case 4: {
+      CCTK_FNAME(deriv_gf_4_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    default:
+      assert(0);
+    }
   }
 } 
diff --git a/src/call_derivs_name.c b/src/call_derivs_name.c
index c696e4e..893aa43 100644
--- a/src/call_derivs_name.c
+++ b/src/call_derivs_name.c
@@ -1,82 +1,87 @@
 #include "cctk.h"
 #include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
 
 #include <assert.h>
 
-void DiffGf ( const CCTK_POINTER cctkGH, const CCTK_INT dir,
+void DiffGf ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
               const char *var_name, const char *dvar_name )
 {
+  cGH const * restrict const cctkGH = cctkGH_;
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+
   CCTK_REAL *var, *dvar;
   CCTK_INT ni, nj, nk, gsize, ic;
   CCTK_REAL delta;
-  CCTK_REAL phys_min[3], phys_max[3];
-  CCTK_REAL int_min[3], int_max[3];
-  CCTK_REAL ext_min[3], ext_max[3];
-  CCTK_REAL spacing[3];
   CCTK_INT ierr;
   CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3];
-  void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(CCTK_REAL *var, const CCTK_INT *ni,
-                                       const CCTK_INT *nj, const CCTK_INT *nk,
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
                                        const CCTK_INT *dir,
                                        const CCTK_INT *bb,
                                        const CCTK_INT *gsize,
                                        const CCTK_REAL *delta,
                                        CCTK_REAL *dvar);
-  void CCTK_FCALL CCTK_FNAME(deriv_gf_4_2)(CCTK_REAL *var, const CCTK_INT *ni,
-                                       const CCTK_INT *nj, const CCTK_INT *nk,
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_4_2)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
                                        const CCTK_INT *dir,
                                        const CCTK_INT *bb,
                                        const CCTK_INT *gsize,
                                        const CCTK_REAL *delta,
                                        CCTK_REAL *dvar);
-  void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(CCTK_REAL *var, const CCTK_INT *ni,
-                                       const CCTK_INT *nj, const CCTK_INT *nk,
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_4_3)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
+                                       const CCTK_INT *dir,
+                                       const CCTK_INT *bb,
+                                       const CCTK_INT *gsize,
+                                       const CCTK_REAL *delta,
+                                       CCTK_REAL *dvar);
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
+                                       const CCTK_INT *dir,
+                                       const CCTK_INT *bb,
+                                       const CCTK_INT *gsize,
+                                       const CCTK_REAL *delta,
+                                       CCTK_REAL *dvar);
+  void CCTK_FCALL CCTK_FNAME(deriv_gf_8_4)(const CCTK_REAL *var,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
                                        const CCTK_INT *dir,
                                        const CCTK_INT *bb,
                                        const CCTK_INT *gsize,
                                        const CCTK_REAL *delta,
                                        CCTK_REAL *dvar);
-  DECLARE_CCTK_PARAMETERS
 
   
-  if (CCTK_GrouplshVN(cctkGH, 3, lsh, var_name) < 0)
-  {
-    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-               "Error getting local size of grid for variable '%s'",var_name);
-  }
-  if (CCTK_GroupbboxVN(cctkGH, 6, bbox, var_name) < 0)
-  {
-    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-               "Error getting bounding box for variable '%s'",var_name);
-  }
-  if (CCTK_GroupnghostzonesVN(cctkGH, 3, nghostzones, var_name) < 0)
-  {
-    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-               "Error getting number of ghostzones for variable '%s'",var_name);
-  }
-
-  ni = lsh[0]; nj = lsh[1]; nk = lsh[2];
-
-  ierr = GetDomainSpecification ( 3, phys_min, phys_max, int_min, int_max,
-                                                  ext_min, ext_max, spacing );
+  ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
   
   switch(dir) {
   case 0: {
-    delta = spacing[0];
-    bb[0] = bbox[0]; bb[1] = bbox[1];
-    gsize = nghostzones[0];
+    delta = CCTK_DELTA_SPACE(0);
+    bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1];
+    gsize = cctk_nghostzones[0];
     break;
   }
   case 1: {
-    delta = spacing[1];
-    bb[0] = bbox[2]; bb[1] = bbox[3];
-    gsize = nghostzones[1];
+    delta = CCTK_DELTA_SPACE(1);
+    bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3];
+    gsize = cctk_nghostzones[1];
     break;
   }
   case 2: {
-    delta = spacing[2];
-    bb[0] = bbox[4]; bb[1] = bbox[5];
-    gsize = nghostzones[2];
+    delta = CCTK_DELTA_SPACE(2);
+    bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5];
+    gsize = cctk_nghostzones[2];
     break;
   }
   default:
@@ -86,20 +91,35 @@ void DiffGf ( const CCTK_POINTER cctkGH, const CCTK_INT dir,
   var = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,var_name));
   dvar = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,dvar_name));
 
-  switch(order) {
-  case 2: {
-    CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  case 4: {
-    CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  case 6: {
-    CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
-    break;
-  }
-  default:
-    assert(0);
+  if ( CCTK_Equals(norm_type,"Diagonal") ) {
+    switch(order) {
+    case 2: {
+      CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 4: {
+      CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 6: {
+      CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    case 8: {
+      CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    default:
+      assert(0);
+    }
+  } else {
+    switch(order) {
+    case 4: {
+      CCTK_FNAME(deriv_gf_4_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+      break;
+    }
+    default:
+      assert(0);
+    }
   }
 } 
diff --git a/src/call_norm_mask.c b/src/call_norm_mask.c
new file mode 100644
index 0000000..5e7996a
--- /dev/null
+++ b/src/call_norm_mask.c
@@ -0,0 +1,24 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include <assert.h>
+
+void SetNMask ( const CCTK_POINTER cctkGH_, CCTK_REAL *mask )
+{
+  cGH const * restrict const cctkGH = cctkGH_;
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+
+  CCTK_INT ni, nj, nk;
+  void CCTK_FCALL CCTK_FNAME(sbp_set_norm_mask)(const CCTK_POINTER cctkGH,
+                                       const CCTK_INT *ni,
+                                       const CCTK_INT *nj,
+                                       const CCTK_INT *nk,
+                                       CCTK_REAL *mask);
+
+  
+  ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
+
+  CCTK_FNAME(sbp_set_norm_mask)(cctkGH_,&ni,&nj,&nk,mask);
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index a5db2d8..893ee3f 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS =  Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c call_derivs_name.c
+SRCS =  Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c Derivatives_4_3.F90 call_derivs_name.c Get_Coeff.F90 call_norm_mask.c set_norm_mask.F90
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/set_norm_mask.F90 b/src/set_norm_mask.F90
new file mode 100644
index 0000000..30b3c63
--- /dev/null
+++ b/src/set_norm_mask.F90
@@ -0,0 +1,124 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+      
+subroutine sbp_set_norm_mask (cctkGH, ni, nj, nk, mask)
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+      
+  CCTK_POINTER, intent(IN) :: cctkGH
+  CCTK_INT, intent(IN) :: ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: mask
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np
+  CCTK_INT, dimension(6) :: symbnd
+  CCTK_POINTER :: psym_name
+  character(len=256) :: symmetry_name
+  CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z
+
+! Note: The first number is twice the value from the paper, since Carpet
+! Multiplies by 1/2 at the boundary when doing the sum reduction.
+  CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 1.0_wp, 1.0_wp /) 
+  CCTK_REAL, dimension(3), parameter :: bmask_3 = (/ 1.0_wp, 1.0_wp, 1.0_wp /) 
+  CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 34.0_wp/48.0_wp, &
+                                                     59.0_wp/48.0_wp, &
+                                                     43.0_wp/48.0_wp, &
+                                                     49.0_wp/48.0_wp /)
+  CCTK_REAL, dimension(6), parameter :: bmask_6 = (/ 27298.0_wp/43200._wp, &
+                                                     12013.0_wp/8640._wp, &
+                                                     2711.0_wp/4320.0_wp, &
+                                                     5359.0_wp/4320.0_wp, &
+                                                     7877.0_wp/8640.0_wp, &
+                                                     43801.0_wp/43200.0_wp /)
+  CCTK_REAL, dimension(8), parameter :: bmask_8 = (/ 2996278.0_wp/5080320.0_wp,&
+                                                     1107307.0_wp/725760.0_wp, &
+                                                     20761.0_wp/80640.0_wp, &
+                                                     1304999.0_wp/725760.0_wp, &
+                                                     299527.0_wp/725760.0_wp, &
+                                                     103097.0_wp/80640.0_wp, &
+                                                     670091.0_wp/725760.0_wp, &
+                                                     5127739.0_wp/5080320.0_wp/)
+                                                     
+  CCTK_REAL, dimension(8) :: bmask
+
+  symtable = SymmetryTableHandleForGrid ( cctkGH )
+  call Util_TableGetIntArray ( n_elements, symtable, 6, &
+                                           symbnd, 'symmetry_handle' )
+  
+  pen_sym_handle = SymmetryHandleOfName ( 'penalty inter-patch' )
+
+  if ( any ( symbnd == pen_sym_handle ) ) then
+    allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
+    mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0
+
+    if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
+      select case (order)
+      case (2)
+        bmask(1:2) = bmask_2
+      case (4)
+        bmask(1:4) = bmask_4
+      case(6)
+        bmask(1:6) = bmask_6
+      case(8)
+        bmask(1:8) = bmask_8
+      end select 
+      np = order
+    else
+      select case (order)
+      case (4)
+        bmask(1:3) = bmask_3
+        np = 3
+      end select
+    end if
+
+    if ( symbnd(1) == pen_sym_handle ) then
+        mask_x(1:np) = bmask(1:np)
+    end if
+    if ( symbnd(2) == pen_sym_handle ) then
+      mask_x(ni:ni-np+1:-1) = bmask(1:np)
+    end if
+    if ( symbnd(3) == pen_sym_handle ) then
+      mask_y(1:np) = bmask(1:np)
+    end if
+    if ( symbnd(4) == pen_sym_handle ) then
+      mask_y(nj:nj-np+1:-1) = bmask(1:np)
+    end if
+    if ( symbnd(5) == pen_sym_handle ) then
+      mask_z(1:np) = bmask(1:np)
+    end if
+    if ( symbnd(6) == pen_sym_handle ) then
+      mask_z(nk:nk-np+1:-1) = bmask(1:np)
+    end if
+
+    mask = outer_prod_3d ( mask_x, mask_y, mask_z )
+  else
+    mask = 1.0_wp
+  end if
+
+contains
+
+  function outer_prod_3d ( ax, ay, az )
+    CCTK_REAL, dimension(:), intent(IN) :: ax
+    CCTK_REAL, dimension(:), intent(IN) :: ay
+    CCTK_REAL, dimension(:), intent(IN) :: az
+    CCTK_REAL, dimension(size(ax),size(ay),size(az)) :: outer_prod_3d
+    CCTK_INT :: nx, ny, nz, i, j, k
+  
+    nx = size(ax); ny = size(ay); nz = size(az)
+    do k = 1, nz
+      do j = 1, ny
+        do i = 1, nx
+          outer_prod_3d(i,j,k) = ax(i) * ay(j) * az(k)
+        end do
+      end do
+    end do
+  end function outer_prod_3d
+
+end subroutine sbp_set_norm_mask