C++ SideIterator类代码示例

void EBPoissonOp::
getOpVoFStencil(VoFStencil&             a_stencil,
                const int&              a_idir,
                const Vector<VolIndex>& a_allMonotoneVoFs,
                const EBISBox&          a_ebisbox,
                const VolIndex&         a_VoF,
                const bool&             a_lowOrder)
{

  for (SideIterator sit; sit.ok(); ++sit)
    {
      Side::LoHiSide side = sit();
      Vector<FaceIndex> faces;

      faces = a_ebisbox.getFaces(a_VoF,a_idir,side);

      for (int iface = 0; iface < faces.size(); iface++)
        {
          FaceIndex face = faces[iface];
          VoFStencil faceStencil;

          getOpFaceStencil(faceStencil,a_allMonotoneVoFs,a_ebisbox,a_VoF,
                           a_idir,side,face,a_lowOrder);

          a_stencil += faceStencil;
        }
    }
}

// -----------------------------------------------------------------------------
// Interpolate ghosts at CF interface using zeros on coarser grids.
// -----------------------------------------------------------------------------
void homogeneousCFInterp (LevelData<FArrayBox>& a_phif,
                          const RealVect&       a_fineDx,
                          const RealVect&       a_crseDx,
                          const CFRegion&       a_cfRegion,
                          const IntVect&        a_applyDirs)
{
    CH_TIME("homogeneousCFInterp (full level)");

    // Loop over grids, directions, and sides and call the worker function.
    DataIterator dit = a_phif.dataIterator();
    for (dit.begin(); dit.ok(); ++dit) {
        if (a_phif[dit].box().isEmpty()) continue;

        for (int dir = 0; dir < SpaceDim; ++dir) {
            if (a_applyDirs[dir] == 0) continue;

            SideIterator sit;
            for (sit.begin(); sit.ok(); sit.next()) {
                homogeneousCFInterp(a_phif,
                                    dit(),
                                    dir,
                                    sit(),
                                    a_fineDx[dir],
                                    a_crseDx[dir],
                                    a_cfRegion);
            }
        }
    }
}

void 
CompGridVTOBC::operator()( FArrayBox&           a_state,
			   const Box&           a_valid,
			   const ProblemDomain& a_domain,
			   Real                 a_dx,
			   bool                 a_homogeneous)
{
  const Box& domainBox = a_domain.domainBox();
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      if (!a_domain.isPeriodic(idir))
        {
          for (SideIterator sit; sit.ok(); ++sit)
            {
              Side::LoHiSide side = sit();              
              if (a_valid.sideEnd(side)[idir] == domainBox.sideEnd(side)[idir])
                {
                  // Dirichlet BC
                  int isign = sign(side);
                  Box toRegion = adjCellBox(a_valid, idir, side, 1);
                  // include corner cells if possible by growing toRegion in transverse direction
                  toRegion.grow(1);
                  toRegion.grow(idir, -1);
                  toRegion &= a_state.box();
                  for (BoxIterator bit(toRegion); bit.ok(); ++bit)
                    {
                      IntVect ivTo = bit();                     
                      IntVect ivClose = ivTo - isign*BASISV(idir);
                      for (int ighost=0;ighost<m_nGhosts[0];ighost++,ivTo += isign*BASISV(idir))
                        {
                          //for (int icomp = 0; icomp < a_state.nComp(); icomp++) a_state(ivTo, icomp) = 0.0;
                          IntVect ivFrom = ivClose;
                          
                          // hardwire to linear BCs for now
                          for (int icomp = 0; icomp < a_state.nComp() ; icomp++)
                            {
                              if (m_bcDiri[idir][side][icomp]) 
                                {
                                  a_state(ivTo, icomp) = (-1.0)*a_state(ivFrom, icomp);
                                }
                              else 
                                {
                                  a_state(ivTo, icomp) = (1.0)*a_state(ivFrom, icomp);
                                }
                            }   
                        }
                    } // end loop over cells
                } // if ends match
            } // end loop over sides
        } // if not periodic in this direction
    } // end loop over directions    
}

