C++ EBISBox::getEBGraph方法代码示例

void
setToExactFlux(EBFluxFAB&             a_flux,
               const EBISBox&         a_ebisBox,
               const Box&             a_region,
               const Real&            a_dx)
{
    IntVectSet ivsregion(a_region);
    for (int faceDir = 0; faceDir < SpaceDim; faceDir++)
    {
        for (FaceIterator faceit(ivsregion, a_ebisBox.getEBGraph(), faceDir, FaceStop::SurroundingWithBoundary); faceit.ok(); ++faceit)
        {
            RealVect xval;
            IntVect iv = faceit().gridIndex(Side::Hi);
            RealVect centroid = a_ebisBox.centroid(faceit());
            for (int idir = 0; idir < SpaceDim; idir++)
            {
                if (idir == faceDir)
                {
                    xval[idir] = (Real(iv[idir]))*a_dx;
                }
                else
                {
                    xval[idir] = (Real(iv[idir]) + 0.5 + centroid[idir])*a_dx;
                }
            }
            Real fluxDir = exactFlux(xval, faceDir);
            a_flux[faceDir](faceit(), 0) = fluxDir;
        }
    } //end loop over face directions

}

bool
EBFluxRegister::
copyBIFFToEBFF(EBFaceFAB&             a_dst,
               const BaseIFFAB<Real>& a_src,
               const Box            & a_box,
               const EBISBox&         a_ebisBox)
{
  IntVectSet ivs = a_src.getIVS();
  ivs &= a_box;

  bool hasCells = !(ivs.isEmpty());
  if (hasCells)
    {
      a_dst.define(a_ebisBox, a_box, a_src.direction(), a_src.nComp());
      a_dst.setVal(0.);
      for (FaceIterator faceit(ivs, a_ebisBox.getEBGraph(),
                              a_src.direction(),
                              FaceStop::SurroundingWithBoundary);
          faceit.ok(); ++faceit)
        {
          for (int icomp = 0; icomp < a_src.nComp(); icomp++)
            {
              a_dst(faceit(), icomp) = a_src(faceit(), icomp);
            }
        }
    }
  return (hasCells);
}

int
testEBFluxFAB(const EBISBox& a_ebisBox, const Box& a_box)
{
  Interval comps(0,0);
  IntVectSet ivs(a_box);
  EBFluxFAB srcFab( a_ebisBox, a_box, 1);
  EBFluxFAB dstFab( a_ebisBox, a_box, 1);
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      //set source fab to right ans
      for (FaceIterator faceit(ivs, a_ebisBox.getEBGraph(), idir, FaceStop::SurroundingWithBoundary);
          faceit.ok(); ++faceit)
        {
          srcFab[idir](faceit(), 0) = rightAns(faceit());
        }
    }

  //linearize the data to dst
  int sizeFab = srcFab.size(a_box, comps);
  unsigned char* buf = new unsigned char[sizeFab];
  srcFab.linearOut(buf, a_box, comps);
  dstFab.linearIn( buf, a_box, comps);
  delete[] buf;

  //check the answer
  int eekflag = 0;
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      Real tolerance = 0.001;
      for (FaceIterator faceit(ivs, a_ebisBox.getEBGraph(), idir, FaceStop::SurroundingWithBoundary);
          faceit.ok(); ++faceit)
        {
          Real correct  = rightAns(faceit());
          if (Abs(dstFab[idir](faceit(), 0) - correct) > tolerance)
            {
              pout() << "ivfab test failed at face "
                     << faceit().gridIndex(Side::Lo)
                     << faceit().gridIndex(Side::Hi) << endl;

              eekflag = -4;
              return eekflag;
            }
        }
    }
  return 0;
}

void
NoFlowAdvectBC::
fluxBC(EBFluxFAB&            a_primGdnv,
       const EBCellFAB&      a_primCenter,
       const EBCellFAB&      a_primExtrap,
       const Side::LoHiSide& a_side,
       const Real&           a_time,
       const EBISBox&        a_ebisBox,
       const DataIndex&      a_dit,
       const Box&            a_box,
       const Box&            a_faceBox,
       const int&            a_dir)
{
  CH_assert(m_isDefined);

  Box FBox = a_faceBox;
  Box cellBox = FBox;
  CH_assert(a_primGdnv[a_dir].nComp()==1);

  // Determine which side and thus shifting directions
  int isign = sign(a_side);
  cellBox.shiftHalf(a_dir,isign);

  // Is there a domain boundary next to this grid
  if (!m_domain.contains(cellBox))
    {
      cellBox &= m_domain;
      // Find the strip of cells next to the domain boundary
      Box bndryBox = adjCellBox(cellBox, a_dir, a_side, 1);

      // Shift things to all line up correctly
      bndryBox.shift(a_dir,-isign);

      IntVectSet ivs(bndryBox);
      for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();

          Vector<FaceIndex> bndryFaces = a_ebisBox.getFaces(vof, a_dir, a_side);
          for (int iface= 0; iface < bndryFaces.size(); iface++)
            {
              const FaceIndex& face = bndryFaces[iface];
              //set all fluxes to zero then fix momentum flux
              //solid wall
              if (a_dir == m_velComp)
                {
                  a_primGdnv[a_dir](face, 0) = 0.0;
                }
              else
                {
                  a_primGdnv[a_dir](face, 0) = a_primExtrap(vof, 0);
                }
            }
        }
    }
}

