/**
 * Mesh.cpp
 *
 *  Created on: May 10, 2017
 *      Author: vadims
 */

#include <istream>
#include <algorithm>
#include <optional>
#include "Mesh.h"
#include <TopoDS_Edge.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopExp.hxx>
#include <BRep_Tool.hxx>
#include <GeomLProp_CLProps.hxx>
#include <BRep_Tool.hxx>
#include "cenos_exception.h"

#include "Mesh_IO.h"
#include "Gmsh_IO.h"
#include "OFMesh_IO.h"
#include "ElmerMesh_IO.h"

Mesh::Mesh() 
{ 
	meshScale = 1; 
	meshDimension = Element::Dim::ONEDIM;
}

Mesh::~Mesh() 
{
}

void Mesh::clear()
{
	// do not clear mesh dimension and mesh scale!!!
	nodes.clear();
	elements.clear();
	entities.clear();
	filename.clear();
}


void Mesh::setMeshScale(double muf)
{
	meshScale = 1/ muf;
}

void Mesh::setWDir(std::string w)
{
	wDir = w;
}

void Mesh::setFileName(std::string fname)
{
	filename = fname;
}

void Mesh::reserveElementsSize(int s)
{
	elements.reserve(s);
}

void Mesh::reserveNodesSize(int s)
{
	nodes.reserve(s);
}

size_t Mesh::getSizeNodes()
{
	return nodes.size();
}

size_t Mesh::getSizeEntities()
{
	return entities.size();
}

size_t Mesh::getSizeElements()
{
	return elements.size();
}

void Mesh::addNode(std::shared_ptr<Node> n0)
{
	nodes.push_back(n0);
}

void Mesh::addEntity(std::shared_ptr<MeshEntity> e0)
{
	entities.push_back(e0);
}

void Mesh::addElement(std::shared_ptr<Element> e0)
{
	elements.push_back(e0);
	if (static_cast<int>(meshDimension) < static_cast<int>(e0->getDimension()))
		meshDimension = e0->getDimension();

	//add element and node ownership
	for (auto& nd : e0->getNodes())
	{
		nd->addOwnerElement(e0);
		e0->getOwnerEntity()->addNode(nd);
	}
	e0->getOwnerEntity()->addElement(e0);
	e0->setMidsideNodes();
}

NetgenMesh_ Mesh::getNetgenMesh()
{
	NetgenMesh_ ng_mesh = std::make_shared< NetgenMesh>();

	for (auto nd : this->nodes)
	{
		double p[3] = { nd->x(), nd->y(), nd->z() };
		ng_mesh->AddPoint(p);
	}

	//edge elements
	int ned = 0;
	std::map<int, int> old_to_new_edgeids;
	for (auto ent : entities)
	{
		if (ent->getDimension() == 1)
		{
			ned++;
			old_to_new_edgeids[ent->getId()] = ned;
		}
	}

	for (auto f : old_to_new_edgeids)
		std::cout << f.first << " " << f.second << std::endl;
	for (auto el : elements)
	{
		if (el->getDimension() != Element::Dim::ONEDIM)
			continue;
		ng_mesh->AddSegment(el->getNodes()[0]->getId(), el->getNodes()[1]->getId(), old_to_new_edgeids[el->getOwnerEntity()->getId()]);
	}


	//face elements
	int nfd = 0;
	std::map<int, int> old_to_new_faceids;
	for (auto ent : entities)
	{
		if (ent->getDimension() == 2)
		{
			nfd++;
			ng_mesh->AddFaceDescriptor(nfd, 0, 0, 0);
			old_to_new_faceids[ent->getId()] = nfd;
		}
	}

	for (auto el : elements)
	{
		if (el->getDimension() != Element::Dim::TWODIM)
			continue;
		// TO DO Fix for quads
		int nds[3] = { el->getNodes()[0]->getId(), el->getNodes()[1]->getId(), el->getNodes()[2]->getId() };
		ng_mesh->AddSurfaceElement(ElementTypes::TRI, nds, old_to_new_faceids[el->getOwnerEntity()->getId()]);
	}

	// volume elements
	for (auto el : elements)
	{
		if (el->getDimension() != Element::Dim::THREEDIM)
			continue;

		int nds[8];
		for (int n = 0; n < el->getNodes().size(); n++)
			nds[n] = el->getNodes()[n]->getId();
		ng_mesh->AddVolumeElement(ElementTypes::TETRA, nds, el->getOwnerEntity()->getId());
	}

	return ng_mesh;
}


