#include "AutoMeshSize.h"
#include "misc/miscFunctions.hpp"
#include "topology/ShapeAlgorithms.h"
#include "setup/physics.hpp"
#include "cenos_exception.h"

AutoMeshSize::AutoMeshSize(std::shared_ptr<GeometryData> gd, std::shared_ptr<MeshGenerator> mg, json physics_json, json mesh_params) : geom_data(gd), mesh_generator(mg)
{
	// Read meshing config template.
	auto configDir = misc::getCenosConfigDir();
	if (!configDir.has_value())
	{
		throw cenos_exception("Meshing config file not found");
	}
	std::string meshing_config_filename = configDir.value() + "/jsons/_generated/meshing.json";

	std::ifstream conf_stream(meshing_config_filename);
	meshing_config = json::parse(conf_stream);

	std::string md = mesh_params["meshDensity"].get<std::string>();
	if (md == "rough")
		mesh_density_factor = 1.5;
	else if (md == "fine")
		mesh_density_factor = 0.55;
	else
		mesh_density_factor = 1.0;

	max_mesh_size = 1e10;  // initial value, later is set by analyzePhysics
	analyzePhysics(physics_json);

}

json cascadeSize(json params2D, json params1D) {
	// This prevents scenarios where edge size is greater than face size
	// or the face size is greater than the volume size.
	// Essentialy lower dimension sizing can only be finer, never coarser.
	// 
	// If the 2D max size is 5.0 and 1D max size is 7.5 it becomes 5.0
	// If the 2D max size is 5.0 and 1D max size is 3.5 it stays at 3.5

	double max2D = params2D["maxh"].get<double>();
	double max1D = params1D["maxh"].get<double>();
	double min1D = params1D["minh"].get<double>();

	if (max1D > max2D) {
		max1D = max2D;

		if (min1D > max1D) {
			min1D = max1D * 0.7;
		}
	}
	params1D["maxh"] = max1D;
	params1D["minh"] = min1D;
	return params1D;
}

