OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#6106 Apply clang-format on AppFwk

Browse Source
This commit is contained in:
Magne Sjaastad
2020-06-19 07:53:59 +02:00
parent e2ef6a910b
commit bdc536dfc4
407 changed files with 88976 additions and 42309 deletions
Download Patch File Download Diff File Expand all files Collapse all files

110
Fwk/AppFwk/cafTensor/cafTensor3.cpp
View File

@@ -2,39 +2,39 @@
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "cvfBase.h"
#include "cafTensor3.h"
#include "cvfBase.h"
#include "cvfVector3.h"
#include "cvfMatrix3.h"
namespace caf {
namespace caf
{
//--------------------------------------------------------------------------------------------------
/// Compute the cofactors of the 3x3 matrix
/// Cofactor:
/// The determinant obtained by deleting the row and column of a given element and
/// Cofactor:
/// The determinant obtained by deleting the row and column of a given element and
/// preceded by a + or <20> sign depending whether the element is in a + or <20> position as follows:
///
/// + - +
///
/// + - +
/// - + -
/// + - +
//--------------------------------------------------------------------------------------------------
cvf::Mat3d cofactor3(const cvf::Mat3d& mx)
cvf::Mat3d cofactor3( const cvf::Mat3d& mx )
{
int detIdxi[2];
int detIdxj[2];
@@ -42,81 +42,81 @@ cvf::Mat3d cofactor3(const cvf::Mat3d& mx)
cvf::Mat3d cof;
double sign;
for (int i = 0; i < 3; i++)
for ( int i = 0; i < 3; i++ )
{
detIdxi[0] = (i == 0) ? 1 : 0;
detIdxi[1] = (i == 2) ? 1 : 2;
detIdxi[0] = ( i == 0 ) ? 1 : 0;
detIdxi[1] = ( i == 2 ) ? 1 : 2;
for (int j = 0; j < 3; j++)
for ( int j = 0; j < 3; j++ )
{
detIdxj[0] = (j == 0) ? 1 : 0;
detIdxj[1] = (j == 2) ? 1 : 2;
sign = (abs(j - i) == 1) ? -1 : 1;
detIdxj[0] = ( j == 0 ) ? 1 : 0;
detIdxj[1] = ( j == 2 ) ? 1 : 2;
sign = ( abs( j - i ) == 1 ) ? -1 : 1;
cof(i, j) = sign * ( mx(detIdxi[0], detIdxj[0]) * mx(detIdxi[1], detIdxj[1])
- mx(detIdxi[0], detIdxj[1]) * mx(detIdxi[1], detIdxj[0]));
cof( i, j ) = sign * ( mx( detIdxi[0], detIdxj[0] ) * mx( detIdxi[1], detIdxj[1] ) -
mx( detIdxi[0], detIdxj[1] ) * mx( detIdxi[1], detIdxj[0] ) );
}
}
return cof;
}
//--------------------------------------------------------------------------------------------------
/// Compute the eigenvector of the matrix corresponding to the provided eigenvalue
/// Compute the eigenvector of the matrix corresponding to the provided eigenvalue
/// The provided eigenvalue must be an actual eigenvalue of the matrix
//--------------------------------------------------------------------------------------------------
cvf::Vec3d eigenVector3(const cvf::Mat3d& mx, double eigenValue, bool* computedOk)
cvf::Vec3d eigenVector3( const cvf::Mat3d& mx, double eigenValue, bool* computedOk )
{
const double doubleThreshold = 1.0e-60;
if (computedOk) (*computedOk) = false;
if ( computedOk ) ( *computedOk ) = false;
cvf::Mat3d mxMinusEigv = mx;
for (int i = 0; i < 3; i++)
for ( int i = 0; i < 3; i++ )
{
mxMinusEigv(i, i) -= eigenValue;
mxMinusEigv( i, i ) -= eigenValue;
}
cvf::Mat3d cof = cofactor3(mxMinusEigv);
cvf::Mat3d cof = cofactor3( mxMinusEigv );
// Find largest absolute cofactor
int largestCof_i = -1;
int largestCof_j = -1;
double largestCof = 0.0;
for (int i = 0; i < 3; i++)
int largestCof_i = -1;
int largestCof_j = -1;
double largestCof = 0.0;
for ( int i = 0; i < 3; i++ )
{
for (int j = 0; j < 3; j++)
for ( int j = 0; j < 3; j++ )
{
double absCof = fabs(cof(i,j));
double absCof = fabs( cof( i, j ) );
if (absCof > largestCof)
if ( absCof > largestCof )
{
largestCof = absCof;
largestCof = absCof;
largestCof_i = i;
largestCof_j = j;
}
}
}
if (fabs(largestCof) < doubleThreshold) return cvf::Vec3d::ZERO;
if ( fabs( largestCof ) < doubleThreshold ) return cvf::Vec3d::ZERO;
// Find largest matrix element not in the max cofactor row/col
int largestMxElm_i = -1;
int largestMxElm_j = -1;
double largestMxElm = 0.0;
for (int i = 0; i < 3; i++)
int largestMxElm_i = -1;
int largestMxElm_j = -1;
double largestMxElm = 0.0;
for ( int i = 0; i < 3; i++ )
{
if (i != largestCof_i)
if ( i != largestCof_i )
{
for (int j = 0; j < 3; j++)
for ( int j = 0; j < 3; j++ )
{
if (j != largestCof_j)
if ( j != largestCof_j )
{
double absMxElm = fabs(mxMinusEigv(i,j));
double absMxElm = fabs( mxMinusEigv( i, j ) );
if (absMxElm > largestMxElm)
if ( absMxElm > largestMxElm )
{
largestMxElm = absMxElm;
largestMxElm = absMxElm;
largestMxElm_i = i;
largestMxElm_j = j;
}
@@ -126,25 +126,25 @@ cvf::Vec3d eigenVector3(const cvf::Mat3d& mx, double eigenValue, bool* computedO
}
// Check if largest coefficient is zero
if (fabs(largestMxElm) < doubleThreshold) return cvf::Vec3d::ZERO;
if ( fabs( largestMxElm ) < doubleThreshold ) return cvf::Vec3d::ZERO;
// Find last component index
int lastComp_j = 0;
for (int i = 0; i < 3; i++)
for ( int i = 0; i < 3; i++ )
{
if ((i != largestCof_j) && (i != largestMxElm_j)) lastComp_j = i;
if ( ( i != largestCof_j ) && ( i != largestMxElm_j ) ) lastComp_j = i;
}
cvf::Vec3d eigenVector;
eigenVector[largestCof_j] = 1.0;
eigenVector[lastComp_j] = cof(largestCof_i, lastComp_j) / cof(largestCof_i, largestCof_j);
eigenVector[largestMxElm_j] = (-mxMinusEigv(largestMxElm_i, largestCof_j) - mxMinusEigv(largestMxElm_i, lastComp_j)*eigenVector[lastComp_j] )
/ mxMinusEigv(largestMxElm_i, largestMxElm_j);
eigenVector[lastComp_j] = cof( largestCof_i, lastComp_j ) / cof( largestCof_i, largestCof_j );
eigenVector[largestMxElm_j] = ( -mxMinusEigv( largestMxElm_i, largestCof_j ) -
mxMinusEigv( largestMxElm_i, lastComp_j ) * eigenVector[lastComp_j] ) /
mxMinusEigv( largestMxElm_i, largestMxElm_j );
if (computedOk) (*computedOk) = true;
if ( computedOk ) ( *computedOk ) = true;
return eigenVector;
}
}
} // namespace caf

60
Fwk/AppFwk/cafTensor/cafTensor3.h
View File

@@ -2,17 +2,17 @@
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
@@ -24,43 +24,51 @@
namespace caf
{
template< typename S>
template <typename S>
class Tensor3
{
S m_tensor[6]; // SXX, SYY, SZZ, SXY, SYZ, SZX
public:
Tensor3();
Tensor3(S sxx, S syy, S szz, S sxy, S syz, S szx);
Tensor3(const Tensor3& other);
template<typename T>
explicit Tensor3(const Tensor3<T>& other);
Tensor3( S sxx, S syy, S szz, S sxy, S syz, S szx );
Tensor3( const Tensor3& other );
template <typename T>
explicit Tensor3( const Tensor3<T>& other );
static Tensor3 invalid();
static Tensor3 invalid();
inline Tensor3& operator=(const Tensor3& rhs);
inline Tensor3 operator+(const Tensor3& rhs) const;
inline Tensor3 operator*(S scale) const;
inline Tensor3& operator=( const Tensor3& rhs );
inline Tensor3 operator+( const Tensor3& rhs ) const;
inline Tensor3 operator*( S scale ) const;
bool equals(const Tensor3& mat) const;
bool operator==(const Tensor3& rhs) const;
bool operator!=(const Tensor3& rhs) const;
bool equals( const Tensor3& mat ) const;
bool operator==( const Tensor3& rhs ) const;
bool operator!=( const Tensor3& rhs ) const;
enum TensorComponentEnum { SXX, SYY, SZZ, SXY, SYZ, SZX };
inline S& operator[](TensorComponentEnum comp);
inline S operator[](TensorComponentEnum comp) const;
enum TensorComponentEnum
{
SXX,
SYY,
SZZ,
SXY,
SYZ,
SZX
};
inline S& operator[]( TensorComponentEnum comp );
inline S operator[]( TensorComponentEnum comp ) const;
void setFromInternalLayout(S* tensorData);
void setFromAbaqusLayout(S* tensorData);
void setFromInternalLayout( S* tensorData );
void setFromAbaqusLayout( S* tensorData );
cvf::Vec3f calculatePrincipals(cvf::Vec3f principalDirections[3]) const;
float calculateVonMises() const;
cvf::Vec3f calculatePrincipals( cvf::Vec3f principalDirections[3] ) const;
float calculateVonMises() const;
Tensor3 rotated(const cvf::Matrix3<S>& rotMx) const;
Tensor3 rotated( const cvf::Matrix3<S>& rotMx ) const;
};
typedef Tensor3<float> Ten3f;
typedef Tensor3<float> Ten3f;
typedef Tensor3<double> Ten3d;
}
} // namespace caf
#include "cafTensor3.inl"

262
Fwk/AppFwk/cafTensor/cafTensor3.inl
View File

@@ -2,52 +2,51 @@
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "cvfAssert.h"
#include "cvfMath.h"
#include "cvfMatrix3.h"
#include "cvfSystem.h"
#include <algorithm>
#include "cvfMatrix3.h"
#include <cmath>
namespace caf {
namespace caf
{
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
template< typename S>
template <typename S>
caf::Tensor3<S>::Tensor3()
{
m_tensor[0] = (S) 0.0;
m_tensor[1] = (S) 0.0;
m_tensor[2] = (S) 0.0;
m_tensor[3] = (S) 0.0;
m_tensor[4] = (S) 0.0;
m_tensor[5] = (S) 0.0;
m_tensor[0] = (S)0.0;
m_tensor[1] = (S)0.0;
m_tensor[2] = (S)0.0;
m_tensor[3] = (S)0.0;
m_tensor[4] = (S)0.0;
m_tensor[5] = (S)0.0;
}
//----------------------------------------------------------------------------------------------------
/// Copy constructor
//----------------------------------------------------------------------------------------------------
template <typename S>
inline Tensor3<S>::Tensor3(const Tensor3& other)
inline Tensor3<S>::Tensor3( const Tensor3& other )
{
cvf::System::memcpy(m_tensor, sizeof(m_tensor), other.m_tensor, sizeof(other.m_tensor));
cvf::System::memcpy( m_tensor, sizeof( m_tensor ), other.m_tensor, sizeof( other.m_tensor ) );
}
//----------------------------------------------------------------------------------------------------
@@ -55,22 +54,22 @@ inline Tensor3<S>::Tensor3(const Tensor3& other)
//----------------------------------------------------------------------------------------------------
template <typename S>
template <typename T>
Tensor3<S>::Tensor3(const Tensor3<T>& other)
Tensor3<S>::Tensor3( const Tensor3<T>& other )
{
m_tensor[SXX] = other[Tensor3<T>::SXX];
m_tensor[SYY] = other[Tensor3<T>::SYY];
m_tensor[SZZ] = other[Tensor3<T>::SZZ];
m_tensor[SXY] = other[Tensor3<T>::SXY];
m_tensor[SYZ] = other[Tensor3<T>::SYZ];
m_tensor[SZX] = other[Tensor3<T>::SZX];
m_tensor[SZX] = other[Tensor3<T>::SZX];
}
//----------------------------------------------------------------------------------------------------
/// Constructor with explicit initialization of all tensor elements.
///
///
//----------------------------------------------------------------------------------------------------
template <typename S>
Tensor3<S>::Tensor3(S sxx, S syy, S szz, S sxy, S syz, S szx)
Tensor3<S>::Tensor3( S sxx, S syy, S szz, S sxy, S syz, S szx )
{
m_tensor[0] = sxx;
m_tensor[1] = syy;
@@ -86,28 +85,31 @@ Tensor3<S>::Tensor3(S sxx, S syy, S szz, S sxy, S syz, S szx)
template <typename S>
Tensor3<S> caf::Tensor3<S>::invalid()
{
return caf::Tensor3<S>(std::numeric_limits<S>::infinity(), std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(), std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(), std::numeric_limits<S>::infinity());
return caf::Tensor3<S>( std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity(),
std::numeric_limits<S>::infinity() );
}
//----------------------------------------------------------------------------------------------------
/// Assignment operator
//----------------------------------------------------------------------------------------------------
template <typename S>
inline Tensor3<S>& Tensor3<S>::operator=(const Tensor3& obj)
inline Tensor3<S>& Tensor3<S>::operator=( const Tensor3& obj )
{
cvf::System::memcpy(m_tensor, sizeof(m_tensor), obj.m_tensor, sizeof(obj.m_tensor));
cvf::System::memcpy( m_tensor, sizeof( m_tensor ), obj.m_tensor, sizeof( obj.m_tensor ) );
return *this;
}
//--------------------------------------------------------------------------------------------------
/// Component-wise addition
//--------------------------------------------------------------------------------------------------
template< typename S>
Tensor3<S> caf::Tensor3<S>::operator+(const Tensor3& rhs) const
template <typename S>
Tensor3<S> caf::Tensor3<S>::operator+( const Tensor3& rhs ) const
{
Tensor3<S> result(*this);
Tensor3<S> result( *this );
result.m_tensor[0] += rhs.m_tensor[0];
result.m_tensor[1] += rhs.m_tensor[1];
@@ -120,12 +122,12 @@ Tensor3<S> caf::Tensor3<S>::operator+(const Tensor3& rhs) const
}
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
template< typename S>
Tensor3<S> caf::Tensor3<S>::operator*(S scale) const
template <typename S>
Tensor3<S> caf::Tensor3<S>::operator*( S scale ) const
{
Tensor3<S> result(*this);
Tensor3<S> result( *this );
result.m_tensor[0] *= scale;
result.m_tensor[1] *= scale;
@@ -137,42 +139,39 @@ Tensor3<S> caf::Tensor3<S>::operator*(S scale) const
return result;
}
//----------------------------------------------------------------------------------------------------
/// Check if matrices are equal using exact comparisons.
//----------------------------------------------------------------------------------------------------
template<typename S>
bool Tensor3<S>::equals(const Tensor3& ten) const
template <typename S>
bool Tensor3<S>::equals( const Tensor3& ten ) const
{
for (int i = 0; i < 6; i++)
for ( int i = 0; i < 6; i++ )
{
if (m_tensor[i] != ten.m_tensor[i]) return false;
if ( m_tensor[i] != ten.m_tensor[i] ) return false;
}
return true;
}
//----------------------------------------------------------------------------------------------------
/// Comparison operator. Checks for equality using exact comparisons.
//----------------------------------------------------------------------------------------------------
template <typename S>
bool Tensor3<S>::operator==(const Tensor3& rhs) const
bool Tensor3<S>::operator==( const Tensor3& rhs ) const
{
return this->equals(rhs);
return this->equals( rhs );
}
//----------------------------------------------------------------------------------------------------
/// Comparison operator. Checks for not equal using exact comparisons.
//----------------------------------------------------------------------------------------------------
template <typename S>
bool Tensor3<S>::operator!=(const Tensor3& rhs) const
bool Tensor3<S>::operator!=( const Tensor3& rhs ) const
{
int i;
for (i = 0; i < 6; i++)
for ( i = 0; i < 6; i++ )
{
if (m_tensor[i] != rhs.m_tensor[i]) return true;
if ( m_tensor[i] != rhs.m_tensor[i] ) return true;
}
return false;
@@ -181,34 +180,32 @@ bool Tensor3<S>::operator!=(const Tensor3& rhs) const
//--------------------------------------------------------------------------------------------------
/// Get modifiable component 0,1,2. E.g. x = v[0];
//--------------------------------------------------------------------------------------------------
template<typename S>
inline S Tensor3<S>::operator[](TensorComponentEnum index) const
template <typename S>
inline S Tensor3<S>::operator[]( TensorComponentEnum index ) const
{
CVF_TIGHT_ASSERT(index >= 0);
CVF_TIGHT_ASSERT(index < 6);
CVF_TIGHT_ASSERT( index >= 0 );
CVF_TIGHT_ASSERT( index < 6 );
return m_tensor[index];
}
//--------------------------------------------------------------------------------------------------
/// Get const component 0,1,2. E.g. x = v[0];
//--------------------------------------------------------------------------------------------------
template<typename S>
inline S& Tensor3<S>::operator[](TensorComponentEnum index)
template <typename S>
inline S& Tensor3<S>::operator[]( TensorComponentEnum index )
{
CVF_TIGHT_ASSERT(index >= 0);
CVF_TIGHT_ASSERT(index < 6);
CVF_TIGHT_ASSERT( index >= 0 );
CVF_TIGHT_ASSERT( index < 6 );
return m_tensor[index];
}
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
template< typename S>
void Tensor3<S>::setFromInternalLayout(S* tensorData)
template <typename S>
void Tensor3<S>::setFromInternalLayout( S* tensorData )
{
m_tensor[0] = tensorData[0];
m_tensor[1] = tensorData[1];
@@ -218,12 +215,11 @@ void Tensor3<S>::setFromInternalLayout(S* tensorData)
m_tensor[5] = tensorData[5];
}
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
template< typename S>
void Tensor3<S>::setFromAbaqusLayout(S* tensorData)
template <typename S>
void Tensor3<S>::setFromAbaqusLayout( S* tensorData )
{
m_tensor[0] = tensorData[0];
m_tensor[1] = tensorData[1];
@@ -233,32 +229,30 @@ void Tensor3<S>::setFromAbaqusLayout(S* tensorData)
m_tensor[5] = tensorData[4];
}
cvf::Mat3d cofactor3(const cvf::Mat3d& mx);
cvf::Vec3d eigenVector3(const cvf::Mat3d& mx, double eigenValue, bool* computedOk);
cvf::Mat3d cofactor3( const cvf::Mat3d& mx );
cvf::Vec3d eigenVector3( const cvf::Mat3d& mx, double eigenValue, bool* computedOk );
//--------------------------------------------------------------------------------------------------
/// Compute principal values and optionally the principal directions
/// The tensor must be laid out as follows: SXX, SYY, SZZ, SXY, SYZ, SZX
//--------------------------------------------------------------------------------------------------
template<typename S>
cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) const
template <typename S>
cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3] ) const
{
CVF_TIGHT_ASSERT(m_tensor);
CVF_TIGHT_ASSERT( m_tensor );
const float floatThreshold = 1.0e-30f;
const float floatThreshold = 1.0e-30f;
const double doubleThreshold = 1.0e-60;
cvf::Vec3f principalValues;
// Init return arrays to invalid
// Init return arrays to invalid
principalValues[0] = std::numeric_limits<float>::infinity();
principalValues[1] = std::numeric_limits<float>::infinity();
principalValues[2] = std::numeric_limits<float>::infinity();
if (principalDirections)
if ( principalDirections )
{
principalDirections[0] = cvf::Vec3f::ZERO;
principalDirections[1] = cvf::Vec3f::ZERO;
@@ -266,29 +260,29 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
}
// Return if we have an undefined component
int i;
for (i = 0; i < 6; i++)
for ( i = 0; i < 6; i++ )
{
if (m_tensor[i] == std::numeric_limits<S>::infinity())
if ( m_tensor[i] == std::numeric_limits<S>::infinity() )
{
return principalValues;
}
}
// Return 0, 0, 0 if all components are zero
bool isAllTensCompsZero = true;
for (i = 0; i < 6; i++)
for ( i = 0; i < 6; i++ )
{
if (!(fabs(m_tensor[i]) < floatThreshold))
if ( !( fabs( m_tensor[i] ) < floatThreshold ) )
{
isAllTensCompsZero = false;
break;
}
}
if (isAllTensCompsZero)
if ( isAllTensCompsZero )
{
return cvf::Vec3f::ZERO;
}
@@ -296,7 +290,7 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
double SXX = m_tensor[0], SYY = m_tensor[1], SZZ = m_tensor[2];
double SXY = m_tensor[3], SYZ = m_tensor[4], SZX = m_tensor[5];
double pressure = -(SXX + SYY + SZZ)/3.0;
double pressure = -( SXX + SYY + SZZ ) / 3.0;
// Normally we would solve the eigenvalues by solving the 3'rd degree equation:
// -sigma^3 + A*sigma^2 - B*sigma + C = 0
@@ -314,15 +308,11 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
double S1, S2, S3;
double AA, BB, CC, DD, angleP;
AA = SXY*SXY + SYZ*SYZ + SZX*SZX - SXX*SYY - SYY*SZZ - SXX*SZZ;
AA = SXY * SXY + SYZ * SYZ + SZX * SZX - SXX * SYY - SYY * SZZ - SXX * SZZ;
BB = SXX * SYZ * SYZ
+ SYY * SZX * SZX
+ SZZ * SXY * SXY
- SXX * SYY * SZZ
- 2.0 * SXY * SYZ * SZX;
BB = SXX * SYZ * SYZ + SYY * SZX * SZX + SZZ * SXY * SXY - SXX * SYY * SZZ - 2.0 * SXY * SYZ * SZX;
if (fabs(AA) < doubleThreshold)
if ( fabs( AA ) < doubleThreshold )
{
S1 = 0.0;
S2 = 0.0;
@@ -330,16 +320,18 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
}
else
{
CC = -sqrt(27.0/AA) * BB * 0.5 / AA;
CC = -sqrt( 27.0 / AA ) * BB * 0.5 / AA;
if (CC > 1.0) CC = 1.0;
else if (CC < -1.0) CC = -1.0;
if ( CC > 1.0 )
CC = 1.0;
else if ( CC < -1.0 )
CC = -1.0;
angleP = acos(CC)/3.0;
DD = 2.0*sqrt(AA/3.0);
S1 = DD*cos(angleP);
S2 = DD*cos(angleP + 4.0*cvf::PI_D/3.0);
S3 = DD*cos(angleP + 2.0*cvf::PI_D/3.0);
angleP = acos( CC ) / 3.0;
DD = 2.0 * sqrt( AA / 3.0 );
S1 = DD * cos( angleP );
S2 = DD * cos( angleP + 4.0 * cvf::PI_D / 3.0 );
S3 = DD * cos( angleP + 2.0 * cvf::PI_D / 3.0 );
}
int idxPMin = 2;
@@ -347,37 +339,43 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
int idxPMax = 0;
double principalsd[3];
principalsd[idxPMax] = (S1 - pressure);
principalsd[idxPMid] = (S2 - pressure);
principalsd[idxPMin] = (S3 - pressure);
principalsd[idxPMax] = ( S1 - pressure );
principalsd[idxPMid] = ( S2 - pressure );
principalsd[idxPMin] = ( S3 - pressure );
// Sort the principals if we have no Z component in the tensor at all
if ((m_tensor[2] == 0.0f) && (m_tensor[4] == 0.0f) && (m_tensor[5] == 0.0f))
// Sort the principals if we have no Z component in the tensor at all
if ( ( m_tensor[2] == 0.0f ) && ( m_tensor[4] == 0.0f ) && ( m_tensor[5] == 0.0f ) )
{
if (fabs(principalsd[idxPMin]) > fabs(principalsd[idxPMid])) std::swap(idxPMin, idxPMid);
if (fabs(principalsd[idxPMin]) > fabs(principalsd[idxPMax])) std::swap(idxPMin, idxPMax);
if (principalsd[idxPMax] < principalsd[idxPMid]) std::swap(idxPMax, idxPMid);
if ( fabs( principalsd[idxPMin] ) > fabs( principalsd[idxPMid] ) ) std::swap( idxPMin, idxPMid );
if ( fabs( principalsd[idxPMin] ) > fabs( principalsd[idxPMax] ) ) std::swap( idxPMin, idxPMax );
if ( principalsd[idxPMax] < principalsd[idxPMid] ) std::swap( idxPMax, idxPMid );
principalsd[idxPMin] = 0;
}
// Calculate the principal directions if needed
if (principalDirections)
if ( principalDirections )
{
cvf::Mat3d T;
T(0,0) = m_tensor[0]; T(0,1) = m_tensor[3]; T(0,2) = m_tensor[5];
T(1,0) = m_tensor[3]; T(1,1) = m_tensor[1]; T(1,2) = m_tensor[4];
T(2,0) = m_tensor[5]; T(2,1) = m_tensor[4]; T(2,2) = m_tensor[2];
T( 0, 0 ) = m_tensor[0];
T( 0, 1 ) = m_tensor[3];
T( 0, 2 ) = m_tensor[5];
T( 1, 0 ) = m_tensor[3];
T( 1, 1 ) = m_tensor[1];
T( 1, 2 ) = m_tensor[4];
T( 2, 0 ) = m_tensor[5];
T( 2, 1 ) = m_tensor[4];
T( 2, 2 ) = m_tensor[2];
principalDirections[0] = cvf::Vec3f(eigenVector3(T, principalsd[idxPMax], NULL));
principalDirections[0] = cvf::Vec3f( eigenVector3( T, principalsd[idxPMax], NULL ) );
principalDirections[0].normalize();
principalDirections[1] = cvf::Vec3f(eigenVector3(T, principalsd[idxPMid], NULL));
principalDirections[1] = cvf::Vec3f( eigenVector3( T, principalsd[idxPMid], NULL ) );
principalDirections[1].normalize();
principalDirections[2] = cvf::Vec3f(eigenVector3(T, principalsd[idxPMin], NULL));
principalDirections[2] = cvf::Vec3f( eigenVector3( T, principalsd[idxPMin], NULL ) );
principalDirections[2].normalize();
}
principalValues[0] = (float)principalsd[idxPMax];
principalValues[1] = (float)principalsd[idxPMid];
principalValues[2] = (float)principalsd[idxPMin];
@@ -386,34 +384,42 @@ cvf::Vec3f Tensor3<S>::calculatePrincipals( cvf::Vec3f principalDirections[3]) c
}
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
template< typename S>
template <typename S>
float caf::Tensor3<S>::calculateVonMises() const
{
return (float) sqrt( ( (m_tensor[0]*m_tensor[0] + m_tensor[1]*m_tensor[1] + m_tensor[2]*m_tensor[2]) ) +
( -(m_tensor[0]*m_tensor[1] + m_tensor[1]*m_tensor[2] + m_tensor[0]*m_tensor[2]) ) +
( 3*(m_tensor[3]*m_tensor[3] + m_tensor[4]*m_tensor[4] + m_tensor[5]*m_tensor[5]) ) );
return (float)sqrt( ( ( m_tensor[0] * m_tensor[0] + m_tensor[1] * m_tensor[1] + m_tensor[2] * m_tensor[2] ) ) +
( -( m_tensor[0] * m_tensor[1] + m_tensor[1] * m_tensor[2] + m_tensor[0] * m_tensor[2] ) ) +
( 3 * ( m_tensor[3] * m_tensor[3] + m_tensor[4] * m_tensor[4] + m_tensor[5] * m_tensor[5] ) ) );
}
//--------------------------------------------------------------------------------------------------
/// Calculates Trot = rotMx*T*transpose(rotMx)
//--------------------------------------------------------------------------------------------------
template< typename S>
Tensor3<S> caf::Tensor3<S>::rotated(const cvf::Matrix3<S>& rotMx) const
template <typename S>
Tensor3<S> caf::Tensor3<S>::rotated( const cvf::Matrix3<S>& rotMx ) const
{
cvf::Matrix3<S> tensor(m_tensor[SXX], m_tensor[SXY], m_tensor[SZX],
m_tensor[SXY], m_tensor[SYY], m_tensor[SYZ],
m_tensor[SZX], m_tensor[SYZ], m_tensor[SZZ]);
cvf::Matrix3<S> tensor( m_tensor[SXX],
m_tensor[SXY],
m_tensor[SZX],
m_tensor[SXY],
m_tensor[SYY],
m_tensor[SYZ],
m_tensor[SZX],
m_tensor[SYZ],
m_tensor[SZZ] );
cvf::Matrix3<S> transposedRotMx = rotMx;
transposedRotMx.transpose();
cvf::Matrix3<S> rotatedTensor = rotMx * tensor * transposedRotMx;
return Tensor3(rotatedTensor(0,0), rotatedTensor(1,1), rotatedTensor(2,2),
rotatedTensor(1,0), rotatedTensor(1,2), rotatedTensor(0,2));
}
cvf::Matrix3<S> rotatedTensor = rotMx * tensor * transposedRotMx;
return Tensor3( rotatedTensor( 0, 0 ),
rotatedTensor( 1, 1 ),
rotatedTensor( 2, 2 ),
rotatedTensor( 1, 0 ),
rotatedTensor( 1, 2 ),
rotatedTensor( 0, 2 ) );
}
} // namespace caf