Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Attributes | Friends | List of all members
moab::Matrix3 Class Reference

#include <Matrix3.hpp>

Public Member Functions

 Matrix3 ()
 
 Matrix3 (double diagonal)
 
 Matrix3 (const CartVect &diagonal)
 
 Matrix3 (const std::vector< double > &diagonal)
 
 Matrix3 (double v00, double v01, double v02, double v10, double v11, double v12, double v20, double v21, double v22)
 
 Matrix3 (const Matrix3 &f)
 
template<typename Vector >
 Matrix3 (const Vector &row0, const Vector &row1, const Vector &row2, const bool isRow)
 
 Matrix3 (const double v[size])
 
void copyto (double v[Matrix3::size])
 
Matrix3 & operator= (const Matrix3 &m)
 
Matrix3 & operator= (const double v[size])
 
double * operator[] (unsigned i)
 
const double * operator[] (unsigned i) const
 
double & operator() (unsigned r, unsigned c)
 
double operator() (unsigned r, unsigned c) const
 
double & operator() (unsigned i)
 
double operator() (unsigned i) const
 
double * array ()
 
const double * array () const
 
Matrix3 & operator+= (const Matrix3 &m)
 
Matrix3 & operator-= (const Matrix3 &m)
 
Matrix3 & operator*= (double s)
 
Matrix3 & operator/= (double s)
 
Matrix3 & operator*= (const Matrix3 &m)
 
bool is_symmetric ()
 
bool is_positive_definite ()
 
template<typename Vector >
ErrorCode eigen_decomposition (Vector &evals, Matrix3 &evecs)
 
template<typename Vector >
ErrorCode eigen_decomposition_native (Vector &eigenvalues, Matrix3 &eigenvectors)
 
void transpose_inplace ()
 
Matrix3 transpose () const
 
template<typename Vector >
void copycol (int index, Vector &vol)
 
void swapcol (int srcindex, int destindex)
 
template<typename Vector >
Vector vcol (int index) const
 
void colscale (int index, double scale)
 
void rowscale (int index, double scale)
 
CartVect col (int index) const
 
CartVect row (int index) const
 
CartVect diagonal () const
 
double trace () const
 
double determinant () const
 
Matrix3 inverse () const
 
bool invert ()
 
double subdet (int r, int c) const
 
void print (std::ostream &s) const
 

Static Public Member Functions

static Matrix3 identity ()
 

Static Public Attributes

static const int size = 9
 

Private Attributes

double _mat [size]
 

Friends

Matrix3 operator+ (const Matrix3 &a, const Matrix3 &b)
 
Matrix3 operator- (const Matrix3 &a, const Matrix3 &b)
 
Matrix3 operator* (const Matrix3 &a, const Matrix3 &b)
 

Detailed Description

Definition at line 173 of file Matrix3.hpp.

Constructor & Destructor Documentation

◆ Matrix3() [1/8]

moab::Matrix3::Matrix3 ( )
inline

Definition at line 187 of file Matrix3.hpp.

188  {
189 #ifdef MOAB_HAVE_EIGEN3
190  _mat.fill( 0.0 );
191 #else
192  MOAB_DMEMZERO( _mat, Matrix3::size );
193 #endif
194  }

References _mat, MOAB_DMEMZERO, and size.

Referenced by identity(), inverse(), and transpose().

◆ Matrix3() [2/8]

moab::Matrix3::Matrix3 ( double  diagonal)
inline

Definition at line 202 of file Matrix3.hpp.

203  {
204 #ifdef MOAB_HAVE_EIGEN3
205  _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal;
206 #else
207  MOAB_DMEMZERO( _mat, Matrix3::size );
208  _mat[0] = _mat[4] = _mat[8] = diagonal;
209 #endif
210  }

References _mat, diagonal(), MOAB_DMEMZERO, and size.

◆ Matrix3() [3/8]

moab::Matrix3::Matrix3 ( const CartVect &  diagonal)
inline

Definition at line 212 of file Matrix3.hpp.

213  {
214 #ifdef MOAB_HAVE_EIGEN3
215  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
216 #else
217  MOAB_DMEMZERO( _mat, Matrix3::size );
218  _mat[0] = diagonal[0];
219  _mat[4] = diagonal[1];
220  _mat[8] = diagonal[2];
221 #endif
222  }

References _mat, diagonal(), MOAB_DMEMZERO, and size.

◆ Matrix3() [4/8]

moab::Matrix3::Matrix3 ( const std::vector< double > &  diagonal)
inline

Definition at line 229 of file Matrix3.hpp.

230  {
231 #ifdef MOAB_HAVE_EIGEN3
232  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
233 #else
234  MOAB_DMEMZERO( _mat, Matrix3::size );
235  _mat[0] = diagonal[0];
236  _mat[4] = diagonal[1];
237  _mat[8] = diagonal[2];
238 #endif
239  }

References _mat, diagonal(), MOAB_DMEMZERO, and size.

◆ Matrix3() [5/8]

moab::Matrix3::Matrix3 ( double  v00,
double  v01,
double  v02,
double  v10,
double  v11,
double  v12,
double  v20,
double  v21,
double  v22 
)
inline

Definition at line 241 of file Matrix3.hpp.

250  {
251 #ifdef MOAB_HAVE_EIGEN3
252  _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22;
253 #else
254  MOAB_DMEMZERO( _mat, Matrix3::size );
255  _mat[0] = v00;
256  _mat[1] = v01;
257  _mat[2] = v02;
258  _mat[3] = v10;
259  _mat[4] = v11;
260  _mat[5] = v12;
261  _mat[6] = v20;
262  _mat[7] = v21;
263  _mat[8] = v22;
264 #endif
265  }

