vectorial.hh

Go to the documentation of this file.

 
#ifndef vectorialThings_hh
#define vectorialThings_hh
 
#include <memory>
#include "basics/cloneable.hh"
#include "toolbox/array.hh"
#include "toolbox/scannerConnectors.hh"
#include "toolbox/sharedPointer.hh"
#include "space/tmatrix.hh"
#include "space/element.hh"
#include "space/space.hh"
#include "operator/bilinearForm.hh"
#include "vectorial/graphics.hh"
#include "geometry/integral.hh"
 
// debugging
#include "basics/debug.hh"
 
#define VectorialConstr_D 0
#define VectorialPut_D 0
 
#define TMatrixOffsetConstr_D 0
 
#define TMatrixConstr_D 0
 
#define ElementConstr_D 0
 
#define BilinearFormConstr_D 0
 
#define LinearFormConstr_D   0
 
#define SpaceConstr_D 0
 
namespace concepts {
  // forward declaration
  template<typename F>
  class ElementGraphics;
 
  template<typename F, typename G>
  class ElementFunction;
}
 
namespace vectorial {
 
  using concepts::Real;
 
  // ************************************************************* Vectorial **
 
  template<class F>
  class Vectorial {
  public:
    Vectorial(uint vdim, uint arrayWidth) : 
      vdim_(vdim), idx_(0), vdata_(arrayWidth, 0) {
      DEBUGL(VectorialConstr_D,
             "vdim = " << vdim_ << ", arraywidth = " << arrayWidth);
    }
    virtual ~Vectorial();
    virtual void insert(F& vdata, const int a = 0, const int b = 0);
    virtual const F* get(const int a, const int b = 0) const
    { return vdata_[a]; }
    uint vdim() const { return vdim_; }
  protected:
    uint vdim_;
    uint idx_;
    mutable concepts::Array<F*> vdata_;
  private:
    Vectorial(const Vectorial<F>& c);
  }; 
 
  template<class F>
  Vectorial<F>::~Vectorial() {}
 
  template<class F>
  void Vectorial<F>::insert(F& vdata, int, int) {
    vdata_[idx_++] = &vdata;
    DEBUGL(VectorialConstr_D, "idx = " << idx_);
  }
 
  // ********************************************************* TMatrixOffset **
 
  template<class F>
  class TMatrix;
 
  template<class F>
  class TMatrixOffset : public concepts::TMatrixBase<F> {
    friend class TMatrix<F>;
  public:
    TMatrixOffset(const concepts::TMatrixBase<F> &T, int offsetRow,
                  int offsetColumn) :
      concepts::TMatrixBase<F>(0),
      offsetRow_(offsetRow), offsetColumn_(offsetColumn), T_(T) {
      DEBUGL(TMatrixOffsetConstr_D,
             "offsetRow = " << offsetRow_ << ", " <<
             "offsetColumn = " << offsetColumn_ << ", T = " << T_);
    }
    virtual ~TMatrixOffset();
    virtual void operator ()(const concepts::ElementMatrix<F>& A,
                             concepts::ElementMatrix<F>& B) const
    { return T_(A,B); }
    virtual void operator ()
      (const concepts::ElementMatrix<std::complex<F> >& A,
       concepts::ElementMatrix<std::complex<F> >& B) const
    { return T_(A,B); }
 
    virtual uint index(unsigned int i) const { return T_.index(i); }
 
    inline uint m() const { return T_.m(); }
 
    inline uint n() const { return T_.n(); }
 
    virtual void usedIdx(concepts::TColumn<bool>& c) const {}
    virtual void extract(const concepts::Vector<F>& solution,
                         concepts::Array<F>& coeff) const;
    virtual void extract(const concepts::Vector<std::complex<F> >& solution,
                         concepts::Array<std::complex<F> >& coeff) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    //private:
    uint offsetRow_, offsetColumn_;
    const concepts::TMatrixBase<F> &T_;
  };
 
  // *************************************************************** TMatrix **
 
