LSProblem.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 _LSPROBLEM_H_
#define _LSPROBLEM_H_

#include <map>
using std::map;

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

#include "NamespaceHeader.H"

template <int dim> class CutCellMoments;

template <int dim> class LSProblem
{
  typedef IndexTM<int,dim>IvDim;
  typedef map<IvDim,int,LexLT <IvDim> > PthMomentLoc;
  typedef map<int,IvDim> LocPthMoment;

public:

  //destructor will free m_matrix and free matrixes in m_rhs
 ~LSProblem();

  //constructor for just accessing fill map function
  LSProblem(const int  & a_degreeP,
            const bool & a_useConstraints);

  //constructor
  LSProblem(const int                & a_orderAccuracy,
            const int                & a_degreeP,
            const bool               & a_useConstraints,
            const IndexTM<Real,dim>  & a_normal);

  //Solves over determined system matrix*unknown = rhs
  void invertNormalEq(const Vector<Real>& a_rhs,Vector<Real>& a_residual);

  //accessor functions
  int getDegree(){ return m_degreeP; }
  int getOrderAccuracy(){ return m_order;}

  void getMatrix(Real ** a_matrix);
  void getRhs(Vector<Real>& a_rhs);

  //getUnknowns will return the geometry info at the current level
  void getUnknowns(Vector<Real>& a_unknown);


  void computeBounds(const IndexTM<Real,dim> & a_dx,
                     const CutCellMoments<dim> & a_ccm);




  void print(ostream& out) const;

  void outputMatrix()const;
  void outputRhs()const;
  void outputUnknowns()const;

  const LocPthMoment& getLocMonomialMapDegreeP()const{return (const LocPthMoment&)m_locMonoP;}
  const PthMomentLoc& getMonomialLocMapDegreeP()const{return (const PthMomentLoc&)m_monoLocP;}
  const PthMomentLoc& getMonomialLocMapDegreePLess1()const{ return (const PthMomentLoc&)m_monoLocPLess1;}
  Real getUnknown(int loc){ return m_unknowns[loc]; }
  int getNumberDegP(){return m_numP;}
  int getNumberDegPLess1(){return m_numPLess1;}


private:
  void setRhs(const Vector<Real>& a_rhs); //deprecated

  //constructs matrix of overdetermined system
  void setMatrix();



  void momentBounds(Real & lobnd,
                    Real & hibnd,
                    const IvDim & mono,
                    const IndexTM<Real,dim>& a_dx);





  //\partial/\partial xidir (mono) = coeff*Dmono
  void differentiate(int        & a_coeff,
                     IvDim      & a_Dmono,
                     int        & a_idir,
                     const IvDim& a_mono);



  int nChooseR(int n, int r);

  //makes a list of monomials of degreeP
  void fillMap(PthMomentLoc& m_monoLoc,
               LocPthMoment& m_locMono,
               const int& a_degree);

  //computes number of monomials of degree a_monoDegree
  int numMonomials(const int& a_monoDegree);
  int factorial(const int& a_n,const int &a_m=0);

  //array tools
  void allocArray(const int& rows,
                  const int& cols,
                  Real**& A);

  void freeArray(const int& rows,
                 const int& cols,
                 Real**& A);


  //order of accuracy
  int m_order;

  //degree of monomial
  int m_degreeP;

  //use constrained LS
  bool m_useConstraints;

  //normal
  IndexTM<Real,dim>m_normal;

  //monomials of degree P
  PthMomentLoc m_monoLocP;
  LocPthMoment m_locMonoP;
  //number of monomials of degreeP
  int  m_numP;

  //monomials of degree P-1
  PthMomentLoc m_monoLocPLess1;
  LocPthMoment m_locMonoPLess1;
  //number of monomials of degreeP
  int  m_numPLess1;

  //matrix
  Real**  m_matrix;
  Vector<Real> m_unknowns;
  Vector<Real> m_rhs;



  Real *  m_constraintVec;
  Real ** m_constraintMat;
  int m_nInequality;

};

template<>class LSProblem<1>
{
public:

  //number of monomials of degree P
  int m_numP;

  //number of monomials of degree P-1
  int m_numPLess1;

  //monomials of degree P
  Real m_pMoments;

  //used in filling the map
  int m_ithMoment;



  //Destructor
  ~LSProblem();

  //copy constructor
  LSProblem(const LSProblem<1>& a_lsProblem);


  //empty Constructor
  LSProblem();

  //count num monomials
  int recursiveCount(const int& a_degreeP);

  void setNumMonomials();

  //print
  void print(ostream& a_out)const;

  //equals operator
  void operator = (const LSProblem & a_lSProblem);

};

#include "NamespaceFooter.H"

#include "LSProblemImplem.H"

#endif

---