C++ ProblemDomain类代码示例

NewPoissonOp* NewPoissonOpFactory::MGnewOp(const ProblemDomain& a_FineindexSpace,
                                           int   a_depth,
                                           bool  a_homoOnly)
{
  CH_assert(a_depth >= 0 );
  CH_assert(m_bc != NULL);
  NewPoissonOp* newOp = new NewPoissonOp();
  RealVect dx = m_dx;
  ProblemDomain domain = a_FineindexSpace;
  for (int i=0; i<a_depth; i++)
    {
      Box d = domain.domainBox();
      d.coarsen(8);
      d.refine(8);
      if (domain.domainBox() == d)
        {
          dx*=2;
          domain.coarsen(2);
        }
      else
        {
          return NULL;
        }
    }
  newOp->define(dx, domain, m_bc);
  return newOp;
}

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 mtac::LiveVariableAnalysisProblem::meet(ProblemDomain& out, const ProblemDomain& in){
    if(out.top()){
        out = in;
        return;
    } else if(in.top()){
        //out does not change
        return;
    }

    for(auto& value : in.values()){
        out.values().insert(value);
    }
}

// Define new AMR level
void AMRLevelPluto::define(AMRLevel*            a_coarserLevelPtr,
                           const ProblemDomain& a_problemDomain,
                           int                  a_level,
                           int                  a_refRatio)
{
  if (s_verbosity >= 3)
  {
    pout() << "AMRLevelPluto::define " << a_level << endl;
  }

  // Call inherited define
  AMRLevel::define(a_coarserLevelPtr,
                   a_problemDomain,
                   a_level,
                   a_refRatio);

  // Get setup information from the next coarser level
  if (a_coarserLevelPtr != NULL)
  {
    AMRLevelPluto* amrGodPtr = dynamic_cast<AMRLevelPluto*>(a_coarserLevelPtr);

    if (amrGodPtr != NULL)
    {
      m_cfl = amrGodPtr->m_cfl;
      m_domainLength = amrGodPtr->m_domainLength;
      m_refineThresh = amrGodPtr->m_refineThresh;
      m_tagBufferSize = amrGodPtr->m_tagBufferSize;
    }
    else
    {
      MayDay::Error("AMRLevelPluto::define: a_coarserLevelPtr is not castable to AMRLevelPluto*");
    }
  }

  // Compute the grid spacing
  m_dx = m_domainLength / (a_problemDomain.domainBox().bigEnd(0)-a_problemDomain.domainBox().smallEnd(0)+1.);

  // Nominally, one layer of ghost cells is maintained permanently and
  // individual computations may create local data with more
  m_numGhost = 1;

  CH_assert(m_patchPluto != NULL);
  CH_assert(isDefined());
  m_patchPluto->define(m_problem_domain,m_dx,m_level,m_numGhost);

  // Get additional information from the patch integrator
  m_numStates      = m_patchPluto->numConserved();
  m_ConsStateNames = m_patchPluto->ConsStateNames();
  m_PrimStateNames = m_patchPluto->PrimStateNames();
}

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);,

/** Initializes the DisjointBoxLayouts, each step. */
void
EBRestart::makeHierarchy(Vector<DisjointBoxLayout>& a_dbl,
                         const ProblemDomain&       a_baseDomain,
                         const IntVectSet&          a_baseTags,
                         const Vector<int>&         a_refRatios,
                         const InputParams&         a_inputs)
{
  a_dbl.resize(a_inputs.nlevs);
  Real fillRatio     = 0.85;
  BRMeshRefine regridder(a_baseDomain, a_refRatios, fillRatio,
                         a_inputs.blocking_factor,
                         a_inputs.buffer_size, a_inputs.max_size);

  Vector<Vector<Box> > oldGrids(a_inputs.nlevs,1), newGrids(a_inputs.nlevs);
  oldGrids[0][0]=a_baseDomain.domainBox();
  for ( int l=1; l<a_inputs.nlevs; ++l )
  {
    oldGrids[l][0]=coarsen(oldGrids[l-1][0], a_refRatios[l-1]);
    // FIXME: is a_refRatios[l-1] correct??
  }

  regridder.regrid(newGrids, a_baseTags, 0, 1, oldGrids);
  // 0 is baselevel, 1 is toplevel.

  Vector<int> procs;
  for (int l=0; l<a_inputs.nlevs; l++)
  {
    newGrids[l].sort();
    LoadBalance(procs, newGrids[l]);
    a_dbl[l] = DisjointBoxLayout(newGrids[l], procs);
  }
}

void mtac::LiveVariableAnalysisProblem::transfer(mtac::basic_block_p/* basic_block*/, mtac::Quadruple& quadruple, ProblemDomain& in){
    if(in.top()){
        ProblemDomain::Values values;
        in.int_values = values;
    }

    if(quadruple.op != mtac::Operator::NOP){
        if(mtac::erase_result(quadruple.op)){
            in.values().erase(quadruple.result);
        } else {
            in.values().insert(quadruple.result);
        }

        if_init<std::shared_ptr<Variable>>(quadruple.arg1, [&in](std::shared_ptr<Variable>& var){in.values().insert(var);});
        if_init<std::shared_ptr<Variable>>(quadruple.arg2, [&in](std::shared_ptr<Variable>& var){in.values().insert(var);});
    }
}