json AutoMeshSize::getRefinmentJson(std::shared_ptr<Group> gr) {
	json params3D, params2D, params1D;
	// We will only add these params to the refinemnt json if they are not empty
	params3D = params2D = params1D = nullptr;

	double size_max = 0; // maximal allowable size for this volume <- what does this mean?
	double size_min = 0;

	if (gr->getEntities()[0]->getType() == ShapeType::SOLID) {
		double volume = ShapeAlgorithms::getVolume(gr->getShape());
		size_max = std::pow(volume / 1000.0, 1.0 / 3);
		if (size_max > max_mesh_size)
			size_max = max_mesh_size;
		size_min = std::pow(volume / 30000.0, 1.0 / 3);
	}
	else if (gr->getEntities()[0]->getType() == ShapeType::FACE) {
		double area = ShapeAlgorithms::getArea(gr->getShape());
		size_max = std::pow(area / 100.0, 1.0 / 2);
		if (size_max > max_mesh_size)
			size_max = max_mesh_size;
		size_min = std::pow(area / 1000.0, 1.0 / 2);
	}

	for (const auto& item : meshing_config.at("meshing").items()) {
		if (item.value().at("key").get<std::string>() == gr->getRole()) {
			json mesh_defaults = item.value().at("meshDefaults");

			// Add params3D not needed for 2D groups.
			if (gr->getEntities()[0]->getType() == ShapeType::SOLID) {

				if (mesh_defaults.find("mesh3D") == mesh_defaults.end()) {
					throw cenos_exception("mesh3D not found in mesh config for role: " + gr->getRole());
				}
				json conf_params_3D = mesh_defaults.at("mesh3D");
				std::string type_3D = conf_params_3D.at("type").get<std::string>();
				if (type_3D != "DEFAULT") {
					double sf = this->getScalingFactor(gr, type_3D);
					double maxh = std::min(max_mesh_size, conf_params_3D.at("maxSize").get<double>() * sf);
					double minh = conf_params_3D.at("minSize").get<double>() * sf;
					maxh = std::max(maxh, size_max);
					if (minh > maxh / 5)
						minh = maxh / 5;

					params3D = {
						{"maxh", 0.0},
						{ "minh", 0.0 },
						{ "grading", 0.0 },
						{ "entities", json::array()},
					};
					params3D["maxh"] = maxh;
					params3D["minh"] = minh;
					params3D["grading"] = conf_params_3D.at("grading").get<double>();
					
					json entities = json::array();
					for (auto ent : gr->getEntities()) {
						entities.push_back(ent->getName());
					}
					params3D["entities"] = entities;

					// domain_skinparams were defined in analyzePhysics
					//if (domain_skinparams.find(gr->getName()) != domain_skinparams.end()) {
					//	params3D["skin_layer"] = { 
					//		{"skin_min", domain_skinparams[gr->getName()].first},
					//		{"skin_max", domain_skinparams[gr->getName()].second}
					//	};
					//}
				}
			}


			if (mesh_defaults.find("mesh2D") == mesh_defaults.end()) {
				throw cenos_exception("mesh2D not found in mesh config for role: " + gr->getRole());
			}
			json conf_params_2D = mesh_defaults.at("mesh2D");
			std::string type_2D = conf_params_2D.at("type").get<std::string>();

			if (type_2D != "DEFAULT") {
				double sf = this->getScalingFactor(gr, type_2D);
				double maxh = std::min(max_mesh_size, conf_params_2D.at("maxSize").get<double>() * sf);
				double minh = conf_params_2D.at("minSize").get<double>() * sf;
				maxh = std::max(maxh, size_max);
				if (minh > maxh / 5)
					minh = maxh / 5;

				params2D["maxh"] = maxh;
				params2D["minh"] = minh;
				params2D["grading"] = conf_params_2D.at("grading").get<double>();
				params2D["entities"] = json::array();

				if (gr->getEntities()[0]->getType() == ShapeType::SOLID) {
					// Add the subshapes of solids if group is solid
					for (auto ent : geom_data->getEntities()) {
						if (ent->getType() == ShapeType::FACE and ent->isSubShapeOf(gr->getShape())) {
							params2D["entities"].push_back(ent->getName());
						}
					}
				} else {
					// otherwise just add the faces of the group, as it is 2D
					for (auto ent : gr->getEntities()) {
						params2D["entities"].push_back(ent->getName());
					}
					// domain_skinparams were defined in analyzePhysics
					//if (domain_skinparams.find(gr->getName()) != domain_skinparams.end()) {
					//	params2D["skin_layer"] = {
					//		{"skin_min", domain_skinparams[gr->getName()].first},
					//		{"skin_max", domain_skinparams[gr->getName()].second}
					//	};
					//}
				}
			}
			

			if (mesh_defaults.find("mesh1D") == mesh_defaults.end()) {
				throw cenos_exception("mesh1D not found in mesh config for role: " + gr->getRole());
			}

			json conf_params_1D = mesh_defaults.at("mesh1D");
			std::string type_1D = conf_params_1D.at("type").get<std::string>();
			if (type_1D != "REFINE") {
				double sf = this->getScalingFactor(gr, type_1D);
				double maxh = std::min(max_mesh_size, conf_params_1D.at("maxSize").get<double>() * sf);
				double minh = conf_params_1D.at("minSize").get<double>() * sf;
				maxh = std::max(maxh, size_max);
				if (minh > maxh / 5)
					minh = maxh / 5;

				params1D["maxh"] = maxh;
				params1D["minh"] = minh;
				params1D["grading"] = conf_params_1D.at("grading").get<double>();
				params1D["entities"] = json::array();
			}
			else {
				params1D["maxh"] = params2D["maxh"] * conf_params_1D.at("scale").get<double>();
				params1D["minh"] = params2D["minh"] * conf_params_1D.at("scale").get<double>();
				params1D["grading"] = 0.4;
				params1D["entities"] = json::array();
			}

			for (auto ent : geom_data->getEntities()) {
				if (ent->getType() == ShapeType::EDGE and ent->isSubShapeOf(gr->getShape())) {
					params1D["entities"].push_back(ent->getName());
				}
			}
			
		}
	}

	json refinment = nullptr;
	if (params3D != nullptr or params2D != nullptr or params1D != nullptr) {
		refinment = {
			{"key", gr->getName()},
			{"label", gr->getLabel()},
			{"source", gr->getName()},
			{"sizeMultiplier", mesh_density_factor},
		};
	}

	if (params3D != nullptr and params2D != nullptr) {
		params2D = cascadeSize(params3D, params2D);
	}
	if (params2D != nullptr and params1D != nullptr) {
		params1D = cascadeSize(params2D, params1D);
	}

	if (params3D != nullptr){
		refinment["params3D"] = params3D;
	}
	if (params2D != nullptr) {
		refinment["params2D"] = params2D;
	}
	if (params1D != nullptr) {
		refinment["params1D"] = params1D;
	}

	return refinment;
}

