Chombo + EB: NWOViscousTensorOp.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 _NWOVISCOUSTENSOROP_H_
00012 #define _NWOVISCOUSTENSOROP_H_
00013 
00014 #include "AMRMultiGrid.H"
00015 #include "REAL.H"
00016 #include "Box.H"
00017 #include "LevelDataOps.H"
00018 #include "BCFunc.H"
00019 #include "FArrayBox.H"
00020 #include "FluxBox.H"
00021 #include "CFIVS.H"
00022 #include "AMRTGA.H"
00023 #include "NWOQuadCFInterp.H"
00024 #include "FourthOrderCFInterp.H"
00025 #include "CoarseAverage.H"
00026 #include "LevelFluxRegister.H"
00027 #include "AMRIO.H"
00028 #include "NamespaceHeader.H"
00029 
00030 #define VTOP_DEFAULT_SAFETY 0.9
00031 ///
00032 /**
00033    operator is alpha a I + (divF) = alpha*acoef I + beta*div(eta(grad B + grad BT) + lambda I div B )
00034    beta, lambda,  and eta are incorporated into the flux F
00035 */
00036 class NWOViscousTensorOp : public LevelTGAHelmOp<LevelData<FArrayBox>, FluxBox>
00037 {
00038 public:
00039 
00040   enum prolongationType
00041     {
00042       piecewiseConstant = 0,
00043       linearInterp,
00044       NUM_INTERP_TYPE
00045     };
00046 
00047   virtual void diagonalScale(LevelData<FArrayBox>& a_rhs, bool a_kappaWeighted)
00048   {
00049     diagonalScale(a_rhs);
00050   }
00051 
00052   virtual void setAlphaAndBeta(const Real& a_alpha,
00053                                const Real& a_beta)
00054   {
00055     m_alpha = a_alpha;
00056     m_beta = a_beta;
00057     defineRelCoef();
00058   }
00059 
00060   virtual void diagonalScale(LevelData<FArrayBox>& a_rhs)
00061   {
00062     DisjointBoxLayout grids = a_rhs.disjointBoxLayout();
00063     for (DataIterator dit = grids.dataIterator(); dit.ok(); ++dit)
00064       {
00065         for (int idir = 0; idir < SpaceDim; idir++)
00066           {
00067             int isrc = 0; int idst = idir; int inco = 1;
00068             a_rhs[dit()].mult((*m_acoef)[dit()], isrc, idst, inco);
00069           }
00070       }
00071   }
00072 
00073   virtual void divideByIdentityCoef(LevelData<FArrayBox>& a_rhs)
00074   {
00075     DisjointBoxLayout grids = a_rhs.disjointBoxLayout();
00076     for (DataIterator dit = grids.dataIterator(); dit.ok(); ++dit)
00077       {
00078         for (int idir = 0; idir < SpaceDim; idir++)
00079           {
00080             int isrc = 0; int idst = idir; int inco = 1;
00081             a_rhs[dit()].divide((*m_acoef)[dit()], isrc, idst, inco);
00082           }
00083       }
00084   }
00085 
00086   ///
00087   /**
00088    */
00089   virtual ~NWOViscousTensorOp()
00090   {
00091   }
00092 
00093   ///
00094   /**
00095    */
00096   NWOViscousTensorOp(const DisjointBoxLayout&                    a_grids,
00097                      const DisjointBoxLayout&                    a_gridsFine,
00098                      const DisjointBoxLayout&                    a_gridsCoar,
00099                      const RefCountedPtr<LevelData<FluxBox> >&   a_eta,
00100                      const RefCountedPtr<LevelData<FluxBox> >&   a_lambda,
00101                      const RefCountedPtr<LevelData<FArrayBox> >& a_acoef,
00102                      Real                                        a_alpha,
00103                      Real                                        a_beta,
00104                      int                                         a_refToFine,
00105                      int                                         a_refToCoar,
00106                      const ProblemDomain&                        a_domain,
00107                      const Real&                                 a_dxLevel,
00108                      const Real&                                 a_dxCoar,
00109                      BCFunc                                      a_bc,
00110                      Real                                        a_safety = VTOP_DEFAULT_SAFETY,
00111                      Real                                        a_relaxTolerance = 0.0,
00112                      int                                         a_relaxMinIter = 1000
00113                      );
00114 
00115   ///
00116   /**
00117    */
00118   NWOViscousTensorOp(const DisjointBoxLayout&                    a_grids,
00119                      const DisjointBoxLayout&                    a_gridsFine,
00120                      const DisjointBoxLayout&                    a_gridsCoar,
00121                      const RefCountedPtr<LevelData<FluxBox> >&   a_eta,
00122                      const RefCountedPtr<LevelData<FluxBox> >&   a_lambda,
00123                      const RefCountedPtr<LevelData<FArrayBox> >& a_acoef,
00124                      Real                                        a_alpha,
00125                      Real                                        a_beta,
00126                      int                                         a_refToFine,
00127                      int                                         a_refToCoar,
00128                      const ProblemDomain&                        a_domain,
00129                      const Real&                                 a_dxLevel,
00130                      const Real&                                 a_dxCoar,
00131                      BCHolder&                                   a_bc,
00132                      Real                                        a_safety = VTOP_DEFAULT_SAFETY,
00133                      Real                                        a_relaxTolerance = 0.0,
00134                      int                                         a_relaxMinIter = 1000
00135                      );
00136 
00137   ///
00138   /**
00139      Fillgrad here is a placeholder that does not do anything.  The whole fillGrad paradigm is an artifact 
00140      of the how old ViscousTensorOp worked.  Sometimes performance tuning eventually forces you to do the
00141      right thing.
00142    */
00143   void fillGrad(const LevelData<FArrayBox>& a_phiFine)
00144   {
00145   }
00146 
00147   static void loHiCenterFace(Box&                 a_loBox,
00148                              int&                 a_hasLo,
00149                              Box&                 a_hiBox,
00150                              int&                 a_hasHi,
00151                              Box&                 a_centerBox,
00152                              const ProblemDomain& a_eblg,
00153                              const Box&           a_inBox,
00154                              const int&           a_dir);
00155   void cellGrad(FArrayBox&             a_gradPhi,
00156                 const  FArrayBox&      a_phi,
00157                 const Box&             a_grid);
00158 
00159   static void getFaceDivAndGrad(FArrayBox&             a_faceDiv,
00160                                 FArrayBox&             a_faceGrad,
00161                                 const FArrayBox&       a_data,
00162                                 const FArrayBox&       a_gradData,
00163                                 const ProblemDomain&   a_domain,
00164                                 const Box&             a_faceBox,
00165                                 const int&             a_faceDir,
00166                                 const Real             a_dx);
00167   ///
00168   /**
00169    */
00170   void define(const DisjointBoxLayout&                    a_grids,
00171               const DisjointBoxLayout&                    a_gridsFine,
00172               const DisjointBoxLayout&                    a_gridsCoar,
00173               const RefCountedPtr<LevelData<FluxBox> >&   a_eta,
00174               const RefCountedPtr<LevelData<FluxBox> >&   a_lambda,
00175               const RefCountedPtr<LevelData<FArrayBox> >& a_acoef,
00176               Real                                        a_alpha,
00177               Real                                        a_beta,
00178               int                                         a_refToFine,
00179               int                                         a_refToCoar,
00180               const ProblemDomain&                        a_domain,
00181               const Real&                                 a_dxLevel,
00182               const Real&                                 a_dxCoar,
00183               BCHolder&                                   a_bc,
00184               Real                                        a_safety = VTOP_DEFAULT_SAFETY,
00185               Real                                        a_relaxTolerance = 0.0,
00186               int                                         a_relaxMinIter = 1000
00187               );
00188 
00189   virtual void residual(  LevelData<FArrayBox>& a_lhs,
00190                           const LevelData<FArrayBox>& a_phi,
00191                           const LevelData<FArrayBox>& a_rhs,
00192                           bool a_homogeneous = false);
00193 
00194   virtual void preCond(   LevelData<FArrayBox>& a_correction,
00195                           const LevelData<FArrayBox>& a_residual);
00196 
00197   virtual void applyOp(   LevelData<FArrayBox>& a_lhs,
00198                           const LevelData<FArrayBox>& a_phi,
00199                           bool a_homogeneous = false);
00200   virtual void create(    LevelData<FArrayBox>& a_lhs,
00201                           const LevelData<FArrayBox>& a_rhs);
00202   virtual void createCoarsened(    LevelData<FArrayBox>& a_lhs,
00203                                    const LevelData<FArrayBox>& a_rhs,
00204                                    const int& a_refRat);
00205 
00206   void reflux(const LevelData<FArrayBox>& a_phiFine,
00207               const LevelData<FArrayBox>& a_phi,
00208               LevelData<FArrayBox>& residual,
00209               AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00210 
00211   //use at thy peril
00212   virtual void getFlux(FluxBox&             a_flux,
00213                        const LevelData<FArrayBox>&  a_data,
00214                        const Box&       a_grid,
00215                        const DataIndex& a_dit,
00216                        Real             a_scale)
00217   {
00218     const FArrayBox& data = a_data[a_dit];
00219     a_flux.define(a_grid, SpaceDim);
00220     for (int idir = 0; idir < SpaceDim; idir++)
00221       {
00222         const FArrayBox& etaFace     =    (*m_eta)[a_dit][idir];
00223         const FArrayBox& lambdaFace  = (*m_lambda)[a_dit][idir];
00224         Box faceBox = a_flux[idir].box();
00225         Box domFaceBox = surroundingNodes(m_domain.domainBox(), idir);
00226         faceBox &= domFaceBox;
00227         getFlux(a_flux[idir], data, etaFace, lambdaFace, faceBox, idir, 1);
00228         a_flux[idir] *= a_scale;
00229       }
00230   }
00231   void applyOpNoBoundary( LevelData<FArrayBox>& a_lhs,
00232                           const LevelData<FArrayBox>& a_phi)
00233   {
00234     computeOperatorNoBCs(a_lhs, a_phi);
00235   }
00236 
00237   void getFlux(FArrayBox&       a_flux,
00238                const FArrayBox& a_data,
00239                const FArrayBox& a_etaFace,
00240                const FArrayBox& a_lambdaFace,
00241                const Box&       a_facebox,
00242                int              a_dir,
00243                int ref = 1);
00244 
00245   static void
00246   getFlux(FArrayBox&       a_flux,
00247           const FArrayBox& a_data,
00248           const FArrayBox& a_etaFace,
00249           const FArrayBox& a_lamFace,
00250           const Box&       a_faceBox,
00251           const ProblemDomain& a_domain,
00252           const Real         & a_dx,
00253           const Real         & a_beta,
00254           int              a_dir,
00255           int a_ref);
00256 
00257   /// utility function which computes operator after all bc's have been set
00258   void computeOperatorNoBCs(LevelData<FArrayBox>& a_lhs,
00259                             const LevelData<FArrayBox>& a_phi);
00260 
00261   virtual void assign(    LevelData<FArrayBox>& a_lhs,
00262                           const LevelData<FArrayBox>& a_rhs) ;
00263   virtual Real dotProduct(const LevelData<FArrayBox>& a_1,
00264                           const LevelData<FArrayBox>& a_2) ;
00265   virtual void incr( LevelData<FArrayBox>& a_lhs,
00266                      const LevelData<FArrayBox>& a_x,
00267                      Real a_scale) ;
00268   virtual void axby( LevelData<FArrayBox>& a_lhs,
00269                      const LevelData<FArrayBox>& a_x,
00270                      const LevelData<FArrayBox>& a_y,
00271                      Real a, Real b) ;
00272 
00273   virtual void scale(LevelData<FArrayBox>& a_lhs, const Real& a_scale) ;
00274 
00275   virtual Real norm(const LevelData<FArrayBox>& a_x, int a_ord);
00276 
00277   virtual Real dx() const
00278   {
00279     return m_dx;
00280   }
00281 
00282   virtual Real dxCrse() const
00283   {
00284     return m_dxCrse;
00285   }
00286 
00287   virtual void setToZero( LevelData<FArrayBox>& a_x);
00288   /*@}*/
00289 
00290   /**
00291      \name MGLevelOp functions */
00292   /*@{*/
00293 
00294   virtual void relax(LevelData<FArrayBox>& a_e,
00295                      const LevelData<FArrayBox>& a_residual,
00296                      int iterations);
00297 
00298   virtual void createCoarser(LevelData<FArrayBox>& a_coarse,
00299                              const LevelData<FArrayBox>& a_fine,
00300                              bool ghosted);
00301   /**
00302      calculate restricted residual
00303      a_resCoarse[2h] = I[h->2h] (rhsFine[h] - L[h](phiFine[h])
00304   */
00305   virtual void restrictResidual(LevelData<FArrayBox>& a_resCoarse,
00306                                 LevelData<FArrayBox>& a_phiFine,
00307                                 const LevelData<FArrayBox>& a_rhsFine);
00308 
00309   /**
00310      correct the fine solution based on coarse correction
00311      a_phiThisLevel += I[2h->h](a_correctCoarse)
00312   */
00313   virtual void prolongIncrement(LevelData<FArrayBox>& a_phiThisLevel,
00314                                 const LevelData<FArrayBox>& a_correctCoarse);
00315 
00316   /*@}*/
00317 
00318   /**
00319      \name AMRLevelOp functions */
00320   /*@{*/
00321 
00322   /** returns 1 when there are no coarser AMRLevelOp objects */
00323   virtual int refToCoarser()
00324   {
00325     return m_refToCoar;
00326   }
00327 
00328   /** a_residual = a_rhs - L(a_phi, a_phiFine, a_phiCoarse) */
00329   virtual void AMRResidual(LevelData<FArrayBox>& a_residual,
00330                            const LevelData<FArrayBox>& a_phiFine,
00331                            const LevelData<FArrayBox>& a_phi,
00332                            const LevelData<FArrayBox>& a_phiCoarse,
00333                            const LevelData<FArrayBox>& a_rhs,
00334                            bool a_homogeneousPhysBC,
00335                            AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00336 
00337   /** residual assuming no more coarser AMR levels */
00338 
00339   virtual void AMRResidualNC(LevelData<FArrayBox>& a_residual,
00340                              const LevelData<FArrayBox>& a_phiFine,
00341                              const LevelData<FArrayBox>& a_phi,
00342                              const LevelData<FArrayBox>& a_rhs,
00343                              bool a_homogeneousPhysBC,
00344                              AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00345 
00346   /** a_residual = a_rhs - L(a_phi, a_phiCoarse)  */
00347   virtual void AMRResidualNF(LevelData<FArrayBox>& a_residual,
00348                              const LevelData<FArrayBox>& a_phi,
00349                              const LevelData<FArrayBox>& a_phiCoarse,
00350                              const LevelData<FArrayBox>& a_rhs,
00351                              bool a_homogeneousPhysBC);
00352 
00353   /** a_resCoarse = I[h-2h]( a_residual - L(a_correction, a_coarseCorrection))
00354       it is assumed that a_resCoarse has already been filled in with the coarse
00355       version of AMRResidualNF and that this operation is free to overwrite
00356       in the overlap regions.
00357   */
00358 
00359   virtual void AMRRestrict(LevelData<FArrayBox>& a_resCoarse,
00360                            const LevelData<FArrayBox>& a_residual,
00361                            const LevelData<FArrayBox>& a_correction,
00362                            const LevelData<FArrayBox>& a_coarseCorrection, 
00363                            bool a_skip_res = false );
00364 
00365   /** a_correction += I[h->h](a_coarseCorrection) */
00366   virtual void AMRProlong(LevelData<FArrayBox>& a_correction,
00367                           const LevelData<FArrayBox>& a_coarseCorrection);
00368 
00369   /** a_residual = a_residual - L(a_correction, a_coarseCorrection) */
00370   virtual void AMRUpdateResidual(LevelData<FArrayBox>& a_residual,
00371                                  const LevelData<FArrayBox>& a_correction,
00372                                  const LevelData<FArrayBox>& a_coarseCorrection);
00373 
00374   ///
00375   /**
00376      compute norm over all cells on coarse not covered by finer
00377   */
00378   virtual Real AMRNorm(const LevelData<FArrayBox>& a_coarseResid,
00379                        const LevelData<FArrayBox>& a_fineResid,
00380                        const int&                  a_refRat,
00381                        const int&                  a_ord);
00382 
00383   virtual void outputLevel(LevelData<FArrayBox>& a_rhs,
00384                            string& a_name)
00385   {
00386     string fname(a_name);
00387     fname.append(".hdf5");
00388 #ifdef CH_HDF5
00389     writeLevelname(&a_rhs, fname.c_str());
00390 #endif
00391   }
00392 
00393   //debugging hook
00394   virtual void outputAMR(Vector<LevelData<FArrayBox>* >& a_rhs, string& a_name)
00395   {
00396     // for now, klugs is to use outputLevel on level 0
00397     outputLevel(*a_rhs[0], a_name);
00398   }
00399 
00400   void homogeneousCFInterp(LevelData<FArrayBox>& a_phif);
00401   ///
00402   /**
00403      does homogeneous coarse/fine interpolation for phi
00404   */
00405   void
00406   homogeneousCFInterpPhi(LevelData<FArrayBox>& a_phif,
00407                          const DataIndex& a_datInd,
00408                          int a_idir,
00409                          Side::LoHiSide a_hiorlo);
00410 
00411   ///
00412   /** does homogeneous coarse/fine interpolation for tangential gradient
00413       (needs phi ghost cells to be filled in first, so should call
00414       homogeneousCFInterpPhi first)
00415   */
00416   void homogeneousCFInterpTanGrad(LevelData<FArrayBox>& a_tanGrad,
00417                                   const LevelData<FArrayBox>& a_phi,
00418                                   const DataIndex& a_datInd,
00419                                   int a_idir,
00420                                   Side::LoHiSide a_hiorlo);
00421 
00422   void interpOnIVSHomo(LevelData<FArrayBox>& a_phif,
00423                        const DataIndex& a_datInd,
00424                        const int a_idir,
00425                        const Side::LoHiSide a_hiorlo,
00426                        const IntVectSet& a_interpIVS);
00427   ///
00428   /**
00429      Apply the AMR operator, including coarse-fine matching
00430   */
00431   void AMROperator(      LevelData<FArrayBox>& a_LofPhi,
00432                          const LevelData<FArrayBox>& a_phiFine,
00433                          const LevelData<FArrayBox>& a_phi,
00434                          const LevelData<FArrayBox>& a_phiCoarse,
00435                          bool a_homogeneousDomBC,
00436                          AMRLevelOp<LevelData<FArrayBox> >*  a_finerOp);
00437 
00438   void AMROperatorNF(      LevelData<FArrayBox>& a_LofPhi,
00439                            const LevelData<FArrayBox>& a_phi,
00440                            const LevelData<FArrayBox>& a_phiCoarse,
00441                            bool a_homogeneousBC);
00442 
00443   virtual void AMROperatorNC(      LevelData<FArrayBox>& a_LofPhi,
00444                                    const LevelData<FArrayBox>& a_phiFine,
00445                                    const LevelData<FArrayBox>& a_phi,
00446                                    bool a_homogeneousBC,
00447                                    AMRLevelOp<LevelData<FArrayBox> >* a_finerOp);
00448 
00449   void cfinterp(const LevelData<FArrayBox>& a_phiFine,
00450                 const LevelData<FArrayBox>& a_phiCoarse);
00451 
00452   /// prolongation type
00453   /** make this public/static to make it easy to change, while
00454       also ensuring that it remains constent across all instances
00455   */
00456   static int s_prolongType;
00457 
00458   /// access function
00459   
00460   RefCountedPtr<LevelData<FluxBox>   > getEta()    const {return m_eta;}
00461   RefCountedPtr<LevelData<FluxBox>   > getLambda() const {return m_lambda;}
00462   RefCountedPtr<LevelData<FArrayBox> > getACoef()  const {return m_acoef;}
00463   Real getAlpha() const {return m_alpha;}
00464   Real getBeta()  const {return m_beta;}
00465 
00466 protected:
00467   void defineRelCoef();
00468   RefCountedPtr<LevelData<FluxBox> >         m_eta;
00469   RefCountedPtr<LevelData<FluxBox> >         m_lambda;
00470   RefCountedPtr<LevelData<FArrayBox> >       m_acoef;
00471   Real                                       m_alpha;
00472   Real                                       m_beta;
00473   int                                        m_refToCoar;
00474   int                                        m_refToFine;
00475   BCHolder                                   m_bc;
00476 
00477   DisjointBoxLayout m_grids, m_gridsCoar, m_gridsFine;
00478   
00479   //b is a scalar in 2d, vector in 3d
00480   int m_ncomp;
00481   Real                    m_dx;
00482   Real                    m_dxCrse;
00483   ProblemDomain           m_domain;
00484   //relaxation coef
00485 public: // I want to get this from an owner of this -mfa
00486   LevelData<FArrayBox>    m_relaxCoef;
00487 protected:
00488   // underrelaxation parameter
00489   Real                    m_safety;
00490   // stopping tolerance for relaxation
00491   Real                    m_relaxTolerance;
00492   // minimum number of smooths before tolerance check
00493   int                     m_relaxMinIter;
00494 
00495   LevelDataOps<FArrayBox> m_levelOps;
00496   Copier                  m_exchangeCopier;
00497   NWOQuadCFInterp         m_interpWithCoarser;
00498   //FourthOrderCFInterp         m_interpWithCoarser;
00499 
00500 
00501   Vector<IntVect> m_colors;
00502 
00503 private:
00504   ///weak construction is bad
00505   NWOViscousTensorOp()
00506   {
00507     MayDay::Error("invalid operator");
00508   }
00509   //copy constructor and operator= disallowed for all the usual reasons
00510   NWOViscousTensorOp(const NWOViscousTensorOp& a_opin)
00511   {
00512     MayDay::Error("invalid operator");
00513   }
00514   void operator=(const NWOViscousTensorOp& a_opin)
00515   {
00516     MayDay::Error("invalid operator");
00517   }
00518 
00519 };
00520 
00521 ///
00522 /**
00523    Factory to create NWOViscousTensorOps
00524 */
00525 class NWOViscousTensorOpFactory: public AMRLevelOpFactory<LevelData<FArrayBox> >
00526 {
00527 public:
00528   virtual ~NWOViscousTensorOpFactory()
00529   {
00530 #if 0
00531     for (int ivec = 0; ivec < m_eta.size(); ivec++)
00532       {
00533         m_eta   [ivec] = RefCountedPtr<LevelData<FluxBox> >();
00534         m_lambda[ivec] = RefCountedPtr<LevelData<FluxBox> >();
00535         m_acoef [ivec] = RefCountedPtr<LevelData<FArrayBox> >();
00536       }
00537 #endif
00538 
00539   }
00540 
00541   NWOViscousTensorOpFactory(const Vector<DisjointBoxLayout>&                     a_grids,
00542                             const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_eta,
00543                             const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_lambda,
00544                             const Vector<RefCountedPtr<LevelData<FArrayBox> > >& a_acoef,
00545                             Real                                                 a_alpha,
00546                             Real                                                 a_beta,
00547                             const Vector<int>&                                   a_refRatios,
00548                             const ProblemDomain&                                 a_domainCoar,
00549                             const Real&                                          a_dxCoar,
00550                             BCFunc                                               a_bc,
00551                             Real                                                 a_safety = VTOP_DEFAULT_SAFETY,
00552                             Real                                                 a_relaxTolerance = 0.0,
00553                             int                                                  a_relaxMinIter = 1000
00554                             );
00555 
00556   NWOViscousTensorOpFactory(const Vector<DisjointBoxLayout>&                     a_grids,
00557                             const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_eta,
00558                             const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_lambda,
00559                             const Vector<RefCountedPtr<LevelData<FArrayBox> > >& a_acoef,
00560                             Real                                                 a_alpha,
00561                             Real                                                 a_beta,
00562                             const Vector<int>&                                   a_refRatios,
00563                             const ProblemDomain&                                 a_domainCoar,
00564                             const Real&                                          a_dxCoar,
00565                             BCHolder                                             a_bc,
00566                             Real                                                 a_safety = VTOP_DEFAULT_SAFETY,
00567                             Real                                                 a_relaxTolerance = 0.0,
00568                             int                                                  a_relaxMinIter = 1000
00569                             );
00570 
00571   virtual void define(const Vector<DisjointBoxLayout>&                     a_grids,
00572                       const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_eta,
00573                       const Vector<RefCountedPtr<LevelData<FluxBox> > >&   a_lambda,
00574                       const Vector<RefCountedPtr<LevelData<FArrayBox> > >& a_acoef,
00575                       Real                                                 a_alpha,
00576                       Real                                                 a_beta,
00577                       const Vector<int>&                                   a_refRatios,
00578                       const ProblemDomain&                                 a_domainCoar,
00579                       const Real&                                          a_dxCoar,
00580                       BCHolder                                            a_bc,
00581                       Real                                                 a_safety = VTOP_DEFAULT_SAFETY,
00582                       Real                                                 a_relaxTolerance = 0.0,
00583                       int                                                  a_relaxMinIter = 1000
00584                       );
00585 
00586   ///
00587   virtual NWOViscousTensorOp*
00588   MGnewOp(const ProblemDomain& a_FineindexSpace,
00589           int depth,
00590           bool homoOnly = true);
00591 
00592   ///
00593   virtual  NWOViscousTensorOp* AMRnewOp(const ProblemDomain& a_indexSpace);
00594 
00595   ///
00596   virtual int refToFiner(const ProblemDomain&) const;
00597 
00598   /// set any relevant default values
00599   void setDefaults();
00600 
00601   // 0 = arithmetic, 1 = harmonic
00602   static int s_coefficientAverageType;
00603 
00604 
00605 private:
00606   Vector<RefCountedPtr<LevelData<FluxBox> > >    m_eta;
00607   Vector<RefCountedPtr<LevelData<FluxBox> > >    m_lambda;
00608   Vector<RefCountedPtr<LevelData<FArrayBox> > >  m_acoef;
00609   LevelData<FArrayBox>                           m_phic;
00610   Vector<ProblemDomain>                          m_domains;
00611   Vector<DisjointBoxLayout>                      m_boxes;
00612   Vector<Real>                                   m_dx;
00613   // refinement to next coarser level
00614   Vector<int>                                    m_refRatios;
00615   BCHolder  m_bc;
00616   Real m_alpha;
00617   Real m_beta;
00618   Real m_safety;
00619   Real m_relaxTolerance;
00620   int                     m_relaxMinIter;
00621 
00622   ///weak construction is bad
00623   NWOViscousTensorOpFactory()
00624   {
00625     MayDay::Error("invalid operator");
00626   }
00627   //copy constructor and operator= disallowed for all the usual reasons
00628   NWOViscousTensorOpFactory(const NWOViscousTensorOpFactory& a_opin)
00629   {
00630     MayDay::Error("invalid operator");
00631   }
00632 
00633   void operator=(const NWOViscousTensorOpFactory& a_opin)
00634   {
00635     MayDay::Error("invalid operator");
00636   }
00637 };
00638 
00639 extern void
00640 NWOcoarsenStuff(LevelData<FluxBox> &          a_etaCoar,
00641                 LevelData<FluxBox> &          a_lambdaCoar,
00642                 const LevelData<FluxBox> &    a_etaFine,
00643                 const LevelData<FluxBox> &    a_lambdaFine,
00644                 const int &                   a_refToDepth);
00645 
00646 #include "NamespaceFooter.H"
00647 #endif