C++ FArrayBox::nComp方法代码示例

void ConstBCFunction::operator()(FArrayBox&           a_state,
                                 const Box&           a_valid,
                                 const ProblemDomain& a_domain,
                                 Real                 a_dx,
                                 bool                 a_homogeneous)
{
  const Box& domainBox = a_domain.domainBox();

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    if (!a_domain.isPeriodic(idir))
    {
      for (SideIterator sit; sit.ok(); ++sit)
      {
        Side::LoHiSide side = sit();

        if (a_valid.sideEnd(side)[idir] == domainBox.sideEnd(side)[idir])
        {
          int  bcType;
          Real bcValue;

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

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

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

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

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

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

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

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

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

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

              Real inhomogVal = 0.0;

              if (!a_homogeneous)
              {
                inhomogVal = bcValue;
              }

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

                Real ghostVal = linearInterp(inhomogVal, nearVal);

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

void
LevelFluxRegisterEdge::incrementFine(
                                     FArrayBox& a_fineFlux,
                                     Real a_scale,
                                     const DataIndex& a_fineDataIndex,
                                     const Interval& a_srcInterval,
                                     const Interval& a_dstInterval,
                                     int a_dir,
                                     Side::LoHiSide a_sd)
{
  CH_assert(isDefined());
  CH_assert(!a_fineFlux.box().isEmpty());
  CH_assert(a_srcInterval.size() == a_dstInterval.size());
  CH_assert(a_srcInterval.begin() >= 0);
  CH_assert(a_srcInterval.end() < a_fineFlux.nComp());
  CH_assert(a_dstInterval.begin() >= 0);
  CH_assert(a_dstInterval.end() < m_nComp);

  CH_assert(a_dir >= 0);
  CH_assert(a_dir < SpaceDim);
  CH_assert((a_sd == Side::Lo)||(a_sd == Side::Hi));
  //
  //
  //denom is the number of fine faces per coarse face
  //this is intrinsically dimension-dependent
#if (CH_SPACEDIM == 2)
  Real denom = 1;
#elif (CH_SPACEDIM == 3)
  Real denom = m_nRefine;
#else
  // This code doesn't make any sense in 1D, and hasn't been implemented
  // for DIM > 3
  Real denom = -1.0;
  MayDay::Error("LevelFluxRegisterEdge -- bad SpaceDim");
#endif

  Real scale = a_scale/denom;

  // need which fluxbox face we're doing this for
  Box thisBox = a_fineFlux.box();
  int fluxComp = -1;
  for (int sideDir=0; sideDir<SpaceDim; sideDir++)
  {
    // we do nothing in the direction normal to face
    if (sideDir != a_dir)
    {
      if (thisBox.type(sideDir) == IndexType::CELL)
      {
        fluxComp = sideDir;
      }
    }
  }
  CH_assert (fluxComp >= 0);
  int regcomp = getRegComp(a_dir, fluxComp);

  FluxBox& thisReg = m_fabFine[index(a_dir, a_sd)][a_fineDataIndex];
  FArrayBox& reg = thisReg[regcomp];

  a_fineFlux.shiftHalf(a_dir, sign(a_sd));

  // this is a way of geting a face-centered domain
  // box which we can then use to intersect with things
  // to screen out cells outside the physical domain
  // (nothing is screened out in periodic case)
  Box shiftedValidDomain = m_domainCoarse.domainBox();
  shiftedValidDomain.grow(2);
  shiftedValidDomain &= m_domainCoarse;
  shiftedValidDomain.surroundingNodes(regcomp);

  BoxIterator regIt(reg.box() & shiftedValidDomain);
  for (regIt.begin(); regIt.ok(); ++regIt)
    {
      const IntVect& coarseIndex = regIt();
      // create a cell-centered box, then shift back to face-centered
      Box box(coarseIndex, coarseIndex);
      box.shiftHalf(regcomp,-1);
      // to avoid adding in edges which do not overlie coarse-grid
      // edges, will refine only in non-fluxComp directions to
      // determine box from which to grab fluxes.
      IntVect refineVect(m_nRefine*IntVect::Unit);
      //refineVect.setVal(fluxComp,1);
      box.refine(refineVect);
      if (a_sd == Side::Lo) box.growLo(a_dir,-(m_nRefine-1));
      else                 box.growHi(a_dir,-(m_nRefine-1));
      BoxIterator fluxIt(box);
      for (fluxIt.begin(); fluxIt.ok(); ++fluxIt)
        {
          int src = a_srcInterval.begin();
          int dest =  a_dstInterval.begin();
          for ( ; src <=a_srcInterval.end(); ++src,++dest)
            reg(coarseIndex, dest) += scale*a_fineFlux(fluxIt(), src);
        }
    }
  a_fineFlux.shiftHalf(a_dir, -sign(a_sd));
}

//.........这里部分代码省略.........

  // Compute flattening once for all slopes if needed
  FArrayBox flattening(slopeBox, 1); // cell-centered
  if (m_useFlattening)
    {
      Interval velInt    = m_gdnvPhysics->velocityInterval();
      int pressureIndex  = m_gdnvPhysics->pressureIndex();
      Real smallPressure = m_gdnvPhysics->smallPressure();
      int bulkIndex      = m_gdnvPhysics->bulkModulusIndex();

      m_util.computeFlattening(flattening,
                               W,
                               velInt,
                               pressureIndex,
                               smallPressure,
                               bulkIndex,
                               slopeBox);
    }

  // Intermediate, extrapolated primitive variables
  FArrayBox WMinus[SpaceDim];
  FArrayBox WPlus [SpaceDim];

  // Initial fluxes
  FArrayBox WHalf1[SpaceDim];

  // Source term computed from current state
  FArrayBox localSource;

  // If the source term is valid, make a local copy, increment it, and scale
  // it by 1/2 the timestep
  if (!a_S.box().isEmpty())
    {
      localSource.define(a_S.box(), a_S.nComp());
      localSource.copy(a_S);

      m_gdnvPhysics->incrementSource(localSource,W,slopeBox);

      localSource *= 0.5 * a_dt;
    }

  // Compute initial fluxes
  for (int dir1 = 0; dir1 < SpaceDim; dir1++)
    {
      // Size the intermediate, extrapolated primitive variables
      WMinus[dir1].resize(slopeBox,numPrim); // cell-centered
      WPlus [dir1].resize(slopeBox,numPrim); // cell-centered

      // Compute predictor step to obtain extrapolated primitive variables
      if (m_normalPredOrder == 0)
        {
          CTUNormalPred(WMinus[dir1],WPlus[dir1],a_dt,m_dx,W,
                        flattening,dir1,slopeBox);
        }
      else if (m_normalPredOrder == 1)
        {
          PLMNormalPred(WMinus[dir1],WPlus[dir1],a_dt,m_dx,W,
                        flattening,dir1,slopeBox);
        }
      else if (m_normalPredOrder == 2)
        {
          PPMNormalPred(WMinus[dir1],WPlus[dir1],a_dt,m_dx,W,
                        flattening,dir1,slopeBox);
        }
      else
        {