  template<class F>
  class TMatrix : public concepts::TMatrixBase<F>,
                  public Vectorial<TMatrixOffset<F> > {
  public:
    TMatrix(uint vdim, uint arrayWidth = 0) :
      concepts::TMatrixBase<F>(0),
      Vectorial<TMatrixOffset<F> >(vdim,vdim),
      offsetRow_(0), offsetColumn_(0) {
      DEBUGL(TMatrixConstr_D, "done.");
    }
    virtual ~TMatrix();
 
    virtual void operator()(const concepts::ElementMatrix<F>& A, 
                            concepts::ElementMatrix<F>& B) const;
    virtual void operator()
      (const concepts::ElementMatrix<std::complex<F> >& A, 
       concepts::ElementMatrix<std::complex<F> >& B) const;
    void put(const concepts::TMatrixBase<F> &T, const int dim);
    
    virtual uint index(unsigned int i) const;
 
    // not needed in Vectorial problems
    virtual void usedIdx(concepts::TColumn<bool>& c) const {}
    virtual void extract(const concepts::Vector<F>& solution,
                         concepts::Array<F>& coeff) const;
    virtual void extract(const concepts::Vector<std::complex<F> >& solution,
                         concepts::Array<std::complex<F> >& coeff) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uint offsetRow_, offsetColumn_;
  };
 
 
  // ******************************************************* ElementNotBuilt **
 
  class ElementNotBuilt : public concepts::ExceptionBase {};
 
 
  // *********************************************************** ElementBase **
 
  // forward declaration
  template<class F>
  class Graphics;
 
  template<class ElementType>
  class ElementBase : public Vectorial<ElementType>,
                      public ElementType {
  public:
    typedef typename ElementType::type F;
 
    ElementBase(uint vdim, uint arrayWidth, std::string name = "ElementBase") :
      Vectorial<ElementType>(vdim, vdim), ElementType(), 
      T_(new TMatrix<F>(vdim)), name_(name) {
      DEBUGL(ElementConstr_D, "done.");
    }
    virtual ~ElementBase();
 
    virtual void put(ElementType& em, const int dim, const int b = 0);
 
    virtual vectorial::TMatrix<F>& T() const { return *T_; }
    virtual const Graphics<F>* graphics() const {
      if (! graphics_.get())
        graphics_.reset(new Graphics<F>());
      return graphics_.get();
    }
    static Graphics<F>* vecGraphics() {
      if (! graphics_.get())
        graphics_.reset(new Graphics<F>());
      return graphics_.get();
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    std::unique_ptr<vectorial::TMatrix<F> > T_;
  private:
    static std::unique_ptr<Graphics<F> > graphics_;
    std::string name_;
  };
 
  // *************************************************************** Element **
 
  template<class F>
  class Element : public ElementBase<concepts::Element<F> > {
  public:
    Element(uint vdim, uint arrayWidth=0) :
      ElementBase<concepts::Element<F> >(vdim, arrayWidth,
                                         "vectorial::Element") {}
    virtual ~Element();
  };
 
  // ******************************************************* ElementWithCell **
 
  template<class F>
  class ElementWithCell : public ElementBase<concepts::ElementWithCell<F> >,
                          public concepts::IntegrationCell {
  public:
    ElementWithCell(uint vdim, uint arrayWidth=0) :
      ElementBase<concepts::ElementWithCell<F> >(vdim, arrayWidth,
                                                 "vectorial::ElementWithCell"),
      cell_(0), intCell_(0) {}
    virtual ~ElementWithCell();
    
    virtual void put(concepts::ElementWithCell<F>& em,
                     const int dim, const int b = 0);
 
    virtual bool quadraturePoint(uint i, intPoint &p, intFormType form,
                                 bool localCoord) const;
    
    virtual const concepts::Cell& cell() const;
  private:
    const concepts::Cell* cell_;
    const concepts::IntegrationCell* intCell_;
  };
 
  // ********************************************************** BilinearForm **
 
  template<class F, class G = typename concepts::Realtype<F>::type>
  class BilinearForm : public concepts::BilinearForm<F,G>,
                       public Vectorial<concepts::BilinearForm<F,G> > {
  public:
    BilinearForm(const uint vdim1, const uint vdim2 = 0) : 
      Vectorial<concepts::BilinearForm<F,G> >(vdim1,
                                              vdim1*(vdim2 ? vdim2 : vdim1)),
      vdim2(vdim2 ? vdim2 : vdim1), deepCopies_(0)
    {
      DEBUGL(BilinearFormConstr_D, "done.");
    }
    BilinearForm(const BilinearForm& b);
    virtual ~BilinearForm();
    virtual BilinearForm* clone() const;
 
    virtual void put(concepts::BilinearForm<F,G> &bf, const int i, const int j);
    void putStore(concepts::BilinearForm<F,G>* bf, const int i, const int j);
    void putStore(concepts::RCP<concepts::BilinearForm<F,G> > bf,
                  const int i, const int j);
 
    void putData(concepts::Cloneable* data) { sharedData_.reset(data); }
 
    BilinearForm<F,G>& reset();
 
    virtual void operator ()(const concepts::Element<G>& elmX,
                             const concepts::Element<G>& elmY,
                             concepts::ElementMatrix<F>& em) const;
    virtual void operator ()(const concepts::Element<G>& elmV,
                             const concepts::Element<G>& elmU,
                             concepts::ElementMatrix<F>& em,
                             const concepts::ElementPair<G>& elmp) const;
 
    uint vdim2;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    concepts::Array<concepts::RCP<concepts::BilinearForm<F,G> > > deepCopies_;
    std::unique_ptr<concepts::Cloneable> sharedData_;
  };
 
} // namespace vectorial
 
namespace concepts {
 
