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

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

}

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

int checkCoarseAssortment(const Box& a_domain)
{
  int retval = 0;
  const EBIndexSpace* const ebisPtr = Chombo_EBIS::instance();
  CH_assert(ebisPtr->isDefined());
  Box fineDomain = a_domain;

  int numLevels = ebisPtr->numLevels();
  for (int ilev = 1; ilev < numLevels; ilev++)
    {
      CH_assert(!fineDomain.isEmpty());
      Vector<Box> vbox(1, fineDomain);
      Vector<int> proc(1, 0);
      DisjointBoxLayout fineDBL(vbox, proc);
      EBISLayout fineEBISL;
      int nghost = 4;
      ebisPtr->fillEBISLayout(fineEBISL, fineDBL, fineDomain, nghost);
      Box coarDomain = coarsen(fineDomain, 2);
      DisjointBoxLayout coarDBL;
      coarsen(coarDBL, fineDBL, 2);
      EBISLayout coarEBISL;
      ebisPtr->fillEBISLayout(coarEBISL, coarDBL, coarDomain, nghost);

      for (DataIterator dit = fineDBL.dataIterator(); dit.ok(); ++dit)
        {
          retval = checkEBISBox(coarDBL.get(dit()), coarEBISL[dit()], fineEBISL[dit()]);
          if (retval != 0)
            {
              pout() << "problem in coarsening  " << fineDomain << " to " << coarDomain << endl;
              return retval;
            }
        }
      fineDomain.coarsen(2);
   }
return retval;
}

int
testIFFAB(const DisjointBoxLayout&              a_dbl,
          const EBISLayout       &              a_ebisl,
          const Box              &              a_domain,
          const Real             &              a_dx )
{
  int faceDir = 0;
  int nFlux = 1;
  LayoutData<IntVectSet>        irregSetsGrown;
  LevelData< BaseIFFAB<Real> > fluxInterpolant;

  EBArith::defineFluxInterpolant(fluxInterpolant,
                                 irregSetsGrown,
                                 a_dbl, a_ebisl, a_domain, nFlux, faceDir);

  //set source fab to right ans over set only on grids interior cells
  int ibox = 0;
  for (DataIterator dit = a_dbl.dataIterator(); dit.ok(); ++dit)
    {
      BaseIFFAB<Real>& srcFab = fluxInterpolant[dit()];
      srcFab.setVal(-1.0);
      IntVectSet ivsSmall = irregSetsGrown[dit()];
      const Box& grid = a_dbl.get(dit());
      ivsSmall &=  grid;
      for (FaceIterator faceit(ivsSmall, a_ebisl[dit()].getEBGraph(), faceDir, FaceStop::SurroundingWithBoundary);
          faceit.ok(); ++faceit)
        {
          srcFab(faceit(), 0) = rightAns(faceit());
        }
      ibox++;
    }

  //diagnostics
  if (g_diagnosticMode)
    {
      pout() << " diagnostics for processor " << procID() << endl;
      for (DataIterator dit = a_dbl.dataIterator(); dit.ok(); ++dit)
        {
          const IntVectSet& ivsGrown = irregSetsGrown[dit()];
          const Box& grid = a_dbl.get(dit());
          pout() << "============" << endl;
          pout() << " box = " << grid;
          pout() << ", full ivs = "  ;
          dumpIVS(&ivsGrown);
          pout() << "============" << endl;
          for (LayoutIterator lit = a_dbl.layoutIterator(); lit.ok(); ++lit)
            {
              const Box& grid2 = a_dbl.get(lit());
              IntVectSet ivsIntersect = ivsGrown;
              ivsIntersect &= grid2;
              pout() << "intersection with box " << grid2 << " = ";
              dumpIVS(&ivsIntersect);
              pout() << "============" << endl;
            }
        }

    }

  BaseIFFAB<Real>::setVerbose(true);
  //do the evil exchange
  Interval interv(0, nFlux-1);
  fluxInterpolant.exchange(interv);
  ibox = 0;
  //check the answer over grown set
  Real tolerance = 0.001;
  for (DataIterator dit = a_dbl.dataIterator(); dit.ok(); ++dit)
    {
      const BaseIFFAB<Real>& srcFab = fluxInterpolant[dit()];
      const IntVectSet& ivsGrown = irregSetsGrown[dit()];
      for (FaceIterator faceit(ivsGrown, a_ebisl[dit()].getEBGraph(), faceDir, FaceStop::SurroundingWithBoundary);
          faceit.ok(); ++faceit)
        {
          Real correct = rightAns(faceit());
          Real fabAns  =   srcFab(faceit(), 0);
          if (Abs(correct - fabAns)  > tolerance)
            {
              pout() << "iffab test failed at face "
                     << faceit().gridIndex(Side::Lo)
                     << faceit().gridIndex(Side::Hi) << endl;
              pout() << " right ans = " << correct  << endl;
              pout() << " data holds= " << fabAns << endl;

              int eekflag = -3;
              return eekflag;
            }

        }
      ibox++;
    }
  return 0;
}

