OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed: Made the Vec3Oper unit tests a little bit more challenging.

Browse Source 
Changed: The global *-operator for symmetric tensors is renamed to ddot,
while adding a new *-operator multiplying two symmetric tensors.
Added: New method Tensor::outerProd and an overloaded method Tensor::principal.
Added: Missing default: in switch statements to satisfy picky compilers.
Added: class Tensor4 and a default constructor for SymmTensor4.
Added: Some math utilities in the utl namespace.
This commit is contained in:
Knut Morten Okstad
2015-10-22 17:04:33 +02:00
parent 5c54db94eb
commit d8a4c5501d
11 changed files with 642 additions and 205 deletions
Download Patch File Download Diff File Expand all files Collapse all files

65
src/Utility/Math.C Normal file
View File

@@ -0,0 +1,65 @@
// $Id$
//==============================================================================
//!
//! \file Math.C
//!
//! \date Oct 20 2015
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Various math utility methods.
//!
//==============================================================================
#include "Math.h"
size_t utl::Pascal (int p, unsigned short int nsd)
{
size_t nM = 1;
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
if (nsd == 2)
nM++;
else for (int j = q; j >= 0; j--)
if (i+j <= q) nM++;
return nM;
}
void utl::Pascal (int p, Real x, Real y, std::vector<Real>& phi)
{
phi.clear();
phi.reserve(Pascal(p,2));
phi.push_back(Real(1));
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
{
int k, j = q-i;
Real a = Real(1);
for (k = 0; k < i; k++) a *= x;
for (k = 0; k < j; k++) a *= y;
phi.push_back(a);
}
}
void utl::Pascal (int p, Real x, Real y, Real z, std::vector<Real>& phi)
{
phi.clear();
phi.reserve(Pascal(p,3));
phi.push_back(Real(1));
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
for (int j = q; j >= 0; j--)
if (i+j <= q)
{
int l, k = q-i-j;
Real a = Real(1);
for (l = 0; l < i; l++) a *= x;
for (l = 0; l < j; l++) a *= y;
for (l = 0; l < k; l++) a *= z;
phi.push_back(a);
}
}

41
src/Utility/Math.h Normal file
View File

@@ -0,0 +1,41 @@
// $Id$
//==============================================================================
//!
//! \file Math.h
//!
//! \date Oct 20 2015
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Various math utility methods.
//!
//==============================================================================
#ifndef UTL_MATH_H
#define UTL_MATH_H
#include <cmath>
#include <cstddef>
#include <vector>
namespace utl
{
//! \brief Returns the number of monomials in Pascal's triangle.
//! \param[in] p Polynomial order (>= 0)
//! \param[in] nsd Number of spatial dimensions (2 or 3)
size_t Pascal(int p, unsigned short int nsd);
//! \brief Evaluates the monomials of Pascal's triangle in 2D for order \a p.
void Pascal(int p, Real x, Real y, std::vector<Real>& phi);
//! \brief Evaluates the monomials of Pascal's triangle in 3D for order \a p.
void Pascal(int p, Real x, Real y, Real z, std::vector<Real>& phi);
//! \brief Evaluates the positive ramp function \f$(x+|x|)/2\f$.
inline Real Pos(Real x) { return x > Real(0) ? x : Real(0); }
//! \brief Evaluates the negative ramp function \f$(x-|x|)/2\f$.
inline Real Neg(Real x) { return x < Real(0) ? x : Real(0); }
//! \brief Converts from degrees to radians.
inline Real Rad(Real x) { return x*M_PI/Real(180); }
}
#endif

164
src/Utility/Tensor.C
View File