  // ****************************************************************** Scan **
 
  template<class F>
  class Scan<vectorial::Element<F> > :
    public concepts::Scan<concepts::Element<F> > {
  public:  
    virtual vectorial::Element<F>& operator++(int) = 0;
  };
 
  template<class F>
  class Scan<vectorial::ElementWithCell<F> > :
    public concepts::Scan<concepts::ElementWithCell<F> > {
  public:  
    virtual vectorial::ElementWithCell<F>& operator++(int) = 0;
  };
 
} // namespace concepts
 
namespace vectorial {
 
  template<class F>
  class VectorialElementType {};
 
  template<class F>
  class VectorialElementType<concepts::Element<F> > {
  public:
    typedef Element<F> type;
  };
 
  template<class F>
  class VectorialElementType<concepts::ElementWithCell<F> > {
  public:
    typedef ElementWithCell<F> type;
  };
 
  // ***************************************************************** SpaceBase **
 
  template<class SpaceType>
  class SpaceBase : public SpaceType,
                    public Vectorial<SpaceType> {
  public:
    typedef 
    typename VectorialElementType<typename SpaceType::type>::type ElementType;
 
    SpaceBase(const uint vdim, const uint arrayWidth = 0,
          std::string name = "SpaceBase")
      : SpaceType()
      , Vectorial<SpaceType>(vdim, vdim)
      , dim_(0), nelm_(0), elm_(0), name_(name)  
    {
      DEBUGL(SpaceConstr_D, "done.");
    }
    SpaceBase(const SpaceBase & spcb)
      : Vectorial<SpaceType>(spcb.vdim(), spcb.vdim())
      , dim_(spcb.dim_) , nelm_(spcb.nelm_)
      , elm_(spcb.elm_) , name_(spcb.name_)
    {
    }
 
    virtual ~SpaceBase();
 
    virtual uint dim() const { return dim_; }
    virtual uint nelm() const { return nelm_; }
    typedef concepts::Scan<ElementType> Scanner;
    virtual Scanner* scan() const;
 
    virtual void put(SpaceType &spc, const int i = 0, const int j = 0);
 
    void rebuild();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uint dim_, nelm_;
    concepts::Joiner<ElementType*, 1>* elm_;
    std::string name_;
    void buildElm_(const SpaceType& spc, uint idx);
  };
 
  // ***************************************************************** Space **
 
  template<class F>
  class Space : public SpaceBase<concepts::Space<F> > {
  public:
    Space(const uint vdim, const uint arrayWidth=0)
      : SpaceBase<concepts::Space<F> >(vdim, arrayWidth, "Space") {}
    virtual ~Space();
  };
 
  // ********************************************************** SpaceOnCells **
 
  template<class F>
  class SpaceOnCells : public SpaceBase<concepts::SpaceOnCells<F> > {
  public:
    SpaceOnCells(const uint vdim, const uint arrayWidth=0)
      : SpaceBase<concepts::SpaceOnCells<F> >(vdim, arrayWidth, "SpaceOnCells") {}
    virtual ~SpaceOnCells();
  };
 
  // *********************************************************** ElementPair **
  
  template<class F>
  class ElementPair 
    : public concepts::ElementPair<F>,
      public Vectorial<concepts::ElementPair<F> > {
  public:
    ElementPair(const Element<F>& elm1, const Element<F>& elm2) 
      : Vectorial<concepts::ElementPair<F> >(elm2.vdim(),
                                             elm1.vdim()*elm2.vdim()), 
        elm1_(elm1), elm2_(elm2) { }
    virtual ~ElementPair() { }
    virtual void put(concepts::ElementPair<F>& ep,
                     const int iV, const int iU);
    virtual const concepts::Element<F>& elm1() const;
    virtual const concepts::Element<F>& elm2() const;
  private:
    const Element<F>& elm1_;
    const Element<F>& elm2_;
  };
 
} // namespace vectorial
 
#endif // vectorialThings_hh