diff options
Diffstat (limited to 'src/Cartoon2DBC.c')
| -rw-r--r-- | src/Cartoon2DBC.c | 217 |
1 files changed, 212 insertions, 5 deletions
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c index f74c30f..7653424 100644 --- a/src/Cartoon2DBC.c +++ b/src/Cartoon2DBC.c @@ -24,6 +24,7 @@ #include "cctk_FortranString.h" #include "Cartoon2D_tensors.h" +#include "SpaceMask.h" static const char *rcsid = "$Id$"; @@ -38,9 +39,11 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, CCTK_REAL *f, int i, int di, int ijk, CCTK_REAL x); CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], CCTK_REAL x); + CCTK_REAL y[], CCTK_REAL x, + CCTK_INT excised_points[]); CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], CCTK_REAL x); + CCTK_REAL y[], CCTK_REAL x, + CCTK_INT excised_points[]); void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, int *tensortype, int *vi); void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH, @@ -49,6 +52,9 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN) (int *retval, const cGH **GH, int *tensortype, ONE_FORTSTRING_ARG); +void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS); +void SetupExcisionVars(CCTK_ARGUMENTS); + /* set boundaries of a grid tensor assuming axisymmetry - handles lower boundary in x - does not touch other boundaries @@ -90,7 +96,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi) CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz; CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz; - gi = CCTK_GroupIndexFromVarI (vi); assert (gi >= 0); n_vars = CCTK_NumVarsInGroupI (gi); @@ -342,6 +347,59 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, int lnx; int n, offset; + CCTK_REAL *old_mask; + CCTK_INT *new_mask; + + CCTK_INT excised_points[6]; + + if (*excision_active < 0) + { + SetupExcisionVars((cGH*)cctkGH); + } + + if (old_excision) + { + old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi); + } + else + { + old_mask = NULL; + } + + if (new_excision) + { + new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi); + } + else + { + new_mask = NULL; + } + + /* + If this boundary point is excised then we return + without doing anything + */ + + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + if (old_mask[ijk] < 0.75) + { + return f[ijk]; + } + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + if (SpaceMask_CheckStateBits(new_mask, (ijk), + *new_excision_field, + *new_excision_descriptor)) + { + return f[ijk]; + } + } + } + lnx = cctkGH->cctk_lsh[0]; dx0 = CCTK_DELTA_SPACE(0); @@ -369,10 +427,25 @@ static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, x0 = lx0 + dx0 * (i-offset); for (n = 0; n <= order; n++) { y[n] = f[ijk-offset+n]; + excised_points[n] = 0; + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1; + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + excised_points[n] |= + (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n), + *new_excision_field, + *new_excision_descriptor)); + } + } } /* call interpolator */ - r = interpolate_local(order, x0, dx0, y, x); + r = interpolate_local(order, x0, dx0, y, x, excised_points); return r; } @@ -391,6 +464,59 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, int lnx; int n, offset; int tmp_order = 2; + + CCTK_REAL *old_mask; + CCTK_INT *new_mask; + + CCTK_INT excised_points[11]; + + if (*excision_active < 0) + { + SetupExcisionVars((cGH*)cctkGH); + } + + if (old_excision) + { + old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi); + } + else + { + old_mask = NULL; + } + + if (new_excision) + { + new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi); + } + else + { + new_mask = NULL; + } + + /* + If this boundary point is excised then we return + without doing anything + */ + + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + if (old_mask[ijk] < 0.75) + { + return f[ijk]; + } + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + if (SpaceMask_CheckStateBits(new_mask, (ijk), + *new_excision_field, + *new_excision_descriptor)) + { + return f[ijk]; + } + } + } lnx = cctkGH->cctk_lsh[0]; @@ -424,11 +550,92 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, for (n = 0; n < 2 * eno_order + 1; ++n) { y[n] = f[ijk-offset+n]; + excised_points[n] = 0; + if (*excision_active) + { + if ((*excision_active == 1)||(*excision_active == 3)) + { + excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0; + } + if ((*excision_active == 2)||(*excision_active == 3)) + { + excised_points[n] |= + (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n), + *new_excision_field, + *new_excision_descriptor)); + } + } } /* call interpolator */ /* r = interpolate_eno(tmp_order, x0, dx0, y, x); */ - r = interpolate_eno(eno_order, x0, dx0, y, x); + r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points); return r; } + +void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + *excision_active = -1; + *old_mask_vi = -1; + *new_mask_vi = -1; + *new_excision_field = -1; + *new_excision_descriptor = -1; + +} + +void SetupExcisionVars(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (old_excision) + { + *old_mask_vi = CCTK_VarIndex(old_style_excision_var); + } + else + { + *old_mask_vi = -1; + } + if (new_excision) + { + *new_mask_vi = CCTK_VarIndex(new_style_excision_var); + } + else + { + *new_mask_vi = -1; + } + + if (*old_mask_vi > -1) + { + if (*new_mask_vi > -1) + { + *excision_active = 3; + *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name); + *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name, + new_mask_excised_name); + } + else + { + *excision_active = 1; + } + } + else + { + if (*new_mask_vi > -1) + { + *excision_active = 2; + *new_excision_field = SpaceMask_GetTypeBits(new_mask_field_name); + *new_excision_descriptor = SpaceMask_GetStateBits(new_mask_field_name, + new_mask_excised_name); + } + else + { + *excision_active = 0; + } + } + +} |