C++ BoxArray::intersections方法代码示例

int
iMultiFab::norm0 (int comp, const BoxArray& ba, bool local) const
{
    int nm0 = -std::numeric_limits<int>::max();

#ifdef _OPENMP
#pragma omp parallel
#endif
    {
	std::vector< std::pair<int,Box> > isects;
	int priv_nm0 = -std::numeric_limits<int>::max();

	for (MFIter mfi(*this); mfi.isValid(); ++mfi)
	{
	    ba.intersections(mfi.validbox(),isects);
	    
	    for (int i = 0, N = isects.size(); i < N; i++)
	    {
		priv_nm0 = std::max(priv_nm0, get(mfi).norm(isects[i].second, 0, comp, 1));
	    }
	}
#ifdef _OPENMP
#pragma omp critical (imultifab_norm0_ba)
#endif
	{
	    nm0 = std::max(nm0, priv_nm0);
	}
    }
 
    if (!local)
	ParallelDescriptor::ReduceIntMax(nm0);
 
    return nm0;
}

MultiFab*
Nyx::build_fine_mask()
{
    BL_ASSERT(level > 0); // because we are building a mask for the coarser level

    if (fine_mask != 0) return fine_mask;

    BoxArray baf = parent->boxArray(level);
    baf.coarsen(crse_ratio);

    const BoxArray& bac = parent->boxArray(level-1);
    fine_mask = new MultiFab(bac,parent->DistributionMap(level-1), 1,0);
    fine_mask->setVal(1.0);

#ifdef _OPENMP
#pragma omp parallel
#endif
    for (MFIter mfi(*fine_mask); mfi.isValid(); ++mfi)
    {
        FArrayBox& fab = (*fine_mask)[mfi];

        std::vector< std::pair<int,Box> > isects = baf.intersections(fab.box());

        for (int ii = 0; ii < isects.size(); ii++)
        {
            fab.setVal(0.0,isects[ii].second,0);
        }
    }

    return fine_mask;
}