json AutoMeshSize::generateMeshSetupJson() {
	// do meshing on sorted domains
	// start with lowest priority, because sizing gets overwritten

	// This is needed because "DomainGroup" is needed instead of "Group"
	std::vector<std::shared_ptr<Group>> domains;
	for (auto dom : geom_data->getDomains()) {
		domains.push_back(dom);
	}
	std::vector<std::shared_ptr<Group>> sorted_domains = getSortedGroups(domains);
	
	json refinments = json::array();

	for (std::shared_ptr<Group> dom : sorted_domains) {
		json ref = getRefinmentJson(dom);
		if (ref != nullptr) {
			refinments.push_back(ref);
		}
	}
	for (auto bnd : geom_data->getBoundaries()) {
		json ref = getRefinmentJson(bnd);
		if (ref != nullptr) {
			refinments.push_back(ref);
		}
	}

	json meshSetup = {
	{"element_order", 1},
	{"refinments", refinments},
	};
	return meshSetup;
}


void AutoMeshSize::applyMeshSetupJson(json meshSetup) {
	for (json::iterator r = meshSetup["refinments"].begin(); r != meshSetup["refinments"].end(); ++r) {
		json ref = r.value();

		for (std::string params : {"params3D", "params2D", "params1D"}) {
			// 2D group refinents wont have "params3D"
			if (ref.contains(params)) {

				double max_size = ref[params]["maxh"].get<double>();
				double min_size = ref[params]["minh"].get<double>();
				double grading = ref[params]["grading"].get<double>();
				std::vector<std::string> entitity_names = ref[params]["entities"].get<std::vector<std::string>>();

				if (!entitity_names.empty() and max_size != 0.0) {
					std::shared_ptr<MeshParameters> mp = std::make_shared<MeshParameters>();
					if (max_size == 0.0 or grading == 0.0) {
						throw cenos_exception("One of the meshing parameters is 0 for refinemnt: " + ref["key"]);
					}
					mp->setMaxH(max_size);
					mp->setMinH(min_size);
					mp->setGrading(grading);
					if (ref[params]["skin_layer"] != nullptr) {
						mp->calculateLayers(ref[params]["skin_layer"]["skin_min"], ref[params]["skin_layer"]["skin_max"]);
					}

					for (std::string ent_name : entitity_names) {
						TopoEntity_ ent = geom_data->getEntityByName(ent_name);
						mesh_generator->setParameters(mp, ent);
					}
				}
			}
		}
	}
}

json AutoMeshSize::applyMultiplier(json meshSetup) {
	int i = 0;
	for (json::iterator r = meshSetup["refinments"].begin(); r != meshSetup["refinments"].end(); ++r) {
		json ref = r.value();
		double sizeMultiplier = ref["sizeMultiplier"];
		ref.erase("sizeMultiplier");

		for (std::string params : {"params3D", "params2D", "params1D"}) {
			if (ref.contains(params)) {
				ref[params]["maxh"] = ref[params]["maxh"].get<double>() * sizeMultiplier;
				ref[params]["minh"] = ref[params]["minh"].get<double>() * sizeMultiplier;
			}
		}
		meshSetup["refinments"][i] = ref;
		i++;
	}
	return meshSetup;
}

