C++ ProblemDomain::domainBox方法代码示例

NewPoissonOp* NewPoissonOpFactory::MGnewOp(const ProblemDomain& a_FineindexSpace,
                                           int   a_depth,
                                           bool  a_homoOnly)
{
  CH_assert(a_depth >= 0 );
  CH_assert(m_bc != NULL);
  NewPoissonOp* newOp = new NewPoissonOp();
  RealVect dx = m_dx;
  ProblemDomain domain = a_FineindexSpace;
  for (int i=0; i<a_depth; i++)
    {
      Box d = domain.domainBox();
      d.coarsen(8);
      d.refine(8);
      if (domain.domainBox() == d)
        {
          dx*=2;
          domain.coarsen(2);
        }
      else
        {
          return NULL;
        }
    }
  newOp->define(dx, domain, m_bc);
  return newOp;
}

// Define new AMR level
void AMRLevelPluto::define(AMRLevel*            a_coarserLevelPtr,
                           const ProblemDomain& a_problemDomain,
                           int                  a_level,
                           int                  a_refRatio)
{
  if (s_verbosity >= 3)
  {
    pout() << "AMRLevelPluto::define " << a_level << endl;
  }

  // Call inherited define
  AMRLevel::define(a_coarserLevelPtr,
                   a_problemDomain,
                   a_level,
                   a_refRatio);

  // Get setup information from the next coarser level
  if (a_coarserLevelPtr != NULL)
  {
    AMRLevelPluto* amrGodPtr = dynamic_cast<AMRLevelPluto*>(a_coarserLevelPtr);

    if (amrGodPtr != NULL)
    {
      m_cfl = amrGodPtr->m_cfl;
      m_domainLength = amrGodPtr->m_domainLength;
      m_refineThresh = amrGodPtr->m_refineThresh;
      m_tagBufferSize = amrGodPtr->m_tagBufferSize;
    }
    else
    {
      MayDay::Error("AMRLevelPluto::define: a_coarserLevelPtr is not castable to AMRLevelPluto*");
    }
  }

  // Compute the grid spacing
  m_dx = m_domainLength / (a_problemDomain.domainBox().bigEnd(0)-a_problemDomain.domainBox().smallEnd(0)+1.);

  // Nominally, one layer of ghost cells is maintained permanently and
  // individual computations may create local data with more
  m_numGhost = 1;

  CH_assert(m_patchPluto != NULL);
  CH_assert(isDefined());
  m_patchPluto->define(m_problem_domain,m_dx,m_level,m_numGhost);

  // Get additional information from the patch integrator
  m_numStates      = m_patchPluto->numConserved();
  m_ConsStateNames = m_patchPluto->ConsStateNames();
  m_PrimStateNames = m_patchPluto->PrimStateNames();
}

/** Initializes the DisjointBoxLayouts, each step. */
void
EBRestart::makeHierarchy(Vector<DisjointBoxLayout>& a_dbl,
                         const ProblemDomain&       a_baseDomain,
                         const IntVectSet&          a_baseTags,
                         const Vector<int>&         a_refRatios,
                         const InputParams&         a_inputs)
{
  a_dbl.resize(a_inputs.nlevs);
  Real fillRatio     = 0.85;
  BRMeshRefine regridder(a_baseDomain, a_refRatios, fillRatio,
                         a_inputs.blocking_factor,
                         a_inputs.buffer_size, a_inputs.max_size);

  Vector<Vector<Box> > oldGrids(a_inputs.nlevs,1), newGrids(a_inputs.nlevs);
  oldGrids[0][0]=a_baseDomain.domainBox();
  for ( int l=1; l<a_inputs.nlevs; ++l )
  {
    oldGrids[l][0]=coarsen(oldGrids[l-1][0], a_refRatios[l-1]);
    // FIXME: is a_refRatios[l-1] correct??
  }

  regridder.regrid(newGrids, a_baseTags, 0, 1, oldGrids);
  // 0 is baselevel, 1 is toplevel.

  Vector<int> procs;
  for (int l=0; l<a_inputs.nlevs; l++)
  {
    newGrids[l].sort();
    LoadBalance(procs, newGrids[l]);
    a_dbl[l] = DisjointBoxLayout(newGrids[l], procs);
  }
}

