C++ MultiFab::define方法代码示例

void
MCLinOp::makeCoefficients (MultiFab&       cs,
                           const MultiFab& fn,
                           int             level)
{
    const int nc = fn.nComp();
    //
    // Determine index type of incoming MultiFab.
    //
    const IndexType iType(fn.boxArray().ixType());
    const IndexType cType(D_DECL(IndexType::CELL, IndexType::CELL, IndexType::CELL));
    const IndexType xType(D_DECL(IndexType::NODE, IndexType::CELL, IndexType::CELL));
    const IndexType yType(D_DECL(IndexType::CELL, IndexType::NODE, IndexType::CELL));
#if (BL_SPACEDIM == 3)    
    const IndexType zType(D_DECL(IndexType::CELL, IndexType::CELL, IndexType::NODE));
#endif
    int cdir;
    if (iType == cType)
    {
        cdir = -1;
    }
    else if (iType == xType)
    {
        cdir = 0;
    }
    else if (iType == yType)
    {
        cdir = 1;
    }
#if (BL_SPACEDIM == 3)
    else if (iType == zType)
    {
        cdir = 2;
    }
#endif
    else
        BoxLib::Abort("MCLinOp::makeCoeffients(): Bad index type");
    
    BoxArray d(gbox[level]);
    if (cdir >= 0)
	d.surroundingNodes(cdir);

    int nGrow=0;
    cs.define(d, nc, nGrow, Fab_allocate);
    cs.setVal(0.0);

    const BoxArray& grids = gbox[level];

    for (MFIter csmfi(cs); csmfi.isValid(); ++csmfi)
    {
        const Box&       grd   = grids[csmfi.index()];
        FArrayBox&       csfab = cs[csmfi];
        const FArrayBox& fnfab = fn[csmfi];

	switch(cdir)
        {
	case -1:
	    FORT_AVERAGECC(
		csfab.dataPtr(),
                ARLIM(csfab.loVect()), ARLIM(csfab.hiVect()),
		fnfab.dataPtr(),
                ARLIM(fnfab.loVect()), ARLIM(fnfab.hiVect()),
		grd.loVect(),
                grd.hiVect(), &nc);
	    break;
	case 0:
	case 1:
	case 2:
	    if ( harmavg )
            {
		FORT_HARMONIC_AVERAGEEC(
		    csfab.dataPtr(), 
                    ARLIM(csfab.loVect()), ARLIM(csfab.hiVect()),
		    fnfab.dataPtr(), 
                    ARLIM(fnfab.loVect()), ARLIM(fnfab.hiVect()),
		    grd.loVect(),
                    grd.hiVect(), &nc, &cdir);
	    }
            else
            {
		FORT_AVERAGEEC(
		    csfab.dataPtr(), 
                    ARLIM(csfab.loVect()), ARLIM(csfab.hiVect()),
		    fnfab.dataPtr(), 
                    ARLIM(fnfab.loVect()), ARLIM(fnfab.hiVect()),
		    grd.loVect(),
                    grd.hiVect(), &nc, &cdir);
	    }
	    break;
	default:
	    BoxLib::Error("MCLinOp::makeCoeffients(): bad coefficient coarsening direction!");
	}
    }
}

