element.hh

Go to the documentation of this file.

 
#ifndef spcElement_hh
#define spcElement_hh
 
#if __GNUC__ == 2
#  include <float.h>
#  define EPS DBL_EPSILON
#else
#  include <limits>
#  define EPS std::numeric_limits<double>::epsilon()
#endif
 
#include <cstring>
#include "basics/typedefs.hh"
#include "basics/outputOperator.hh"
#include "basics/exceptions.hh"
#include "basics/vectorsMatricesForward.hh"
#include "toolbox/array.hh"
#include "geometry/cell.hh"
 
#include "basics/debug.hh"
 
// debugging
#define ElementMatrixAppl_D 0
 
namespace concepts {
 
  // forward declaration
  template<class F>
  class ElementGraphics;
 
  template<class F>
  class TMatrixBase;
 
  class Cell;
 
  // *************************************************************** Element **
 
  template<class F>
  class Element : public virtual OutputOperator {
  public:
    typedef F type;
 
    Element() : tag_(0) {}
    virtual ~Element() {}
 
    virtual const TMatrixBase<F>& T() const = 0;
 
    virtual const ElementGraphics<F>* graphics() const { return 0; }
 
    inline uint& tag() { return tag_; }
  private:
    uint tag_;
  };
 
  // ******************************************************* ElementWithCell **
 
  template<typename F>
  class ElementWithCell : public Element<F> {
  public:
    typedef F type;
 
    virtual const Cell& cell() const = 0;
 
    Real3d elemMap(const Real coord_local) const {
      return cell().elemMap(coord_local);
    }
 
    Real3d elemMap(const Real2d& coord_local) const {
      return cell().elemMap(coord_local);
    }
 
    Real3d elemMap(const Real3d& coord_local) const {
      return cell().elemMap(coord_local);
    }
 
    virtual const TMatrixBase<F>& T() const = 0;
 
  };
 
 
  // ****************************************************** ElementAndFacette **
 
  template<class F>
  class ElementAndFacette : public virtual OutputOperator {
  public:
    ElementAndFacette(const F* elm, const uint k)
      : elm(elm), k(k) {}
    const F* elm;
    uint k;
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      conceptsAssert(elm, concepts::Assertion());
      return os << concepts::typeOf(*this)<<"(" << *elm << ", k = " << k << ")";
    }
  };
 
 
  // ***************************************************** ElementMatrixBase **
 
  template<class F>
  class ElementMatrixBase {
  public:
    typedef F  value_type;
    
    ElementMatrixBase() : data_(0), m_(0), n_(0), t_(false) {}
 
    uint m() const { return t_ ? n_ : m_; }
    uint n() const { return t_ ? m_ : n_; }
 
    inline const F& operator()(const uint i, const uint j) const {
      DEBUGL(ElementMatrixAppl_D, '(' << (t_ ? n_ : m_) << 'x'
       << (t_ ? m_ : n_) << ") - " << (t_ ? j : i) << ", "
       << (t_ ? i : j));
      conceptsAssert(t_ ? j < m_ : i < m_, Assertion());
      conceptsAssert(t_ ? i < n_ : j < n_, Assertion());
      return t_ ? data_[j * n_ + i] : data_[i * n_ + j]; 
    }
 
    uint size() const { return data_.size(); }
 
    operator const F*() const { return (const F*)data_; }
 
    const Array<F>& getData() const { return data_; }
 
    bool isTranspose() const { return t_; }
 
    bool storeMatlab(const std::string filename, std::string name = "",
                     bool append = false) const;
  protected:
    Array<F> data_;
    uint m_;
    uint n_;
    bool t_;
  };
 
  // ********************************************************* ElementMatrix **
 
  template<class F>
  class ElementMatrix : public ElementMatrixBase<F> {
  public:
    typedef F  value_type;
 
    ElementMatrix(const uint m = 0, const uint n = 0);
 
    ElementMatrix(const uint m, const uint n, const F* data);
 
    virtual ~ElementMatrix() {}
 
    ElementMatrix<F>& operator=(const ElementMatrixBase<F>& A);
    void add(const ElementMatrix<F>& A, uint offm = 0, uint offn = 0);
 
    inline const F& operator()(const uint i, const uint j) const {
      return ElementMatrixBase<F>::operator()(i, j);
    }
    inline F& operator()(const uint i, const uint j) {
      DEBUGL(ElementMatrixAppl_D, '(' << this->m() << 'x' << this->n() << ") - " 
             << (this->t_ ? j : i) << ", " << (this->t_ ? i : j));
      conceptsAssert(this->t_ ? j < this->m_ : i < this->m_, Assertion());
      conceptsAssert(this->t_ ? i < this->n_ : j < this->n_, Assertion());
      return this->t_ ? this->data_[j * this->n_ + i] 
        : this->data_[i * this->n_ + j];
    }
 
    template<class G>
    ElementMatrix<F>& operator=(const ElementMatrix<G>& other) 
    {
      this->data_ = other.getData();
      i_ = other.getIndex();
      this->m_ = other.isTranspose() ? other.n() : other.m();
      this->n_ = other.isTranspose() ? other.m() : other.n();
      this->t_ = other.isTranspose();
      return *this; 
    }
 
    virtual ElementMatrix<F>& operator=(const F& v) 
    {
      this->data_ = v;
      return *this; 
    }
 
    virtual void resize(uint m, uint n);
 
    template<uint dim>
    void writeColumn(uint col, const Point<F,dim>& val);
 
    void writeColumn(uint col, const Array<F>& val);
 
    void writeColumn(uint col, const F val);
 
    bool transpose() { return this->t_ ^= 1; }
    
    void setTranspose(bool t) { this->t_ = t; }
 
    int getIndex() const { return i_; }
 
    void setIndex(int i) { i_ = i; }
 
    void zeros() { std::memset((F*)this->data_, 0, this->n_*this->m_*sizeof(F)); }
 
    operator F*() { return (F*)this->data_; }
 
    ElementMatrix<F>& operator*=(const F n) {
      this->data_ *= n;
      return *this;
    }
 
    uint row() const { return this->t_ ? i_ % this->n_ : i_ / this->n_; }
    uint col() const { return this->t_ ? i_ / this->n_ : i_ % this->n_; }
 
    void compress(const Real threshold = EPS);
 
    std::ostream& info(std::ostream& os) const;
  private:
    uint i_;
  };
 
  template<class F>
  template<uint dim>
  void ElementMatrix<F>::writeColumn(uint col, const Point<F,dim>& val)
  {
    conceptsAssert(this->m() >= dim, Assertion());
    for(uint i = 0; i < dim; ++i)
      (*this)(i, col) = val[i];
  }
 
  template<class F>
  std::ostream& operator<<(std::ostream& os, const ElementMatrix<F>& o) {
#ifdef DEBUG
    os << std::flush;
#endif
    return o.info(os);
  }
 
} // namespace concepts
 
#endif // spcElement_hh