double AutoMeshSize::getScalingFactor(std::shared_ptr<Group> gr, std::string type)
{
	if (type == "MIN_EDGE_LENGTH")
	{
		return ShapeAlgorithms::getMinEdgeLength(gr->getShape());
	}
	else if (type == "DIAGONAL")
	{
		return ShapeAlgorithms::getDiagonal(gr->getShape());
	}
}

// returns a vector or groups reordered according to priority of roles in meshing
// list starts with smallest priority, ends with largest
std::vector<std::shared_ptr<Group>> AutoMeshSize::getSortedGroups(std::vector<std::shared_ptr<Group>> groups)
{
	std::vector<std::shared_ptr<Group>> sorted_groups;
	int current_priority = 0;

	while (!groups.empty())
	{
		auto it = groups.begin();
		while (it != groups.end())
		{
			bool f = false;
			for (const auto& item : meshing_config.at("meshing").items())
			{
				if (item.value().at("key").get<std::string>() == (*it)->getRole() and
					item.value().at("priority").get<int>() == current_priority)
				{
					sorted_groups.push_back(*it);
					it = groups.erase(it);
					f = true;
					break;
				}
			}
			if (!f)
				++it;
		}
		current_priority++;
		
		// 5 is the highest priority in meshing config
		if (current_priority > 5)
			break;
	}

	return sorted_groups;
}

// Sets max_mesh_size and domain_skinparams based on 
// the physics (frequency, wavelength, material properties)
void AutoMeshSize::analyzePhysics(json physics_json)
{
	for (auto ph : physics_json["physicsList"])
	{
		Physics phys;
		phys.readInput(ph);
		phys.validateInput(); // throws exception if something is not ok

		if (phys.getPhysicsKey() == "physicsMicrowaves")
		{
			double frequency = phys.getTimeParameters()["frequency"].getDoubleValue() * phys.getTimeParameters()["frequency"].getFactor();
			double wavelength = 299792458.0 / frequency / lengthToSI(geom_data->getLengthUnit()); //  speed of light over frequency = wavelength
			max_mesh_size = wavelength / 5;
		}
		else if (phys.getSolverKey() == "getdpEMGeneral2D" or phys.getSolverKey() == "getdpEMGeneral3D")
		{
			ShapeStats stat = geom_data->getShapeStatistics();
			max_mesh_size = sqrt(std::pow(stat.maxX - stat.minX, 2) + std::pow(stat.maxY - stat.minY, 2) + std::pow(stat.maxZ - stat.minZ, 2)) / 10;

			float frequency_min = phys.getTimeParameters()["frequency"].getMinValue();
			float frequency_max = phys.getTimeParameters()["frequency"].getMaxValue();

			for (auto dom : phys.getDomainsList())
			{
				if (dom->getTypeKey() == "em_cc" or dom->getTypeKey() == "source_dom")
				{
					auto mat = phys.getMaterialById(dom->getMaterialId());
					auto skin_layers = mat.getSkinLayerThicknessRange(frequency_min, frequency_max);
					double skin_min = skin_layers.first / lengthToSI(geom_data->getLengthUnit());
					double skin_max = skin_layers.second / lengthToSI(geom_data->getLengthUnit());
					domain_skinparams[dom->getName()] = std::make_pair(skin_min, skin_max);
				}
			}
		}
		else if (phys.getSolverKey() == "elmerRotatingMagnets")
		{
			ShapeStats stat = geom_data->getShapeStatistics();
			max_mesh_size = sqrt(std::pow(stat.maxX - stat.minX, 2) + std::pow(stat.maxY - stat.minY, 2) + std::pow(stat.maxZ - stat.minZ, 2)) / 10;
		}
	}
}