CutCellMoments.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef _CUTCELLMOMENTS_H_
#define _CUTCELLMOMENTS_H_

#if defined(CH_Darwin) && defined(__GNUC__) && ( __GNUC__ == 3 )
// deal with the broken isnan()/isinf() in GCC on MacOS
#include <unistd.h>
#define _GLIBCPP_USE_C99 1
#endif

#include <map>
using std::map;

#include "Vector.H"
#include "REAL.H"
#include "IndexTM.H"

#include "Notation.H"
#include "LSquares.H"
#include "IFData.H"
#include "RefinementCriterion.H"

#include "NamespaceHeader.H"

template <int dim> class LSProblem;

template <int dim> class CutCellMoments
{
public:
  typedef IndexTM<int,dim>  IvDim;
  typedef IndexTM<Real,dim> RvDim;

  typedef map<IvDim,Real,LexLT <IvDim > >  PthMoment;

  typedef map<IndexTM<int,dim-1>,Real,LexLT < IndexTM<int,dim-1> > > PthMomentLesserDimension;

  typedef map<IndexTM<int,1>,Real >     OneDMoments;
  typedef map<int,IvDim>                LocPthMoment;
  typedef map<IvDim,int,LexLT <IvDim> > PthMomentLoc;

  typedef map<Iv2,CutCellMoments<dim-1>, LexLT<Iv2> > BdCutCellMoments;

  // Constructors
  inline CutCellMoments();
  inline CutCellMoments(const CutCellMoments<dim>& a_cutCellMoments);

  // This is used to build bd CutCellMoments
  inline CutCellMoments(const IFData<dim>& a_info);

  // Destructor
  inline ~CutCellMoments();

  // Returns a boundary element of one dim less
  inline const CutCellMoments<dim - 1> getBdCutCellMoments(const Iv2& a_bdId) const;

  // Solve
  inline void computeMoments(const int           & a_order,
                             const int           & a_degreeP,
                             const bool          & a_useConstraints,
                             RefinementCriterion & a_refinementCriterion);

  // Integrates a monomial over a full cell
  inline Real fullCellQuadrature(const IvDim& a_mono);

  inline Vector<Real> computeRhs(LSProblem<dim> & a_lsp,
                                 const int      & a_order);

  // Solve the problem on a refined cell
  inline Vector< CutCellMoments<dim> > refine(const int           & a_order,
                                              const int           & a_degreeP,
                                              const bool          & a_useConstraints,
                                              RefinementCriterion & a_refinementCriterion,
                                              const Vector<int>   & a_refineInDir);

  inline void addMomentMaps(const Vector<CutCellMoments<dim> > & a_refinedCutCellVector,
                            const int                          & a_degreeP,
                            const bool                         & a_useConstraints);

  inline void addMoments(PthMoment               & a_momentMap,
                         PthMoment               & a_refinedMomentMap,
                         const IndexTM<Real,dim> & a_refinedCenterDelta);

  inline void addBdMoments(CutCellMoments<dim>     & a_coarseCutCell,
                           const IFData<dim+1>     & a_iFData,
                           const int               & a_degreeP,
                           const bool              & a_useConstraints,
                           const IndexTM<Real,dim> & a_refinedCenterDelta,
                           const IndexTM<int,dim>  & a_localHilo);

  inline Real changeMomentCoordinates(PthMoment               & a_refinedMomentMap,
                                      const IndexTM<int,dim>  & a_monomial,
                                      const IndexTM<Real,dim> & a_refinedCenterDelta);

  inline void initialize(CutCellMoments<dim> & a_refinedCutCell);

  inline void initializeMap(PthMoment & a_map1,
                     PthMoment & a_map2);

  inline void initializeMap(PthMomentLesserDimension & a_map1,
                     PthMomentLesserDimension & a_map2);

  inline void initializeMap(OneDMoments & a_map1,
                     OneDMoments & a_map2);

  // Get the value of the moment depending on the values of the booleans
  // m_allVerticesIn/Out these functions are used under refinement, when
  // moments on the faces are needed even on a covered/regular cell
  inline Real getBdMoment(const IvDim             & a_mono,
                          const IFData<dim+1>     & a_iFData,
                          const IndexTM<Real,dim> & a_refinedCenterDelta,
                          PthMoment                 a_fullCellMap = PthMoment());

  inline Real getBdEBMoment(const IvDim             & a_mono,
                            const IFData<dim+1>     & a_iFData,
                            const IndexTM<Real,dim> & a_refinedCenterDelta);

  inline void computeResiduals(const int  & a_order,
                               const int  & a_degreeP,
                               const bool & a_useConstraints);

  inline void computeResiduals(const Vector<CutCellMoments<dim> > & a_refinedCCMoms,
                               const int & a_degreeP);

  // Output methods that check a_mono is in the map
  inline Real getMoment(const IvDim   & a_mono,
                        const EBorVol & a_EBorVOL) const;

  // Methods for reading geom data that do sanity checks on results
  inline Real getVol(const EBorVol & a_EBorVol) const;

  inline RvDim getCentroid(const EBorVol & a_EBorVOL) const;

  inline Real getResidual(const int & a_iDegree,
                          const int & a_normJ) const;

  inline bool isCovered() const;

  inline bool isRegular() const;

  // Called via: pout() foo
  inline void print(ostream& out) const;

  inline void dump() const;

  inline void printMoments();

  inline void printBdMoments();

  // Operators
  inline void operator=(const CutCellMoments<dim>& a_cutCellMoments);

  // Member data
  PthMoment         m_moments;
  PthMoment         m_EBmoments;
  BdCutCellMoments  m_bdCutCellMoments;
  IFData<dim>       m_iFData;
  bool              m_bdCCOn;

  // Least Square residual
  Vector<Vector<Real> > m_residual;
};

