C++ ProblemDomain::isPeriodic方法代码示例

void DenseIntVectSet::nestingRegion(int radius, const ProblemDomain& a_domain)
{
  CH_assert(radius >= 0);
  if (radius == 0) return;

  DenseIntVectSet tmp;
  Box region = a_domain.domainBox();
  if (!a_domain.isPeriodic())
    {
      tmp = *this;
    }
  else
    {
      D_TERM6(if (a_domain.isPeriodic(0)) region.grow(0, radius);,
              if (a_domain.isPeriodic(1)) region.grow(1, radius);,
              if (a_domain.isPeriodic(2)) region.grow(2, radius);,

void CFStencil::buildPeriodicVector(Vector<Box>& a_periodicVector,
                                    const ProblemDomain& a_fineDomain,
                                    const DisjointBoxLayout& a_fineBoxes)
{
  Box periodicTestBox(a_fineDomain.domainBox());
  if (a_fineDomain.isPeriodic())
    {
      for (int idir=0; idir<SpaceDim; idir++)
        {
          if (a_fineDomain.isPeriodic(idir))
            {
              periodicTestBox.grow(idir,-1);
            }
        }
    }
  a_periodicVector.clear();
  a_periodicVector.reserve(a_fineBoxes.size());

  LayoutIterator lit = a_fineBoxes.layoutIterator();
  for (lit.reset(); lit.ok(); ++lit)
    {
      const Box& box = a_fineBoxes[lit()];
      a_periodicVector.push_back(box);
      // if periodic, also need to add periodic images
      // only do this IF we're periodic and box
      // adjacent to  the domain box boundary somewhere
      if (a_fineDomain.isPeriodic()
          && !periodicTestBox.contains(box))
        {
          ShiftIterator shiftIt = a_fineDomain.shiftIterator();
          IntVect shiftMult(a_fineDomain.domainBox().size());
          Box shiftedBox(box);
          for (shiftIt.begin(); shiftIt.ok(); ++shiftIt)
            {
              IntVect shiftVect = shiftMult*shiftIt();
              shiftedBox.shift(shiftVect);
              a_periodicVector.push_back(shiftedBox);
              shiftedBox.shift(-shiftVect);
            } // end loop over periodic shift directions
        } // end if periodic
    }
  a_periodicVector.sort();
}

void 
CompGridVTOBC::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])
                {
                  // Dirichlet BC
                  int isign = sign(side);
                  Box toRegion = adjCellBox(a_valid, idir, side, 1);
                  // include corner cells if possible by growing toRegion in transverse direction
                  toRegion.grow(1);
                  toRegion.grow(idir, -1);
                  toRegion &= a_state.box();
                  for (BoxIterator bit(toRegion); bit.ok(); ++bit)
                    {
                      IntVect ivTo = bit();                     
                      IntVect ivClose = ivTo - isign*BASISV(idir);
                      for (int ighost=0;ighost<m_nGhosts[0];ighost++,ivTo += isign*BASISV(idir))
                        {
                          //for (int icomp = 0; icomp < a_state.nComp(); icomp++) a_state(ivTo, icomp) = 0.0;
                          IntVect ivFrom = ivClose;
                          
                          // hardwire to linear BCs for now
                          for (int icomp = 0; icomp < a_state.nComp() ; icomp++)
                            {
                              if (m_bcDiri[idir][side][icomp]) 
                                {
                                  a_state(ivTo, icomp) = (-1.0)*a_state(ivFrom, icomp);
                                }
                              else 
                                {
                                  a_state(ivTo, icomp) = (1.0)*a_state(ivFrom, icomp);
                                }
                            }   
                        }
                    } // end loop over cells
                } // if ends match
            } // end loop over sides
        } // if not periodic in this direction
    } // end loop over directions    
}

