Chombo + EB: EBArith.H Source File

00001 #ifdef CH_LANG_CC
00002 /*
00003  *      _______              __
00004  *     / ___/ /  ___  __ _  / /  ___
00005  *    / /__/ _ \/ _ \/  V \/ _ \/ _ \
00006  *    \___/_//_/\___/_/_/_/_.__/\___/
00007  *    Please refer to Copyright.txt, in Chombo's root directory.
00008  */
00009 #endif
00010 
00011 #ifndef _EBARITH_H_
00012 #define _EBARITH_H_
00013 
00014 #include <fstream>
00015 
00016 #include "EBCellFAB.H"
00017 #include "EBISBox.H"
00018 #include "EBISLayout.H"
00019 #include "EBFaceFAB.H"
00020 #include "BoxLayout.H"
00021 #include "BoxLayoutData.H"
00022 #include "Stencils.H"
00023 #include "LayoutData.H"
00024 #include "Tuple.H"
00025 #include "EBFluxFAB.H"
00026 #include "EBCellFAB.H"
00027 #include "NamespaceHeader.H"
00028 ///
00029 /**
00030    enumeration class to ennumerate the different
00031    extents of norm covereage. sometimes one
00032    wants to only take a norm over the irregular
00033    cells.  sometimes one only wants to take a norm
00034    over the regular cells, and sometimes one wants
00035    both.
00036  */
00037 class EBNormType
00038 {
00039 public:
00040 
00041   ///
00042   /**
00043      OverOnlyRegular means take the norm only over
00044      regular cells. \\
00045      OverOnlyirregular means take the norm only over
00046      irregular cells. \\
00047      OverBoth means use all uncovered cells. \\
00048      All norms exclude covered cells.
00049    */
00050   enum  NormMode
00051   {
00052     OverOnlyRegular=0,
00053     OverOnlyIrregular,
00054     OverBoth
00055   };
00056 };
00057 
00058 ///
00059 /**
00060    class to encapsulate the common arithmatic operations
00061    for ebfabs
00062  */
00063 class EBArith
00064 {
00065 public:
00066 
00067   ///
00068   static void convertToITM(IndexTM<Real, SpaceDim>& a_diffrv,  const RealVect& a_rv);
00069 
00070 
00071   ///
00072   static void convertToITM(IndexTM<int, SpaceDim>& a_diffrv,  const IntVect& a_rv);
00073 
00074   ///
00075   static void convertToITM(IndexTM<Real, SpaceDim-1>& a_diffrv,  const RealVect& a_rv, const int& a_ignoreIndex);
00076 
00077   ///
00078   static void convertToITM(IndexTM<int, SpaceDim-1>& a_diffrv,  const IntVect& a_rv, const int& a_ignoreIndex);
00079 
00080   ///
00081   /**
00082      Gets the stencil to take the first derivative of  cell centered data.
00083      no high order one sided stuff (to keep stencil close)
00084   */
00085   static int
00086   getFirstDerivStencilWidthOne(VoFStencil&      a_sten,
00087                                const VolIndex&  a_vof,
00088                                const EBISBox&   a_ebisBox,
00089                                const int&       a_idir,
00090                                const Real&      a_dx,
00091                                IntVectSet*    a_cfivsPtr,
00092                                int ivar);
00093 
00094   ///
00095   /**
00096      Least square sten using all cells in a monotone path of radius a_rad
00097    */
00098   static void
00099   getLeastSquaresGradStenAllVoFsRad(VoFStencil&          a_stencil,
00100                                     Real&                a_weight,
00101                                     const RealVect&      a_normal,
00102                                     const RealVect&      a_centroid,
00103                                     const VolIndex&      a_vof,
00104                                     const EBISBox&       a_ebisBox,
00105                                     const RealVect&      a_dx,
00106                                     const ProblemDomain& a_domain,
00107                                     int                  a_ivar,
00108                                     int                  a_rad);
00109 
00110   static int
00111   getFirstOrderExtrapolationStencil(VoFStencil&     a_stencil,
00112                                     const RealVect& a_dist,
00113                                     const RealVect& a_dx,
00114                                     const VolIndex& a_startVoF,
00115                                     const EBISBox&  a_ebisBox,
00116                                     int a_noExtrapThisDir,
00117                                     IntVectSet*    a_cfivsPtr,
00118                                     int ivar);
00119   ///
00120   /**
00121    */
00122   static int orderScript(int icomp, int inorm, int ncomp);
00123 
00124   ///
00125   static void
00126   compareError(Vector<Real>&            a_orders,
00127                const  LevelData<EBCellFAB>&    a_errorFine,
00128                const  LevelData<EBCellFAB>&    a_errorCoar,
00129                const  DisjointBoxLayout&       a_gridsFine,
00130                const  DisjointBoxLayout&       a_gridsCoar,
00131                const  EBISLayout&              a_ebislFine,
00132                const  EBISLayout&              a_ebislCoar,
00133                const Box&                      a_coarseDom,
00134                int                             a_testverbosity,
00135                fstream*                        a_fout,
00136                Vector<string>                  a_names,
00137                string  a_prefix)          ;
00138   ///
00139   /**o
00140    */
00141   static void
00142   compareError(Vector<Real>&                            a_orders,
00143                const Vector< LevelData<EBCellFAB>* >&   a_errorFine,
00144                const Vector< LevelData<EBCellFAB>* >&   a_errorCoar,
00145                const Vector< DisjointBoxLayout >&       a_gridsFine,
00146                const Vector< DisjointBoxLayout >&       a_gridsCoar,
00147                const Vector< EBISLayout >&              a_ebislFine,
00148                const Vector< EBISLayout >&              a_ebislCoar,
00149                const Vector<int>&                       a_refRat,
00150                const Box&                               a_coarseDom,
00151                int                                      a_testverbosity,
00152                fstream*                                 a_fout = NULL,
00153                Vector<string> names = Vector<string>(),
00154                string prefix = string());
00155 
00156   static Real
00157   irregNorm(const LevelData<BaseIVFAB<Real> >    & a_dataOne,
00158             const DisjointBoxLayout              & a_layout,
00159             const EBISLayout                     & a_ebisl,
00160             const int                            & a_comp,
00161             const int                            & a_p);
00162   static void 
00163   irregNorm(Real& a_ebIrregNorm,
00164             const BaseIVFAB<Real>& a_ebiError,
00165             const IntVectSet& a_ivsIrreg,
00166             const EBISBox& a_ebisBox,
00167             const int& a_comp,
00168             const int& a_normtype);
00169   static void
00170   compareIrregError(Vector<Real>&                          a_orders,
00171                     const  LevelData<BaseIVFAB<Real> >&    a_errorFine,
00172                     const  LevelData<BaseIVFAB<Real> >&    a_errorCoar,
00173                     const  DisjointBoxLayout&              a_gridsFine,
00174                     const  DisjointBoxLayout&              a_gridsCoar,
00175                     const  EBISLayout&                     a_ebislFine,
00176                     const  EBISLayout&                     a_ebislCoar,
00177                     const Box&                             a_coarseDom,
00178                     const string &                         a_prefix,
00179                     Vector<string>                         a_names) ;         
00180   ///
00181   static void shrinkIVS(IntVectSet& a_ivs, const int& a_numShrink);
00182 
00183   static void
00184   timeInterpolate(LevelData<EBCellFAB>&       a_U,
00185                   const LevelData<EBCellFAB>& a_UOld,
00186                   const LevelData<EBCellFAB>& a_UNew,
00187                   const DisjointBoxLayout&    a_grids,
00188                   const Real&                 a_time,
00189                   const Real&                 a_told,
00190                   const Real&                 a_tnew);
00191   ///
00192   static Real
00193   getDiagWeight(  VoFStencil&     a_vofStencil,
00194                   const VolIndex& a_vof,
00195                   int a_var = 0);
00196 
00197   ///
00198   static void getMultiColors(Vector<IntVect>& a_colors);
00199 
00200   ///
00201   static void
00202   interpolateCFH(EBCellFAB&                  a_phi,
00203                  const int  &                a_idir,
00204                  const Side::LoHiSide&       a_hiorlo,
00205                  const EBISBox&              a_ebisBox,
00206                  const Real&                 a_dxfine,
00207                  const Real&                 a_dxcoar,
00208                  const IntVectSet&           a_interpIVS);
00209 
00210   ///
00211   /**
00212      returns the order of the extrapolation.
00213      the reason for the last argument is that you might not want the stencil
00214      to leak over in the noExtrap direction even though you have set a_dist to zero.
00215      This happens in CF interpolation where you have to really worry about stencil width.
00216    **/
00217   static int
00218   getExtrapolationStencil(VoFStencil&     a_stencil,
00219                           const RealVect& a_dist,
00220                           const RealVect& a_dx,
00221                           const VolIndex& a_startVoF,
00222                           const EBISBox&  a_ebisBox,
00223                           int a_noExtrapThisDirection = -1,
00224                           IntVectSet*    a_cfivsPtr = NULL,
00225                           int ivar = 0);
00226 
00227   ///
00228   /**
00229      returns the order of the extrapolation.
00230      the reason for the last argument is that you might not want the stencil
00231      to leak over in the noExtrap direction even though you have set a_dist to zero.
00232      This happens in CF interpolation where you have to really worry about stencil width.
00233    **/
00234   static int
00235   get1stOrderExtrapolationStencil(VoFStencil&     a_stencil,
00236                                   const RealVect& a_dist,
00237                                   const RealVect& a_dx,
00238                                   const VolIndex& a_startVoF,
00239                                   const EBISBox&  a_ebisBox,
00240                                   int a_noExtrapThisDirection = -1,
00241                                   IntVectSet*    a_cfivsPtr = NULL,
00242                                   int ivar = 0);
00243 
00244   ///
00245   /**
00246    **/
00247   static void
00248   getDir1Dir2(int& a_dir1, int& a_dir2, const int& a_dir);
00249 
00250   ///
00251   /**
00252      send null ivs if you do not want lohicenter to care about cfivs
00253    **/
00254   static void
00255   loHiCenter(Box&                 a_loBox,
00256              int&                 a_hasLo,
00257              Box&                 a_hiBox,
00258              int&                 a_hasHi,
00259              Box&                 a_centerBox,
00260              const ProblemDomain  &   a_eblg,
00261              const Box&           a_inBox,
00262              const int&           a_dir,
00263              IntVectSet*    a_cfivsPtr = NULL);
00264 
00265   ///
00266   /**
00267    **/
00268   static void
00269   loHi(Box&                 a_loBox,
00270        int&                 a_hasLo,
00271        Box&                 a_hiBox,
00272        int&                 a_hasHi,
00273        const ProblemDomain& a_eblg,
00274        const Box&           a_inBox,
00275        const int&           a_dir);
00276 
00277   ///
00278   /**
00279      testRef is the size of the coarsest domain allowed in multigrid. If testRef=2,
00280      then the coarsest domain in multigrid will be 2x2(x2)
00281    **/
00282   static bool getCoarserLayouts(DisjointBoxLayout&       a_dblCoar,
00283                                 ProblemDomain&           a_domainCoar,
00284                                 const DisjointBoxLayout& a_dblFine,
00285                                 const ProblemDomain&     a_domainFine,
00286                                 int                      a_refToCoar,
00287                                 int                      a_maxBoxSize,
00288                                 bool&                    a_layoutChanged,
00289                                 int                      a_testRef = 2);
00290 
00291   ///
00292   /**
00293      Gets the stencil to take the first derivative in the given direction
00294      of  cell centered data.
00295      Returns the expected order of the derivative.
00296      When we need them, we prefer
00297      first derivataves to be second order if at all possible, so we take some pains
00298      to achieve that
00299   */
00300   static int
00301   getFirstDerivStencil(VoFStencil&      a_sten,
00302                        const VolIndex&  a_vof,
00303                        const EBISBox&   a_ebisBox,
00304                        const int&       a_idir,
00305                        const Real&      a_dx,
00306                        IntVectSet*    a_cfivsPtr = NULL,
00307                        int ivar = 0);
00308 
00309 
00310   ///
00311   /**
00312      Gets the stencil to take the mixed (in the given directions)
00313      derivative of  cell centered data.
00314      Returns the expected order of the derivative.
00315      When we need them, we usually only need second derivatives to O(h), so
00316      this just shifts stencil when it has to.
00317   */
00318   static int
00319   getMixedDerivStencil(VoFStencil&      a_sten,
00320                        const VolIndex&  a_vof,
00321                        const EBISBox&   a_ebisBox,
00322                        const int&       a_dir1,
00323                        const int&       a_dir2,
00324                        const Real&      a_dx1,
00325                        const Real&      a_dx2,
00326                        IntVectSet*    a_cfivsPtr = NULL,
00327                        int ivar = 0);
00328 
00329   ///
00330   /**
00331      Gets the stencil to take the second derivative in the given direction)
00332      of  cell centered data.
00333      Returns the expected order of the derivative.
00334      When we need them, we usually only need second derivatives to O(h), so
00335      this just shifts stencil when it has to.
00336   */
00337   static int
00338   getSecondDerivStencil(VoFStencil&      a_sten,
00339                         const VolIndex&  a_vof,
00340                         const EBISBox&   a_ebisBox,
00341                         const int&       a_idir,
00342                         const Real&      a_dx,
00343                         IntVectSet*    a_cfivsPtr = NULL,
00344                         int ivar = 0);
00345 
00346   ///
00347   /**
00348      a function that could possibly take over the world.
00349   */
00350   static void
00351   getVoFsDir(bool& a_hasClose, VolIndex& a_closeVoF,
00352              bool& a_hasFar,   VolIndex& a_farVoF,
00353              const EBISBox& a_ebisBox,
00354              const VolIndex& a_vof,
00355              int a_idir, Side::LoHiSide a_sd,
00356              IntVectSet*    a_cfivsPtr);
00357 
00358   ///
00359   /**
00360      define an intvectset of ghost cells which live on
00361      the coarse fine interface.  includes corner cells.
00362    */
00363   static void
00364   defineCFIVS(LayoutData<IntVectSet>&   a_cfivs,
00365               const DisjointBoxLayout&  a_grids,
00366               const ProblemDomain&      a_probDom);
00367 
00368   ///
00369   static Real
00370   extrapFaceGradToOutflow(const FaceIndex&      a_bndryFace,
00371                           const Side::LoHiSide& a_side,
00372                           const int&            a_idir,
00373                           const EBGraph&        a_ebGraph,
00374                           const EBFaceFAB&      a_faceData,
00375                           const int&            a_comp);
00376 
00377   ///
00378   static Real
00379   extrapFaceValueToDomain(const FaceIndex&      a_bndryFace,
00380                           const Side::LoHiSide& a_side,
00381                           const int&            a_idir,
00382                           const EBGraph&        a_ebGraph,
00383                           const EBFaceFAB&      a_faceData,
00384                           const int&            a_comp,
00385                           const Real&           a_dropOrderValue = 0.);
00386 
00387   ///
00388   static Real
00389   extrapFaceVelToOutflow(const FaceIndex&      a_bndryFace,
00390                          const Side::LoHiSide& a_side,
00391                          const int&            a_idir,
00392                          const EBGraph&        a_ebGraph,
00393                          const EBFaceFAB&      a_faceData,
00394                          const int&            a_comp);
00395 
00396   ///
00397   static Real
00398   interpolateVel(const EBFaceFAB& a_vel,
00399                  const EBISBox&   a_ebisBox,
00400                  const FaceIndex& a_face);
00401 
00402   ///
00403   static Real
00404   deInterpolateVel(const Real&      a_centroidVel,
00405                    const EBFaceFAB& a_vel,
00406                    const EBISBox&   a_ebisBox,
00407                    const FaceIndex& a_face);
00408 
00409   ///
00410   static Real
00411   getFaceVelForDivFreeCell(const FaceIndex&      a_face,
00412                            const VolIndex&       a_vof,
00413                            const int&            a_idir,
00414                            const Side::LoHiSide& a_side,
00415                            const EBFluxFAB&      a_vel,
00416                            const RealVect&       a_dx,
00417                            const EBISBox&        a_ebisBox);
00418   ///
00419   /**
00420      evaluates integral(a_divu dV)/integral(dV) over a hierarchy.
00421      pval = -1--->assumes volfrac already multiplied in
00422      pval = -2--->multiply by volume fraction
00423    */
00424   static void
00425   meanOverHierarchy(Real&                                   a_mean,
00426                     const Vector<LevelData<EBCellFAB> *>&   a_divu,
00427                     const Vector<DisjointBoxLayout>&        a_grids,
00428                     const Vector<EBISLayout>&               a_ebisl,
00429                     const Vector<int>&                      a_refRat,
00430                     const ProblemDomain&                    a_domainCoarsest,
00431                     const RealVect&                         a_dxCoarsest,
00432                     const int& a_comp,
00433                     int pval= -1);
00434 
00435   ///
00436   /**
00437    //gets the weight of a vof to interpolate (linear in 2d, bilinear in 3d) to an interpolation point in the plane
00438    //this is used to get one of the points along the johansen ray
00439    */
00440   static Real
00441   getJohanVoFWeight(const VolIndex& a_whichVoF,     //vof for which we are trying to find weight
00442                     const RealVect& a_intersectLoc, //location of intersection in real space
00443                     const IntVect&  a_intersectIV,  //intvect through which ray passes
00444                     const int&      a_nMaxDir,      //direction of largest normal
00445                     const RealVect& a_dx);           //grid spacing but you knew that.
00446 
00447   ///version which uses incoming normal, boundary centroid
00448   /**
00449      splitting up stuff this way to facillitate multifluid
00450      which can have multiple normals and boundary centroids per cell.
00451    */
00452   static void
00453   johanStencil(bool&               a_dropOrder,
00454                Vector<VoFStencil>& a_pointStencils,
00455                Vector<Real>&       a_distanceAlongLine,
00456                const RealVect&     a_normal,
00457                const RealVect&     a_bndryCentroid,
00458                const VolIndex&     a_vof,
00459                const EBISBox&      a_ebisBox,
00460                const RealVect&     a_dx,
00461                const IntVectSet&   a_cfivs,
00462                int a_ivar = 0);
00463 
00464   ///
00465   /**
00466      Gets the stencils to get the data points along the Johansen ray.   If a_dropOrder
00467      returns true then the stencil does not exist.  Also returns the distances along the ray
00468      to each point.
00469    */
00470   static void
00471   johanStencil(bool&               a_dropOrder,
00472                Vector<VoFStencil>& a_pointStencils,
00473                Vector<Real>&       a_distanceAlongLine,
00474                const VolIndex&     a_vof,
00475                const EBISBox&      a_ebisBox,
00476                const RealVect&     a_dx,
00477                const IntVectSet&   a_cfivs,
00478                int a_ivar = 0);
00479 
00480   ///standard interface
00481   /**
00482      gets the normal and calls the other version
00483    */
00484   static void
00485   getLeastSquaresGradSten(VoFStencil&     a_stencil,
00486                           Real&           a_weight,
00487                           const VolIndex& a_vof,
00488                           const EBISBox&  a_ebisBox,
00489                           const RealVect& a_dx,
00490                           const ProblemDomain& a_domain,
00491                           int a_ivar = 0);
00492 
00493   ///version which uses incoming normal, boundary centroid
00494   /**
00495      splitting up stuff this way to facillitate multifluid
00496      which can have multiple normals and boundary centroids per cell.
00497    */
00498   static void
00499   getLeastSquaresGradSten(VoFStencil&     a_stencil,
00500                           Real&           a_weight,
00501                           const RealVect& a_normal  ,
00502                           const RealVect& a_centroid,
00503                           const VolIndex& a_vof,
00504                           const EBISBox&  a_ebisBox,
00505                           const RealVect& a_dx,
00506                           const ProblemDomain& a_domain,
00507                           int a_ivar = 0);
00508 
00509 
00510   ///version which uses incoming normal, boundary centroid
00511   /**
00512      splitting up stuff this way to facillitate multifluid
00513      which can have multiple normals and boundary centroids per cell.
00514    */
00515   static void
00516   getLeastSquaresGradSten(VoFStencil&     a_stencil,
00517                           Real&           a_weight,
00518                           const RealVect& a_normal  ,
00519                           const RealVect& a_centroid,
00520                           const IntVect&  a_quadrant,
00521                           const VolIndex& a_vof,
00522                           const EBISBox&  a_ebisBox,
00523                           const RealVect& a_dx,
00524                           const ProblemDomain& a_domain,
00525                           int a_ivar = 0);
00526 
00527   ///version which uses incoming normal, boundary centroid, and
00528   /// all available VoFs.
00529   /**
00530      splitting up stuff this way to facillitate multifluid
00531      which can have multiple normals and boundary centroids per cell.
00532    */
00533   static void
00534   getLeastSquaresGradStenAllVoFs(VoFStencil&          a_stencil,
00535                                  Real&                a_weight,
00536                                  const RealVect&      a_normal,
00537                                  const RealVect&      a_centroid,
00538                                  const VolIndex&      a_vof,
00539                                  const EBISBox&       a_ebisBox,
00540                                  const RealVect&      a_dx,
00541                                  const ProblemDomain& a_domain,
00542                                  int                  a_ivar);
00543 
00544   ///version which uses incoming normal, boundary centroid, and
00545   /// all quadrants.
00546   /**
00547      splitting up stuff this way to facillitate multifluid
00548      which can have multiple normals and boundary centroids per cell.
00549    */
00550   static void
00551   getLeastSquaresGradStenAllQuad(VoFStencil&          a_stencil,
00552                                  Real&                a_weight,
00553                                  const RealVect&      a_normal,
00554                                  const RealVect&      a_centroid,
00555                                  const VolIndex&      a_vof,
00556                                  const EBISBox&       a_ebisBox,
00557                                  const RealVect&      a_dx,
00558                                  const ProblemDomain& a_domain,
00559                                  int                  a_ivar,
00560                                  bool                 a_doSymmetric = false);
00561 
00562   /**
00563      split matrix computations for least squares into separate
00564      routine that is independent of stencil size, so that the
00565      code can be reused
00566    */
00567   static void
00568   calculateWeightingMatrix(RealVect           x0,
00569                            Vector<RealVect>&  xp,
00570                            Vector<RealVect>&  weightMatrix,
00571                            bool&              detZero);
00572 
00573   static void
00574   calculateWeightingMatrixRed(RealVect           x0,
00575                               Vector<RealVect>&  xp,
00576                               IntVect            dimm,
00577                               Vector<RealVect>&  weightMatrix,
00578                               bool&              deadRed);
00579 
00580   ///gets the normal of a domain face
00581   static RealVect  getDomainNormal(int a_idir, Side::LoHiSide a_side);
00582 
00583   ///gets the location in real space of a face center
00584   static RealVect  getFaceLocation(const FaceIndex& a_face,
00585                                    const RealVect&  a_dx,
00586                                    const RealVect&  a_probLo);
00587 
00588 
00589   ///gets the location in real space of a face center
00590   static RealVect  getFaceLocation(const FaceIndex& a_face,
00591                                    const Real&  a_dx,
00592                                    const RealVect&  a_probLo)
00593     {
00594       return getFaceLocation(a_face, a_dx*RealVect::Unit, a_probLo);
00595     }
00596 
00597   ///get stencils for points along a line (generalized johansen stencil)
00598   static void dataRayCast(bool&               a_dropOrder,
00599                           Vector<VoFStencil>& a_pointStencils,
00600                           Vector<Real>&       a_distanceAlongLine,
00601                           const RealVect&     a_normal,
00602                           const RealVect&     a_bndryCentroid,
00603                           const VolIndex&     a_vof,
00604                           const EBISBox&      a_ebisBox,
00605                           const RealVect&     a_dx,
00606                           const IntVectSet&   a_cfivs,
00607                           int a_ivar,
00608                           int a_numPoints);
00609 
00610   ///gets the location in real space of a cell center
00611   static RealVect  getVofLocation(const VolIndex&  a_vof,
00612                                   const RealVect&  a_dx,
00613                                   const RealVect&  a_probLo);
00614 
00615   ///I have misspelled this one time too many
00616   static RealVect  getVoFLocation(const VolIndex&  a_vof,
00617                                   const RealVect&  a_dx,
00618                                   const RealVect&  a_probLo)
00619   {
00620     return getVofLocation(a_vof, a_dx, a_probLo);
00621   }
00622 
00623   ///I have misspelled this one time too many
00624   static RealVect  getVoFLocation(const VolIndex&  a_vof,
00625                                   const Real&  a_dx,
00626                                   const RealVect&  a_probLo)
00627   {
00628     RealVect rvdx = a_dx*RealVect::Unit;
00629     return getVofLocation(a_vof, rvdx, a_probLo);
00630   }
00631 
00632 
00633   ///I have misspelled this one time too many
00634   static RealVect  getVofLocation(const VolIndex&  a_vof,
00635                                   const Real&  a_dx,
00636                                   const RealVect&  a_probLo)
00637   {
00638     RealVect rvdx = a_dx*RealVect::Unit;
00639     return getVofLocation(a_vof, rvdx, a_probLo);
00640   }
00641 
00642 
00643   ///gets the location in real space of a cell center
00644   static RealVect  getIVLocation(const IntVect&   a_iv,
00645                                  const RealVect&  a_dx,
00646                                  const RealVect&  a_probLo);
00647 
00648   ///
00649   /**
00650      A very easy to screw up piece of code that was in four places.
00651    */
00652   static void
00653   defineFluxInterpolant(LevelData<BaseIFFAB<Real> >& a_fluxInterpolant,
00654                         LayoutData<IntVectSet>     & a_irregSetsGrown,
00655                         const DisjointBoxLayout    & a_dbl,
00656                         const EBISLayout           & a_ebisl,
00657                         const ProblemDomain        & a_domain,
00658                         const int                  & a_ncomp,
00659                         const int                  & a_faceDir);
00660 
00661   ///
00662   static void
00663   interpolateFluxToCentroids(LevelData<EBFluxFAB>&       a_centroidFlux,
00664                              const LevelData<EBFluxFAB>& a_faceCentFlux,
00665                              const DisjointBoxLayout&    a_grids,
00666                              const EBISLayout&           a_ebisl,
00667                              const ProblemDomain&        a_domain);
00668 
00669   ///
00670   static void
00671   interpolateFluxToCentroids(EBFaceFAB&                  a_centroidFlux,
00672                              const EBFaceFAB&            a_faceCentFlux,
00673                              const Box&                  a_box,
00674                              const EBISBox&              a_ebisBox,
00675                              const ProblemDomain&        a_domain,
00676                              const int&                  a_idir);
00677 
00678 
00679   ///
00680   static FaceStencil getInterpStencil(const FaceIndex&     a_face,
00681                                       const IntVectSet&    a_coarseFineIVS,
00682                                       const EBISBox&       a_ebisBox,
00683                                       const ProblemDomain& a_domainBox);
00684 
00685   ///
00686   static void getInterpStencil2D(FaceStencil&          a_sten,
00687                                  const FaceIndex&      a_face,
00688                                  const IntVectSet&     a_coarseFineIVS,
00689                                  const EBISBox&        a_ebisBox,
00690                                  const ProblemDomain&  a_domainBox) ;
00691 
00692   ///
00693   static void getInterpStencil3D(FaceStencil&         a_sten,
00694                                  const FaceIndex&     a_face,
00695                                  const IntVectSet&    a_coarseFineIVS,
00696                                  const EBISBox&       a_ebisBox,
00697                                  const ProblemDomain& a_domainBox);
00698 
00699   ///
00700   /**
00701      Return true if there is unique vof associated with a_cell
00702      that lies in the a_vofsStencil. If so, a_vof = vof. Else return false.
00703   */
00704   static bool
00705   isVoFHere(VolIndex& a_vof,
00706             const Vector<VolIndex>& a_vofsStencil,
00707             const IntVect& a_cell );
00708 
00709   ///
00710   static bool
00711   isVoFHere(VolIndex& a_vof, int& a_whichVoF,
00712             const Vector<VolIndex>& a_vofsStencil,
00713             const IntVect& a_cell );
00714 
00715   ///
00716   /** Returns all vofs that can be reached from a_vof via a
00717       monotone path of length <= than radius
00718   */
00719   static void
00720   getAllVoFsInMonotonePath(Vector<VolIndex>& a_vofList,
00721                            const VolIndex&   a_vof,
00722                            const EBISBox&    a_ebisBox,
00723                            const int&        a_redistRad);
00724 
00725   /*******/
00726   static void
00727   getAllVoFsWithinRadius(Vector<VolIndex>& a_vofList,
00728                          Vector<IntVect>    & a_dist,
00729                          const VolIndex&   a_vof,
00730                          const EBISBox&    a_ebisBox,
00731                          const int&        a_redistRad);
00732 
00733   ///get volume of a vof in rz coords
00734   static void
00735   getKVolRZ(Real&           a_kvol,
00736             Real&           a_cellVol,
00737             const EBISBox&  a_ebisBox,
00738             const Real&     a_dx,
00739             const VolIndex& a_vof);
00740 
00741   ///get volume of a vof in rz coords
00742   static void
00743   getKVolRZNoDx(Real&           a_kvol,
00744                 Real&           a_cellVolDx3,
00745                 const EBISBox&  a_ebisBox,
00746                 const VolIndex& a_vof);
00747 
00748   ///
00749   /**
00750      Return true if there is a unique adjacent face in the
00751      given direction.  False if either covered or multivalued.
00752    */
00753   static
00754   bool
00755   getAdjacentFace(FaceIndex& a_adjacentFace,
00756                   const FaceIndex& a_face,
00757                   const EBISBox& a_ebisBox,
00758                   const ProblemDomain& a_domain,
00759                   const int& a_idir,
00760                   const Side::LoHiSide& a_side);
00761 
00762   ///
00763   /**
00764      compute a stencil for one of our fancy interpolated
00765      gradients.
00766    */
00767   static
00768   void computeGradFluxStencil(VoFStencil& a_thisStencil,
00769                               const FaceIndex& a_thisFace,
00770                               const EBISBox& a_ebisBox,
00771                               const ProblemDomain& a_domain,
00772                               const int& a_dir);
00773 
00774   ///
00775   /**
00776      compute a stencil for one of our fancy interpolated
00777      fluxes.
00778    */
00779   static
00780   void computeInterpStencil(FaceStencil& a_thisStencil,
00781                             const FaceIndex& a_thisFace,
00782                             const EBISBox& a_ebisBox,
00783                             const ProblemDomain& a_domain,
00784                             const int& a_dir);
00785 
00786   ///
00787   /**
00788      return lp-norm of component comp of a_src,
00789      weighted by local volume fraction.
00790      not normalized by number of points or anything like that.
00791      if p==0, volume returned is one
00792      and norm returned is Max(abs(a_src)) over
00793      uncovered regions.  otherwise,
00794      returns sum(volfrac*a_src(iv,comp)**p)  of component comp of a_src
00795      weighted by local volume fraction and
00796      also returns volume of uncovered regions.
00797      Only uncovered regions count here.
00798    */
00799   static
00800   void
00801   volWeightedSum(Real&              a_norm,
00802                  Real&              a_volume,
00803                  const EBCellFAB&   a_src,
00804                  const Box&         a_region,
00805                  const EBISBox&     a_ebisBox,
00806                  const RealVect&    a_dx,
00807                  const IntVectSet&  a_ivsExclude,
00808                  const int&         a_comp,
00809                  const int&         a_p,
00810                  EBNormType::NormMode = EBNormType::OverBoth);
00811 
00812 
00813   static void
00814   getCompVolRZ(Real&           a_compVol,
00815                const EBISBox&  a_ebisBox,
00816                const Real&     a_dx,
00817                const VolIndex& a_vof,
00818                bool a_verbose = false);
00819   ///
00820   /**
00821      return volume-weighted sum of component comp of a_src,
00822      weighted by local volume fraction.
00823      and norm returned is Max(abs(a_src)) over
00824      uncovered regions.  otherwise,
00825      returns sum(volfrac*a_src(iv,comp)**p)  of component comp of a_src
00826      weighted by local volume fraction and
00827      also returns volume of uncovered regions.
00828      Only uncovered regions count here.
00829    */
00830   static
00831   void
00832   computeSum(Real&              a_norm,
00833              Real&              a_volume,
00834              const EBCellFAB&   a_src,
00835              const Box&         a_region,
00836              const EBISBox&     a_ebisBox,
00837              const RealVect&    a_dx,
00838              const IntVectSet&  a_ivsExclude,
00839              const int&         a_comp,
00840              EBNormType::NormMode = EBNormType::OverBoth);
00841 
00842 
00843   ///
00844   /**
00845      return volume-weighted sum of component comp of a_src,
00846      weighted by local volume fraction.
00847      and norm returned is Max(abs(a_src)) over
00848      uncovered regions.  otherwise,
00849      returns sum(volfrac*a_src(iv,comp)**p)  of component comp of a_src
00850      weighted by local volume fraction and
00851      also returns volume of uncovered regions.
00852      Only uncovered regions count here.
00853    */
00854   static
00855   void
00856   computeUnweightedSum(Real&             a_norm,
00857                        Real&             a_volume,
00858                        const EBCellFAB&  a_src,
00859                        const Box&        a_region,
00860                        const EBISBox&    a_ebisBox,
00861                        const RealVect&   a_dx,
00862                        const IntVectSet& a_ivsExclude,
00863                        const int&        a_comp,
00864                        EBNormType::NormMode = EBNormType::OverBoth);
00865 
00866   ///
00867   /** return l-p norm of a_src.
00868       if p==0, v
00869       norm returned is Max(abs(a_src)) over
00870       uncovered regions.  otherwise,
00871       returns 1/vol(sum(volfrac*a_src(iv,comp)**p)^(1/p))
00872       of component comp of a_src
00873       weighted by local volume fraction and
00874       also returns volume of uncovered regions.
00875       Only uncovered regions count here.
00876       The data must have the same layout as a_layout with
00877       the possible exception of ghost cells.
00878 
00879       int pmode = -2;  //norm = (1/v)(sum(phi dv)) ---no absolute values and multiply kappa as you go
00880       int pmode = -1;  //norm = (1/v)(sum(phi dv)) ---no absolute values and assume kappa already multiplied in
00881 
00882    */
00883   static
00884   Real
00885   norm(const BoxLayoutData<EBCellFAB >& a_dataOne,
00886        const BoxLayout& a_layout,
00887        const EBISLayout& a_ebisl,
00888        const int& comp,
00889        const int&  p,
00890        EBNormType::NormMode = EBNormType::OverBoth);
00891 
00892   static
00893   Real
00894   norm(const Vector<LevelData<EBCellFAB>* >& a_data,
00895        const Vector< DisjointBoxLayout > &   a_layout,
00896        const Vector< EBISLayout >&           a_ebisl,
00897        const Vector< int >&                  a_refRatio,
00898        const int& comp,
00899        const int&  p,
00900        EBNormType::NormMode = EBNormType::OverBoth);
00901 
00902   ///
00903   static
00904   void
00905   computeCoveredFaces(Vector<VolIndex>&     a_coveredFace,
00906                       IntVectSet&           a_coveredSets,
00907                       IntVectSet&           a_irregIVS,
00908                       const int&            a_idir,
00909                       const Side::LoHiSide& a_sd,
00910                       const EBISBox&        a_ebisBox,
00911                       const Box&            a_region);
00912 
00913   //testing interface
00914   static
00915   Real
00916   norm(Real& volume, const BoxLayoutData<EBCellFAB >& a_dataOne,
00917        const BoxLayout& a_layout,
00918        const EBISLayout& a_ebisl,
00919        const int& comp,
00920        const int&  p,
00921        EBNormType::NormMode = EBNormType::OverBoth);
00922 
00923   ///
00924   /** return l-p norm of a_src.
00925       if p==0, v
00926       norm returned is Max(abs(a_src)) over
00927       uncovered regions.  otherwise,
00928       returns 1/vol(sum(volfrac*a_src(iv,comp)**p)^(1/p))
00929       of component comp of a_src
00930       weighted by local volume fraction and
00931       also returns volume of uncovered regions.
00932       Only uncovered regions count here.
00933       The data must have the same layout as a_layout with
00934       the possible exception of ghost cells.
00935    */
00936   static
00937   Real
00938   norm(const EBCellFAB& a_dataOne,
00939        const Box&       a_grid,
00940        const EBISBox&   a_ebisl,
00941        const int&       a_comp,
00942        const int&       a_p,
00943        EBNormType::NormMode = EBNormType::OverBoth);
00944 
00945   static
00946   Real
00947   norm(Real& volume, const EBCellFAB& a_dataOne,
00948        const Box&       a_grid,
00949        const EBISBox&   a_ebisl,
00950        const int&       a_comp,
00951        const int&       a_p,
00952        EBNormType::NormMode = EBNormType::OverBoth);
00953 
00954   ///
00955   /** return l-p norm of a_src.
00956       if p==0, v
00957       norm returned is Max(abs(a_src)) over
00958       uncovered regions.  otherwise,
00959       returns sum(volfrac*a_src(iv,comp)**p)  of component comp of a_src
00960       weighted by local volume fraction and
00961       also returns volume of uncovered regions.
00962       Only uncovered regions count here.
00963       The data must have the same layout as a_layout with
00964       the possible exception of ghost cells.
00965    */
00966   static
00967   void
00968   volWeightedSum(Real&                            a_sum,
00969                  Real&                            a_volume,
00970                  const BoxLayoutData<EBCellFAB >& a_dataOne,
00971                  const BoxLayout&                 a_layout,
00972                  const EBISLayout&                a_ebisl,
00973                  const RealVect&                  a_dx,
00974                  const IntVectSet&                a_ivsExclude,
00975                  const int&                       a_comp,
00976                  const int&                       a_p,
00977                  EBNormType::NormMode = EBNormType::OverBoth);
00978 
00979   ///
00980   /** return the sum of all irregular boundary areas
00981    */
00982   static
00983   void
00984   sumBndryArea(Real&               a_area,
00985                const BoxLayout&    a_region,
00986                const EBISLayout&   a_ebisl);
00987 
00988   ///
00989   /** return the dotproduct of two leveldatas of ebfabs,
00990       Only uncovered regions count here.
00991    */
00992   static
00993   Real
00994   dotProduct(const BoxLayoutData<EBCellFAB >& a_dataOne,
00995              const BoxLayoutData<EBCellFAB >& a_dataTwo,
00996              const BoxLayout& a_layout,
00997              const EBISLayout& a_ebisl,
00998              const int& a_comp);
00999 
01000   ///
01001   /** return the dotproduct of two  ebfabs,
01002       Only uncovered regions count here.
01003    */
01004   static
01005   Real
01006   dotProduct(const EBCellFAB& a_dataOne,
01007              const EBCellFAB& a_dataTwo,
01008              const Box& a_layout,
01009              const EBISBox& a_ebisBox,
01010              const int& a_comp);
01011 
01012   ///
01013   /** Given a VoF, a_vof1, and a cell, a_cell2, determine if there is a vector of
01014       VoFs in a_cell2 that connects to a_vof1 via a monotone path.
01015       A vector of VolIndices is returned, even if empty.
01016    */
01017   static
01018   void
01019   monotonePathVoFToCellMultiVoFs(Vector<VolIndex>& a_vofsInCell,
01020                                  const Vector<VolIndex>& a_vofsInPath,
01021                                  const IntVect& a_cell2,
01022                                  const EBISBox& a_ebisBox);
01023   ///
01024   /** Given a VoF, a_vof1, and a cell, a_cell2, determine if there is a single
01025       VoF in a_cell2 that connects to a_vof1 via a monotone path.  If there is
01026       one such VoF then TRUE is returned (and the VolIndex is returned,
01027       a_vof2).  If there is no such VoF or if there are more than one such
01028       VoF then FALSE is returned (and a_vof2 is unchanged).
01029    */
01030   static
01031   bool
01032   monotonePathVoFToCellVoF(VolIndex& a_vof2,
01033                            const VolIndex& a_vof1,
01034                            const IntVect& a_cell2,
01035                            const EBISBox& a_ebisBox);
01036   /*******/
01037   static void
01038   getAllVoFsWithinRadius(Vector<VolIndex>& a_vofList,
01039                          const VolIndex&   a_vof,
01040                          const EBISBox&    a_ebisBox,
01041                          const int&        a_redistRad);
01042 
01043   /*******/
01044   static void
01045   getAllVoFsWithinRadius(Vector<VolIndex>& a_vofList,
01046                          const IntVect&    a_timesMoved,
01047                          const VolIndex&   a_vof,
01048                          const EBISBox&    a_ebisBox,
01049                          const int&        a_redistRad);
01050 
01051   ///
01052   /** Returns all vofs that can be reached from a_vof via a
01053       monotone path of length <= than radius
01054   */
01055   static
01056   void
01057   getAllVoFsInMonotonePath(Vector<VolIndex>& a_vofs,
01058                            const IntVect&    a_timesMoved,
01059                            const IntVect&    a_pathSign,
01060                            const VolIndex&   a_vof,
01061                            const EBISBox&    a_ebisBox,
01062                            const int&        a_radius);
01063 
01064   static VolIndex& getVoFMax()
01065   {
01066     return s_vofMax;
01067   }
01068   static Real&     getValMax()
01069   {
01070     return s_valMax;
01071   }
01072 
01073   static void setMinVolumeFraction(Real a_minVolFrac)
01074   {
01075     s_minVolFrac = a_minVolFrac;
01076   }
01077 
01078 private:
01079   static VolIndex s_vofMax;
01080   static Real     s_valMax;
01081   static Real     s_minVolFrac;
01082 
01083 };
01084 
01085 extern IntVect ebcoarsen (const IntVect& b,
01086                           int        refinement_ratio);
01087 
01088 extern
01089 Box ebrefine (const Box& b,
01090               int        refinement_ratio);
01091 extern
01092 ProblemDomain ebrefine (const ProblemDomain& b,
01093                         int        refinement_ratio);
01094 
01095 extern Box ebcoarsen (const Box& b,
01096                       int        refinement_ratio);
01097 extern ProblemDomain ebcoarsen (const ProblemDomain& b,
01098                                 int        refinement_ratio);
01099 
01100 extern void ebrefine(DisjointBoxLayout& output,
01101                      const DisjointBoxLayout& input,
01102                      int refinement);
01103 
01104 extern void ebcoarsen(DisjointBoxLayout& output,
01105                       const DisjointBoxLayout& input,
01106                       int refinement);
01107 
01108 #include "NamespaceFooter.H"
01109 #endif