References _mat, MOAB_DMEMZERO, and size.

◆ Matrix3() [6/8]

moab::Matrix3::Matrix3 ( const Matrix3 &  f)
inline

Definition at line 268 of file Matrix3.hpp.

269  {
270 #ifdef MOAB_HAVE_EIGEN3
271  _mat = f._mat;
272 #else
273  memcpy( _mat, f._mat, size * sizeof( double ) );
274 #endif
275  }

References _mat, and size.

◆ Matrix3() [7/8]

template<typename Vector >
moab::Matrix3::Matrix3 ( const Vector &  row0,
const Vector &  row1,
const Vector &  row2,
const bool  isRow 
)
inline

Definition at line 279 of file Matrix3.hpp.

280  {
281 #ifdef MOAB_HAVE_EIGEN3
282  if( isRow )
283  {
284  _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2];
285  }
286  else
287  {
288  _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2];
289  }
290 #else
291  MOAB_DMEMZERO( _mat, Matrix3::size );
292  if( isRow )
293  {
294  _mat[0] = row0[0];
295  _mat[1] = row0[1];
296  _mat[2] = row0[2];
297  _mat[3] = row1[0];
298  _mat[4] = row1[1];
299  _mat[5] = row1[2];
300  _mat[6] = row2[0];
301  _mat[7] = row2[1];
302  _mat[8] = row2[2];
303  }
304  else
305  {
306  _mat[0] = row0[0];
307  _mat[1] = row1[0];
308  _mat[2] = row2[0];
309  _mat[3] = row0[1];
310  _mat[4] = row1[1];
311  _mat[5] = row2[1];
312  _mat[6] = row0[2];
313  _mat[7] = row1[2];
314  _mat[8] = row2[2];
315  }
316 #endif
317  }

References _mat, MOAB_DMEMZERO, and size.

◆ Matrix3() [8/8]

moab::Matrix3::Matrix3 ( const double  v[size])
inline

Definition at line 323 of file Matrix3.hpp.

324  {
325 #ifdef MOAB_HAVE_EIGEN3
326  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
327 #else
328  memcpy( _mat, v, size * sizeof( double ) );
329 #endif
330  }

References _mat, and size.

Member Function Documentation

◆ array() [1/2]

double* moab::Matrix3::array ( )
inline

Definition at line 416 of file Matrix3.hpp.

417  {
418 #ifdef MOAB_HAVE_EIGEN3
419  return _mat.data();
420 #else
421  return _mat;
422 #endif
423  }

References _mat.

Referenced by moab::EvalSet::evaluate_reverse(), and moab::LinearTet::evaluate_reverse().

◆ array() [2/2]

const double* moab::Matrix3::array ( ) const
inline

Definition at line 425 of file Matrix3.hpp.

426  {
427 #ifdef MOAB_HAVE_EIGEN3
428  return _mat.data();
429 #else
430  return _mat;
431 #endif
432  }

References _mat.

◆ col()

CartVect moab::Matrix3::col ( int  index) const
inline

Definition at line 971 of file Matrix3.hpp.

972  {
973  assert( index < Matrix3::size );
974 #ifdef MOAB_HAVE_EIGEN3
975  Eigen::Vector3d mvec = _mat.col( index );
976  return CartVect( mvec[0], mvec[1], mvec[2] );
977 #else
978  switch( index )
979  {
980  case 0:
981  return CartVect( _mat[0], _mat[3], _mat[6] );
982  case 1:
983  return CartVect( _mat[1], _mat[4], _mat[7] );
984  case 2:
985  return CartVect( _mat[2], _mat[5], _mat[8] );
986  }
987  return CartVect( 0.0 );
988 #endif
989  }

References _mat, and size.

Referenced by moab::OrientedBox::area(), moab::OrientedBox::axis(), moab::box_from_axes(), moab::OrientedBox::closest_location_in_box(), moab::OrientedBox::contained(), moab::OrientedBox::dimensions(), eigen_decomposition_native(), moab::OrientedBox::inner_radius(), moab::OrientedBox::inner_radius_squared(), moab::OrientedBox::make_hex(), moab::operator<<(), moab::OrientedBox::OrientedBox(), moab::OrientedBox::outer_radius(), moab::OrientedBox::outer_radius_squared(), moab::OrientedBox::scaled_axis(), and moab::OrientedBox::volume().

◆ colscale()

void moab::Matrix3::colscale ( int  index,
double  scale 
)
inline

Definition at line 917 of file Matrix3.hpp.

918  {
919  assert( index < Matrix3::size );
920 #ifdef MOAB_HAVE_EIGEN3
921  _mat.col( index ) *= scale;
922 #else
923  switch( index )
924  {
925  case 0:
926  _mat[0] *= scale;
927  _mat[3] *= scale;
928  _mat[6] *= scale;
929  break;
930  case 1:
931  _mat[1] *= scale;
932  _mat[4] *= scale;
933  _mat[7] *= scale;
934  break;
935  case 2:
936  _mat[2] *= scale;
937  _mat[5] *= scale;
938  _mat[8] *= scale;
939  break;
940  }
941 #endif
942  }

References _mat, and size.

Referenced by moab::box_from_axes(), and moab::OrientedBox::order_axes_by_length().

◆ copycol()

template<typename Vector >
void moab::Matrix3::copycol ( int  index,
Vector &  vol 
)
inline

