C++ SideIterator::begin方法代码示例

// -----------------------------------------------------------------------------
// 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
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
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);

    }
  }
}

void
LevelFluxRegisterEdge::refluxCurl(LevelData<FluxBox>& a_uCoarse,
                                  Real a_scale)
{
  CH_assert(isDefined());
  CH_assert(a_uCoarse.nComp() == m_nComp);

  SideIterator side;
  // idir is the normal direction to the coarse-fine interface
  for (int idir=0 ; idir<SpaceDim; ++idir)
  {
    for (side.begin(); side.ok(); ++side)
    {
      LevelData<FluxBox>& fineReg = m_fabFine[index(idir, side())];

      // first, create temp LevelData<FluxBox> to hold "coarse flux"
      const DisjointBoxLayout coarseBoxes = m_regCoarse.getBoxes();

      // this fills the place of what used to be m_fabCoarse in the old
      // implementation
      LevelData<FluxBox> coarReg(coarseBoxes, m_nComp, IntVect::Unit);


      // now fill the coarReg with the curl of the stored coarse-level
      // edge-centered flux
      DataIterator crseDit = coarseBoxes.dataIterator();
      for (crseDit.begin(); crseDit.ok(); ++crseDit)
        {
          FluxBox& thisCoarReg = coarReg[crseDit];
          thisCoarReg.setVal(0.0);

          EdgeDataBox& thisEdgeData = m_regCoarse[crseDit];

          for (int edgeDir=0; edgeDir<SpaceDim; edgeDir++)
            {
              if (idir != edgeDir)
                {
                  FArrayBox& crseEdgeDataDir = thisEdgeData[edgeDir];
                  for (int faceDir = 0; faceDir<SpaceDim; faceDir++)
                    {
                      if (faceDir != edgeDir)
                        {
                          FArrayBox& faceData = thisCoarReg[faceDir];
                          int shiftDir = -1;
                          for (int i=0; i<SpaceDim; i++)
                            {
                              if ((i != faceDir) && (i != edgeDir) )
                                {
                                  shiftDir = i;
                                }
                            }
                          CH_assert(shiftDir >= 0);
                          crseEdgeDataDir.shiftHalf(shiftDir, sign(side()));
                          // scaling already taken care of in incrementCrse
                          Real scale = 1.0;
                          faceData.plus(crseEdgeDataDir, scale, 0, 0, faceData.nComp());
                          crseEdgeDataDir.shiftHalf(shiftDir, -sign(side()));
                        } // end if not normal direction
                    } // end loop over face directions
                } // end if edgeDir != idir
            } // end loop over edge directions
        } // end loop over crse boxes


      // first, we need to create a temp LevelData<FluxBox>
      // to make a local copy in the coarse layout space of
      // the fine register increments

      LevelData<FluxBox> fineRegLocal(coarReg.getBoxes(), m_nComp, IntVect::Unit);

      fineReg.copyTo(fineReg.interval(), fineRegLocal,
                     fineRegLocal.interval(),
                     m_crseCopiers[index(idir,side())]);

      for (DataIterator it = a_uCoarse.dataIterator(); it.ok(); ++it)
        {
          // loop over flux components here
          for (int fluxComp=0; fluxComp < SpaceDim; fluxComp++)
            {
              // we don't do anything in the normal direction
              if (fluxComp != idir)
                {
                  // fluxDir is the direction of the face-centered flux
                  FArrayBox& U = a_uCoarse[it()][fluxComp];
                  // set up IntVectSet to avoid double counting of updates
                  Box coarseGridBox = U.box();
                  // transfer to Cell-centered, then create IVS
                  coarseGridBox.shiftHalf(fluxComp,1);
                  IntVectSet nonUpdatedEdges(coarseGridBox);

                  // remember, we want to take the curl here
                  // also recall that fluxComp is the component
                  // of the face-centered curl (not the edge-centered
                  // vector field that we're refluxing, which is why
                  // the sign may seem like it's the opposite of what
                  // you might expect!
                  Real local_scale = -sign(side())*a_scale;
                  //int testDir = (fluxComp+1)%(SpaceDim);
                  if (((fluxComp+1)%(SpaceDim)) == idir)  local_scale *= -1;
                  Vector<IntVectSet>& ivsV =
//.........这里部分代码省略.........

// new define
void
LevelFluxRegisterEdge::define(
                              const DisjointBoxLayout& a_dbl,
                              const DisjointBoxLayout& a_dblCoarse,
                              const ProblemDomain& a_dProblem,
                              int a_nRefine,
                              int a_nComp)
{
  m_isDefined = true;
  CH_assert(a_nRefine > 0);
  CH_assert(a_nComp > 0);
  CH_assert(!a_dProblem.isEmpty());
  m_nComp = a_nComp;
  m_nRefine = a_nRefine;
  m_domainCoarse = coarsen(a_dProblem, a_nRefine);
  CH_assert (a_dblCoarse.checkPeriodic(m_domainCoarse));

  // allocate copiers
  m_crseCopiers.resize(SpaceDim*2);

  SideIterator side;

  // create a Vector<Box> of the fine boxes which also includes periodic images,
  // since we don't really care about the processor layouts, etc
  Vector<Box> periodicFineBoxes;
  CFStencil::buildPeriodicVector(periodicFineBoxes, a_dProblem, a_dbl);
  // now coarsen these boxes...
  for (int i=0; i<periodicFineBoxes.size(); i++)
    {
      periodicFineBoxes[i].coarsen(m_nRefine);
    }

  for (int idir=0 ; idir<SpaceDim; ++idir)
  {
    for (side.begin(); side.ok(); ++side)
    {
      // step one, build fineBoxes, flux register boxes
      // indexed by the fine level but in the coarse index
      // space
      DisjointBoxLayout fineBoxes,tmp;
      // first create coarsened dbl, then compute flux register boxes
      // adjacent to coarsened fine boxes
      coarsen(tmp, a_dbl, m_nRefine);
      if (side() == Side::Lo)
        {
          adjCellLo(fineBoxes, tmp, idir,1);
        }
      else
        {
          adjCellHi(fineBoxes, tmp, idir,1);
        }

      // now define the FluxBoxes of fabFine on this DisjointBoxLayout
      m_fabFine[index(idir, side())].define(fineBoxes, a_nComp);



      LayoutData<Vector<Vector<IntVectSet> > >& ivsetsVect
        = m_refluxLocations[index(idir, side())];
      ivsetsVect.define(a_dblCoarse);

      LayoutData<Vector<DataIndex> >& mapsV =
        m_coarToCoarMap[index(idir, side())];
      mapsV.define(a_dblCoarse);

      DisjointBoxLayout coarseBoxes = a_dblCoarse;
      DataIterator dit = a_dblCoarse.dataIterator();
      for (dit.begin(); dit.ok(); ++dit)
        {
          unsigned int thisproc = a_dblCoarse.procID(dit());
          if (thisproc == procID())
          {
            ivsetsVect[DataIndex(dit())].resize(SpaceDim);
          }
          const Box& coarseBox = a_dblCoarse[dit()];
          int count = 0;
          for (int i=0; i<periodicFineBoxes.size(); i++)
            {
              Box regBox;
              if (side() == Side::Lo)
                {
                  regBox = adjCellLo(periodicFineBoxes[i], idir, 1);
                }
              else
                {
                  regBox = adjCellHi(periodicFineBoxes[i], idir, 1);
                }

              // do this little dance in order to ensure that
              // we catch corner cells which might be in different
              // boxes.
              Box testBox(regBox);
              testBox.grow(1);
              testBox.grow(idir,-1);
              if (testBox.intersectsNotEmpty(coarseBox))
                {
                  testBox &= coarseBox;
                  ++count;
                  unsigned int proc = a_dblCoarse.procID(dit());
//.........这里部分代码省略.........

void VCAMRPoissonOp2::reflux(const LevelData<FArrayBox>&        a_phiFine,
                            const LevelData<FArrayBox>&        a_phi,
                            LevelData<FArrayBox>&              a_residual,
                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp)
{
  CH_TIME("VCAMRPoissonOp2::reflux");

  int ncomp = 1;
  ProblemDomain fineDomain = refine(m_domain, m_refToFiner);
  LevelFluxRegister levfluxreg(a_phiFine.disjointBoxLayout(),
                               a_phi.disjointBoxLayout(),
                               fineDomain,
                               m_refToFiner,
                               ncomp);

  levfluxreg.setToZero();
  Interval interv(0,a_phi.nComp()-1);

  DataIterator dit = a_phi.dataIterator();
  for (dit.reset(); dit.ok(); ++dit)
    {
      const FArrayBox& coarfab = a_phi[dit];
      const FluxBox& coarBCoef  = (*m_bCoef)[dit];
      const Box& gridBox = a_phi.getBoxes()[dit];

      for (int idir = 0; idir < SpaceDim; idir++)
        {
          FArrayBox coarflux;
          Box faceBox = surroundingNodes(gridBox, idir);
          getFlux(coarflux, coarfab, coarBCoef , faceBox, idir);

          Real scale = 1.0;
          levfluxreg.incrementCoarse(coarflux, scale,dit(),
                                     interv,interv,idir);
        }
    }
  LevelData<FArrayBox>& p = ( LevelData<FArrayBox>&)a_phiFine;

  // has to be its own object because the finer operator
  // owns an interpolator and we have no way of getting to it
  VCAMRPoissonOp2* finerAMRPOp = (VCAMRPoissonOp2*) a_finerOp;
  QuadCFInterp& quadCFI = finerAMRPOp->m_interpWithCoarser;

  quadCFI.coarseFineInterp(p, a_phi);
  // p.exchange(a_phiFine.interval()); // BVS is pretty sure this is not necesary.
  IntVect phiGhost = p.ghostVect();

  DataIterator ditf = a_phiFine.dataIterator();
  const  DisjointBoxLayout& dblFine = a_phiFine.disjointBoxLayout();
  for (ditf.reset(); ditf.ok(); ++ditf)
    {
      const FArrayBox& phifFab = a_phiFine[ditf];
      const FluxBox& fineBCoef  = (*(finerAMRPOp->m_bCoef))[ditf];
      const Box& gridbox = dblFine.get(ditf());
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          int normalGhost = phiGhost[idir];
          SideIterator sit;
          for (sit.begin(); sit.ok(); sit.next())
            {
              Side::LoHiSide hiorlo = sit();
              Box fabbox;
              Box facebox;

              // assumption here that the stencil required
              // to compute the flux in the normal direction
              // is 2* the number of ghost cells for phi
              // (which is a reasonable assumption, and probably
              // better than just assuming you need one cell on
              // either side of the interface
              // (dfm 8-4-06)
              if (sit() == Side::Lo)
                {
                  fabbox = adjCellLo(gridbox,idir, 2*normalGhost);
                  fabbox.shift(idir, 1);
                  facebox = bdryLo(gridbox, idir,1);
                }
              else
                {
                  fabbox = adjCellHi(gridbox,idir, 2*normalGhost);
                  fabbox.shift(idir, -1);
                  facebox = bdryHi(gridbox, idir, 1);
                }

              // just in case we need ghost cells in the transverse direction
              // (dfm 8-4-06)
              for (int otherDir=0; otherDir<SpaceDim; ++otherDir)
                {
                  if (otherDir != idir)
                    {
                      fabbox.grow(otherDir, phiGhost[otherDir]);
                    }
                }
              CH_assert(!fabbox.isEmpty());

              FArrayBox phifab(fabbox, a_phi.nComp());
              phifab.copy(phifFab);

              FArrayBox fineflux;
              getFlux(fineflux, phifab, fineBCoef, facebox, idir,
//.........这里部分代码省略.........

//
// VCAMRPoissonOp2::reflux()
//   There are currently the new version (first) and the old version (second)
//   in this file.  Brian asked to preserve the old version in this way for
//   now. - TJL (12/10/2007)
//
void VCAMRPoissonOp2::reflux(const LevelData<FArrayBox>&        a_phiFine,
                            const LevelData<FArrayBox>&        a_phi,
                            LevelData<FArrayBox>&              a_residual,
                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp)
{
  CH_TIMERS("VCAMRPoissonOp2::reflux");

  m_levfluxreg.setToZero();
  Interval interv(0,a_phi.nComp()-1);

  CH_TIMER("VCAMRPoissonOp2::reflux::incrementCoarse", t2);
  CH_START(t2);

  DataIterator dit = a_phi.dataIterator();
  for (dit.reset(); dit.ok(); ++dit)
  {
    const FArrayBox& coarfab   = a_phi[dit];
    const FluxBox&   coarBCoef = (*m_bCoef)[dit];
    const Box&       gridBox   = a_phi.getBoxes()[dit];

    if (m_levfluxreg.hasCF(dit()))
    {
      for (int idir = 0; idir < SpaceDim; idir++)
      {
        FArrayBox coarflux;
        Box faceBox = surroundingNodes(gridBox, idir);

        getFlux(coarflux, coarfab, coarBCoef, faceBox, idir);

        Real scale = 1.0;
        m_levfluxreg.incrementCoarse(coarflux, scale,dit(),
            interv, interv, idir);
      }
    }
  }

  CH_STOP(t2);

  // const cast:  OK because we're changing ghost cells only
  LevelData<FArrayBox>& phiFineRef = ( LevelData<FArrayBox>&)a_phiFine;

  VCAMRPoissonOp2* finerAMRPOp = (VCAMRPoissonOp2*) a_finerOp;
  QuadCFInterp& quadCFI = finerAMRPOp->m_interpWithCoarser;

  quadCFI.coarseFineInterp(phiFineRef, a_phi);
  // I'm pretty sure this is not necessary. bvs -- flux calculations use
  // outer ghost cells, but not inner ones
  // phiFineRef.exchange(a_phiFine.interval());
  IntVect phiGhost = phiFineRef.ghostVect();
  int ncomps = a_phiFine.nComp();

  CH_TIMER("VCAMRPoissonOp2::reflux::incrementFine", t3);
  CH_START(t3);

  DataIterator ditf = a_phiFine.dataIterator();
  const DisjointBoxLayout& dblFine = a_phiFine.disjointBoxLayout();
  for (ditf.reset(); ditf.ok(); ++ditf)
    {
      const FArrayBox& phifFab   = a_phiFine[ditf];
      const FluxBox&   fineBCoef = (*(finerAMRPOp->m_bCoef))[ditf];
      const Box&       gridbox   = dblFine.get(ditf());

      for (int idir = 0; idir < SpaceDim; idir++)
        {
          //int normalGhost = phiGhost[idir];
          SideIterator sit;
          for (sit.begin(); sit.ok(); sit.next())
            {
              if (m_levfluxreg.hasCF(ditf(), sit()))
                {
                  Side::LoHiSide hiorlo = sit();
                  Box fluxBox = bdryBox(gridbox,idir,hiorlo,1);

                  FArrayBox fineflux(fluxBox,ncomps);
                  getFlux(fineflux, phifFab, fineBCoef, fluxBox, idir,
                          m_refToFiner);

                  Real scale = 1.0;
                  m_levfluxreg.incrementFine(fineflux, scale, ditf(),
                                             interv, interv, idir, hiorlo);
                }
            }
        }
    }

  CH_STOP(t3);

  Real scale = 1.0/m_dx;
  m_levfluxreg.reflux(a_residual, scale);
}