space/element.hh

Go to the documentation of this file.

/* Elements
 * Elements to build a space
 */

#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 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:
    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);
    }

  };


  // ****************************************************** ElementAndFacette **

  template<class F>
  class ElementAndFacette : public concepts::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 << "ElementAndFacette(" << *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 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);

    ElementMatrix<F>& operator=(const ElementMatrixBase<F>& A);
    void add(const ElementMatrix<F>& A, uint offm = 0, uint offn = 0);

    inline 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; 
    }

    void resize(uint m, uint n);

    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>
  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