Definition at line 823 of file Matrix3.hpp.

824  {
825 #ifdef MOAB_HAVE_EIGEN3
826  _mat.col( index ).swap( vol );
827 #else
828  switch( index )
829  {
830  case 0:
831  _mat[0] = vol[0];
832  _mat[3] = vol[1];
833  _mat[6] = vol[2];
834  break;
835  case 1:
836  _mat[1] = vol[0];
837  _mat[4] = vol[1];
838  _mat[7] = vol[2];
839  break;
840  case 2:
841  _mat[2] = vol[0];
842  _mat[5] = vol[1];
843  _mat[8] = vol[2];
844  break;
845  }
846 #endif
847  }

References _mat.

◆ copyto()

void moab::Matrix3::copyto ( double  v[Matrix3::size])
inline

Definition at line 332 of file Matrix3.hpp.

333  {
334 #ifdef MOAB_HAVE_EIGEN3
335  std::copy( _mat.data(), _mat.data() + size, v );
336 #else
337  memcpy( v, _mat, size * sizeof( double ) );
338 #endif
339  }

References _mat, and size.

Referenced by moab::LinearTet::initFcn(), and moab::LinearTri::initFcn().

◆ determinant()

double moab::Matrix3::determinant ( ) const
inline

Definition at line 1033 of file Matrix3.hpp.

1034  {
1035 #ifdef MOAB_HAVE_EIGEN3
1036  return _mat.determinant();
1037 #else
1038  return ( _mat[0] * _mat[4] * _mat[8] + _mat[1] * _mat[5] * _mat[6] + _mat[2] * _mat[3] * _mat[7] -
1039  _mat[0] * _mat[5] * _mat[7] - _mat[1] * _mat[3] * _mat[8] - _mat[2] * _mat[4] * _mat[6] );
1040 #endif
1041  }

References _mat.

Referenced by moab::Element::Map::det_ijacobian(), moab::Element::Map::det_jacobian(), moab::EvalSet::evaluate_reverse(), moab::LinearTet::evaluate_reverse(), moab::LinearTri::evaluate_reverse(), moab::Element::Map::ievaluate(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), moab::Element::SpectralHex::integrate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::integrate_scalar_field(), moab::ElemUtil::integrate_trilinear_hex(), moab::LinearHex::integrateFcn(), moab::LinearQuad::integrateFcn(), inverse(), invert(), moab::AffineXform::reflection(), moab::AffineXform::scale(), moab::Element::LinearTet::set_vertices(), moab::Element::LinearTri::set_vertices(), moab::GeomUtil::VolMap::solve_inverse(), moab::ElemUtil::VolMap::solve_inverse(), and test_spectral_hex().

◆ diagonal()

CartVect moab::Matrix3::diagonal ( ) const
inline

Definition at line 1011 of file Matrix3.hpp.

1012  {
1013 #ifdef MOAB_HAVE_EIGEN3
1014  return CartVect( _mat( 0, 0 ), _mat( 1, 1 ), _mat( 2, 2 ) );
1015 #else
1016  return CartVect( _mat[0], _mat[4], _mat[8] );
1017 #endif
1018  }

References _mat.

Referenced by Matrix3().

◆ eigen_decomposition()

template<typename Vector >
ErrorCode moab::Matrix3::eigen_decomposition ( Vector &  evals,
Matrix3 &  evecs 
)
inline

Definition at line 635 of file Matrix3.hpp.

636  {
637  const bool bisSymmetric = this->is_symmetric();
638  if( !bisSymmetric )
639  {
640  MB_CHK_SET_ERR( MB_FAILURE, "Unsymmetric matrix implementation with Matrix3 are currently not supported." );
641  }
642 #if defined( MOAB_HAVE_EIGEN3 )
643  Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > eigensolver( this->_mat );
644  if( eigensolver.info() != Eigen::Success ) return MB_FAILURE;
645  const Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d >::RealVectorType& e3evals = eigensolver.eigenvalues();
646  evals[0] = e3evals( 0 );
647  evals[1] = e3evals( 1 );
648  evals[2] = e3evals( 2 );
649  evecs._mat = eigensolver.eigenvectors(); //.col(1)
650  return MB_SUCCESS;
651 
652 #elif defined( MOAB_HAVE_LAPACK )
653  return eigen_decomposition_lapack( bisSymmetric, evals, evecs );
654 
655  // Vector evals_native;
656  // Matrix3 evecs_native;
657  // eigen_decomposition_native( evals_native, evecs_native );
658  // // std::cout << "LAPACK: " << evals << " NATIVE: " << evals_native << " Error: " << evals-evals_native << std::endl;
659  // return MB_SUCCESS;
660 #else
661  return eigen_decomposition_native( evals, evecs );
662 #endif
663  }

References _mat, eigen_decomposition_native(), is_symmetric(), MB_CHK_SET_ERR, and MB_SUCCESS.

Referenced by moab::OrientedBox::compute_from_covariance_data(), and moab::OrientedBox::compute_from_vertices().

◆ eigen_decomposition_native()

template<typename Vector >
ErrorCode moab::Matrix3::eigen_decomposition_native ( Vector &  eigenvalues,
Matrix3 &  eigenvectors 
)
inline

Definition at line 666 of file Matrix3.hpp.

