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