int checkEBISL(const EBISLayout&        a_ebisl,
               const DisjointBoxLayout& a_grids,
               const Box&               a_domain,
               const Real&              a_dx,
               const RealVect&          a_origin,
               const Real&              a_radius,
               const RealVect&          a_center,
               const bool&              a_insideRegular)
{
  //check to see that the sum of the volume fractions
  //comes to close to exactly the total volume
  int eekflag =  0;
  //First: calculate volume of domain

  const IntVect& ivSize = a_domain.size();
  RealVect hiCorn;
  RealVect domLen;
  Real cellVolume = 1.0;
  Real totalVolume = 1.0;

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    hiCorn[idir] = a_origin[idir] + a_dx*Real(ivSize[idir]);
    cellVolume *= a_dx;
    domLen[idir] = hiCorn[idir] - a_origin[idir];
    totalVolume *= domLen[idir];
  }

  // Calculate the exact volume of the regular region
  Real volExact;
  Real volError;
  Real bndryExact;
  Real bndryError;

  if (SpaceDim == 2)
  {
    volExact = acos(-1.0) * a_radius*a_radius;
    volError = 0.0001851;

    bndryExact = 2*acos(-1.0) * a_radius;
    bndryError = 0.000047;
  }

  if (SpaceDim == 3)
  {
    volExact = 4.0/3.0 * acos(-1.0) * a_radius*a_radius*a_radius;
    volError = 0.000585;

    bndryExact = 4.0 * acos(-1.0) * a_radius*a_radius;
    bndryError = 0.000287;
  }

  if (a_insideRegular == false)
  {
    volExact = totalVolume - volExact;
  }

  // Now calculate the volume of approximate sphere
  Real tolerance = 0.001;
  Real volTotal = 0.0;
  Real bndryTotal = 0.0;


  for (DataIterator dit = a_grids.dataIterator(); dit.ok(); ++dit)
  {
    const EBISBox& ebisBox = a_ebisl[dit()];
    for (BoxIterator bit(a_grids.get(dit())); bit.ok(); ++bit)
    {
      const IntVect& iv = bit();
      Vector<VolIndex> vofs = ebisBox.getVoFs(iv);

      if (vofs.size() > 1)
      {
        eekflag = 2;
        return eekflag;
      }

      for (int ivof = 0; ivof < vofs.size(); ivof++)
      {
        Real volFrac = ebisBox.volFrac(vofs[ivof]);
        if (ebisBox.isIrregular(iv))
          {
            pout() << "iv = " << iv << ", volFrac = " << volFrac << endl;
          }
        volTotal += volFrac * cellVolume;

        Real bndryArea = ebisBox.bndryArea(vofs[ivof]);
        bndryTotal += bndryArea * cellVolume / a_dx;

        if (volFrac > tolerance && volFrac < 1.0 - tolerance)
        {
          if (!ebisBox.isIrregular(iv))
          {
            eekflag = 4;
            return eekflag;
          }
        }
      }
    }
  }