667  {
668  constexpr int n = 3;
669  constexpr double EPSILON = 1e-14; // Tolerance for stopping the iteration
670  constexpr int maxiters = 500;
671  Matrix3 matrix = *this;
672 
673  // Initialize the eigenvector matrix as the identity matrix
674  eigenvectors = moab::Matrix3::identity();
675 
676  // Function to perform a Jacobi rotation
677  auto jacobiRotate = []( moab::Matrix3& A, moab::Matrix3& V, int p, int q ) {
678  double theta, t, c, s;
679  theta = ( A( q, q ) - A( p, p ) ) / ( 2.0 * A( p, q ) );
680  t = ( theta >= 0 ) ? 1.0 / ( theta + std::sqrt( 1.0 + theta * theta ) )
681  : 1.0 / ( theta - std::sqrt( 1.0 + theta * theta ) );
682  c = 1.0 / std::sqrt( 1.0 + t * t );
683  s = t * c;
684 
685  double app = A( p, p ), aqq = A( q, q ), apq = A( p, q );
686  A( p, p ) = c * c * app - 2.0 * c * s * apq + s * s * aqq;
687  A( q, q ) = s * s * app + 2.0 * c * s * apq + c * c * aqq;
688  A( p, q ) = A( q, p ) = 0.0; // Zero the off-diagonal element
689 
690  // Update other elements
691  for( int i = 0; i < 3; i++ )
692  {
693  if( i != p && i != q )
694  {
695  double aip = A( i, p ), aiq = A( i, q );
696  A( i, p ) = A( p, i ) = c * aip - s * aiq;
697  A( i, q ) = A( q, i ) = s * aip + c * aiq;
698  }
699  }
700 
701  // Update the eigenvector matrix
702  for( int i = 0; i < 3; i++ )
703  {
704  double vip = V( i, p ), viq = V( i, q );
705  V( i, p ) = c * vip - s * viq;
706  V( i, q ) = s * vip + c * viq;
707  }
708  };
709 
710  // Iterate to apply Jacobi rotations to find
711  // eigenvalues and eigenvectors of a symmetric 3x3 matrix
712  bool converged = false;
713  for( int iter = 0; iter < maxiters; ++iter )
714  {
715  // Find the largest off-diagonal element
716  int p = 0, q = 1;
717  double maxOffDiag = std::abs( matrix( p, q ) );
718  for( int i = 0; i < n; ++i )
719  {
720  for( int j = i + 1; j < n; ++j )
721  {
722  if( std::abs( matrix( i, j ) ) > maxOffDiag )
723  {
724  maxOffDiag = std::abs( matrix( i, j ) );
725  p = i;
726  q = j;
727  }
728  }
729  }
730 
731  // If the largest off-diagonal element is smaller than tolerance, stop
732  if( maxOffDiag < EPSILON )
733  {
734  converged = true;
735  break;
736  }
737 
738  // Apply Jacobi rotation to zero out matrix[p][q]
739  if( fabs( matrix( p, q ) ) > EPSILON ) jacobiRotate( matrix, eigenvectors, p, q );
740  }
741 
742  // The diagonal elements of A are the eigenvalues
743  for( int i = 0; i < n; ++i )
744  {
745  // Extract eigenvalues (diagonal elements of the matrix)
746  eigenvalues[i] = matrix( i, i );
747  eigenvectors.col( i ).normalize();
748  }
749 
750  if( !converged )
751  {
752  // if iterative scheme did not converge, show warning and return with failure
753  std::cerr << "Jacobi method did not converge within " << maxiters << " iterations\n";
754  }
755  else
756  {
757  // iterative method converged; let us sort the eigenvalues and eigenvectors in increasing order
758  auto sortIndices = []( const CartVect& vec ) -> std::array< int, 3 > {
759  // Initialize the index array with values 0, 1, 2
760  std::array< int, 3 > indices = { 0, 1, 2 };
761 
762  // Sort the indices based on the values in the original vector
763  std::sort( indices.begin(), indices.end(), [&]( int i, int j ) { return vec[i] < vec[j]; } );
764 
765  return indices;
766  };
767 
768  auto newIndices = sortIndices( eigenvalues );
769 
770  // eigenvectors.transpose_inplace();
771  CartVect teigenvalues = eigenvalues;
772  Matrix3 teigenvectors = eigenvectors;
773  for( int i = 0; i < 3; i++ )
774  {
775  eigenvalues[i] = teigenvalues[newIndices[i]];
776  eigenvectors( i, 0 ) = teigenvectors( i, newIndices[0] );
777  eigenvectors( i, 1 ) = teigenvectors( i, newIndices[1] );
778  eigenvectors( i, 2 ) = teigenvectors( i, newIndices[2] );
779  }
780  }
781 
782  return ( converged ? MB_SUCCESS : MB_FAILURE );
783  }

References col(), EPSILON, identity(), MB_SUCCESS, and moab::CartVect::normalize().

Referenced by eigen_decomposition().

◆ identity()

static Matrix3 moab::Matrix3::identity ( )
inlinestatic

Definition at line 785 of file Matrix3.hpp.

786  {
787  // Identity matrix
788  return Matrix3( 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 );
789  }

References Matrix3().

Referenced by eigen_decomposition_native().

◆ inverse()

Matrix3 moab::Matrix3::inverse ( ) const
inline

Definition at line 1043 of file Matrix3.hpp.

