[cig-commits] commit: Make P_Refine work in 3D

[Mercurial]

changeset:   163:e29160299bda
user:        Walter Landry <wlandry at caltech.edu>
date:        Fri Apr 15 20:31:02 2011 -0700
files:       P_Refine.C
description:
Make P_Refine work in 3D


diff -r d5b47f287653 -r e29160299bda P_Refine.C
--- a/P_Refine.C	Fri Apr 15 20:30:39 2011 -0700
+++ b/P_Refine.C	Fri Apr 15 20:31:02 2011 -0700
@@ -72,78 +72,48 @@ void SAMRAI::geom::P_Refine::refine(
    tbox::Pointer<geom::CartesianPatchGeometry>
      geom = coarse_patch.getPatchGeometry();
 
-   for(int j=fine_box.lower(1); j<=fine_box.upper(1); ++j)
-     for(int i=fine_box.lower(0); i<=fine_box.upper(0); ++i)
-       {
-         pdat::CellIndex fine(tbox::Dimension(2));
-         fine[0]=i;
-         fine[1]=j;
+   for(pdat::CellIterator ci(fine_box); ci; ci++)
+     {
+       pdat::CellIndex fine(*ci);
 
-         pdat::CellIndex ip(fine), jp(fine);
-         ip[0]=1;
-         ip[1]=0;
-         jp[0]=0;
-         jp[1]=1;
-         pdat::CellIndex center(hier::Index::coarsen(fine,hier::Index(2,2)));
-         double dp_dx,dp_dy;
+       pdat::CellIndex
+         center(hier::Index::coarsen(fine,hier::Index::getOneIndex(dim)*2));
 
-         /* Pressure is cell-centered, so prolongation is a
-            linear interpolation from nearby cells. */
+       /* Pressure is cell-centered, so prolongation is a
+          linear interpolation from nearby cells. */
 
-         /* This assumes that the levels are always properly nested,
-            so that we always have an extra grid space for
-            interpolation.  So we only have to have a special case for
-            physical boundaries, where we do not have an extra grid
-            space. */
+       /* This assumes that the levels are always properly nested,
+          so that we always have an extra grid space for
+          interpolation.  So we only have to have a special case for
+          physical boundaries, where we do not have an extra grid
+          space. */
 
-         /* We could, in theory, use a refine_patch_strategy to fill
-            in the ghost zones with extrapolations, and then use
-            simple differences here. */
+       /* We could, in theory, use a refine_patch_strategy to fill
+          in the ghost zones with extrapolations, and then use
+          simple differences here. */
 
-         if(center[0]==coarse_box.lower(0)
-            && geom->getTouchesRegularBoundary(0,0))
-           {
-             dp_dx=((*p)(center+ip)-(*p)(center))/4;
-           }
-         else if(center[0]==coarse_box.upper(0)
-                 && geom->getTouchesRegularBoundary(0,1))
-           {
-             dp_dx=((*p)(center)-(*p)(center-ip))/4;
-           }
-         else
-           {
-             dp_dx=((*p)(center+ip)-(*p)(center-ip))/8;
-           }
-
-         if(center[1]==coarse_box.lower(1)
-            && geom->getTouchesRegularBoundary(1,0))
-           {
-             dp_dy=((*p)(center+jp)-(*p)(center))/4;
-           }
-         else if(center[1]==coarse_box.upper(1)
-                 && geom->getTouchesRegularBoundary(1,1))
-           {
-             dp_dy=((*p)(center)-(*p)(center-jp))/4;
-           }
-         else
-           {
-             dp_dy=((*p)(center+jp)-(*p)(center-jp))/8;
-           }
-
-         (*p_fine)(fine)=(*p)(center)
-           + ((i%2==0) ? (-dp_dx) : dp_dx)
-           + ((j%2==0) ? (-dp_dy) : dp_dy);
-
-         // tbox::plog << "P_Refine "
-         //            << fine_patch.getPatchLevelNumber() << " "
-         //            << i << " "
-         //            << j << " "
-         //            << (*p_fine)(fine) << " "
-         //            << (*p)(center) << " "
-         //            << dp_dx << " "
-         //            << dp_dy << " "
-         //            << "\n";
-
+       (*p_fine)(fine)=(*p)(center);
+       for(int d=0; d<dim.getValue(); ++d)
+         {
+           hier::Index ip(hier::Index::getZeroIndex(dim));
+           ip[d]=1;
+           double dp;
+           if(center[d]==coarse_box.lower(d)
+              && geom->getTouchesRegularBoundary(d,0))
+             {
+               dp=((*p)(center+ip)-(*p)(center))/4;
+             }
+           else if(center[d]==coarse_box.upper(d)
+                   && geom->getTouchesRegularBoundary(d,1))
+             {
+               dp=((*p)(center)-(*p)(center-ip))/4;
+             }
+           else
+             {
+               dp=((*p)(center+ip)-(*p)(center-ip))/8;
+             }
+           (*p_fine)(fine)+=((fine[d]%2==0) ? (-dp) : dp);
+         }           
        }
 }
 #endif