// One dimensional cutCellMoments
template <> class CutCellMoments<1>
{
public:
  typedef map<IndexTM<int,1>,Real> OneDMoments;

  // Constructors
  inline CutCellMoments();
  inline CutCellMoments(const CutCellMoments<1> & a_cutCellMoments);
  inline CutCellMoments(const IFData<1>& a_info);

  // Destructor
  inline ~CutCellMoments();

  inline void computeMoments(const int           & a_order,
                             const int           & a_degree,
                             const bool          & a_useConstraints,
                             RefinementCriterion & a_refinementCriterion);

  inline void simpleComputeMoments(const Real & a_loPt,
                                   const Real & a_hiPt,
                                   const int  & a_degree);

  inline void computeMomentsUsingBinomial(const Real & a_loPt,
                                          const Real & a_hiPt,
                                          const int  & a_loSign,
                                          const int  & a_hiSign,
                                          const int  & a_degree);

  inline void addBdMoments(CutCellMoments<1>     & a_coarseCutCell,
                           const IFData<2>       & a_iFData,
                           const int             & a_degreeP,
                           const bool            & a_useConstraints,
                           const IndexTM<Real,1> & a_refinedCenterDelta,
                           const IndexTM<int,1>  & a_localHilo);

  inline Real changeMomentCoordinates(OneDMoments           & a_refinedMap,
                                      const IndexTM<int,1>  & a_monomial,
                                      const IndexTM<Real,1> & a_refinedCenterDelta);

  inline void initialize(CutCellMoments<1> & a_refinedCutCell);

  inline void initializeMap(OneDMoments & a_map1,
                            OneDMoments & a_map2);

  inline Real getBdMoment(const IndexTM<int,1>  & a_mono,
                          const IFData<2>       & a_iFData,
                          const IndexTM<Real,1> & a_refinedCenterDelta,
                          OneDMoments             a_fullCellMap = OneDMoments());

  inline Real getBdEBMoment(const IndexTM<int,1>  & a_mono,
                            const IFData<2>       & a_iFData,
                            const IndexTM<Real,1> & a_refinedCenterDelta);

  // Output method that check a_mono is in the map
  inline Real getMoment(const IndexTM<int,1> & a_mono,
                        const EBorVol        & a_EBorVOL) const;

  inline Real getMoment(const IndexTM<int,1>& a_mono) const;

  // Methods for reading geom data that do sanity checks on results
  inline Real getVol(const EBorVol& a_EBorVol) const;

  inline IndexTM<Real,1> getCentroid(const EBorVol& a_EBorVol) const;

  inline bool isCovered() const;

  inline bool isRegular() const;

  // Output
  inline void print(ostream& out) const;

  inline void dump() const;

  // Operators
  inline void operator=(const CutCellMoments<1>& a_cutCellMoments);

  // Member data
  OneDMoments  m_moments;
  IFData<1>    m_iFData;
  bool         m_bdCCOn;

  // This is a dummy quantity: no EB in 1D
  OneDMoments  m_EBmoments;
};

#include "NamespaceFooter.H"

#include "CutCellMomentsImplem.H"

#endif

---