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

/* ************************************************************************* */
void PatchPluto::computeRefGradient(FArrayBox& gFab, FArrayBox& UFab, const Box& b)
/*!
 * Compute numerical gradient of the solution in order to tag
 * zones for refinement.
 * The gradient is computed by standard finite differences using
 *
 * - REF_CRIT equal to 1 --> compute (normalized) gradient using 1st
 *                           derivative of the solution;
 * - REF_CRIT equal to 2 --> compute (normalized) gradient using 2nd
 *                           derivative of the solution (default);
 *                           This approach is based on Lohner (1987).
 *
 * Zones will be flagged for refinement whenever grad[k][j][i] exceeds
 * the threshold value specified by the 'Refine_thresh' parameter read in
 * pluto.ini.
 *
 * Derivatives are computed using the conserved variable U[REF_VAR]
 * where REF_VAR is taken to be energy density (default).
 * However, by setting REF_VAR = -1, you can provide your own
 * physical variable through the function computeRefVar().
 *
 * \authors C. Zanni   ([email protected])\n
 *          A. Mignone ([email protected])
 * \date    Oct 11, 2012
 *************************************************************************** */
{
    CH_assert(m_isDefined);

    int nv, i, j, k;

    double rp, rm, r, tp, tm;
    double x, dqx_p, dqx_m, dqx, d2qx, den_x;
    double y, dqy_p, dqy_m, dqy, d2qy, den_y;
    double z, dqz_p, dqz_m, dqz, d2qz, den_z;

    double gr1, gr2, eps = 0.01;
    double ***UU[NVAR], ***q, ***grad;
    RBox  Ubox, Gbox;

    rp = rm = r = 1.0;
    tp = tm = 1.0;

    /* -- check ref criterion -- */

#if REF_CRIT != 1 && REF_CRIT != 2
    print ("! TagCells.cpp: Refinement criterion not valid\n");
    QUIT_PLUTO(1);
#endif

    /* -----------------------------------------------
        The solution array U is defined on the box
        [Uib, Uie] x [Ujb, Uje] x [Ukb, Uke], which
        differs from that of gFab ([Gib,...Gke]),
        typically one point larger in each direction.
       ----------------------------------------------- */

    Ubox.jb = Ubox.je = Ubox.kb = Ubox.ke = 0;
    Gbox.jb = Gbox.je = Gbox.kb = Gbox.ke = 0;

    D_EXPAND(Ubox.ib = UFab.loVect()[IDIR]; Ubox.ie = UFab.hiVect()[IDIR]; ,

void
CellBilinear::interp (const FArrayBox&  crse,
                      int               crse_comp,
                      FArrayBox&        fine,
                      int               fine_comp,
                      int               ncomp,
                      const Box&        fine_region,
                      const IntVect &   ratio,
                      const Geometry& /*crse_geom*/,
                      const Geometry& /*fine_geom*/,
                      Array<BCRec>&   /*bcr*/,
                      int               actual_comp,
                      int               actual_state)
{
    BL_PROFILE("CellBilinear::interp()");
#if (BL_SPACEDIM == 3)
    BoxLib::Error("interp: not implemented");
#endif
    //
    // Set up to call FORTRAN.
    //
    const int* clo = crse.box().loVect();
    const int* chi = crse.box().hiVect();
    const int* flo = fine.loVect();
    const int* fhi = fine.hiVect();
    const int* lo  = fine_region.loVect();
    const int* hi  = fine_region.hiVect();
    int num_slope  = D_TERM(2,*2,*2)-1;
    int len0       = crse.box().length(0);
    int slp_len    = num_slope*len0;

    Array<Real> slope(slp_len);

    int strp_len = len0*ratio[0];

    Array<Real> strip(strp_len);

    int strip_lo = ratio[0] * clo[0];
    int strip_hi = ratio[0] * chi[0];

    const Real* cdat  = crse.dataPtr(crse_comp);
    Real*       fdat  = fine.dataPtr(fine_comp);
    const int* ratioV = ratio.getVect();

    FORT_CBINTERP (cdat,ARLIM(clo),ARLIM(chi),ARLIM(clo),ARLIM(chi),
                   fdat,ARLIM(flo),ARLIM(fhi),ARLIM(lo),ARLIM(hi),
                   D_DECL(&ratioV[0],&ratioV[1],&ratioV[2]),&ncomp,
                   slope.dataPtr(),&num_slope,strip.dataPtr(),&strip_lo,&strip_hi,
                   &actual_comp,&actual_state);
}

/* ************************************************************************* */
void PatchPluto::computeRefGradient(FArrayBox& gFab, FArrayBox& UFab, const Box& b)
/*
 *
 * PURPOSE
 *
 *   Tag cells for refinement by computing grad[k][j][i]. 
 *   By default a convex combination of the first and second
 *   derivative of the total energy density is used.
 *   alpha = 0 --> triggers refinement towards the 2nd derivative
 *   alpha = 1 --> triggers refinement towards the 1st derivative
 *
 * 
 *
 *************************************************************************** */
{
  CH_assert(m_isDefined);

  int nv, i, j, k;
  int Uib, Uie, Ujb=0, Uje=0, Ukb=0, Uke=0;
  int Gib, Gie, Gjb=0, Gje=0, Gkb=0, Gke=0;

  double rp, rm, r;
  double x, dqx_p, dqx_m, d2qx, den_x;
  double y, dqy_p, dqy_m, d2qy, den_y;
  double z, dqz_p, dqz_m, d2qz, den_z;

  double alpha, qref, gr1, gr2;

  double eps = 0.01;
  double ***UU[NVAR], ***q, ***grad;

  double us[NVAR], vs[NVAR], mu;
  static double **T;

  rp = rm = r = 1.0;

/* -----------------------------------------------
    The solution array U is defined on the box 
    [Uib, Uie] x [Ujb, Uje] x [Ukb, Uke], which 
    differs from that of gFab ([Gib,...Gke]), 
    typically one point larger in each direction. 
   ----------------------------------------------- */
    
  D_EXPAND(Uib = UFab.loVect()[IDIR]; Uie = UFab.hiVect()[IDIR]; ,

void
NodeBilinear::interp (const FArrayBox&  crse,
                      int               crse_comp,
                      FArrayBox&        fine,
                      int               fine_comp,
                      int               ncomp,
                      const Box&        fine_region,
                      const IntVect&    ratio,
                      const Geometry& /*crse_geom */,
                      const Geometry& /*fine_geom */,
                      Array<BCRec>&   /*bcr*/,
                      int               actual_comp,
                      int               actual_state)
{
    BL_PROFILE("NodeBilinear::interp()");
    //
    // Set up to call FORTRAN.
    //
    const int* clo = crse.box().loVect();
    const int* chi = crse.box().hiVect();
    const int* flo = fine.loVect();
    const int* fhi = fine.hiVect();
    const int* lo  = fine_region.loVect();
    const int* hi  = fine_region.hiVect();
    int num_slope  = D_TERM(2,*2,*2)-1;
    int len0       = crse.box().length(0);
    int slp_len    = num_slope*len0;

    Array<Real> strip(slp_len);

    const Real* cdat  = crse.dataPtr(crse_comp);
    Real*       fdat  = fine.dataPtr(fine_comp);
    const int* ratioV = ratio.getVect();

    FORT_NBINTERP (cdat,ARLIM(clo),ARLIM(chi),ARLIM(clo),ARLIM(chi),
                   fdat,ARLIM(flo),ARLIM(fhi),ARLIM(lo),ARLIM(hi),
                   D_DECL(&ratioV[0],&ratioV[1],&ratioV[2]),&ncomp,
                   strip.dataPtr(),&num_slope,&actual_comp,&actual_state);
}

void
CellConservativeLinear::interp (const FArrayBox& crse,
                                int              crse_comp,
                                FArrayBox&       fine,
                                int              fine_comp,
                                int              ncomp,
                                const Box&       fine_region,
                                const IntVect&   ratio,
                                const Geometry&  crse_geom,
                                const Geometry&  fine_geom,
                                Array<BCRec>&    bcr,
                                int              actual_comp,
                                int              actual_state)
{
    BL_PROFILE("CellConservativeLinear::interp()");
    BL_ASSERT(bcr.size() >= ncomp);

    //
    // Make box which is intersection of fine_region and domain of fine.
    //
    Box target_fine_region = fine_region & fine.box();
    //
    // crse_bx is coarsening of target_fine_region, grown by 1.
    //
    Box crse_bx = CoarseBox(target_fine_region,ratio);
    //
    // Slopes are needed only on coarsening of target_fine_region.
    //
    Box cslope_bx(crse_bx);
    cslope_bx.grow(-1);
    //
    // Make a refinement of cslope_bx
    //
    Box fine_version_of_cslope_bx = BoxLib::refine(cslope_bx,ratio);
    //
    // Get coarse and fine edge-centered volume coordinates.
    //
    Array<Real> fvc[BL_SPACEDIM];
    Array<Real> cvc[BL_SPACEDIM];
    int dir;
    for (dir = 0; dir < BL_SPACEDIM; dir++)
    {
        fine_geom.GetEdgeVolCoord(fvc[dir],fine_version_of_cslope_bx,dir);
        crse_geom.GetEdgeVolCoord(cvc[dir],crse_bx,dir);
    }
    //
    // alloc tmp space for slope calc.
    //
    // In ucc_slopes and lcc_slopes , there is a slight abuse of
    // the number of compenents argument
    // --> there is a slope for each component in each coordinate
    //     direction
    //
    FArrayBox ucc_slopes(cslope_bx,ncomp*BL_SPACEDIM);
    FArrayBox lcc_slopes(cslope_bx,ncomp*BL_SPACEDIM);
    FArrayBox slope_factors(cslope_bx,BL_SPACEDIM);

    FArrayBox  cmax(cslope_bx,ncomp);
    FArrayBox  cmin(cslope_bx,ncomp);
    FArrayBox alpha(cslope_bx,ncomp);

    Real* fdat       = fine.dataPtr(fine_comp);
    const Real* cdat = crse.dataPtr(crse_comp);
    Real* ucc_xsldat = ucc_slopes.dataPtr(0);
    Real* lcc_xsldat = lcc_slopes.dataPtr(0);
    Real* xslfac_dat = slope_factors.dataPtr(0);
#if (BL_SPACEDIM>=2)
    Real* ucc_ysldat = ucc_slopes.dataPtr(ncomp);
    Real* lcc_ysldat = lcc_slopes.dataPtr(ncomp);
    Real* yslfac_dat = slope_factors.dataPtr(1);
#endif
#if (BL_SPACEDIM==3)
    Real* ucc_zsldat = ucc_slopes.dataPtr(2*ncomp);
    Real* lcc_zsldat = lcc_slopes.dataPtr(2*ncomp);
    Real* zslfac_dat = slope_factors.dataPtr(2);
#endif

    const int* flo    = fine.loVect();
    const int* fhi    = fine.hiVect();
    const int* clo    = crse.loVect();
    const int* chi    = crse.hiVect();
    const int* fblo   = target_fine_region.loVect();
    const int* fbhi   = target_fine_region.hiVect();
    const int* csbhi  = cslope_bx.hiVect();
    const int* csblo  = cslope_bx.loVect();
    int lin_limit     = (do_linear_limiting ? 1 : 0);
    const int* cvcblo = crse_bx.loVect();
    const int* fvcblo = fine_version_of_cslope_bx.loVect();
    int slope_flag    = 1;

    int cvcbhi[BL_SPACEDIM];
    int fvcbhi[BL_SPACEDIM];

    for (dir=0; dir<BL_SPACEDIM; dir++)
    {
        cvcbhi[dir] = cvcblo[dir] + cvc[dir].size() - 1;
        fvcbhi[dir] = fvcblo[dir] + fvc[dir].size() - 1;
    }

    D_TERM(Real* voffx = new Real[fvc[0].size()]; ,