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