1044  {
1045 #ifdef MOAB_HAVE_EIGEN3
1046  return Matrix3( _mat.inverse() );
1047 #else
1048  // return Matrix::compute_inverse( *this, this->determinant() );
1049  Matrix3 m( 0.0 );
1050  const double d_determinant = 1.0 / this->determinant();
1051  m._mat[0] = d_determinant * ( _mat[4] * _mat[8] - _mat[5] * _mat[7] );
1052  m._mat[1] = d_determinant * ( _mat[2] * _mat[7] - _mat[8] * _mat[1] );
1053  m._mat[2] = d_determinant * ( _mat[1] * _mat[5] - _mat[4] * _mat[2] );
1054  m._mat[3] = d_determinant * ( _mat[5] * _mat[6] - _mat[8] * _mat[3] );
1055  m._mat[4] = d_determinant * ( _mat[0] * _mat[8] - _mat[6] * _mat[2] );
1056  m._mat[5] = d_determinant * ( _mat[2] * _mat[3] - _mat[5] * _mat[0] );
1057  m._mat[6] = d_determinant * ( _mat[3] * _mat[7] - _mat[6] * _mat[4] );
1058  m._mat[7] = d_determinant * ( _mat[1] * _mat[6] - _mat[7] * _mat[0] );
1059  m._mat[8] = d_determinant * ( _mat[0] * _mat[4] - _mat[3] * _mat[1] );
1060  return m;
1061 #endif
1062  }

References _mat, determinant(), and Matrix3().

Referenced by moab::Element::Map::det_ijacobian(), moab::EvalSet::evaluate_reverse(), moab::LinearTet::evaluate_reverse(), moab::LinearTri::evaluate_reverse(), moab::Element::Map::ievaluate(), moab::Element::Map::ijacobian(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), moab::AffineXform::inverse(), moab::Element::LinearTet::set_vertices(), moab::Element::LinearTri::set_vertices(), moab::element_utility::Linear_hex_map::solve_inverse(), moab::GeomUtil::VolMap::solve_inverse(), and moab::ElemUtil::VolMap::solve_inverse().

◆ invert()

bool moab::Matrix3::invert ( )
inline

Definition at line 1064 of file Matrix3.hpp.

1065  {
1066  bool invertible = false;
1067  double d_determinant;
1068 #ifdef MOAB_HAVE_EIGEN3
1069  Eigen::Matrix3d invMat;
1070  _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible );
1071  if( !Util::is_finite( d_determinant ) ) return false;
1072  _mat = invMat;
1073  return invertible;
1074 #else
1075  d_determinant = this->determinant();
1076  if( d_determinant > 1e-13 ) invertible = true;
1077  d_determinant = 1.0 / d_determinant; // invert the determinant
1078  std::vector< double > _m;
1079  _m.assign( _mat, _mat + size );
1080  _mat[0] = d_determinant * ( _m[4] * _m[8] - _m[5] * _m[7] );
1081  _mat[1] = d_determinant * ( _m[2] * _m[7] - _m[8] * _m[1] );
1082  _mat[2] = d_determinant * ( _m[1] * _m[5] - _m[4] * _m[2] );
1083  _mat[3] = d_determinant * ( _m[5] * _m[6] - _m[8] * _m[3] );
1084  _mat[4] = d_determinant * ( _m[0] * _m[8] - _m[6] * _m[2] );
1085  _mat[5] = d_determinant * ( _m[2] * _m[3] - _m[5] * _m[0] );
1086  _mat[6] = d_determinant * ( _m[3] * _m[7] - _m[6] * _m[4] );
1087  _mat[7] = d_determinant * ( _m[1] * _m[6] - _m[7] * _m[0] );
1088  _mat[8] = d_determinant * ( _m[0] * _m[4] - _m[3] * _m[1] );
1089 #endif
1090  return invertible;
1091  }

References _mat, determinant(), moab::Util::is_finite(), and size.

◆ is_positive_definite()

bool moab::Matrix3::is_positive_definite ( )
inline

Definition at line 516 of file Matrix3.hpp.

517  {
518 #ifdef MOAB_HAVE_EIGEN3
519  double subdet6 = _mat( 1 ) * _mat( 5 ) - _mat( 2 ) * _mat( 4 );
520  double subdet7 = _mat( 2 ) * _mat( 3 ) - _mat( 0 ) * _mat( 5 );
521  double subdet8 = _mat( 0 ) * _mat( 4 ) - _mat( 1 ) * _mat( 3 );
522  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
523  const double det = _mat( 6 ) * subdet6 + _mat( 7 ) * subdet7 + _mat( 8 ) * subdet8;
524  return _mat( 0 ) > 0 && subdet8 > 0 && det > 0;
525 #else
526  double subdet6 = _mat[1] * _mat[5] - _mat[2] * _mat[4];
527  double subdet7 = _mat[2] * _mat[3] - _mat[0] * _mat[5];
528  double subdet8 = _mat[0] * _mat[4] - _mat[1] * _mat[3];
529  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
530  const double det = _mat[6] * subdet6 + _mat[7] * subdet7 + _mat[8] * subdet8;
531  return _mat[0] > 0 && subdet8 > 0 && det > 0;
532 #endif
533  }

References _mat.

◆ is_symmetric()

bool moab::Matrix3::is_symmetric ( )
inline

Definition at line 500 of file Matrix3.hpp.

501  {
502  const double EPS = 1e-14;
503 #ifdef MOAB_HAVE_EIGEN3
504  if( ( fabs( _mat( 1 ) - _mat( 3 ) ) < EPS ) && ( fabs( _mat( 2 ) - _mat( 6 ) ) < EPS ) &&
505  ( fabs( _mat( 5 ) - _mat( 7 ) ) < EPS ) )
506  return true;
507 #else
508  if( ( fabs( _mat[1] - _mat[3] ) < EPS ) && ( fabs( _mat[2] - _mat[6] ) < EPS ) &&
509  ( fabs( _mat[5] - _mat[7] ) < EPS ) )
510  return true;
511 #endif
512  else
513  return false;
514  }

