/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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