C++ Box类代码示例

void VCAMRPoissonOp2::reflux(const LevelData<FArrayBox>&        a_phiFine,
                            const LevelData<FArrayBox>&        a_phi,
                            LevelData<FArrayBox>&              a_residual,
                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp)
{
  CH_TIME("VCAMRPoissonOp2::reflux");

  int ncomp = 1;
  ProblemDomain fineDomain = refine(m_domain, m_refToFiner);
  LevelFluxRegister levfluxreg(a_phiFine.disjointBoxLayout(),
                               a_phi.disjointBoxLayout(),
                               fineDomain,
                               m_refToFiner,
                               ncomp);

  levfluxreg.setToZero();
  Interval interv(0,a_phi.nComp()-1);

  DataIterator dit = a_phi.dataIterator();
  for (dit.reset(); dit.ok(); ++dit)
    {
      const FArrayBox& coarfab = a_phi[dit];
      const FluxBox& coarBCoef  = (*m_bCoef)[dit];
      const Box& gridBox = a_phi.getBoxes()[dit];

      for (int idir = 0; idir < SpaceDim; idir++)
        {
          FArrayBox coarflux;
          Box faceBox = surroundingNodes(gridBox, idir);
          getFlux(coarflux, coarfab, coarBCoef , faceBox, idir);

          Real scale = 1.0;
          levfluxreg.incrementCoarse(coarflux, scale,dit(),
                                     interv,interv,idir);
        }
    }
  LevelData<FArrayBox>& p = ( LevelData<FArrayBox>&)a_phiFine;

  // has to be its own object because the finer operator
  // owns an interpolator and we have no way of getting to it
  VCAMRPoissonOp2* finerAMRPOp = (VCAMRPoissonOp2*) a_finerOp;
  QuadCFInterp& quadCFI = finerAMRPOp->m_interpWithCoarser;

  quadCFI.coarseFineInterp(p, a_phi);
  // p.exchange(a_phiFine.interval()); // BVS is pretty sure this is not necesary.
  IntVect phiGhost = p.ghostVect();

  DataIterator ditf = a_phiFine.dataIterator();
  const  DisjointBoxLayout& dblFine = a_phiFine.disjointBoxLayout();
  for (ditf.reset(); ditf.ok(); ++ditf)
    {
      const FArrayBox& phifFab = a_phiFine[ditf];
      const FluxBox& fineBCoef  = (*(finerAMRPOp->m_bCoef))[ditf];
      const Box& gridbox = dblFine.get(ditf());
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          int normalGhost = phiGhost[idir];
          SideIterator sit;
          for (sit.begin(); sit.ok(); sit.next())
            {
              Side::LoHiSide hiorlo = sit();
              Box fabbox;
              Box facebox;

              // assumption here that the stencil required
              // to compute the flux in the normal direction
              // is 2* the number of ghost cells for phi
              // (which is a reasonable assumption, and probably
              // better than just assuming you need one cell on
              // either side of the interface
              // (dfm 8-4-06)
              if (sit() == Side::Lo)
                {
                  fabbox = adjCellLo(gridbox,idir, 2*normalGhost);
                  fabbox.shift(idir, 1);
                  facebox = bdryLo(gridbox, idir,1);
                }
              else
                {
                  fabbox = adjCellHi(gridbox,idir, 2*normalGhost);
                  fabbox.shift(idir, -1);
                  facebox = bdryHi(gridbox, idir, 1);
                }

              // just in case we need ghost cells in the transverse direction
              // (dfm 8-4-06)
              for (int otherDir=0; otherDir<SpaceDim; ++otherDir)
                {
                  if (otherDir != idir)
                    {
                      fabbox.grow(otherDir, phiGhost[otherDir]);
                    }
                }
              CH_assert(!fabbox.isEmpty());

              FArrayBox phifab(fabbox, a_phi.nComp());
              phifab.copy(phifFab);

              FArrayBox fineflux;
              getFlux(fineflux, phifab, fineBCoef, facebox, idir,
//.........这里部分代码省略.........

void CollisionTest::display()
{
	const static GLfloat lightPosition[] = {1, -1, 0, 0};
	const static GLfloat lightPositionMirror[] = { 1, 1, 0, 0 };

	// Update the transform matrices of each box in turn.
	for (Box *box = boxData; box < boxData + boxes; box++)
	{
		box->calculateInternals();
		box->isOverlapping = false;
	}

	// Update the transform matrices of each ball in turn.
	for (Ball *ball = ballData; ball < ballData + balls; ball++)
	{
		// Run the physics.
		ball->calculateInternals();
	}

	// Clear the viewport and set the camera direction.
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	gluLookAt(18.0f, 0, 0, 0, 0, 0, 0, 1.0f, 0);
	glRotatef(-phi, 0, 0, 1);
	glRotatef(theta, 0, 1, 0);
	glTranslatef(0, -5.0f, 0);

	// Render each element in turn as a shadow.
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_LIGHTING);
	glLightfv(GL_LIGHT0, GL_POSITION, lightPosition);
	glEnable(GL_LIGHT0);
	glPushMatrix();
	glMultMatrixf(floorMirror);
	glColorMaterial(GL_FRONT_AND_BACK, GL_DIFFUSE);
	glEnable(GL_COLOR_MATERIAL);
	for (Box *box = boxData; box < boxData + boxes; box++)
	{
		box->render();
	}
	for (Ball *ball = ballData; ball < ballData + balls; ball++)
	{
		ball->render();
	}
	glPopMatrix();
	glDisable(GL_COLOR_MATERIAL);
	glDisable(GL_LIGHTING);
	glDisable(GL_LIGHT0);

	DrawGrid(Vector3(35, 0, 30), 60);

	// Render each shadow in turn
	glEnable(GL_BLEND);
	glColor4f(0, 0, 0, 0.1f);
	glDisable(GL_DEPTH_TEST);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	for (Box *box = boxData; box < boxData + boxes; box++)
	{
		box->renderShadow();
	}
	for (Ball *ball = ballData; ball < ballData + balls; ball++)
	{
		ball->renderShadow();
	}
	glDisable(GL_BLEND);
	
	// Render the boxes themselves
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_LIGHTING);
	glLightfv(GL_LIGHT0, GL_POSITION, lightPositionMirror);
	glColorMaterial(GL_FRONT_AND_BACK, GL_DIFFUSE);
	glEnable(GL_LIGHT0);
	glEnable(GL_COLOR_MATERIAL);
	for (Box *box = boxData; box < boxData + boxes; box++)
	{
		box->render();
	}
	for (Ball *ball = ballData; ball < ballData + balls; ball++)
	{
		ball->render();
	}
	glDisable(GL_COLOR_MATERIAL);
	glDisable(GL_LIGHTING);
	glDisable(GL_LIGHT0);
	
	drawDebug();
}

void
LevelFluxRegisterEdge::refluxCurl(LevelData<FluxBox>& a_uCoarse,
                                  Real a_scale)
{
  CH_assert(isDefined());
  CH_assert(a_uCoarse.nComp() == m_nComp);

  SideIterator side;
  // idir is the normal direction to the coarse-fine interface
  for (int idir=0 ; idir<SpaceDim; ++idir)
  {
    for (side.begin(); side.ok(); ++side)
    {
      LevelData<FluxBox>& fineReg = m_fabFine[index(idir, side())];

      // first, create temp LevelData<FluxBox> to hold "coarse flux"
      const DisjointBoxLayout coarseBoxes = m_regCoarse.getBoxes();

      // this fills the place of what used to be m_fabCoarse in the old
      // implementation
      LevelData<FluxBox> coarReg(coarseBoxes, m_nComp, IntVect::Unit);


      // now fill the coarReg with the curl of the stored coarse-level
      // edge-centered flux
      DataIterator crseDit = coarseBoxes.dataIterator();
      for (crseDit.begin(); crseDit.ok(); ++crseDit)
        {
          FluxBox& thisCoarReg = coarReg[crseDit];
          thisCoarReg.setVal(0.0);

          EdgeDataBox& thisEdgeData = m_regCoarse[crseDit];

          for (int edgeDir=0; edgeDir<SpaceDim; edgeDir++)
            {
              if (idir != edgeDir)
                {
                  FArrayBox& crseEdgeDataDir = thisEdgeData[edgeDir];
                  for (int faceDir = 0; faceDir<SpaceDim; faceDir++)
                    {
                      if (faceDir != edgeDir)
                        {
                          FArrayBox& faceData = thisCoarReg[faceDir];
                          int shiftDir = -1;
                          for (int i=0; i<SpaceDim; i++)
                            {
                              if ((i != faceDir) && (i != edgeDir) )
                                {
                                  shiftDir = i;
                                }
                            }
                          CH_assert(shiftDir >= 0);
                          crseEdgeDataDir.shiftHalf(shiftDir, sign(side()));
                          // scaling already taken care of in incrementCrse
                          Real scale = 1.0;
                          faceData.plus(crseEdgeDataDir, scale, 0, 0, faceData.nComp());
                          crseEdgeDataDir.shiftHalf(shiftDir, -sign(side()));
                        } // end if not normal direction
                    } // end loop over face directions
                } // end if edgeDir != idir
            } // end loop over edge directions
        } // end loop over crse boxes


      // first, we need to create a temp LevelData<FluxBox>
      // to make a local copy in the coarse layout space of
      // the fine register increments

      LevelData<FluxBox> fineRegLocal(coarReg.getBoxes(), m_nComp, IntVect::Unit);

      fineReg.copyTo(fineReg.interval(), fineRegLocal,
                     fineRegLocal.interval(),
                     m_crseCopiers[index(idir,side())]);

      for (DataIterator it = a_uCoarse.dataIterator(); it.ok(); ++it)
        {
          // loop over flux components here
          for (int fluxComp=0; fluxComp < SpaceDim; fluxComp++)
            {
              // we don't do anything in the normal direction
              if (fluxComp != idir)
                {
                  // fluxDir is the direction of the face-centered flux
                  FArrayBox& U = a_uCoarse[it()][fluxComp];
                  // set up IntVectSet to avoid double counting of updates
                  Box coarseGridBox = U.box();
                  // transfer to Cell-centered, then create IVS
                  coarseGridBox.shiftHalf(fluxComp,1);
                  IntVectSet nonUpdatedEdges(coarseGridBox);

                  // remember, we want to take the curl here
                  // also recall that fluxComp is the component
                  // of the face-centered curl (not the edge-centered
                  // vector field that we're refluxing, which is why
                  // the sign may seem like it's the opposite of what
                  // you might expect!
                  Real local_scale = -sign(side())*a_scale;
                  //int testDir = (fluxComp+1)%(SpaceDim);
                  if (((fluxComp+1)%(SpaceDim)) == idir)  local_scale *= -1;
                  Vector<IntVectSet>& ivsV =
//.........这里部分代码省略.........

  // Computes   omega = A*psi
  void execute( const vector<complex>& psi, vector<complex>& omega )
  {
#if STAGE_TIMES
    StopWatch timer;
    timer.start();
#endif

    int nLevels = tree.numLevels();

    // Compute lowest level multipole
    Level<DIM>& leafLevel = getLevel(nLevels);
    leafLevel.zeroFields();
    for( BoxIter bi = leafLevel.boxbegin(); bi != leafLevel.boxend(); ++bi ) {
      Box* b = *bi;

      const Vec3& center = b->center;

      const Box::pointIter piEnd = b->pointIndex.end();
      for( Box::pointIter pi = b->pointIndex.begin(); pi != piEnd; ++pi ) {
        int index = *pi;
        const Vec3 r = center - sourcefield[index];
        cout << "Accumulating " << index << " into Box " << b->n << endl;
        Translation_Function::add( r, psi[index], b->getMultipole() );
      }
    }

#if STAGE_TIMES
    cout << "Lowest Level Multipoles Done " << timer.stop() << endl;
#endif

    // Upward Pass - Compute the multipoles
    for( int L = nLevels; L >= 3; --L ) {

      Level<DIM>& level = getLevel(L);
      Level<DIM>& pLevel = getLevel(L-1);
      pLevel.zeroFields();
      NFunction& Ms = pLevel.getScratch();

      // For all boxes at level L
      for( BoxIter bi = level.boxbegin(); bi != level.boxend(); ++bi ) {
        Box* b = *bi;

        // Interpolate to the parent level L-1
        level.getInterp().apply( b->getMultipole(), Ms );

        // Multiply by the transfer function and accumulate into parent
        b->parent->getMultipole().addProduct( Ms, pLevel.getTranslationUp(b) );
      }
    }

#if STAGE_TIMES
    cout << "Upward Pass Done " << timer.stop() << endl;
#endif

    // Do the Transfers
    for( int L = nLevels; L >= 2; --L ) {
      getLevel(L).applyTransferFunctions();
    }

#if STAGE_TIMES
    cout << "Transfer List Applied " << timer.stop() << endl;
#endif

    for( int L = 3; L <= nLevels; ++L ) {

      Level<DIM>& level = getLevel(L);
      Level<DIM>& pLevel = getLevel(L-1);
      NFunction& Lb = level.getScratch();
      NFunction& LB = pLevel.getScratch();

      // For all boxes at level L
      for( BoxIter bi = level.boxbegin(); bi != level.boxend(); ++bi ) {
        Box* b = *bi;

        // Multiply parent's local expansion by translation function
        LB.setProduct( b->parent->getLocal(), pLevel.getTranslationDown(b) );

        // Anterpolate
        pLevel.getAnterp().apply( LB, Lb );

        // Accumulate into the box's local expansion
        b->getLocal().add( Lb );
      }
    }

#if STAGE_TIMES
    cout << "Downward Pass Done " << timer.stop() << endl;
#endif

    // Zero omega
    omega.assign( omega.size(), 0 );

    // Integrate and add close interactions

    // For all the boxes at the lowest level, integrate local
    NFunction& LE = leafLevel.getScratch();           // Scratch space
    for( BoxIter bi = leafLevel.boxbegin(); bi != leafLevel.boxend(); ++bi ) {
      Box* b = *bi;
      Vec3& center = b->center;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	glPopMatrix();
	glPushMatrix();
	u3.drawstrikerbox(0.58f,7.0f,-4.0f,0.0f,-8.0f,0.0f,0.0f,0.0f);
	glPopMatrix();
	//stricker box holes
	glPushMatrix();
	hole[0].drawhole(3.5f,4.0f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[1].drawhole(-4.0f,3.5f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[2].drawhole(-3.5f,-4.0f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[3].drawhole(4.0f,-3.5f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[4].drawhole(3.5f,-4.0f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[5].drawhole(4.0f,3.5f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[6].drawhole(-3.5f,4.0f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	glPushMatrix();
	hole[7].drawhole(-4.0f,-3.5f,-8.0f,1.0f,0.0f,0.0f,0.3f);
	glPopMatrix();
	//coins
	glPushMatrix();
	coin[0].drawcoin(coin[0].a,coin[0].c,-8.0f,0.0f,0.0f,0.0f,coin[0].rad);
	glPopMatrix();
	glPushMatrix();
	coin[1].drawcoin(coin[1].a,coin[1].c,-8.0f,0.0f,0.0f,0.0f,coin[1].rad);
	glPopMatrix();
	glPushMatrix();
	coin[2].drawcoin(coin[2].a,coin[2].c,-8.0f,0.0f,0.0f,0.0f,coin[2].rad);
	glPopMatrix();
	glPushMatrix();
	coin[3].drawcoin(coin[3].a,coin[3].c,-8.0f,0.0f,0.0f,0.0f,coin[3].rad);
	glPopMatrix();
	glPushMatrix();
	coin[4].drawcoin(coin[4].a,coin[4].c,-8.0f,0.8f,0.4f,0.0f,coin[4].rad);
	glPopMatrix();
	glPushMatrix();
	coin[5].drawcoin(coin[5].a,coin[5].c,-8.0f,0.8f,0.4f,0.0f,coin[5].rad);
	glPopMatrix();
	glPushMatrix();
	coin[6].drawcoin(coin[6].a,coin[6].c,-8.0f,0.8f,0.4f,0.0f,coin[6].rad);
	glPopMatrix();
	glPushMatrix();
	coin[7].drawcoin(coin[7].a,coin[7].c,-8.0f,0.8f,0.4f,0.0f,coin[7].rad);
	glPopMatrix();
	glPushMatrix();
	coin[8].drawcoin(coin[8].a,coin[8].c,-8.0f,1.0f,0.0f,0.0f,coin[8].rad);
	glPopMatrix();
	// striker
	glPushMatrix();
	s[0].drawstriker(s[0].a,s[0].c,-8.0f,0.0f,1.0f,0.0f,s[0].rad);
	glPopMatrix();
	//striker pointer
	glPushMatrix();
	glTranslatef(s[0].a,s[0].c,-8.0f);
	glColor3f(0.0f,0.0f,0.0f);
	glRotatef(theta,0.0f,0.0f,-1.0f);
	glBegin(GL_LINES);
	glVertex2f(0.0f,0.0f);
	glVertex2f(lx,ly);
	glEnd();
	glPopMatrix();
	// speedometer
	glPushMatrix();
	glBegin(GL_QUADS);
	glColor3f(0.0,1.0,0.0);
	glVertex2f(3.8,-3.0);
	glVertex2f(4.0,-3.0);
	glVertex2f(4.0,pointer);
	glVertex2f(3.8,pointer);
	glColor3f(1.0,0.0,0.1);
	glVertex2f(4.0,2.0);
	glVertex2f(3.8,2.0);
	glVertex2f(3.8,pointer);
	glVertex2f(4.0,pointer);
	glEnd();
	glColor3f(0.4,0.3,0.1);
	glVertex2f(0.48,pointer);
	glVertex2f(0.55,pointer);
	glPopMatrix();
	//score-board
	scoreBoard.drawBox();
	glColor3f(0.0f,0.0f,1.0f);
	char string[] = "Score";
	scoreBoard.drawText(string,strlen(string),-0.7f,0.2f);
	ostringstream sc;
	sc<<scoreBoard.score1;
	score1 = sc.str();
	scoreBoard.drawText(score1.data(),score1.size(),-0.2f,-0.4f); 
	glutSwapBuffers();
}

KOKKOS_INLINE_FUNCTION
bool isValid( Box const &b )
{
    return isValid( b.minCorner() ) && isValid( b.maxCorner() );
}

		Mesh* Sphere::GenerateMesh()
		{
			Box box;
			return box.GenerateMesh();
		}

  void MeshFromSpline2D (SplineGeometry2d & geometry,
			 Mesh *& mesh, 
			 MeshingParameters & mp)
  {
    PrintMessage (1, "Generate Mesh from spline geometry");

    double h = mp.maxh;

    Box<2> bbox = geometry.GetBoundingBox ();

    if (bbox.Diam() < h) 
      {
	h = bbox.Diam();
	mp.maxh = h;
      }

    mesh = new Mesh;
    mesh->SetDimension (2);

    geometry.PartitionBoundary (h, *mesh);

    // marks mesh points for hp-refinement
    for (int i = 0; i < geometry.GetNP(); i++)
      if (geometry.GetPoint(i).hpref)
	{
	  double mindist = 1e99;
	  PointIndex mpi(0);
	  Point<2> gp = geometry.GetPoint(i);
	  Point<3> gp3(gp(0), gp(1), 0);
	  for (PointIndex pi = PointIndex::BASE; 
	       pi < mesh->GetNP()+PointIndex::BASE; pi++)
	    if (Dist2(gp3, (*mesh)[pi]) < mindist)
	      {
		mpi = pi;
		mindist = Dist2(gp3, (*mesh)[pi]);
	      }
	  (*mesh)[mpi].Singularity(1.);
	}


    int maxdomnr = 0;
    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      {
	if ( (*mesh)[si].domin > maxdomnr) maxdomnr = (*mesh)[si].domin;
	if ( (*mesh)[si].domout > maxdomnr) maxdomnr = (*mesh)[si].domout;
      }

    mesh->ClearFaceDescriptors();
    for (int i = 1; i <= maxdomnr; i++)
      mesh->AddFaceDescriptor (FaceDescriptor (i, 0, 0, i));

    // set Array<string*> bcnames... 
    // number of bcnames
    int maxsegmentindex = 0;
    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      {
	if ( (*mesh)[si].si > maxsegmentindex) maxsegmentindex = (*mesh)[si].si;
      }

    mesh->SetNBCNames(maxsegmentindex);

    for ( int sindex = 0; sindex < maxsegmentindex; sindex++ )
      mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) );

    for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
      (*mesh)[si].SetBCName ( (*mesh).GetBCNamePtr( (*mesh)[si].si-1 ) );

    Point3d pmin(bbox.PMin()(0), bbox.PMin()(1), -bbox.Diam());
    Point3d pmax(bbox.PMax()(0), bbox.PMax()(1), bbox.Diam());

    mesh->SetLocalH (pmin, pmax, mparam.grading);
    mesh->SetGlobalH (h);
  
    mesh->CalcLocalH();

    int bnp = mesh->GetNP(); // boundary points

    int hquad = mparam.quad;


    for (int domnr = 1; domnr <= maxdomnr; domnr++)
      if (geometry.GetDomainTensorMeshing (domnr))
        { // tensor product mesh
          
          Array<PointIndex, PointIndex::BASE> nextpi(bnp);
          Array<int, PointIndex::BASE> si1(bnp), si2(bnp);
          PointIndex firstpi;
          
          nextpi = -1;
          si1 = -1;
          si2 = -1;
          for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++)
            {
              int p1 = -1, p2 = -2;

              if ( (*mesh)[si].domin == domnr)
                { p1 = (*mesh)[si][0]; p2 = (*mesh)[si][1]; }
              if ( (*mesh)[si].domout == domnr)
                { p1 = (*mesh)[si][1]; p2 = (*mesh)[si][0]; }
              
//.........这里部分代码省略.........

void
Moof::ParseTraf(Box& aBox, Trex& aTrex, Mvhd& aMvhd, Mdhd& aMdhd, Edts& aEdts, Sinf& aSinf, uint64_t* aDecodeTime, bool aIsAudio)
{
  MOZ_ASSERT(aDecodeTime);
  Tfhd tfhd(aTrex);
  Tfdt tfdt;
  for (Box box = aBox.FirstChild(); box.IsAvailable(); box = box.Next()) {
    if (box.IsType("tfhd")) {
      tfhd = Tfhd(box, aTrex);
    } else if (!aTrex.mTrackId || tfhd.mTrackId == aTrex.mTrackId) {
      if (box.IsType("tfdt")) {
        tfdt = Tfdt(box);
      } else if (box.IsType("saiz")) {
        mSaizs.AppendElement(Saiz(box, aSinf.mDefaultEncryptionType));
      } else if (box.IsType("saio")) {
        mSaios.AppendElement(Saio(box, aSinf.mDefaultEncryptionType));
      }
    }
  }
  if (aTrex.mTrackId && tfhd.mTrackId != aTrex.mTrackId) {
    return;
  }
  // Now search for TRUN boxes.
  uint64_t decodeTime =
    tfdt.IsValid() ? tfdt.mBaseMediaDecodeTime : *aDecodeTime;
  for (Box box = aBox.FirstChild(); box.IsAvailable(); box = box.Next()) {
    if (box.IsType("trun")) {
      if (ParseTrun(box, tfhd, aMvhd, aMdhd, aEdts, &decodeTime, aIsAudio)) {
        mValid = true;
      } else {
        mValid = false;
        break;
      }
    }
  }
  *aDecodeTime = decodeTime;
}

//.........这里部分代码省略.........
    exactStep = nStep*crseMult;

    ostrstream exactFile;
    ostrstream computedFile;
    ostrstream errorFile;

    exactFile.fill('0');
    computedFile.fill('0');
    errorFile.fill('0');

    if (isTimeDep)
    {
      constructPlotFileName(exactFile, exactRoot, intFieldSize, exactStep);
      constructPlotFileName(computedFile, computedRoot, intFieldSize, nStep);
      constructPlotFileName(errorFile, errorRoot, intFieldSize, nStep);

    }
    else
    {
      // if not time dependent, file roots are really filenames
      exactFile << exactRoot << ends;
      computedFile << computedRoot << ends;
      errorFile << errorRoot << ends;
    }

    pout() << "exact Filename = " << exactFile.str() << endl;
    pout() << "computed Filename = " << computedFile.str() << endl;
    if (doPlots)
    {
      pout() << "error Filename = " << errorFile.str() << endl;
    }

    // declare memory
    Vector<LevelData<FArrayBox>* > exactSoln;
    Vector<string> exactVars; // exact solution variable names
    Vector<DisjointBoxLayout> exactGrids;
    Box exactDomain;
    Real exactDx, exactDt, exactTime;
    Vector<int> exactRefRatio;
    int exactNumLevels;
    IntVect ghostVect = IntVect::Unit;
    string exactFileName(exactFile.str());

    // get exact solution
    if (verbose)
    {
      pout() << "read exact solution..." << endl;
    }

    ReadAMRHierarchyHDF5(exactFileName,
                         exactGrids,
                         exactSoln,
                         exactVars,
                         exactDomain,
                         exactDx,
                         exactDt,
                         exactTime,
                         exactRefRatio,
                         exactNumLevels);

    if (verbose)
    {
      pout () << "done reading exact soln" << endl;
    }

    // we assume that exact soln is single-grid if we're doing averaging

// this function averages down the fine solution to the valid
// regions of the computed solution, then subtracts ir from
// the computed solution.  (error is exact-computed)
void computeAMRError(Vector<LevelData<FArrayBox>* >&       a_error,
                     const Vector<string>&                 a_errorVars,
                     const Vector<LevelData<FArrayBox>* >& a_computedSoln,
                     const Vector<string>&                 a_computedVars,
                     const Vector<DisjointBoxLayout>&      a_computedGrids,
                     const Real                            a_computedDx,
                     const Vector<int>&                    a_computedRefRatio,
                     const Vector<LevelData<FArrayBox>* >& a_exactSoln,
                     const Vector<string>&                 a_exactVars,
                     const Real                            a_exactDx,
                     Real                                  a_bogus_value,
                     bool                                  a_HOaverage,
                     bool                                  a_computeRelativeError)
{
  int numLevels = a_computedSoln.size();

  CH_assert(a_exactSoln.size() == 1);
  CH_assert(a_error.size() == numLevels);
  CH_assert(a_exactDx <= a_computedDx);
  CH_assert(a_computedRefRatio.size() >= numLevels - 1);

  if (a_exactDx == a_computedDx)
  {
    cerr << "Exact dx and computed dx are equal." << endl;
  }

  // check whether input file selects "sum all variables"
  bool sumAll = false;
  ParmParse pp;
  pp.query("sumAll",sumAll);
  
  // const DisjointBoxLayout& exactGrids = a_exactSoln[0]->getBoxes();

  Real dxLevel = a_computedDx;

  // do a bit of sleight-of-hand in the case where there are no
  // ghost cells in the exact solution -- allocate a temporary which
  // _has_ ghost cells, and do a copyTo
  LevelData<FArrayBox>* exactSolnPtr = NULL;
  bool allocatedMemory = false;
  if (a_exactSoln[0]->ghostVect() == IntVect::Zero)
    {
      exactSolnPtr = new LevelData<FArrayBox>(a_exactSoln[0]->getBoxes(),
                                              a_exactSoln[0]->nComp(),
                                              IntVect::Unit);
      a_exactSoln[0]->copyTo(*exactSolnPtr);

      allocatedMemory = true;
     }
   else
     {
       // if there are ghost cells, we can use the exactSoln as-is
       exactSolnPtr = a_exactSoln[0];
     }
   LevelData<FArrayBox>& exactSolnRef = *exactSolnPtr;

   // first need to set boundary conditions on exactsoln
   // this is for the Laplacian which is needed in AverageHO
   DataIterator ditFine = exactSolnRef.dataIterator();
   DomainGhostBC exactBC;
   Interval exactComps(0, a_exactVars.size() - 1);
   for (int dir = 0; dir < SpaceDim; dir++)
   {
     SideIterator sit;
     for (sit.reset(); sit.ok(); ++sit)
     {
       // use HO extrapolation at physical boundaries
       HOExtrapBC thisBC(dir, sit(), exactComps);
       exactBC.setBoxGhostBC(thisBC);
     }
   }

   for (ditFine.begin(); ditFine.ok(); ++ditFine)
   {
     FArrayBox& thisFineSoln = exactSolnRef[ditFine()];
     const Box& fineBox = exactSolnRef.getBoxes()[ditFine()];
     exactBC.applyInhomogeneousBCs(thisFineSoln, fineBox, a_exactDx);
   }
   exactSolnRef.exchange(exactComps);

   // outer loop is over levels
   for (int level = 0; level < numLevels; level++)
   {
     LevelData<FArrayBox>& thisLevelError = *a_error[level];
     LevelData<FArrayBox>& thisLevelComputed = *a_computedSoln[level];

     // compute refinement ratio between solution at this level
     // and exact solution
     Real nRefTemp = (dxLevel / a_exactDx);
     int nRefExact = (int) nRefTemp;

     // this is to do rounding properly if necessary
     if (nRefTemp - nRefExact > 0.5) nRefExact += 1;

     // make sure it's not zero
     if (nRefExact == 0) nRefExact =1;

//.........这里部分代码省略.........

void
CartCellDoubleBoundsPreservingConservativeLinearRefine::refine(Patch<NDIM>& fine,
                                                               const Patch<NDIM>& coarse,
                                                               const int dst_component,
                                                               const int src_component,
                                                               const Box<NDIM>& fine_box,
                                                               const IntVector<NDIM>& ratio)
    const
{
    // Determine the box over which we can apply the bounds-preserving
    // correction, and construct a list of boxes that will not be corrected.
    bool empty_correction_box = false;
    Box<NDIM> correction_box = Box<NDIM>::refine(Box<NDIM>::coarsen(fine_box, ratio), ratio);
    for (unsigned int axis = 0; axis < NDIM; ++axis)
    {
        int& lower = correction_box.lower()(axis);
        while (lower < fine_box.lower()(axis))
        {
            lower += ratio(axis);
        }

        int& upper = correction_box.upper()(axis);
        while (upper > fine_box.upper()(axis))
        {
            upper -= ratio(axis);
        }

        if (lower >= upper)
        {
            empty_correction_box = true;
        }
    }
    const Box<NDIM> coarse_correction_box = Box<NDIM>::coarsen(correction_box, ratio);

    BoxList<NDIM> uncorrected_boxes(fine_box);
    if (!empty_correction_box)
    {
        uncorrected_boxes.removeIntersections(correction_box);
    }

    // Employ limited conservative interpolation to prolong data on the
    // correction box.
    d_conservative_linear_refine_op.refine(
        fine, coarse, dst_component, src_component, correction_box, ratio);

    // Employ constant interpolation to prolong data on the rest of the fine
    // box.
    for (BoxList<NDIM>::Iterator b(uncorrected_boxes); b; b++)
    {
        d_constant_refine_op.refine(fine, coarse, dst_component, src_component, b(), ratio);
    }

    // There is nothing left to do if the correction box is empty.
    if (empty_correction_box) return;

    // Correct the data within the correction box.
    Pointer<CellData<NDIM, double> > fdata = fine.getPatchData(dst_component);
    Pointer<CellData<NDIM, double> > cdata = coarse.getPatchData(src_component);
#if !defined(NDEBUG)
    TBOX_ASSERT(fdata);
    TBOX_ASSERT(cdata);
    TBOX_ASSERT(fdata->getDepth() == cdata->getDepth());
#endif
    const int data_depth = fdata->getDepth();
    const Box<NDIM>& patch_box_crse = coarse.getBox();
    const Index<NDIM>& patch_lower_crse = patch_box_crse.lower();
    const Index<NDIM>& patch_upper_crse = patch_box_crse.upper();
    Pointer<CartesianPatchGeometry<NDIM> > pgeom_crse = coarse.getPatchGeometry();
    for (int depth = 0; depth < data_depth; ++depth)
    {
        for (Box<NDIM>::Iterator b(coarse_correction_box); b; b++)
        {
            const Index<NDIM>& i_crse = b();
            const Index<NDIM> i_fine = i_crse * ratio;

            // Determine the lower/upper bounds.
            Box<NDIM> stencil_box_crse(i_crse, i_crse);
            for (unsigned int axis = 0; axis < NDIM; ++axis)
            {
                if (i_crse(axis) > patch_lower_crse(axis) ||
                    !pgeom_crse->getTouchesRegularBoundary(axis, 0))
                {
                    stencil_box_crse.growLower(axis, 1);
                }
                if (i_crse(axis) < patch_upper_crse(axis) ||
                    !pgeom_crse->getTouchesRegularBoundary(axis, 1))
                {
                    stencil_box_crse.growUpper(axis, 1);
                }
            }

            double l = std::numeric_limits<double>::max();
            double u = -(l - std::numeric_limits<double>::epsilon());
            for (Box<NDIM>::Iterator b(stencil_box_crse); b; b++)
            {
                const double& m = (*cdata)(b(), depth);
                l = std::min(l, m);
                u = std::max(u, m);
            }

//.........这里部分代码省略.........

DTboolean     PrimitiveCollisions::ray_intersect_box (  const Vector3 &from, const Vector3 &direction,
                                                        const Box &b, DTfloat &t)
{
    DTfloat tmin = -std::numeric_limits<DTfloat>::infinity();
    DTfloat tmax = std::numeric_limits<DTfloat>::infinity();
    
    DTfloat t1,t2;
    
    // X axis
    if (direction.x > EPSILON || direction.x < -EPSILON) {
        DTfloat direction_inverse = 1.0F / direction.x;
        t1 = (b.minus_x() - from.x) * direction_inverse;
        t2 = (b.plus_x() - from.x) * direction_inverse;
        
        tmin = MoreMath::min(tmin,t1,t2);
        tmax = MoreMath::max(tmin,t1,t2);
        
        if (tmin > tmax)    return false;
        if (tmax < 0.0F)    return false;
    } else if (from.x < b.minus_x() || from.x > b.plus_x()) {
        return false;
    }   
    
    // Y axis
    if (direction.y > EPSILON || direction.y < -EPSILON) {
        DTfloat direction_inverse = 1.0F / direction.y;
        t1 = (b.minus_y() - from.y) * direction_inverse;
        t2 = (b.plus_y() - from.y) * direction_inverse;
        
        tmin = MoreMath::min(tmin,t1,t2);
        tmax = MoreMath::max(tmin,t1,t2);
        
        if (tmin > tmax)    return false;
        if (tmax < 0.0F)    return false;
    } else if (from.y < b.minus_y() || from.y > b.plus_y()) {
        return false;
    }
    
    // Z axis
    if (direction.z > EPSILON || direction.z < -EPSILON) {
        DTfloat direction_inverse = 1.0F / direction.z;
        t1 = (b.minus_z() - from.z) * direction_inverse;
        t2 = (b.plus_z() - from.z) * direction_inverse;
        
        tmin = MoreMath::min(tmin,t1,t2);
        tmax = MoreMath::max(tmin,t1,t2);
        
        if (tmin > tmax)    return false;
        if (tmax < 0.0F)    return false;
    } else if (from.z < b.minus_z() || from.z > b.plus_z()) {
        return false;
    }
    
    if (tmin <= tmax) {
        t = tmin;
        return true;
    }

    return false;
}

// Set boundary fluxes
void VelIBC::primBC(FArrayBox&            a_WGdnv,
                    const FArrayBox&      a_Wextrap,
                    const FArrayBox&      a_W,
                    const int&            a_dir,
                    const Side::LoHiSide& a_side,
                    const Real&           a_time)
{
  CH_assert(m_isDefined == true);

  // In periodic case, this doesn't do anything
  if (!m_domain.isPeriodic(a_dir))
    {
      // This needs to be fixed
      // CH_assert(m_isBCvalSet);

      int lohisign;
      Box tmp = a_WGdnv.box() & m_domain;
      Real bcVal;

      // Determine which side and thus shifting directions
      if (a_side == Side::Lo)
        {
          lohisign = -1;
          bcVal = m_bcVal[a_dir][0];
        }
      else
        {
          lohisign = 1;
          bcVal = m_bcVal[a_dir][1];
        }

      tmp.shiftHalf(a_dir,lohisign);

      // Is there a domain boundary next to this grid
      if (!m_domain.contains(tmp))
        {
          tmp &= m_domain;

          Box boundaryBox;

          // Find the strip of cells next to the domain boundary
          if (a_side == Side::Lo)
            {
              boundaryBox = bdryLo(tmp,a_dir);
            }
          else
            {
              boundaryBox = bdryHi(tmp,a_dir);
            }

          // Set the boundary fluxes
          FORT_SOLIDVELBCF(CHF_FRA(a_WGdnv),
                           CHF_CONST_FRA(a_Wextrap),
                           CHF_CONST_REAL(bcVal),
                           CHF_CONST_INT(lohisign),
                           CHF_CONST_REAL(m_dx),
                           CHF_CONST_INT(a_dir),
                           CHF_BOX(boundaryBox));
        }
    }
}

// calculate the centroid of a box
KOKKOS_INLINE_FUNCTION
void centroid( Box const &box, Point &c )
{
    for ( int d = 0; d < 3; ++d )
        c[d] = 0.5 * ( box.minCorner()[d] + box.maxCorner()[d] );
}