//.........这里部分代码省略.........
    typedef FabArray<NodeFab> MultiNodeFab; 
    PArray<MultiNodeFab> nodes(Nlev,PArrayManage);

    std::cerr << "Before nodes allocated" << endl;
    for (int lev=0; lev<Nlev; ++lev)
        nodes.set(lev,new MultiNodeFab(gridArray[lev],1,nGrow));
    std::cerr << "After nodes allocated" << endl;

    int cnt = 0;
    typedef std::map<Node,int> NodeMap;
    NodeMap nodeMap;
    for (int lev=0; lev<Nlev; ++lev)
    {
        for (MFIter fai(nodes[lev]); fai.isValid(); ++fai)
        {
            NodeFab& ifab = nodes[lev][fai];
            const Box& box = ifab.box() & subboxArray[lev];
            for (IntVect iv=box.smallEnd(); iv<=box.bigEnd(); box.next(iv))
                ifab(iv,0) = Node(iv,lev,fai.index(),Node::VALID);
        }
            
        if (lev != 0)
        {
            const int ref = amrData.RefRatio()[lev-1];
            const Box& rangeBox = Box(IntVect::TheZeroVector(),
                                     (ref-1)*IntVect::TheUnitVector());

            BoxArray bndryCells = GetBndryCells(nodes[lev].boxArray(),ref,geom[lev]);

            for (MFIter fai(nodes[lev]); fai.isValid(); ++fai)
            {
                const Box& box = BoxLib::grow(fai.validbox(),ref) & subboxArray[lev];
                NodeFab& ifab = nodes[lev][fai];
                std::vector< std::pair<int,Box> > isects = bndryCells.intersections(box);
                for (int i = 0; i < isects.size(); i++)
                {
                    Box co = isects[i].second & fai.validbox();
                    if (co.ok())
                        std::cout << "BAD ISECTS: " << co << std::endl;

                    const Box& dstBox = isects[i].second;
                    const Box& srcBox = BoxLib::coarsen(dstBox,ref);

                    NodeFab dst(dstBox,1);
                    for (IntVect iv(srcBox.smallEnd());
                         iv<=srcBox.bigEnd();
                         srcBox.next(iv))
                    {
                        const IntVect& baseIV = ref*iv;
                        for (IntVect ivt(rangeBox.smallEnd());ivt<=rangeBox.bigEnd();rangeBox.next(ivt))
                            dst(baseIV + ivt,0) = Node(iv,lev-1,-1,Node::VALID);
                    }
                    const Box& ovlp = dstBox & ifab.box();

                    Box mo = ovlp & fai.validbox();
                    if (mo.ok())
                    {
                        std::cout << "BAD OVERLAP: " << mo << std::endl;
                        std::cout << "         vb: " << fai.validbox() << std::endl;
                    }
                    if (ovlp.ok())
                        ifab.copy(dst,ovlp,0,ovlp,0,1);
                }
            }
        }

//.........这里部分代码省略.........
        std::cout << "plotVarNames " << plotVarNames[comps[i]] << std::endl;;
        destFillComps[i] = i;
    }

    const int nGrow = 0;

    if (ParallelDescriptor::IOProcessor() && verbose>0)
    {
        cerr << "Will read the following states: ";
        for (int i=0; i<nComp; ++i)
            cerr << " " << amrData.StateNumber(inVarNames[i]) << " (" << inVarNames[i] << ")" ;
        cerr << '\n';
    }

    const Box& probDomain = amrData.ProbDomain()[finestLevel];
    int dir=BL_SPACEDIM-1; pp.query("dir",dir);
    const IntVect lo=probDomain.smallEnd();
    IntVect hi=lo; hi[dir] = probDomain.bigEnd()[dir];
    const Box resBox(lo,hi);
    FArrayBox resFab(resBox,nComp); resFab.setVal(0.);

    int accumFac = 1;
    for (int lev=finestLevel; lev>=0; --lev)
    {
        const BoxArray& ba = amrData.boxArray(lev);
	const DistributionMapping& dm = amrData.DistributionMap(lev);
        MultiFab mf(ba,dm,nComp,nGrow);
        if (ParallelDescriptor::IOProcessor() && verbose>0)
            cerr << "...filling data at level " << lev << endl;
        amrData.FillVar(mf,lev,inVarNames,destFillComps);

        // Zero covered regions
        if (lev < finestLevel)
        {
            const BoxArray baf = BoxArray(amrData.boxArray(lev+1)).coarsen(amrData.RefRatio()[lev]);

	    for (MFIter mfi(mf); mfi.isValid(); ++mfi)
            {
		FArrayBox& myFab = mf[mfi];
                std::vector< std::pair<int,Box> > isects = baf.intersections(ba[mfi.index()]);
                
		for (int ii = 0; ii < isects.size(); ii++)
		    myFab.setVal(0,isects[ii].second,0,nComp);
            }
        }

        // Scale by "volume" (plane, 1 cell thick) of cells, normalized to volume of finest cells
        mf.mult(accumFac*accumFac,0,nComp);

        for (MFIter mfi(mf); mfi.isValid(); ++mfi)
        {
            const FArrayBox& fab = mf[mfi];
            const Box& box = mfi.validbox();
            IntVect ivlo = box.smallEnd();
            IntVect ivhi = box.bigEnd(); ivhi[dir] = ivlo[dir];

            for (int plane=box.smallEnd(dir); plane<=box.bigEnd(dir); ++plane)
            {
                ivlo[dir] = plane;
                ivhi[dir] = plane;

                Box subbox(ivlo,ivhi);
                Vector<Real> thisSum(nComp);
                for (int n=0; n<nComp; ++n)
                    thisSum[n] = fab.sum(subbox,n,1);

                // Now, increment each sum affected by this coarse plane
                for (int r=0; r<accumFac; ++r)
                {
                    const int finePlane = plane*accumFac + r;
                    IntVect ivF=resBox.smallEnd(); ivF[dir] = finePlane;
                    for (int n=0; n<nComp; ++n)
                        resFab(ivF,n) += thisSum[n];
                }
            }
        }

        if (lev>0)
            accumFac *= amrData.RefRatio()[lev-1];
    }

    // Accumulate sums from all processors to IOProc
    ParallelDescriptor::ReduceRealSum(resFab.dataPtr(),nComp*resBox.numPts(),ParallelDescriptor::IOProcessorNumber());

    if (ParallelDescriptor::IOProcessor())
    {
        // Scale result by total "volume" (plane, 1 cell thick) of slab at each plane
        Real volume=1;
        for (int d=0; d<BL_SPACEDIM; ++d)
            if (d!=dir)
                volume*=probDomain.length(d);

        resFab.mult(1/volume);
        std::cout << "RESFAB " << resFab << std::endl;
    }

    
    amrex::Finalize();
    return 0;
}

static
BoxArray
GetBndryCells (const BoxArray& ba,
               int             ngrow,
               const Geometry& geom)
{
    //
    // First get list of all ghost cells.
    //
    BoxList gcells, bcells;

    for (int i = 0; i < ba.size(); ++i)
	gcells.join(BoxLib::boxDiff(BoxLib::grow(ba[i],ngrow),ba[i]));
    //
    // Now strip out intersections with original BoxArray.
    //
    for (BoxList::const_iterator it = gcells.begin(); it != gcells.end(); ++it)
    {
        std::vector< std::pair<int,Box> > isects = ba.intersections(*it);

        if (isects.empty())
            bcells.push_back(*it);
        else
        {
            //
            // Collect all the intersection pieces.
            //
            BoxList pieces;
            for (int i = 0; i < isects.size(); i++)
                pieces.push_back(isects[i].second);
            BoxList leftover = BoxLib::complementIn(*it,pieces);
            bcells.catenate(leftover);
        }
    }
    //
    // Now strip out overlaps.
    //
    gcells.clear();
    gcells = BoxLib::removeOverlap(bcells);
    bcells.clear();

    if (geom.isAnyPeriodic())
    {
        Array<IntVect> pshifts(27);

        const Box& domain = geom.Domain();

        for (BoxList::const_iterator it = gcells.begin(); it != gcells.end(); ++it)
        {
            if (!domain.contains(*it))
            {
                //
                // Add in periodic ghost cells shifted to valid region.
                //
                geom.periodicShift(domain, *it, pshifts);

                for (int i = 0; i < pshifts.size(); i++)
                {
                    const Box& shftbox = *it + pshifts[i];

                    const Box& ovlp = domain & shftbox;
                    BoxList bl = BoxLib::complementIn(ovlp,BoxList(ba));
                    bcells.catenate(bl);
                }
            }
        }

        gcells.catenate(bcells);
    }

    return BoxArray(gcells);
}