void CFStencil::buildPeriodicVector(Vector<Box>& a_periodicVector,
                                    const ProblemDomain& a_fineDomain,
                                    const DisjointBoxLayout& a_fineBoxes)
{
  Box periodicTestBox(a_fineDomain.domainBox());
  if (a_fineDomain.isPeriodic())
    {
      for (int idir=0; idir<SpaceDim; idir++)
        {
          if (a_fineDomain.isPeriodic(idir))
            {
              periodicTestBox.grow(idir,-1);
            }
        }
    }
  a_periodicVector.clear();
  a_periodicVector.reserve(a_fineBoxes.size());

  LayoutIterator lit = a_fineBoxes.layoutIterator();
  for (lit.reset(); lit.ok(); ++lit)
    {
      const Box& box = a_fineBoxes[lit()];
      a_periodicVector.push_back(box);
      // if periodic, also need to add periodic images
      // only do this IF we're periodic and box
      // adjacent to  the domain box boundary somewhere
      if (a_fineDomain.isPeriodic()
          && !periodicTestBox.contains(box))
        {
          ShiftIterator shiftIt = a_fineDomain.shiftIterator();
          IntVect shiftMult(a_fineDomain.domainBox().size());
          Box shiftedBox(box);
          for (shiftIt.begin(); shiftIt.ok(); ++shiftIt)
            {
              IntVect shiftVect = shiftMult*shiftIt();
              shiftedBox.shift(shiftVect);
              a_periodicVector.push_back(shiftedBox);
              shiftedBox.shift(-shiftVect);
            } // end loop over periodic shift directions
        } // end if periodic
    }
  a_periodicVector.sort();
}

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 DenseIntVectSet::nestingRegion(int radius, const ProblemDomain& a_domain)
{
  CH_assert(radius >= 0);
  if (radius == 0) return;

  DenseIntVectSet tmp;
  Box region = a_domain.domainBox();
  if (!a_domain.isPeriodic())
    {
      tmp = *this;
    }
  else
    {
      D_TERM6(if (a_domain.isPeriodic(0)) region.grow(0, radius);,
              if (a_domain.isPeriodic(1)) region.grow(1, radius);,
              if (a_domain.isPeriodic(2)) region.grow(2, radius);,

int VCAMRPoissonOp2Factory::refToFiner(const ProblemDomain& a_domain) const
{
  int retval = -1;
  bool found = false;

  for (int ilev = 0; ilev < m_domains.size(); ilev++)
    {
      if (m_domains[ilev].domainBox() == a_domain.domainBox())
        {
          retval = m_refRatios[ilev];
          found = true;
        }
    }

  if (!found)
    {
      MayDay::Abort("Domain not found in AMR hierarchy");
    }

  return retval;
}

void
CFStencil::define(
                  const ProblemDomain& a_fineDomain,
                  const Box& a_grid,
                  const DisjointBoxLayout& a_fineBoxes,
                  const DisjointBoxLayout& a_coarBoxes,
                  int a_refRatio,
                  int a_direction,
                  Side::LoHiSide a_hiorlo)
{
  m_isDefined = true;
  CH_assert(a_refRatio >= 1);
  CH_assert(a_direction >= 0);
  CH_assert(a_direction < SpaceDim);
  CH_assert((a_hiorlo == Side::Lo) ||
         (a_hiorlo == Side::Hi));
  CH_assert(!a_fineDomain.isEmpty());

  //set internal vars.  most of these are kept around
  //just to keep the class from having an identity crisis.
  m_direction = a_direction;
  m_hiorlo =  a_hiorlo;

  Box finebox = a_grid;


  //compute intvectset of all points on fine grid that
  //need to be interpolated

  //shift direction
  int hilo = sign(a_hiorlo);

  //create fine stencil
  Box edgebox;
  CH_assert((hilo ==1) || (hilo == -1));
  if (hilo == -1)
    {
      edgebox = adjCellLo(finebox,m_direction,1);
    }
  else
    {
      edgebox = adjCellHi(finebox,m_direction,1);
    }
  edgebox = a_fineDomain & edgebox;
  if (!edgebox.isEmpty())
    {
      Box periodicTestBox(a_fineDomain.domainBox());
      if (a_fineDomain.isPeriodic())
        {
          for (int idir=0; idir<SpaceDim; idir++)
            {
              if (a_fineDomain.isPeriodic(idir))
                {
                  periodicTestBox.grow(idir,-1);
                }
            }
        }

      m_fineIVS.define(edgebox);
      LayoutIterator lit = a_fineBoxes.layoutIterator();
      for (lit.reset(); lit.ok(); ++lit)
        {
          m_fineIVS -= a_fineBoxes[lit()];
          // if periodic, also need to subtract periodic images
          // only do this IF we're periodic _and_ both boxes
          // adjoin the domain box boundary somewhere
          if (a_fineDomain.isPeriodic() && !periodicTestBox.contains(edgebox)
              && !periodicTestBox.contains(a_fineBoxes[lit()]))
            {
              ShiftIterator shiftIt = a_fineDomain.shiftIterator();
              IntVect shiftMult(a_fineDomain.domainBox().size());
              Box shiftedBox(a_fineBoxes[lit()]);
              for (shiftIt.begin(); shiftIt.ok(); ++shiftIt)
                {
                  IntVect shiftVect = shiftMult*shiftIt();
                  shiftedBox.shift(shiftVect);
                  m_fineIVS -= shiftedBox;
                  shiftedBox.shift(-shiftVect);
                } // end loop over periodic shift directions
            } // end if periodic
        }
    }

  //ivs where all coarse slopes are defined
  //== coarsened fine ivs
  m_coarIVS.define(m_fineIVS);
  m_coarIVS.coarsen(a_refRatio);
  // this is a trick to get around the lack of a IntVectSet intersection
  // operator which works with a ProblemDomain
  ProblemDomain coardom= coarsen(a_fineDomain, a_refRatio);
  Box domainIntersectBox = m_coarIVS.minBox();
  domainIntersectBox = coardom & domainIntersectBox;
  m_coarIVS &= domainIntersectBox;

  m_packedBox = m_fineIVS.minBox();
  if (m_fineIVS.numPts() == m_packedBox.numPts())
    {
      m_isPacked = true;
    }
  else
//.........这里部分代码省略.........

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 ReductionCopier::define(const BoxLayout& a_level,
                             const BoxLayout& a_dest,
                             const ProblemDomain& a_domain,
                             const IntVect& a_ghost,
                             const Vector<int>& a_transverseDir,
                             bool  a_exchange)
{
#ifdef MULTIDIM_TIMER
  CH_TIME("ReductionCopier::define")
#endif
  m_isDefined = true;
  m_transverseDir = a_transverseDir;

  CH_assert(a_level.isClosed());
  CH_assert(a_dest.isClosed());
  //  CH_assert(a_level.checkPeriodic(a_domain));

  clear();
  buffersAllocated = false;
  //bool self = a_dest == a_level;
  const BoxLayout&         level= a_level;
  const BoxLayout&         dest = a_dest;

  // note that while much of the implementation for the
  // periodic case remains in this function (as copied, more
  // or less, from the original Copier), at least at the moment
  // periodic support is "turned off" for the ReductionCopier
  // because it's not entirely clear it will do the right thing
  // for this case.  I've left the vast majority of the
  // periodic implementation in place, however, in case we
  // change our minds... (DFM 3/10/09)

  // set up vector of dataIndexes to keep track of which
  // "to" boxes are not completely contained within the primary
  // domain.  these boxes are then candidates for filling by
  // periodic images of the "from" data.
  Vector<DataIndex> periodicallyFilledToVect;

  // in order to cull which "from" data may be needed to
  // fill the "to" data, keep track of the radius around the
  // primary domain in which all these cells lie.
  // do this by incrementally growing the domain box and
  // keeping track of what this radius is.
  // just to make things simpler, start off with a radius of one
  Box grownDomainCheckBox = a_domain.domainBox();
  grownDomainCheckBox.grow(1);
  int periodicCheckRadius = 1;

  // since valid regions of the "from" DBL may also be outside
  // the primary domain, need to keep track of whether any of these
  // need to be checked separately.
  Vector<DataIndex> periodicFromVect;
  // use same domain trick here as well
  Box grownFromDomainCheckBox = a_domain.domainBox();
  int periodicFromCheckRadius = 1;

  Box domainBox(a_domain.domainBox());
  // grab index extents of domain in transverse direction
  Vector<int> transverseLo(m_transverseDir.size());
  Vector<int> transverseHi(m_transverseDir.size());

  for (int n=0; n<m_transverseDir.size(); n++)
    {
      // grab index extents of domain in transverse direction
      transverseLo[n] = domainBox.smallEnd(m_transverseDir[n]);
      transverseHi[n] = domainBox.bigEnd(m_transverseDir[n]);
    }

  bool isPeriodic = false;
  // I think the right thing here is that the ReductionCopier
  // completely ignores periodicity
  //if (!domainBox.isEmpty())
  //isPeriodic = a_domain.isPeriodic();

  // (dfm -- 9/13/05) as currently written, the Copier won't correctly
  // handle periodic cases where the number of ghost cells is greater
  // than the width of the domain.  We _should_ do multiple wraparounds,
  // but we don't. So, put in this assertion. We can revisit this if it
  // becomes an issue
  if (isPeriodic)
    {
      for (int dir=0; dir<SpaceDim; dir++)
        {
          if (a_domain.isPeriodic(dir))
            {
              CH_assert (a_ghost[dir] <= domainBox.size(dir));
            }
        }
    }

  unsigned int myprocID = procID();

  // The following 4 for loops are the result of a performance optimization.
  // When increasing the size of the problem, we found that the code was
  // looping over every destination box for every source box which was N1*N2
  // loop iterations (essentially an N-squared approach).
  // The following code attempts to simply reduce N1 and N2 by first separating
  // the boxes (or LayoutIndexes to boxes) that reside on the current processor.
  // Then the loop to determine which boxes of the first list intersect with
  // which boxes of the second list can be done in N1' * N2' iterations,
//.........这里部分代码省略.........

void ConstBCFunction::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])
        {
          int  bcType;
          Real bcValue;

          if (side == Side::Lo)
          {
            bcType  = m_loSideType [idir];
            bcValue = m_loSideValue[idir];
          }
          else
          {
            bcType  = m_hiSideType [idir];
            bcValue = m_hiSideValue[idir];
          }

          if (bcType == 0)
          {
            // Neumann BC
            int isign = sign(side);

            Box toRegion = adjCellBox(a_valid, idir, side, 1);
            toRegion &= a_state.box();

            Box fromRegion = toRegion;
            fromRegion.shift(idir, -isign);

            a_state.copy(a_state, fromRegion, 0, toRegion, 0, a_state.nComp());

            if (!a_homogeneous)
            {
              for (BoxIterator bit(toRegion); bit.ok(); ++bit)
              {
                const IntVect& ivTo = bit();
                // IntVect ivClose = ivTo - isign*BASISV(idir);

                for (int icomp = 0; icomp < a_state.nComp(); icomp++)
                {
                  a_state(ivTo, icomp) += Real(isign)*a_dx*bcValue;
                }
              }
            }
          }
          else if (bcType == 1)
          {
            // Dirichlet BC
            int isign = sign(side);

            Box toRegion = adjCellBox(a_valid, idir, side, 1);
            toRegion &= a_state.box();

            for (BoxIterator bit(toRegion); bit.ok(); ++bit)
            {
              const IntVect& ivTo = bit();

              IntVect ivClose = ivTo -   isign*BASISV(idir);
              // IntVect ivFar   = ivTo - 2*isign*BASISV(idir);

              Real inhomogVal = 0.0;

              if (!a_homogeneous)
              {
                inhomogVal = bcValue;
              }

              for (int icomp = 0; icomp < a_state.nComp(); icomp++)
              {
                Real nearVal = a_state(ivClose, icomp);
                // Real farVal  = a_state(ivFar,   icomp);

                Real ghostVal = linearInterp(inhomogVal, nearVal);

                a_state(ivTo, icomp) = ghostVal;
              }
            }
          }
          else
          {
            MayDay::Abort("ConstBCFunction::operator() - unknown BC type");
          }
        } // if ends match
      } // end loop over sides
    } // if not periodic in this direction
  } // end loop over directions
