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

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

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