NetgenMesh_ Mesh::getNetgenMesh(std::string ent_name, gp_Trsf trans)
{
	auto entity = this->getEntityByName(ent_name);
	auto m_nodes = entity->getNodes();
	auto m_elements = entity->getElements();

	// *****************
	// check if reversed
	// *****************
	bool reversed = false;
	if (entity->getType() == ShapeType::EDGE)
	{
		TopoDS_Edge edge = TopoDS::Edge(entity->getShape());
		TopoDS_Vertex vstart = TopExp::FirstVertex(edge);
		gp_Pnt pstart = BRep_Tool::Pnt(vstart);

		// find element which is attached to vstart
		Element_ start_element;
		for (auto el : m_elements)
		{
			for (auto nd : el->getNodes())
			{
				gp_Pnt n = gp_Pnt(nd->x(), nd->y(), nd->z());
				start_element = el;
				break;
			}
			if (start_element != nullptr)
				break;
		}

		//check if found
		if (start_element == nullptr)
			throw(cenos_exception("No elements attached to staring vertex found for entity " + entity->getName()));

		auto el_nodes = start_element->getNodes();
		gp_Pnt n1 = gp_Pnt(el_nodes[0]->x(), el_nodes[0]->y(), el_nodes[0]->z());
		gp_Pnt n2 = gp_Pnt(el_nodes[1]->x(), el_nodes[1]->y(), el_nodes[1]->z());

		gp_Dir element_dir = gp_Dir(gp_Vec(n1, n2));
		Standard_Real first, last;
		Handle(Geom_Curve) geom_curve = BRep_Tool::Curve(edge, first, last);

		GeomLProp_CLProps prop = GeomLProp_CLProps(geom_curve, 1, 1e-7);
		prop.SetParameter(first);
		gp_Dir dir;
		prop.Tangent(dir);

		if (edge.Orientation() == TopAbs_Orientation::TopAbs_REVERSED)
			dir.Reverse();

		// CP-1288 increased tolerance angle for opposite directions to 0.5 radians which is about 30 degrees
		// should it be more?
		//CP-1443 increased to 45.8 degrees
		if (element_dir.IsOpposite(dir, 0.8))
			reversed = true;

	}
	// *****************


	NetgenMesh_ ng_mesh = std::make_shared<NetgenMesh>();

	std::map<int, int> old_to_new_nodeids;
	int i = 1;
	for (auto nd : m_nodes)
	{
		gp_XYZ xyz(nd->x(), nd->y(), nd->z());
		trans.Transforms(xyz);
		double p[3] = { xyz.X(), xyz.Y(), xyz.Z()};

		ng_mesh->AddPoint( p);
		old_to_new_nodeids[nd->getId()] = i++;
	}


	if (m_elements[0]->getDimension() == Element::Dim::TWODIM)
		ng_mesh->AddFaceDescriptor(1, 0, 0, 0);

	for (auto el : m_elements)
	{
		if (el->getDimension() == Element::Dim::ONEDIM)
		{
			if (reversed)
				ng_mesh->AddSegment(old_to_new_nodeids[el->getNodes()[1]->getId()],
											old_to_new_nodeids[el->getNodes()[0]->getId()], 1);
			else
				ng_mesh->AddSegment(old_to_new_nodeids[el->getNodes()[0]->getId()],
					old_to_new_nodeids[el->getNodes()[1]->getId()], 1);
		}
		else if (el->getDimension() == Element::Dim::TWODIM)
		{
			if (el->getElementType() == ElementTypes::TRI)
			{
				int nds[3] = { old_to_new_nodeids[el->getNodes()[0]->getId()],
				old_to_new_nodeids[el->getNodes()[1]->getId()],
				old_to_new_nodeids[el->getNodes()[2]->getId()] };
				ng_mesh->AddSurfaceElement(ElementTypes::TRI, nds, 1);
			}
			else if (el->getElementType() == ElementTypes::QUAD)
			{
				int nds[4] = { old_to_new_nodeids[el->getNodes()[0]->getId()],
				old_to_new_nodeids[el->getNodes()[1]->getId()],
				old_to_new_nodeids[el->getNodes()[2]->getId()],
				old_to_new_nodeids[el->getNodes()[3]->getId()] };
				ng_mesh->AddSurfaceElement(ElementTypes::QUAD, nds, 1);
			}
		}
		else if (el->getDimension() == Element::Dim::THREEDIM)
		{
			int nds[8];
			for (int n = 0; n < el->getNodes().size(); n++)
				nds[n] = old_to_new_nodeids[el->getNodes()[n]->getId()];
			ng_mesh->AddVolumeElement(static_cast<ElementTypes>(el->getElementType()), nds, 1);		}
	}

	return ng_mesh;
}

