#ifndef _TETRA_VZ_H
#define _TETRA_VZ_H

#include "cmath"

#define ZERO 1e-20
#define G0 6.67408e-11

/**
 * Define the point structure in 3D space
 */

struct point3d
{
	double x, y, z;
	double module() const; // module
	point3d normal() const; // normal vector
};

double point3d::module() const
{
	return sqrt(x*x + y*y + z*z);
}

point3d point3d::normal() const
{
	point3d nor_v;
	double length = module();
	if (length <= ZERO)
	{
		nor_v.x = nor_v.y = nor_v.z = 0.0;
	}
	else
	{
		nor_v.x = x/length;
		nor_v.y = y/length;
		nor_v.z = z/length;
	}
	return nor_v;
}

point3d operator+(const point3d &a, const point3d &b) // Add
{
	point3d m;
	m.x = a.x+b.x;
	m.y = a.y+b.y;
	m.z = a.z+b.z;
	return m;
}

point3d operator-(const point3d &a, const point3d &b) // Minus
{
	point3d m;
	m.x = a.x-b.x;
	m.y = a.y-b.y;
	m.z = a.z-b.z;
	return m;
}

point3d operator*(double scale, const point3d &b) // Multiply
{
	point3d m;
	m.x = scale*b.x;
	m.y = scale*b.y;
	m.z = scale*b.z;
	return m;
}

double dot(const point3d &a, const point3d &b) // Dot
{
	return a.x*b.x+a.y*b.y+a.z*b.z;
}

double distance(const point3d &a, const point3d &b) // Distance
{
	return (a - b).module();
}

point3d cross(const point3d &a, const point3d &b) // Cross
{
	point3d v;
	v.x = a.y*b.z-a.z*b.y;
	v.y = a.z*b.x-a.x*b.z;
	v.z = a.x*b.y-a.y*b.x;
	if (fabs(v.x) <= ZERO) v.x = 0;
	if (fabs(v.y) <= ZERO) v.y = 0;
	if (fabs(v.z) <= ZERO) v.z = 0;
	return v;
}

/**
 * Define the tensor structure
 */

struct tensor
{
	double val[3][3];
};

tensor kron(const point3d &a, const point3d &b) // tensor product
{
	tensor 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;
}

/**
 * Define the tetrahedron structure
 * 
 * regular type of tetrahedron
 *          3
 *         /|\
 *        / | \    y
 *   z   /  |  \  /
 *   |  /   1   \
 *   | /  /   \  \
 *   |/ /       \ \
 *   0-------------2----> x
 * 
 * triangle list (anti-clockwise)
 * 0 1 2
 * 0 2 3
 * 0 3 1
 * 1 3 2
 * 
 * inverse type of tetrahedron
 *          3
 *         /|\
 *        / | \    y
 *   z   /  |  \  /
 *   |  /   2   \
 *   | /  /   \  \
 *   |/ /       \ \
 *   0-------------1----> x
 * 
 * triangle list (anti-clockwise)
 * 0 1 3
 * 0 2 1
 * 0 3 2
 * 1 2 3
 */
static int regular_order[12] = {0,1,2,0,2,3,0,3,1,1,3,2};
static int inverse_order[12] = {0,1,3,0,2,1,0,3,2,1,2,3};

struct tetrahedron
{
	bool regular_type;
	point3d *vec_ptr[4];
	tensor F[4];
	tensor E[12];
	double rho;
	double edglen[12];
	void initialize_tensors();
};

void tetrahedron::initialize_tensors()
{
	point3d v1, v2, v3, nf, ne;

	v1 = *vec_ptr[1] - *vec_ptr[0];
	v2 = *vec_ptr[2] - *vec_ptr[0];
	v3 = *vec_ptr[3] - *vec_ptr[0];

	int *triangle_order;
	if (dot(v3, cross(v1, v2)) <= 0.0)
	{
		regular_type = true;
		triangle_order = regular_order;
	}
	else
	{
		regular_type = false;
		triangle_order = inverse_order;
	}

	for (int i = 0; i < 4; i++)
	{
		v1 = *vec_ptr[triangle_order[1+i*3]] - *vec_ptr[triangle_order[i*3]];
		v2 = *vec_ptr[triangle_order[2+i*3]] - *vec_ptr[triangle_order[i*3]];
		nf = cross(v1, v2).normal();
		F[i] = kron(nf, nf);

		for (int j = 0; j < 3; j++)
		{
			edglen[j+i*3] = distance(*vec_ptr[triangle_order[j+i*3]], *vec_ptr[triangle_order[(j+1)%3+i*3]]);

			v3 = *vec_ptr[triangle_order[(j+1)%3+i*3]] - *vec_ptr[triangle_order[j+i*3]];
			ne = cross(v3, nf).normal();
			E[j+i*3] = kron(nf, ne);
		}
	}
	return;
}

double tetra_vz(tetrahedron *ele_ptr, point3d *op_ptr)
{
	int f,e;
	double Le,wf;
	double dv1,dv2;
	double face_sum,edge_sum;
	point3d re;
	point3d r_ijk[3];
	double L_ijk[3];

	int *triangle_order;
	if (ele_ptr->regular_type)
		triangle_order = regular_order;
	else
		triangle_order = inverse_order;

	face_sum = edge_sum = 0.0;
	for (f = 0; f < 4; f++)
	{
		r_ijk[0] = *ele_ptr->vec_ptr[triangle_order[3*f]] - *op_ptr;
		r_ijk[1] = *ele_ptr->vec_ptr[triangle_order[3*f+1]] - *op_ptr;
		r_ijk[2] = *ele_ptr->vec_ptr[triangle_order[3*f+2]] - *op_ptr;

		L_ijk[0] = r_ijk[0].module();
		L_ijk[1] = r_ijk[1].module();
		L_ijk[2] = r_ijk[2].module();

		wf =2*atan2(dot(r_ijk[0], cross(r_ijk[1],r_ijk[2])),
					L_ijk[0]*L_ijk[1]*L_ijk[2] + L_ijk[0]*dot(r_ijk[1],r_ijk[2]) + 
					L_ijk[1]*dot(r_ijk[2],r_ijk[0]) + L_ijk[2]*dot(r_ijk[0],r_ijk[1]));

		face_sum += (ele_ptr->F[f].val[2][0]*r_ijk[0].x + ele_ptr->F[f].val[2][1]*r_ijk[0].y + ele_ptr->F[f].val[2][2]*r_ijk[0].z) * wf;

		for (e = 0; e < 3; e++)
		{
			dv1 = distance(*ele_ptr->vec_ptr[triangle_order[e+3*f]], *op_ptr);
			dv2 = distance(*ele_ptr->vec_ptr[triangle_order[(e+1)%3+3*f]], *op_ptr);

			re = 0.5*(*ele_ptr->vec_ptr[triangle_order[e+3*f]] + *ele_ptr->vec_ptr[triangle_order[(e+1)%3+3*f]]) - *op_ptr;
			Le = log((dv1+dv2+ele_ptr->edglen[e+3*f])/(dv1+dv2-ele_ptr->edglen[e+3*f]));

			edge_sum += (ele_ptr->E[e+3*f].val[2][0]*re.x + ele_ptr->E[e+3*f].val[2][1]*re.y + ele_ptr->E[e+3*f].val[2][2]*re.z) * Le;
		}
	}
	return 1.0e+8*G0*(ele_ptr->rho)*(edge_sum - face_sum);
}

#endif // _TETRA_VZ_H