References _mat, and EPS.

Referenced by eigen_decomposition().

◆ operator()() [1/4]

double& moab::Matrix3::operator() ( unsigned  i)
inline

Definition at line 397 of file Matrix3.hpp.

398  {
399 #ifdef MOAB_HAVE_EIGEN3
400  return _mat( i );
401 #else
402  return _mat[i];
403 #endif
404  }

References _mat.

◆ operator()() [2/4]

double moab::Matrix3::operator() ( unsigned  i) const
inline

Definition at line 406 of file Matrix3.hpp.

407  {
408 #ifdef MOAB_HAVE_EIGEN3
409  return _mat( i );
410 #else
411  return _mat[i];
412 #endif
413  }

References _mat.

◆ operator()() [3/4]

double& moab::Matrix3::operator() ( unsigned  r,
unsigned  c 
)
inline

Definition at line 379 of file Matrix3.hpp.

380  {
381 #ifdef MOAB_HAVE_EIGEN3
382  return _mat( r, c );
383 #else
384  return _mat[r * 3 + c];
385 #endif
386  }

References _mat.

◆ operator()() [4/4]

double moab::Matrix3::operator() ( unsigned  r,
unsigned  c 
) const
inline

Definition at line 388 of file Matrix3.hpp.

389  {
390 #ifdef MOAB_HAVE_EIGEN3
391  return _mat( r, c );
392 #else
393  return _mat[r * 3 + c];
394 #endif
395  }

References _mat.

◆ operator*=() [1/2]

Matrix3& moab::Matrix3::operator*= ( const Matrix3 &  m)
inline

Definition at line 478 of file Matrix3.hpp.

479  {
480 #ifdef MOAB_HAVE_EIGEN3
481  _mat *= m._mat;
482 #else
483  // Uncomment below if you want point-wise multiplication instead (.*)
484  // for (int i=0; i < Matrix3::size; ++i) _mat[i] *= m._mat[i];
485  std::vector< double > dmat;
486  dmat.assign( _mat, _mat + size );
487  _mat[0] = dmat[0] * m._mat[0] + dmat[1] * m._mat[3] + dmat[2] * m._mat[6];
488  _mat[1] = dmat[0] * m._mat[1] + dmat[1] * m._mat[4] + dmat[2] * m._mat[7];
489  _mat[2] = dmat[0] * m._mat[2] + dmat[1] * m._mat[5] + dmat[2] * m._mat[8];
490  _mat[3] = dmat[3] * m._mat[0] + dmat[4] * m._mat[3] + dmat[5] * m._mat[6];
491  _mat[4] = dmat[3] * m._mat[1] + dmat[4] * m._mat[4] + dmat[5] * m._mat[7];
492  _mat[5] = dmat[3] * m._mat[2] + dmat[4] * m._mat[5] + dmat[5] * m._mat[8];
493  _mat[6] = dmat[6] * m._mat[0] + dmat[7] * m._mat[3] + dmat[8] * m._mat[6];
494  _mat[7] = dmat[6] * m._mat[1] + dmat[7] * m._mat[4] + dmat[8] * m._mat[7];
495  _mat[8] = dmat[6] * m._mat[2] + dmat[7] * m._mat[5] + dmat[8] * m._mat[8];
496 #endif
497  return *this;
498  }

References _mat, and size.

◆ operator*=() [2/2]

Matrix3& moab::Matrix3::operator*= ( double  s)
inline

Definition at line 456 of file Matrix3.hpp.

457  {
458 #ifdef MOAB_HAVE_EIGEN3
459  _mat *= s;
460 #else
461  for( int i = 0; i < Matrix3::size; ++i )
462  _mat[i] *= s;
463 #endif
464  return *this;
465  }

References _mat, and size.

◆ operator+=()

Matrix3& moab::Matrix3::operator+= ( const Matrix3 &  m)
inline

Definition at line 434 of file Matrix3.hpp.

435  {
436 #ifdef MOAB_HAVE_EIGEN3
437  _mat += m._mat;
438 #else
439  for( int i = 0; i < Matrix3::size; ++i )
440  _mat[i] += m._mat[i];
441 #endif
442  return *this;
443  }

References _mat, and size.

◆ operator-=()

Matrix3& moab::Matrix3::operator-= ( const Matrix3 &  m)
inline

Definition at line 445 of file Matrix3.hpp.

446  {
447 #ifdef MOAB_HAVE_EIGEN3
448  _mat -= m._mat;
449 #else
450  for( int i = 0; i < Matrix3::size; ++i )
451  _mat[i] -= m._mat[i];
452 #endif
453  return *this;
454  }

References _mat, and size.

◆ operator/=()

Matrix3& moab::Matrix3::operator/= ( double  s)
inline

Definition at line 467 of file Matrix3.hpp.

468  {
469 #ifdef MOAB_HAVE_EIGEN3
470  _mat /= s;
471 #else
472  for( int i = 0; i < Matrix3::size; ++i )
473  _mat[i] /= s;
474 #endif
475  return *this;
476  }

References _mat, and size.

◆ operator=() [1/2]

Matrix3& moab::Matrix3::operator= ( const double  v[size])
inline

Definition at line 351 of file Matrix3.hpp.

