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

void EBPoissonOp::
getOpVoFStencil(VoFStencil&             a_stencil,
                const int&              a_idir,
                const Vector<VolIndex>& a_allMonotoneVoFs,
                const EBISBox&          a_ebisbox,
                const VolIndex&         a_VoF,
                const bool&             a_lowOrder)
{

  for (SideIterator sit; sit.ok(); ++sit)
    {
      Side::LoHiSide side = sit();
      Vector<FaceIndex> faces;

      faces = a_ebisbox.getFaces(a_VoF,a_idir,side);

      for (int iface = 0; iface < faces.size(); iface++)
        {
          FaceIndex face = faces[iface];
          VoFStencil faceStencil;

          getOpFaceStencil(faceStencil,a_allMonotoneVoFs,a_ebisbox,a_VoF,
                           a_idir,side,face,a_lowOrder);

          a_stencil += faceStencil;
        }
    }
}

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
EBGradDivFilter::
getAreaFracs(Vector<FaceIndex> a_facesLo[SpaceDim],
             Vector<FaceIndex> a_facesHi[SpaceDim],
             bool              a_hasFacesLo[SpaceDim],
             bool              a_hasFacesHi[SpaceDim],
             RealVect&         a_areaFracLo,
             RealVect&         a_areaFracHi,
             const VolIndex&   a_vof,
             const EBISBox&    a_ebisBox)
{
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      Real areaSumLo = 0;
      Real areaSumHi = 0;
      a_facesLo[idir] = a_ebisBox.getFaces(a_vof, idir, Side::Lo);
      a_facesHi[idir] = a_ebisBox.getFaces(a_vof, idir, Side::Hi);
      a_hasFacesLo[idir] = ( ( a_facesLo[idir].size() > 0) &&
                             (!a_facesLo[idir][0].isBoundary()));
      a_hasFacesHi[idir] = ( ( a_facesHi[idir].size() > 0) &&
                             (!a_facesHi[idir][0].isBoundary()));
      //want boundary faces to look covered
      if (a_hasFacesLo[idir])
        {
          for (int iface = 0; iface < a_facesLo[idir].size(); iface++)
            {
              areaSumLo += a_ebisBox.areaFrac(a_facesLo[idir][iface]);
            }
        }
      //want boundary faces to look covered
      if (a_hasFacesHi[idir])
        {
          for (int iface = 0; iface < a_facesHi[idir].size(); iface++)
            {
              areaSumHi += a_ebisBox.areaFrac(a_facesHi[idir][iface]);
            }
        }

      a_areaFracLo[idir] = areaSumLo;
      a_areaFracHi[idir] = areaSumHi;
    }
}

Real
divergence(const EBFluxFAB&  a_func,
           const RealVect&   a_bndryFlux,
           const EBISBox&    a_ebisBox,
           const VolIndex&   a_vof,
           const Real&       a_dx)
{
  Real retval = 0;

  Real bndryArea = a_ebisBox.bndryArea(a_vof);
  RealVect normal = a_ebisBox.normal(a_vof);
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      Real bndryFlux    = a_bndryFlux[idir]; //normal already dealt with
      Real bndryContrib = -bndryFlux*bndryArea*normal[idir];
      Real openContrib = 0;
      if (bndryArea > 1.0e-3)
        {
          openContrib = 0;
        }
      for (SideIterator sit; sit.ok(); ++sit)
        {
          int isign = sign(sit());
          Real rsign = isign;
          Vector<FaceIndex> faces = a_ebisBox.getFaces(a_vof, idir, sit());
          for (int iface = 0; iface < faces.size(); ++iface)
            {
              Real areaFrac = a_ebisBox.areaFrac(faces[iface]);
              Real faceFlux = a_func[idir](faces[iface], 0);
              openContrib += rsign*faceFlux*areaFrac;
            }
        }
      retval += openContrib + bndryContrib;

    }
  retval /= a_dx;
  return retval;
}

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);
//.........这里部分代码省略.........

