moments.hh

Go to the documentation of this file.

 
// TODO Some more documentation what this class is about, <emph>i.e.</emph> <grad_n u_h> as Mean value or for neumannboundary etc
 
#ifndef moments_A0_hh
#define moments_A0_hh
 
#include "geometry.hh"
#include "basics/typedefs.hh"
 
#include "hp2D/space.hh"
 
#include "patches.hh"
//#include "linInnerProd.hh"
#include "innerResidual.hh"
 
//needed spaces
#include "hp2D/neumannTraceSpace.hh"
#include "hp1D.hh"
 
#include "operator.hh"
 
#include "toolbox/hashMap.hh"
#include "formula/exceptions.hh"
#include "basics/exceptions.hh"
 
#include "geometry/cell.hh"
#include "space/element.hh"
#include "space/formula.hh"
 
#include "formula/parsedFormula.hh"
 
#include <limits>
#  define EPS std::numeric_limits<double>::epsilon()
 
namespace hp2D {
 
template<typename F = concepts::Real>
class ApproxMoments: public concepts::OutputOperator{
  //need of access to underlying Element Info
  template<class>
  friend class EquilibratedMomentsAO;
 
public:
  typedef concepts::ElementAndFacette<hp2D::Element<F> > UnderlyingElement;
 
 
  ApproxMoments(const concepts::SpaceOnCells<F>& spc,const concepts::Vector<F>& sol);
 
    void addNeumann(concepts::ParsedFormula<F>& frm,
        concepts::Set<uint>& set) {
      frms_.push_back(frm);
      nAttrbs_.push_back(set);
    }
 
    void addNeumann(concepts::ParsedFormula<F>& frm, uint attrib) {
      frms_.push_back(frm);
      nAttrbs_.push_back(concepts::Set<uint>(attrib));
    }
 
    void addDirichlet(concepts::Set<uint>& set) {
      dAttrbs_ = dAttrbs_||set;
    }
 
    void addDirichlet(uint attrb) {
      dAttrbs_.insert(attrb);
    }
 
 
    void addComplete();
 
 
    //Access routine to write and read moments.
    const concepts::ElementMatrix<F>& operator()(uint edgeKey) const {
      typename concepts::HashMap<concepts::ElementMatrix<F> >::const_iterator iter = hashM_.find(edgeKey);
      conceptsAssert(iter!= hashM_.end(),concepts::Assertion());
      return iter->second;
    }
 
    //@pre the edgeKey exists in the underlying list
    //@pre the InputQuad is a underlying quad of Edge edgekey
    //sets a weight if input Quad is first underlying quad of edge E to w = 1  if it is the second w = -1, else
    Real getWeight(uint quadKey,uint edgeKey) const{
      typename concepts::HashMap<concepts::Sequence<UnderlyingElement> >::const_iterator iter = uelm_.find(edgeKey);
      conceptsAssert(iter!=uelm_.end(),concepts::Assertion());
      if (iter->second[0].elm->support().key().key()==quadKey)
        return 1.0;
      return -1.0;
 
//      conceptsAssert(hashM_.exists(edgeKey),concepts::Assertion());
//      //the first (or only) uelm is the inputquad
//      if(uelm_[edgeKey][0].elm->support().key().key()==quadKey)
//          return 1.0;
//      //the second due to pre have to be the InputQuad
//      return -1.0;
 
    }
 
protected:
 
  virtual std::ostream& info(std::ostream& os) const;
 
 
private:
 
 
 
  //local reference of the space on which the moments depend
  const concepts::SpaceOnCells<Real>& spc_;
 
  //map representing the moments via  hashM[E] = mean(moments) on Edge E corresponding to the first underlying Element (e.g. Quad) K
  concepts::HashMap<concepts::ElementMatrix<F> > hashM_;
 
 
  const concepts::Vector<F>& sol_;
 
  //HashMap from Edge Key to its underlying elements
  // as info is given in TraceSpace and NeumannTracespaces, its just pointed here to it
  //this info might be the same in the Patch map ittself
  concepts::HashMap<concepts::Sequence<UnderlyingElement> > uelm_;
 
 
  void computeID_();
  void computeN_();
 
  //ushort m_;
 
  //Neumann formulas
  concepts::Sequence<concepts::ParsedFormula<F> > frms_;
  //Neumann Attributes
  concepts::Sequence<concepts::Set<uint> > nAttrbs_;
  //Dirichlet Attributes
  concepts::Set<uint> dAttrbs_;
 
};//class
 
 
 
 
//@pre continuous Space, since of resize init routine, ie.
 
template<typename F = concepts::Real>
class EquilibratedMomentsAO: public concepts::OutputOperator{
 
public :
 
  typedef concepts::ElementAndFacette<hp2D::Element<F> > UnderlyingElement;
 
  EquilibratedMomentsAO(const concepts::SpaceOnCells<F>& spc,
            const geometry::VtxToPatchMaps& patchMap,
            const concepts::InnerResidual<F>& innerRes,
            const hp2D::ApproxMoments<F>& apprxMoments);
 
  virtual ~EquilibratedMomentsAO();
 