void
CFStencil::define(
                  const ProblemDomain& a_fineDomain,
                  const Box& a_grid,
                  const DisjointBoxLayout& a_fineBoxes,
                  const DisjointBoxLayout& a_coarBoxes,
                  int a_refRatio,
                  int a_direction,
                  Side::LoHiSide a_hiorlo)
{
  m_isDefined = true;
  CH_assert(a_refRatio >= 1);
  CH_assert(a_direction >= 0);
  CH_assert(a_direction < SpaceDim);
  CH_assert((a_hiorlo == Side::Lo) ||
         (a_hiorlo == Side::Hi));
  CH_assert(!a_fineDomain.isEmpty());

  //set internal vars.  most of these are kept around
  //just to keep the class from having an identity crisis.
  m_direction = a_direction;
  m_hiorlo =  a_hiorlo;

  Box finebox = a_grid;


  //compute intvectset of all points on fine grid that
  //need to be interpolated

  //shift direction
  int hilo = sign(a_hiorlo);

  //create fine stencil
  Box edgebox;
  CH_assert((hilo ==1) || (hilo == -1));
  if (hilo == -1)
    {
      edgebox = adjCellLo(finebox,m_direction,1);
    }
  else
    {
      edgebox = adjCellHi(finebox,m_direction,1);
    }
  edgebox = a_fineDomain & edgebox;
  if (!edgebox.isEmpty())
    {
      Box periodicTestBox(a_fineDomain.domainBox());
      if (a_fineDomain.isPeriodic())
        {
          for (int idir=0; idir<SpaceDim; idir++)
            {
              if (a_fineDomain.isPeriodic(idir))
                {
                  periodicTestBox.grow(idir,-1);
                }
            }
        }

      m_fineIVS.define(edgebox);
      LayoutIterator lit = a_fineBoxes.layoutIterator();
      for (lit.reset(); lit.ok(); ++lit)
        {
          m_fineIVS -= a_fineBoxes[lit()];
          // if periodic, also need to subtract periodic images
          // only do this IF we're periodic _and_ both boxes
          // adjoin the domain box boundary somewhere
          if (a_fineDomain.isPeriodic() && !periodicTestBox.contains(edgebox)
              && !periodicTestBox.contains(a_fineBoxes[lit()]))
            {
              ShiftIterator shiftIt = a_fineDomain.shiftIterator();
              IntVect shiftMult(a_fineDomain.domainBox().size());
              Box shiftedBox(a_fineBoxes[lit()]);
              for (shiftIt.begin(); shiftIt.ok(); ++shiftIt)
                {
                  IntVect shiftVect = shiftMult*shiftIt();
                  shiftedBox.shift(shiftVect);
                  m_fineIVS -= shiftedBox;
                  shiftedBox.shift(-shiftVect);
                } // end loop over periodic shift directions
            } // end if periodic
        }
    }

  //ivs where all coarse slopes are defined
  //== coarsened fine ivs
  m_coarIVS.define(m_fineIVS);
  m_coarIVS.coarsen(a_refRatio);
  // this is a trick to get around the lack of a IntVectSet intersection
  // operator which works with a ProblemDomain
  ProblemDomain coardom= coarsen(a_fineDomain, a_refRatio);
  Box domainIntersectBox = m_coarIVS.minBox();
  domainIntersectBox = coardom & domainIntersectBox;
  m_coarIVS &= domainIntersectBox;

  m_packedBox = m_fineIVS.minBox();
  if (m_fineIVS.numPts() == m_packedBox.numPts())
    {
      m_isPacked = true;
    }
  else
//.........这里部分代码省略.........

void
MappedLevelFluxRegister::define(const DisjointBoxLayout& a_dbl,
                          const DisjointBoxLayout& a_dblCoarse,
                          const ProblemDomain&     a_dProblem,
                          const IntVect&           a_nRefine,
                          int                      a_nComp,
                          bool                     a_scaleFineFluxes)
{
    CH_TIME("MappedLevelFluxRegister::define");
    m_isDefined = FluxRegDefined;  // Basically, define was called
    m_nRefine   = a_nRefine;

    m_scaleFineFluxes = a_scaleFineFluxes;

    DisjointBoxLayout coarsenedFine;
    coarsen(coarsenedFine, a_dbl, a_nRefine);

#ifndef DISABLE_TEMPORARY_FLUX_REGISTER_OPTIMIZATION
    // This doesn't work for multi-block calculations, which are
    // not properly nested. -JNJ

    //begin temporary optimization.  bvs
    int numPts = 0;
    for (LayoutIterator lit = a_dblCoarse.layoutIterator(); lit.ok(); ++lit) {
        numPts += a_dblCoarse[lit].numPts();
    }
    for (LayoutIterator lit = coarsenedFine.layoutIterator(); lit.ok(); ++lit) {
        numPts -= coarsenedFine[lit].numPts();
    }

    if (numPts == 0) {
        m_coarFlux.clear();
        // OK, fine region completely covers coarse region.  no registers.
        return;
    }
#endif

    //end temporary optimization.   bvs
    m_coarFlux.define( a_dblCoarse, a_nComp);
    m_isDefined |= FluxRegCoarseDefined;
    m_domain = a_dProblem;
    ProblemDomain coarsenedDomain;
    coarsen(coarsenedDomain, a_dProblem, a_nRefine);

    m_fineFlux.define( coarsenedFine, a_nComp, IntVect::Unit);
    m_isDefined |= FluxRegFineDefined;

    m_reverseCopier.ghostDefine(coarsenedFine, a_dblCoarse,
                                coarsenedDomain, IntVect::Unit);

    for (int i = 0; i < CH_SPACEDIM; i++) {
        m_coarseLocations[i].define(a_dblCoarse);
        m_coarseLocations[i + CH_SPACEDIM].define(a_dblCoarse);
    }


    DataIterator dC   = a_dblCoarse.dataIterator();
    LayoutIterator dF = coarsenedFine.layoutIterator();

    for (dC.begin(); dC.ok(); ++dC) {
        const Box& cBox = a_dblCoarse.get(dC);

        for (dF.begin(); dF.ok(); ++dF) {
            const Box& fBox  = coarsenedFine.get(dF);

            if (fBox.bigEnd(0) + 1 < cBox.smallEnd(0)) {
                //can skip this box since they cannot intersect, due to sorting
            } else if (fBox.smallEnd(0) - 1 > cBox.bigEnd(0)) {
                //skip to end, since all the rest of boxes will not intersect either
                dF.end();
            } else {

                for (int i = 0; i < CH_SPACEDIM; i++) {
                    Vector<Box>& lo = m_coarseLocations[i][dC];
                    Vector<Box>& hi = m_coarseLocations[i + CH_SPACEDIM][dC];

                    Box loBox = adjCellLo(fBox, i, 1);
                    Box hiBox = adjCellHi(fBox, i, 1);
                    if (cBox.intersectsNotEmpty(loBox)) lo.push_back(loBox & cBox);
                    if (cBox.intersectsNotEmpty(hiBox)) hi.push_back(hiBox & cBox);
                }
            }
        }
    }

    Box domainBox = coarsenedDomain.domainBox();

    if (a_dProblem.isPeriodic()) {
        Vector<Box> periodicBoxes[2 * CH_SPACEDIM];
        for (dF.begin(); dF.ok(); ++dF) {
            const Box& fBox  = coarsenedFine.get(dF);
            for (int i = 0; i < CH_SPACEDIM; i++) {
                if (a_dProblem.isPeriodic(i)) {
                    if (fBox.smallEnd(i) == domainBox.smallEnd(i))
                        periodicBoxes[i].push_back(adjCellLo(fBox, i, 1));
                    if (fBox.bigEnd(i) == domainBox.bigEnd(i))
                        periodicBoxes[i + CH_SPACEDIM].push_back(adjCellHi(fBox, i, 1));
                }
            }
        }
//.........这里部分代码省略.........

void ReductionCopier::define(const BoxLayout& a_level,
                             const BoxLayout& a_dest,
                             const ProblemDomain& a_domain,
                             const IntVect& a_ghost,
                             const Vector<int>& a_transverseDir,
                             bool  a_exchange)
{
#ifdef MULTIDIM_TIMER
  CH_TIME("ReductionCopier::define")
#endif
  m_isDefined = true;
  m_transverseDir = a_transverseDir;

  CH_assert(a_level.isClosed());
  CH_assert(a_dest.isClosed());
  //  CH_assert(a_level.checkPeriodic(a_domain));

  clear();
  buffersAllocated = false;
  //bool self = a_dest == a_level;
  const BoxLayout&         level= a_level;
  const BoxLayout&         dest = a_dest;

  // note that while much of the implementation for the
  // periodic case remains in this function (as copied, more
  // or less, from the original Copier), at least at the moment
  // periodic support is "turned off" for the ReductionCopier
  // because it's not entirely clear it will do the right thing
  // for this case.  I've left the vast majority of the
  // periodic implementation in place, however, in case we
  // change our minds... (DFM 3/10/09)

  // set up vector of dataIndexes to keep track of which
  // "to" boxes are not completely contained within the primary
  // domain.  these boxes are then candidates for filling by
  // periodic images of the "from" data.
  Vector<DataIndex> periodicallyFilledToVect;

  // in order to cull which "from" data may be needed to
  // fill the "to" data, keep track of the radius around the
  // primary domain in which all these cells lie.
  // do this by incrementally growing the domain box and
  // keeping track of what this radius is.
  // just to make things simpler, start off with a radius of one
  Box grownDomainCheckBox = a_domain.domainBox();
  grownDomainCheckBox.grow(1);
  int periodicCheckRadius = 1;

  // since valid regions of the "from" DBL may also be outside
  // the primary domain, need to keep track of whether any of these
  // need to be checked separately.
  Vector<DataIndex> periodicFromVect;
  // use same domain trick here as well
  Box grownFromDomainCheckBox = a_domain.domainBox();
  int periodicFromCheckRadius = 1;

  Box domainBox(a_domain.domainBox());
  // grab index extents of domain in transverse direction
  Vector<int> transverseLo(m_transverseDir.size());
  Vector<int> transverseHi(m_transverseDir.size());

  for (int n=0; n<m_transverseDir.size(); n++)
    {
      // grab index extents of domain in transverse direction
      transverseLo[n] = domainBox.smallEnd(m_transverseDir[n]);
      transverseHi[n] = domainBox.bigEnd(m_transverseDir[n]);
    }

  bool isPeriodic = false;
  // I think the right thing here is that the ReductionCopier
  // completely ignores periodicity
  //if (!domainBox.isEmpty())
  //isPeriodic = a_domain.isPeriodic();

  // (dfm -- 9/13/05) as currently written, the Copier won't correctly
  // handle periodic cases where the number of ghost cells is greater
  // than the width of the domain.  We _should_ do multiple wraparounds,
  // but we don't. So, put in this assertion. We can revisit this if it
  // becomes an issue
  if (isPeriodic)
    {
      for (int dir=0; dir<SpaceDim; dir++)
        {
          if (a_domain.isPeriodic(dir))
            {
              CH_assert (a_ghost[dir] <= domainBox.size(dir));
            }
        }
    }

  unsigned int myprocID = procID();

  // The following 4 for loops are the result of a performance optimization.
  // When increasing the size of the problem, we found that the code was
  // looping over every destination box for every source box which was N1*N2
  // loop iterations (essentially an N-squared approach).
  // The following code attempts to simply reduce N1 and N2 by first separating
  // the boxes (or LayoutIndexes to boxes) that reside on the current processor.
  // Then the loop to determine which boxes of the first list intersect with
  // which boxes of the second list can be done in N1' * N2' iterations,
//.........这里部分代码省略.........

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