void
EBCompositeMACProjector::
correctVelocityComponent(BaseIVFAB<Real>       &  a_coveredVel,
                         const Vector<VolIndex>&  a_coveredFace,
                         const IntVectSet      &  a_coveredSets,
                         const EBFluxFAB       &  a_macGradient,
                         const EBISBox         &  a_ebisBox,
                         int    a_coveredFaceDir,  Side::LoHiSide a_sd,    int a_faceGradComp)
{
  CH_TIME("EBCompositeMACProjector::correctVelocityComponent");
  //if a_coveredFaceDir  == faceGradComp, just linearly extrapolate gradient and subtract.
  //if a_coveredFaceDir != faceGradComp, need to average to cell centers then extrapolate and subtract
  CH_assert(a_coveredVel.nComp() == 1);
  CH_assert(a_macGradient.nComp() == 1);
  const EBFaceFAB& macGradFAB = a_macGradient[a_faceGradComp];
  for (int ivof = 0; ivof < a_coveredFace.size(); ivof++)
    {
      const VolIndex& vof = a_coveredFace[ivof];
      Vector<FaceIndex> nearFaces, farFaces;
      bool hasNearFaces = false;
      bool hasFarFaces  = false;
      VolIndex  nearVoF, farVoF;
      FaceIndex nearFace, farFace;

      nearFaces = a_ebisBox.getFaces(vof, a_coveredFaceDir, flip(a_sd));
      hasNearFaces = (nearFaces.size()==1) && (!nearFaces[0].isBoundary());
      if (hasNearFaces)
        {
          nearFace = nearFaces[0];
          nearVoF = nearFace.getVoF(flip(a_sd));

          farFaces = a_ebisBox.getFaces(nearVoF, a_coveredFaceDir, flip(a_sd));
          hasFarFaces = (farFaces.size()==1)  && (!farFaces[0].isBoundary());
          if (hasFarFaces)
            {
              farFace = farFaces[0];
              farVoF  = farFace.getVoF(flip(a_sd));
            }
        }
      Real extrapGrad;
      if (a_coveredFaceDir == a_faceGradComp)
        {
          //if a_coveredFaceDir  == faceGradComp, just linearly extrapolate gradient and subtract.
          if (hasNearFaces && hasFarFaces)
            {
              Real nearGrad  = macGradFAB(nearFace, 0);
              Real farGrad   = macGradFAB(farFace,  0);
              extrapGrad = 2.0*nearGrad - farGrad;
            }
          else if (hasNearFaces)
            {
              extrapGrad  = macGradFAB(nearFace, 0);
            }
          else
            {
              extrapGrad = 0.0;
            }

        }
      else
        {
          //if a_coveredFaceDir != faceGradComp, need to average to cell centers then extrapolate and subtract
          Real nearGrad = getAverageFaceGrad(vof,     macGradFAB, a_ebisBox, a_faceGradComp);
          if (hasNearFaces)
            {
              Real farGrad =  getAverageFaceGrad(nearVoF, macGradFAB, a_ebisBox, a_faceGradComp);
              extrapGrad = 1.5*nearGrad - 0.5*farGrad;
            }
          else
            {
              extrapGrad = nearGrad;
            }
        }
      a_coveredVel(vof, 0) -= extrapGrad;
    }
}

