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

9 files changed, 558 insertions, 76 deletions

diff --git a/interface.ccl b/interface.ccl
index 9eef8aa..c05a3c1 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -25,6 +25,13 @@ SUBROUTINE Diff_gv ( CCTK_POINTER IN cctkGH, \
                      CCTK_REAL OUT ARRAY dvar )
 PROVIDES FUNCTION Diff_gv WITH DiffGv LANGUAGE C
 
+CCTK_REAL FUNCTION GetScalProdCoeff ()
+PROVIDES FUNCTION GetScalProdCoeff WITH GetCoeff LANGUAGE Fortran
+
+SUBROUTINE SetNormMask ( CCTK_POINTER IN cctkGH, \
+                         CCTK_REAL OUT ARRAY mask )
+PROVIDES FUNCTION SetNormMask WITH SetNMask LANGUAGE C
+
 CCTK_INT FUNCTION GetDomainSpecification \
   (CCTK_INT IN size, \
    CCTK_REAL OUT ARRAY physical_min, \
@@ -35,3 +42,10 @@ CCTK_INT FUNCTION GetDomainSpecification \
    CCTK_REAL OUT ARRAY exterior_max, \
    CCTK_REAL OUT ARRAY spacing)
 USES FUNCTION GetDomainSpecification
+
+CCTK_INT FUNCTION \
+    SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
+REQUIRES FUNCTION SymmetryTableHandleForGrid
+
+CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name)
+REQUIRES FUNCTION SymmetryHandleOfName
diff --git a/param.ccl b/param.ccl
index 4881ebf..285f928 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,6 +1,12 @@
 # Parameter definitions for thorn SummationByParts
 # $Header$
 
+KEYWORD norm_type "Type of norm"
+{
+  "Diagonal" :: "Diagonal norm"
+  "Full restricted" :: "Full resricted norm"
+} "Diagonal"
+
 INT order "Order of accuracy" STEERABLE=always
 {
   2:8:2 :: ""
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