Robin BC in C - not tested yet, not in make file

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@61 6a38eb6e-646e-4a02-a296-d141613ad6c4

1 files changed, 128 insertions, 0 deletions

diff --git a/src/RobinBoundary.c b/src/RobinBoundary.c
new file mode 100644
index 0000000..b7b84ad
--- /dev/null
+++ b/src/RobinBoundary.c
@@ -0,0 +1,128 @@
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "math.h"
+
+void RobinBCi(cGH *GH, int *sw, int *doBC, int *lsh, 
+	      CCTK_REAL *var, CCTK_REAL *r, CCTK_REAL finf, CCTK_INT npow) {
+
+  int xgp1,xgp2,xgp3,xgp4;
+  int ygp1,ygp2,ygp3,ygp4;
+  int zgp1,zgp2,zgp3,zgp4;
+
+  int i,j,k;
+
+  if ((sw[0]!=1) || (sw[1]!=1) || (sw[2]!=1))
+    CCTK_WARN(0,"Robin BC only with stencil_size=1");
+
+  /* Lower x-bound: x(2,:,:) --> x[xgp2] */
+  if (doBC[0]==1) {
+    for (k=0;k<lsh[2];k++) {
+      for (j=0;j<lsh[1];j++) {
+	
+        xgp1 = CCTK_GFINDEX3D(GH,0,j,k);
+	xgp2 = CCTK_GFINDEX3D(GH,1,j,k);
+	xgp3 = CCTK_GFINDEX3D(GH,2,j,k);
+	xgp4 = CCTK_GFINDEX3D(GH,3,j,k);
+
+	var[xgp1] = finf + (
+                pow(r[xgp4],npow)*( var[xgp4] - finf )
+                -3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
+                +3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
+                )/pow(r[xgp1],npow);
+      }
+    }
+  }
+
+  /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
+  if (doBC[1]==1){
+    for (k=0;k<lsh[2];k++) {
+      for (j=0;j<lsh[1];j++) {	
+	xgp1 = CCTK_GFINDEX3D(GH,lsh[0]-1,j,k);
+        xgp2 = CCTK_GFINDEX3D(GH,lsh[0]-2,j,k);
+        xgp3 = CCTK_GFINDEX3D(GH,lsh[0]-3,j,k);
+	xgp4 = CCTK_GFINDEX3D(GH,lsh[0]-4,j,k);
+	
+	var[xgp1] = finf + (
+             pow(r[xgp3],npow)*( var[xgp4] - finf )
+             -3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
+             +3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
+             )/(r[xgp1],npow);
+      }
+    }
+  }
+  
+  /* Lower y-bound */
+  if (doBC[2] == 1) {
+    for (k=0;k<lsh[2];k++) {
+      for (i=0;i<lsh[0];i++) {
+	ygp1 = CCTK_GFINDEX3D(GH,i,0,k);
+	ygp2 = CCTK_GFINDEX3D(GH,i,1,k);
+	ygp3 = CCTK_GFINDEX3D(GH,i,2,k);
+	ygp4 = CCTK_GFINDEX3D(GH,i,3,k);
+	
+	var[ygp1] = finf + (
+                pow(r[ygp4],npow)*( var[ygp4] - finf )
+                -3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
+                +3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
+                )/pow(r[ygp1],npow);
+
+      }
+    }
+  }
+
+  /* Upper y bound */
+  if (doBC[3] == 1) {
+    for (k=0;k<lsh[2];k++) {
+      for (i=0;i<lsh[0];i++) {
+	ygp1 = CCTK_GFINDEX3D(GH,i,lsh[1]-1,k);
+	ygp2 = CCTK_GFINDEX3D(GH,i,lsh[1]-2,k);
+ 	ygp3 = CCTK_GFINDEX3D(GH,i,lsh[1]-3,k);
+	ygp4 = CCTK_GFINDEX3D(GH,i,lsh[1]-4,k);
+
+   	var[ygp1] = finf + (
+             pow(r[ygp3],npow)*( var[ygp4] - finf )
+             -3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
+             +3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
+             )/(r[ygp1],npow);
+	  }
+    }
+  }
+
+  /* Lower z-bound */
+  if (doBC[4]==1) {
+    for (j=0;j<lsh[1];j++) {
+      for (i=0;i<lsh[0];i++) {
+	zgp1 = CCTK_GFINDEX3D(GH,i,j,1);
+        zgp2 = CCTK_GFINDEX3D(GH,i,j,2);
+        zgp3 = CCTK_GFINDEX3D(GH,i,j,3);
+        zgp4 = CCTK_GFINDEX3D(GH,i,j,4);
+    
+	var[zgp1] = finf + (
+                pow(r[zgp4],npow)*( var[zgp4] - finf )
+                -3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
+                +3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
+                )/pow(r[zgp1],npow);
+      }
+    }
+  }
+
+  /* Upper z-bound */
+  if (doBC[5] == 1) {
+    for (j=0;j<lsh[1];j++) {
+      for (i=0;i<lsh[0];i++) {
+	zgp1 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-1);
+	zgp2 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-2);
+	zgp3 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-3);
+	zgp4 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-4);
+
+     	var[zgp1] = finf + (
+             pow(r[zgp3],npow)*( var[zgp4] - finf )
+             -3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
+             +3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
+             )/(r[zgp1],npow);
+
+      }
+    }
+  }
+}