int Mesh::getDimension()
{
	return meshDimension;
}

const std::vector<std::shared_ptr<Node>> Mesh::getNodes()
{
	return nodes;
}

const std::shared_ptr<Node> Mesh::getNode(int i)
{
 	return nodes[i-1];
}

const std::vector<std::shared_ptr<Node>> Mesh::getNodes(std::shared_ptr<MeshEntity> entity)
{
	std::vector<std::shared_ptr<Node>> nodesOfEntity = entity->getNodes();
	sort(nodesOfEntity.begin(), nodesOfEntity.end(), [](std::shared_ptr<Node> n1, std::shared_ptr<Node> n2) -> bool { return n1->getId() < n2->getId(); });
	return nodesOfEntity;
}

const std::vector<std::shared_ptr<Element> > Mesh::getElements()
{
	return elements;
}




double Mesh::getMeshScale()
{
	return meshScale;
}

std::string Mesh::getFileName()
{
	return filename;
}

std::vector<std::shared_ptr<MeshEntity>> Mesh::getEntities()
{
	return entities;
}

double Mesh::calculateArea(std::shared_ptr<MeshEntity>  ent)
{
	double area = 0;
	for (auto el : elements)
	{
		if (el->getOwnerEntity() == ent)
		{
			area = area + el->getArea();
		}
	}
	return area * meshScale * meshScale;
}

std::vector<double> Mesh::calculateCentroid(std::shared_ptr<MeshEntity>  ent) // currently area only
{
	double x = 0, y=0, z=0;
	double area = 0;
	if (ent->getType() == ShapeType::SOLID)
		return { 0,0,0 };

	for (auto el : elements)
	{
		if (el->getOwnerEntity() == ent)
		{
			std::vector<double> c = el->getCentroid();
			double elArea = el->getArea();
			x = x + c[0]* elArea;
			y = y + c[1]* elArea;
			z = z + c[2]* elArea;
			area = area + elArea;
		}
	}
	if (area == 0)
		return { 0, 0, 0 };
	else
		return { x * meshScale / area, y * meshScale / area, z * meshScale / area };
}


double Mesh::calculateVolume(std::shared_ptr<MeshEntity>  ent)
{
	double volume = 0;
	for (auto el : elements)
	{
		if (el->getOwnerEntity() == ent)
		{
			if (el->getDimension() == Element::THREEDIM)
				volume = volume + el->getVolume();
			else if (el->getDimension() == Element::TWODIM)
				volume = volume + el->getArea();
		}
	}
	return volume * meshScale * meshScale * meshScale;
}

double Mesh::calculateAxiVolume(std::shared_ptr<MeshEntity>  ent)
{
	double volume = 0;
	for (auto el : elements)
	{
		if (el->getOwnerEntity() == ent)
		{
			if (el->getDimension() == Element::THREEDIM)
				volume = volume + el->getVolume();
			else if (el->getDimension() == Element::TWODIM)
				volume = volume + el->getArea() * 2 * 3.14159 * el->getCentroid()[0] ; //area * 2 * PI * x
		}
	}
	return volume * meshScale* meshScale* meshScale;
}

void Mesh::readMesh(std::string& filename, std::string& format, char flags)
{
	Mesh_IO* reader;
	if (format == "GMSH")
		reader = new Gmsh_IO();
	else if (format == "OF")
		reader = new Gmsh_IO();
	else
		throw cenos_exception("Mesh format " + format + " not supported");

	reader->flags = flags;
	reader->readMesh(filename, this);
	delete reader;
}