352  {
353 #ifdef MOAB_HAVE_EIGEN3
354  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
355 #else
356  memcpy( _mat, v, size * sizeof( double ) );
357 #endif
358  return *this;
359  }

References _mat, and size.

◆ operator=() [2/2]

Matrix3& moab::Matrix3::operator= ( const Matrix3 &  m)
inline

Definition at line 341 of file Matrix3.hpp.

342  {
343 #ifdef MOAB_HAVE_EIGEN3
344  _mat = m._mat;
345 #else
346  memcpy( _mat, m._mat, size * sizeof( double ) );
347 #endif
348  return *this;
349  }

References _mat, and size.

◆ operator[]() [1/2]

double* moab::Matrix3::operator[] ( unsigned  i)
inline

Definition at line 361 of file Matrix3.hpp.

362  {
363 #ifdef MOAB_HAVE_EIGEN3
364  return _mat.row( i ).data();
365 #else
366  return &_mat[i * 3]; // Row Major
367 #endif
368  }

References _mat.

◆ operator[]() [2/2]

const double* moab::Matrix3::operator[] ( unsigned  i) const
inline

Definition at line 370 of file Matrix3.hpp.

371  {
372 #ifdef MOAB_HAVE_EIGEN3
373  return _mat.row( i ).data();
374 #else
375  return &_mat[i * 3];
376 #endif
377  }

References _mat.

◆ print()

void moab::Matrix3::print ( std::ostream &  s) const
inline

Definition at line 1110 of file Matrix3.hpp.

1111  {
1112 #ifdef MOAB_HAVE_EIGEN3
1113  s << "| " << _mat( 0 ) << " " << _mat( 1 ) << " " << _mat( 2 ) << " | " << _mat( 3 ) << " " << _mat( 4 ) << " "
1114  << _mat( 5 ) << " | " << _mat( 6 ) << " " << _mat( 7 ) << " " << _mat( 8 ) << " |";
1115 #else
1116  s << "| " << _mat[0] << " " << _mat[1] << " " << _mat[2] << " | " << _mat[3] << " " << _mat[4] << " " << _mat[5]
1117  << " | " << _mat[6] << " " << _mat[7] << " " << _mat[8] << " |";
1118 #endif
1119  }

References _mat.

◆ row()

CartVect moab::Matrix3::row ( int  index) const
inline

Definition at line 991 of file Matrix3.hpp.

992  {
993  assert( index < Matrix3::size );
994 #ifdef MOAB_HAVE_EIGEN3
995  Eigen::Vector3d mvec = _mat.row( index );
996  return CartVect( mvec[0], mvec[1], mvec[2] );
997 #else
998  switch( index )
999  {
1000  case 0:
1001  return CartVect( _mat[0], _mat[1], _mat[2] );
1002  case 1:
1003  return CartVect( _mat[3], _mat[4], _mat[5] );
1004  case 2:
1005  return CartVect( _mat[6], _mat[7], _mat[8] );
1006  }
1007  return CartVect( 0.0 );
1008 #endif
1009  }

References _mat, and size.

◆ rowscale()

void moab::Matrix3::rowscale ( int  index,
double  scale 
)
inline

Definition at line 944 of file Matrix3.hpp.

945  {
946  assert( index < Matrix3::size );
947 #ifdef MOAB_HAVE_EIGEN3
948  _mat.row( index ) *= scale;
949 #else
950  switch( index )
951  {
952  case 0:
953  _mat[0] *= scale;
954  _mat[1] *= scale;
955  _mat[2] *= scale;
956  break;
957  case 1:
958  _mat[3] *= scale;
959  _mat[4] *= scale;
960  _mat[5] *= scale;
961  break;
962  case 2:
963  _mat[6] *= scale;
964  _mat[7] *= scale;
965  _mat[8] *= scale;
966  break;
967  }
968 #endif
969  }

References _mat, and size.

◆ subdet()

double moab::Matrix3::subdet ( int  r,
int  c 
) const
inline

Definition at line 1095 of file Matrix3.hpp.

1096  {
1097  assert( r >= 0 && c >= 0 );
1098  if( r < 0 || c < 0 ) return DBL_MAX;
1099 #ifdef MOAB_HAVE_EIGEN3
1100  const int r1 = ( r + 1 ) % 3, r2 = ( r + 2 ) % 3;
1101  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
1102  return _mat( r1, c1 ) * _mat( r2, c2 ) - _mat( r1, c2 ) * _mat( r2, c1 );
1103 #else
1104  const int r1 = Matrix3::size * ( ( r + 1 ) % 3 ), r2 = Matrix3::size * ( ( r + 2 ) % 3 );
1105  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
1106  return _mat[r1 + c1] * _mat[r2 + c2] - _mat[r1 + c2] * _mat[r2 + c1];
1107 #endif
1108  }

References _mat, and size.

◆ swapcol()

void moab::Matrix3::swapcol ( int  srcindex,
int  destindex 
)
inline

Definition at line 849 of file Matrix3.hpp.

