C++ DisjointBoxLayout::dataIterator方法代码示例

void
PoisselleTube::
initializeVelocity(LevelData<EBCellFAB>&    a_velocity,
                   const DisjointBoxLayout& a_grids,
                   const EBISLayout&        a_ebisl,
                   const ProblemDomain&     a_domain,
                   const RealVect&          a_origin,
                   const Real&              a_time,
                   const RealVect&          a_dx) const
{
  for (DataIterator dit = a_grids.dataIterator(); dit.ok(); ++dit)
    {
      IntVectSet ivsBox(a_grids.get(dit()));
      for (VoFIterator vofit(ivsBox, a_ebisl[dit()].getEBGraph()); vofit.ok(); ++vofit)
        {
          RealVect xval = EBArith::getVofLocation(vofit() , a_dx, RealVect::Zero);
          for (int idir = 0; idir < SpaceDim; idir++)
            {
              //bogus values for normal and time because they do not matter
              Real velComp = m_bcval[idir].value(xval, RealVect::Zero, a_time, idir);
              a_velocity[dit()](vofit(), idir) = velComp;
            }
        }
    }
}

void
NoFlowVortex::
initializeVelocity(LevelData<EBCellFAB>&    a_velocity,
                   const DisjointBoxLayout& a_grids,
                   const EBISLayout&        a_ebisl,
                   const ProblemDomain&     a_domain,
                   const RealVect&          a_origin,
                   const Real&              a_time,
                   const RealVect&          a_dx) const
{
  Real pi = 4.0*atan(1.0);
  for (DataIterator dit = a_grids.dataIterator(); dit.ok(); ++dit)
    {
      IntVectSet ivsBox(a_grids.get(dit()));
      for (VoFIterator vofit(ivsBox, a_ebisl[dit()].getEBGraph()); vofit.ok(); ++vofit)
        {
          RealVect xval, velpt;
          getXVal(xval, a_origin,vofit(), a_dx);
          getVelPt(velpt, xval, pi);
          for (int idir = 0; idir < SpaceDim; idir++)
            {
              a_velocity[dit()](vofit(), idir) = velpt[idir];
            }
        }
    }
}

void
getError(LevelData<EBCellFAB>&       a_error,
         const EBISLayout&           a_ebisl,
         const DisjointBoxLayout&    a_dbl,
         const Real&                 a_dx)
{
    EBCellFactory ebcellfact(a_ebisl);
    EBFluxFactory ebfluxfact(a_ebisl);
    a_error.define(a_dbl, 1, IntVect::Zero,   ebcellfact);
    LevelData<EBCellFAB> divFCalc(a_dbl, 1, IntVect::Zero,   ebcellfact);
    LevelData<EBCellFAB> divFExac(a_dbl, 1, IntVect::Zero,   ebcellfact);
    LevelData<EBFluxFAB> faceFlux(a_dbl, 1, IntVect::Zero,   ebfluxfact);

    for (DataIterator dit = a_dbl.dataIterator(); dit.ok(); ++dit)
    {
        a_error[dit()].setVal(0.);
        divFCalc[dit()].setVal(0.);
        divFExac[dit()].setVal(0.);
    }

    setToExactDivFLD(divFExac,  a_ebisl, a_dbl, a_dx);
    setToExactFluxLD(faceFlux,  a_ebisl, a_dbl, a_dx);

    Interval interv(0, 0);
    faceFlux.exchange(interv);

    kappaDivergenceLD(divFCalc, faceFlux, a_ebisl, a_dbl, a_dx);

    for (DataIterator dit = a_dbl.dataIterator(); dit.ok(); ++dit)
    {
        EBCellFAB& errorFAB = a_error[dit()];
        EBCellFAB& exactFAB = divFExac[dit()];
        EBCellFAB& calcuFAB = divFCalc[dit()];

        errorFAB += calcuFAB;
        errorFAB -= exactFAB;
    }
}

void
setToExactFluxLD(LevelData<EBFluxFAB>&       a_flux,
                 const EBISLayout&           a_ebisl,
                 const DisjointBoxLayout&    a_dbl,
                 const Real&                 a_dx)
{
    for (DataIterator dit= a_dbl.dataIterator(); dit.ok(); ++dit)
    {
        setToExactFlux(a_flux[dit()],
                       a_ebisl[dit()],
                       a_dbl.get(dit()),
                       a_dx);
    }
}

