lib/src/AMRTools/TensorCFInterp.cpp
/***********************/
//get extended phi (lives next to interpivs)
/***********************/
get extended phi (lives next to interpivs)
get extended phi ( lives next to interpivs )
void
TensorCFInterp::getPhiStar(BaseFab<Real> & a_phistar,
                           BaseFab<Real> & a_tanGradStar,
                           const BaseFab<Real> & a_phic,
                           const QuadCFStencil& a_qcfs,
                           const Side::LoHiSide a_hiorlo,
                           const int a_idir,
                           const Interval& a_variables) const
{
  CH_assert(isDefined());
  CH_assert(a_qcfs.isDefined());
  Real dxf = m_dxFine;
  Real dxc = m_refRatio*dxf;
  if (!a_qcfs.isEmpty())
    {
      CH_assert(m_level > 0);
      const IntVectSet& interp_ivs = a_qcfs.getFineIVS();
      const IntVectSet& coarsl_ivs = a_qcfs.getCoarIVS();
      Box coarinterpbox = coarsl_ivs.minBox();
      int ncomp = a_phic.nComp();
      CH_assert(ncomp == m_nComp);
      CH_assert(a_phic.box().contains((coarinterpbox)));
      int ihilo = sign(a_hiorlo);
      Box phistarbox = interp_ivs.minBox();
      phistarbox.shift(a_idir, ihilo);
      a_phistar.resize(phistarbox, ncomp);
      a_phistar.setVal(666.666);
      a_tanGradStar.resize(phistarbox, ncomp*SpaceDim);
      a_tanGradStar.setVal(666.666);
      for (int ivar = a_variables.begin(); ivar <= a_variables.end(); ivar++)
        {
          BaseFab<Real>  coarslope(coarinterpbox, SpaceDim);
          BaseFab<Real>  coarcurva(coarinterpbox, SpaceDim);
#if (CH_SPACEDIM == 3)
          Vector<int> xydir;
          for (int isldir = 0; isldir < SpaceDim; isldir++)
            if (isldir != a_idir)
              xydir.push_back(isldir);
          BaseFab<Real>  coarmixed(coarinterpbox, 1);
#endif
          coarslope.setVal(0.);
          IVSIterator coar_ivsit(coarsl_ivs);
          for (coar_ivsit.begin(); coar_ivsit.ok(); ++coar_ivsit)
            {
              IntVect coariv = coar_ivsit ();
              for (int isldir = 0; isldir < SpaceDim; isldir++)
                {
                  if (isldir == a_idir)
                    {
                      coarslope(coariv, isldir) = 0.0;
                      coarcurva(coariv, isldir) = 0.0;
                    }
                  else
                    {
                      coarslope(coariv, isldir) =
                        a_qcfs.computeFirstDerivative(a_phic, isldir,
                                                      ivar, coariv, dxc);
                      coarcurva(coariv, isldir) =
                        a_qcfs.computeSecondDerivative(a_phic, isldir,
                                                       ivar, coariv, dxc);
                    } 
                } 
#if (CH_SPACEDIM == 3)
              coarmixed(coariv) =
                a_qcfs.computeMixedDerivative(a_phic, ivar, coariv, dxc);
#endif
            } 
          IVSIterator fine_ivsit(interp_ivs);
          for (fine_ivsit.begin(); fine_ivsit.ok(); ++fine_ivsit)
            {
              IntVect ivf = fine_ivsit();
              IntVect ivc = coarsen(ivf,m_refRatio);
              IntVect ivstar = ivf;
              ivstar.shift(a_idir, ihilo);
              a_phistar(ivstar, ivar)  = a_phic(ivc,ivar);
              for (int gradDir=0; gradDir<SpaceDim; gradDir++)
                {
                  if (gradDir != a_idir)
                    {
                      int gradComp = gradIndex(ivar, gradDir);
                      a_tanGradStar(ivstar, gradComp) = coarslope(ivc,gradDir);
                    }
                }
              for (int isldir = 0; isldir < SpaceDim; isldir++)
                {
                  if (isldir != a_idir)
                    {
                      int jf = ivf[isldir];
                      int jc = ivc[isldir];
                      Real xf = (jf+0.5)*dxf;
                      Real xc = (jc+0.5)*dxc;
                      Real x = xf-xc;
                      a_phistar(ivstar, ivar) +=
                        (      x*coarslope(ivc, isldir)) +
                        (0.5*x*x*coarcurva(ivc, isldir));
                      int gradComp = gradIndex(ivar, isldir);
                      a_tanGradStar(ivstar, gradComp) +=
                        x*coarcurva(ivc,isldir);
                    }
                } 
#if (CH_SPACEDIM==3)
              int dif = ivf[xydir[0]];
              int dic = ivc[xydir[0]];
              int djf = ivf[xydir[1]];
              int djc = ivc[xydir[1]];
              Real xf = (dif+0.5)*dxf;
              Real xc = (dic+0.5)*dxc;
              Real yf = (djf+0.5)*dxf;
              Real yc = (djc+0.5)*dxc;
              Real x = xf-xc;
              Real y = yf-yc;
              a_phistar(ivstar, ivar) += x*y*coarmixed(ivc);
              CH_assert (SpaceDim <= 3);
              int gradComp = gradIndex(ivar,xydir[0]);
              a_tanGradStar(ivstar, gradComp) += y*coarmixed(ivc);
              gradComp = gradIndex(ivar,xydir[1]);
              a_tanGradStar(ivstar, gradComp) += x*coarmixed(ivc);
#endif
            } 
        } 
    } 
}
void TensorCFInterp :: getPhiStar ( BaseFab < Real > & a_phistar , BaseFab < Real > & a_tanGradStar , const BaseFab < Real > & a_phic , const QuadCFStencil & a_qcfs , const Side :: LoHiSide a_hiorlo , const int a_idir , const Interval & a_variables ) const { CH_assert ( isDefined ( ) ) ; CH_assert ( a_qcfs . isDefined ( ) ) ; Real dxf = m_dxFine ; Real dxc = m_refRatio * dxf ; if ( ! a_qcfs . isEmpty ( ) ) { CH_assert ( m_level > 0 ) ; const IntVectSet & interp_ivs = a_qcfs . getFineIVS ( ) ; const IntVectSet & coarsl_ivs = a_qcfs . getCoarIVS ( ) ; Box coarinterpbox = coarsl_ivs . minBox ( ) ; int ncomp = a_phic . nComp ( ) ; CH_assert ( ncomp == m_nComp ) ; CH_assert ( a_phic . box ( ) . contains ( ( coarinterpbox ) ) ) ; int ihilo = sign ( a_hiorlo ) ; Box phistarbox = interp_ivs . minBox ( ) ; phistarbox . shift ( a_idir , ihilo ) ; a_phistar . resize ( phistarbox , ncomp ) ; a_phistar . setVal ( 666.666 ) ; a_tanGradStar . resize ( phistarbox , ncomp * SpaceDim ) ; a_tanGradStar . setVal ( 666.666 ) ; for ( int ivar = a_variables . begin ( ) ; ivar <= a_variables . end ( ) ; ivar ++ ) { BaseFab < Real > coarslope ( coarinterpbox , SpaceDim ) ; BaseFab < Real > coarcurva ( coarinterpbox , SpaceDim ) ; #if ( CH_SPACEDIM == 3 ) Vector < int > xydir ; for ( int isldir = 0 ; isldir < SpaceDim ; isldir ++ ) if ( isldir != a_idir ) xydir . push_back ( isldir ) ; BaseFab < Real > coarmixed ( coarinterpbox , 1 ) ; #endif coarslope . setVal ( 0. ) ; IVSIterator coar_ivsit ( coarsl_ivs ) ; for ( coar_ivsit . begin ( ) ; coar_ivsit . ok ( ) ; ++ coar_ivsit ) { IntVect coariv = coar_ivsit ( ) ; for ( int isldir = 0 ; isldir < SpaceDim ; isldir ++ ) { if ( isldir == a_idir ) { coarslope ( coariv , isldir ) = 0.0 ; coarcurva ( coariv , isldir ) = 0.0 ; } else { coarslope ( coariv , isldir ) = a_qcfs . computeFirstDerivative ( a_phic , isldir , ivar , coariv , dxc ) ; coarcurva ( coariv , isldir ) = a_qcfs . computeSecondDerivative ( a_phic , isldir , ivar , coariv , dxc ) ; } } #if ( CH_SPACEDIM == 3 ) coarmixed ( coariv ) = a_qcfs . computeMixedDerivative ( a_phic , ivar , coariv , dxc ) ; #endif } IVSIterator fine_ivsit ( interp_ivs ) ; for ( fine_ivsit . begin ( ) ; fine_ivsit . ok ( ) ; ++ fine_ivsit ) { IntVect ivf = fine_ivsit ( ) ; IntVect ivc = coarsen ( ivf , m_refRatio ) ; IntVect ivstar = ivf ; ivstar . shift ( a_idir , ihilo ) ; a_phistar ( ivstar , ivar ) = a_phic ( ivc , ivar ) ; for ( int gradDir = 0 ; gradDir < SpaceDim ; gradDir ++ ) { if ( gradDir != a_idir ) { int gradComp = gradIndex ( ivar , gradDir ) ; a_tanGradStar ( ivstar , gradComp ) = coarslope ( ivc , gradDir ) ; } } for ( int isldir = 0 ; isldir < SpaceDim ; isldir ++ ) { if ( isldir != a_idir ) { int jf = ivf [ isldir ] ; int jc = ivc [ isldir ] ; Real xf = ( jf + 0.5 ) * dxf ; Real xc = ( jc + 0.5 ) * dxc ; Real x = xf - xc ; a_phistar ( ivstar , ivar ) += ( x * coarslope ( ivc , isldir ) ) + ( 0.5 * x * x * coarcurva ( ivc , isldir ) ) ; int gradComp = gradIndex ( ivar , isldir ) ; a_tanGradStar ( ivstar , gradComp ) += x * coarcurva ( ivc , isldir ) ; } } #if ( CH_SPACEDIM == 3 ) int dif = ivf [ xydir [ 0 ] ] ; int dic = ivc [ xydir [ 0 ] ] ; int djf = ivf [ xydir [ 1 ] ] ; int djc = ivc [ xydir [ 1 ] ] ; Real xf = ( dif + 0.5 ) * dxf ; Real xc = ( dic + 0.5 ) * dxc ; Real yf = ( djf + 0.5 ) * dxf ; Real yc = ( djc + 0.5 ) * dxc ; Real x = xf - xc ; Real y = yf - yc ; a_phistar ( ivstar , ivar ) += x * y * coarmixed ( ivc ) ; CH_assert ( SpaceDim <= 3 ) ; int gradComp = gradIndex ( ivar , xydir [ 0 ] ) ; a_tanGradStar ( ivstar , gradComp ) += y * coarmixed ( ivc ) ; gradComp = gradIndex ( ivar , xydir [ 1 ] ) ; a_tanGradStar ( ivstar , gradComp ) += x * coarmixed ( ivc ) ; #endif } } } }
get extended phi (lives next to interpivs)
get extended phi ( lives next to interpivs )
//if cfsten is empty, nothing to interpolate. //if there IS something to interpolate, the level ident //had better be greater than zero.  Otherwise a null //was sent in as coarse grids on construction // allocate phistar here // tanGradStar has same location as phiStar // notice that we're allocating SpaceDim components for // each component in phi (for each coordinate direction). // possible memory optimization would only allocate SpaceDim-1 // components (won't need normal dir) //phi = phino + slope*x + half*x*x*curvature //directions != idir,  needed for mixed derivative in 3d //add to phi x*y*mixed_derivative //first find extended value phistar //get slopes of coarse solution //end if (isldir == idir) // end loop over slope dir //end loop over coarse intvects // initialize tanGrad to derivatives at coarse point // grad will be stored so that dimensions // vary fastest // increment tangential gradient //end if (isldir != idir) //end loop over slope directions //add in mixed derivative component // mixed derivative gets added to each component of // tanGrad.  this is simple enough to do in 3d -- // we can just hardwire it.  if higher dimensional // problem, this needs to be rethought. // first do direction 0 // now do direction 1 //end loop over fine intvects //end loop over variables //end if (level>0 && !hocfs.isempty())