850  {
851  assert( srcindex < Matrix3::size );
852  assert( destindex < Matrix3::size );
853 #ifdef MOAB_HAVE_EIGEN3
854  _mat.col( srcindex ).swap( _mat.col( destindex ) );
855 #else
856  CartVect svol = this->vcol< CartVect >( srcindex );
857  CartVect dvol = this->vcol< CartVect >( destindex );
858  switch( srcindex )
859  {
860  case 0:
861  _mat[0] = dvol[0];
862  _mat[3] = dvol[1];
863  _mat[6] = dvol[2];
864  break;
865  case 1:
866  _mat[1] = dvol[0];
867  _mat[4] = dvol[1];
868  _mat[7] = dvol[2];
869  break;
870  case 2:
871  _mat[2] = dvol[0];
872  _mat[5] = dvol[1];
873  _mat[8] = dvol[2];
874  break;
875  }
876  switch( destindex )
877  {
878  case 0:
879  _mat[0] = svol[0];
880  _mat[3] = svol[1];
881  _mat[6] = svol[2];
882  break;
883  case 1:
884  _mat[1] = svol[0];
885  _mat[4] = svol[1];
886  _mat[7] = svol[2];
887  break;
888  case 2:
889  _mat[2] = svol[0];
890  _mat[5] = svol[1];
891  _mat[8] = svol[2];
892  break;
893  }
894 #endif
895  }

References _mat, and size.

Referenced by moab::box_from_axes(), and moab::OrientedBox::order_axes_by_length().

◆ trace()

double moab::Matrix3::trace ( ) const
inline

Definition at line 1020 of file Matrix3.hpp.

1021  {
1022 #ifdef MOAB_HAVE_EIGEN3
1023  return _mat.trace();
1024 #else
1025  return _mat[0] + _mat[4] + _mat[8];
1026 #endif
1027  }

References _mat.

◆ transpose()

Matrix3 moab::Matrix3::transpose ( ) const
inline

Definition at line 806 of file Matrix3.hpp.

807  {
808 #ifdef MOAB_HAVE_EIGEN3
809  return Matrix3( _mat.transpose() );
810 #else
811  Matrix3 mtmp( *this );
812  mtmp._mat[1] = _mat[3];
813  mtmp._mat[3] = _mat[1];
814  mtmp._mat[2] = _mat[6];
815  mtmp._mat[6] = _mat[2];
816  mtmp._mat[5] = _mat[7];
817  mtmp._mat[7] = _mat[5];
818  return mtmp;
819 #endif
820  }

References _mat, and Matrix3().

Referenced by moab::OrientedBox::intersect_ray().

◆ transpose_inplace()

void moab::Matrix3::transpose_inplace ( )
inline

Definition at line 791 of file Matrix3.hpp.

792  {
793 #ifdef MOAB_HAVE_EIGEN3
794  _mat.transposeInPlace();
795 #else
796  Matrix3 mtmp( *this );
797  _mat[1] = mtmp._mat[3];
798  _mat[3] = mtmp._mat[1];
799  _mat[2] = mtmp._mat[6];
800  _mat[6] = mtmp._mat[2];
801  _mat[5] = mtmp._mat[7];
802  _mat[7] = mtmp._mat[5];
803 #endif
804  }

References _mat.

◆ vcol()

template<typename Vector >
Vector moab::Matrix3::vcol ( int  index) const
inline

Definition at line 898 of file Matrix3.hpp.

899  {
900  assert( index < Matrix3::size );
901 #ifdef MOAB_HAVE_EIGEN3
902  return _mat.col( index );
903 #else
904  switch( index )
905  {
906  case 0:
907  return Vector( _mat[0], _mat[3], _mat[6] );
908  case 1:
909  return Vector( _mat[1], _mat[4], _mat[7] );
910  case 2:
911  return Vector( _mat[2], _mat[5], _mat[8] );
912  }
913  return Vector( 0.0 );
914 #endif
915  }

References _mat, and size.

Friends And Related Function Documentation

◆ operator*

Matrix3 operator* ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1154 of file Matrix3.hpp.

1155 {
1156 #ifdef MOAB_HAVE_EIGEN3
1157  return Matrix3( a._mat * b._mat );
1158 #else
1159  return Matrix::mmult3( a, b );
1160 #endif
1161 }

◆ operator+

Matrix3 operator+ ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1130 of file Matrix3.hpp.

1131 {
1132 #ifdef MOAB_HAVE_EIGEN3
1133  return Matrix3( a._mat + b._mat );
1134 #else
1135  Matrix3 s( a );
1136  for( int i = 0; i < Matrix3::size; ++i )
1137  s( i ) += b._mat[i];
1138  return s;
1139 #endif
1140 }

◆ operator-

Matrix3 operator- ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1142 of file Matrix3.hpp.

1143 {
1144 #ifdef MOAB_HAVE_EIGEN3
1145  return Matrix3( a._mat - b._mat );
1146 #else
1147  Matrix3 s( a );
1148  for( int i = 0; i < Matrix3::size; ++i )
1149  s( i ) -= b._mat[i];
1150  return s;
1151 #endif
1152 }

Member Data Documentation

◆ _mat

double moab::Matrix3::_mat[size]
private

Definition at line 182 of file Matrix3.hpp.

Referenced by array(), col(), colscale(), copycol(), copyto(), determinant(), diagonal(), eigen_decomposition(), inverse(), invert(), is_positive_definite(), is_symmetric(), Matrix3(), operator()(), operator*=(), operator+=(), operator-=(), operator/=(), operator=(), operator[](), print(), row(), rowscale(), subdet(), swapcol(), trace(), transpose(), transpose_inplace(), and vcol().

◆ size

const int moab::Matrix3::size = 9
static

Definition at line 176 of file Matrix3.hpp.

Referenced by col(), colscale(), copyto(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), invert(), Matrix3(), operator*=(), operator+=(), operator-=(), operator/=(), operator=(), row(), rowscale(), subdet(), swapcol(), and vcol().

The documentation for this class was generated from the following file: