#pragma once

#include <algorithm>
#include <ostream>

class TetIntersectMgr
{
private:
	/* Definition of the tetrahedra */
    double p[4][3];          // Vertex coordinate
	double v[4];             // Vertex value

	/* Description of the intersected surface */
	int sp_surf;             // Shape
	double pnt_surf[4][3];   // Vertex coordinate
	int loc_surf[4][2];      // Vertex location

	/* Description of the remaining element below the intersected surface */
	int sp_vol;              // Shape
	int idx_vol;             // Cell index in the mesh
	double pnt_vol[6][3];    // Vertex coordinate
	int face_vol[5][4];      // Face composition of vertexes

public:
	TetIntersectMgr()
    {
		reset();
    }

	double *coordinate(int idx)
	{
		return p[idx];
	}

	double value(int idx) const
	{
		return v[idx];
	}

	void reset()
	{
		for (int i = 0; i < 4; i++)
		{
			v[i] = 0.0;
			for (int j = 0; j < 3; j++)
        		p[i][j] = 0.0;
		}
			
        sp_surf = -1;
		for (int i = 0; i < 4; i++)
		{
			for (int j = 0; j < 2; j++)
				loc_surf[i][j] = 0;

			for (int j = 0; j < 3; j++)
        		pnt_surf[i][j] = 0.0;
		}

        sp_vol = -1;
		for (int i = 0; i < 6; i++)
			for (int j = 0; j < 3; j++)
        		pnt_vol[i][j] = 0.0;
	}

    int set_coordinate(int idx, double x_, double y_, double z_)
    {
        if (idx < 0 || idx >= 4)
            return 1;

        p[idx][0] = x_;
        p[idx][1] = y_;
        p[idx][2] = z_;

        return 0;
    }

	void set_all_coordinate(const double *p0, const double *p1, const double *p2, const double *p3)
	{
		copy_coordinate(p0, p[0]);
		copy_coordinate(p1, p[1]);
		copy_coordinate(p2, p[2]);
		copy_coordinate(p3, p[3]);
	}

    int set_value(int idx, double val_)
    {
        if (idx < 0 || idx > 3)
            return 1;

        v[idx] = val_;

        return 0;
    }

	void set_all_value(double val0, double val1, double val2, double val3)
	{
		v[0] = val0;
		v[1] = val1;
		v[2] = val2;
		v[3] = val3;
	}

	void set_all_value(const double *v_)
	{
		std::copy(v_, v_+4, v);
	}

    int set_coordinate_and_value(int idx, double x_, double y_, double z_, double val_)
    {
        if (idx < 0 || idx > 3)
            return 1;

        p[idx][0] = x_;
        p[idx][1] = y_;
        p[idx][2] = z_;
        v[idx] = val_;
        return 0;
    }

    int intersect()
	{
		int err = 0;

        const int st = value_state();
		switch (st)
		{
		case 0b0000: // -,-,-,- : all negative
        case 0b1111: // +,+,+,+ : all positive
            sp_surf = 0;
			sp_vol = 0;
			break;
		
		case 0b0001: // -,-,-,+ : only one positive, triangle
		case 0b1110: // +,+,+,- : only one negative, triangle
			err = get_tri(3, 2, 1, 0);
			break;

		case 0b0010: // -,-,+,- : only one positive, triangle
		case 0b1101: // +,+,-,+ : only one negative, triangle
			err = get_tri(2, 1, 3, 0);
			break;

		case 0b0100: // -,+,-,- : only one positive, triangle
		case 0b1011: // +,-,+,+ : only one negative, triangle
			err = get_tri(1, 0, 3, 2);
			break;

		case 0b1000: // +,-,-,- : only one positive, triangle		
		case 0b0111: // -,+,+,+ : only one negative, triangle
			err = get_tri(0, 3, 1, 2);
			break;
       
		case 0b0011: // -,-,+,+ : two positive & two negative, quadrilateral
        case 0b1100: // +,+,-,- : two positive & two negative, quadrilateral
			err = get_quad(2, 3, 1, 0);
			break;
		
		case 0b0101: // -,+,-,+ : two positive & two negative, quadrilateral
		case 0b1010: // +,-,+,- : two positive & two negative, quadrilateral
			err = get_quad(1, 3, 0, 2);
			break;

		case 0b0110: // -,+,+,- : two positive & two negative, quadrilateral        
		case 0b1001: // +,-,-,+ : two positive & two negative, quadrilateral
			err = get_quad(1, 2, 3, 0);
			break;
                
		default:
			sp_surf = -1;
			sp_vol = -1;
			err = 1;
			break;
		}

		return err;
	}

	int report(std::ostream &out)
	{
		out << "Tetrahedra:" << std::endl;
		for (int i = 0; i < 4; i++)
			out << "[" << i << "]: p=(" << p[i][0] << ", " << p[i][1] << ", " << p[i][2] << "), v=" << v[i] << std::endl;

		out << "Interface: ";
		if (sp_surf == 3)
			out << "TRI" << std::endl;
		else if (sp_surf == 4)
			out << "QUAD" << std::endl;
		else
		{
			out << "ERR (" << sp_surf << ")" << std::endl;
			return 1;
		}

		for (int i = 0; i < sp_surf; i++)
			out << "[" << i << "]: (" << pnt_surf[i][0] << ", " << pnt_surf[i][1] << ", " << pnt_surf[i][2] << "), {" << loc_surf[i][0] << ", " << loc_surf[i][1] << "}" << std::endl;

		out << "Remaining element: ";
		if (sp_vol == 4)
			out << "TET" << std::endl;
		else if (sp_vol == 6)
			out << "PRISM" << std::endl;
		else
		{
			out << "ERR (" << sp_vol << ")" << std::endl;
			return 2; 
		}
		for (int i = 0; i < sp_vol; i++)
			out << "N[" << i << "]: (" << pnt_vol[i][0] << ", " << pnt_vol[i][1] << ", " << pnt_vol[i][2] << ")" << std::endl;

		return 0;
	}

private:
	static void copy_coordinate(const double *src, double *dst)
	{
		std::copy(src, src+3, dst);
	}

    int value_state() const
    {
        const int sign[4] = {
			v[0] > 0.0 ? 1 : 0,
			v[1] > 0.0 ? 1 : 0,
			v[2] > 0.0 ? 1 : 0,
			v[3] > 0.0 ? 1 : 0
		};

		const int st = (sign[0] << 3) | (sign[1] << 2) | (sign[2] << 1) | (sign[3]);

        return st;
    }

	// Calculate the ratio that gives the linear interpolation value f between f1 and f2
	static double calcInterpParam(double f, double f1, double f2) 
	{
		return (f - f1) / (f2 - f1);
	}

	static void interpVec(const double *pa, const double *pb, double t, double *p)
	{
		p[0] = (1.0 - t) * pa[0] + t * pb[0];
		p[1] = (1.0 - t) * pa[1] + t * pb[1];
		p[2] = (1.0 - t) * pa[2] + t * pb[2];
	}

	int get_tri(int iv0, int iv1, int iv2, int iv3)
	{
		sp_surf = 3;

		// Get the ratio for the intersected values in the edges
		const double t01 = calcInterpParam(0.0, v[iv0], v[iv1]);
		const double t02 = calcInterpParam(0.0, v[iv0], v[iv2]);
		const double t03 = calcInterpParam(0.0, v[iv0], v[iv3]);

		// Interpolate values and get the triangle vertices
		interpVec(p[iv0], p[iv1], t01, pnt_surf[0]); loc_surf[0][0] = iv0; loc_surf[0][1] = iv1;
		interpVec(p[iv0], p[iv2], t02, pnt_surf[1]); loc_surf[1][0] = iv0; loc_surf[1][1] = iv2;
		interpVec(p[iv0], p[iv3], t03, pnt_surf[2]); loc_surf[2][0] = iv0; loc_surf[2][1] = iv3;

		// Determine the shape of the remaining volume
		if (v[iv0] > 0)
		{
			/* Prism */
			sp_vol = 6;
			copy_coordinate(p[iv1], pnt_vol[0]);
			copy_coordinate(p[iv2], pnt_vol[1]);
			copy_coordinate(p[iv3], pnt_vol[2]);
			copy_coordinate(pnt_surf[0], pnt_vol[3]);
			copy_coordinate(pnt_surf[1], pnt_vol[4]);
			copy_coordinate(pnt_surf[2], pnt_vol[5]);
		}
		else
		{
			/* Tet */
			sp_vol = 4;
			copy_coordinate(p[iv0], pnt_vol[0]);
			copy_coordinate(pnt_surf[0], pnt_vol[1]);
			copy_coordinate(pnt_surf[1], pnt_vol[2]);
			copy_coordinate(pnt_surf[2], pnt_vol[3]);
		}

		return 0;
	}

	int get_quad(int iv0, int iv1, int iv2, int iv3)
	{
		sp_surf = 4;

		// Get the ratio for the intersected values in the edges
		const double t02 = calcInterpParam(0.0, v[iv0], v[iv2]);
		const double t21 = calcInterpParam(0.0, v[iv2], v[iv1]);
		const double t13 = calcInterpParam(0.0, v[iv1], v[iv3]);
		const double t30 = calcInterpParam(0.0, v[iv3], v[iv0]);

		// Interpolate values and get the triangle vertices
		interpVec(p[iv0], p[iv2], t02, pnt_surf[0]); loc_surf[0][0] = iv0; loc_surf[0][1] = iv2;
		interpVec(p[iv2], p[iv1], t21, pnt_surf[1]); loc_surf[1][0] = iv2; loc_surf[1][1] = iv1;
		interpVec(p[iv1], p[iv3], t13, pnt_surf[2]); loc_surf[2][0] = iv1; loc_surf[2][1] = iv3;
		interpVec(p[iv3], p[iv0], t30, pnt_surf[3]); loc_surf[3][0] = iv3; loc_surf[3][1] = iv0;

		// Determine the shape of the remaining volume
		sp_vol = 6; /* Prism */
		if (v[iv0] > 0)
		{
			copy_coordinate(pnt_surf[0], pnt_vol[0]);
			copy_coordinate(pnt_surf[1], pnt_vol[1]);
			copy_coordinate(p[iv2], pnt_vol[2]);
			copy_coordinate(pnt_surf[3], pnt_vol[3]);
			copy_coordinate(pnt_surf[2], pnt_vol[4]);
			copy_coordinate(p[iv3], pnt_vol[5]);
		}
		else
		{
			copy_coordinate(pnt_surf[3], pnt_vol[0]);
			copy_coordinate(p[iv0], pnt_vol[1]);
			copy_coordinate(pnt_surf[0], pnt_vol[2]);
			copy_coordinate(pnt_surf[2], pnt_vol[3]);
			copy_coordinate(p[iv1], pnt_vol[5]);
			copy_coordinate(pnt_surf[1], pnt_vol[4]);
		}

		return 0;
	}
};