C++ DisjointBoxLayout类代码示例

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

    }

}

int
makeLayout(DisjointBoxLayout& a_dbl,
           const Box& a_domain)
{
    ParmParse pp;
    int eekflag= 0;
    int ipieces;
    ipieces = Max(ipieces, 1);
    int maxsize;
    pp.get("maxboxsize",maxsize);
    Vector<Box> vbox(1, a_domain);
    domainSplit(a_domain, vbox,  maxsize);
    if (eekflag != 0)
    {
        pout() << "problem in domainsplit" << endl;
        return eekflag;
    }
    Vector<int>  procAssign;
    eekflag = LoadBalance(procAssign,vbox);
    if (eekflag != 0)
    {
        pout() << "problem in loadbalance" << endl;
        return eekflag;
    }
    a_dbl.define(vbox, procAssign);
    return eekflag;
}

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
makeLayout(DisjointBoxLayout& a_dbl,
           const Box& a_domainFine)
{
  //set up mesh refine object
  ParmParse pp;
  int eekflag= 0;
  int maxsize;
  pp.get("maxboxsize",maxsize);
  int bufferSize = 1;
  int blockFactor = 2;
  Real fillrat = 0.75;
  Box domainCoar = coarsen(a_domainFine, 2);
  Vector<int> refRat(2,2);
  BRMeshRefine meshRefObj(domainCoar, refRat, fillrat,
                          blockFactor, bufferSize, maxsize);

  Vector<Vector<Box> > oldMeshes(2);
  oldMeshes[0] = Vector<Box>(1,   domainCoar);
  oldMeshes[1] = Vector<Box>(1, a_domainFine);

  //set up coarse tags
  int nc = domainCoar.size(0);
  int nmi = nc/2;//16
  int nqu = nc/4;//8
  int ntf = (nc*3)/4;  //24
  int nte = (nc*3)/8; //12
  int nfe = (nc*5)/8; //20
#if (CH_SPACEDIM ==2)
  Box boxf1(IntVect(0, nqu), IntVect(nmi-1,ntf-1));
  Box boxf2(IntVect(nmi,nte), IntVect(ntf-1,nfe-1));
  Box boxf3(IntVect(nqu,0  ), IntVect(nfe-1,nqu-1));
  Box boxf4(IntVect(nfe,nqu), IntVect(nc -1,nte-1));
#else
  Box boxf1(IntVect(0, nqu,nqu), IntVect(nmi-1,ntf-1,ntf-1));
  Box boxf2(IntVect(nmi,nte,nte), IntVect(ntf-1,nfe-1,nfe-1));
  Box boxf3(IntVect(nqu,0,0  ), IntVect(nfe-1,nqu-1,nqu-1));
  Box boxf4(IntVect(nfe,nqu,nqu), IntVect(nc -1,nte-1,nte-1));
#endif
  IntVectSet tags;
  tags |= boxf1;
  tags |= boxf2;
  tags |= boxf3;
  tags |= boxf4;

  int baseLevel = 0;
  int topLevel = 0;
  Vector<Vector<Box> > newMeshes;
  meshRefObj.regrid(newMeshes, tags, baseLevel,
                    topLevel, oldMeshes);

  const Vector<Box>& vbox = newMeshes[1];
  Vector<int>  procAssign;
  eekflag = LoadBalance(procAssign,vbox);
  if (eekflag != 0) return eekflag;
  a_dbl.define(vbox, procAssign);
  return eekflag;
}

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

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

int
makeLayoutPart(DisjointBoxLayout& a_dbl,
               const Box& a_domainCoarsest)
{
  //  Vector<int> vecRefRat(2, 2);
  ParmParse pp;
  int eekflag= 0;
  int blockFactor, bufferSize, maxSize;
  Real fillRat;
  pp.get("block_factor", blockFactor);
  pp.get("buffer_size", bufferSize);
  pp.get("maxboxsize", maxSize);
  pp.get("fill_ratio", fillRat);
  int nlevels =  2;
  Vector<int> vecRefRat(nlevels);
  pp.getarr("ref_ratio", vecRefRat,0,nlevels);
  BRMeshRefine mesher(a_domainCoarsest, vecRefRat,
                      fillRat, blockFactor, bufferSize, maxSize);

  int topLevel  =  0;
  int baseLevel =  0;
  //tags at base level
  IntVectSet tags;
  eekflag = makeTags(tags, a_domainCoarsest);
  if (eekflag < 0) return eekflag;

  Vector<Vector<Box> > oldMeshes(nlevels);
  oldMeshes[0] = Vector<Box>(1, a_domainCoarsest);
  Box finerDomain = a_domainCoarsest;
  for (int ilev = 1; ilev < nlevels; ilev++)
    {
      finerDomain.refine(vecRefRat[ilev]);
      oldMeshes[ilev] = Vector<Box>(1, finerDomain);
    }
  Vector<Vector<Box> > newMeshes;
  mesher.regrid(newMeshes, tags, baseLevel, topLevel, oldMeshes);

  Vector<int> procAssign;
  eekflag = LoadBalance(procAssign, newMeshes[1]);
  if (eekflag != 0) return eekflag;
  a_dbl.define(newMeshes[1], procAssign);
  int iverbose;
  pp.get("verbose", iverbose);
  if (iverbose == 1)
    {
      pout() << "the grids coarser domain= " << a_domainCoarsest
             << " = " << a_dbl << endl;
    }

  return eekflag;
}

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
NoFlowVortex::
initializeScalar ( LevelData<EBCellFAB>&    a_scalar,
                   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;
          Real scal;
          getXVal(xval, a_origin,vofit(), a_dx);
          getScalarPt(scal, xval);
          a_scalar[dit()](vofit(), 0) = scal;
        }
    }
}

void
NoFlowVortex::
initializePressure(LevelData<EBCellFAB>&    a_pressure,
                   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, pt;
          getXVal(xval, a_origin,vofit(), a_dx);
          for (int idir = 0; idir < SpaceDim; idir++)
            {
              a_pressure[dit()](vofit(), idir) = 1.e99;
            }
        }
    }
}

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