void NodeNeumBC(NodeFArrayBox&  a_state,
                const Box&      a_valid,
                Real            a_dx,
                bool            a_homogeneous,
                BCValueHolder   a_value)
{
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      for (SideIterator sit; sit.ok(); ++sit)
        {
          NodeNeumBC(a_state, a_valid, a_dx, a_homogeneous, a_value, idir,sit());
        }
    }
}

void
EBCompositeMACProjector::
correctVelocityComponent(Vector<LayoutData< Vector< BaseIVFAB<Real> * > >* >      &  a_coveredVelLo,
                         Vector<LayoutData< Vector< BaseIVFAB<Real> * > >* >      &  a_coveredVelHi,
                         const Vector< LayoutData< Vector< Vector<VolIndex> > >* >&  a_coveredFaceLo,
                         const Vector< LayoutData< Vector< Vector<VolIndex> > >* >&  a_coveredFaceHi,
                         const Vector< LayoutData< Vector< IntVectSet > >* >      &  a_coveredSetsLo,
                         const Vector< LayoutData< Vector< IntVectSet > >* >      &  a_coveredSetsHi,
                         const Vector<LevelData<EBFluxFAB>* >                     &  a_macGradient,
                         int                                                         a_coveredFaceDir,
                         int                                                         a_velComp)
{
  CH_TIME("EBCompositeMACProjector::correctVelocityComponent");
  for (int ilev = 0; ilev < m_numLevels; ilev++)
    {
      for (SideIterator sit; sit.ok(); ++sit)
        {
          LayoutData<Vector<BaseIVFAB<Real>* > >*        velPtr = NULL;
          const LayoutData<Vector<Vector<VolIndex> > >* facePtr = NULL;
          const LayoutData<Vector<IntVectSet>  >*       setsPtr = NULL;
          if (sit() == Side::Lo)
            {
              velPtr  =  a_coveredVelLo[ilev];
              facePtr = a_coveredFaceLo[ilev];
              setsPtr = a_coveredSetsLo[ilev];
            }
          else
            {
              velPtr  =  a_coveredVelHi[ilev];
              facePtr = a_coveredFaceHi[ilev];
              setsPtr = a_coveredSetsHi[ilev];
            }

          const LevelData<EBFluxFAB>& macGradient = *a_macGradient[ilev];
          for (DataIterator dit = m_eblg[ilev].getDBL().dataIterator(); dit.ok(); ++dit)
            {
              const EBFluxFAB       & macGradFAB  = macGradient[dit()];
              const Vector<VolIndex>& coveredFace =  (*facePtr)[dit()][a_coveredFaceDir];
              const IntVectSet      & coveredSets =  (*setsPtr)[dit()][a_coveredFaceDir];
              BaseIVFAB<Real>       & coveredVel  = *((*velPtr)[dit()][a_coveredFaceDir]);
              const EBISBox& ebisBox = m_eblg[ilev].getEBISL()[dit()];

              correctVelocityComponent(coveredVel, coveredFace, coveredSets, macGradFAB, ebisBox,
                                       a_coveredFaceDir, sit(), a_velComp);
            }
        }
    }
}

void
LevelFluxRegisterEdge::incrementFine(
                                     FArrayBox& a_fineFlux,
                                     Real a_scale,
                                     const DataIndex& a_fineDataIndex,
                                     const Interval& a_srcInterval,
                                     const Interval& a_dstInterval)
{
  CH_assert(isDefined());
  CH_assert(!a_fineFlux.box().isEmpty());
  CH_assert(a_srcInterval.size() == a_dstInterval.size());
  CH_assert(a_srcInterval.begin() >= 0);
  CH_assert(a_srcInterval.end() < a_fineFlux.nComp());
  CH_assert(a_dstInterval.begin() >= 0);
  CH_assert(a_dstInterval.end() < m_nComp);

  int edgeDir = -1;
  for (int sideDir = 0; sideDir<SpaceDim; sideDir++)
    {
      if (a_fineFlux.box().type(sideDir) == IndexType::CELL)
        {
          edgeDir = sideDir;
        }
    }
  CH_assert(edgeDir >= 0);
  CH_assert(edgeDir < SpaceDim);

  for (int faceDir=0; faceDir<SpaceDim; faceDir++)
    {
      if (faceDir != edgeDir)
        {

          SideIterator sit;
          for (sit.begin(); sit.ok(); ++sit)
            {
              incrementFine(a_fineFlux,
                            a_scale,
                            a_fineDataIndex,
                            a_srcInterval,
                            a_dstInterval,
                            faceDir,
                            sit());
            }
        }
    }
}