@@ -13,6 +13,7 @@
#include "Tensor.h"
#include "Vec3.h"
#include <array>
#include <algorithm>
#include <cstring>
@@ -49,9 +50,9 @@ std::ostream& Tensor::print (std::ostream& os) const
return os << v[0] <<' '<< v[3] <<' '<< v[6] <<'\n'
<< v[1] <<' '<< v[4] <<' '<< v[7] <<'\n'
<< v[2] <<' '<< v[5] <<' '<< v[8] << std::endl;
default:
return os;
}
return os;
}
@@ -225,16 +226,21 @@ Tensor::Tensor (Real alpha, t_ind axis) : n(3)
v[5] = sa;
v[7] = -sa;
break;
case 2:
v[0] = v[8] = ca;
v[2] = -sa;
v[6] = -sa;
break;
case 3:
v[0] = v[4] = ca;
v[1] = sa;
v[3] = -sa;
break;
default:
break;
}
}
@@ -405,22 +411,27 @@ Tensor& Tensor::postMult (const Tensor& B)
case 1:
v[0] *= B.v[0];
break;
case 2:
{
Tensor A(*this);
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j);
break;
}
break;
case 3:
{
Tensor A(*this);
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j);
break;
}
break;
default:
break;
}
return *this;
@@ -433,28 +444,43 @@ Tensor& Tensor::preMult (const Tensor& A)
case 1:
v[0] *= A.v[0];
break;
case 2:
{
Tensor B(*this);
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j);
break;
}
break;
case 3:
{
Tensor B(*this);
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
v[this->index(i,j)] = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j);
break;
}
break;
default:
break;
}
return *this;
}
Tensor& Tensor::outerProd (const Vec3& a, const Vec3& b)
{
for (t_ind i = 1; i <= n; i++)
for (t_ind j = 1; j <= n; j++)
v[this->index(i,j)] = a[i-1]*b[j-1];
return *this;
}
Real Tensor::innerProd (const Tensor& T) const
{
Real value = Real(0);
@@ -493,10 +519,15 @@ Tensor& Tensor::transpose ()
case 2:
std::swap(v[1],v[2]);
break;
case 3:
std::swap(v[1],v[3]);
std::swap(v[2],v[6]);
std::swap(v[5],v[7]);
break;
default:
break;
}
return *this;
@@ -509,10 +540,15 @@ Tensor& Tensor::symmetrize ()
case 2:
v[1] = v[2] = 0.5*(v[1]+v[2]);
break;
case 3:
v[1] = v[3] = 0.5*(v[1]+v[3]);
v[2] = v[6] = 0.5*(v[2]+v[6]);
v[5] = v[7] = 0.5*(v[5]+v[7]);
break;
default:
break;
}
return *this;
@@ -627,8 +663,9 @@ Vec3 operator* (const Tensor& T, const Vec3& v)
return Vec3(T(1,1)*v.x + T(1,2)*v.y + T(1,3)*v.z,
T(2,1)*v.x + T(2,2)*v.y + T(2,3)*v.z,
T(3,1)*v.x + T(3,2)*v.y + T(3,3)*v.z);
default:
return v;
}
return v;
}
@@ -648,8 +685,9 @@ Vec3 operator* (const Vec3& v, const Tensor& T)
return Vec3(T(1,1)*v.x + T(2,1)*v.y + T(3,1)*v.z,
T(1,2)*v.x + T(2,2)*v.y + T(3,2)*v.z,
T(1,3)*v.x + T(2,3)*v.y + T(3,3)*v.z);
default:
return v;
}
return v;
}
@@ -665,16 +703,21 @@ Tensor operator* (const Tensor& A, const Tensor& B)
case 1:
C.v[0] = A.v[0]*B.v[0];
break;
case 2:
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j);
break;
case 3:
for (int i = 1; i <= 3; i++)
for (int j = 1; j <= 3; j++)
C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j);
break;
default:
break;
}
return C;
@@ -694,9 +737,9 @@ std::ostream& SymmTensor::print (std::ostream& os) const
return os << v[0] <<'\n'
<< v[3] <<' '<< v[1] <<'\n'
<< v[5] <<' '<< v[4] <<' '<< v[2] << std::endl;
default:
return os;
}
return os;
}
@@ -812,6 +855,10 @@ SymmTensor& SymmTensor::transform (const Tensor& T)
v[4] = S23;
v[5] = S31*T(1,1) + S32*T(1,2);
}
break;
default:
break;
}
return *this;
@@ -918,6 +965,40 @@ SymmTensor& SymmTensor::rightCauchyGreen (const Tensor& F)
v[3] = F(1,1)*F(1,2) + F(2,1)*F(2,2) + F(3,1)*F(3,2);
v[4] = F(1,2)*F(1,3) + F(2,2)*F(2,3) + F(3,2)*F(3,3);
v[5] = F(1,1)*F(1,3) + F(2,1)*F(2,3) + F(3,1)*F(3,3);
break;
default:
break;
}
return *this;
}
SymmTensor& SymmTensor::outerProd (const Vec3& u)
{
switch (n) {
case 1:
v[0] = u.x*u.x;
break;
case 2:
v[0] = u.x*u.x;
v[1] = u.y*u.y;
v[2] = u.x*u.y;
break;
case 3:
v[0] = u.x*u.x;
v[1] = u.y*u.y;
v[2] = u.z*u.z;
v[3] = u.x*u.y;
v[4] = u.y*u.z;
v[5] = u.x*u.z;
break;
default:
break;
}
return *this;
@@ -995,7 +1076,7 @@ bool SymmTensor::principal (Vec3& p) const
p.z = Real(0);
return true;
}
}
// Compute mean and deviatoric (upper triangular part) tensors
const Real tol(1.0e-12);
@@ -1142,6 +1223,19 @@ bool SymmTensor::principal (Vec3& p, Vec3* pdir, int ndir) const
}
bool SymmTensor::principal (Vec3& p, SymmTensor* M) const
{
std::array<Vec3,3> pdir;
if (!this->principal(p,pdir.data()))
return false;
for (t_ind a = 0; a < n; a++)
M[a].outerProd(pdir[a]);
return true;
}
/*!
\brief Adding a scaled unit tensor to a symmetric tensor.
*/
@@ -1193,33 +1287,39 @@ SymmTensor operator* (Real a, const SymmTensor& T)
}
SymmTensor4::SymmTensor4 (const std::vector<Real>& x, t_ind nsd)
: n(nsd), m(0), v(x)
/*!
\brief Multiplication between two symmetric tensors.
*/
SymmTensor operator* (const SymmTensor& A, const SymmTensor& B)
{
if (n == 3)
m = 6;
else if (n == 2)
m = 3;
else
std::cerr <<" *** Invalid fourth-order tensor, dim="<< n << std::endl;
SymmTensor C(std::min(A.n,B.n));
if (v.size() < (size_t)m*m)
std::cerr <<" *** Invalid fourth-order tensor,"
<<" matrix represention too small, size="<< v.size() << std::endl;
switch (C.n) {
case 1:
C.v[0] = A.v[0]*B.v[0];
break;
ptr = (Real*)&v.front();
}
case 2:
C.v[0] = A.v[0]*B.v[0] + A.v[2]*B.v[2];
C.v[1] = A.v[2]*B.v[2] + A.v[1]*B.v[1];
C.v[2] = A.v[0]*B.v[2] + A.v[2]*B.v[1];
break;
case 3:
C.v[0] = A.v[0]*B.v[0] + A.v[3]*B.v[3] + A.v[5]*B.v[5];
C.v[1] = A.v[3]*B.v[3] + A.v[1]*B.v[1] + A.v[4]*B.v[4];
C.v[2] = A.v[5]*B.v[5] + A.v[4]*B.v[4] + A.v[2]*B.v[2];
C.v[3] = A.v[0]*B.v[3] + A.v[3]*B.v[1] + A.v[5]*B.v[4];
C.v[4] = A.v[3]*B.v[5] + A.v[1]*B.v[4] + A.v[4]*B.v[2];
C.v[5] = A.v[0]*B.v[5] + A.v[3]*B.v[4] + A.v[5]*B.v[2];
break;
const Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l) const
{
return v[this->index(i,j)*m+this->index(k,l)];
}
default:
break;
}
Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l)
{
return ptr[this->index(i,j)*m+this->index(k,l)];
return C;
}