int main(int argc, char* argv[])
{
  BoxLib::Initialize(argc,argv);

  BL_PROFILE_VAR("main()", pmain);

  std::cout << std::setprecision(15);

  solver_type = BoxLib_C;
  bc_type = Periodic;

  Real     a = 0.0;
  Real     b = 1.0;

  // ---- First use the number of processors to decide how many grids you have.
  // ---- We arbitrarily decide to have one grid per MPI process in a uniform
  // ---- cubic domain, so we require that the number of processors be N^3.
  // ---- This requirement is somewhat arbitrary, but convenient for now.

  int nprocs = ParallelDescriptor::NProcs();

  // N is the cube root of the number of processors
  int N(0);
  for(int i(1); i*i*i <= nprocs; ++i) {
    if(i*i*i == nprocs) {
      N = i;
    }
  }

  if(N == 0) {  // not a cube
    if(ParallelDescriptor::IOProcessor()) {
      std::cerr << "**** Error:  nprocs = " << nprocs << " is not currently supported." << std::endl;
    }
    BoxLib::Error("We require that the number of processors be a perfect cube");
  }
  if(ParallelDescriptor::IOProcessor()) {
    std::cout << "N = " << N << std::endl;
  }


  // ---- make a box, then a boxarray with maxSize
  int domain_hi = (N*maxGrid) - 1;
  Box domain(IntVect(0,0,0), IntVect(domain_hi,domain_hi,domain_hi));
  BoxArray bs(domain);
  bs.maxSize(maxGrid);

  // This defines the physical size of the box.  Right now the box is [0,1] in each direction.
  RealBox real_box;
  for (int n = 0; n < BL_SPACEDIM; n++) {
    real_box.setLo(n, 0.0);
    real_box.setHi(n, 1.0);
  }
 
  // This says we are using Cartesian coordinates
  int coord = 0;
  
  // This sets the boundary conditions to be periodic or not
  int is_per[BL_SPACEDIM];
  
  if (bc_type == Dirichlet || bc_type == Neumann) {
    for (int n = 0; n < BL_SPACEDIM; n++) is_per[n] = 0;
  } 
  else {
    for (int n = 0; n < BL_SPACEDIM; n++) is_per[n] = 1;
  }
 
  // This defines a Geometry object which is useful for writing the plotfiles
  Geometry geom(domain,&real_box,coord,is_per);

  for ( int n=0; n<BL_SPACEDIM; n++ ) {
    dx[n] = ( geom.ProbHi(n) - geom.ProbLo(n) )/domain.length(n);
  }

  if (ParallelDescriptor::IOProcessor()) {
     std::cout << "Domain size     : " << N << std::endl;
     std::cout << "Max_grid_size   : " << maxGrid << std::endl;
     std::cout << "Number of grids : " << bs.size() << std::endl;
  }

  // Allocate and define the right hand side.
  MultiFab rhs(bs, Ncomp, 0, Fab_allocate); 
  setup_rhs(rhs, geom, a, b);

  MultiFab alpha(bs, Ncomp, 0, Fab_allocate);
  MultiFab beta[BL_SPACEDIM];
  for ( int n=0; n<BL_SPACEDIM; ++n ) {
    BoxArray bx(bs);
    beta[n].define(bx.surroundingNodes(n), Ncomp, 1, Fab_allocate);
  }

  setup_coeffs(bs, alpha, beta, geom);

  MultiFab anaSoln;
  if (comp_norm) {
    anaSoln.define(bs, Ncomp, 0, Fab_allocate);
    compute_analyticSolution(anaSoln);
  }

  // Allocate the solution array 
  // Set the number of ghost cells in the solution array.
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    std::cout << "Grid resolution : " << n_cell << " (cells)" << std::endl;
    std::cout << "Domain size     : " << real_box.hi(0) - real_box.lo(0) << " (length unit) " << std::endl;
    std::cout << "Max_grid_size   : " << max_grid_size << " (cells)" << std::endl;
    std::cout << "Number of grids : " << bs.size() << std::endl;
  }

  // Allocate and define the right hand side.
  bool do_4th = (solver_type==BoxLib_C4 || solver_type==All);
  int ngr = (do_4th ? 1 : 0);
  MultiFab rhs(bs, Ncomp, ngr); 
  setup_rhs(rhs, geom);

  // Set up the Helmholtz operator coefficients.
  MultiFab alpha(bs, Ncomp, 0);
  PArray<MultiFab> beta(BL_SPACEDIM, PArrayManage);
  for ( int n=0; n<BL_SPACEDIM; ++n ) {
    BoxArray bx(bs);
    beta.set(n, new MultiFab(bx.surroundingNodes(n), Ncomp, 0, Fab_allocate));
  }

  // The way HPGMG stores face-centered data is completely different than the
  // way BoxLib does it, and translating between the two directly via indexing
  // magic is a nightmare. Happily, the way this tutorial calculates
  // face-centered values is by first calculating cell-centered values and then
  // interpolating to the cell faces. HPGMG can do the same thing, so rather
  // than converting directly from BoxLib's face-centered data to HPGMG's, just
  // give HPGMG the cell-centered data and let it interpolate itself.

  MultiFab beta_cc(bs,Ncomp,1); // cell-centered beta
  setup_coeffs(bs, alpha, beta, geom, beta_cc);

  MultiFab alpha4, beta4;
  if (do_4th) {
    alpha4.define(bs, Ncomp, 4, Fab_allocate);
    beta4.define(bs, Ncomp, 3, Fab_allocate);
    setup_coeffs4(bs, alpha4, beta4, geom);
  }

  MultiFab anaSoln;
  if (comp_norm || plot_err || plot_asol) {
    anaSoln.define(bs, Ncomp, 0, Fab_allocate);
    compute_analyticSolution(anaSoln,Array<Real>(BL_SPACEDIM,0.5));
    
    if (plot_asol) {
      writePlotFile("ASOL", anaSoln, geom);
    }
  }

  // Allocate the solution array 
  // Set the number of ghost cells in the solution array.
  MultiFab soln(bs, Ncomp, 1);
  MultiFab soln4;
  if (do_4th) {
    soln4.define(bs, Ncomp, 3, Fab_allocate);
  }
  MultiFab gphi(bs, BL_SPACEDIM, 0);