  const concepts::Vector<F>& operator()(uint edgeKey) const {
    typename concepts::HashMap<concepts::Vector<F> >::const_iterator iter = hashM_.find(edgeKey);
    conceptsAssert(iter!= hashM_.end(),concepts::Assertion());
    return iter->second;
  }
 
 
protected:
  virtual std::ostream& info(std::ostream& os) const;
 
private:
 
  //local reference of the space on which the moments depend
  const concepts::SpaceOnCells<Real>& spc_;
 
  //map representing the moments via  hashM[E] = moments on Edge E corresponding to first UnderlyingElement (e.g. Quad) K
  //to get moments for the second underlying element K', moments just have to be multiplicated with (-1)
  concepts::HashMap<concepts::Vector<F> > hashM_;
 
 
 
  const concepts::InnerResidual<F>& innerRes_;
  const hp2D::ApproxMoments<F>& apprxMoments_;
 
  void initMoments_(const geometry::VtxToPatchMaps& patchMap);
 
  //add routine for higher Moments, @pre Edge E is the k-th Edge of Quad K with Quad->p()={px,py}, and K is the first uelm of E.
  void addHigherMoments_(const uint K, const uint E, const uint k,const uint px, const uint py);
 
  void buildRhs_(const geometry::VtxToPatchMaps& patchMap,
         const concepts::InnerResidual<F>& innerRes, const  hp2D::ApproxMoments<F>& apprxMoments, concepts::HashMap<concepts::Vector<Real> >& rhs);
 
  //Routine that solves the first order considered Patchproblems
  //vec respresenting the right hand sides as input and representing the solutions as output
  void solve_(const geometry::VtxToPatchMaps& patchMap, concepts::HashMap<concepts::Vector<Real> >& vec);
 
  //routine to precompute all needed inverse matrices due to patchsizes and appearing patchtypes
  void precomputeInverse_(const geometry::VtxToPatchMaps& patchMap, concepts::HashMap<concepts::Vector<Real> >& rhs);
 
  void computeMoments_(const geometry::VtxToPatchMaps& patchMap,
      const concepts::InnerResidual<F>& innerRes,
      const hp2D::ApproxMoments<F>& apprxMoments,
      concepts::HashMap<concepts::Vector<Real> >& sigma);
 
  //Method to solve a inner Patch problem with
  //@ param patchsize size of the Patch
  //@ param patchRhs  given Rhs of the Patch computed with InnerResidual and apprxMoments
  // the patchRhs holds the Solution later, coz of reference the rhs can not be used anymore, indeed this
  // is no necessary anyway in the algorithm
  void solveInner_(concepts::Vector<Real>& patchRhs);
 
  //solves inner Patch problem for given Patchsize and given rhs
  // uses precomputed T+D, since T is singular  the rhs gets corrected
  // @ param N  size of the patch
  // @ param patchRhs right handed side of the Patchproblem
  void solveNN_(concepts::Vector<Real>& patchRhs);
 
 
  //solves a Dirichlet_Neumann patch problem for given patchsize and given rhs
  // uses precomputed inverse of T , since T is regular on those problems
  // @ param N  size of the patch
  // @ param patchRhs right handed side of the Patchproblem
  void solveDN_(concepts::Vector<Real>& patchRhs);
 
 
  //solves a Neumann-Dirichlet patch problem for given patchsize and given rhs
  // uses precomputed inverse of the T_{DN} matrices due to symmetrie to the DN problem.
  void solveND_(concepts::Vector<Real>& patchRhs);
 
  //solves a local Dirichlet-Dirichlet patch problem for given patchsize and given rhs-
  // uses precomputed inverses of the T_{DD} matrices in ddInvM_
  void solveDD_(concepts::Vector<Real>& patchRhs);
 
  //Precomputes the inner Matrix  inv(0.5*( T+ones() ) for requested sizes
  void precomputeInnerInverse_(const concepts::Set<uint>& dSet);
 
  //precomputes the NN- Matrix inv(0.5*(T+ones() ) )
  void precomputeNNInverse_(const concepts::Set<uint>& dSet);
 
  //precomputes the DN - Matrix inverses  inv (0.5*T) since here T is invertable
  //@ param nr = number of precomputed inverse Matrices
  void precomputeDNInverse_(const concepts::Set<uint>& dSet);
 
  //precomputed DD - matrix inverse for some request of demanded patchsizes
  // computes : (1/2*T_{DD})^{-1}
  void precomputeDDInverse_(const concepts::Set<uint>& dSet);
 
  //InnerPatch Circulant Matrix inverses :
  concepts::HashMap<concepts::Array<Real> > innerCircM_;
 
  //Neumann_Neumann Patch Inverse Matrix (T+D)^(-1) computed up to ..
  concepts::HashMap<concepts::DenseMatrix<Real> > nnInvM_;
 
  //Dirichlet_Neumann Patch Inverse Matrix  T^(-1) computed up to wanted number
  concepts::HashMap<concepts::DenseMatrix<Real> > dnInvM_;
 
  //Dirichlet_Dirichlecht Patch Inverse Matrix T^(-1) computed up to given dim
  concepts::HashMap<concepts::DenseMatrix<Real> > ddInvM_;
 
  //local reference to the underlying Elements information in the approximated Moments
  const concepts::HashMap<concepts::Sequence<UnderlyingElement> >& uelm_;
 
};
 
 
 
 
} //namespace
 
#endif // moments_A0_hh