(Changes from Ian Hawke.)

Possibly use ENO interpolation. Useful e.g. for hydro where shocks may appear.

To make Cartoon use the ENO interpolator you need an additional tags table entry. The same tags are used as for Carpet. Either

tags='Prolongation="TVD"'

or

tags='Prolongation="ENO"'

will work.

The default (i.e., with no tag) is to use the standard Lagrange polynomials. The tags "Lagrange" and "None" will also do this. Also, any direct call to the Cartoon functions (BndCartoon2DVI etc.) will use Lagrange interpolation.

Code from Burkhard Zink for interpolation; tags table stuff from me.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@76 eec4d7dc-71c2-46d6-addf-10296150bf52

1 files changed, 55 insertions, 0 deletions

diff --git a/src/interpolate.c b/src/interpolate.c
index 8a39c44..d48da6e 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -16,6 +16,7 @@
 */
 
 #include "cctk.h"
+#include "math.h"
 
 #define PRIVATE 
 #define PUBLIC 
@@ -40,6 +41,9 @@ PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
 					CCTK_REAL y[],
 					CCTK_REAL x);
 
+PRIVATE CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+                                  CCTK_REAL y[],
+                                  CCTK_REAL x);
 
 /*****************************************************************************/
 
@@ -225,3 +229,54 @@ CCTK_REAL xr = (x - xc) / dx;
 return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
 }
 
+
+/*****************************************************************************/
+
+CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
+{
+
+  CCTK_REAL diff[4];
+  CCTK_INT seed = 0;
+  CCTK_INT j,r;  
+
+  /* Find seed index */
+  while (x > x0 + dx * ((float)seed+0.5)) seed++;
+
+  if (seed!=2) {
+    /* Not enough stencil, only perform linear interpolation */   
+    seed = 0;
+    while (x > x0 + dx * ((float)(seed+1))) seed++;
+    if (seed==4) seed=3;
+    return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(float)seed*dx);
+  }
+
+  /*  Calculate the undivided differences */
+  for (j=0; j<=3; j++)
+    diff[j] = y[j+1] - y[j];
+
+  /*  Find the stencil */
+  if ( fabs(diff[1]) < fabs(diff[2]) ) {
+    if ( fabs(diff[1]-diff[0]) < fabs(diff[2]-diff[1]) )
+      r = 0;
+    else
+      r = 1;
+  } else {
+    if ( fabs(diff[2]-diff[1]) < fabs(diff[3]-diff[2]) )
+      r = 1;
+    else
+      r = 2;
+  }
+
+  /* Interpolate second order */
+  double c0 = y[r+1];
+  double c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2];
+  double c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2];
+
+  double xc = x0 + dx * (double)(r+1);
+  double xr = (x - xc) / dx;
+
+  return(  c0 + xr*(c1 + xr*c2)  );  
+
+}