You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
263 lines
6.1 KiB
C
263 lines
6.1 KiB
C
1 year ago
|
/***************************************************************************
|
||
|
* Copyright (C) gempa GmbH *
|
||
|
* All rights reserved. *
|
||
|
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||
|
* *
|
||
|
* GNU Affero General Public License Usage *
|
||
|
* This file may be used under the terms of the GNU Affero *
|
||
|
* Public License version 3.0 as published by the Free Software Foundation *
|
||
|
* and appearing in the file LICENSE included in the packaging of this *
|
||
|
* file. Please review the following information to ensure the GNU Affero *
|
||
|
* Public License version 3.0 requirements will be met: *
|
||
|
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||
|
* *
|
||
|
* Other Usage *
|
||
|
* Alternatively, this file may be used in accordance with the terms and *
|
||
|
* conditions contained in a signed written agreement between you and *
|
||
|
* gempa GmbH. *
|
||
|
***************************************************************************/
|
||
|
|
||
|
|
||
|
#ifndef SEISCOMP_MATH_TENSOR_H
|
||
|
#define SEISCOMP_MATH_TENSOR_H
|
||
|
|
||
|
|
||
|
#include <seiscomp/math/matrix3.h>
|
||
|
#include <stdio.h>
|
||
|
|
||
|
|
||
|
namespace Seiscomp {
|
||
|
namespace Math {
|
||
|
|
||
|
|
||
|
template <typename T>
|
||
|
class Tensor2S;
|
||
|
|
||
|
|
||
|
/**
|
||
|
* Non-Symmetric second rank tensor (gradient & rotation tensors)
|
||
|
*/
|
||
|
template <typename T>
|
||
|
class Tensor2N {
|
||
|
public:
|
||
|
Tensor2N() { assign(0.0); }
|
||
|
|
||
|
//! Equality operator
|
||
|
void assign(Tensor2N<T> &other);
|
||
|
//! Initialize
|
||
|
void assign(T k);
|
||
|
|
||
|
//! Make unit tensor
|
||
|
void unit();
|
||
|
|
||
|
//! Add scaled tensor
|
||
|
void sum(const Tensor2N<T> &other, T k);
|
||
|
|
||
|
//! Multiply by a scalar
|
||
|
void scale(T k);
|
||
|
|
||
|
//! Determinant of gradient
|
||
|
T det() const;
|
||
|
|
||
|
//! Make inverse tensor
|
||
|
void inverse(const Tensor2N<T>& src);
|
||
|
|
||
|
void derive(T *dX1,
|
||
|
T *dX2,
|
||
|
T *dX3,
|
||
|
T *disp,
|
||
|
int * inodes,
|
||
|
int n);
|
||
|
|
||
|
//! Spin tensor
|
||
|
void spin (Tensor2N<T> &U);
|
||
|
|
||
|
//! FD = [det(F)^-1/3] * F
|
||
|
void deviator(Tensor2N<T> &F);
|
||
|
|
||
|
//! Product of non-symmetric
|
||
|
//! tensor with symmetric tensor
|
||
|
void product(const Tensor2N<T> &N,
|
||
|
const Tensor2S<T> &S);
|
||
|
|
||
|
//! Product of two non-symmetric
|
||
|
//! tensors
|
||
|
void product(const Tensor2N<T> &A,
|
||
|
const Tensor2N<T> &B);
|
||
|
|
||
|
//! Performs polar decomposition of the deformation gradient
|
||
|
//! tensor into the right stretch tensor and the rotation tensor
|
||
|
void polarDecomp(Tensor2N<T> &F);
|
||
|
|
||
|
//! Hughes-Winget rotation increment
|
||
|
void HughesWinget(Tensor2N<T> &U);
|
||
|
|
||
|
|
||
|
public:
|
||
|
T _11, _12, _13;
|
||
|
T _21, _22, _23;
|
||
|
T _31, _32, _33;
|
||
|
};
|
||
|
|
||
|
|
||
|
typedef Tensor2N<float> Tensor2Nf;
|
||
|
typedef Tensor2N<double> Tensor2Nd;
|
||
|
|
||
|
|
||
|
|
||
|
/**
|
||
|
* Symmetric second rank tensor (stress & strain tensors)
|
||
|
*/
|
||
|
template <typename T>
|
||
|
class Tensor2S {
|
||
|
public:
|
||
|
Tensor2S() { assign(0.0); }
|
||
|
Tensor2S(T __11, T __12, T __13,
|
||
|
T __22, T __23,
|
||
|
T __33);
|
||
|
|
||
|
//! Equality operator
|
||
|
void assign(const Tensor2S<T> &A);
|
||
|
//! Initialize
|
||
|
void assign(T k = 0.0);
|
||
|
|
||
|
//! Make spheric tensor
|
||
|
void spheric(T k = 1.0);
|
||
|
|
||
|
//! Add scaled tensor
|
||
|
void sum(const Tensor2S<T> &A, T k = 1.0);
|
||
|
|
||
|
//! Multiply by a scalar
|
||
|
void scale(T k);
|
||
|
|
||
|
//! Strain tensor
|
||
|
void strain(const Tensor2N<T> &U);
|
||
|
|
||
|
//! Strain deviator
|
||
|
void deviator(const Tensor2N<T> &U);
|
||
|
|
||
|
//! Green-Lagrange strain tensor
|
||
|
void GLStrain(const Tensor2N<T> &F);
|
||
|
|
||
|
//! Left Cauchy-Green deformation tensor
|
||
|
void leftCG(const Tensor2N<T> &F);
|
||
|
|
||
|
//! Right Cauchy-Green deformation tensor
|
||
|
void rightCG(const Tensor2N<T> &F);
|
||
|
|
||
|
//! Square of symmetric tensor C^2 = C^t*C
|
||
|
void square(const Tensor2S<T> &C);
|
||
|
|
||
|
//! Stress Push-Forward (PF)
|
||
|
void pshFrwd(const Tensor2N<T> &F,
|
||
|
const Tensor2S<T> &S);
|
||
|
|
||
|
//! Spatial to unrotated tensor transf.
|
||
|
void unrotate(const Matrix3<T> &R);
|
||
|
|
||
|
//! Spatial to unrotated tensor transf.
|
||
|
void unrotate(const Matrix3<T> &R,
|
||
|
const Tensor2S<T> &S);
|
||
|
|
||
|
//! Unrotated to spatial transf.
|
||
|
void rotate(const Matrix3<T> &R);
|
||
|
|
||
|
//! Unrotated to spatial tensor transf.
|
||
|
void rotate(const Matrix3<T> &R,
|
||
|
const Tensor2S<T> &U);
|
||
|
|
||
|
//! Trace of tensor
|
||
|
T I1() const;
|
||
|
|
||
|
//! Second invariant of tensor
|
||
|
T I2() const;
|
||
|
|
||
|
//! Mean stress
|
||
|
T mean() const;
|
||
|
|
||
|
//! Euclidean norm of tensor
|
||
|
T norm() const;
|
||
|
|
||
|
//! Euclidean norm of deviator
|
||
|
T devnorm() const;
|
||
|
|
||
|
//! Make deviator & return trace
|
||
|
T dtrace();
|
||
|
|
||
|
//! Make deviator & return mean value
|
||
|
T dmean();
|
||
|
|
||
|
//! Augment deviator with spheric tensor
|
||
|
void augment(T p);
|
||
|
|
||
|
//! Update deviatoric stress
|
||
|
void update(Tensor2S<T> &d,
|
||
|
T G);
|
||
|
|
||
|
//! Prints tensor to file
|
||
|
void print(FILE * fl);
|
||
|
|
||
|
//! apply Jacoby rotation
|
||
|
void jacoby(T &pp,
|
||
|
T &qq,
|
||
|
T &pq,
|
||
|
T &rp,
|
||
|
T &rq);
|
||
|
|
||
|
//! compute eigenvalues
|
||
|
bool eigenval(T tol = 1e-15,
|
||
|
int itmax = 30);
|
||
|
|
||
|
public:
|
||
|
T _11, _12, _13;
|
||
|
T _22, _23;
|
||
|
T _33;
|
||
|
};
|
||
|
|
||
|
|
||
|
typedef Tensor2S<float> Tensor2Sf;
|
||
|
typedef Tensor2S<double> Tensor2Sd;
|
||
|
|
||
|
|
||
|
/**
|
||
|
* Spectral representation of symmetric tensor
|
||
|
*/
|
||
|
template <typename T>
|
||
|
class Spectral2S {
|
||
|
public:
|
||
|
Spectral2S() { assign(0.0); }
|
||
|
|
||
|
//! initialize
|
||
|
void assign(T k);
|
||
|
|
||
|
//! spectral decomposition
|
||
|
bool spect(const Tensor2S<T> &A,
|
||
|
T atol = 1e-12,
|
||
|
int itmax = 50);
|
||
|
//! sort principal values in decending order & permute principal directions
|
||
|
void sort();
|
||
|
|
||
|
//! sort absolute principal values in decending order & permute principal directions
|
||
|
void absSort();
|
||
|
|
||
|
//! compose symmetric tensor from spectral representation
|
||
|
void compose(Tensor2S<T> &A);
|
||
|
|
||
|
//! compute norm of spectrally decomposed tensor
|
||
|
T norm();
|
||
|
|
||
|
public:
|
||
|
T a1, a2, a3; // principal values of tensor
|
||
|
Vector3<T> n1, n2, n3; // principal directions of tensor
|
||
|
};
|
||
|
|
||
|
|
||
|
typedef Spectral2S<float> Spectral2Sf;
|
||
|
typedef Spectral2S<double> Spectral2Sd;
|
||
|
|
||
|
|
||
|
}
|
||
|
}
|
||
|
|
||
|
#endif
|