innerResidual.hh

Go to the documentation of this file.

 
#ifndef innerResidual_own_hh
#define innerResidual_own_hh
 
//TODO: Includes
#include "basics/exceptions.hh"
 
#include "space.hh"
#include "hp2D.hh"
#include "geometry.hh"
#include "toolbox.hh"
#include "operator.hh"
 
#include "basics.hh"
#include "function.hh"
 
#define  fillEntrieAddM 0
#define  ADD_D 0
 
namespace concepts{
 
// forward declaration
template<class F, class G>
class LinearForm;
 
template<class F>
class Space;
 
using concepts::Real;
 
/* TODO:
 * - Documenntation beispielklassen nennen die nur auf randbasisfunktionen aufbauen
 * - Schreibe Operator und PatchRHS operator
 * - dont make ElementMatrix use Array
 *
 * -Der Space muss ein hpAdaptiveSpaceH1 erstmal sein sonst nicht implementiert
 * -Die Linearform muss Elemente aus dem entsprechenden Raum nur annehmen können!
 * -Die Linearform darf nur eine Linearform_0 oder _1 sein da diese nur die enstsprechenden
 * -Shapefunktion integrale berechnet ! Alles andere wäre falsch
 * -Vermeide doppelte Matrix resizing durch fillEntries bei add und Constructor
 *
 * -Deconstructor schreiben
 * -Add routine muss linearForm bekommen die auf den Space zugreift
 *
 * TODO: Const routine umschreibung, todo: documentation
 */
 
template<class F>
class InnerResidual: public OutputOperator { //Function<F> {
public:
 
  template<class G>
  InnerResidual(const SpaceOnCells<G>& spc);
 
 
  template<class G>
  InnerResidual(const SpaceOnCells<G>& spc,LinearForm<F, G>& lform, const F a = 1);
 
 
  virtual ~InnerResidual() {};
 
  inline const concepts::ElementMatrix<F>& operator[](uint quadKey){
    conceptsAssert(hashM_.exists(quadKey),Assertion());
    return hashM_[quadKey];
  }
 
 
  inline const concepts::ElementMatrix<F>& operator()(uint quadKey) const{
    typename concepts::HashMap<concepts::ElementMatrix<F> >::const_iterator iter = hashM_.find(quadKey);
    assert(iter!=hashM_.end());
    return iter->second;
  }
 
  const F operator()(const hp2D::Quad<F>& quad,const uint locNode) {
    uint K = quad.support().key();
    conceptsAssert(hashM_.exists(K),Assertion());
    return hashM_[K].getData()[getPos_(locNode,quad.p())];
  }
 
 
  template<class G>
  void add(LinearForm<F, G>& lform, const F a = 1);
 
 
 
  /*Operator gets a Patch an computes right hand side for the
   * corresponding linear system due to the topology of the
   * Patch
   */
 
  /* template<class G>
   Vector<F> operator()(ElementPatch<G> Patch){
   }*/
 
 
 
protected:
  inline std::ostream& info(std::ostream& os) const {
    return os << concepts::typeOf(*this) <<"( \n" << hashM_ << "\n )";
  }
 
 
 
 
private:
 
  //TODO: dont make ElementMatrix use Array
  concepts::HashMap<concepts::ElementMatrix<F> > hashM_;
  //local reference to the current space
  const Space<F>& spc_;
 
 
  //maps from ElementKey to its belonging polynomial degree.
  concepts::HashMap<const ushort*> degP_;
 
  //all appearing polynomial degrees in x and y direction
  //concepts::Sequence<const ushort*> degP_;
 
 
 
  template<class G>
  void fillEntries_(LinearForm<F, G>& lform, const F a = 1);
 
 
 
//    Method is used to read the correct innerResidual due to first order Basisfunctions.
//    Let N0,N1, E1,...,EP-1 be the Basis functions on [-1,1] ( [0,1] )
//    Then Innerresidual is computed due to first x then y tensor basis ordering :
//    Elementmatrix of size 2(p[0]+p[1]) iterative filled with :
//
//    N0_x*N0_y , N1_x * N0_y , E1_x *N0_y,... EP-1_x * N0_y    first p[0]+1 entries
//    increasing y-directed basis:
//    N0_x*N1_y , N1_x * N1_y , E1_x *N1_y,... EP-1_x * N1_y       ...
//    increasing y-directed basis:
//    ...
//
//    So the necassary nodal basis functions are
//    N0_x*N0_y       at positon 0
//    N1_x * N0_y    at position 1
//    N0_x*N1_y       at position p[0]+1
//    N1_x * N1_y     at position p[0]+2
//
//    where p[0] is degree of basis functions in x direction
//
//
//    i=0..3  is the i-th (local numbered) node of quad, this local numbering is counter clockwise.
//
//    But the polynomials are given first x then y order,
//    means for first order basis functioN:
//
//         0 node is 0th basis function
//          1 node is 1th basis function
//          2 node is 3th basis function
//          3 node is 2th basis function
//
//    For higher order including edge functions therefore  0-> 0 , 1 -> 1 , 2->p[0]+1+1 , 3->p[0]+1
//    For p[0]=1 this coincides with the shown example.
 
  //TODO: This will be be used in new class anymore
  uint getPos_(const uint i, const ushort* p) {
    conceptsAssert(i<4, concepts::Assertion());
    uint pos = 0;
    switch (i) {
    case (0):
      return pos = 0;
      break;
    case (1):
      return pos = 1;
      break;
    case (2):
      return pos = p[0] + 2;
      break;
    case (3):
      return pos = p[0] + 1;
      break;
    default:
      return -1;
    }
  }
 
};
 
}//namespace
#endif