54
src/Utility/Tensor.h
View File

@@ -106,6 +106,9 @@ public:
//! \brief Pre-multiplication with another Tensor.
Tensor& preMult(const Tensor& A);
//! \brief Dyadic (outer) product between two vectors.
Tensor& outerProd(const Vec3& a, const Vec3& b);
//! \brief Returns the inner-product of \a *this and the given tensor.
Real innerProd(const Tensor& T) const;
@@ -236,6 +239,9 @@ public:
//! \brief Constructs the right Cauchy-Green tensor from a deformation tensor.
SymmTensor& rightCauchyGreen(const Tensor& F);
//! \brief Dyadic (outer) product between two identical vectors.
SymmTensor& outerProd(const Vec3& u);
//! \brief Returns the inner-product (L2-norm) of the symmetric tensor.
Real L2norm(bool doSqrt = true) const;
//! \brief Returns the von Mises value of the symmetric tensor.
@@ -244,6 +250,8 @@ public:
bool principal(Vec3& p) const;
//! \brief Computes the principal values and associated principal directions.
bool principal(Vec3& p, Vec3* pdir, int ndir = 0) const;
//! \brief Computes the principal values and associated principal directions.
bool principal(Vec3& p, SymmTensor* M) const;
// Global operators
@@ -254,6 +262,8 @@ public:
//! \brief Multiplication between a scalar and a symmetric tensor.
friend SymmTensor operator*(Real a, const SymmTensor& T);
//! \brief Multiplication between two symmetric tensors.
friend SymmTensor operator*(const SymmTensor& A, const SymmTensor& B);
};
@@ -269,53 +279,13 @@ inline SymmTensor operator-(const SymmTensor& A, const SymmTensor& B)
SymmTensor C(A); C -= B; return C;
}
//! \brief Inner-product of two symmetric tensors.
inline Real operator*(const SymmTensor& A, const SymmTensor& B)
//! \brief Inner-product (:-operator) of two symmetric tensors.
inline Real ddot(const SymmTensor& A, const SymmTensor& B)
{
return A.innerProd(B);
}
/*!
\brief Simple class for representing a symmetric fourth-order tensor.
*/
class SymmTensor4
{
typedef unsigned short int t_ind; //!< Tensor index type
t_ind n; //!< Number of spatial dimensions for the tensor
t_ind m; //!< Dimension of the matrix representation
const std::vector<Real>& v; //!< The actual tensor component values
Real* ptr; //!< Non-const pointer to tensor component values
//! \brief Returns a 0-based array index for the given row and column indices.
//! \details Symmetric second-order tensors in 3D are assumed stored with the
//! following order in a one-dimensional array: s11, s22, s33, s12, s23, s13.
inline t_ind index(t_ind i, t_ind j) const
{
if (i == j)
return i-1; // diagonal term
else if (n == 2)
return 2; // off-diagonal term (2D)
if (i == j+1 || i+2 == j) std::swap(i,j);
return i+2; // upper triangular term (3D)
}
public:
//! \brief The constructor creates a tensor from a vector of components.
//! \details The provided vector is assumed to contain the components of the
//! matrix representation of the tensor, stored in a one-dimensional array.
SymmTensor4(const std::vector<Real>& x, t_ind nsd = 3);
//! \brief Index-1 based component reference.
const Real& operator()(t_ind i, t_ind j, t_ind d, t_ind l) const;
//! \brief Index-1 based component access.
Real& operator()(t_ind i, t_ind j, t_ind k, t_ind l);
};
/*!
\brief Abstract interface to problem-specific local coordinate systems.
*/

151
src/Utility/Tensor4.C Normal file
View File

@@ -0,0 +1,151 @@
// $Id$
//==============================================================================
//!
//! \file Tensor4.C
//!
//! \date Oct 27 2015
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of fourth-order tensors with some basic operations.
//!
//==============================================================================
#include "Tensor4.h"
#include "matrix.h"
Tensor4::Tensor4 (t_ind nd, Real scale, bool makeJ)
{
this->redim(nd);
v.resize(m*m,Real(0));
if (scale != Real(0))
{
if (makeJ) // Define the "identity" tensor, J
for (t_ind i = 1; i <= n; i++)
for (t_ind j = 1; j <= n; j++)
v[this->index(i,i,j,j)] = scale;
else // Define the identity tensor, I
for (t_ind i = 1; i <= n; i++)
v[this->index(i,i,i,i)] = scale;
}
}
Tensor4::Tensor4 (const std::vector<Real>& x, t_ind nd) : v(x)
{
this->redim(nd);
if (v.size() < (size_t)m*m)
std::cerr <<" *** Invalid fourth-order tensor,"
<<" matrix represention too small, size="<< v.size() << std::endl;
else if (v.size() > (size_t)m*m)
v.resize(m*m);
}
std::ostream& Tensor4::print (std::ostream& os) const
{
utl::matrix<Real> A(m,m);
A.fill(&v.front());
return os << A;;
}
const Real& Tensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l) const
{
return v[this->index(i,j,k,l)];
}
Real& Tensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l)
{
return v[this->index(i,j,k,l)];
}
Tensor4& Tensor4::operator= (const Tensor4& T)
{
this->redim(T.n);
v.resize(m*m,Real(0));
std::copy(T.v.begin(),T.v.end(),v.begin());
return *this;
}
Tensor4& Tensor4::operator= (Real val)
{
this->zero();
for (t_ind i = 1; i <= n; i++)
v[this->index(i,i,i,i)] = val;
return *this;
}
Tensor4& Tensor4::operator+= (Real val)
{
for (t_ind i = 1; i <= n; i++)
v[this->index(i,i,i,i)] += val;
return *this;
}
SymmTensor4::SymmTensor4 (t_ind nd, bool makeJ) : Tensor4(nd,Real(0))
{
if (makeJ) // Define the "identity" tensor, J
for (t_ind i = 0; i < n; i++)
for (t_ind j = 0; j < n; j++)
v[i*m+j] = Real(1);
else // Define the identity tensor, I
for (t_ind i = 0; i < n; i++)
v[i*m+i] = Real(1);
}
SymmTensor4::SymmTensor4 (const std::vector<Real>& x, t_ind nd) : Tensor4(x,nd)
{
}
std::ostream& SymmTensor4::print (std::ostream& os) const
{
utl::matrix<Real> A(m,m);
for (t_ind l = 1; l <= n; l++)
for (t_ind k = 1; k <= n; k++)
for (t_ind j = 1; j <= n; j++)
for (t_ind i = 1; i <= n; i++)
A(i+n*(j-1),k+l*(n-1)) = v[this->index(i,j)*m+this->index(k,l)];
return os << A;;
}
void SymmTensor4::redim (t_ind ndim)
{
switch (n = ndim) {
case 1: m = 1; break;
case 2: m = 3; break;
case 3: m = 6; break;
default: m = 0;
std::cerr <<" *** Invalid fourth-order tensor, dim="<< n << std::endl;
}
}
const Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l) const
{
return v[this->index(i,j)*m+this->index(k,l)];
}
Real& SymmTensor4::operator() (t_ind i, t_ind j, t_ind k, t_ind l)
{
return v[this->index(i,j)*m+this->index(k,l)];
}