#ifdef USEHYPRE
  if (solver_type == Hypre || solver_type == All) {
    if (ParallelDescriptor::IOProcessor()) {
      std::cout << "----------------------------------------" << std::endl;
      std::cout << "Solving with Hypre " << std::endl;
    }

    solve(soln, anaSoln, gphi, a, b, alpha, beta, beta_cc, rhs, bs, geom, Hypre);
  }

void
LinOp::makeCoefficients (MultiFab&       cs,
                         const MultiFab& fn,
                         int             level)
{
    BL_PROFILE("LinOp::makeCoefficients()");

    int nc = 1;
    //
    // Determine index type of incoming MultiFab.
    //
    const IndexType iType(fn.boxArray().ixType());
    const IndexType cType(D_DECL(IndexType::CELL, IndexType::CELL, IndexType::CELL));
    const IndexType xType(D_DECL(IndexType::NODE, IndexType::CELL, IndexType::CELL));
    const IndexType yType(D_DECL(IndexType::CELL, IndexType::NODE, IndexType::CELL));
#if (BL_SPACEDIM == 3)    
    const IndexType zType(D_DECL(IndexType::CELL, IndexType::CELL, IndexType::NODE));
#endif

    int cdir;
    if (iType == cType)
    {
        cdir = -1;
    }
    else if (iType == xType)
    {
        cdir = 0;
    }
    else if (iType == yType)
    {
        cdir = 1;
#if (BL_SPACEDIM == 3)
    }
    else if (iType == zType)
    {
        cdir = 2;
#endif    
    }
    else
    {
        BoxLib::Error("LinOp::makeCoeffients: Bad index type");
    }

    BoxArray d(gbox[level]);
    if (cdir >= 0)
        d.surroundingNodes(cdir);
    //
    // Only single-component solves supported (verified) by this class.
    //
    const int nComp=1;
    const int nGrow=0;
    cs.define(d, nComp, nGrow, Fab_allocate);

    const bool tiling = true;

    switch (cdir)
    {
    case -1:
#ifdef _OPENMP
#pragma omp parallel
#endif
        for (MFIter csmfi(cs,tiling); csmfi.isValid(); ++csmfi)
        {
            const Box& tbx = csmfi.tilebox();
            FArrayBox&       csfab = cs[csmfi];
            const FArrayBox& fnfab = fn[csmfi];

            FORT_AVERAGECC(csfab.dataPtr(), ARLIM(csfab.loVect()),
                           ARLIM(csfab.hiVect()),fnfab.dataPtr(),
                           ARLIM(fnfab.loVect()),ARLIM(fnfab.hiVect()),
                           tbx.loVect(),tbx.hiVect(), &nc);
        }
        break;
    case 0:
    case 1:
    case 2:
        if (harmavg)
        {
#ifdef _OPENMP
#pragma omp parallel
#endif
  	    for (MFIter csmfi(cs,tiling); csmfi.isValid(); ++csmfi)
            {
	        const Box& tbx = csmfi.tilebox();
                FArrayBox&       csfab = cs[csmfi];
                const FArrayBox& fnfab = fn[csmfi];

                FORT_HARMONIC_AVERAGEEC(csfab.dataPtr(),
                                        ARLIM(csfab.loVect()),
                                        ARLIM(csfab.hiVect()),
                                        fnfab.dataPtr(),
                                        ARLIM(fnfab.loVect()),
                                        ARLIM(fnfab.hiVect()),
                                        tbx.loVect(),tbx.hiVect(),
                                        &nc,&cdir);
            }
        }
        else
        {
#ifdef _OPENMP
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    // Create coarse boundary register, fill w/data from coarse FAB
    int bndry_InRad=0;
    int bndry_OutRad=1;
    int bndry_Extent=1;
    BoxArray cbs = BoxArray(bs).coarsen(ratio);
    BndryRegister cbr(cbs,bndry_InRad,bndry_OutRad,bndry_Extent,Ncomp);
    for (OrientationIter face; face; ++face)
    {
	Orientation f = face();
	FabSet& bnd_fs(cbr[f]);
	bnd_fs.copyFrom(crse_mf, 0, 0, 0, Ncomp);
    }
  
    // Interpolate crse data to fine boundary, where applicable
    int cbr_Nstart=0;
    int fine_Nstart=0;
    int bndry_Nstart=0;
    vbd.setBndryValues(cbr,cbr_Nstart,fine,fine_Nstart,
		       bndry_Nstart,Ncomp,ratio,pbcarray);
  
    Nghost = 1; // other variables don't need extra space
    
    DivVis lp(vbd,H);
    
    Real a = 0.0;
    Real b[BL_SPACEDIM];
    b[0] = 1.0;
    b[1] = 1.0;
#if BL_SPACEDIM>2
    b[2] = 1.0;
#endif
    MultiFab  acoefs;
    int NcompA = (BL_SPACEDIM == 2  ?  2  :  1);
    acoefs.define(bs, NcompA, Nghost, Fab_allocate);
    acoefs.setVal(a);
    MultiFab bcoefs[BL_SPACEDIM];
    for (n=0; n<BL_SPACEDIM; ++n)
    {
	BoxArray bsC(bs);
	bcoefs[n].define(bsC.surroundingNodes(n), 1,
			 Nghost, Fab_allocate);
#if 1
	for(MFIter bmfi(bcoefs[n]); bmfi.isValid(); ++bmfi)
	{
	    FORT_MAKEMU(bcoefs[n][bmfi].dataPtr(),
			ARLIM(bcoefs[n][bmfi].loVect()),ARLIM(bcoefs[n][bmfi].hiVect()),H,n);
	}
#else
	bcoefs[n].setVal(b[n]);
#endif
    } // -->> over dimension
    lp.setCoefficients(acoefs, bcoefs);
#if 1
    lp.maxOrder(4);
#endif
    
    Nghost = 1;
    MultiFab tsoln(bs, Ncomp, Nghost, Fab_allocate); 
    tsoln.setVal(0.0);
#if 1
    tsoln.copy(fine);
#endif
#if 0
    // testing apply
    lp.apply(out,tsoln);
    Box subbox = out[0].box();