void
setToExactDivFLD(LevelData<EBCellFAB>&       a_soln,
                 const EBISLayout&           a_ebisl,
                 const DisjointBoxLayout&    a_dbl,
                 const Real&                 a_dx)
{
    for (DataIterator dit= a_dbl.dataIterator(); dit.ok(); ++dit)
    {
        setToExactDivF(a_soln[dit()],
                       a_ebisl[dit()],
                       a_dbl.get(dit()),
                       a_dx);
    }
}

// Set up initial conditions
void AdvectTestIBC::initialize(LevelData<FArrayBox>& a_U)
{
  DisjointBoxLayout grids = a_U.disjointBoxLayout();
  for (DataIterator dit = grids.dataIterator(); dit.ok(); ++dit)
    {
      const Box& grid = grids.get(dit());
      FORT_ADVECTINITF(CHF_FRA1(a_U[dit()],0),
                       CHF_CONST_REALVECT(m_center),
                       CHF_CONST_REAL(m_size),
                       CHF_CONST_REAL(m_dx),
                       CHF_BOX(grid));

    }

}

void
kappaDivergenceLD(LevelData<EBCellFAB>&       a_divF,
                  const LevelData<EBFluxFAB>& a_flux,
                  const EBISLayout&           a_ebisl,
                  const DisjointBoxLayout&    a_dbl,
                  const Real&                 a_dx)
{
    for (DataIterator dit= a_dbl.dataIterator(); dit.ok(); ++dit)
    {
        kappaDivergence(a_divF[dit()],
                        a_flux[dit()],
                        a_ebisl[dit()],
                        a_dbl.get(dit()),
                        a_dx);
    }

}

// ------------------------------------------------------------
// version of 27 March 2003
Real
norm(const LevelData<NodeFArrayBox>& a_phi,
     const ProblemDomain& a_domain,
     const DisjointBoxLayout& a_finerGridsCoarsened,
     const LayoutData< Vector<Box> >& a_IVSVext,
     const LayoutData< Vector<Box> >& a_IVSVintFinerCoarsened,
     const int a_nRefFine,
     const Real a_dx,
     const Interval& a_comps,
     const int a_p,
     bool a_verbose)
{
  // Idea:  copy a_phi to temp, then zero out temp on:
  // - exterior nodes of grids at this level;
  // - projections of interior nodes of the finer grids.
  int ncomps = a_comps.size();
  const DisjointBoxLayout& grids = a_phi.getBoxes();
  LevelData<NodeFArrayBox> temp(grids, ncomps);

  // Copy a_phi to temp.
  Interval newcomps(0, ncomps-1);
  for (DataIterator dit(grids.dataIterator()); dit.ok(); ++dit)
    {
      const NodeFArrayBox& nfab = a_phi[dit()];
      const Box& bx = grids.get(dit());
      temp[dit()].copy(bx, newcomps, bx, nfab, a_comps);
    }

  // Zero out temp on exterior nodes.
  zeroBoundaryNodes(temp, a_IVSVext);

  // Define zeroCoarsened to be all zero on the coarsened finer grids.
  LevelData<NodeFArrayBox>
    zeroCoarsened(a_finerGridsCoarsened, ncomps, IntVect::Zero);
  for (DataIterator dit(a_finerGridsCoarsened.dataIterator()); dit.ok(); ++dit)
    zeroCoarsened[dit()].getFab().setVal(0.);

  // Set temp to zero on interior nodes of coarsened finer grids.
  copyInteriorNodes(temp, zeroCoarsened, a_IVSVintFinerCoarsened);

  Real normLevel = norm(temp, a_dx, a_p, newcomps, a_verbose);

  return normLevel;
}