131
src/Utility/Tensor4.h Normal file
View File

@@ -0,0 +1,131 @@
// $Id$
//==============================================================================
//!
//! \file Tensor4.h
//!
//! \date Oct 27 2015
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of fourth-order tensors with some basic operations.
//!
//==============================================================================
#ifndef _TENSOR4_H
#define _TENSOR4_H
#include <vector>
#include <iostream>
/*!
\brief Simple class for representing a non-symmetric fourth-order tensor.
*/
class Tensor4
{
protected:
typedef unsigned short int t_ind; //!< Tensor index type (for convenience)
t_ind n; //!< Number of spatial dimensions for the tensor
t_ind m; //!< Dimension of the matrix representation
std::vector<Real> v; //!< The actual tensor component values
//! \brief Auxilliary method used by the constructors.
virtual void redim(t_ind nsd) { n = nsd; m = n*n; }
//! \brief Prints out the tensor to an output stream.
virtual std::ostream& print(std::ostream& os) const;
private:
//! \brief Returns a 0-based array index for the given tensor indices.
inline t_ind index(t_ind i, t_ind j, t_ind k, t_ind l) const
{
return i-1 + n*(j-1 + n*(k-1 + n*(l-1)));
}
public:
//! \brief The default constructor creates a (scaled) identity tensor.
Tensor4(t_ind nsd = 3, Real scale = Real(1), bool makeJ = false);
//! \brief Constructor creating a tensor from a vector of components.
//! \details The provided vector is assumed to contain the components of the
//! matrix representation of the tensor, stored in a one-dimensional array.
Tensor4(const std::vector<Real>& x, t_ind nsd = 3);
//! \brief Sets \a this to the 0-tensor.
void zero() { std::fill(v.begin(),v.end(),Real(0)); }
//! \brief Type casting to a one-dimensional vector, for referencing.
operator const std::vector<Real>&() const { return v; }
//! \brief Type casting to a one-dimensional vector, for assignment.
operator std::vector<Real>&() { return v; }
//! \brief Reference through a pointer.
const Real* ptr() const { return &v.front(); }
//! \brief Index-1 based component reference.
const Real& operator()(t_ind i, t_ind j, t_ind d, t_ind l) const;
//! \brief Index-1 based component access.
Real& operator()(t_ind i, t_ind j, t_ind k, t_ind l);
//! \brief Assignment operator.
Tensor4& operator=(const Tensor4& T);
//! \brief Overloaded assignment operator.
Tensor4& operator=(Real val);
//! \brief Incrementation operator.
Tensor4& operator+=(Real val);
//! \brief Output stream operator.
friend std::ostream& operator<<(std::ostream& os, const Tensor4& T)
{
return T.print(os);
}
};
/*!
\brief Simple class for representing a symmetric fourth-order tensor.
*/
class SymmTensor4 : public Tensor4
{
protected:
//! \brief Auxilliary method used by the constructors.
virtual void redim(t_ind nsd);
//! \brief Prints out the tensor to an output stream.
virtual std::ostream& print(std::ostream& os) const;
private:
//! \brief Returns a 0-based array index for the given row and column indices.
//! \details Symmetric second-order tensors in 3D are assumed stored with the
//! following order in a one-dimensional array: s11, s22, s33, s12, s23, s13.
inline t_ind index(t_ind i, t_ind j) const
{
if (i == j)
return i-1; // diagonal term
else if (n == 2)
return 2; // off-diagonal term (2D)
else if (n == 1)
return 0;
if (i == j+1 || i+2 == j) std::swap(i,j);
return i+2; // upper triangular term (3D)
}
public:
//! \brief The default constructor creates an identity tensor.
SymmTensor4(t_ind nsd = 3, bool makeJ = false);
//! \brief Constructor creating a tensor from a vector of components.
//! \details The provided vector is assumed to contain the components of the
//! matrix representation of the tensor, stored in a one-dimensional array.
SymmTensor4(const std::vector<Real>& x, t_ind nsd = 3);
//! \brief Index-1 based component reference.
const Real& operator()(t_ind i, t_ind j, t_ind d, t_ind l) const;
//! \brief Index-1 based component access.
Real& operator()(t_ind i, t_ind j, t_ind k, t_ind l);
};
#endif