//.........这里部分代码省略.........

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

// ---------------------------------------------------------
// 7 Dec 2005
Real
norm(const BoxLayoutData<NodeFArrayBox>& a_layout,
     const LevelData<NodeFArrayBox>& a_mask,
     const ProblemDomain& a_domain,
     const Real a_dx,
     const int a_p,
     const Interval& a_interval,
     bool a_verbose)
{
  if (a_p == 0)
    return maxnorm(a_layout, a_mask, a_domain, a_interval, a_verbose);

  Real normTotal = 0.;
  int ncomp = a_interval.size();
  Box domBox = a_domain.domainBox();
  for (DataIterator it = a_layout.dataIterator(); it.ok(); ++it)
    {
      const NodeFArrayBox& thisNfab = a_layout[it()];
      const FArrayBox& dataFab = thisNfab.getFab();
      const FArrayBox& maskFab = a_mask[it()].getFab();
      const Box& thisBox(a_layout.box(it())); // CELL-centered
      NodeFArrayBox dataMasked(thisBox, ncomp);
      FArrayBox& dataMaskedFab = dataMasked.getFab();
      dataMaskedFab.copy(dataFab);
      // dataMaskedFab *= maskFab;
      for (int comp = a_interval.begin(); comp <= a_interval.end(); comp++)
        {
          // Set dataMaskedFab[comp] *= maskFab[0].
          dataMaskedFab.mult(maskFab, 0, comp);
        }

      Real thisNfabNorm = 0.;
      if (thisBox.intersects(domBox))
        {
          Box thisBoxInDomain = thisBox & domBox;
          thisNfabNorm = norm(dataMasked, a_dx, thisBoxInDomain, a_p,
                              a_interval.begin(), a_interval.size());
        }

      if (a_verbose)
        cout << a_p << "norm(" << thisBox << ") = " << thisNfabNorm << endl;

      if (a_p == 1)
        {
          normTotal += thisNfabNorm;
        }
      else if (a_p == 2)
        {
          normTotal += thisNfabNorm * thisNfabNorm;
        }
      else
        {
          normTotal += pow(thisNfabNorm, Real(a_p));
        }
    }
# ifdef CH_MPI
  Real recv;
  // add up (a_p is not 0)
  int result = MPI_Allreduce(&normTotal, &recv, 1, MPI_CH_REAL,
                             MPI_SUM, Chombo_MPI::comm);
  if (result != MPI_SUCCESS)
    { //bark!!!
      MayDay::Error("sorry, but I had a communication error on norm");
    }
  normTotal = recv;
# endif
  // now do sqrt, etc
  if (a_p == 2)
    normTotal = sqrt(normTotal);
  else
    if ((a_p != 0) && (a_p != 1))
      normTotal = pow(normTotal, (Real)1.0/Real(a_p));

  return normTotal;
}