void Mesh::writeMesh(std::string& filename, std::string& format, char flags)
{
	Mesh_IO* writer;
	if (format == "GMSH")
		writer = new Gmsh_IO();
	else if (format == "OF")
		writer = new OFMesh_IO();
	else if (format == "ELMER")
		writer = new ElmerMesh_IO();
	else
		throw cenos_exception("Mesh format " + format + " not supported");

	std::vector<std::string> selected_entities;
	for (auto ent : entities)
		if (ent->getDimension() == 3)
		{
			selected_entities.push_back(ent->getName());
			break;
		}

	writer->restricted_entity_names = selected_entities;
	writer->flags = flags;
	writer->writeMesh(filename, this);
	delete writer;
}

//calculates cohomology
void Mesh::addCut(std::shared_ptr<MeshEntity> ent)
{
	std::shared_ptr<MeshEntity>  ent0 = std::make_shared<MeshEntity>(ent->getShape(), ent->getName() + "_CH");
	ent0->setEnabled(false);
	entities.push_back(ent0);

	std::vector<std::shared_ptr<Node> > nodeList = getNodes(ent);

	for (auto& node : nodeList)
	{
		std::set<std::shared_ptr<Node>> nodeNeighbourSet;
		// nodes on surface can belong to several entities
		// to surface entity and volume entity

		for (auto& el : node->getOwnerElements())
		{

			if (ent->isSubShapeOf(el->getOwnerEntity()))
			{
				std::set<std::shared_ptr<Node> > connectedNodes = getConnectedNodes(node, el, nodeList);
				std::copy(connectedNodes.begin(), connectedNodes.end(), std::inserter(nodeNeighbourSet, nodeNeighbourSet.end()));
			}
		}

		for (auto& nnode : nodeNeighbourSet)
		{
			std::shared_ptr<Element>  e0 = Element::createElement(ElementTypes::LINE);
			e0->setId(getSizeElements() + 1);
			std::vector<std::shared_ptr<Node>> nodeList;
			nodeList.push_back(node);
			nodeList.push_back(nnode);
			e0->setNodes(nodeList);
			e0->setEntity(ent0);
			addElement(e0);
		}
	}
}

std::shared_ptr<MeshEntity> Mesh::getEntityById(int id)
{
	for (auto ent : entities)
	{
		if (ent->getId() == id)
			return ent;
	}
	return nullptr;
}

std::shared_ptr<MeshEntity> Mesh::getEntityByName(std::string n)
{
	for (auto ent : entities)
	{
		if (ent->getName() == n)
			return ent;
	}
	return nullptr;
}

MeshStats Mesh::getMeshStats()
{
	MeshEntityStats total_stats;
	total_stats.element_count = elements.size();
	total_stats.node_count = nodes.size();

	std::map<ElementTypes, std::vector<std::shared_ptr<Element>>> elementTypeMap;
	for (auto el : elements)
	{
		elementTypeMap[static_cast<ElementTypes>(el->getElementType())].push_back(el);
	}

	for (auto eltype : elementTypeMap)
	{
		total_stats.element_type_count[eltype.first] = eltype.second.size();
	}



	MeshStats stat;
	stat.total_stats = total_stats;

	for (auto ent : entities)
	{
		MeshEntityStats entity_stats;
		entity_stats.element_count = ent->getElements().size();
		entity_stats.node_count = ent->getNodes().size();

		std::map<ElementTypes, std::vector<std::shared_ptr<Element>>> ent_elementTypeMap;
		for (auto el : ent->getElements())
		{
			ent_elementTypeMap[static_cast<ElementTypes>(el->getElementType())].push_back(el);
		}

		for (auto eltype : ent_elementTypeMap)
		{
			entity_stats.element_type_count[eltype.first] = eltype.second.size();
		}
		stat.entity_stats[ent->getName()] = entity_stats;
	}

	return stat;
}


std::set<std::shared_ptr<Node>> Mesh::getConnectedNodes(std::shared_ptr<Node> node, std::shared_ptr<Element> element, std::vector<std::shared_ptr<Node>> surfaceNodeList)
{
	std::set<std::shared_ptr<Node>> retVec;

	for (auto& nd : element->getConnectedNodes(node) )
	{
		if (std::find(surfaceNodeList.begin(), surfaceNodeList.end(), nd) == surfaceNodeList.end())
			retVec.insert(nd);
	}
	return retVec;

}