402 lines
9.7 KiB
C++
402 lines
9.7 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _GCTL_TENSOR_H
|
|
#define _GCTL_TENSOR_H
|
|
|
|
#include "../core.h"
|
|
|
|
#include "point3c.h"
|
|
|
|
namespace gctl
|
|
{
|
|
template <typename T> struct type_tensor;
|
|
|
|
typedef type_tensor<double> tensor;
|
|
|
|
#ifndef IO_PSN
|
|
#define IO_PSN
|
|
// static variable for controlling the IO process
|
|
static int io_psn = 6;
|
|
#endif // IO_PSN
|
|
|
|
/**
|
|
* @brief Structure of a 3*3 tensor product.
|
|
*
|
|
* 3x3 的张量
|
|
*/
|
|
template <typename T>
|
|
struct type_tensor
|
|
{
|
|
T val[3][3];
|
|
// 构造函数与拷贝构造函数
|
|
type_tensor();
|
|
type_tensor(T a);
|
|
type_tensor(T a0, T a1, T a2, T a3, T a4, T a5, T a6, T a7, T a8);
|
|
type_tensor(const array<T> &p); // 数组的长度必须是9
|
|
type_tensor(const matrix<T> &p); // 二维数组的大小必须是3*3
|
|
type_tensor(const type_tensor<T> &b);
|
|
// 析构函数
|
|
virtual ~type_tensor(){}
|
|
// 结构体方法
|
|
void set(T a);
|
|
void set(T a0, T a1, T a2, T a3, T a4, T a5, T a6, T a7, T a8);
|
|
void set(const array<T> &p);
|
|
void set(const matrix<T> &p);
|
|
void set(const type_tensor<T> &b);
|
|
type_tensor<T> transpose() const;
|
|
// 元素获取方法
|
|
T &at(unsigned int row_index, unsigned int col_index);
|
|
T at(unsigned int row_index, unsigned int col_index) const;
|
|
T *operator[](unsigned int index);
|
|
};
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor()
|
|
{
|
|
for (int i = 0; i < 3; ++i)
|
|
{
|
|
for (int j = 0; j < 3; ++j)
|
|
{
|
|
val[i][j] = NAN;
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor(T a)
|
|
{
|
|
set(a);
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor(T a0, T a1, T a2, T a3, T a4, T a5, T a6, T a7, T a8)
|
|
{
|
|
set(a0, a1, a2, a3, a4, a5, a6, a7, a8);
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor(const array<T> &p)
|
|
{
|
|
set(p);
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor(const matrix<T> &p)
|
|
{
|
|
set(p);
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T>::type_tensor(const type_tensor<T> &b)
|
|
{
|
|
set(b);
|
|
}
|
|
|
|
template <typename T>
|
|
void type_tensor<T>::set(T a)
|
|
{
|
|
for (int i = 0; i < 3; ++i)
|
|
{
|
|
for (int j = 0; j < 3; ++j)
|
|
{
|
|
val[i][j] = a;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void type_tensor<T>::set(T a0, T a1, T a2, T a3, T a4, T a5, T a6, T a7, T a8)
|
|
{
|
|
val[0][0] = a0; val[0][1] = a1; val[0][2] = a2;
|
|
val[1][0] = a3; val[1][1] = a4; val[1][2] = a5;
|
|
val[2][0] = a6; val[2][1] = a7; val[2][2] = a8;
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void type_tensor<T>::set(const array<T> &p)
|
|
{
|
|
if (p.size() != 9)
|
|
throw runtime_error("Incompatible array size. From type_tensor::set(...)");
|
|
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
val[i][j] = p[3*i + j];
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void type_tensor<T>::set(const matrix<T> &p)
|
|
{
|
|
if (p.row_size() != 3 || p.col_size() != 3)
|
|
throw runtime_error("Incompatible matrix size. From type_tensor::set(...)");
|
|
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
val[i][j] = p[i][j];
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
void type_tensor<T>::set(const type_tensor<T> &b)
|
|
{
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
val[i][j] = b.val[i][j];
|
|
return;
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T> type_tensor<T>::transpose() const
|
|
{
|
|
type_tensor<T> m;
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
for(int j=0;j<3;j++)
|
|
{
|
|
m.val[i][j] = val[j][i];
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
T &type_tensor<T>::at(unsigned int row_index, unsigned int col_index)
|
|
{
|
|
return val[row_index][col_index];
|
|
}
|
|
|
|
template <typename T>
|
|
T type_tensor<T>::at(unsigned int row_index, unsigned int col_index) const
|
|
{
|
|
return val[row_index][col_index];
|
|
}
|
|
|
|
template <typename T>
|
|
T *type_tensor<T>::operator[](unsigned int index)
|
|
{
|
|
return val[index];
|
|
}
|
|
|
|
/**
|
|
* @brief 重载逻辑操作符, 执行两个 gctl::tensor 类型的加法。
|
|
*
|
|
* @param[in] a gctl::tensor 类型值，加数。
|
|
* @param[in] b gctl::tensor 类型值，加数。
|
|
*
|
|
* @return gctl::tensor 类型值，和。
|
|
*/
|
|
template <typename T>
|
|
type_tensor<T> operator +(const type_tensor<T> &a, const type_tensor<T> &b)
|
|
{
|
|
type_tensor<T> m;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
m[i][j] = a.val[i][j] + b.val[i][j];
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载逻辑操作符, 执行两个 gctl::tensor 类型的减法。
|
|
*
|
|
* @param[in] a gctl::tensor 类型值，减数。
|
|
* @param[in] b gctl::tensor 类型值，被减数。
|
|
*
|
|
* @return gctl::tensor 类型值，差。
|
|
*/
|
|
template <typename T>
|
|
type_tensor<T> operator -(const type_tensor<T> &a, const type_tensor<T> &b)
|
|
{
|
|
type_tensor<T> m;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
m[i][j] = a.val[i][j] - b.val[i][j];
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载运算操作符, 执行两个 gctl::tensor 类型的乘法。
|
|
*
|
|
* @param[in] a gctl::tensor 类型值，乘数。
|
|
* @param[in] b gctl::tensor 类型值，被乘数。
|
|
*
|
|
* @return gctl::tensor 类型值，积。
|
|
*/
|
|
template <typename T>
|
|
type_tensor<T> operator *(const type_tensor<T> &a, const type_tensor<T> &b)
|
|
{
|
|
type_tensor<T> m;
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
for (int j = 0; j < 3; j++)
|
|
{
|
|
m[i][j] = 0.0;
|
|
for (int k = 0; k < 3; k++)
|
|
{
|
|
m[i][j] += a.val[i][k] * b.val[k][j];
|
|
}
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载运算操作符, 张量乘向量。
|
|
*
|
|
* @param[in] m 张量
|
|
* @param[in] a 向量
|
|
*
|
|
* @return 积向量
|
|
*/
|
|
template <typename T>
|
|
point3c<T> operator *(const type_tensor<T> &m, const point3c<T> &a)
|
|
{
|
|
point3c<T> v;
|
|
v.x = a.x*m.val[0][0] + a.y*m.val[0][1] + a.z*m.val[0][2];
|
|
v.y = a.x*m.val[1][0] + a.y*m.val[1][1] + a.z*m.val[1][2];
|
|
v.z = a.x*m.val[2][0] + a.y*m.val[2][1] + a.z*m.val[2][2];
|
|
return v;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载运算操作符, 向量乘张量。
|
|
*
|
|
* @param[in] m 张量
|
|
* @param[in] a 向量
|
|
*
|
|
* @return 积向量
|
|
*/
|
|
template <typename T>
|
|
point3c<T> operator *(const point3c<T> &a, const type_tensor<T> &m)
|
|
{
|
|
point3c<T> v;
|
|
v.x = a.x*m.val[0][0] + a.y*m.val[1][0] + a.z*m.val[2][0];
|
|
v.y = a.x*m.val[0][1] + a.y*m.val[1][1] + a.z*m.val[2][1];
|
|
v.z = a.x*m.val[0][2] + a.y*m.val[1][2] + a.z*m.val[2][2];
|
|
return v;
|
|
}
|
|
|
|
/**
|
|
* @brief 标量成张量
|
|
*
|
|
* @param[in] d 标量
|
|
* @param[in] a 张量
|
|
*
|
|
* @return 积张量
|
|
*/
|
|
template <typename T>
|
|
type_tensor<T> operator *(double d, const type_tensor<T> &a)
|
|
{
|
|
type_tensor<T> m;
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
for(int j=0;j<3;j++)
|
|
{
|
|
m.val[i][j] = a.val[i][j]*d;
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
template <typename T>
|
|
type_tensor<T> operator *(const type_tensor<T> &a, double d)
|
|
{
|
|
type_tensor<T> m;
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
for(int j=0;j<3;j++)
|
|
{
|
|
m.val[i][j] = a.val[i][j]*d;
|
|
}
|
|
}
|
|
return m;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载输出流
|
|
*
|
|
* @param os 输出流
|
|
* @param[in] a 张量对象
|
|
*
|
|
* @return 输出流
|
|
*/
|
|
template <typename T>
|
|
std::ostream &operator <<(std::ostream & os, const type_tensor<T> &a)
|
|
{
|
|
os << std::setprecision(io_psn) << a.val[0][0] << " " << a.val[0][1] << " "
|
|
<< a.val[0][2] << " " << a.val[1][0] << " " << a.val[1][1] << " "
|
|
<< a.val[1][2] << " " << a.val[2][0] << " " << a.val[2][1] << " " << a.val[2][2];
|
|
return os;
|
|
}
|
|
|
|
/**
|
|
* @brief 重载输入流
|
|
*
|
|
* @param os 输入流
|
|
* @param a 张量对象
|
|
*
|
|
* @return 输入流
|
|
*/
|
|
template <typename T>
|
|
std::istream &operator >>(std::istream & os, type_tensor<T> &a)
|
|
{
|
|
os >> a.val[0][0] >> a.val[0][1] >> a.val[0][2]
|
|
>> a.val[1][0] >> a.val[1][1] >> a.val[1][2]
|
|
>> a.val[2][0] >> a.val[2][1] >> a.val[2][2];
|
|
return os;
|
|
}
|
|
|
|
/**
|
|
* @brief 两个矢量的张量积
|
|
*
|
|
* @param[in] a 输出的矢量a
|
|
* @param[in] b 输出的矢量b
|
|
*
|
|
* @return 返回的张量。
|
|
*/
|
|
template <typename T>
|
|
type_tensor<T> kron(const point3c<T> &a, const point3c<T> &b)
|
|
{
|
|
type_tensor<T> t;
|
|
t.val[0][0]=a.x*b.x; t.val[0][1]=a.x*b.y; t.val[0][2]=a.x*b.z;
|
|
t.val[1][0]=a.y*b.x; t.val[1][1]=a.y*b.y; t.val[1][2]=a.y*b.z;
|
|
t.val[2][0]=a.z*b.x; t.val[2][1]=a.z*b.y; t.val[2][2]=a.z*b.z;
|
|
return t;
|
|
}
|
|
}
|
|
|
|
#endif // _GCTL_TENSOR_H
|