Chombo + EB + MF: NewFourthOrderCoordSys.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 _NEWFOURTHORDERCOORDSYS_H_
00012 #define _NEWFOURTHORDERCOORDSYS_H_
00013 
00014 #include "REAL.H"
00015 #include "RealVect.H"
00016 #include "FArrayBox.H"
00017 #include "FluxBox.H"
00018 #include "EdgeQuadrature.H"
00019 #include "NewCoordSys.H"
00020 
00021 #include "NamespaceHeader.H"
00022 
00023 /// Virtual base class encapsulating mapped-grid coordinate systems
00024 /**
00025    The NewFourthOrderCoordSys class is a pure-virtual base class
00026    an API for performing operations in a mapped grid coordinate space.
00027    Implementations of specific coordinate systems will be provided by
00028    derived classes.
00029 
00030    The following pure virtual functions of NewCoordSys are defined in
00031    NewFourthOrderCoordSys.cpp :
00032 
00033 
00034   virtual void getNodeRealCoordinates(FArrayBox& a_nodeCoords,
00035                                       const Box& a_box) const;
00036   // calls virtual RealVect realCoord(const RealVect& a_Xi) const = 0;
00037 
00038   virtual int getNumN() const
00039   {
00040     return SpaceDim*SpaceDim;
00041   }
00042 
00043   virtual int getNcomponent(const int a_row, const int a_col) const
00044   {return (SpaceDim*a_col) + a_row;}
00045 
00046   virtual void cellVol(FArrayBox& a_vol,
00047                        const FluxBox& a_N,
00048                        const Box& a_box) const;
00049   // calls virtual RealVect realCoord(const RealVect& a_Xi) const = 0
00050   // on corners, then calls fourthOrderAverageFace and computeDivergence.
00051 
00052   virtual void getN(FluxBox& a_N, const Box& a_box) const;
00053 
00054   virtual void getAvgJ(FArrayBox& a_avgJ, const Box& a_box) const;
00055   // calls cellVol.
00056 
00057   virtual void getAvgJinverse(FluxBox& a_avgJinverse,
00058                               const Box& a_box) const;
00059 
00060   virtual void pointwiseJ(FArrayBox& a_J,
00061                           const FArrayBox& a_Xi,
00062                           const Box& a_box) const;
00063   // calls dXdXi.
00064 
00065   virtual void computeMetricTermProductAverage(FluxBox& a_product,
00066                                                const FluxBox& a_F,
00067                                                const FluxBox& a_N,
00068                                                const FluxBox& a_FforGrad,
00069                                                const Box& a_box,
00070                                                bool a_fourthOrder=true) const;
00071 
00072   virtual const RealVect& dx() const
00073   {
00074     return m_dx;
00075   }
00076 */
00077 class NewFourthOrderCoordSys : public NewCoordSys
00078 {
00079 public:
00080   /// default constructor
00081   NewFourthOrderCoordSys();
00082 
00083   /**
00084      Destructor.
00085    */
00086   virtual ~NewFourthOrderCoordSys();
00087 
00088   /// given coordinate in mapped space, return its location in real space
00089   virtual RealVect realCoord(const RealVect& a_Xi) const = 0;
00090 
00091   /// given coordinate in real space, return its location in the mapped space
00092   virtual RealVect mappedCoord(const RealVect& a_x) const = 0;
00093 
00094   /// given coordinates in mapped space, return locations in real space
00095   virtual void realCoord(FArrayBox& a_x, const FArrayBox& a_Xi,
00096                          const Box& a_box) const;
00097 
00098   /// given coordinate in real space, return its location in the mapped space
00099   virtual void mappedCoord(FArrayBox& a_Xi, const FArrayBox& a_x,
00100                            const Box& a_box) const;
00101 
00102   /// return Cartesian XYZ locations of nodes
00103   /** nodeCoords should have dimension returned by dimension()
00104    */
00105   virtual void getNodeRealCoordinates(FArrayBox& a_nodeCoords,
00106                                       const Box& a_box) const;
00107 
00108   /// note that a_Xi is in mapped space.
00109   virtual Real  dXdXi(const RealVect& a_Xi,
00110                       int a_dirX,
00111                       int a_dirXi) const = 0;
00112 
00113   /// note that a_X is in mapped space.
00114   /** fills the destComp component of a_dxdXi
00115       with the derivative of x w/ respect to Xi
00116       in the dirX direction
00117       (default implementation is there in case derived class doesn't
00118       provide it)
00119   */
00120   virtual void dXdXi(FArrayBox& a_dxdXi,
00121                      const FArrayBox& a_Xi,
00122                      int a_destComp,
00123                      int a_dirX,
00124                      int a_dirXi,
00125                      const Box& a_box) const;
00126 
00127 
00128   /// returns number of components in the metric term matrix (N)
00129   virtual int getNumN() const
00130   {
00131     return SpaceDim*SpaceDim;
00132   }
00133 
00134   /// index function into face-centered metric terms
00135   /** returns which component of the face-centered metric term in
00136       which N^row_col is stored.
00137       \note
00138       <ul>
00139         <li> \f$N\f$ is stored in column-major order (fortran).
00140         <li> Then using the same format, \f$N^T\f$ is stored in
00141              row-major order (C).
00142         <li> so (0,x) gives the start of a column of \f$N\f$ and
00143              a row of \f$N^T\f$.
00144         <li> If you want the computational flux in direction 'd',
00145              from the physical flux, 'f', take the scalar product
00146              of N(getNcomponent(0, d)) dot f
00147       </ul>
00148   */
00149   virtual int getNcomponent(const int a_row, const int a_col) const
00150   {return (SpaceDim*a_col) + a_row;}
00151 
00152 
00153   /// computes the volume flux on the faces
00154   virtual void volFlux(FluxBox& a_volFlux,
00155                        const FluxBox& a_Nt,
00156                        const Box& a_box) const;
00157 
00158   /// computes cell volumes
00159   /**
00160      Return cell volume in a_vol on cells of a_box,
00161      using a_N from getN() defined on a_box grown by 1.
00162    */
00163   virtual void cellVol(FArrayBox& a_vol,
00164                        const FluxBox& a_N,
00165                        const Box& a_box) const;
00166 
00167   /// computes integral of N over each face of a_box
00168   virtual void getN(FluxBox& a_N, const Box& a_box) const;
00169 
00170   /// computes cell-averaged J
00171   virtual void getAvgJ(FArrayBox& a_avgJ,
00172                        const FluxBox& a_volFlux,
00173                        const Box& a_box) const;
00174 
00175   /// computes cell-averaged J
00176   virtual void getAvgJ(FArrayBox& a_avgJ, const Box& a_box) const;
00177 
00178   /// computes cell-averaged 1/J
00179   virtual void getAvgJinverse(FluxBox& a_avgJinverse,
00180                               const Box& a_box) const;
00181 
00182   /// computes cell-averaged 1/J, using pre-computed cell-avg J
00183   virtual void getAvgJinverse(FluxBox& a_avgJinverse,
00184                               const FArrayBox& a_avgJ,
00185                               const Box& a_box) const;
00186 
00187   /// Jacobian evaluated at location Xi in mapped space
00188   virtual Real pointwiseJ(const RealVect& a_Xi) const;
00189 
00190   /// Jacobian evaluated at locations Xi in mapped space
00191   virtual void pointwiseJ(FArrayBox& a_J,
00192                           const FArrayBox& a_Xi,
00193                           const Box& a_box) const;
00194 
00195   /// Computes 4th-order average of product = N^T*F.  F can be a flux dyad
00196   virtual void computeMetricTermProductAverage(
00197     FluxBox&                   a_product,
00198     const FluxBox&             a_F,
00199     const FluxBox&             a_N,
00200     const int                  a_NNumComp,
00201     const FluxBox&             a_FforGrad,
00202     const Box&                 a_box,
00203     bool                       a_fourthOrder = true,
00204     Interval                   a_varIntervalProduct = Interval(),
00205     Interval                   a_varIntervalF = Interval(),
00206     int                        a_fluxSpaceContiguous = 1,
00207     const ProblemDomain *const a_problemDomainPtr = NULL) const;
00208 
00209   /// computes 4th-order average of product = N^T*F.  F is vector of a scalar
00210   /** if a_fourthOrder is false, then only do a second-order dot product
00211       aFforGrad is the F used to compute grad_perp(F);
00212       returns results in a_product on a_box;
00213       both a_FforGrad and a_N (from getN()) must be defined on a_box grown by 1.
00214    */
00215   virtual void computeMetricTermProductAverage(FluxBox& a_product,
00216                                                const FluxBox& a_F,
00217                                                const FluxBox& a_N,
00218                                                const FluxBox& a_FforGrad,
00219                                                const Box& a_box,
00220                                                bool a_fourthOrder=true) const;
00221 
00222 
00223   /// computes 4th-order average of product = N^T*F.  F is vector of a scalar
00224   /** if a_fourthOrder is false, then only do a second-order dot product;
00225       returns results in a_product on a_box;
00226       both a_F and a_N (from getN()) must be defined on a_box grown by 1.
00227    */
00228   virtual void computeMetricTermProductAverage(FluxBox& a_product,
00229                                                const FluxBox& a_F,
00230                                                const FluxBox& a_N,
00231                                                const Box& a_box,
00232                                                bool a_fourthOrder=true) const;
00233 
00234   /// Computes magnitude of <N^T> on each face
00235   //  This is not virtual because I don't think it should even be in this class
00236   //  -- see comments in NewFourthOrderCoordSys.cpp.
00237   void magnitudeN(FluxBox&       a_NMag,
00238                   const FluxBox& a_N,
00239                   const Box&     a_box) const;
00240 
00241   /// given N and Jinverse, computes N/J
00242   /** provided for convenience, since the multicomponent N/J product
00243       requires a bit of processing to work with the existing
00244       fourth-order multiplication function
00245   */
00246   virtual void computeNJinverse(FluxBox& a_NJinverse,
00247                                 const FluxBox& a_Jinverse,
00248                                 const FluxBox& a_N,
00249                                 const Box& a_cellBox) const;
00250 
00251 
00252   /// access function to simplify things -- returns mapped-space cell spacing
00253   virtual const RealVect& dx() const
00254   {
00255     return m_dx;
00256   }
00257 
00258   /// this evaluates the script N values from equation 12 in Phil's notes
00259   /** note that a_Xi is in mapped space.
00260    */
00261   virtual Real getN(const RealVect& a_Xi, int a_s, int a_d, int a_d1) const;
00262 
00263   /// note that a_Xi is in mapped space.
00264   virtual Real getNMatrixEntry(const RealVect& a_Xi,
00265                                int a_s, int a_d, int a_d1,
00266                                int a_row, int a_column) const;
00267 
00268   ///  Integrates \f$\mathcal{N}^s\f$ on a specific set of codimension 2
00269   ///  hyperedges that are orthogonal to the given directions.
00270   virtual void integrateScriptN(FArrayBox& a_scrN,
00271                                 const int a_dir0,
00272                                 const int a_dir1,
00273                                 const Box& a_box) const;
00274 
00275   virtual void incrementFaceMetricWithEdgeTerm(FArrayBox& a_faceMetrics,
00276                                                int a_faceDir,
00277                                                int a_edgeDir,
00278                                                const Box& a_box,
00279                                                const Side::LoHiSide& a_side) const;
00280 
00281   virtual void computeTransverseFaceMetric(FArrayBox& a_faceMetrics,
00282                                            const Box& a_box,
00283                                            int a_faceDir,
00284                                            int a_dDir) const;
00285 
00286 
00287   virtual void computeTangentialGrad(FluxBox& a_gradPhi,
00288                                      const FluxBox& a_phiFace,
00289                                      const Box& a_box) const;
00290 
00291 
00292   virtual int tanGradComp(const int a_faceDir, const int a_tanDir,
00293                           const int a_comp) const;
00294 
00295   /// return row a_idir of contravariant metric
00296   virtual void contravariantMetric(FArrayBox& a_metric,
00297                                    int a_dir) const;
00298 
00299   /// convert vector components (in a_vectorIntv) of a_flux to orthonormal frame to that of this block
00300   virtual void orthonormalize(FluxBox& a_flux,
00301                               const Interval& a_vectorIntv) const;
00302 
00303   /// overwrite vector components (in a_vectorIntv) of a_fluxFab with orthonormalized
00304   virtual void orthonormalize(FArrayBox& a_fluxFab,
00305                               const Interval& a_vectorIntv) const;
00306 
00307   /// convert components in a_csComps of vector a_csFab in coordinate-system basis to components in a_orthoComps of vector a_orthoFab in orthonormal basis, at points with indices in a_box, on lines where direction a_idir is constant
00308   virtual void orthonormalize(
00309                               const FArrayBox& a_csFab,
00310                               FArrayBox& a_orthoFab,
00311                               const Box& a_box,
00312                               int a_idir,
00313                               const IntVect& a_csComps,
00314                               const IntVect& a_orthoComps) const;
00315 
00316   /// convert vector components (in a_vectorIntv) of a_flux from orthonormal frame to that of this block
00317   virtual void deorthonormalize(FluxBox& a_flux,
00318                                 const Interval& a_vectorIntv) const;
00319 
00320   /// overwrite vector components (in a_vectorIntv) of a_fluxFab with deorthonormalized
00321   virtual void deorthonormalize(FArrayBox& a_fluxFab,
00322                                 const Interval& a_vectorIntv) const;
00323 
00324   /// convert components in a_orthoComps of vector a_orthoFab in orthonormal basis to components in a_csComps of vector a_csFab in coordinate-system basis, at points with indices in a_box, on lines where direction a_idir is constant
00325   virtual void deorthonormalize(
00326                                 const FArrayBox& a_orthoFab,
00327                                 FArrayBox& a_csFab,
00328                                 const Box& a_box,
00329                                 int a_idir,
00330                                 const IntVect& a_orthoComps,
00331                                 const IntVect& a_csComps) const;
00332 protected:
00333 
00334   RealVect m_dx;
00335 
00336   // which component directions to use when computing cell volumes
00337   Interval m_volInterval;
00338 
00339   EdgeQuadrature* m_quadraturePtr;
00340 };
00341 
00342 
00343 
00344 #include "NamespaceFooter.H"
00345 
00346 #endif