ProblemDomain mtac::LiveVariableAnalysisProblem::meet(ProblemDomain& out, ProblemDomain& in){
    if(out.top()){
        return in;
    } else if(in.top()){
        return out;
    }

    typename ProblemDomain::Values values;
    ProblemDomain result(values);

    for(auto& value : in.values()){
        result.values().insert(value);
    }
    
    for(auto& value : out.values()){
        result.values().insert(value);
    }

    return result;
}

void mtac::global_cse::meet(ProblemDomain& in, const ProblemDomain& out){
    eddic_assert(!in.top() || !out.top(), "At least one lattice should not be a top element");

    if(in.top()){
        in = out;
    } else if(out.top()){
        //in does not change
    } else {
        auto& first = in.values();
        auto& second = out.values();
        
        std::set<expression> intersection;

        std::set_intersection(first.begin(), first.end(), second.begin(), second.end(), std::inserter(intersection, intersection.begin()));

        in.values() = std::move(intersection);
    }
}

int VCAMRPoissonOp2Factory::refToFiner(const ProblemDomain& a_domain) const
{
  int retval = -1;
  bool found = false;

  for (int ilev = 0; ilev < m_domains.size(); ilev++)
    {
      if (m_domains[ilev].domainBox() == a_domain.domainBox())
        {
          retval = m_refRatios[ilev];
          found = true;
        }
    }

  if (!found)
    {
      MayDay::Abort("Domain not found in AMR hierarchy");
    }

  return retval;
}

void mtac::global_cse::transfer(mtac::basic_block_p basic_block, ProblemDomain& out){
    auto& out_values = out.values();
    auto it = out_values.begin();

    //Compute AEin - Kill(i)

    while(it != out_values.end()){
        if(Kill[basic_block].find(*it) != Kill[basic_block].end()){
            it = out_values.erase(it);
            continue;
        }

        ++it;
    }

    //Compute Eval(i) U (AEin - Kill(i))

    for(auto& expression : Eval[basic_block]){
        out_values.insert(expression);
    }
}

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
CFStencil::define(
                  const ProblemDomain& a_fineDomain,
                  const Box& a_grid,
                  const Vector<Box>& a_periodicVector,
                  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()) return;


  int w1 = edgebox.smallEnd()[0];
  int w2 = edgebox.bigEnd()[0];
  m_fineIVS.define(edgebox);
  // moving window loop in i-direction (bvs)
  for (int i=0; i<a_periodicVector.size(); ++i)
    {
      const Box& b = a_periodicVector[i];
      if (b.bigEnd()[0] >= w1)
        {
          m_fineIVS -= b;
          if (b.smallEnd()[0] > w2)
          {
            i=a_periodicVector.size();
          }
        }

    }


  //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
    {
      m_isPacked = false;
      m_packedBox = Box();
    }
}

// ---------------------------------------------------------
// 7 Dec 2005
Real
norm(const BoxLayoutData<NodeFArrayBox>& a_layout,
     const LevelData<NodeFArrayBox>& a_mask,
     const ProblemDomain& a_domain,
     const Real a_dx,
     const int a_p,
     const Interval& a_interval,
     bool a_verbose)
{
  if (a_p == 0)
    return maxnorm(a_layout, a_mask, a_domain, a_interval, a_verbose);

  Real normTotal = 0.;
  int ncomp = a_interval.size();
  Box domBox = a_domain.domainBox();
  for (DataIterator it = a_layout.dataIterator(); it.ok(); ++it)
    {
      const NodeFArrayBox& thisNfab = a_layout[it()];
      const FArrayBox& dataFab = thisNfab.getFab();
      const FArrayBox& maskFab = a_mask[it()].getFab();
      const Box& thisBox(a_layout.box(it())); // CELL-centered
      NodeFArrayBox dataMasked(thisBox, ncomp);
      FArrayBox& dataMaskedFab = dataMasked.getFab();
      dataMaskedFab.copy(dataFab);
      // dataMaskedFab *= maskFab;
      for (int comp = a_interval.begin(); comp <= a_interval.end(); comp++)
        {
          // Set dataMaskedFab[comp] *= maskFab[0].
          dataMaskedFab.mult(maskFab, 0, comp);
        }

      Real thisNfabNorm = 0.;
      if (thisBox.intersects(domBox))
        {
          Box thisBoxInDomain = thisBox & domBox;
          thisNfabNorm = norm(dataMasked, a_dx, thisBoxInDomain, a_p,
                              a_interval.begin(), a_interval.size());
        }

      if (a_verbose)
        cout << a_p << "norm(" << thisBox << ") = " << thisNfabNorm << endl;

      if (a_p == 1)
        {
          normTotal += thisNfabNorm;
        }
      else if (a_p == 2)
        {
          normTotal += thisNfabNorm * thisNfabNorm;
        }
      else
        {
          normTotal += pow(thisNfabNorm, Real(a_p));
        }
    }
# ifdef CH_MPI
  Real recv;
  // add up (a_p is not 0)
  int result = MPI_Allreduce(&normTotal, &recv, 1, MPI_CH_REAL,
                             MPI_SUM, Chombo_MPI::comm);
  if (result != MPI_SUCCESS)
    { //bark!!!
      MayDay::Error("sorry, but I had a communication error on norm");
    }
  normTotal = recv;
# endif
  // now do sqrt, etc
  if (a_p == 2)
    normTotal = sqrt(normTotal);
  else
    if ((a_p != 0) && (a_p != 1))
      normTotal = pow(normTotal, (Real)1.0/Real(a_p));

  return normTotal;
}

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