void setBorkedFlux(EBFluxFAB&          a_flux,
                   const EBISBox&      a_ebisBox,
                   const Box&          a_box,
                   const RealVect&     a_fluxVal,
                   const BaseFab<int>& a_map)
{
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      Real bogval = 0;
      //the bogus value will be zero so it will not
      //do to have the correct value be zero
      if (Abs(a_fluxVal[idir]) < 1.0e-3) bogval = 1.0;

      a_flux[idir].setVal(a_fluxVal[idir]);
      //if i am in a 1d box and the box next to me
      //is 2d, set the border flux to
      for (SideIterator sit; sit.ok(); ++sit)
        {
          Box borderBox;
          if (sit() == Side::Lo)
            {
              borderBox = adjCellLo(a_box, idir,  -1);
            }
          else
            {
              borderBox = adjCellHi(a_box, idir,  -1);
            }
          for (BoxIterator bit(borderBox); bit.ok(); ++bit)
            {
              const IntVect& thisIV = bit();
              IntVect thatIV = thisIV + sign(sit())*BASISV(idir);
              if ((a_map(thisIV, 0) == 1) && (a_map(thatIV, 0) == 2))
                {
                  Vector<FaceIndex> faces = a_ebisBox.getAllFaces(thisIV, idir, sit());
                  for (int iface = 0; iface < faces.size(); iface++)
                    {
                      a_flux[idir](faces[iface],0) = bogval;
                    }
                }
            }
        }
    }

}

void
EBFluxRegister::
incrementCoarseIrregular(const BaseIFFAB<Real>& a_coarFlux,
                         const Real&            a_scale,
                         const DataIndex&       a_coarDatInd,
                         const Interval&        a_variables,
                         const int&             a_dir)
{
  if (m_hasEBCF)
    {
      EBFaceFAB facefab;
      bool hasCells = copyBIFFToEBFF(facefab, a_coarFlux, m_eblgCoar.getDBL()[a_coarDatInd], m_eblgCoar.getEBISL()[a_coarDatInd]);
      if (hasCells)
        {
          for (SideIterator sit; sit.ok(); ++sit)
            {
              incrementCoarIrreg(facefab, a_scale, a_coarDatInd, a_variables, a_dir, sit());
            }
        }
    }
}

void
LevelFluxRegisterEdge::setToZero()
{

  for (DataIterator dit = m_regCoarse.dataIterator(); dit.ok(); ++dit)
    m_regCoarse[dit()].setVal(0.0);


  SideIterator side;
  for (int idir=0 ; idir<SpaceDim; ++idir)
  {
    for (side.begin(); side.ok(); ++side)
    {
      LevelData<FluxBox>& fineReg = m_fabFine[index(idir, side())];


      for (DataIterator dit = fineReg.dataIterator(); dit.ok(); ++dit)
        fineReg[dit()].setVal(0.0);

    }
  }
}

Real
divergence(const EBFluxFAB&  a_func,
           const RealVect&   a_bndryFlux,
           const EBISBox&    a_ebisBox,
           const VolIndex&   a_vof,
           const Real&       a_dx)
{
  Real retval = 0;

  Real bndryArea = a_ebisBox.bndryArea(a_vof);
  RealVect normal = a_ebisBox.normal(a_vof);
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      Real bndryFlux    = a_bndryFlux[idir]; //normal already dealt with
      Real bndryContrib = -bndryFlux*bndryArea*normal[idir];
      Real openContrib = 0;
      if (bndryArea > 1.0e-3)
        {
          openContrib = 0;
        }
      for (SideIterator sit; sit.ok(); ++sit)
        {
          int isign = sign(sit());
          Real rsign = isign;
          Vector<FaceIndex> faces = a_ebisBox.getFaces(a_vof, idir, sit());
          for (int iface = 0; iface < faces.size(); ++iface)
            {
              Real areaFrac = a_ebisBox.areaFrac(faces[iface]);
              Real faceFlux = a_func[idir](faces[iface], 0);
              openContrib += rsign*faceFlux*areaFrac;
            }
        }
      retval += openContrib + bndryContrib;

    }
  retval /= a_dx;
  return retval;
}