//.........这里部分代码省略.........

void
MappedLevelFluxRegister::define(const DisjointBoxLayout& a_dbl,
                          const DisjointBoxLayout& a_dblCoarse,
                          const ProblemDomain&     a_dProblem,
                          const IntVect&           a_nRefine,
                          int                      a_nComp,
                          bool                     a_scaleFineFluxes)
{
    CH_TIME("MappedLevelFluxRegister::define");
    m_isDefined = FluxRegDefined;  // Basically, define was called
    m_nRefine   = a_nRefine;

    m_scaleFineFluxes = a_scaleFineFluxes;

    DisjointBoxLayout coarsenedFine;
    coarsen(coarsenedFine, a_dbl, a_nRefine);

#ifndef DISABLE_TEMPORARY_FLUX_REGISTER_OPTIMIZATION
    // This doesn't work for multi-block calculations, which are
    // not properly nested. -JNJ

    //begin temporary optimization.  bvs
    int numPts = 0;
    for (LayoutIterator lit = a_dblCoarse.layoutIterator(); lit.ok(); ++lit) {
        numPts += a_dblCoarse[lit].numPts();
    }
    for (LayoutIterator lit = coarsenedFine.layoutIterator(); lit.ok(); ++lit) {
        numPts -= coarsenedFine[lit].numPts();
    }

    if (numPts == 0) {
        m_coarFlux.clear();
        // OK, fine region completely covers coarse region.  no registers.
        return;
    }
#endif

    //end temporary optimization.   bvs
    m_coarFlux.define( a_dblCoarse, a_nComp);
    m_isDefined |= FluxRegCoarseDefined;
    m_domain = a_dProblem;
    ProblemDomain coarsenedDomain;
    coarsen(coarsenedDomain, a_dProblem, a_nRefine);

    m_fineFlux.define( coarsenedFine, a_nComp, IntVect::Unit);
    m_isDefined |= FluxRegFineDefined;

    m_reverseCopier.ghostDefine(coarsenedFine, a_dblCoarse,
                                coarsenedDomain, IntVect::Unit);

    for (int i = 0; i < CH_SPACEDIM; i++) {
        m_coarseLocations[i].define(a_dblCoarse);
        m_coarseLocations[i + CH_SPACEDIM].define(a_dblCoarse);
    }


    DataIterator dC   = a_dblCoarse.dataIterator();
    LayoutIterator dF = coarsenedFine.layoutIterator();

    for (dC.begin(); dC.ok(); ++dC) {
        const Box& cBox = a_dblCoarse.get(dC);

        for (dF.begin(); dF.ok(); ++dF) {
            const Box& fBox  = coarsenedFine.get(dF);

            if (fBox.bigEnd(0) + 1 < cBox.smallEnd(0)) {
                //can skip this box since they cannot intersect, due to sorting
            } else if (fBox.smallEnd(0) - 1 > cBox.bigEnd(0)) {
                //skip to end, since all the rest of boxes will not intersect either
                dF.end();
            } else {

                for (int i = 0; i < CH_SPACEDIM; i++) {
                    Vector<Box>& lo = m_coarseLocations[i][dC];
                    Vector<Box>& hi = m_coarseLocations[i + CH_SPACEDIM][dC];

                    Box loBox = adjCellLo(fBox, i, 1);
                    Box hiBox = adjCellHi(fBox, i, 1);
                    if (cBox.intersectsNotEmpty(loBox)) lo.push_back(loBox & cBox);
                    if (cBox.intersectsNotEmpty(hiBox)) hi.push_back(hiBox & cBox);
                }
            }
        }
    }

    Box domainBox = coarsenedDomain.domainBox();

    if (a_dProblem.isPeriodic()) {
        Vector<Box> periodicBoxes[2 * CH_SPACEDIM];
        for (dF.begin(); dF.ok(); ++dF) {
            const Box& fBox  = coarsenedFine.get(dF);
            for (int i = 0; i < CH_SPACEDIM; i++) {
                if (a_dProblem.isPeriodic(i)) {
                    if (fBox.smallEnd(i) == domainBox.smallEnd(i))
                        periodicBoxes[i].push_back(adjCellLo(fBox, i, 1));
                    if (fBox.bigEnd(i) == domainBox.bigEnd(i))
                        periodicBoxes[i + CH_SPACEDIM].push_back(adjCellHi(fBox, i, 1));
                }
            }
        }
//.........这里部分代码省略.........

// ---------------------------------------------------------
void
NodeQCFI::define(const DisjointBoxLayout& a_grids,
                 Real a_dx,
                 const ProblemDomain& a_domain,
                 const LayoutData<NodeCFIVS>* const a_loCFIVS,
                 const LayoutData<NodeCFIVS>* const a_hiCFIVS,
                 int a_refToCoarse,
                 NodeBCFunc               a_bc,
                 int a_interpolationDegree,
                 int a_ncomp,
                 bool a_verbose)
{
  CH_assert(a_refToCoarse >= 2);
  m_bc = a_bc;
  m_grids = a_grids;
  m_dx = a_dx;
  m_verbose = a_verbose;

  // m_coarsenings == log2(a_refToCoarse);
  m_coarsenings = 0;
  for (int interRatio = a_refToCoarse; interRatio >= 2; interRatio /= 2)
    m_coarsenings++;

  m_qcfi2.resize(m_coarsenings);
  m_inter.resize(m_coarsenings-1);
  m_loCFIVScoarser.resize(m_coarsenings-1);
  m_hiCFIVScoarser.resize(m_coarsenings-1);
  // for i=1:m_c-2, m_qcfi2[i] interpolates m_inter[i] from m_inter[i-1].

  if (m_coarsenings == 1)
    {
      m_qcfi2[0] = new NodeQuadCFInterp2(a_grids, a_domain,
                                         a_loCFIVS, a_hiCFIVS,
                                         false, a_interpolationDegree, a_ncomp);
    }
  else
    {
      // if m_coarsenings == 2:
      // m_qcfi2[1] = new NodeQuadCFInterp2(a_grids, true);
      // m_qcfi2[0] = new NodeQuadCFInterp2(a_grids.coarsen(2), false);
      // m_inter[0] = new LevelData(a_grids.coarsen(2));

      ProblemDomain domainLevel(a_domain);
      Real dxLevel = m_dx;
      DisjointBoxLayout gridsInterFine(a_grids);
      DisjointBoxLayout gridsInterCoarse;
      // m_qcfi2[m_coarsenings-1] interpolates from a 2-coarsening of
      // the fine grids onto the fine grids (a_grids).
      // Interpolation is from the interface only.
      m_qcfi2[m_coarsenings-1] =
        new NodeQuadCFInterp2(a_grids, domainLevel,
                              a_loCFIVS, a_hiCFIVS,
                              true, a_interpolationDegree, a_ncomp);
      for (int interlev = m_coarsenings - 2; interlev >= 0; interlev--)
        {
          // m_coarsenings >= 2, so this will be executed at least once.
          coarsen(gridsInterCoarse, gridsInterFine, 2);
          domainLevel.coarsen(2);
          dxLevel *= 2.;
          bool interfaceOnly = (interlev > 0);
          // m_qcfi2[interlev] interpolates from some coarsening of
          // the fine grids onto gridsInterCoarse.
          // Except for the first time (interlev == 0),
          // interpolation is from the interface only.

          // Define objects containing the nodes of the coarse/fine interfaces.
          m_loCFIVScoarser[interlev] = new LayoutData<NodeCFIVS>[SpaceDim];
          m_hiCFIVScoarser[interlev] = new LayoutData<NodeCFIVS>[SpaceDim];
          for (int idir = 0; idir < SpaceDim; idir++)
            {
              LayoutData<NodeCFIVS>& loCFIVS =
                m_loCFIVScoarser[interlev][idir];
              LayoutData<NodeCFIVS>& hiCFIVS =
                m_hiCFIVScoarser[interlev][idir];

              loCFIVS.define(gridsInterCoarse);
              hiCFIVS.define(gridsInterCoarse);
              for (DataIterator dit = gridsInterCoarse.dataIterator(); dit.ok(); ++dit)
                {
                  const Box& bx = gridsInterCoarse.get(dit());
                  loCFIVS[dit()].define(domainLevel, bx, gridsInterCoarse,
                                        idir, Side::Lo);
                  hiCFIVS[dit()].define(domainLevel, bx, gridsInterCoarse,
                                        idir, Side::Hi);
                }
            }

          m_qcfi2[interlev] =
            new NodeQuadCFInterp2(gridsInterCoarse, domainLevel,
                                  m_loCFIVScoarser[interlev],
                                  m_hiCFIVScoarser[interlev],
                                  interfaceOnly, a_interpolationDegree, a_ncomp);

          m_inter[interlev] =
            new LevelData<NodeFArrayBox>(gridsInterCoarse, a_ncomp, IntVect::Zero);
          gridsInterFine = gridsInterCoarse;
        }
      m_domainPenultimate = domainLevel;
      m_dxPenultimate = dxLevel;
//.........这里部分代码省略.........

void
EBCoarseAverage::define(const DisjointBoxLayout& a_dblFine,
                        const DisjointBoxLayout& a_dblCoar,
                        const EBISLayout& a_ebislFine,
                        const EBISLayout& a_ebislCoar,
                        const ProblemDomain& a_domainCoar,
                        const int& a_nref,
                        const int& a_nvar,
                        const EBIndexSpace* ebisPtr)
{
  CH_TIME("EBCoarseAverage::define");
  CH_assert(ebisPtr->isDefined());

  ProblemDomain domainFine = a_domainCoar;
  domainFine.refine(a_nref);
  EBLevelGrid eblgFine;
  EBLevelGrid eblgCoar = EBLevelGrid(a_dblCoar, a_ebislCoar, a_domainCoar);
  EBLevelGrid eblgCoFi;

  //check to see if the input layout is coarsenable.
  //if so, proceed with ordinary drill
  //otherwise, see if the layout covers the domain.
  //if it does, we can use domainsplit
  if (a_dblFine.coarsenable(a_nref))
    {
      eblgFine = EBLevelGrid(a_dblFine, a_ebislFine,   domainFine);
      m_useFineBuffer = false;
    }
  else
    {
      Box fineDomBox = refine(a_domainCoar.domainBox(), a_nref);
      int numPtsDom = fineDomBox.numPts();
      //no need for gathers here because the meta data is global
      int numPtsLayout = 0;
      for (LayoutIterator lit = a_dblFine.layoutIterator(); lit.ok(); ++lit)
        {
          numPtsLayout += a_dblFine.get(lit()).numPts();
        }
      bool coveringDomain = (numPtsDom == numPtsLayout);
      if (coveringDomain)
        {
          m_useFineBuffer = true;
          int maxBoxSize = 4*a_nref;
          Vector<Box> boxes;
          Vector<int> procs;
          domainSplit(fineDomBox, boxes, maxBoxSize);
          mortonOrdering(boxes);
          LoadBalance(procs, boxes);
          DisjointBoxLayout dblBufFine(boxes, procs);

          eblgFine = EBLevelGrid(dblBufFine, domainFine, 2, eblgCoar.getEBIS());
        }
      else
        {
          pout() << "EBCoarseAverage::input layout is not coarsenable and does not cover the domain--bailing out" << endl;
          MayDay::Error();
        }
    }

  coarsen(eblgCoFi, eblgFine, a_nref);
  define(eblgFine, eblgCoar, eblgCoFi, a_nref, a_nvar);
}

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