releasedExamples/AMRINS/util/Gradient.cpp
// utility function to do the actual computin' (to reduce code duplication)
utility function to do the actual computin' (to reduce code duplication)
utility function to do the actual computin ' ( to reduce code duplication )
void
Gradient::singleBoxMacGrad(FArrayBox& a_gradFab,
                           const FArrayBox& a_phiFab,
                           int a_gradComp,
                           int a_phiComp,
                           int a_numComp,
                           const Box& a_edgeBox,
                           Real a_dx,
                           int a_dir,
                           int a_edgeDir,
                           const IntVectSet& a_gridIVS)
{
  Box BLcorner(a_edgeBox);
  Box BRcorner(a_edgeBox);
  Box ULcorner(a_edgeBox);
  Box URcorner(a_edgeBox);
  if (a_edgeDir != a_dir)
    {
      for (int tempDir=0; tempDir<SpaceDim; ++tempDir)
        {
          if ((tempDir==a_edgeDir))
            {
              BLcorner.setBig(tempDir,a_edgeBox.smallEnd(tempDir));
              ULcorner.setBig(tempDir,a_edgeBox.smallEnd(tempDir));
              BRcorner.setSmall(tempDir,a_edgeBox.bigEnd(tempDir));
              URcorner.setSmall(tempDir,a_edgeBox.bigEnd(tempDir));
            }
          else if (tempDir == a_dir)
            {
              BLcorner.setBig(tempDir,a_edgeBox.smallEnd(tempDir));
              BRcorner.setBig(tempDir,a_edgeBox.smallEnd(tempDir));
              ULcorner.setSmall(tempDir,a_edgeBox.bigEnd(tempDir));
              URcorner.setSmall(tempDir,a_edgeBox.bigEnd(tempDir));
            }
        }
    }
  int phiComp = a_phiComp;
  for (int comp=a_gradComp; comp< a_gradComp+a_numComp; comp++)
    {
      FORT_NEWMACGRAD(CHF_FRA1(a_gradFab,comp),
                      CHF_CONST_FRA1(a_phiFab,phiComp),
                      CHF_BOX(a_edgeBox),
                      CHF_CONST_REAL(a_dx),
                      CHF_INT(a_dir),
                      CHF_INT(a_edgeDir));
      if (a_dir != a_edgeDir)
        {
          Real factor = 0.25/a_dx;
          IntVect dirBasis = BASISV(a_dir);
          IntVect edgeBasis = BASISV(a_edgeDir);
          {
            BoxIterator bit(BLcorner);
            for (bit.begin(); bit.ok(); ++bit)
              {
                IntVect thisIV = bit();
                IntVect cornerVect = thisIV -dirBasis - edgeBasis;
                if (!a_gridIVS.contains(cornerVect))
                  {
                    a_gradFab(thisIV, comp) = factor*(3.0*(a_phiFab(thisIV+dirBasis, phiComp)
                                                          -a_phiFab(thisIV-dirBasis, phiComp))
                                                      -( a_phiFab(thisIV+dirBasis+edgeBasis, phiComp)
                                                         -a_phiFab(thisIV-dirBasis+edgeBasis, phiComp)));
                  }
              } 
          }  
          {
            BoxIterator bit(BRcorner);
            for (bit.begin(); bit.ok(); ++bit)
              {
                IntVect thisIV = bit();
                IntVect cornerVect = thisIV -dirBasis;
                if (!a_gridIVS.contains(cornerVect))
                  {
                    a_gradFab(thisIV, comp) = factor*(3.0*(a_phiFab(thisIV+dirBasis-edgeBasis, phiComp)
                                                          -a_phiFab(thisIV-dirBasis-edgeBasis, phiComp))
                                                      -( a_phiFab(thisIV+dirBasis-2*edgeBasis, phiComp)
                                                         -a_phiFab(thisIV-dirBasis-2*edgeBasis, phiComp)));
                  }
              } 
          }  
          {
            BoxIterator bit(ULcorner);
            for (bit.begin(); bit.ok(); ++bit)
              {
                IntVect thisIV = bit();
                IntVect cornerVect = thisIV +dirBasis - edgeBasis;
                if (!a_gridIVS.contains(cornerVect))
                  {
                    a_gradFab(thisIV, comp) = factor*(3.0*(a_phiFab(thisIV+dirBasis, phiComp)
                                                          -a_phiFab(thisIV-dirBasis, phiComp))
                                                      -( a_phiFab(thisIV+dirBasis+edgeBasis, phiComp)
                                                       -a_phiFab(thisIV-dirBasis+edgeBasis, phiComp)));
                  }
              } 
          }  
          {
            BoxIterator bit(URcorner);
            for (bit.begin(); bit.ok(); ++bit)
              {
                IntVect thisIV = bit();
                IntVect cornerVect = thisIV +dirBasis;
                if (!a_gridIVS.contains(cornerVect))
                  {
                    a_gradFab(thisIV, comp) = factor*(3.0*(a_phiFab(thisIV+dirBasis-edgeBasis, phiComp)
                                                          -a_phiFab(thisIV-dirBasis-edgeBasis, phiComp))
                                                      -( a_phiFab(thisIV+dirBasis-2*edgeBasis, phiComp)
                                                         -a_phiFab(thisIV-dirBasis-2*edgeBasis, phiComp)));
                  }
              } 
          }  
        }
      phiComp++;
    }
}
void Gradient :: singleBoxMacGrad ( FArrayBox & a_gradFab , const FArrayBox & a_phiFab , int a_gradComp , int a_phiComp , int a_numComp , const Box & a_edgeBox , Real a_dx , int a_dir , int a_edgeDir , const IntVectSet & a_gridIVS ) { Box BLcorner ( a_edgeBox ) ; Box BRcorner ( a_edgeBox ) ; Box ULcorner ( a_edgeBox ) ; Box URcorner ( a_edgeBox ) ; if ( a_edgeDir != a_dir ) { for ( int tempDir = 0 ; tempDir < SpaceDim ; ++ tempDir ) { if ( ( tempDir == a_edgeDir ) ) { BLcorner . setBig ( tempDir , a_edgeBox . smallEnd ( tempDir ) ) ; ULcorner . setBig ( tempDir , a_edgeBox . smallEnd ( tempDir ) ) ; BRcorner . setSmall ( tempDir , a_edgeBox . bigEnd ( tempDir ) ) ; URcorner . setSmall ( tempDir , a_edgeBox . bigEnd ( tempDir ) ) ; } else if ( tempDir == a_dir ) { BLcorner . setBig ( tempDir , a_edgeBox . smallEnd ( tempDir ) ) ; BRcorner . setBig ( tempDir , a_edgeBox . smallEnd ( tempDir ) ) ; ULcorner . setSmall ( tempDir , a_edgeBox . bigEnd ( tempDir ) ) ; URcorner . setSmall ( tempDir , a_edgeBox . bigEnd ( tempDir ) ) ; } } } int phiComp = a_phiComp ; for ( int comp = a_gradComp ; comp < a_gradComp + a_numComp ; comp ++ ) { FORT_NEWMACGRAD ( CHF_FRA1 ( a_gradFab , comp ) , CHF_CONST_FRA1 ( a_phiFab , phiComp ) , CHF_BOX ( a_edgeBox ) , CHF_CONST_REAL ( a_dx ) , CHF_INT ( a_dir ) , CHF_INT ( a_edgeDir ) ) ; if ( a_dir != a_edgeDir ) { Real factor = 0.25 / a_dx ; IntVect dirBasis = BASISV ( a_dir ) ; IntVect edgeBasis = BASISV ( a_edgeDir ) ; { BoxIterator bit ( BLcorner ) ; for ( bit . begin ( ) ; bit . ok ( ) ; ++ bit ) { IntVect thisIV = bit ( ) ; IntVect cornerVect = thisIV - dirBasis - edgeBasis ; if ( ! a_gridIVS . contains ( cornerVect ) ) { a_gradFab ( thisIV , comp ) = factor * ( 3.0 * ( a_phiFab ( thisIV + dirBasis , phiComp ) - a_phiFab ( thisIV - dirBasis , phiComp ) ) - ( a_phiFab ( thisIV + dirBasis + edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis + edgeBasis , phiComp ) ) ) ; } } } { BoxIterator bit ( BRcorner ) ; for ( bit . begin ( ) ; bit . ok ( ) ; ++ bit ) { IntVect thisIV = bit ( ) ; IntVect cornerVect = thisIV - dirBasis ; if ( ! a_gridIVS . contains ( cornerVect ) ) { a_gradFab ( thisIV , comp ) = factor * ( 3.0 * ( a_phiFab ( thisIV + dirBasis - edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis - edgeBasis , phiComp ) ) - ( a_phiFab ( thisIV + dirBasis - 2 * edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis - 2 * edgeBasis , phiComp ) ) ) ; } } } { BoxIterator bit ( ULcorner ) ; for ( bit . begin ( ) ; bit . ok ( ) ; ++ bit ) { IntVect thisIV = bit ( ) ; IntVect cornerVect = thisIV + dirBasis - edgeBasis ; if ( ! a_gridIVS . contains ( cornerVect ) ) { a_gradFab ( thisIV , comp ) = factor * ( 3.0 * ( a_phiFab ( thisIV + dirBasis , phiComp ) - a_phiFab ( thisIV - dirBasis , phiComp ) ) - ( a_phiFab ( thisIV + dirBasis + edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis + edgeBasis , phiComp ) ) ) ; } } } { BoxIterator bit ( URcorner ) ; for ( bit . begin ( ) ; bit . ok ( ) ; ++ bit ) { IntVect thisIV = bit ( ) ; IntVect cornerVect = thisIV + dirBasis ; if ( ! a_gridIVS . contains ( cornerVect ) ) { a_gradFab ( thisIV , comp ) = factor * ( 3.0 * ( a_phiFab ( thisIV + dirBasis - edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis - edgeBasis , phiComp ) ) - ( a_phiFab ( thisIV + dirBasis - 2 * edgeBasis , phiComp ) - a_phiFab ( thisIV - dirBasis - 2 * edgeBasis , phiComp ) ) ) ; } } } } phiComp ++ ; } }
utility function to do the actual computin' (to reduce code duplication)
utility function to do the actual computin ' ( to reduce code duplication )
// corner boxes identify which cells need special // stencil for transverse derivatives // "bottom", "top", "left", and "right" are relative // to edgeDir=0, dir=1 // now loop over directions and set box edge extents // loop over components // now compute gradient in direction dir on this grid // now do corner cells specially (replace gradient values computed // in fortran with one-sided gradients where necessary) // "lower left" // end loop over cells in BLcorner // end BLcorner context // "lower right" // end loop over cells in BRcorner // end BRcorner context // "upper left" // end loop over cells in ULcorner // end ULcorner context // "upper right" // end loop over cells in URcorner // end URcorner context