7
src/Utility/Test/TestTensor.C
View File

@@ -12,6 +12,7 @@
#include "Tensor.h"
#include "Vec3.h"
#include <array>
#include "gtest/gtest.h"
@@ -128,7 +129,7 @@ TEST(TestTensor, Principal)
ASSERT_FLOAT_EQ(p.y, 2.0);
ASSERT_FLOAT_EQ(p.z, 1.0);
Vec3Vec dir(3);
std::array<Vec3,3> dir;
T1.principal(p,dir.data());
ASSERT_FLOAT_EQ(p.x, 3.0);
ASSERT_FLOAT_EQ(p.y, 2.0);
@@ -171,4 +172,8 @@ TEST(TestTensor, Principal)
ASSERT_NEAR(dir[1].x, Trans3(1,1), 1.0e-15);
ASSERT_NEAR(dir[1].y, Trans3(2,1), 1.0e-15);
ASSERT_NEAR(dir[1].z, Trans3(3,1), 1.0e-15);
ASSERT_FLOAT_EQ(T1.vonMises(), sqrt(3.0));
ASSERT_FLOAT_EQ(T2.vonMises(), sqrt(3.0));
ASSERT_FLOAT_EQ(T3.vonMises(), sqrt(3.0));
}

53
src/Utility/Test/TestTensor4.C Normal file
View File

