#include "NodeCentricRepMatrix.h"

using namespace OBJREC;
using namespace std;

/** Member-Methodes **/
NodeCentricRepMatrix::NodeCentricRepMatrix ( const uint m, const uint n, const uint e ) : std::valarray<double> ( m * n * e ) {
  this->m = m;
  this->n = n;
  this->e = e;
}
NodeCentricRepMatrix::NodeCentricRepMatrix ( const NodeCentricRepMatrix& matrix ) : std::valarray<double> ( matrix ) {
  this->m = matrix.getM();
  this->n = matrix.getN();
  this->e = matrix.getE();
}


NodeCentricRepMatrix::~NodeCentricRepMatrix() {}

double NodeCentricRepMatrix::at ( uint m, uint n, uint e ) const {
  checkRange ( m, n, e );
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}

double& NodeCentricRepMatrix::at ( uint m, uint n, uint e ) {
  checkRange ( m, n, e );
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}

double NodeCentricRepMatrix::atQuick ( uint m, uint n, uint e ) const
{
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}

double& NodeCentricRepMatrix::atQuick ( uint m, uint n, uint e )
{
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}


void NodeCentricRepMatrix::checkRange ( uint m, uint n, uint e ) const {
  if ( m >= ( *this ).m || n >= ( *this ).n || e >= ( *this ).e )
    throw range_error ( "NodeCentricRepMatrix: You try to access an element outside the matrix!" );
}

NodeCentricRepMatrix& NodeCentricRepMatrix::operator= ( const NodeCentricRepMatrix & matrix ) {
  // Selbstkopie verhindern
  if ( this == &matrix ) return ( *this );

  if ( ( *this ).getM() != matrix.getM() ||
       ( *this ).getN() != matrix.getN() ||
       ( *this ).getE() != matrix.getE() ) {
    throw range_error ( "NodeCentricRepMatrix: Assigment not possible. Matrix-dim does not match." );
  }

  for ( uint i = 0; i < ( *this ).getN(); i++ ) {
    for ( uint j = 0; j < ( *this ).getM(); j++ ) {
      for ( uint e = 0; e < ( *this ).getE(); e++ ) {
        ( *this ) ( j, i, e ) = matrix ( j, i, e );
      }
    }
  }
  return ( *this );
}

NodeCentricRepMatrix& NodeCentricRepMatrix::operator-= ( NodeCentricRepMatrix const & rhs )
{
  if ( ( *this ).getM() != rhs.getM() ||
       ( *this ).getN() != rhs.getN() ||
       ( *this ).getE() != rhs.getE() )
  {
    throw range_error ( "NodeCentricRepMatrix: Minus not possible. Matrix-dim does not match." );
  }
#pragma omp parallel for
  for ( uint i = 0; i < ( *this ).getN(); i++ )
  {
    for ( uint j = 0; j < ( *this ).getM(); j++ )
    {
      for ( uint e = 0; e < ( *this ).getE(); e++ )
      {
        ( *this ) ( j, i, e ) -= rhs ( j, i, e );
      }
    }
  }

  return ( *this );
}

NodeCentricRepMatrix NodeCentricRepMatrix::operator- ( NodeCentricRepMatrix const& rhs )
{
  NodeCentricRepMatrix tmp ( ( *this ) );
  tmp -= rhs;
  return tmp;
}

double NodeCentricRepMatrix::operator() ( const uint m, const uint n, const uint e ) const
{
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}

double& NodeCentricRepMatrix::operator() ( const uint m, const uint n, const uint e )
{
  return valarray<double>::operator[] ( m * ( ( *this ).n * ( *this ).e ) + n * ( *this ).e + e );
}

bool NodeCentricRepMatrix::changes ( const NodeCentricRepMatrix& matrix, double epsilon ) const {
  if ( ( *this ).getM() != matrix.getM() ||
       ( *this ).getN() != matrix.getN() ||
       ( *this ).getE() != matrix.getE() ) {
    throw range_error ( "NodeCentricRepMatrix: Calculate difference not possible. Matrix-dim does not match." );
  }

  bool retValue = false;

//#pragma omp parallel for
  for ( uint i = 0; i < ( *this ).getN(); i++ ) {
    for ( uint j = 0; j < ( *this ).getM(); j++ ) {
      for ( uint e = 0; e < ( *this ).getE(); e++ ) {
        if ( abs ( ( *this ) ( j, i, e ) - matrix ( j, i, e ) ) > epsilon )
        {
          return true;
        }

      }
    }
  }
  return retValue;
}