int 
EBAMRTestCommon::
checkForZero(const LevelData<EBCellFAB>& a_errorVelo,
             const DisjointBoxLayout&    a_gridsFine,
             const EBISLayout&           a_ebislFine,
             const Box&                  a_domainFine,
             string a_funcname)
{
  int eekflag = 0;
  Real eps = 1.0e-8;
#ifdef CH_USE_FLOAT
  eps = 1.0e-3;
#endif

  for (DataIterator dit = a_gridsFine.dataIterator(); dit.ok(); ++dit)
    {
      Box grownBox = a_gridsFine.get(dit());
      grownBox.grow(1);
      grownBox &= a_domainFine;
      IntVectSet ivsBox(grownBox);
      for (VoFIterator vofit(ivsBox, a_ebislFine[dit()].getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();
          int ihere = 0;
          if (vof.gridIndex() == EBDebugPoint::s_ivd)
            {
              ihere = 1;
            }
          for (int ivar = 0; ivar < a_errorVelo.nComp(); ivar++)
            {
              Real errorIn = a_errorVelo[dit()](vof, ivar);
              if (Abs(errorIn) > eps)
                {
                  pout() << "check for zero error too big for test " << a_funcname << endl;
                  pout() << "ivar = " << ivar <<  endl;
                  pout() << "vof = " << vof.gridIndex() << " error = " << errorIn << endl;
                  return -1;
                }
            }
        }

    }
  return eekflag;
}

// 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
EBMGAverage::
defineStencils()
{
  CH_TIME("EBMGInterp::defineStencils");

  DisjointBoxLayout gridsStenCoar;
  EBISLayout        ebislStenCoar;
  DisjointBoxLayout gridsStenFine;
  EBISLayout        ebislStenFine;
  if (m_layoutChanged)
    {
      if (m_coarsenable)
        {
          gridsStenCoar = m_buffGrids;
          ebislStenCoar = m_buffEBISL;
          gridsStenFine = m_fineGrids;
          ebislStenFine = m_fineEBISL;
        }
      else
        {
          gridsStenCoar = m_coarGrids;
          ebislStenCoar = m_coarEBISL;
          gridsStenFine = m_buffGrids;
          ebislStenFine = m_buffEBISL;
        }

    }
  else
    {
      gridsStenCoar = m_coarGrids;
      ebislStenCoar = m_coarEBISL;
      gridsStenFine = m_fineGrids;
      ebislStenFine = m_fineEBISL;
    }

  {
    CH_TIME("graph walking");

    LayoutData<Vector<VoFStencil> > averageStencil;
    LayoutData<VoFIterator > vofItIrregCoar;
    LayoutData<VoFIterator > vofItIrregFine;
    m_averageEBStencil.define(gridsStenCoar);
    averageStencil.define(    gridsStenCoar);
    vofItIrregFine.define(    gridsStenCoar); //not wrong.  trust me.
    vofItIrregCoar.define(    gridsStenCoar);

    for (DataIterator dit = gridsStenCoar.dataIterator(); dit.ok(); ++dit)
      {
        Box refBox(IntVect::Zero, IntVect::Zero);
        refBox.refine(m_refRat);
        int numFinePerCoar = refBox.numPts();

        const     Box&     boxFine = gridsStenFine[dit()];
        const EBISBox& ebisBoxFine = ebislStenFine[dit()];
        const EBGraph& ebGraphFine = ebisBoxFine.getEBGraph();

        const         Box& boxCoar = gridsStenCoar[dit()];
        const EBISBox& ebisBoxCoar = ebislStenCoar[dit()];
        const EBGraph& ebGraphCoar = ebisBoxCoar.getEBGraph();

        IntVectSet notRegularCoar = ebisBoxCoar.getIrregIVS(boxCoar);
        vofItIrregCoar[dit()].define(notRegularCoar, ebGraphCoar);

        IntVectSet ivsFine = refine(notRegularCoar, m_refRat);
        vofItIrregFine[dit()].define(ivsFine, ebGraphFine);

        const Vector<VolIndex>& allCoarVofs = vofItIrregCoar[dit()].getVector();

        Vector<VoFStencil>& averageStencils = averageStencil[dit()];

        averageStencils.resize(allCoarVofs.size());

        for (int icoar = 0; icoar < allCoarVofs.size(); icoar++)
          {
            Vector<VolIndex> fineVofs;
            if (m_refRat > 2)
              {
                fineVofs = ebislStenCoar.refine(allCoarVofs[icoar], m_refRat, dit());
              }
            else
              {
                fineVofs = ebislStenCoar[dit()].refine(allCoarVofs[icoar]);
              }

            VoFStencil& averageStencil = averageStencils[icoar];

            for (int ifine = 0; ifine < fineVofs.size(); ifine++)
              {
                averageStencil.add(fineVofs[ifine], 1./numFinePerCoar);
              }
          }

        m_averageEBStencil[dit()] = RefCountedPtr<EBStencil>(new EBStencil(allCoarVofs, averageStencil[dit()], boxCoar,  boxFine, ebisBoxCoar,  ebisBoxFine,  m_ghost, m_ghost));
      }
  }
}

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;
                    }
                }
            }
        }
    }
}