Real
getAverageFaceGrad(const VolIndex&   a_vof,
                   const EBFaceFAB&  a_macGradient,
                   const EBISBox&    a_ebisBox,
                   int               a_faceDir)
{
  CH_TIME("EBCompositeMACProjector::getAverageFaceGrad");
  Vector<FaceIndex> hiFaces = a_ebisBox.getFaces(a_vof, a_faceDir, Side::Hi);
  Vector<FaceIndex> loFaces = a_ebisBox.getFaces(a_vof, a_faceDir, Side::Lo);
  bool hasHiFaces = (hiFaces.size() > 0);
  bool hasLoFaces = (loFaces.size() > 0);
  Real hiVal = 0.0;
  Real loVal = 0.0;
  if (hasHiFaces)
    {
      for (int iface = 0; iface < hiFaces.size(); iface++)
        {
          hiVal += a_macGradient(hiFaces[iface], 0);
        }
      hiVal /= hiFaces.size();

    }
  if (hasLoFaces)
    {
      for (int iface = 0; iface < loFaces.size(); iface++)
        {
          loVal += a_macGradient(loFaces[iface], 0);
        }
      loVal /= loFaces.size();
    }
  if (!hasLoFaces && hasHiFaces)
    {
      //only want to do the wacky extrapolation thing if  not multivalued
      if ((hiFaces.size() == 1) && (!hiFaces[0].isBoundary()))
        {
          const VolIndex& hiVoF = hiFaces[0].getVoF(Side::Hi);
          Vector<FaceIndex> farFaces = a_ebisBox.getFaces(hiVoF, a_faceDir, Side::Hi);
          if (farFaces.size() > 0)
            {
              Real farVal = 0.0;
              for (int iface = 0; iface < farFaces.size(); iface++)
                {
                  farVal += a_macGradient(farFaces[iface], 0);
                }
              farVal /= farFaces.size();
              loVal = 2*hiVal - farVal;
            }
          else
            {
              loVal = hiVal;
            }
        }
      else
        {
          loVal = hiVal;
        }
    }


  if (!hasHiFaces && hasLoFaces)
    {
      //only want to do the wacky extrapolation thing if not multivalued
      if ((loFaces.size() == 1) && (!loFaces[0].isBoundary()))
        {
          const VolIndex& hiVoF = loFaces[0].getVoF(Side::Lo);
          Vector<FaceIndex> farFaces = a_ebisBox.getFaces(hiVoF, a_faceDir, Side::Lo);
          if (farFaces.size() > 0)
            {
              Real farVal = 0.0;
              for (int iface = 0; iface < farFaces.size(); iface++)
                {
                  farVal += a_macGradient(farFaces[iface], 0);
                }
              farVal /= farFaces.size();
              hiVal = 2*loVal - farVal;
            }
          else
            {
              hiVal = loVal;
            }
        }
      else
        {
          hiVal = loVal;
        }
    }

  Real retval = 0.5*(hiVal + loVal);
  return retval;
}

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
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
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
                    {
//.........这里部分代码省略.........

Real
EBGradDivFilter::
getDomainDivergence(const EBCellFAB&       a_gradVel,
                    const EBCellFAB&       a_vel,
                    const EBFluxFAB&       a_fluxVel,
                    const Box&             a_grid,
                    const EBISBox&         a_ebisBox,
                    const int&             a_faceDir,
                    const FaceIndex&       a_face,
                    const Side::LoHiSide&  a_side)
{
  CH_TIME("EBGradDivFilter::getDomainDivergence");
  Real divVel;
  CH_assert(a_fluxVel.nComp() == 1);
  Real velB = a_fluxVel[a_faceDir](a_face, 0);
  //compute normal derivative.
  //need second-order derivative at the boundary so
  //given u_b, u_i, u_i+1,
  //ux_b = (9u_i - u_i+1 - 8u_b)/(6 dx)
  VolIndex  vofNear, vofStart;
  Real valNear=0;
  Real valStart =0;
  bool hasVoFNear;
  vofStart = a_face.getVoF(flip(a_side));
  valStart = a_vel(vofStart, a_faceDir);
  int iflipsign = sign(flip(a_side));
  Real rsign = Real(iflipsign);

  Vector<FaceIndex> facesNear = a_ebisBox.getFaces(vofStart, a_faceDir, flip(a_side));
  hasVoFNear = (facesNear.size() == 1);
  if (hasVoFNear)
    {
      vofNear = facesNear[0].getVoF(flip(a_side));
      valNear = a_vel(vofNear, a_faceDir);
    }

  Real hsign = rsign*m_dxFine[a_faceDir];
  if (hasVoFNear)
    {
      divVel  = (2.25*(valStart-velB)  - 0.25*(valNear-velB))/(0.75*hsign);
    }
  else
    {
      //if there is no farther away value, drop order of derivative
      divVel = (valStart - velB)/(0.5*hsign);
    }

  //now add in tangential derivatives of tangential velocities
  //using extrapolated gradients
  for (int tanDir = 0; tanDir< SpaceDim; tanDir++)
    {
      if (tanDir != a_faceDir)
        {
          Real gradVel = 0.0;
          int gradcomp = getGradComp(tanDir, tanDir);
          if (hasVoFNear)
            {
              Real gradNear = a_gradVel(vofNear,   gradcomp);
              Real gradStart = a_gradVel(vofStart, gradcomp);
              gradVel = 1.5*gradStart - 0.5*gradNear;
            }
          else
            {
              Real gradStart = a_gradVel(vofStart, gradcomp);
              gradVel = gradStart;
            }
          divVel += gradVel;
        }
    }

  return divVel;

}