void
kappaDivergence(EBCellFAB&             a_divF,
                const EBFluxFAB&       a_flux,
                const EBISBox&         a_ebisBox,
                const Box&             a_box,
                const Real&            a_dx)
{
    //set the divergence initially to zero
    //then loop through directions and increment the divergence
    //with each directions flux difference.
    a_divF.setVal(0.0);
    BaseFab<Real>&       regDivF = a_divF.getSingleValuedFAB();
    regDivF.setVal(0.);
    for (int idir = 0; idir < SpaceDim; idir++)
    {
        //update for the regular vofs in the nonconservative
        //case  works for all single valued vofs.
        /* do the regular vofs */
        /**/
        const EBFaceFAB& fluxDir = a_flux[idir];
        const BaseFab<Real>& regFluxDir = fluxDir.getSingleValuedFAB();
        int ncons = 1;
        FORT_DIVERGEF( CHF_BOX(a_box),
                       CHF_FRA(regDivF),
                       CHF_CONST_FRA(regFluxDir),
                       CHF_CONST_INT(idir),
                       CHF_CONST_INT(ncons),
                       CHF_CONST_REAL(a_dx));
        /**/
    }
    //update the irregular vofs using conservative diff
    IntVectSet ivsIrreg = a_ebisBox.getIrregIVS(a_box);
    for (VoFIterator vofit(ivsIrreg, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
    {
        const VolIndex& vof = vofit();
        //divergence was set in regular update.  we reset it
        // to zero and recalc.
        Real update = 0.;
        for ( int idir = 0; idir < SpaceDim; idir++)
        {
            const EBFaceFAB& fluxDir = a_flux[idir];
            for (SideIterator sit; sit.ok(); ++sit)
            {
                int isign = sign(sit());
                Vector<FaceIndex> faces =
                    a_ebisBox.getFaces(vof, idir, sit());
                for (int iface = 0; iface < faces.size(); iface++)
                {
                    const FaceIndex& face = faces[iface];
                    Real areaFrac = a_ebisBox.areaFrac(face);
                    Real faceFlux =fluxDir(face, 0);
                    update += isign*areaFrac*faceFlux;

                }
            }
        }
        //add EB boundary condtions in divergence
        const IntVect& iv = vof.gridIndex();
        Real bndryArea = a_ebisBox.bndryArea(vof);
        RealVect bndryCent = a_ebisBox.bndryCentroid(vof);
        RealVect normal = a_ebisBox.normal(vof);
        RealVect bndryLoc;
        RealVect exactF;
        for (int idir = 0; idir < SpaceDim; idir++)
        {
            bndryLoc[idir] = a_dx*(iv[idir] + 0.5 + bndryCent[idir]);
        }
        for (int idir = 0; idir < SpaceDim; idir++)
        {
            exactF[idir] = exactFlux(bndryLoc, idir);
        }
        Real bndryFlux = PolyGeom::dot(exactF, normal);

        update -= bndryFlux*bndryArea;
        update /= a_dx; //note NOT divided by volfrac

        a_divF(vof, 0) = update;
    }
}

// this function averages down the fine solution to the valid
// regions of the computed solution, then subtracts ir from
// the computed solution.  (error is exact-computed)
void computeAMRError(Vector<LevelData<FArrayBox>* >&       a_error,
                     const Vector<string>&                 a_errorVars,
                     const Vector<LevelData<FArrayBox>* >& a_computedSoln,
                     const Vector<string>&                 a_computedVars,
                     const Vector<DisjointBoxLayout>&      a_computedGrids,
                     const Real                            a_computedDx,
                     const Vector<int>&                    a_computedRefRatio,
                     const Vector<LevelData<FArrayBox>* >& a_exactSoln,
                     const Vector<string>&                 a_exactVars,
                     const Real                            a_exactDx,
                     Real                                  a_bogus_value,
                     bool                                  a_HOaverage,
                     bool                                  a_computeRelativeError)
{
  int numLevels = a_computedSoln.size();

  CH_assert(a_exactSoln.size() == 1);
  CH_assert(a_error.size() == numLevels);
  CH_assert(a_exactDx <= a_computedDx);
  CH_assert(a_computedRefRatio.size() >= numLevels - 1);

  if (a_exactDx == a_computedDx)
  {
    cerr << "Exact dx and computed dx are equal." << endl;
  }

  // check whether input file selects "sum all variables"
  bool sumAll = false;
  ParmParse pp;
  pp.query("sumAll",sumAll);
  
  // const DisjointBoxLayout& exactGrids = a_exactSoln[0]->getBoxes();

  Real dxLevel = a_computedDx;

  // do a bit of sleight-of-hand in the case where there are no
  // ghost cells in the exact solution -- allocate a temporary which
  // _has_ ghost cells, and do a copyTo
  LevelData<FArrayBox>* exactSolnPtr = NULL;
  bool allocatedMemory = false;
  if (a_exactSoln[0]->ghostVect() == IntVect::Zero)
    {
      exactSolnPtr = new LevelData<FArrayBox>(a_exactSoln[0]->getBoxes(),
                                              a_exactSoln[0]->nComp(),
                                              IntVect::Unit);
      a_exactSoln[0]->copyTo(*exactSolnPtr);

      allocatedMemory = true;
     }
   else
     {
       // if there are ghost cells, we can use the exactSoln as-is
       exactSolnPtr = a_exactSoln[0];
     }
   LevelData<FArrayBox>& exactSolnRef = *exactSolnPtr;

   // first need to set boundary conditions on exactsoln
   // this is for the Laplacian which is needed in AverageHO
   DataIterator ditFine = exactSolnRef.dataIterator();
   DomainGhostBC exactBC;
   Interval exactComps(0, a_exactVars.size() - 1);
   for (int dir = 0; dir < SpaceDim; dir++)
   {
     SideIterator sit;
     for (sit.reset(); sit.ok(); ++sit)
     {
       // use HO extrapolation at physical boundaries
       HOExtrapBC thisBC(dir, sit(), exactComps);
       exactBC.setBoxGhostBC(thisBC);
     }
   }

   for (ditFine.begin(); ditFine.ok(); ++ditFine)
   {
     FArrayBox& thisFineSoln = exactSolnRef[ditFine()];
     const Box& fineBox = exactSolnRef.getBoxes()[ditFine()];
     exactBC.applyInhomogeneousBCs(thisFineSoln, fineBox, a_exactDx);
   }
   exactSolnRef.exchange(exactComps);

   // outer loop is over levels
   for (int level = 0; level < numLevels; level++)
   {
     LevelData<FArrayBox>& thisLevelError = *a_error[level];
     LevelData<FArrayBox>& thisLevelComputed = *a_computedSoln[level];

     // compute refinement ratio between solution at this level
     // and exact solution
     Real nRefTemp = (dxLevel / a_exactDx);
     int nRefExact = (int) nRefTemp;

     // this is to do rounding properly if necessary
     if (nRefTemp - nRefExact > 0.5) nRefExact += 1;

     // make sure it's not zero
     if (nRefExact == 0) nRefExact =1;

//.........这里部分代码省略.........

void
EBLevelTGA::
setSourceGhostCells(LevelData<EBCellFAB>&    a_src,
                    const DisjointBoxLayout& a_grids,
                    int a_lev)
{
  int ncomp = a_src.nComp();
  for (DataIterator dit = a_grids.dataIterator(); dit.ok(); ++dit)
    {
      const Box& grid   = a_grids.get(dit());
      const Box& srcBox = a_src[dit()].box();
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          for (SideIterator sit; sit.ok(); ++sit)
            {
              int iside = sign(sit());
              Box bc_box = adjCellBox(grid, idir, sit(), 1);

              for (int jdir = 0; jdir < SpaceDim; jdir++)
                {
                  //want corners too
                  if (jdir != idir)
                    {
                      bc_box.grow(jdir, 1);
                    }
                }

              //if fails might not have a ghost cell.
              bc_box &= m_eblg[a_lev].getDomain().domainBox();

              CH_assert(srcBox.contains(bc_box));
              if (grid.size(idir) >= 4)
                {
                  FORT_HORESGHOSTBC(CHF_FRA(a_src[dit()].getSingleValuedFAB()),
                                    CHF_BOX(bc_box),
                                    CHF_CONST_INT(idir),
                                    CHF_CONST_INT(iside),
                                    CHF_CONST_INT(ncomp));
                }
              else
                {
                  // valid region not wide enough to apply HOExtrap -- drop
                  // to linear extrap
                  FORT_RESGHOSTBC(CHF_FRA(a_src[dit()].getSingleValuedFAB()),
                                  CHF_BOX(bc_box),
                                  CHF_CONST_INT(idir),
                                  CHF_CONST_INT(iside),
                                  CHF_CONST_INT(ncomp));
                }
              IntVectSet ivs = m_eblg[a_lev].getEBISL()[dit()].getIrregIVS(bc_box);
              for (VoFIterator vofit(ivs, m_eblg[a_lev].getEBISL()[dit()].getEBGraph()); vofit.ok(); ++vofit)
                {
                  for (int icomp = 0; icomp < ncomp; icomp++)
                    {
                      Real valNeigh = 0;
                      Vector<FaceIndex> faces = m_eblg[a_lev].getEBISL()[dit()].getFaces(vofit(), idir, flip(sit()));
                      for (int iface = 0; iface < faces.size(); iface++)
                        {
                          VolIndex vofNeigh = faces[iface].getVoF(flip(sit()));
                          valNeigh += a_src[dit()](vofNeigh, icomp);
                        }
                      if (faces.size() > 1) valNeigh /= faces.size();
                      a_src[dit()](vofit(), icomp) = valNeigh;
                    }
                }
            }
        }
    }
}

void
EBCompositeMACProjector::
kappaDivergence(Vector<LevelData<EBCellFAB>* >&              a_divu,
                Vector<LevelData<EBFluxFAB>* >&              a_velo,
                const Vector<LevelData<BaseIVFAB<Real> >* >* a_boundaryVelo)
{
  CH_TIME("EBCompositeMACProjector::kappaDivergence");
  for (int ilev = 0; ilev < m_numLevels; ilev++)
    {
      EBLevelMACProjector::setCurLevel(ilev);

      EBFluxFactory fluxfact(m_eblg[ilev].getEBISL());
      //need one ghost cell so that exchange is meaningful
      LevelData<EBFluxFAB> centroidVelocity(m_eblg[ilev].getDBL(), 1, IntVect::Unit, fluxfact);

      Interval interv(0, 0);
      LevelData<EBFluxFAB>& velExch = (LevelData<EBFluxFAB>&)(*a_velo[ilev]);
      velExch.exchange(interv);

      //interpolate velocity to face centroids
      EBArith::interpolateFluxToCentroids(centroidVelocity, *a_velo[ilev],
                                          m_eblg[ilev].getDBL(), m_eblg[ilev].getEBISL(), m_eblg[ilev].getDomain());

      centroidVelocity.exchange(interv);
      LevelData<BaseIVFAB<Real> >*  levelBoundaryVel = NULL;
      if (a_boundaryVelo != NULL)
        {
          levelBoundaryVel = (*a_boundaryVelo)[ilev];
        }
      //compute the divergence igoring other levels
      macKappaDivergence(*a_divu[ilev], centroidVelocity,
                         m_eblg[ilev].getDBL(), m_eblg[ilev].getEBISL(),
                         m_eblg[ilev].getDomain(),   m_dx[ilev],
                         levelBoundaryVel);

      //use flux registers to correct divergence
      //at coarse-fine interface with velocity from finer
      //level
      if (ilev < (m_numLevels-1))
        {
          Real incrScale = 1.0;
          m_fluxReg[ilev]->setToZero();
          //increment with coarse values
          for (DataIterator dit = m_eblg[ilev].getDBL().dataIterator();dit.ok(); ++dit)
            {
              const EBFluxFAB& veloFlux = (*a_velo[ilev])[dit()];
              for (int idir = 0; idir < SpaceDim; idir++)
                {
                  // This assumes that embedded boundaries and coarse-fine boundaries do not cross.
                  // To remove this assumption use incrementCoarseBoth.
                  m_fluxReg[ilev]->incrementCoarseRegular(veloFlux[idir], incrScale, dit(), interv, idir);
                }
            }
          //increment with fine velocities
          for (DataIterator dit = m_eblg[ilev+1].getDBL().dataIterator();dit.ok(); ++dit)
            {
              const EBFluxFAB& veloFlux = (*a_velo[ilev+1])[dit()];
              for (int idir = 0; idir < SpaceDim; idir++)
                {
                  for (SideIterator sit; sit.ok(); ++sit)
                    {
                      // This assumes that embedded boundaries and coarse-fine boundaries do not cross.
                      // To remove this assumption use incrementFineBoth.
                      m_fluxReg[ilev]->incrementFineRegular(veloFlux[idir], incrScale, dit(), interv, idir, sit());
                    }
                }
            }
          //reflux
          Real reflScale = 1.0/m_dx[ilev][0];
          m_fluxReg[ilev]->reflux(*a_divu[ilev], interv, reflScale);
        }
    }
}

void
EBFluxRegister::
incrementRedistRegister(EBCoarToCoarRedist& a_register,
                        const Interval&     a_variables,
                        const Real&         a_scale)
{
  if (m_hasEBCF)
    {
      LevelData<BaseIVFAB<Real> >& registerMass = a_register.m_regsCoar;
      LayoutData<IntVectSet>&      registerSets = a_register.m_setsCoar;
      //reflux into an empty LevelData<EBCellFAB>
      //Multiply this by (kappa)(1-kappa)
      //add result into register mass
      EBCellFactory ebcfCoar(m_eblgCoar.getEBISL());
      LevelData<EBCellFAB> increment(m_eblgCoar.getDBL(), m_nComp, m_saveCoar.ghostVect(), ebcfCoar);
      EBLevelDataOps::clone (increment, m_saveCoar);
      EBLevelDataOps::setVal(increment, 0.0);

      reflux(increment, a_variables, a_scale, true);

      for (DataIterator dit = m_eblgCoar.getDBL().dataIterator(); dit.ok(); ++dit)
        {
          for (int idir = 0; idir < SpaceDim; idir++)
            {
              for (SideIterator sit; sit.ok(); ++sit)
                {
                  int iindex = index(idir, sit());
                  Vector<IntVectSet> setsCoar = (m_setsCoar[iindex])[dit()];
                  for (int iset = 0; iset < setsCoar.size(); iset++)
                    {
                      const IntVectSet& setCoa  = registerSets[dit()];
                      const IntVectSet& setReg  = setsCoar[iset];
                      IntVectSet set = setCoa;
                      set &= setReg;
                      const EBISBox& ebisBox =m_eblgCoar.getEBISL()[dit()];
                      for (VoFIterator vofit(set, ebisBox.getEBGraph()); vofit.ok(); ++vofit)
                        {
                          const VolIndex& vof = vofit();

                          int ibleck = 0;
                          if ((vof.gridIndex() == ivdebugfr) && EBFastFR::s_verbose)
                            {
                              ibleck = 1;
                              pout()    << setprecision(10)
                                        << setiosflags(ios::showpoint)
                                        << setiosflags(ios::scientific);
                              pout() << "incrcotoco:" << endl;

                            }

                          for (int icomp = a_variables.begin(); icomp <= a_variables.end(); icomp++)
                            {
                              Real extraMass = increment[dit()](vofit(), icomp);
                              Real oldMass   =  registerMass[dit()](vofit(), icomp);
                              Real newMass   = oldMass + extraMass;
                              Real diff = extraMass;
                              if (ibleck == 1)
                                {
                                  pout() << "( " << oldMass << ", " << newMass << ", "  << diff << ")";
                                }

                              registerMass[dit()](vofit(), icomp) += extraMass;

                              //set increment to zero in case it gets
                              //hit twice (more than one direction or
                              //whatever
                              increment[dit()](vof, icomp) = 0;

                            } //loop over comps
                          if (ibleck == 1)
                            {
                              ibleck = 0;
                              pout() << endl;
                            }
                        }//loop over vofs in the set
                    }//loop over sets in this coarse box
                } //loop over sides
            } //loop over directions
        } //dataiterator loop
    } //You are using Bonetti's defense against me, uh?
} //I thought it fitting, considering the rocky terrain.