int getError(LevelData<EBCellFAB>& a_errorFine,
             const EBISLayout& a_ebislFine,
             const DisjointBoxLayout& a_gridsFine,
             const Box& a_domainFine,
             const Real& a_dxFine,
             const EBISLayout& a_ebislCoar,
             const DisjointBoxLayout& a_gridsCoar,
             const Box& a_domainCoar,
             const Real& a_dxCoar,
             const int& a_refRatio)
{
  int eekflag = 0;
  int nvar = 1;

  EBCellFactory ebcellfactFine(a_ebislFine);
  EBCellFactory ebcellfactCoar(a_ebislCoar);
  LevelData<EBCellFAB> phiFine(a_gridsFine, nvar,
                               IntVect::Zero,ebcellfactFine);
  LevelData<EBCellFAB> phiCoar(a_gridsCoar, nvar,
                               IntVect::Zero,ebcellfactCoar);

  //fill phi fine and phiCExact
  //put phiFexact into a_errorFine
  for (DataIterator dit = a_gridsFine.dataIterator();
      dit.ok(); ++dit)
    {
      IntVectSet ivsBox(a_gridsFine.get(dit()));
      phiFine[dit()].setVal(0.0);
      EBCellFAB& phiFineFAB = a_errorFine[dit()];
      phiFineFAB.setCoveredCellVal(0.0,0);
      for (VoFIterator vofit(ivsBox, a_ebislFine[dit()].getEBGraph());
          vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();
          Real rightAns = exactFunc(vof.gridIndex(), a_dxFine);
          phiFineFAB(vof, 0) = rightAns;
        }
    }

  for (DataIterator dit = a_gridsCoar.dataIterator();
      dit.ok(); ++dit)
    {
      IntVectSet ivsBox(a_gridsCoar.get(dit()));
      EBCellFAB& phiCoarFAB = phiCoar[dit()];
      phiCoarFAB.setCoveredCellVal(0.0,0);
      for (VoFIterator vofit(ivsBox, a_ebislCoar[dit()].getEBGraph());
          vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();
          Real rightAns = exactFunc(vof.gridIndex(), a_dxCoar);
          phiCoarFAB(vof, 0) = rightAns;
        }
    }

  EBPWLFineInterp interpOp(a_gridsFine, a_gridsCoar,
                           a_ebislFine, a_ebislCoar,
                           a_domainCoar, a_refRatio, nvar);
  Interval zeroiv(0,0);
  interpOp.interpolate(phiFine, phiCoar, zeroiv);
  //error = phiC - phiCExact
  for (DataIterator dit = a_gridsFine.dataIterator();
      dit.ok(); ++dit)
    {
      EBCellFAB& errorFAB = a_errorFine[dit()];
      EBCellFAB& phiFineFAB = phiFine[dit()];

      errorFAB -= phiFineFAB;
    }

  return eekflag;
}

//----------------------------------------------------------------------------
void
EBNormalizeByVolumeFraction::
operator()(LevelData<EBCellFAB>& a_Q,
           const Interval& a_compInterval) const
{
  CH_TIME("EBNormalizer::operator()");
   // Endpoints of the given interval.
   int begin = a_compInterval.begin(), end = a_compInterval.end(),
       length = a_compInterval.size();

   // Loop over the EBISBoxes within our grid. The EB data structures are
   // indexed in the same manner as the non-EB data structures, so we piggy-
   // back the former on the latter.
   a_Q.exchange();
   EBISLayout ebisLayout = m_levelGrid.getEBISL();
   DisjointBoxLayout layout = m_levelGrid.getDBL();
   for (DataIterator dit = layout.dataIterator(); dit.ok(); ++dit)
   {
      const EBISBox& box = ebisLayout[dit()];
      EBCellFAB& QFAB = a_Q[dit()];

      // Go over the irregular cells in this box.
      const IntVectSet& irregCells = box.getIrregIVS(layout[dit()]);

      // The average has to be computed from the uncorrected data from all
      // the neighbors, so we can't apply the corrections in place. For now,
      // we stash them in a map.
      map<VolIndex, vector<Real> > correctedValues;
      for (VoFIterator vit(irregCells, box.getEBGraph()); vit.ok(); ++vit)
        {
          Real kappajSum = 0.0;
          vector<Real> kappajQjSum(length, 0.0);

          // Get all of the indices of the VoFs within a monotone path
          // radius of 1.
          VolIndex vofi = vit();
          Vector<VolIndex> vofjs;
          EBArith::getAllVoFsInMonotonePath(vofjs, vofi, box, 1);

          // Accumulate the contributions from the neighboring cells.
          for (unsigned int j = 0; j < vofjs.size(); ++j)
            {
              VolIndex vofj = vofjs[j];

              Real kappaj = box.volFrac(vofj);
              for (int icomp = begin; icomp <= end; ++icomp)
                {
                  kappajQjSum[icomp] += QFAB(vofj, icomp);
                }

              // Add this volume fraction to the sum.
              kappajSum += kappaj;
            }

          if (kappajSum > 0.)
            {
              // Normalize the quantity and stow it.
              vector<Real> correctedValue(length);
              //         Real kappai = box.volFrac(vofi);  //unused dtg
              for (int icomp = begin; icomp <= end; ++icomp)
                {
                  // correctedValue[icomp - begin] =
                  //    QFAB(vofi, icomp) + (1.0 - kappai) * kappajQjSum[icomp] / kappajSum;
                  correctedValue[icomp - begin] = kappajQjSum[icomp] / kappajSum;
                }
              correctedValues[vofi] = correctedValue;
            }
        }

      // Apply the corrections.
      for (map<VolIndex, vector<Real> >::const_iterator
           cit = correctedValues.begin(); cit != correctedValues.end(); ++cit)
      {
         for (int icomp = begin; icomp <= end; ++icomp)
         {
            QFAB(cit->first, icomp) = cit->second[icomp-begin];
         }
      }
   }
}

