1 files changed, 293 insertions, 0 deletions
diff --git a/Carpet/CarpetRegrid2/src/boundary.cc b/Carpet/CarpetRegrid2/src/boundary.cc
new file mode 100644
index 000000000..79e595784
--- /dev/null
+++ b/Carpet/CarpetRegrid2/src/boundary.cc
@@ -0,0 +1,293 @@
+#include <cctk.h>
+#include <cctk_Parameters.h>
+
+#include <carpet.hh>
+
+#include "boundary.hh"
+
+
+
+namespace CarpetRegrid2 {
+  
+  using namespace Carpet;
+  
+  
+  
+  // Convert a coordinate location to an index location. For cell
+  // centring, shift upwards.
+  ivect
+  rpos2ipos (rvect const & rpos,
+             rvect const & origin, rvect const & scale,
+             gh const & hh, int const rl)
+  {
+    ivect const istride  = hh.baseextent(0,rl).stride();
+    ivect const bistride = hh.baseextent(0,0).stride();
+    
+    if (hh.refcent == cell_centered) {
+      assert (all (istride % 2 == 0));
+    }
+    
+#if 1
+    ivect const ipos =
+      hh.refcent == vertex_centered
+      ? ivect (floor (((rpos - origin) * scale                    ) / rvect(istride) + rvect(0.5))) * istride
+      : ivect (floor (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride)             )) * istride + istride/2 + bistride/2;
+#else
+    ivect const ipos =
+      hh.refcent == vertex_centered
+      ? ivect (floor (((rpos - origin) * scale                    ) / rvect(istride) + rvect(0.5))) * istride
+      : ivect (floor (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride) + rvect(0.5))) * istride + istride/2 + bistride/2;
+#endif
+    
+    return ipos;
+  }
+  
+  // Convert a coordinate location to an index location, rounding in
+  // the opposite manner as rpos2ipos. For cell centring, shift
+  // downwards instead of upwards.
+  ivect
+  rpos2ipos1 (rvect const & rpos,
+              rvect const & origin, rvect const & scale,
+              gh const & hh, int const rl)
+  {
+    ivect const istride  = hh.baseextent(0,rl).stride();
+    ivect const bistride = hh.baseextent(0,0).stride();
+    
+    if (hh.refcent == cell_centered) {
+      assert (all (istride % 2 == 0));
+    }
+    
+#if 1
+    ivect const ipos =
+      hh.refcent == vertex_centered
+      ? ivect (ceil (((rpos - origin) * scale                    ) / rvect(istride) - rvect(0.5))) * istride
+      : ivect (ceil (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride)             )) * istride - istride/2 + bistride/2;
+#else
+    ivect const ipos =
+      hh.refcent == vertex_centered
+      ? ivect (ceil (((rpos - origin) * scale                    ) / rvect(istride) - rvect(0.5))) * istride
+      : ivect (ceil (((rpos - origin) * scale - rvect(bistride/2)) / rvect(istride) - rvect(0.5))) * istride - istride/2 + bistride/2;
+#endif
+    
+    return ipos;
+  }
+  
+  // Convert an index location to a coordinate location
+  rvect
+  ipos2rpos (ivect const & ipos,
+             rvect const & origin, rvect const & scale,
+             gh const & hh, int const rl)
+  {
+    return rvect(ipos) / scale + origin;
+  }
+  
+  // Convert an index bbox to a coordinate bbox
+  rbbox
+  ibbox2rbbox (ibbox const & ib,
+               rvect const & origin, rvect const & scale,
+               gh const & hh, int const rl)
+  {
+    rvect const zero (0);
+    return rbbox (ipos2rpos (ib.lower() , origin, scale, hh, rl),
+                  ipos2rpos (ib.upper() , origin, scale, hh, rl),
+                  ipos2rpos (ib.stride(), zero  , scale, hh, rl));
+  }
+  
+  
+  
+  void
+  get_boundary_specification (jjvect & nboundaryzones,
+                              jjvect & is_internal,
+                              jjvect & is_staggered,
+                              jjvect & shiftout)
+  {
+    if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_GetBoundarySpecification
+        (Carpet::map, 2*dim,
+         & nboundaryzones[0][0],
+         & is_internal[0][0],
+         & is_staggered[0][0],
+         & shiftout[0][0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("GetBoundarySpecification")) {
+      CCTK_INT const ierr = GetBoundarySpecification
+        (2*dim,
+         & nboundaryzones[0][0],
+         & is_internal[0][0],
+         & is_staggered[0][0],
+         & shiftout[0][0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  void
+  get_physical_boundary (rvect & physical_lower,
+                         rvect & physical_upper,
+                         rvect & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    rvect exterior_lower, exterior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_GetDomainSpecification
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("GetDomainSpecification")) {
+      CCTK_INT const ierr = GetDomainSpecification
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  void
+  calculate_exterior_boundary (rvect const & physical_lower,
+                               rvect const & physical_upper,
+                               rvect       & exterior_lower,
+                               rvect       & exterior_upper,
+                               rvect const & spacing)
+  {
+    rvect interior_lower, interior_upper;
+    if (CCTK_IsFunctionAliased ("MultiPatch_ConvertFromPhysicalBoundary")) {
+      assert (Carpet::map >= 0);
+      CCTK_INT const ierr = MultiPatch_ConvertFromPhysicalBoundary
+        (Carpet::map, dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else if (CCTK_IsFunctionAliased ("ConvertFromPhysicalBoundary")) {
+      CCTK_INT const ierr =
+        ConvertFromPhysicalBoundary
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    } else {
+      assert (0);
+    }
+  }
+  
+  
+  
+  // Location and description of the outer boundary
+  domain_boundary::
+  domain_boundary (gh const& hh, dh const& dd)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (veryverbose) {
+      cout << "Determining domain outer boundary...\n";
+    }
+    
+    get_boundary_specification (nboundaryzones, is_internal,
+                                is_staggered, shiftout);
+    
+    boundary_staggering_mismatch =
+      xpose ((hh.refcent == vertex_centered) != (is_staggered == 0));
+#warning "TODO: This is too strict"
+    assert (all (all (not boundary_staggering_mismatch)));
+    
+    get_physical_boundary (physical_lower, physical_upper, spacing);
+    
+    // Adapt spacing for convergence level
+    spacing *= ipow (CCTK_REAL(mgfact), basemglevel);
+    
+    calculate_exterior_boundary (physical_lower, physical_upper,
+                                 exterior_lower, exterior_upper,
+                                 spacing);
+    
+    // The physical boundary
+    origin = exterior_lower;
+    scale = rvect (hh.baseextent(0,0).stride()) / spacing;
+    
+    // The location of the outermost grid points. For cell centring,
+    // these are 1/2 grid spacing inside of the boundary.
+    physical_ilower = rpos2ipos (physical_lower, origin, scale, hh, 0);
+    physical_iupper = rpos2ipos1 (physical_upper, origin, scale, hh, 0);
+  }
+  
+  level_boundary::
+  level_boundary (gh const& hh, dh const& dd, int const rl)
+    : domain_boundary (hh, dd)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": determining outer boundary...\n";
+    }
+    
+    level_physical_lower = physical_lower;
+    level_physical_upper = physical_upper;
+    level_spacing = spacing / rvect (hh.reffacts.at(rl));
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": physical coordinate boundary is at " << r2vect(level_physical_lower, level_physical_upper) << "\n";
+      cout << "Refinement level " << rl << ": spacing is " << level_spacing << "\n";
+    }
+    
+    calculate_exterior_boundary (level_physical_lower, level_physical_upper,
+                                 level_exterior_lower, level_exterior_upper,
+                                 level_spacing);
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": exterior coordinate boundary is at " << r2vect(level_exterior_lower, level_exterior_upper) << "\n";
+    }
+    
+    ibbox const& baseextent = hh.baseextent(0,rl);
+    ivect const istride = baseextent.stride();
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": stride is " << istride << "\n";
+    }
+    
+    // This is the location of the outermost grid points. For cell
+    // centring, these are 1/2 grid spacing inside of the boundary.
+    level_physical_ilower = rpos2ipos (physical_lower, origin, scale, hh, rl);
+    level_physical_iupper = rpos2ipos1 (physical_upper, origin, scale, hh, rl);
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": physical boundary is at " << i2vect(level_physical_ilower, level_physical_iupper) << "\n";
+      cout << "Refinement level " << rl << ": reconstructed physical coordinate boundary is at "
+           << r2vect(ipos2rpos(level_physical_ilower - (hh.refcent == cell_centered ? istride/2 : 0), origin, scale, hh, rl),
+                     ipos2rpos(level_physical_iupper + (hh.refcent == cell_centered ? istride/2 : 0), origin, scale, hh, rl)) << "\n";
+    }
+    
+    level_exterior_ilower =
+      rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
+    level_exterior_iupper =
+      rpos2ipos1 (level_exterior_upper, origin, scale, hh, rl);
+    assert (all (level_exterior_ilower >= baseextent.lower()));
+    assert (all (level_exterior_iupper <= baseextent.upper()));
+    if (veryverbose) {
+      cout << "Refinement level " << rl << ": exterior boundary is at " << i2vect(level_exterior_ilower, level_exterior_iupper) << "\n";
+      cout << "Refinement level " << rl << ": reconstructed exterior coordinate boundary is at "
+           << r2vect(ipos2rpos(level_exterior_ilower, origin, scale, hh, rl),
+                     ipos2rpos(level_exterior_iupper, origin, scale, hh, rl)) << "\n";
+    }
+    
+    // Find the minimum necessary distance away from the outer
+    // boundary due to buffer and ghost zones. This is e.g. the
+    // distance that the lower boundary of a bbox has to have from the
+    // lower boundary. This is in terms of grid points.
+    min_bnd_dist_away = dd.ghost_widths.at(rl);
+    // Find the minimum necessary distance from the outer boundary due
+    // to buffer and ghost zones. This is e.g. the distance that the
+    // upper boundary of a bbox has to have from the lower boundary.
+    // This is in terms of grid points.
+    min_bnd_dist_incl = dd.ghost_widths.at(rl);
+    // TODO: The above is required only near symmetry boundaries.
+  }
+  
+} // namespace CarpetRegrid2