@@ -0,0 +1,53 @@
//==============================================================================
//!
//! \file TestTensor4.C
//!
//! \date Oct 29 2015
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Tests for fourth-order tensors with some basic operations.
//!
//==============================================================================
#include "Tensor4.h"
#include "gtest/gtest.h"
TEST(TestTensor4, Constructor)
{
unsigned short int i, j, k, l, n;
for (n = 1; n <= 3; n++)
{
SymmTensor4 I(n), J(n,true);
Tensor4 Is(n), Js(n,1.2,true);
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
for (k = 1; k <= n; k++)
for (l = 1; l <= n; l++)
{
if (i == j && j == k && k == l)
{
ASSERT_FLOAT_EQ( I(i,j,k,l), 1.0);
ASSERT_FLOAT_EQ(Is(i,j,k,l), 1.0);
}
else
{
ASSERT_FLOAT_EQ( I(i,j,k,l), 0.0);
ASSERT_FLOAT_EQ(Is(i,j,k,l), 0.0);
}
if (i == j && k == l)
{
ASSERT_FLOAT_EQ( J(i,j,k,l), 1.0);
ASSERT_FLOAT_EQ(Js(i,j,k,l), 1.2);
}
else
{
ASSERT_FLOAT_EQ( J(i,j,k,l), 0.0);
ASSERT_FLOAT_EQ(Js(i,j,k,l), 0.0);
}
}
}
}

121
src/Utility/Test/TestVec3Oper.C
View File

@@ -17,7 +17,7 @@
TEST(TestVec3Oper, GetAndSet)
{
Vec3 a(1,2,3);
Vec3 a(1.0,2.0,3.0);
EXPECT_TRUE(a[0] == 1.0);
EXPECT_TRUE(a[1] == 2.0);
EXPECT_TRUE(a[2] == 3.0);
@@ -38,141 +38,122 @@ TEST(TestVec3Oper, GetAndSet)
TEST(TestVec3Oper, MxV)
{
utl::matrix<Real> A(3,3);
A.fill(1.0);
std::vector<Real> x(3, 1.0);
Vec3 result = A*x;
EXPECT_TRUE(result[0] == 3.0);
EXPECT_TRUE(result[1] == 3.0);
EXPECT_TRUE(result[2] == 3.0);
A.fill(std::vector<Real>({1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0}).data());
std::vector<Real> x({1.0,2.0,3.0});
Vec3 result = A * x;
EXPECT_TRUE(result.x == 30.0);
EXPECT_TRUE(result.y == 36.0);
EXPECT_TRUE(result.z == 42.0);
}
TEST(TestVec3Oper, MultiplyScalar)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 result = a*10.0;
Vec3 result2 = 10.0*a;
EXPECT_TRUE(result == result2);
EXPECT_TRUE(result[0] == 10.0);
EXPECT_TRUE(result[1] == 10.0);
EXPECT_TRUE(result[2] == 10.0);
EXPECT_TRUE(result.x == 10.0);
EXPECT_TRUE(result.y == 20.0);
EXPECT_TRUE(result.z == 30.0);
}
TEST(TestVec3Oper, DivideScalar)
{
Vec3 a;
a[0] = a[1] = a[2] = 10.0;
Vec3 a(1.0,2.0,3.0);
Vec3 result = a/10.0;
EXPECT_TRUE(result[0] == 1.0);
EXPECT_TRUE(result[1] == 1.0);
EXPECT_TRUE(result[2] == 1.0);
EXPECT_TRUE(result.x == 0.1);
EXPECT_TRUE(result.y == 0.2);
EXPECT_TRUE(result.z == 0.3);
}
TEST(TestVec3Oper, Dot)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
Real result = a * b;
Real result = a*b;
EXPECT_TRUE(result == 3.0);
EXPECT_TRUE(result == 32.0);
}
TEST(TestVec3Oper, Addition)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
Vec3 result = a + b;
Vec3 result = a+b;
EXPECT_TRUE(result[0] == 2.0);
EXPECT_TRUE(result[1] == 2.0);
EXPECT_TRUE(result[2] == 2.0);
EXPECT_TRUE(result.x == 5.0);
EXPECT_TRUE(result.y == 7.0);
EXPECT_TRUE(result.z == 9.0);
}
TEST(TestVec3Oper, Subtraction)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
Vec3 result = a - b;
Vec3 result = a-b;
EXPECT_TRUE(result[0] == 0.0);
EXPECT_TRUE(result[1] == 0.0);
EXPECT_TRUE(result[2] == 0.0);
EXPECT_TRUE(result.x == -3.0);
EXPECT_TRUE(result.y == -3.0);
EXPECT_TRUE(result.z == -3.0);
}
TEST(TestVec3Oper, Length)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
double sqrt3 = sqrt(3);
Vec3 a(1.0,1.0,1.0);
double sqrt3 = sqrt(3.0);
EXPECT_FLOAT_EQ(a.length(), sqrt3);
}
TEST(TestVec3Oper, Equality)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 2.0;
Vec3 c;
c[0] = c[1] = c[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
Vec3 c(1.0,2.0,3.0);
EXPECT_TRUE( a == c);
EXPECT_TRUE (a == c);
EXPECT_FALSE(a == b);
}
TEST(TestVec3Oper, InEquality)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 2.0;
Vec3 c;
c[0] = c[1] = c[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
Vec3 c(1.0,2.0,3.0);
EXPECT_TRUE( a != b);
EXPECT_TRUE (a != b);
EXPECT_FALSE(a != c);
}
TEST(TestVec3Oper, Less)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 b;
b[0] = b[1] = b[2] = 2.0;
Vec3 a(1.0,2.0,3.0);
Vec3 b(4.0,5.0,6.0);
EXPECT_TRUE(a < b);
EXPECT_FALSE(b < a);
}
TEST(TestVec3Oper, StreamOut)
{
Vec3 a;
a[0] = a[1] = a[2] = 1.0;
Vec3 a(1.0,2.0,3.0);
std::stringstream str;
str << a;
EXPECT_STREQ(str.str().c_str(), "1 1 1");
EXPECT_STREQ(str.str().c_str(), "1 2 3");
}
TEST(TestVec3Oper, StreamIn)
{
std::stringstream str;
str << "1.0 1.0 1.0";
str << "1.0 2.0 3.0";
Vec3 result;
str >> result;
EXPECT_TRUE(result[0] == 1.0);
EXPECT_TRUE(result[1] == 1.0);
EXPECT_TRUE(result[2] == 1.0);
EXPECT_TRUE(result.x == 1.0);
EXPECT_TRUE(result.y == 2.0);
EXPECT_TRUE(result.z == 3.0);
}