void divergence(EBCellFAB&                a_divF,
                const EBFluxFAB&          a_flux,
                const EBISBox&            a_ebisBox,
                const Box&                a_box,
                const RealVect&           a_fluxVal,
                const Real&               a_dx)
{
  IntVectSet ivs(a_box);
  for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
    {
      Real divval = divergence(a_flux, a_fluxVal, a_ebisBox, vofit(), a_dx);
      a_divF(vofit(), 0) = divval;
    }

}

void
EBPlanarShockIBC::
setBndrySlopes(EBCellFAB&       a_deltaPrim,
               const EBCellFAB& a_primState,
               const EBISBox&   a_ebisBox,
               const Box&       a_box,
               const int&       a_dir)
{
  // don't want to deal with the periodic case
 CH_assert(!m_domain.isPeriodic(a_dir));
 CH_assert(m_isDefined);
 CH_assert(m_isFortranCommonSet);

  Box loBox,hiBox,centerBox,domain;
  int hasLo,hasHi;
  Box slopeBox = a_deltaPrim.getRegion()&m_domain;
  int numSlop = a_deltaPrim.nComp();
  // Generate the domain boundary boxes, loBox and hiBox, if there are
  // domain boundarys there
  eblohicenter(loBox,hasLo,hiBox,hasHi,centerBox,domain,
             slopeBox,m_domain,a_dir);

  // Set the boundary slopes if necessary
  if ((hasLo != 0) || (hasHi != 0))
    {
      BaseFab<Real>& regDeltaPrim = a_deltaPrim.getSingleValuedFAB();
      const BaseFab<Real>& regPrimState = a_primState.getSingleValuedFAB();
      /**/
      FORT_SLOPEBCS(CHF_FRA(regDeltaPrim),
                    CHF_CONST_FRA(regPrimState),
                    CHF_CONST_REAL(m_dx),
                    CHF_CONST_INT(a_dir),
                    CHF_BOX(loBox),
                    CHF_CONST_INT(hasLo),
                    CHF_BOX(hiBox),
                    CHF_CONST_INT(hasHi));
      /**/
    }
  for(SideIterator sit; sit.ok(); ++sit)
    {
      bool doThisSide;
      Box thisBox;
      if(sit() == Side::Lo)
        {
          doThisSide = (hasLo != 0);
          thisBox = loBox;
        }
      else
        {
          doThisSide = (hasHi != 0);
          thisBox = hiBox;
        }
      if(doThisSide)
        {

          // the cells for which the regular calculation
          //are incorrect are the grown set of the multivalued
          //cells intersected with the appropriate bndry box
          //and added to the irregular cells on the boundary.
          Box boxGrown = thisBox;
          boxGrown.grow(a_dir, 1);
          boxGrown &= m_domain;
          IntVectSet ivs = a_ebisBox.getMultiCells(boxGrown);
          ivs &= loBox;
          IntVectSet ivsIrreg = a_ebisBox.getIrregIVS(thisBox);
          ivs |= ivsIrreg;
          for(VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
            {
              const VolIndex& vof = vofit();
              //all slopes at boundary get set to zero
              //except normal velocity.  just following the
              //fortran here
              int inormVelVar = a_dir + QVELX;
              for(int ivar = 0; ivar < numSlop; ivar++)
                {
                  if(ivar != inormVelVar)
                    {
                      a_deltaPrim(vof, ivar) = 0.0;
                    }
                }

              //for normal velocity
              //do strangely limited slope
              //just lifted from the fortran.
              Side::LoHiSide otherSide = flip(sit());
              Vector<FaceIndex> faces =
                a_ebisBox.getFaces(vof, a_dir, otherSide);

              Real thisVel = a_primState(vof, inormVelVar);
              Real otherVel = 0.;
              for(int iface = 0; iface < faces.size(); iface++)
                {
                  VolIndex otherVoF = faces[iface].getVoF(otherSide);
                  otherVel += a_primState(otherVoF,inormVelVar);
                }
              if(faces.size() > 0)
                {
                  otherVel /= Real(faces.size());
                  Real slope;
                  if(sit() == Side::Lo)
//.........这里部分代码省略.........

void
EBPlanarShockIBC::
fluxBC(EBFluxFAB&            a_flux,
       const EBCellFAB&      a_primCenter,
       const EBCellFAB&      a_primExtrap,
       const Side::LoHiSide& a_side,
       const Real&           a_time,
       const EBISBox&        a_ebisBox,
       const DataIndex&      a_dit,
       const Box&            a_box,
       const Box&            a_faceBox,
       const int&            a_dir)
{
 CH_assert(m_isDefined);
 CH_assert(m_isFortranCommonSet);
 CH_assert(!m_domain.isPeriodic(a_dir));

  Box FBox = a_flux[a_dir].getSingleValuedFAB().box();
  Box cellBox = FBox;
  int numFlux = a_flux[a_dir].nComp();

  // Determine which side and thus shifting directions
  int isign = sign(a_side);
  cellBox.shiftHalf(a_dir,isign);

  // Is there a domain boundary next to this grid
  if (!m_domain.contains(cellBox))
    {
      cellBox &= m_domain;
      // Find the strip of cells next to the domain boundary
      Box boundaryBox = bdryBox(cellBox, a_dir, a_side, 1);

      // Shift things to all line up correctly
      boundaryBox.shiftHalf(a_dir,-isign);
      BaseFab<Real>& regFlux = a_flux[a_dir].getSingleValuedFAB();
      const BaseFab<Real>& regPrimExtrap = a_primExtrap.getSingleValuedFAB();

      // Set the boundary fluxes
      bool inbackofshock =
        (!m_shockbackward && a_side == Side::Lo) ||
        ( m_shockbackward && a_side == Side::Hi);
      if(a_dir == m_shocknorm && inbackofshock)
        {
          regFlux.shiftHalf(a_dir,-isign);

          // Set the boundary fluxes
          /**/
          FORT_EXTRAPBC(CHF_FRA(regFlux),
                        CHF_CONST_FRA(regPrimExtrap),
                        CHF_CONST_REAL(m_dx),
                        CHF_CONST_INT(a_dir),
                        CHF_BOX(boundaryBox));
          /**/

          // Shift returned fluxes to be face centered
          regFlux.shiftHalf(a_dir,isign);
          //now for the multivalued cells.  Since it is pointwise,
          //the regular calc is correct for all single-valued cells.
          IntVectSet ivs = a_ebisBox.getMultiCells(boundaryBox);
          for(VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
            {
              const VolIndex& vof = vofit();
              Vector<Real> qgdnv(QNUM);
              Vector<Real> fluxv(FNUM);
              for(int ivar = 0; ivar < QNUM; ivar++)
                {
                  qgdnv[ivar] = a_primExtrap(vof, ivar);
                }
              Vector<FaceIndex> bndryFaces =
                a_ebisBox.getFaces(vof, a_dir, a_side);
              for(int iface= 0; iface < bndryFaces.size(); iface++)
                {
                  /**/
                  FORT_POINTGETFLUX(CHF_VR(fluxv),
                                    CHF_VR(qgdnv),
                                    CHF_CONST_INT(a_dir));
                  /**/
                  const FaceIndex& face = bndryFaces[iface];
                  for(int ivar = 0; ivar < FNUM; ivar++)
                    {
                      a_flux[a_dir](face, ivar) = fluxv[ivar];
                    }
                }
            }
        } //end if on the back side of the shock
      else
        {
          regFlux.shiftHalf(a_dir,-isign);
          //solid wall bcs
          FORT_SLIPWALLSOLIDBC(CHF_FRA(regFlux),
                               CHF_CONST_FRA(regPrimExtrap),
                               CHF_CONST_INT(isign),
                               CHF_CONST_REAL(m_dx),
                               CHF_CONST_INT(a_dir),
                               CHF_BOX(boundaryBox));

          // Shift returned fluxes to be face centered
          regFlux.shiftHalf(a_dir,isign);
          int inormMomVar = CMOMX + a_dir;
          int inormVelVar = QVELX + a_dir;
//.........这里部分代码省略.........

void
kappaDivergence(EBCellFAB&             a_divF,
                const EBFluxFAB&       a_flux,
                const EBISBox&         a_ebisBox,
                const Box&             a_box,
                const Real&            a_dx)
{
    //set the divergence initially to zero
    //then loop through directions and increment the divergence
    //with each directions flux difference.
    a_divF.setVal(0.0);
    BaseFab<Real>&       regDivF = a_divF.getSingleValuedFAB();
    regDivF.setVal(0.);
    for (int idir = 0; idir < SpaceDim; idir++)
    {
        //update for the regular vofs in the nonconservative
        //case  works for all single valued vofs.
        /* do the regular vofs */
        /**/
        const EBFaceFAB& fluxDir = a_flux[idir];
        const BaseFab<Real>& regFluxDir = fluxDir.getSingleValuedFAB();
        int ncons = 1;
        FORT_DIVERGEF( CHF_BOX(a_box),
                       CHF_FRA(regDivF),
                       CHF_CONST_FRA(regFluxDir),
                       CHF_CONST_INT(idir),
                       CHF_CONST_INT(ncons),
                       CHF_CONST_REAL(a_dx));
        /**/
    }
    //update the irregular vofs using conservative diff
    IntVectSet ivsIrreg = a_ebisBox.getIrregIVS(a_box);
    for (VoFIterator vofit(ivsIrreg, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
    {
        const VolIndex& vof = vofit();
        //divergence was set in regular update.  we reset it
        // to zero and recalc.
        Real update = 0.;
        for ( int idir = 0; idir < SpaceDim; idir++)
        {
            const EBFaceFAB& fluxDir = a_flux[idir];
            for (SideIterator sit; sit.ok(); ++sit)
            {
                int isign = sign(sit());
                Vector<FaceIndex> faces =
                    a_ebisBox.getFaces(vof, idir, sit());
                for (int iface = 0; iface < faces.size(); iface++)
                {
                    const FaceIndex& face = faces[iface];
                    Real areaFrac = a_ebisBox.areaFrac(face);
                    Real faceFlux =fluxDir(face, 0);
                    update += isign*areaFrac*faceFlux;

                }
            }
        }
        //add EB boundary condtions in divergence
        const IntVect& iv = vof.gridIndex();
        Real bndryArea = a_ebisBox.bndryArea(vof);
        RealVect bndryCent = a_ebisBox.bndryCentroid(vof);
        RealVect normal = a_ebisBox.normal(vof);
        RealVect bndryLoc;
        RealVect exactF;
        for (int idir = 0; idir < SpaceDim; idir++)
        {
            bndryLoc[idir] = a_dx*(iv[idir] + 0.5 + bndryCent[idir]);
        }
        for (int idir = 0; idir < SpaceDim; idir++)
        {
            exactF[idir] = exactFlux(bndryLoc, idir);
        }
        Real bndryFlux = PolyGeom::dot(exactF, normal);

        update -= bndryFlux*bndryArea;
        update /= a_dx; //note NOT divided by volfrac

        a_divF(vof, 0) = update;
    }
}

void
EBScalarAdvectBC::
fluxBC(EBFluxFAB&            a_facePrim,
       const EBCellFAB&      a_Wcenter,
       const EBCellFAB&      a_Wextrap,
       const Side::LoHiSide& a_side,
       const Real&           a_time,
       const EBISBox&        a_ebisBox,
       const DataIndex&      a_dit,
       const Box&            a_box,
       const Box&            a_faceBox,
       const int&            a_faceDir)
{
//  a_facePrim.setVal(1.0); 
  CH_assert(m_isDefined);
  CH_assert(!m_domain.isPeriodic(a_faceDir));
  CH_assert(m_injectionBCFunc != NULL);

  //gets face centered region of data
  Box FBox = a_facePrim[a_faceDir].getRegion();
  Box cellBox = a_faceBox;
  CH_assert(a_facePrim[a_faceDir].nComp() == 1);  

  // Determine which side and thus shifting directions
  int isign = sign(a_side);
  cellBox.shiftHalf(a_faceDir,isign);

  //figure out if we are on a wall and if its an inflow and the inflow value
  bool isInflow   = ((a_side == Side::Lo) && (a_faceDir==m_inflowDir));
  bool isOutflow  = ((a_side == Side::Hi) && (a_faceDir==m_inflowDir));

    // Is there a domain boundary next to this grid
  if (!m_domain.contains(cellBox))
    {
      cellBox &= m_domain;
      // Find the strip of cells next to the domain boundary
      Box bndryBox = adjCellBox(cellBox, a_faceDir, a_side, 1);

      // Shift things to all line up correctly
      bndryBox.shift(a_faceDir,-isign);

      IntVectSet ivs(bndryBox);
      Real scalVal, velBC;
      for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();

          Vector<FaceIndex> bndryFaces = a_ebisBox.getFaces(vof, a_faceDir, a_side);
          for (int iface= 0; iface < bndryFaces.size(); iface++)
            {
              //this all will work for the case of spatially varying inflow
              //for inhomogeneous other faces, this is wrong.
              const FaceIndex& face = bndryFaces[iface];
              if (isOutflow) // assuming injection is from the outflow face (injectibe into counter flow)
                {
                  // figure out if you are inside or outside the tube
                  // inside the tube scalar injection val = 1.0; outside, extrap
 
                  RealVect prob_lo = RealVect::Zero;
                  const RealVect tubeCenter = m_injectionBCFunc->getTubeCenter();
                  Real tubeRadius = m_injectionBCFunc->getTubeRadius();
                  bool isInsideTube = false;
                  const RealVect loc  = EBArith::getFaceLocation(face, m_dx, prob_lo);
                  CH_assert(loc[m_inflowDir] == tubeCenter[m_inflowDir]);
                  Real radius = m_injectionBCFunc->getRadius(loc);
// start hardwire
                  Real wallThickness = 0.075;
                  ParmParse pp;
                  pp.get("wall_thickness",wallThickness);
                  Real bcFactor = 0.5; // % of outflow face to be set to inflow condition
                  tubeRadius += wallThickness * bcFactor /2.0;

                  if (radius <= tubeRadius)                 
                    {
                      isInsideTube = true;
                    }

                  isInsideTube ? scalVal = m_scalarInflowValue : scalVal = a_Wextrap(vof, 0);
                } // end if outflow
              else if (isInflow)
                {
                  scalVal = 0.0;
                }
              else  // solid wall
                {
                  scalVal = 0.0;  
                }
               a_facePrim[a_faceDir](face,0) = scalVal;
             } // end face loop
          }  // end vofit loop
       }          
}

void
EBGradDivFilter::
gradVel(EBCellFAB&             a_gradVel,
        const  EBCellFAB&      a_vel,
        const Box&             a_grid,
        const EBISBox&         a_ebisBox,
        bool a_lowOrderOneSide)
{
  CH_TIME("EBGradDivFilter::gradVel");
  CH_assert(a_gradVel.nComp() == SpaceDim*SpaceDim);
  CH_assert(a_vel.nComp() == SpaceDim);

  int iLowOrderOneSide = 0;
  if (a_lowOrderOneSide)
    {
      iLowOrderOneSide = 1;
    }

  for (int derivDir = 0; derivDir < SpaceDim; derivDir++)
    {
      Box loBox, hiBox, centerBox;
      int hasLo, hasHi;
      EBArith::loHiCenter(loBox, hasLo, hiBox, hasHi, centerBox, m_domainFine, a_grid, derivDir);
      for (int velDir = 0; velDir < SpaceDim; velDir++)
        {
          BaseFab<Real>&   regGradVel  = a_gradVel.getSingleValuedFAB();
          const BaseFab<Real>& regVel  =     a_vel.getSingleValuedFAB();
          int gradcomp = getGradComp(velDir, derivDir);

          FORT_EBGDFGRADVEL(CHF_FRA1(regGradVel, gradcomp),
                            CHF_CONST_FRA1(regVel, velDir),
                            CHF_BOX(loBox),
                            CHF_INT(hasLo),
                            CHF_BOX(hiBox),
                            CHF_INT(hasHi),
                            CHF_BOX(centerBox),
                            CHF_CONST_REAL(m_dxFine[derivDir]),
                            CHF_CONST_INT(derivDir),
                            CHF_CONST_INT(iLowOrderOneSide));

          IntVectSet ivsIrreg = a_ebisBox.getIrregIVS(a_grid);
          if (!a_lowOrderOneSide)
            {
              ivsIrreg.grow(1);
              ivsIrreg &= a_grid;
            }
          for (VoFIterator vofit(ivsIrreg, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
            {
              const VolIndex& vof = vofit();
              bool hasHi, hasLo;
              VolIndex vofHi, vofLo;
              Vector<FaceIndex> facesHi = a_ebisBox.getFaces(vof, derivDir, Side::Hi);
              Vector<FaceIndex> facesLo = a_ebisBox.getFaces(vof, derivDir, Side::Lo);

              hasHi = (facesHi.size() == 1) && (!facesHi[0].isBoundary());
              hasLo = (facesLo.size() == 1) && (!facesLo[0].isBoundary());
              if (hasLo)
                {
                  vofLo = facesLo[0].getVoF(Side::Lo);
                }
              if (hasHi)
                {
                  vofHi = facesHi[0].getVoF(Side::Hi);
                }

              Real gradVal;
              if (hasHi && hasLo)
                {
                  gradVal = (a_vel(vofHi, velDir) - a_vel(vofLo, velDir))/(2.*m_dxFine[derivDir]);
                }
              else if (hasHi || hasLo)
                {
                  //do one-sided diff
                  CH_assert(!(hasHi && hasLo));
                  Side::LoHiSide side;
                  VolIndex vofNear;
                  if (hasLo)
                    {
                      side = Side::Lo;
                      vofNear = vofLo;
                    }
                  else
                    {
                      side = Side::Hi;
                      vofNear = vofHi;
                    }
                  int isign = sign(side);

                  Vector<FaceIndex> facesFar = a_ebisBox.getFaces(vofNear, derivDir, side);
                  bool hasVoFFar = (facesFar.size()==1) && (!facesFar[0].isBoundary());
                  Real valStart = a_vel(vof,     velDir);
                  Real valNear  = a_vel(vofNear, velDir);
                  if (hasVoFFar && !a_lowOrderOneSide)
                    {
                      VolIndex vofFar = facesFar[0].getVoF(side);
                      Real valFar = a_vel(vofFar, velDir);
                      gradVal = (4*(valNear - valStart)- (valFar-valStart))/(2.*isign*m_dxFine[derivDir]);
                    }
                  else
                    {
//.........这里部分代码省略.........

void
InflowOutflowAdvectBC::
fluxBC(EBFluxFAB&            a_primGdnv,
       const EBCellFAB&      a_primCenter,
       const EBCellFAB&      a_primExtrap,
       const Side::LoHiSide& a_side,
       const Real&           a_time,
       const EBISBox&        a_ebisBox,
       const DataIndex&      a_dit,
       const Box&            a_box,
       const Box&            a_faceBox,
       const int&            a_dir)
{
  CH_assert(m_isDefined);
  CH_assert(!m_domain.isPeriodic(a_dir));

  //gets face centered region of data
  Box FBox = a_primGdnv[a_dir].getRegion();
  Box cellBox = a_faceBox;
  CH_assert(a_primGdnv[a_dir].nComp() == 1);

  // Determine which side and thus shifting directions
  int isign = sign(a_side);
  cellBox.shiftHalf(a_dir,isign);

  //figure out if we are on a wall and if its an inflow and the inflow value
  bool isInflow   = ((a_side == Side::Lo) && (a_dir==m_flowDir));
  bool isOutflow  = ((a_side == Side::Hi) && (a_dir==m_flowDir));

  // Is there a domain boundary next to this grid
  if (!m_domain.contains(cellBox))
    {
      cellBox &= m_domain;
      // Find the strip of cells next to the domain boundary
      Box bndryBox = adjCellBox(cellBox, a_dir, a_side, 1);

      // Shift things to all line up correctly
      bndryBox.shift(a_dir,-isign);

      IntVectSet ivs(bndryBox);
      for (VoFIterator vofit(ivs, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();

          Vector<FaceIndex> bndryFaces = a_ebisBox.getFaces(vof, a_dir, a_side);
          for (int iface= 0; iface < bndryFaces.size(); iface++)
            {
              //this all will work for the case of spatially varying inflow
              //for inhomogeneous other faces, this is wrong.
              const FaceIndex& face = bndryFaces[iface];
              Real velComp = 1.e99;
              if (isInflow && (m_velComp == m_flowDir))
                {
                  if (!m_doJet1PoiseInflow)
                    {
                      velComp = m_jet1inflowVel;
                    }
                  else if (m_doJet1PoiseInflow)
                    {
                      //set coordinates based on slip walls
                      RealVect prob_lo = RealVect::Zero;
                      //assumes inflow boundary is always on lo side of domain
                      const RealVect loc  = EBArith::getFaceLocation(face, m_dx, prob_lo);
                      Real radius = m_jet1PoiseInflowFunc->getRadius(loc);
                      velComp = m_jet1PoiseInflowFunc->getVel(radius)[m_flowDir];
                    }
                }
              else if (isInflow && (m_velComp != m_flowDir))
                {
                  velComp = 0.0;
                }
              else if (isOutflow)
                {
                   if (!m_doJet2)
                     {
                       velComp = a_primExtrap(vof, 0);
                     }
                   else if (m_doJet2)
                     {
                       RealVect prob_lo = RealVect::Zero;
                       const RealVect tubeCenter = m_jet2PoiseInflowFunc->getTubeCenter();
                       Real tubeRadius = m_jet2PoiseInflowFunc->getTubeRadius();
                       bool isInsideTube = false;
                       const RealVect loc  = EBArith::getFaceLocation(face, m_dx, prob_lo);                        
                       CH_assert(loc[m_flowDir] == tubeCenter[m_flowDir]);

                       Real radius = m_jet2PoiseInflowFunc->getRadius(loc);
// start hardwire
                       Real wallThickness = 0.075;
                       ParmParse pp;
                       pp.get("wall_thickness",wallThickness);
                       Real bcFactor = 0.5; // % of outflow face to be set to inflow condition
                       tubeRadius += wallThickness * bcFactor /2.0;
// end hardwire
                       if (radius <= tubeRadius)
                         {
                           isInsideTube = true;
                         }
                       if (isInsideTube && (m_velComp == m_flowDir))
                         {
//.........这里部分代码省略.........

void
EBGradDivFilter::
cellGradient(EBCellFAB&             a_gradDiv,
             const EBCellFAB&       a_gradVel,
             const  EBCellFAB&      a_vel,
             const EBFluxFAB&       a_fluxVel,
             const EBFluxFAB&       a_div,
             const Box&             a_grid,
             const EBISBox&         a_ebisBox,
             const DataIndex&       a_dataIndex,
             bool a_multiplyByLambda,
             bool a_noExtrapToCovered)
{
  CH_TIME("EBGradDivFilter::cellGradient");
  int icomp = 0;
  EBCellFAB& lambdaFAB = m_lambda[a_dataIndex];
  IntVectSet irregIVS = a_ebisBox.getIrregIVS(a_grid);

  for (int faceDir = 0; faceDir < SpaceDim; faceDir++)
    {
      const EBFaceFAB&  faceDiv       =  a_div[faceDir];
      const BaseFab<Real>& regFaceDiv =   faceDiv.getSingleValuedFAB();
      const BaseFab<Real>& regLambda =   lambdaFAB.getSingleValuedFAB();
      BaseFab<Real>& regGradDiv       = a_gradDiv.getSingleValuedFAB();

      int imultByLambda = 0;
      if (a_multiplyByLambda)
        {
          imultByLambda = 1;
        }
      FORT_EBGDFCELLGRAD(CHF_FRA1(regGradDiv, faceDir),
                         CHF_CONST_FRA1(regFaceDiv, icomp),
                         CHF_CONST_FRA1(regLambda, icomp),
                         CHF_BOX(a_grid),
                         CHF_CONST_REAL(m_dxFine[faceDir]),
                         CHF_CONST_INT(faceDir),
                         CHF_CONST_INT(imultByLambda));

      //nonconservative gradient
      for (VoFIterator vofit(irregIVS, a_ebisBox.getEBGraph()); vofit.ok(); ++vofit)
        {
          const VolIndex& vof = vofit();
          Real gradVal = 0.0;

          Vector<FaceIndex> facesHi = a_ebisBox.getFaces(vof, faceDir, Side::Hi);
          Vector<FaceIndex> facesLo = a_ebisBox.getFaces(vof, faceDir, Side::Lo);

          bool zeroGrad = a_noExtrapToCovered && ((facesHi.size() == 0) || (facesLo.size() == 0));
          if (zeroGrad)
            {
              gradVal = 0;
            }
          else
            {
              Real valHi= 0;
              if (facesHi.size() > 0)
                {
                  for (int iface = 0; iface < facesHi.size(); iface++)
                    {
                      valHi += faceDiv(facesHi[iface], icomp);
                    }
                  valHi /= facesHi.size();
                }
              else
                {
                  valHi = ccpGetCoveredExtrapValue(vof, faceDir, Side::Hi,
                                                   faceDiv, a_ebisBox, a_grid,
                                                   m_domainFine, m_dxFine, icomp);
                }

              Real valLo= 0;
              if (facesLo.size() > 0)
                {
                  for (int iface = 0; iface < facesLo.size(); iface++)
                    {
                      valLo += faceDiv(facesLo[iface], 0);
                    }
                  valLo /= facesLo.size();
                }
              else
                {
                  valLo = ccpGetCoveredExtrapValue(vof, faceDir, Side::Lo ,
                                                   faceDiv, a_ebisBox, a_grid,
                                                   m_domainFine, m_dxFine, icomp);

                }
              gradVal = (valHi-valLo)/(m_dxFine[faceDir]);

              //multiply the lambda in since we have the area fractions lying about
              if (a_multiplyByLambda)
                {
                  Real lambda = lambdaFAB(vof, 0);
                  gradVal *= lambda;
                }
            }

          a_gradDiv(vof, faceDir) = gradVal;
        }
    }
}

void
EBCompositeMACProjector::
correctTangentialVelocity(EBFaceFAB&        a_velocity,
                          const EBFaceFAB&  a_gradient,
                          const Box&        a_grid,
                          const EBISBox&    a_ebisBox,
                          const IntVectSet& a_cfivs)
{
  CH_TIME("EBCompositeMACProjector::correctTangentialVelocity");
  int velDir = a_velocity.direction();
  int gradDir = a_gradient.direction();
  //veldir is the face on which the velocity lives
  //graddir is the direction of the component
  //and the face on which the mac gradient lives
  CH_assert(velDir != gradDir);
  CH_assert(a_velocity.nComp() == 1);
  CH_assert(a_gradient.nComp() == 1);

  //updating in place in fortran so we have to save a acopy
  EBFaceFAB velSave(a_ebisBox, a_velocity.getCellRegion(),  velDir, 1);
  velSave.copy(a_velocity);

  //interior box is the box where all of the stencil can be reached
  //without going out of the domain
  ProblemDomain domainBox = a_ebisBox.getDomain();
  Box interiorBox = a_grid;
  interiorBox.grow(1);
  interiorBox &= domainBox;
  interiorBox.grow(-1);
  Box interiorFaceBox = surroundingNodes(interiorBox, velDir);

  if (!interiorBox.isEmpty())
    {
      BaseFab<Real>&       regVel  = a_velocity.getSingleValuedFAB();
      const BaseFab<Real>& regGrad = a_gradient.getSingleValuedFAB();
      FORT_REGCORRECTTANVEL(CHF_FRA1(regVel,0),
                            CHF_CONST_FRA1(regGrad,0),
                            CHF_BOX(interiorFaceBox),
                            CHF_INT(velDir),
                            CHF_INT(gradDir));
    }

  //do only irregular and boundary cells pointwise
  IntVectSet ivsIrreg = a_ebisBox.getIrregIVS(a_grid);
  IntVectSet ivsGrid(a_grid);
  ivsGrid -= interiorBox;
  ivsGrid |= ivsIrreg;
  FaceIterator faceit(ivsGrid, a_ebisBox.getEBGraph(), velDir, FaceStop::SurroundingWithBoundary);
  for (faceit.reset(); faceit.ok(); ++faceit)
    {
      //average neighboring grads in grad dir direction
      int numGrads = 0;
      Real gradAve = 0.0;
      const FaceIndex& velFace = faceit();
      for (SideIterator sitVel; sitVel.ok(); ++sitVel)
        {
          const VolIndex& vofSide = velFace.getVoF(sitVel());
          const IntVect&   ivSide = vofSide.gridIndex();
          //do not include stuff over coarse-fine interface and inside the domain
          //cfivs includes cells just outside domain
          if (!a_cfivs.contains(ivSide) && domainBox.contains(ivSide))
            {
              for (SideIterator sitGrad; sitGrad.ok(); ++sitGrad)
                {
                  Vector<FaceIndex> gradFaces = a_ebisBox.getFaces(vofSide, gradDir, sitGrad());
                  for (int iface = 0; iface < gradFaces.size(); iface++)
                    {
                      if (!gradFaces[iface].isBoundary())
                        {
                          numGrads++;
                          gradAve += a_gradient(gradFaces[iface], 0);
                        }
                    }
                }//end loop over gradient sides
            }//end cfivs check
          else
            {
              // inside coarse/fine interface or at domain boundary. extrapolate from neighboring vofs to get the gradient
              const Side::LoHiSide inSide    = flip(sitVel());
              const VolIndex&      inSideVof = velFace.getVoF(inSide);
              const IntVect&       inSideIV  = inSideVof.gridIndex();

              IntVect inSideFarIV = inSideIV;
              inSideFarIV[velDir] += sign(inSide);
              if (domainBox.contains(inSideIV) && domainBox.contains(inSideFarIV))
                {
                  Vector<VolIndex> farVofs = a_ebisBox.getVoFs(inSideFarIV);
                  if (farVofs.size()  == 1)
                    {
                      const VolIndex& inSideFarVof = farVofs[0];
                      for (SideIterator sitGrad; sitGrad.ok(); ++sitGrad)
                        {
                          //get the grad for the face adjoining inSideVof on the sitGrad side, in the gradDir direction
                          Vector<FaceIndex> gradFaces    = a_ebisBox.getFaces(inSideVof   , gradDir, sitGrad());
                          Vector<FaceIndex> gradFarFaces = a_ebisBox.getFaces(inSideFarVof, gradDir, sitGrad());
                          if ( (gradFaces.size() == 1) && (gradFarFaces.size() == 1) )
                            {
                              // if ( (!gradFaces[0].isBoundary()) && (!gradFarFaces[0].isBoundary()) )
                              //   {
                                  const Real& inSideGrad    = a_gradient(gradFaces[0], 0);
//.........这里部分代码省略.........

int checkEBISBox(const Box& a_gridCoar, const EBISBox& a_ebisBoxCoar, const EBISBox& a_ebisBoxFine)
{
  IntVectSet ivs = a_ebisBoxCoar.getIrregIVS(a_gridCoar);
  Real dxCoar = 2;  Real dxFine = 1;

#if CH_SPACEDIM==2
  Real areaFineCell = dxFine;
  Real areaCoarCell = dxCoar;
  Real voluFineCell = dxFine*dxFine;
  Real voluCoarCell = dxCoar*dxCoar;
#elif CH_SPACEDIM==3
  Real areaFineCell = dxFine*dxFine;
  Real areaCoarCell = dxCoar*dxCoar;
  Real voluFineCell = dxFine*dxFine*dxFine;
  Real voluCoarCell = dxCoar*dxCoar*dxCoar;
#else
  MayDay::Error();
#endif
  int retval = 0;
  for (VoFIterator vofit(ivs, a_ebisBoxCoar.getEBGraph()); vofit.ok(); ++vofit)
    {
      const VolIndex&  vofCoar = vofit();
      Vector<VolIndex> vofsFine = a_ebisBoxCoar.refine(vofCoar);

      //check the easy bits
      Real volumCoar = a_ebisBoxCoar.volFrac(  vofCoar);
      RealVect areaCritCoar = a_ebisBoxCoar.bndryArea(vofCoar)*
                              a_ebisBoxCoar.normal(vofCoar);
      Real volumFine = 0;
      RealVect areaCritFine = RealVect::Zero;
      for (int ivof = 0; ivof < vofsFine.size(); ivof++)
        {
          volumFine += a_ebisBoxFine.volFrac(  vofsFine[ivof]);
          areaCritFine += a_ebisBoxFine.bndryArea(vofsFine[ivof])*
                          a_ebisBoxFine.normal(vofsFine[ivof]);
        }
      volumFine *= voluFineCell;
      areaCritFine *= areaFineCell;
      volumCoar *= voluCoarCell;
      areaCritCoar *= areaCoarCell;
      Real tolerance = 1.0e-10;
      if (Abs(volumFine -volumCoar) > tolerance*volumCoar)
        {
          pout() << "volume problem in coar cell " << vofCoar.gridIndex() << endl;
          retval = -1;
        }

      // Real maxCc  = 0.0;
      // Real maxDev = 0.0;
      for (int idir=0; idir<SpaceDim; idir++)
        {
          if (Abs(areaCritFine[idir]-areaCritCoar[idir]) >
             tolerance*Abs(areaCritCoar[idir]))
            {
              pout() << "bndry area problem in coar cell " << vofCoar.gridIndex() << endl;
              retval = -2;
            }
        }
      //centroids are a bit uglier to test
//      RealVect bndryCentroidCoar = a_ebisBoxCoar.bndryCentroid(  vofCoar);
//      RealVect bndryCentroidFine = RealVect::Zero;

      //the areas are somewhat more painful to test
//      for (int idir = 0; idir < SpaceDim; idir++)
//        {
//
//        }
    }

  return retval;
}

void
EBGradDivFilter::
faceDivergence(EBFaceFAB&             a_divVel,
               const EBCellFAB&       a_gradVel,
               const EBCellFAB&       a_vel,
               const EBFluxFAB&       a_fluxVel,
               const Box&             a_grid,
               const EBISBox&         a_ebisBox,
               const int&             a_faceDir)
{
  CH_TIME("EBGradDivFilter::faceDivergence");
  CH_assert(a_divVel.nComp() == 1);
  CH_assert(a_vel.nComp() == SpaceDim);

  a_divVel.setVal(0.0);
  BaseFab<Real>& regDivVel       =  a_divVel.getSingleValuedFAB();
  const BaseFab<Real>& regVel    =     a_vel.getSingleValuedFAB();
  const BaseFab<Real>& regGradVel= a_gradVel.getSingleValuedFAB();

  Box interiorFaceBox = a_grid;
  interiorFaceBox.grow(a_faceDir, 1);
  interiorFaceBox &= m_domainFine;
  interiorFaceBox.grow(a_faceDir, -1);
  interiorFaceBox.surroundingNodes(a_faceDir);

  for (int divDir = 0; divDir < SpaceDim; divDir++)
    {
      FORT_EBGDFFACEDIVINCR(CHF_FRA1(regDivVel, 0),
                            CHF_FRA(regVel),
                            CHF_FRA(regGradVel),
                            CHF_BOX(interiorFaceBox),
                            CHF_REAL(m_dxFine[divDir]),
                            CHF_INT(a_faceDir),
                            CHF_INT(divDir));
    }

  IntVectSet irregIVS = a_ebisBox.getIrregIVS(a_grid);
  FaceStop::WhichFaces stopCrit;
  if (m_domainFine.isPeriodic(a_faceDir))
    {
      stopCrit = FaceStop::SurroundingWithBoundary;
    }
  else
    {
      stopCrit = FaceStop::SurroundingNoBoundary;
    }
  for (FaceIterator faceit(irregIVS, a_ebisBox.getEBGraph(), a_faceDir,
                          stopCrit);
      faceit.ok(); ++faceit)
    {
      Real divVal = 0.0;
      for (int divDir=0; divDir<SpaceDim; divDir++)
        {
          if (divDir == a_faceDir)
            {
              divVal += (a_vel(faceit().getVoF(Side::Hi), a_faceDir) -
                         a_vel(faceit().getVoF(Side::Lo), a_faceDir))/m_dxFine[divDir];
            }
          else
            {
              // take average of cell-centered tangential derivatives
              // and increment div
              //so this is partial (vel_divdir)/partial (x_divdir)
              int velcomp = divDir;
              int gradcomp = getGradComp(velcomp, divDir);

              divVal += 0.5*(a_gradVel(faceit().getVoF(Side::Hi), gradcomp) +
                             a_gradVel(faceit().getVoF(Side::Lo), gradcomp));
            }
        }

      a_divVel(faceit(), 0) = divVal;
    } // end loop over interior irregular faces

  if (!m_domainFine.isPeriodic(a_faceDir))
    {
      //set the boundary face divergence to an extrapolation of the neighboring face divergences
      for (SideIterator sit; sit.ok(); ++sit)
        {
          Box bndryBox = adjCellBox(a_grid, a_faceDir, sit(), 1);
          int ishift = -sign(sit());
          bndryBox.shift(a_faceDir, ishift);
          IntVectSet bndryIVS(bndryBox);
          for (FaceIterator faceit(bndryIVS, a_ebisBox.getEBGraph(), a_faceDir,
                                  FaceStop::AllBoundaryOnly);
              faceit.ok(); ++faceit)
            {
              Real faceDiv = getDomainDivergence(a_gradVel,
                                                 a_vel,
                                                 a_fluxVel,
                                                 a_grid,
                                                 a_ebisBox,
                                                 a_faceDir,
                                                 faceit(),
                                                 sit());
              a_divVel(faceit(), 0) = faceDiv;
            }
        }
    }
}