Provide a coordinate-thorn and symmetry-thorn independent way of

specifying the shape of the physical domain and the location of the
boundary points.  The advantage of that is that it is now much easier
to change the resolution.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/CoordBase/trunk@7 0337457d-221f-4ee6-a5f0-14255d5370d8

1 files changed, 338 insertions, 0 deletions

diff --git a/src/Domain.c b/src/Domain.c
new file mode 100644
index 0000000..ab917a9
--- /dev/null
+++ b/src/Domain.c
@@ -0,0 +1,338 @@
+/* $Header$ */
+
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+
+
+CCTK_INT CoordBase_GetBoundarySpecification
+  (CCTK_INT const size,
+   CCTK_INT * const nboundaryzones,
+   CCTK_INT * const is_internal,
+   CCTK_INT * const is_staggered,
+   CCTK_INT * const shiftout)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  assert (size>=0);
+  assert (nboundaryzones);
+  assert (is_internal);
+  assert (is_staggered);
+  assert (shiftout);
+  
+  assert (size == 6);
+  
+  nboundaryzones[0] = boundary_size_x_lower;
+  nboundaryzones[1] = boundary_size_x_upper;
+  nboundaryzones[2] = boundary_size_y_lower;
+  nboundaryzones[3] = boundary_size_y_upper;
+  nboundaryzones[4] = boundary_size_z_lower;
+  nboundaryzones[5] = boundary_size_z_upper;
+  
+  is_internal[0] = boundary_internal_x_lower;
+  is_internal[1] = boundary_internal_x_upper;
+  is_internal[2] = boundary_internal_y_lower;
+  is_internal[3] = boundary_internal_y_upper;
+  is_internal[4] = boundary_internal_z_lower;
+  is_internal[5] = boundary_internal_z_upper;
+  
+  is_staggered[0] = boundary_staggered_x_lower;
+  is_staggered[1] = boundary_staggered_x_upper;
+  is_staggered[2] = boundary_staggered_y_lower;
+  is_staggered[3] = boundary_staggered_y_upper;
+  is_staggered[4] = boundary_staggered_z_lower;
+  is_staggered[5] = boundary_staggered_z_upper;
+  
+  shiftout[0] = boundary_shiftout_x_lower;
+  shiftout[1] = boundary_shiftout_x_upper;
+  shiftout[2] = boundary_shiftout_y_lower;
+  shiftout[3] = boundary_shiftout_y_upper;
+  shiftout[4] = boundary_shiftout_z_lower;
+  shiftout[5] = boundary_shiftout_z_upper;
+  
+  return 0;
+}
+
+
+
+CCTK_INT CoordBase_GetDomainSpecification
+  (CCTK_INT const size,
+   CCTK_REAL * const physical_min,
+   CCTK_REAL * const physical_max,
+   CCTK_REAL * const interior_min,
+   CCTK_REAL * const interior_max,
+   CCTK_REAL * const exterior_min,
+   CCTK_REAL * const exterior_max,
+   CCTK_REAL * const thespacing)
+{
+  DECLARE_CCTK_PARAMETERS;
+  
+  int d;
+  int ierr;
+  
+  assert (size>=0);
+  assert (physical_min);
+  assert (physical_max);
+  assert (interior_min);
+  assert (interior_max);
+  assert (exterior_min);
+  assert (exterior_max);
+  assert (thespacing);
+  
+  assert (size == 3);
+  
+  
+  
+  if (CCTK_EQUALS (domainsize, "minmax")) {
+    
+    physical_min[0] = xmin;
+    physical_min[1] = ymin;
+    physical_min[2] = zmin;
+    physical_max[0] = xmax;
+    physical_max[1] = ymax;
+    physical_max[2] = zmax;
+    if (CCTK_EQUALS (spacing, "gridspacing")) {
+      thespacing[0] = dx;
+      thespacing[1] = dy;
+      thespacing[2] = dz;
+    } else if (CCTK_EQUALS (spacing, "numcells")) {
+      thespacing[0] = (physical_max[0] - physical_min[0]) / ncells_x;
+      thespacing[1] = (physical_max[1] - physical_min[1]) / ncells_y;
+      thespacing[2] = (physical_max[2] - physical_min[2]) / ncells_z;
+    }
+    
+  } else if (CCTK_EQUALS (domainsize, "extent")) {
+    
+    if (zero_origin_x) {
+      physical_min[0] = xmin;
+      physical_max[0] = xmin + xextent;
+    } else {
+      physical_min[0] = - xextent / 2;
+      physical_max[0] = + xextent / 2;
+    }
+    if (zero_origin_y) {
+      physical_min[1] = ymin;
+      physical_max[1] = ymin + yextent;
+    } else {
+      physical_min[1] = - yextent / 2;
+      physical_max[1] = + yextent / 2;
+    }
+    if (zero_origin_z) {
+      physical_min[2] = zmin;
+      physical_max[2] = zmin + zextent;
+    } else {
+      physical_min[2] = - zextent / 2;
+      physical_max[2] = + zextent / 2;
+    }
+    if (CCTK_EQUALS (spacing, "gridspacing")) {
+      thespacing[0] = dx;
+      thespacing[1] = dy;
+      thespacing[2] = dz;
+    } else if (CCTK_EQUALS (spacing, "numcells")) {
+      thespacing[0] = (physical_max[0] - physical_min[0]) / ncells_x;
+      thespacing[1] = (physical_max[1] - physical_min[1]) / ncells_y;
+      thespacing[2] = (physical_max[2] - physical_min[2]) / ncells_z;
+    }
+    
+  } else if (CCTK_EQUALS (domainsize, "spacing")) {
+    
+    if (zero_origin_x) {
+      physical_min[0] = xmin;
+      physical_max[0] = xmin + dx * ncells_x;
+    } else {
+      physical_min[0] = - dx * ncells_x / 2;
+      physical_max[0] = + dx * ncells_x / 2;
+    }
+    if (zero_origin_y) {
+      physical_min[1] = ymin;
+      physical_max[1] = ymin + dy * ncells_y;
+    } else {
+      physical_min[1] = - dy * ncells_y / 2;
+      physical_max[1] = + dy * ncells_y / 2;
+    }
+    if (zero_origin_z) {
+      physical_min[2] = - dz * ncells_z / 2;
+      physical_max[2] = + dz * ncells_z / 2;
+    } else {
+      physical_min[2] = zmin;
+      physical_max[2] = zmin + dz * ncells_z;
+    }
+    thespacing[0] = dx;
+    thespacing[1] = dy;
+    thespacing[2] = dz;
+    
+  } else {
+    assert (0);
+  }
+  
+  ierr = ConvertFromPhysicalBoundary
+    (size, physical_min, physical_max, interior_min, interior_max,
+     exterior_min, exterior_max, thespacing);
+  assert (!ierr);
+  
+  return 0;
+}
+
+
+
+CCTK_INT CoordBase_ConvertFromPhysicalBoundary
+  (CCTK_INT const size,
+   CCTK_REAL const * const physical_min,
+   CCTK_REAL const * const physical_max,
+   CCTK_REAL * const interior_min,
+   CCTK_REAL * const interior_max,
+   CCTK_REAL * const exterior_min,
+   CCTK_REAL * const exterior_max,
+   CCTK_REAL const * const thespacing)
+{
+  CCTK_INT nboundaryzones[6];
+  CCTK_INT is_internal[6];
+  CCTK_INT is_staggered[6];
+  CCTK_INT shiftout[6];
+  
+  int d;
+  int ierr;
+  
+  assert (size>=0);
+  assert (physical_min);
+  assert (physical_max);
+  assert (interior_min);
+  assert (interior_max);
+  assert (exterior_min);
+  assert (exterior_max);
+  assert (thespacing);
+  
+  assert (size == 3);
+  
+  ierr = CoordBase_GetBoundarySpecification
+    (6, nboundaryzones, is_internal, is_staggered, shiftout);
+  assert (!ierr);
+  
+  for (d=0; d<3; ++d) {
+    
+    exterior_min[d] = physical_min[d] - thespacing[d] *
+      (+ (is_internal[2*d] ? 0 : nboundaryzones[2*d] - 1)
+       + (is_staggered[2*d] ? 0.5 : 0.0)
+       + shiftout[2*d]);
+    exterior_max[d] = physical_max[d] + thespacing[d] *
+      (+ (is_internal[2*d+1] ? 0 : nboundaryzones[2*d+1] - 1)
+       + (is_staggered[2*d+1] ? 0.5 : 0.0)
+       + shiftout[2*d+1]);
+  
+    interior_min[d] = exterior_min[d] + thespacing[d] * nboundaryzones[2*d];
+    interior_max[d] = exterior_max[d] - thespacing[d] * nboundaryzones[2*d+1];
+
+  }
+  
+  return 0;
+}
+
+
+
+CCTK_INT CoordBase_ConvertFromInteriorBoundary
+  (CCTK_INT const size,
+   CCTK_REAL * const physical_min,
+   CCTK_REAL * const physical_max,
+   CCTK_REAL const * const interior_min,
+   CCTK_REAL const * const interior_max,
+   CCTK_REAL * const exterior_min,
+   CCTK_REAL * const exterior_max,
+   CCTK_REAL const * const thespacing)
+{
+  CCTK_INT nboundaryzones[6];
+  CCTK_INT is_internal[6];
+  CCTK_INT is_staggered[6];
+  CCTK_INT shiftout[6];
+  
+  int d;
+  int ierr;
+  
+  assert (size>=0);
+  assert (physical_min);
+  assert (physical_max);
+  assert (interior_min);
+  assert (interior_max);
+  assert (exterior_min);
+  assert (exterior_max);
+  assert (thespacing);
+  
+  assert (size == 3);
+  
+  ierr = CoordBase_GetBoundarySpecification
+    (6, nboundaryzones, is_internal, is_staggered, shiftout);
+  assert (!ierr);
+  
+  for (d=0; d<3; ++d) {
+    
+    exterior_min[d] = interior_min[d] - thespacing[d] * nboundaryzones[2*d];
+    exterior_max[d] = interior_max[d] + thespacing[d] * nboundaryzones[2*d+1];
+
+    physical_min[d] = exterior_min[d] + thespacing[d] *
+      (+ (is_internal[2*d] ? 0 : nboundaryzones[2*d] - 1)
+       + (is_staggered[2*d] ? 0.5 : 0.0)
+       + shiftout[2*d]);
+    physical_max[d] = exterior_max[d] - thespacing[d] *
+      (+ (is_internal[2*d+1] ? 0 : nboundaryzones[2*d+1] - 1)
+       + (is_staggered[2*d+1] ? 0.5 : 0.0)
+       + shiftout[2*d+1]);
+  
+  }
+  
+  return 0;
+}
+
+
+
+CCTK_INT CoordBase_ConvertFromExteriorBoundary
+  (CCTK_INT const size,
+   CCTK_REAL * const physical_min,
+   CCTK_REAL * const physical_max,
+   CCTK_REAL * const interior_min,
+   CCTK_REAL * const interior_max,
+   CCTK_REAL const * const exterior_min,
+   CCTK_REAL const * const exterior_max,
+   CCTK_REAL const * const thespacing)
+{
+  CCTK_INT nboundaryzones[6];
+  CCTK_INT is_internal[6];
+  CCTK_INT is_staggered[6];
+  CCTK_INT shiftout[6];
+  
+  int d;
+  int ierr;
+  
+  assert (size>=0);
+  assert (physical_min);
+  assert (physical_max);
+  assert (interior_min);
+  assert (interior_max);
+  assert (exterior_min);
+  assert (exterior_max);
+  assert (thespacing);
+  
+  assert (size == 3);
+  
+  ierr = CoordBase_GetBoundarySpecification
+    (6, nboundaryzones, is_internal, is_staggered, shiftout);
+  assert (!ierr);
+  
+  for (d=0; d<3; ++d) {
+    
+    physical_min[d] = exterior_min[d] + thespacing[d] *
+      (+ (is_internal[2*d] ? 0 : nboundaryzones[2*d] - 1)
+       + (is_staggered[2*d] ? 0.5 : 0.0)
+       + shiftout[2*d]);
+    physical_max[d] = exterior_max[d] - thespacing[d] *
+      (+ (is_internal[2*d+1] ? 0 : nboundaryzones[2*d+1] - 1)
+       + (is_staggered[2*d+1] ? 0.5 : 0.0)
+       + shiftout[2*d+1]);
+  
+    interior_min[d] = exterior_min[d] + thespacing[d] * nboundaryzones[2*d];
+    interior_max[d] = exterior_max[d] - thespacing[d] * nboundaryzones[2*d+1];
+
+  }
+  
+  return 0;
+}