51
src/Utility/Utilities.C
View File

@@ -305,57 +305,6 @@ int utl::gather (const std::vector<int>& index, size_t ir, size_t nr,
}
size_t utl::Pascal (int p, unsigned short int nsd)
{
size_t nM = 1;
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
if (nsd == 2)
nM++;
else for (int j = q; j >= 0; j--)
if (i+j <= q) nM++;
return nM;
}
void utl::Pascal (int p, Real x, Real y, std::vector<Real>& phi)
{
phi.clear();
phi.reserve(Pascal(p,2));
phi.push_back(Real(1));
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
{
int k, j = q-i;
Real a = Real(1);
for (k = 0; k < i; k++) a *= x;
for (k = 0; k < j; k++) a *= y;
phi.push_back(a);
}
}
void utl::Pascal (int p, Real x, Real y, Real z, std::vector<Real>& phi)
{
phi.clear();
phi.reserve(Pascal(p,3));
phi.push_back(Real(1));
for (int q = 1; q <= p; q++)
for (int i = q; i >= 0; i--)
for (int j = q; j >= 0; j--)
if (i+j <= q)
{
int l, k = q-i-j;
Real a = Real(1);
for (l = 0; l < i; l++) a *= x;
for (l = 0; l < j; l++) a *= y;
for (l = 0; l < k; l++) a *= z;
phi.push_back(a);
}
}
size_t utl::find_closest (const std::vector<Real>& a, Real v)
{
// The lower_bound function uses binary search to find the index of the value.

9
src/Utility/Utilities.h
View File

@@ -152,15 +152,6 @@ namespace utl
const std::vector<Real>& in, std::vector<Real>& out,
size_t offset_in = 0, int shift_idx = 0);
//! \brief Returns the number of monomials in Pascal's triangle.
//! \param[in] p Polynomial order (>= 0)
//! \param[in] nsd Number of spatial dimensions (2 or 3)
size_t Pascal(int p, unsigned short int nsd);
//! \brief Evaluates the monomials of Pascal's triangle in 2D for order \a p.
void Pascal(int p, Real x, Real y, std::vector<Real>& phi);
//! \brief Evaluates the monomials of Pascal's triangle in 3D for order \a p.
void Pascal(int p, Real x, Real y, Real z, std::vector<Real>& phi);
//! \brief Searches for a real value in an ordered array of reals.
//! \param[in] a The array of real values
//! \param[in] v The value to search for