Aliased routine for returning the diagonal norm coefficients.

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

3 files changed, 100 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
index 007620b..0c85802 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -26,6 +26,13 @@ SUBROUTINE Diff_coeff ( CCTK_POINTER_TO_CONST IN cctkGH, \
                         CCTK_INT IN table_handle )
 PROVIDES FUNCTION Diff_coeff WITH DiffCoeff LANGUAGE C
 
+#Routine to return the coefficients for the diagonal norms.
+SUBROUTINE GetScalProdDiag ( CCTK_POINTER_TO_CONST IN cctkGH, \
+                             CCTK_INT IN dir, \
+                             CCTK_INT IN nsize, \
+                             CCTK_REAL OUT ARRAY sigmad )
+PROVIDES FUNCTION GetScalProdDiag WITH SBP_GetScalProdDiag LANGUAGE Fortran
+
 CCTK_REAL FUNCTION GetScalProdCoeff ()
 PROVIDES FUNCTION GetScalProdCoeff WITH GetCoeff LANGUAGE Fortran
 
diff --git a/src/GetScalProdDiag.F90 b/src/GetScalProdDiag.F90
new file mode 100644
index 0000000..e2a7e99
--- /dev/null
+++ b/src/GetScalProdDiag.F90
@@ -0,0 +1,92 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+      
+subroutine SBP_GetScalProdDiag ( cctkGH, dir, nsize, sigmad )
+
+  implicit none
+
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_POINTER_TO_CONST, intent(IN) :: cctkGH
+  CCTK_INT, intent(IN) :: dir
+  CCTK_INT, intent(IN) :: nsize
+  CCTK_REAL, dimension(nsize), intent(OUT) :: sigmad
+
+  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
+  integer :: status
+  integer, dimension(6) :: bbox
+
+  CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 0.5_wp, 1.0_wp /) 
+  CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 17.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 = (/ 13649.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 = (/ 1498139.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/)
+                                                     
+  if ( dir < 0 .or. dir > 2 ) then
+    call CCTK_WARN(0, 'Error: dir is outside the legal range')
+  end if
+
+  call CCTK_GroupbboxGN ( status, cctkGH, 6, bbox, &
+                                  'SummationByParts::normmask' )
+  if ( status < 0 ) then
+    call CCTK_WARN(0,'Error: unable to get bounding box information')
+  end if
+
+  symtable = SymmetryTableHandleForGrid ( cctkGH )
+  call Util_TableGetIntArray ( n_elements, symtable, 6, &
+                                           symbnd, 'symmetry_handle' )
+  
+  pen_sym_handle = SymmetryHandleOfName ( 'multipatch' )
+
+  sigmad = 1.0_wp
+
+  if ( symbnd(dir+1) == pen_sym_handle .and. bbox(dir*2+1) == 1 ) then
+    select case (order)
+    case (2)
+      sigmad(1:2) = bmask_2
+    case (4)
+      sigmad(1:4) = bmask_4
+    case (6)
+      sigmad(1:6) = bmask_6
+    case (8)
+      sigmad(1:8) = bmask_8
+    end select
+  end if
+  if ( symbnd(dir+2) == pen_sym_handle .and. bbox(dir*2+2) == 1 ) then
+    select case (order)
+    case (2)
+      sigmad(nsize:nsize-1:-1) = bmask_2
+    case (4)
+      sigmad(nsize:nsize-3:-1) = bmask_4
+    case (6)
+      sigmad(nsize:nsize-5:-1) = bmask_6
+    case (8)
+      sigmad(nsize:nsize-7:-1) = bmask_8
+    end select
+  end if
+
+end subroutine SBP_GetScalProdDiag
diff --git a/src/make.code.defn b/src/make.code.defn
index 44c5741..1519773 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -16,6 +16,7 @@ SRCS =  call_derivs.c \
         Derivatives_6_5_min_err_coeff.F90 \
         Get_Coeff.F90 \
         set_norm_mask.F90 \
+        GetScalProdDiag.F90 \
         CheckGridSizes.F90 \
         dissipation.c \
         Dissipation_2_1.F90 \