// this function works on two solutions on equivalent grids.
// It subtracts the computed solution from the exact solution
// if a_doGhostCells == true, then does box-by-box comparison,
// including ghost cells (boxes must be the same for each).
// Otherwise, only does this for valid cells, but boxes don't
// need to be the same.
void computeSameSizeError(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_dx,
                          const Vector<int>&                    a_refRatio,
                          const Vector<LevelData<FArrayBox>* >& a_exactSoln,
                          const Vector<string>&                 a_exactVars,
                          Real                                  a_bogus_value,
                          bool                                  a_computeRelativeError,
                          bool                                  a_doGhostCells)

{
  int numLevels = a_computedSoln.size();

  CH_assert(a_exactSoln.size() == numLevels);
  CH_assert(a_error.size() == numLevels);
  CH_assert(a_refRatio.size() >= numLevels - 1);

  Real dxLevel = a_dx;

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

    const DisjointBoxLayout levelGrids = thisLevelComputed.getBoxes();
    const DisjointBoxLayout exactGrids = thisLevelExact.getBoxes();

    DataIterator levelDit = levelGrids.dataIterator();
    for (levelDit.begin(); levelDit.ok(); ++levelDit)
    {
      // initialize error to a bogus value
      thisLevelError[levelDit()].setVal(a_bogus_value);
    }

    // loop over variables
    for (int nErr = 0; nErr < a_errorVars.size(); nErr++)
    {
      string thisErrVar = a_errorVars[nErr];
      bool done = false;

      // this is where things differ between the ghost-cell
      // and non-ghost-cell approach.
      if (a_doGhostCells)
      {
        // this is the older approach to things --
        // do everything grid-by-grid

        // first loop over exact variables
        for (int exactComp = 0; exactComp < a_exactVars.size(); exactComp++)
        {
          string thisExactVar = a_exactVars[exactComp];
          // check if this exact variable is "the one"
          if (thisExactVar == thisErrVar)
          {
            int computedComp = 0;

            // now loop over computed variables
            while (!done && (computedComp < a_computedVars.size()))
            {
              if (a_computedVars[computedComp] == thisErrVar)
              {
                // copy exact solution -> error
                // and then subtract computed solution
                DataIterator exactDit = thisLevelExact.dataIterator();
                for (levelDit.reset(); levelDit.ok(); ++levelDit)
                {
                  FArrayBox& thisComputed = thisLevelComputed[levelDit()];
                  FArrayBox& thisError = thisLevelError[levelDit()];
                  const Box& thisBox = levelGrids[levelDit()];

                  for (exactDit.begin(); exactDit.ok(); ++exactDit)
                  {
                    if (thisBox.contains(exactGrids[exactDit()]))
                    {
                      thisError.copy(thisLevelExact[exactDit()],
                                     exactComp, nErr, 1);
                    } // end if exact and computed boxes match
                  } // end loop over exact grids

                  if (a_computeRelativeError)
                  {
                    // do this a little strangely -- relative
                    // error is one - computed/exact.
                    thisError.divide(thisComputed, computedComp, nErr, 1);
                    thisError.invert(-1.0, nErr, 1);
                    thisError.plus(1.0, nErr, 1);
                  }
                  else
                  {
//.........这里部分代码省略.........