#include "GetdpSetup.hpp"
#include "GetdpMaterialWriter.hpp"
#include "GetdpResolutionFns.hpp"
#include <iostream> 
#include <boost/filesystem.hpp>
#include <boost/range/adaptor/indexed.hpp>
#include <fstream> 
#include "misc/Timer.hpp"
#include "misc/miscFunctions.hpp"
#include "GetdpPostFns.hpp"
#include "vtkForCenos/vtkForCenos.hpp"
#include "vtkForCenos/nearToFarField.h"
#include "vtkForCenos/circulation.h"
#include "mesh/meshToVtkWriter.hpp"
#include "module/SalomeModule.hpp"
#include "GetdpSpecialInputs.h"
#include <ryml.hpp>
#include <windows.h> 
#include <sstream> 
#include "setup/MaterialFunctions.h"
#include "SteelPhases/SteelPhases.h"
#include "SteelPhases/SPData.h"
#include "SteelPhases/SHData.h"
#include "BRepBuilderAPI_MakeVertex.hxx"
#include <cstdio>
#include "mesh/MeshGenerator.h"
#include "topology/geometryMover.h"
#include "misc/batchRunner.hpp"
#include <future>
#include "cenos_exception.h"
#include <filesystem>
#include "mesh/Mesh_IO.h"

std::string GetdpSetup::isA()
{
	return std::string("Getdp setup");
}


Numerics::ALGO GetdpSetup::getCouplingType()
{
	Numerics::ALGO coupling;

	// CP-195 returns Fast algorithm if single physics used
	int phys_count = 0;
	for (auto phys : casePhysics)
	{
		for (auto dom : phys.getDomainsList())
		{
			if (dom->isEnabled())
			{
				phys_count++;
				break;
			}
		}
	}
	if (phys_count == 1)
		return Numerics::ALGO::Fast;

	if (caseNumerics.getAlgorithm() == Numerics::ALGO::Fast)
		return Numerics::ALGO::Fast;
	else if (caseNumerics.getAlgorithm() == Numerics::ALGO::Accurate)
		return Numerics::ALGO::Accurate;
	else
	{
		// *************************************
		// 	   This part to be removed when power control with accurate is working
		bool power_control = false;
		for (auto phys : casePhysics)
		{
			for (auto dom : phys.getDomainsList())
			{
				for (auto prop : dom->getProperties())
				{
					if (prop.first == "_power")
						power_control = true;
				}
			}

			for (auto bnd : phys.getBoundariesList())
			{
				if (bnd->getTypeKey() == "power")
					power_control = true;
			}
		}
		if (power_control)
			return Numerics::ALGO::Fast;

		// *************************************

		for (auto& phys : casePhysics)
		{
			for (auto fl : phys.getFlags())
			{
				if ( (fl.isNLFlag()) and (fl.getValue()) )
					return Numerics::ALGO::Accurate;
			}
		}

		return Numerics::ALGO::Fast;
	}
}

std::string GetdpSetup::getDomainGroup(Domains_ dom)
{
	for (auto& domGroup : domainGroups)
	{
		for (auto subgr : domGroup.subgroups)
		{
			if ( (subgr.second.size() == 1) && (subgr.second[0] == "this") && (domGroup.key == dom->getTypeKey()) )
			{
				return subgr.first;
			}
		}
	}

	return "none";
}


std::string GetdpSetup::getBoundaryGroup(Bcs_ bc)
{
	for (auto& bcGroup : boundaryGroups)
	{
		std::string k = bcGroup.key;
		std::string m = bc->getTypeKey();

		if ((bcGroup.key == bc->getTypeKey()) )
		{
			return bcGroup.groupName;
		}

	}

	return "none";
}


// read functions
void GetdpSetup::readSolverConfig()
{
	Timer t;
	t.setMessage("readSolverConfig finished in: ");
	json physicsConfig = getPhysicsConfig();

	for (auto& phys : casePhysics)
	{
		json solverConfig = getSolverConfig(phys.getSolverKey(), phys.getPhysicsKey(), phys.getSubPhysicsKey(), phys.getBasisKey());

		json subPhysicsConfig = getSubPhysicsConfig(physicsConfig, phys.getPhysicsKey(), phys.getSubPhysicsKey(), phys.getBasisKey(), phys.getMultiSubPhysicsKey());

		phys.readPhysicsConfig(subPhysicsConfig);
		readConstraints(solverConfig);
		readBoundaryGroups(solverConfig);
		readBoundaryCohomology(solverConfig);
		readNetworks(solverConfig);
		readDomainGroups(solverConfig);
		readFlags(solverConfig);
		phys.readPostValues(solverConfig, postProcConfig);
	}
	Physics::readConstants(physicsConfig);
}

void GetdpSetup::readConfigurations()
{
	readSolverConfig();
}

void GetdpSetup::convertResults()
{
	//fixed to ENDTIME for now as runtime is removed!
	PostType postType = PostType::ENDTIME;
	vtkForCenos vtuG;
	size_t portNumber = 1;

	vtuG.setOutput(wDir + "/results", "resFile");

	//get boundaries with defined variables
	std::vector<int> bcsWithPost;
	
	for (auto& phys : casePhysics)
	{
		if (phys.getSolverKey() == "getdpMicrowaves") {
			portNumber = std::stod(phys.getLastPort());
		}
		
		for (auto& postVal : phys.getEnabledPostValueList())
		{
			for (auto& bc : phys.getBoundariesList())
			{
				if ((bc->getSolverKey() == postVal->getSolverKey()) and (postVal->isDefinedInDomain(getBoundaryGroup(bc))))
				{
					if (std::find(bcsWithPost.begin(), bcsWithPost.end(), bc->getId()) == bcsWithPost.end())
					{
						bcsWithPost.push_back(bc->getId());
					}
				}
			}
		}
	}

	std::vector<double> tvec = getTimeList(PostType::ENDTIME);

	if (tvec.size() > 1)
	{
		vtuG.setTransientGeometry(true);
		// set continuation of results
		if (hasRestart())
		{
			vtuG.setContinuation();
		}
	}

	MeshToVtkWriter meshToVtk(mesh);
	meshToVtk.createMultiblockMesh(vtuG, bcsWithPost, geom_data);

	//skipping post of zeroth time step if hasRestart
	int startTstep = 0;
	if (tvec.size() > 1)
	{
		if (hasRestart())
		{
			startTstep = 1;
		}
	}


	//post processing for steel phases
	for (auto& phys : casePhysics)
	{
		for (auto& postVal : phys.getEnabledPostValueList())
		{

			//skip if it is not external
			if (!postVal->isExternal())
				continue;

			if (postVal->getExternalData().at("calculator").get<std::string>() != "steel_phases")
				continue;

			// flag to check if steel phases are calculated in at leasts one domain
			bool steel_phases_needed = false;

			for (auto& domain : phys.getDomainsList())
			{
				// if material of this domain has flag _phaseCalculation
				// write steel phase file phase_properties.yaml
				// then calculate phases
				Material m = phys.getMaterialById(domain->getMaterialId());

				if (m.getMaterialDataFlag("_phaseCalculation"))
				{
					steel_phases_needed = true;
					//skip if it is not main field of steel_phases
					if (!postVal->getExternalData().at("mainField").get<bool>())
						continue;

					writeSteelPhases(m, wDir + "/calculation/phase_properties.yaml");

					try
					{
						SPData spd(wDir + "/getdpresults/spdata_" + domain->getName() + ".spd");
						spd.setName(domain->getName());
						for (auto ent : mesh->getEntities())
						{
							for (auto de : domain->getEntities())
							{
								if (ent->getName() == de->getName())
									spd.addNodes(mesh->getNodes(ent), tvec[0]);
							}
						}

						calculateSteelPhases(postVal->getStringId(), wDir, spd, tvec);
					}
					catch (std::string& msg)
					{
						std::string message = "Could not calculate post processing value " + postVal->getName() + "::" + msg;
						std::cout << "WARNING: " << message << std::endl;
						LOG(WARNING) << message;
						postVal->setEnabled(false);
					}


				}
				else
				{
					continue;
				}
			}

			if (!steel_phases_needed)
			{
				postVal->setEnabled(false);
			}
		}

	}
	for (auto& phys : casePhysics)
	{
		for (auto& postVal : phys.getEnabledPostValueList())
		{

			//skip if it is not external
			if (!postVal->isExternal())
				continue;

			if (postVal->getExternalData().at("calculator").get<std::string>() != "steel_hardness")
				continue;

			// flag to check if steel hardness are calculated in at leasts one domain
			bool steel_hardness_needed = false;

			for (auto& domain : phys.getDomainsList())
			{
				// if material of this domain has flag _hardnessCalculation
				// write steel hardness file hardness_properties.yaml
				// then calculate hardness
				Material m = phys.getMaterialById(domain->getMaterialId());

				if (m.getMaterialDataFlag("_maxTempCalculation"))
				{
					try
					{
						SHData shd(wDir + "/getdpresults/shdata_" + domain->getName() + ".spd");
						shd.setName(domain->getName());
						for (auto ent : mesh->getEntities())
						{
							for (auto de : domain->getEntities())
							{
								if (ent->getName() == de->getName())
									shd.addNodes(mesh->getNodes(ent), tvec[0]);
							}
						}

						calculateSteelHardness(postVal->getStringId(), wDir, shd, tvec);
					}
					catch (std::string& msg)
					{
						std::string message = "Could not calculate post processing value " + postVal->getName() + "::" + msg;
						std::cout << "WARNING: " << message << std::endl;
						LOG(WARNING) << message;
						postVal->setEnabled(false);
					}
				}
				else
				{
					continue;
				}
				}

				//if (!steel_hardness_needed)
				//{
				//	postVal->setEnabled(false);
				//}
			}

		}

	//for(auto& tstep: tvec | boost::adaptors::indexed(startTstep))
	
	for (int i = startTstep; i < tvec.size(); i++)
	{
		double tstepValue = tvec[i];
		//check if tstep is already written
		// if it is - add to vtk object.
		// if it is runtime and not the last step - skip rest of the loop
		if (std::find(writtenTimeSteps.begin(), writtenTimeSteps.end(), tstepValue) != writtenTimeSteps.end())
		{
			vtuG.addTimeStepValue(tstepValue);
			if ((postType == PostType::RUNTIME) and (i != tvec.size() - 1))
			{
				continue;
			}
		}
		
		std::vector<std::string> postValuesForVTK;
		
		for (auto& phys : casePhysics)
		{
			for (auto& postVal : phys.getEnabledPostValueList())
			{
				if (postVal->needsExportToVTK())
				{
					postValuesForVTK.push_back(postVal->getName());
				}

				for (auto& domain : phys.getDomainsList())
				{
					if ( postVal->isDefinedInDomain(getDomainGroup(domain)) )
					{
						std::ifstream resStreamFile;
						resStreamFile.open(wDir + "/getdpresults/" + postVal->getStringId() + "_" + domain->getName() + "_" + std::to_string(i) + ".dat");
						if (postVal->getType() == PostValues::NODE_FIELD)
						{
							if (postVal->getNumComponents() == 1)
								vtuG.addCurrentPointFieldToBlock(resStreamFile, postVal->getName(), domain->getId());
							else if (postVal->getNumComponents() == 3)
								vtuG.addCurrentPointVectorFieldToBlock(resStreamFile, postVal->getName(), domain->getId());
						}
						resStreamFile.close();
					}
				}

				//boundaries
				for (auto& bc : phys.getBoundariesList())
				{
					if ( postVal->isDefinedInDomain(getBoundaryGroup(bc)) )
					{

						if (postVal->getType() == PostValues::BOUNDARY_NODE_FIELD)
						{
							std::ifstream resStreamFile;
							resStreamFile.open(wDir + "/getdpresults/" + postVal->getStringId() + "_" + bc->getName() + "_" + std::to_string(i) + ".dat");
							if (postVal->getNumComponents() == 1)
								vtuG.addCurrentPointFieldToBlock(resStreamFile, postVal->getName(), bc->getId());
							else if (postVal->getNumComponents() == 3)
								vtuG.addCurrentPointVectorFieldToBlock(resStreamFile, postVal->getName(), bc->getId());
							resStreamFile.close();
						}
						else if (postVal->getType() == PostValues::BOUNDARY_CELL_FIELD)
						{
							std::ifstream resStreamFile;
							resStreamFile.open(wDir + "/getdpresults/" + postVal->getStringId() + "_" + bc->getName() + "_" + std::to_string(i) + ".dat");
							if (postVal->getNumComponents() == 1)
								vtuG.addCurrentCellFieldToBlock(resStreamFile, postVal->getName(), bc->getId());
							//else if (postVal.getNumComponents() == 3)
								//vtuG.addCurrentCellVectorFieldToBlock(resStreamFile, postVal.getName(), bc.getId());
							resStreamFile.close();
						}
					}
				}
			}

			// write far field
			for (auto& postVal : phys.getEnabledPostValueList())
			{
				if (postVal->getType() == PostValues::FAR_FIELD)
				{
					Timer t;
					t.setMessage("NTFF duration");
					nearToFarField ntff;
					if (phys.getTimeType() == Physics::HARMONIC or phys.getTimeType() == Physics::HARMONIC_WITH_TIME)
					{
						ntff.setWaveNumber(2 * 3.1415 * std::stof(phys.getTimeParameters()["frequency"].getValue() ) * phys.getTimeParameters()["frequency"].getFactor() 
							/ std::stof(Physics::getConstants()["c0"].getValue()));
						ntff.setFrequency(std::stof(phys.getTimeParameters()["frequency"].getValue()));
					}
					if (phys.getTimeType() == Physics::MULTIHARMONIC)
					{
						ntff.setWaveNumber(2 * 3.1415 * tstepValue * phys.getTimeParameters()["frequency"].getFactor()  / std::stof(Physics::getConstants()["c0"].getValue()));
						ntff.setFrequency(tstepValue);
					}
					ntff.setPath(wDir + "/results");
					
					for (auto& domain : phys.getDomainsList())
					{
						if (postVal->isDefinedInDomain(getDomainGroup(domain)))
						{
							ntff.addEnclosingGrid(vtuG.getGrid(domain->getId()));
						}
						if (domain->getTypeKey() == "em_open")
							ntff.setSourceGrid(vtuG.getGrid(domain->getId()));
					}
					for (auto& boundary : phys.getBoundariesList())
					{
						if (postVal->isDefinedInDomain(getBoundaryGroup(boundary)))
						{
							ntff.addEnclosingGrid(vtuG.getGrid(boundary->getId()));
						}
					}


					
					//near to far field results
					vtkForCenos ntfVtk;
					ntfVtk.setOutput(wDir + "/results", "farfield");
					if (i != 0)
						ntfVtk.setContinuation(tvec[i - 1]);

					std::vector<std::string> blockNames;
					std::vector<int> ids;
					ids.push_back(1);
					blockNames.push_back("far_field");
					ids.push_back(2);
					blockNames.push_back("radiation_intensity");
					ids.push_back(3);
					blockNames.push_back("directivity");

					ntfVtk.setBlockNames(ids, blockNames);
					if (i == 0) {
						ntff.setFirstFreq(true);
					}
					else {
						ntff.setFirstFreq(false);
					}
					
					auto ntf_result = ntff.execute();
					for (int j = 0; j < ntf_result.size(); j++)
					{
						ntfVtk.addUGrid(ntf_result[j], j+1);
					}
					ntfVtk.writeTimeStep(tstepValue);
					ntfVtk.finalize();
					if (i == 0)
					{
						std::ofstream ofs(wDir + "/getdpResults/globals/radiatedPower_globals.dat", std::ofstream::out);
						ofs << tstepValue << " " << ntff.getRadiatedPower() / 2 << "\n";
						ofs.close();
					}
					else
					{
						std::ofstream ofs(wDir + "/getdpResults/globals/radiatedPower_globals.dat", std::ofstream::app);
						ofs << tstepValue << " " << ntff.getRadiatedPower() / 2 << "\n";
						ofs.close();
					}
				}

				else if (postVal->getType() == PostValues::CIRCULATION)
				{

					for (auto& boundary : phys.getBoundariesList())
					{
						if (postVal->isDefinedInDomain(getBoundaryGroup(boundary)))
						{
							circulation circ;
							circ.setPath(wDir + "/results");

							circ.setCirculationGrid(vtuG.getGrid(boundary->getId()));
							circ.setDirectionVector(std::stof(boundary->getBcPropertyValue("eDirX")), std::stof(boundary->getBcPropertyValue("eDirY")), std::stof(boundary->getBcPropertyValue("eDirZ")));
							for (auto& domain : phys.getDomainsList())
							{
								if (domain->isEnabled())
									circ.addSourceGrid(vtuG.getGrid(domain->getId()));
							}
							std::complex<double> circulationValue = circ.execute();
							std::ofstream outFile(wDir + "/getdpresults/globals/" + postVal->getStringId() + "_" + boundary->getName() + "_globals.dat");
							outFile << tstepValue << " " << real(circulationValue) << " " << imag(circulationValue);

						}
					}

				}
			}

		}
	
		vtuG.writeTimeStep(tstepValue);
		vtuG.writeLegacyVTK(tstepValue, postValuesForVTK);
		vtuG.newTimeStep();
	}


	if (postType == PostType::ENDTIME)
	{
		for (auto& phys : casePhysics)
		{
			for (auto& postVal : phys.getEnabledPostValueList())
			{
				
				vtuG.enterDataSet(postVal->getName(), postVal->getNumComponents());
			}
		}
		vtuG.finalize();
	}
}

void GetdpSetup::readConstraints(json solverConfig)
{

	json constraintJson = solverConfig.at("constraints").get<json>();
	for (json::iterator it = constraintJson.begin(); it != constraintJson.end(); ++it)
	{
		if (it.value().is_string())
		{
			constraints.push_back(constraint(it.value().get<std::string>(), solverConfig.at("key").get<std::string>()));
		}
		if (it.value().is_object())
		{
			json constraintObj = it.value();
			if (constraintObj.find("region") != constraintObj.end())
			{
				GetdpSetup::constraint c(constraintObj["name"].get<std::string>(), solverConfig.at("key").get<std::string>());
				c.region = constraintObj["region"].get<std::string>();
				constraints.push_back(c);
			}
			if (constraintObj.find("limitedTo") != constraintObj.end())
			{
				GetdpSetup::constraint c(constraintObj["name"].get<std::string>(), solverConfig.at("key").get<std::string>());
				c.limitedTo = constraintObj["limitedTo"].get<std::string>();
				constraints.push_back(c);
			}
			if (constraintObj.find("assignFromResolution") != constraintObj.end())
			{
				GetdpSetup::constraint c(constraintObj["name"].get<std::string>(), solverConfig.at("key").get<std::string>());
				c.assignFromResolution = constraintObj["assignFromResolution"].get<std::string>();
				constraints.push_back(c);
			}
		}
	}

	if (solverConfig.find("case_constants") != solverConfig.end())
	{
		json case_constant_json = solverConfig.at("case_constants").get<json>();
		for (json::iterator it = case_constant_json.begin(); it != case_constant_json.end(); ++it)
		{
			case_constants.push_back(it.value().get<std::string>());
		}
	}
}

void GetdpSetup::readBoundaryGroups(json solverConfig)
{
	json boundaryGroupJson = solverConfig.at("bc_groups").get<json>();
	for (json::iterator it = boundaryGroupJson.begin(); it != boundaryGroupJson.end(); ++it)
	{
		boundaryGroup bgr(it.key(), solverConfig.at("key").get<std::string>());
		bgr.groupName = it.value().get<std::string>();
		boundaryGroups.push_back(bgr);
	}
}

void GetdpSetup::readBoundaryCohomology(json solverConfig)
{
	if (solverConfig.find("needs_cohomology") != solverConfig.end()) {
		json cohomology = solverConfig.at("needs_cohomology").get<json>();

		for (auto& boundaryType : Bcs::existingBoundaryTypes)
		{
			auto it = std::find(cohomology.begin(), cohomology.end(), boundaryType.key);
			boundaryType.needs_cohomology = (it != cohomology.end());
		}
	}
}

void GetdpSetup::readNetworks(json solverConfig)
{
	if (solverConfig.find("networks") != solverConfig.end()) {
		json netwrk_js = solverConfig.at("networks").get<json>();
		for (json::iterator it = netwrk_js.begin(); it != netwrk_js.end(); ++it)
		{
			networks.push_back(GetdpNetwork(it.value()["name"].get<std::string>(),
				solverConfig.at("key").get<std::string>(), it.value()["typeKey"].get<std::string>()));
		}
	}
}

void GetdpSetup::readDomainGroups(json solverConfig)
{
	json domainGroupJson = solverConfig.at("domain_groups").get<json>();
	for (json::iterator it = domainGroupJson.begin(); it != domainGroupJson.end(); ++it)
	{
		domainGroup domGroup(it.key(), solverConfig.at("key").get<std::string>());
		json groupJson = it.value();
		std::string gprStr = groupJson.dump();

		for (auto kt = groupJson.begin(); kt != groupJson.end(); ++kt)
		{
			if (kt.value().is_object())
				domGroup.oneSideOf = std::pair<std::string, std::vector<std::string>>(kt.key(), kt.value().at("oneSideOf").get<std::vector<std::string> >());
			else
			{
				domGroup.subgroups[kt.key()] = kt.value().get<std::vector<std::string> >();

			}
		}


		domainGroups.push_back(domGroup);
	}
}

void GetdpSetup::readFlags(json solverConfig)
{
	if (solverConfig.find("flags") != solverConfig.end())
	{
		json flagsJson = solverConfig.at("flags").get<json>();
		if (flagsJson.find("boundary") != flagsJson.end())
		{
			json bndFlagsJson = flagsJson.at("boundary").get<json>();
			std::string str(bndFlagsJson.dump());

			for (json::iterator it = bndFlagsJson.begin(); it != bndFlagsJson.end(); ++it)
			{
				Flag fl(it.key(), it.value());
				fl.setSolverKey(solverConfig.at("key").get<std::string>());
				for (auto& phys : casePhysics)
				{
					if (fl.getSolverKey() == phys.getSolverKey())
						phys.addFlag(fl);
				}
			}
		}

		if (flagsJson.find("material") != flagsJson.end())
		{
			json matFlagsJson = flagsJson.at("material").get<json>();
			for (json::iterator it = matFlagsJson.begin(); it != matFlagsJson.end(); ++it)
			{
				Flag fl(it.key(), it.value());
				fl.setSolverKey(solverConfig.at("key").get<std::string>());
				for (auto& phys : casePhysics)
				{
					if (fl.getSolverKey() == phys.getSolverKey())
						phys.addFlag(fl);
				}
			}
		}
	}
}

bool GetdpSetup::isConstraint(std::string name)
{
	for (auto& cn : constraints)
	{
		if (name == cn.name)
			return true;
	}
	return false;
}

bool GetdpSetup::isCaseConstant(std::string name)
{
	return std::find(case_constants.begin(), case_constants.end(), name) != case_constants.end();
}



// write setup
void GetdpSetup::writeSetup()
{
	//FIX for correct initialization, need to think of better way
	hasTransientConstraints = true;
    // replace mesh if simulation is restarted
	bool needs_last_mesh = false;
	for (auto& phys : casePhysics)
	{
		for (auto& tp : phys.getTimeParameters())
		{
			if (tp.first == "continueFromLast")
			{
				if (tp.second.getValue() == "1")
					needs_last_mesh = true;
			}
		}
	}

	if (needs_last_mesh)
	{
		try 
		{
			mesh->clear();
			mesh->readMesh(wDir + "/geometry/lastMesh.msh", std::string("GMSH"));
		}
		catch (cenos_exception& e)
		{
			throw(cenos_exception("Could not open last mesh file " + wDir + "/geometry/lastMesh.msh. Cannot continue calculation!"));
		}

	}

	updateMeshEntities();

	std::vector<std::shared_ptr<MeshEntity>> cutList;
	for (auto& phys : casePhysics)
	{
		for (auto& boundary : phys.getBoundariesList())
		{
			for (auto& boundaryType : Bcs::existingBoundaryTypes)
			{
				if ((boundaryType.needs_cohomology) and (boundaryType.key == boundary->getTypeKey()))
				{
					json bnd_json;
					bnd_json["typeKey"] = boundary->getTypeKey() + "_CH";
					bnd_json["typeId"] = boundary->getTypeId();
					bnd_json["domainId"] = boundary->getDomainId();

					Bcs_ bnd = std::make_shared<Bcs>(boundary->getName() + "_CH", bnd_json, boundary->getSolverKey());
					for (auto ent : boundary->getEntities())
					{
						TopoEntity_ e_new = std::make_shared<TopoEntity>(ent->getShape(), ent->getName() + "_CH");
						e_new->setEnabled(false);
						bnd->addEntity(e_new);
						geom_data->addEntity(e_new);    // addEntity will assign correct id to e_new
						mesh->addCut(mesh->getEntityByName(ent->getName()));
						mesh->getEntityByName(ent->getName() + "_CH")->setId(e_new->getId());  // set correct id to mesh group for cut
					}
					bnd->update();
					bnd->setLabel(boundary->getLabel());
					phys.addBoundary(bnd);
					bnd->setRole("_dummy_");
					geom_data->addBoundaryGroup(bnd);
				}
			}
		}
	}

	updateBoundaryConditions();

	//update physics
	initial_step = true;
	current_time = 0;
	current_timestep = 0;
	end_time = 0;
	for (auto& phys : casePhysics)
	{
		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			auto timeParams = phys.getTimeParameters();
			end_time = timeParams["tend"].getDoubleValue();
			current_time = timeParams["tstart"].getDoubleValue();

			// timeParams["isAdaptive"] has scalar. 1 - true, 0 - false
			bool is_adaptive = timeParams["isAdaptive"].getValue() == "1";
			if (!is_adaptive)
			{
				if (timeParams["tstep"].getTypeId() == InputValue::SCALAR)
					current_timestep = timeParams["tstep"].getDoubleValue();
			}

			// if this is continuation of previous simulation
			// set current time to time from results
			if (phys.getTimeParameters().count("continueFromLast"))
			{
				if (std::stof(phys.getTimeParameters()["continueFromLast"].getValue()) == 1)
				{
					vtkForCenos vtuG;
					vtuG.setOutput(wDir + "/results", "resFile");
					current_time = vtuG.getLastTimeInSet();
				}
			}

			// change timeParams if motion is present
			bool has_complex_motion = false;
			for (auto mt : caseMotions)
			{
				// skip if disabled
				if (!mt.isEnabled())
					continue;

				//skip if simple motion
				if (mt.getType() == Motion::SIMPLE)
					continue;

				has_complex_motion = true;
			}

			if (has_complex_motion)
			{
				if (is_adaptive)
					throw(cenos_exception("Adaptive time step is currently not supported for motion"));
				if (current_timestep == 0)
					throw(cenos_exception("Time step as table is currently not supported for motion"));
				phys.setTimeParameter("tstart", current_time);
				if (current_time > 0)
					phys.setTimeParameter("tend", current_time + current_timestep);
				else
					phys.setTimeParameter("tend", current_time);

			}
		}

		//copy simple motion parameters to thermal physics
		if (phys.getSolverKey() == "getdpThermal")
		{
			for (auto motion : caseMotions)
			{
				// skip if disabled
				if (!motion.isEnabled())
					continue;

				//skip if NOT simple motion
				if (motion.getType() != Motion::SIMPLE)
					continue;

				for (auto& dom : phys.getDomainsList())
				{
					auto domnames = motion.getDomainNames();
					if (std::find(domnames.begin(), domnames.end(), dom->getName()) != domnames.end())
					{
						dom->addProperty("_velocityEffects", InputValue(1.0));
						auto v = motion.getVelocity();
						dom->addProperty("Velocity_x", InputValue(v[0]));
						dom->addProperty("Velocity_y", InputValue(v[1]));
						dom->addProperty("Velocity_z", InputValue(v[2]));
						dom->addProperty("Omega_Z", InputValue(motion.getAngularVelocity()));
					}
				}
			}
		}

	}

	//update post values
	for (auto phys : casePhysics)
	{
		for (auto& postVal : phys.getPostValues())
		{
			auto enabling_field = postVal->getEnablingField();
			if (enabling_field == "")
				continue;

			bool found_enabling_field = false;
			for (auto& dom : phys.getDomainsList())
			{
				auto props = dom->getProperties();
				if (props.find(enabling_field) != props.end())
					found_enabling_field = true;
			}

			for (auto& bnd : phys.getBoundariesList())
			{
				auto props = bnd->getProperties();
				std::vector<std::string> prop_names;
				for (auto pr : props)
					prop_names.push_back(pr.key);
				if (std::find(prop_names.begin(), prop_names.end(), enabling_field) != prop_names.end())
					found_enabling_field = true;
			}

			if (!found_enabling_field)
				postVal->setEnabled(false);
		}

	}

	setFlags();
	addNetworkDomain();
	std::filesystem::remove_all(wDir + "/getdpresults");//needed for multifreq compatible with multiport
	misc::createFolder(wDir + "/calculation");
//	std::string times_dat = wDir + "/getdpresults/times.dat";
//	std::remove(times_dat.c_str());

	writeMain(PostType::SOLUTION);
	writeJakobian();
	writeIntegration();
	writeGroup();
	writeFunctionSpace();
	writeFormulation();
	writeUserFunction();
	writeFunction();
	WriteConstraint();
	writeResolution();
	WritePostProc();
	WritePostOp(PostType::SOLUTION);

	//write calculation mesh ascii
	mesh->writeMesh(wDir + "/calculation/calcMesh.msh", std::string("GMSH"));	
}


// write functions
void GetdpSetup::writePostSetup()
{
	writeMain(PostType::ENDTIME);
	WritePostOp(PostType::ENDTIME);
}


void GetdpSetup::writeMain(PostType postType)
{
	std::ofstream mainFile;
	mainFile.open(wDir + "/calculation/main.pro");
	std::string outString = "Include \"group.pro\" \n"
		"Include \"userfunction.pro\" \n"
		"Include \"function.pro\" \n"
		"Include \"constraint.pro\" \n"
		"Include \"functionspace.pro\" \n"
		"Include \"jakobian.pro\" \n"
		"Include \"integration.pro\" \n"
		"Include \"formulation.pro\" \n"
		"Include \"resolution.pro\" \n"
		"Include \"postproc.pro\" \n";
	if (postType == PostType::SOLUTION)
		outString.append("Include \"postop_SOLUTION.pro\" \n");
	else if (postType == PostType::ENDTIME)
		outString.append("Include \"postop_ENDTIME.pro\" \n");
	else
		outString.append("Include \"postop_RUNTIME.pro\" \n");

	outString.append(" \n");
	mainFile << outString;
	mainFile.close();
}

void GetdpSetup::writeJakobian()
{
	std::string volJacobian, surfaceJacobian;

	//check only zeroth physics, as they should be equal in both physics!
	if (casePhysics[0].getSymmetryType() == Physics::symmetryType::AXISYMMETRIC_2D)
	{
		volJacobian = "VolAxi";
		surfaceJacobian = "SurAxi";
	}
	else
	{
		volJacobian = "Vol";
		surfaceJacobian = "Sur";
	}

	std::ofstream jakobFile;
	jakobFile.open(wDir + "/calculation/jakobian.pro");
	std::string outString = str( boost::format("Jacobian { \n"
		"  { Name JVol ; \n"
		"    Case { \n"
		"      { Region All ; Jacobian %1%; } \n"
		"    } \n"
		"  } \n"
		"  { Name JSur ; \n"
		"    Case { \n"
		"      { Region All ; Jacobian %2%; } \n"
		"    } \n"
		"  } \n"
		"} \n") % volJacobian % surfaceJacobian );
	jakobFile << outString;
	jakobFile.close();
}

void GetdpSetup::writeIntegration()
{
	char const* tmp = getenv("CENOS_CONFIG");
	if (!tmp) {
		throw cenos_exception("variable CENOS_CONFIG not set ");
	}
	std::string cenos_env(tmp);
	std::string solverID;
	std::string str;
	std::ifstream inFile;

	inFile.open(cenos_env + "\\jsons\\integration.pro");
	if (inFile.is_open()) {
		std::stringstream strStream;
		strStream << inFile.rdbuf();
		str.append(strStream.str());
	}
	else {
		throw cenos_exception("No integration configuration file found!");
	}
	inFile.close();

	std::ofstream integrFile;
	integrFile.open(wDir + "/calculation/integration.pro");
	integrFile << str;
	integrFile.close();
}

void GetdpSetup::writeGroup()
{
	std::ofstream groupFile;
	groupFile.open(wDir + "/calculation/group.pro");

	std::string outString = "Group{ \n";

	for (auto dom : geom_data->getDomains())
	{
		outString.append("  " + dom->getName() + " = Region[{");
		
		int i = 0;
		for (auto ent : dom->getEntities())
		{
			if (i > 0)
				outString.append(", ");

			auto meshent = mesh->getEntityByName(ent->getName());
			outString.append(std::to_string(meshent->getId()));
			i++;
		}
		outString.append("}]; \n");
	}

	//add new boundary group for each port type unifPortN	
	
	
	for (auto& bgr : boundaryGroups)
	{
	
		
		if (bgr.groupName == "Sur_UNIFPORT_N")
		{
			std::vector<string> ports;//vector to store portNumbers and check if it's boundaryGroup is already created
			std::string newName;//the name to update boundaries and surfaces
			boundaryGroup groupN = bgr;; //aux group that will be edited and pushed if a new port is detected
			for (auto& boundary : getPhysics(bgr.solverKey).getBoundariesList())
				//this is not efficient because it loops arond all boundaries
				//it should loop only on the boundary group, the same happens in the loop below
			{
				// to identify the port number
				auto val = (boundary->getBcPropertyValue("numberPort"));
				if (val != "0")
				{
					newName = "Sur_UNIFPORT_N_" + val;
					// to keep a group with all uniformPortN we will create a new virtual boundary
					auto boundaryN = boundary;
					boundaryN->setTypeKey(newName);
					getPhysics("getdpMicrowaves").addBoundary(boundaryN);
					if (!(std::find(ports.begin(), ports.end(), val) != ports.end())) {
						ports.push_back(val);
						groupN.groupName = newName;
						groupN.key = newName;
						boundaryGroups.push_back(groupN);
					}
				}
			}
		}
	}

	for (auto bnd : geom_data->getBoundaries())
	{
		outString.append("  " + bnd->getName() + " = Region[{");

		int i = 0;
		for (auto ent : bnd->getEntities())
		{
			if (i > 0)
				outString.append(", ");

			auto meshent = mesh->getEntityByName(ent->getName());
			outString.append(std::to_string(meshent->getId()));
			i++;
		}
		outString.append("}]; \n");
	}
	//turn multiInput true for multiple port, maybe this should be a function separated 
	for (auto& phys : casePhysics)
	{
		phys.multiInput = false;
		if (phys.getSolverKey() == "getdpMicrowaves") {
			for (auto& bgr : boundaryGroups)
			{
				for (auto& boundary : getPhysics(bgr.solverKey).getBoundariesList())
					//this is not efficient because it loops arond all boundaries
					//it should loop only on the boundary group, the same happens in the loop below
				{
					// to identify the port number
					auto val = (boundary->getBcPropertyValue("numberPort"));
					if (val != "0")
					{
						phys.multiInput = true;
					}
				}
			}
		}
	}

	


	// write boundary groups
	// three structures exist:
	// boundaries (Bcs) in Physics (std::vector<Bcs>)
	// boundaryGroups in GetdpSetup
	//loop over existing boundary types

	for (auto& bgr : boundaryGroups)
	{
		outString.append("  " + bgr.groupName + " = Region[{");
		int k = 0;
		//loop over boundaries in physics with solver key of bgr
		for (auto& boundary : getPhysics(bgr.solverKey).getBoundariesList())
		{
			if (boundary->getTypeKey() == bgr.key)
			{
				if (k != 0) outString.append(", ");
				outString.append(boundary->getName());
				k++;
			}
			// generate a general group to gather all unifPort
			std::string unifPortN = "Sur_UNIFPORT_N";
			if (bgr.groupName == unifPortN)
			{
				if ((boundary->getTypeKey()).find(unifPortN) != std::string::npos)
				{
					if (k != 0) outString.append(", ");
					outString.append(boundary->getName());
					k++;
				}
			}

		}
		outString.append("}]; \n");

		outString.append("  " + bgr.groupName + "_List() = GetRegions[" + bgr.groupName + "];\n");

	}
	//write a group with all uniforms ports, used to get number of ports in other .pro files
	// 	   IMPORTANT
	// this should be used only for get number of ports as it is filled with numbers that do not refer to faces!!!!!!!!!!!!!!!!!!!!!!!
	{


		outString.append("  Sur_UNIFPORT_N_N = Region[{");
		int k = 1;
		//quick fix. All simulations need a port. So if there are not multiple ports it means there is only one. 
		// So the list is filled always with Sur_UNIFPORT_N_N = Region[{1}]; 
		std::string k_string = std::to_string(k);
		outString.append(k_string);
		for (auto& bgr : boundaryGroups)
		{

			//check for multiple ports
			if ((bgr.groupName).find("Sur_UNIFPORT_N_") != std::string::npos) {
				if (k != 1) {
					outString.append(", ");
					std::string k_string = std::to_string(k);
					outString.append(k_string);
				}
				k++;
			}
		}

		outString.append("}]; \n");
		outString.append("  Sur_UNIFPORT_N_N_List() = GetRegions[Sur_UNIFPORT_N_N];\n");
	}


	// write domain groups
	for (auto& domgroup : domainGroups)
	{
		// loop over subgroups which define if domain has to be in certain groups
		for (auto& subgroupsMap : domgroup.subgroups)
		{
			outString.append("  " + subgroupsMap.first + " = Region[{");
			int k = 0;
			for (auto& subgroup : subgroupsMap.second)
			{
				for (auto& domain : getPhysics(domgroup.solverKey).getDomainsList())
				{
					if ((subgroup == "this") && (domain->getTypeKey() == domgroup.key) )
					{
						if (k != 0) outString.append(", ");
						outString.append(domain->getName());
						k++;
					}
				}

				for (auto& boundary : getPhysics(domgroup.solverKey).getBoundariesList())
				{
					if (subgroup == boundary->getTypeKey())
					{
						if (k != 0) outString.append(", ");
						outString.append(boundary->getName());
						k++;
					}
				}

			}
			outString.append("}]; \n");
		}

		//loop over oneSideOf
		if (domgroup.oneSideOf.second.size() > 0)
		{
			std::vector<std::string> oneSideGroup;
			for (auto& oneSideOf : domgroup.oneSideOf.second)
			{
				for (auto& boundary : getPhysics(domgroup.solverKey).getBoundariesList())
				{
					if (oneSideOf == boundary->getTypeKey())
					{

						outString.append("  " + boundary->getName() + "_oneSide = ElementsOf[{" + std::to_string(boundary->getDomainId() ) + "}, OnOneSideOf " + boundary->getName() + "]; \n");
						oneSideGroup.push_back(boundary->getName() + "_oneSide");
					}
				}
			}

			outString.append("  " + domgroup.oneSideOf.first + " = Region[{");
			int k = 0;
			for (auto gr : oneSideGroup)
			{
				if (k > 0)
					outString.append(", ");
				outString.append(gr);
				k++;
			}
			outString.append("}]; \n");
			outString.append("  " + domgroup.oneSideOf.first + "_List() = GetRegions[" + domgroup.oneSideOf.first + "];\n");

		}
	}

	//write nonlinear domains
	outString = outString + "  NL_domain = Region[{";
	int k = 0;
	bool has_nl = false;
	for (auto& phys : casePhysics)
	{
		for (auto& domain : phys.getDomainsList()) 
		{
			if (domain->isNLDomain())
			{
				for (auto ent : domain->getEntities())
				{
					outString.append(std::to_string(ent->getId()) + ",");
					has_nl = true;
				}
			}
		}
	}
	if (has_nl)
		outString.pop_back();
	outString.append("}]; \n");

	outString = outString + "  moved_domains = Region[{";
	bool has_motions = false;
	for (auto mt: caseMotions)
	{
		if (mt.getType() == Motion::COMPLEX)
		{
			for (auto dom_nm : mt.getDomainNames())
			{
				auto dom = geom_data->getDomainByName(dom_nm);
				outString.append(" " + dom->getName() + ",");
				has_motions = true;
			}
		}
	}
	if (has_motions)
		outString.pop_back();
	outString.append("}]; \n");


	outString = outString + "}\n";
	groupFile << outString;
	groupFile.close();
}

void GetdpSetup::writeFormulation()
{
	char const* tmp = getenv("CENOS_CONFIG");
	if (!tmp) {
		std::cout << "variable CENOS_CONFIG not set " << std::endl;
	}
	std::string cenos_env(tmp);
	std::string solverID;
	std::string str;
	std::ifstream inFile;


	for (auto& phys : casePhysics)
	{
		inFile.open(cenos_env + "\\jsons\\" + phys.getSolverKey() + "_formulation.pro");
		if (inFile.is_open()) {
			std::stringstream strStream;
			strStream << inFile.rdbuf();
			std::string strToAdd = strStream.str();
			std::string from = "cFormulation";
			std::string to;
			if (phys.getSolverKey() == "getdpMicrowaves") {
				to = phys.getSolverKey() + "_formulation~{i}";
			}
			else {
			    to = phys.getSolverKey() + "_formulation";
			}

			size_t start_pos = strToAdd.find(from);
			strToAdd.replace(start_pos, from.length(), to);
			str.append(strToAdd);
		}
		else {
			std::cout << "No function space configuration file found!" << std::endl;
		}
		inFile.close();
	}

	std::ofstream formulationFile;
	formulationFile.open(wDir + "/calculation/formulation.pro");
	formulationFile << str;
	formulationFile.close();

}

void GetdpSetup::writeFunctionSpace()
{
	char const* tmp = getenv("CENOS_CONFIG");
	if (!tmp) {
		std::cout << "variable CENOS_CONFIG not set " << std::endl;
	}
	std::string cenos_env(tmp);
	std::string solverID;
	std::string str;
	std::ifstream inFile;

	for (auto& phys : casePhysics)
	{
		inFile.open(cenos_env + "\\jsons\\" + phys.getSolverKey() + "_fspace.pro");
		if (inFile.is_open()) {
			std::stringstream strStream;
			strStream << inFile.rdbuf();
			str.append(strStream.str());
		}
		else {
			str = "No function space configuration file found!";
			return;
		}
		inFile.close();
	}

	std::ofstream funcspaceFile;
	funcspaceFile.open(wDir + "/calculation/functionspace.pro");
	funcspaceFile << str;
	funcspaceFile.close();
}

void GetdpSetup::writeUserFunction()
{
	//user function section
	std::ofstream user_function_file;
	user_function_file.open(wDir + "/calculation/userfunction.pro");
	std::string user_outstring = "Function{";

	user_outstring.append(solver_script);
	user_outstring.append("\n");

	user_outstring.append("}\n");
	user_function_file << user_outstring;
	user_function_file.close();
}


void GetdpSetup::writeFunction()
{

	std::ofstream functionFile;
	functionFile.open(wDir + "/calculation/function.pro");
	std::string outString = "Function{ \n";
	outString.append("  DefineFunction[");

	//write boundary properties in definition
	std::vector<std::string> definedFunctions;
	for (auto& boundaryType : Bcs::existingBoundaryTypes) 
	{
		// loop over boundary properties (boundary conditions)
		for (auto& property : boundaryType.properties) 
		{
			if (std::find(definedFunctions.begin(), definedFunctions.end(), property) == definedFunctions.end())
			{
				outString.append(property);
				outString.append(", ");
				definedFunctions.push_back(property);
			}
		}

		// loop over domain material properties
		for (auto& matProperty : boundaryType.materialProps)
		{
			if (std::find(definedFunctions.begin(), definedFunctions.end(), matProperty) == definedFunctions.end())
			{
				outString.append(matProperty);
				outString.append(", ");
				definedFunctions.push_back(matProperty);
			}
		}
	}

	//write domain properties in definition
	for (auto& domainType : Domains::existingDomainTypes)
	{
		// loop over domain properties (sources)
		for (auto& property : domainType.properties)
		{
			if (std::find(definedFunctions.begin(), definedFunctions.end(), property) == definedFunctions.end())
			{
				outString.append(property);
				outString.append(", ");
				definedFunctions.push_back(property);
			}
		}

		// loop over domain material properties
		for (auto& matProperty : domainType.materialProps)
		{
			if (std::find(definedFunctions.begin(), definedFunctions.end(), matProperty) == definedFunctions.end())
			{
				outString.append(matProperty);
				outString.append(", ");
				definedFunctions.push_back(matProperty);
			}
		}
	}

	//write constants in definition
	for (auto& constant : Physics::getConstants()) {
		if (std::find(definedFunctions.begin(), definedFunctions.end(), constant.first) == definedFunctions.end())
		{
			outString.append(constant.first);
			outString.append(", ");
			definedFunctions.push_back(constant.first);
		}
	}
	outString.append(" new_xyz];\n");

	// write constants
	for (auto& constant : Physics::getConstants()) {

		if (constant.second.isTable())
			outString.append("  " + constant.first + "[] = " + constant.second.getValue());
		else
			outString.append("  " + constant.first + " = " + constant.second.getValue());
		outString.append(";  \n");
	}

	outString.append("  GeomCoeff[] = " + std::to_string(casePhysics[0].getGeometricFactor()) + "; \n");

	if ( (casePhysics[0].getSymmetryType() == Physics::symmetryType::FULL_3D) or
		(casePhysics[0].getSymmetryType() == Physics::symmetryType::SECTOR_3D))
	{

		outString.append("  symmFactor[] = " + std::to_string(casePhysics[0].getSymmFactor()) + "; \n");
	}

	bool isHarmonic = false;

	// Frequency treatment: harmonic; time dependent frequency.
	bool frequencySet = false;
	for (auto phys : casePhysics)
	{
		if (phys.getTimeType() == Physics::HARMONIC or phys.getTimeType() == Physics::HARMONIC_WITH_TIME)
		{
			auto timeParams = phys.getTimeParameters();
			if (timeParams["frequency"].getTypeId() == InputValue::TABLE)
			{
				outString.append("  FreqFunc[] = " + timeParams["frequency"].getValue() + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
				outString.append("  Freq = " + timeParams["frequency"].getFirstPair().second + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
				outString.append("  FreqFactor = " + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
			}
			else
			{

				// TEMPORAL FIX FOR FreqFuncProblem
				// TODO Check why definition of FreqFunc is necessary in getdp!
				outString.append("  FreqFunc[] = " + timeParams["frequency"].getValue() + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
				outString.append("  Freq = " + timeParams["frequency"].getValue() + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
				outString.append("  FreqFactor = " + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
			}
			outString.append("  rmsCoeff[] = Sqrt[2]; \n");
			frequencySet = true;
		}
		else if  (phys.getTimeType() == Physics::MULTIHARMONIC)
		{
			auto timeParams = phys.getTimeParameters();
			outString.append("  FreqFunc[] = " + timeParams["frequency"].getValue() + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
			outString.append("  Freq = " + timeParams["frequency"].getValue() + "*" + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
			outString.append("  FreqFactor = " + std::to_string(timeParams["frequency"].getFactor()) + "; \n");
			outString.append("  rmsCoeff[] = Sqrt[2]; \n");
			frequencySet = true;
		}
		else if (phys.getTimeType() == Physics::TRANSIENT)
		{
			auto timeParams = phys.getTimeParameters();
			if (timeParams.find("tstep") != timeParams.end())
			{
				if (timeParams["tstep"].getTypeId() == InputValue::TABLE)
				{
					outString.append("  tsteps[] = " + timeParams["tstep"].getValue() + "; \n");
				}
			}
		}
	}

	if (!frequencySet)
	{
		outString.append("  Freq = 0; \n");
		outString.append("  FreqFactor = 0; \n");
		outString.append("  rmsCoeff[] = 1; \n");
	}

	//writing boundary functions
	std::vector<std::string> case_constants;
	for (auto& phys : casePhysics)
	{
		for (auto& boundary : phys.getBoundariesList())
		{
			//workaround for current vaue in stranded coil 3D
			if ((boundary->getTypeKey() == "current_stranded"))
			{
				for (auto domain : phys.getDomainsList())
				{
					if (boundary->getDomainId() == domain->getId())
						outString.append("  I_str[" + boundary->getName() + "_CH] = " + domain->getPropertyValue("I_str") + "; \n");
				}
			}

			if ((boundary->getTypeKey() == "uniformPort"))
			{
				if (boundary->getBcPropertyValue("autoFeedParams") == "1" )
					outString.append(uniformPort(boundary->getName(), phys.getBoundariesList(), phys.getDomainsList(), this->geom_data));
			}
			if ((boundary->getTypeKey()).find("uniformPortN"))
			{

				if (boundary->getBcPropertyValue("autoFeedParams") == "1") {
					//TODO
					outString.append(uniformPort(boundary->getName(), phys.getBoundariesList(), phys.getDomainsList(), this->geom_data));
				}
			}
			if ((boundary->getTypeKey() == "coaxialPort"))
			{
				if (boundary->getBcPropertyValue("autoFeedParams") == "1")
					outString.append(coaxialPort(boundary->getName(), phys.getBoundariesList(), phys.getDomainsList(), this->geom_data));
			}

			for (auto& prop : boundary->getProperties())
			{
				if (isCaseConstant(prop.key))
				{
					if (std::find(case_constants.begin(), case_constants.end(), prop.key) != case_constants.end())
					    continue;
					case_constants.push_back(prop.key);
					outString.append("  " + prop.key + "[] = " + prop.value.getValue() + "; \n");
				}
				else if (!isConstraint(prop.key) )
					outString.append("  " + prop.key + "[" + boundary->getName() + "] = " + prop.value.getValue() + "; \n");
			}
		}
	}

	//writing domain functions
	for (auto& phys : casePhysics)
	{
		for (auto& domain : phys.getDomainsList())
		{
			for (auto& prop : domain->getProperties())
			{
				if (isCaseConstant(prop.first) and (!domain->isPointDomain()))
				{
					if (std::find(case_constants.begin(), case_constants.end(), prop.first) != case_constants.end())
					    continue;
					case_constants.push_back(prop.first);
					outString.append("  " + prop.first + "[] = " + prop.second.getValue() + "; \n");
				}
				else if (!isConstraint(prop.first))
				{
					outString.append("  " + prop.first + "[" + domain->getName() + "] = " + prop.second.getValue() + "; \n");
				}
					
			}
		}
	}

	// write entity areas and volumes
	double gd_scale = lengthToSI(geom_data->getLengthUnit());
	for (auto dom : geom_data->getDomains())
	{

		ShapeStats stat = dom->getShapeStatistics();

		if (mesh->getDimension() == 2)
		{
			std::ostringstream areaObj;
			areaObj << stat.area * gd_scale * gd_scale;
			outString.append("  area[" + dom->getName() + "] = " + areaObj.str() + "; \n");
			std::ostringstream volObj;
			if (casePhysics[0].getSymmetryType() == Physics::symmetryType::AXISYMMETRIC_2D)
			{
				double vol = 0;
				for (auto ent : dom->getEntities())
				{
					auto mesh_ent = mesh->getEntityByName(ent->getName());
					vol = vol + mesh->calculateAxiVolume(mesh_ent);
				}
				volObj << vol;
				outString.append("  volume[" + dom->getName() + "] = " + volObj.str() + "; \n");
			}
			else if (casePhysics[0].getSymmetryType() == Physics::symmetryType::PLANAR_2D)
			{
				double vol_2d = stat.area * gd_scale * gd_scale * casePhysics[0].getGeometricFactor();
				volObj << vol_2d;
				outString.append("  volume[" + dom->getName() + "] = " + volObj.str() + "; \n");

			}
		}
		else if (mesh->getDimension() == 3)
		{
			std::ostringstream volObj;
			volObj << stat.volume * gd_scale * gd_scale * gd_scale;
			outString.append("  volume[" + dom->getName() + "] = " + volObj.str() + "; \n");
		}

		std::ostringstream centroidx, centroidy, centroidz;
		centroidx << stat.centerX * gd_scale;
		centroidy << stat.centerY * gd_scale;
		centroidz << stat.centerZ * gd_scale;
		outString.append("  centroidX[" + dom->getName() + "] = " + centroidx.str() + "; \n");
		outString.append("  centroidY[" + dom->getName() + "] = " + centroidy.str() + "; \n");
		outString.append("  centroidZ[" + dom->getName() + "] = " + centroidz.str() + "; \n");
	}


	for (auto bnd : geom_data->getBoundaries())
	{

		ShapeStats stat = bnd->getShapeStatistics();

		if (mesh->getDimension() == 3)
		{
			std::ostringstream areaObj;
			areaObj << stat.area * gd_scale * gd_scale;
			outString.append("  area[" + bnd->getName() + "] = " + areaObj.str() + "; \n");
		}

		std::ostringstream centroidx, centroidy, centroidz;
		centroidx << stat.centerX * gd_scale;
		centroidy << stat.centerY * gd_scale;
		centroidz << stat.centerZ * gd_scale;
		outString.append("  centroidX[" + bnd->getName() + "] = " + centroidx.str() + "; \n");
		outString.append("  centroidY[" + bnd->getName() + "] = " + centroidy.str() + "; \n");
		outString.append("  centroidZ[" + bnd->getName() + "] = " + centroidz.str() + "; \n");
	}
	
	// write material data
	for (auto& phys : casePhysics)
	{
		for (auto& domain : phys.getDomainsList())
		{
			if (!domain->isPointDomain())
			{
				Material mat = phys.getMaterialById(domain->getMaterialId());
				GetdpMaterialWriter writer(mat, domain->getName());
				auto ph = this->getPhysics(domain->getSolverKey());
				writer.setHarmonic((ph.getTimeType() == Physics::HARMONIC or ph.getTimeType() == Physics::HARMONIC_WITH_TIME));
				outString.append(writer.writeData());
			}
		}

		for (auto& bnd : phys.getBoundariesList())
		{
			if (bnd->hasMaterial())
			{
				Material mat = phys.getMaterialById(bnd->getMaterialId());
				GetdpMaterialWriter writer(mat, bnd->getName());
				auto ph = this->getPhysics(bnd->getSolverKey());
				writer.setHarmonic((ph.getTimeType() == Physics::HARMONIC or ph.getTimeType() == Physics::HARMONIC_WITH_TIME));
				outString.append(writer.writeData());
			}
		}

		for (auto& bc : phys.getBoundariesList())
		{
			if ( (bc->getTypeKey() == "radiationBoundary") or (bc->getTypeKey() == "sibc"))
			{
				for (auto& domain : phys.getDomainsList())
				{
					if (domain->getId() == bc->getDomainId())
					{
						Material mat = phys.getMaterialById(domain->getMaterialId());
						GetdpMaterialWriter writer(mat, bc->getName());
						auto ph = this->getPhysics(bc->getSolverKey());
						writer.setHarmonic((ph.getTimeType() == Physics::HARMONIC or ph.getTimeType() == Physics::HARMONIC_WITH_TIME));
						outString.append(writer.writeData());
					}
				}

			}
		}
	}


	// write flags
	for (auto& phys : casePhysics)
	{
		for (auto fl : phys.getFlags())
		{
			outString.append("  " + fl.getName() + " = " + fl.getStringValue() + "; \n");
			if (fl.getStringValue() == "1")
			{
				for (auto& phys : casePhysics)
				{
					if (phys.getSolverKey() == fl.getSolverKey())
					{
						//TODO
						//phys.nonlinearityList.at(i) = true;
						break;
					}
				}
			}
		}
	}


	//write moved domains transformations
	for (auto mt : caseMotions)
	{
		if (mt.getType() == Motion::COMPLEX)
		{
			auto linear_velocity = mt.getVelocity();
			double rotation_angle = -mt.getAngularVelocity() * current_timestep;
			for (auto dom_nm : mt.getDomainNames())
			{
				auto dom = geom_data->getDomainByName(dom_nm);
				outString.append("  new_xyz[" + dom->getName() + "] = " +
					"Vector[X[] * Cos[" + std::to_string(rotation_angle) + "] - Y[] *Sin[" + std::to_string(rotation_angle) + "] - " 
					+ std::to_string(linear_velocity[0] * current_timestep) + ", "
					+ "X[] * Sin[" + std::to_string(rotation_angle) + "] + Y[] *Cos[" + std::to_string(rotation_angle) + "] - " 
					+ std::to_string(linear_velocity[1] * current_timestep) + " , "
					+ "Z[] - "
					+ std::to_string(linear_velocity[2] * current_timestep) + "];\n");
			}
		}
	}

	// TO DO
	// Make flexible coupling.
	// currently coupling in thermal equation is done if Flag_MagTherm is 2 or 3.
	if (casePhysics.size()>1)
	{
		if (casePhysics[0].getBasisKey() == "mesh2D")
			outString.append("  Flag_MagTherm = 2;\n");
		if (casePhysics[0].getBasisKey() == "mesh3D")
			outString.append("  Flag_MagTherm = 3;\n");
	}
	else
	{
		outString.append("  Flag_MagTherm = 0;\n");
	}

	if (casePhysics[0].getBasisKey() == "mesh2D")
		outString.append("  meshDim = 2;\n");
	if (casePhysics[0].getBasisKey() == "mesh3D")
		outString.append("  meshDim = 3;\n");


	//TODO add real breakpoints
	outString.append("  Breakpoints = {};\n");


	outString.append("}\n");
	functionFile << outString;
	functionFile.close();
}


void GetdpSetup::WriteConstraint()
{

	// Change initialization if restart
	for (auto& phys: casePhysics)
	{
		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			if (phys.getTimeParameters().count("continueFromLast"))
			{
				if (std::stof(phys.getTimeParameters()["continueFromLast"].getValue()) == 1)
				{
					phys.changeInitializationForRestart();
				}
			}
		}
	}

	std::ofstream constrFile;
	constrFile.open(wDir + "/calculation/constraint.pro");
	std::string outString = "Constraint { \n";

	for (auto& constr : constraints)
	{
		outString.append("{ Name " + constr.name + "; \n");
		outString.append("    Case { \n");
		for (auto phys : casePhysics)
		{
			std::string k = constr.region;
			if (constr.region.length()<1)
			{

				for (auto& boundary : phys.getBoundariesList()) 
				{

					for (auto& prop : boundary->getProperties())
					{

						if (prop.key == constr.name)
						{
							std::string bName = boundary->getName();
							if (Bcs::boundaryTypeNeedsCohomology(boundary->getTypeKey()))
								bName.append("_CH");
							if (prop.value.getTypeId() == InputValue::TABLE)
							{
								hasTransientConstraints = true;
								outString.append("       { Region " + bName + "; Type Assign;  Value 1; TimeFunction " + prop.value.getValue() + "; } \n");
							}
							else if (constr.assignFromResolution.length() > 0)
							{
								outString.append("       { Region " + boundary->getName() + "; Type AssignFromResolution; NameOfResolution " + constr.assignFromResolution + "; } \n");
							}
							else
							{
								outString.append("       { Region " + bName + "; Type Assign;  Value " + prop.value.getValue() + "; } \n");
							}

						}
					}
				}

				if ((phys.getTimeType() == Physics::TRANSIENT)) 
				{
					for (auto& domain : phys.getDomainsList()) 
					{
						if (!domain->isEnabled())
							continue;

						for (auto& ic : domain->getInitialConditions())
						{
							if ((ic.first == constr.name) and (domain->getSolverKey() == constr.solverKey))
							{
								if (ic.second.getValue() == "InitFromResolution")
								{
									hasTransientConstraints = true;
									outString.append("       { Region " + domain->getName() + "; Type InitFromResolution;  NameOfResolution  Initialization; } \n");
								}
								else
								{
									outString.append("       { Region " + domain->getName() + "; Type Init;  Value " + ic.second.getValue() + "; } \n");
								}
							}
						}
					}
				}

				for (auto& domain : phys.getDomainsList()) 
				{
					if (!domain->isEnabled())
						continue;

					for (auto& prop : domain->getProperties())
					{
						//check if the constraint is limited to this type of domains
						if ((constr.limitedTo.length() > 0) and (constr.limitedTo != domain->getTypeKey()) )
							continue;

						if (prop.first == constr.name) {
							if (prop.second.getTypeId() == InputValue::TABLE)
							{
								hasTransientConstraints = true;
								outString.append("       { Region " + domain->getName() + "; Type Assign;  Value 1; TimeFunction " + prop.second.getValue() + "; } \n");
							}
							else
							{
								outString.append("       { Region " + domain->getName() + "; Type Assign;  Value " + prop.second.getValue() + "; } \n");
							}
						}
					}
				}
			}
			else
			{
				std::string subRegString;
				subRegString = "       { Region " + constr.region + "; SubRegion Region[{";
				int k = 0;
				for (auto& boundary : phys.getBoundariesList()) 
				{
					for (auto& prop : boundary->getProperties())
					{
						if (prop.key == constr.name) {
							if (k != 0) subRegString.append(", ");
							subRegString.append(boundary->getName());
							k++;
						}
					}
				}
				subRegString.append("}]; Type Assign;  Value 0; } \n");
				if (k != 0)
					outString.append(subRegString);
			}
		}
		outString.append("    } \n");
		outString.append(" } \n");
	}

	for (auto& network : networks)
	{
		std::map<std::string, std::vector<Domains_> > networkDomains; //group; vector of domains
		for (auto& phys : casePhysics)
		{
			for (auto& dom : phys.getDomainsList()) {
				for (auto& src : dom->getProperties())
				{
					if (src.first == network.getName())
					{
						std::string grName;

						if (dom->getGroupKey() == "none")
						{
							grName = dom->getName();
							dom->setGroupKey(grName);
						}
						else
						{
							grName = dom->getGroupKey();
						}
						networkDomains[grName].push_back(dom);
					}
				}
			}
		}


		if (networkDomains.size() > 0)
			outString.append("{ Name " + network.getName() + "; Type Network ; \n");
		int crc = 0;
		for (const auto& nd : networkDomains)
		{
			outString.append("    Case Circuit" + std::to_string(crc) + " {\n");

			int dmint = 1;
			for (auto dm : nd.second)
			{
				if (dmint == 1)
				{
					outString.append("        { Region " + dm->getName() + ";		Branch{" + std::to_string(nd.second.size()) + "," + std::to_string(dmint) + "};}\n");
					dmint++;
				}
				else
				{
					outString.append("        { Region " + dm->getName() + ";		Branch{" + std::to_string(dmint - 1) + "," + std::to_string(dmint) + "};}\n");
					dmint++;
				}
			}
			crc++;
			outString.append("   }\n");
		}
		outString.append(" }\n");
	}
		

	outString.append("}\n");
	constrFile << outString;
	constrFile.close();
}


void GetdpSetup::writeResolution()
{

	auto coupling = getCouplingType();

	LOG(INFO) << "coupling type " << coupling;
	misc::createFolder(wDir + "\\getdpresults\\");
	std::ofstream resolFile;
	resolFile.open(wDir + "/calculation/resolution.pro");
	std::string outString = "";


	bool startTimeNotZero = false;
	float lastTimeStep = 0;
	for (auto& phys : casePhysics)
	{
		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			//if start time greater than zero, set startTimeNotZero = true
			if (std::stof(phys.getTimeParameters()["tstart"].getValue()) > 0)
				startTimeNotZero = true;

			//check if "continueFromLast" exists 
			if (phys.getTimeParameters().count("continueFromLast"))
			{
				if (std::stof(phys.getTimeParameters()["continueFromLast"].getValue()) == 1)
				{
					startTimeNotZero = true;

					 vtkForCenos vtuG;
					 vtuG.setOutput(wDir + "/results", "resFile");
					 lastTimeStep = vtuG.getLastTimeInSet(); 

					if (std::stof(phys.getTimeParameters()["tend"].getValue()) <= lastTimeStep)
					{
						resolFile << outString;
						resolFile.close();
						throw cenos_exception("End time smaller or equal to restart time.");
						return;
					}

					// add trestart time parameter
					// adding this may cause resolution to fail if only post processing is done!
					phys.setTimeParameter("trestart", lastTimeStep);
				}
			}

		}
	}


	// prepare lists for resolutions
	std::vector<std::string> solverKeyList;
	std::vector<std::map<std::string, InputValue>> timeParsList;
	std::vector<std::map<std::string, InputValue>> solverParsList;
	std::vector<Physics::timeType> timeTypeIdList;
	std::vector<std::string> preResolKeys; // preResolution keys
	std::vector<Physics::timeType> preResolTypeIdList;
	std::vector<bool> nonlinearityList;
	std::vector<std::shared_ptr<PostValues>> enabledPostValueList;
	std::vector<std::string> initResolutions;
	for (auto& phys : casePhysics)
	{
		solverKeyList.push_back(phys.getSolverKey());
		timeTypeIdList.push_back(phys.getTimeType());
		timeParsList.push_back(phys.getTimeParameters());
		solverParsList.push_back(phys.getSolverParameters());
		nonlinearityList.push_back(phys.hasNonlinearProps());
		auto pv = phys.getEnabledPostValueList();
		enabledPostValueList.insert(enabledPostValueList.begin(), pv.begin(), pv.end());

		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			if ( (startTimeNotZero) and (phys.getSolverKey() == "getdpThermal") )
				initResolutions.push_back(phys.getSolverKey());
		}

		// add stranded sigma formulation
		for (auto domain : phys.getDomainsList())
		{
			if ( (domain->getTypeKey() == "stranded_source_dom") and (phys.getBasisKey() == "mesh3D") )
			{
				preResolTypeIdList.push_back(Physics::STEADY); 
				preResolKeys.push_back("stranded_sigma");
				break;
			}
		}
	}
	outString.append(resolutionHeader(solverKeyList, timeTypeIdList, initResolutions, preResolKeys, preResolTypeIdList));
	outString.append(initilizationOperations(timeParsList, timeTypeIdList, solverKeyList, preResolKeys));
	outString.append(timeLoopOpen(timeParsList, timeTypeIdList, solverKeyList));

	//************** Calculation *************
	outString.append(outerLoopAccurateOpen(nonlinearityList, coupling, caseNumerics.getMaxIterations(), caseNumerics.getRelaxationFactor(), caseNumerics.getTolerance()));
		outString.append(solutionLoops(nonlinearityList, coupling, timeParsList, 
			caseNumerics.getMaxIterations(), caseNumerics.getRelaxationFactor(), caseNumerics.getTolerance(), 
			enabledPostValueList, solverKeyList, solverParsList, hasTransientConstraints));
	outString.append(outerLoopAccurateClose(nonlinearityList, coupling));

	//************** END Calculation *************

	outString.append(timeAdaptiveSplit(timeParsList, timeTypeIdList));
	outString.append(finalizeTimeStep(timeParsList, solverKeyList, timeTypeIdList));

	outString.append(timeLoopClose(timeTypeIdList, solverKeyList));
	outString.append(closeOperations());
	outString.append(closeResolution());

	resolFile << outString;
	resolFile.close();
}

void GetdpSetup::WritePostProc()
{

	std::ofstream postprocFile;
	postprocFile.open(wDir + "/calculation/postproc.pro");
	// postprocFile.open(wDir + "/postprocRunTime.pro");

	std::string outString;
	for (auto phys : casePhysics)
	{
		std::string slvKey = phys.getSolverKey();
		outString.append("PostProcessing { \n");
		std::string multiPort = "";
		if (slvKey == "getdpMicrowaves") {
			outString.append("For i In {1:#Sur_UNIFPORT_N_N_List()}\n");
			multiPort = "~{i}";
		}

		outString.append("  { Name postProc_" + slvKey + multiPort + "; NameOfFormulation " + slvKey + "_formulation" + multiPort + "; NameOfSystem resol_" + slvKey + multiPort + "; \n");

		outString.append("    PostQuantity { \n");

		for (auto& postValue : phys.getEnabledPostValueList()) {
			if (postValue->getSolverKey() == slvKey)
			{
				if (postValue->needLoop()) {
					if (postValue->getStringId() == "acS~{j}~{k}") {
						outString.append("	  For k In {1:#Sur_UNIFPORT_N_N_List()}\n");	
					}
					outString.append("	  For j In {1:#Sur_UNIFPORT_N_N_List()}\n");
					outString.append(getPostProcessing(postValue));

					if (postValue->getStringId() == "acS~{j}~{k}") {
						outString.append("	     EndFor \n");
					}
					outString.append("	  EndFor \n");
				}
				else {
					outString.append(getPostProcessing(postValue));
				}

			}
		}

		outString.append("    } \n");
		outString.append("  } \n");
		if (slvKey == "getdpMicrowaves") {
			outString.append("EndFor\n");
		}
		outString.append("} \n");
	}

	postprocFile << outString;
	postprocFile.close();
}


void GetdpSetup::WritePostOp(PostType postType)
{
	std::ofstream postOpFile;
	//std::ofstream runTimePostOpFile;
	if (postType == PostType::SOLUTION)
		postOpFile.open(wDir + "/calculation/postop_SOLUTION.pro");
	else if (postType == PostType::ENDTIME)
		postOpFile.open(wDir + "/calculation/postop_ENDTIME.pro");
	else if (postType == PostType::RUNTIME)
		postOpFile.open(wDir + "/calculation/postop_RUNTIME.pro");

	std::string outString;

	//std::filesystem::remove_all(wDir + "\\getdpresults\\");
	misc::createFolder(wDir + "\\getdpresults\\");
	misc::createFolder(wDir + "\\getdpresults\\globals\\");

	for (auto phys : casePhysics)
	{
		
		std::string slvKey = phys.getSolverKey();
		std::string multPort = "";
		if (slvKey == "getdpMicrowaves") {
			multPort = "~{i}";
		}
		// Scaling part
		if (postType == PostType::SOLUTION)
		{
			outString.append("PostOperation{ ");
			if (slvKey == "getdpMicrowaves") {
				multPort = "~{i}";
				outString.append("\nFor i In {1:#Sur_UNIFPORT_N_N_List()}\n");
			}
			outString.append("  {Name scalingPostOp_" + slvKey + multPort + "; NameOfPostProcessing postProc_" + slvKey + multPort + "; \n");
			outString.append("    Operation{ \n");

			for (auto& postValue : phys.getEnabledPostValueList())
			{
				if (postValue->getSolverKey() == slvKey)
				{
					outString.append(getPostOperations(postValue, "", 0.0, PostValues::SCALING_VALUE));
				}
			}

			outString.append("    } \n");
			outString.append("  } \n");
			if (slvKey == "getdpMicrowaves") {

				outString.append(" EndFor\n");
			}
			outString.append("} \n");
		}

		// in transient cases, transient cases have zeroth step initalization, other physics do not have this step
		bool initializationSkipped = false;
		if ( (hasTransientPhysics()) and (phys.getTimeType() != Physics::TRANSIENT) )
		{
			initializationSkipped = true;
		}

		outString.append("PostOperation{ ");
		if (slvKey == "getdpMicrowaves") {
			multPort = "~{i}";
			outString.append("\nFor i In {1:#Sur_UNIFPORT_N_N_List()}\n");
		}

		outString.append("  {Name postOp_" + slvKey + multPort + "; NameOfPostProcessing postProc_" + slvKey + multPort + "; \n");
		outString.append("    Operation{ \n");




		if (postType == PostType::SOLUTION)
		{
			for (auto& postValue : phys.getEnabledPostValueList())
			{
				if (postValue->getType() == PostValues::RESIDUAL_FOR_ADAPTIVE)
					continue;  
				if (postValue->getType() == PostValues::SCALING_VALUE)
					continue;

				outString.append(getPostOperations(postValue, "", 0.0));
			}
		}

		if (postType == PostType::ENDTIME)
		{
			std::vector<double> timeVec = getTimeList(postType);
			//to avoid repeating times step, usefull for multiport
			std::string portNumber = phys.getLastPort();
			
			for (size_t ts = 0; ts < timeVec.size(); ts++)
			{
				
				bool lastStep = false;
				if ( (ts == timeVec.size() - 1) and (timeVec.size() > 1) )
					lastStep = true;

				for (auto& postValue : phys.getEnabledPostValueList())
				{
					if (postValue->getType() == PostValues::RESIDUAL_FOR_ADAPTIVE)
						continue;

					if (postValue->getSolverKey() == slvKey)
					{

						// skip if it is last step for non-transient phys in transient case
						if ((postValue->getType() == PostValues::TOTAL_VALUE) or
							(postValue->getType() == PostValues::INTEGRAL) or
							(postValue->getType() == PostValues::BOUNDARY_INTEGRAL) or
							(postValue->getType() == PostValues::GLOBAL) or
							(postValue->getType() == PostValues::BOUNDARY_GLOBAL))
						{
							if (lastStep and initializationSkipped)
								continue;
						}



						if (postValue->getType() == PostValues::TOTAL_VALUE)
						{

							std::string domList = "Region[{";
							int g = 0;
							for (auto& domain : phys.getDomainsList())
							{
								if ((domain->getSolverKey() == postValue->getSolverKey()) and (postValue->isDefinedInDomain(getDomainGroup(domain))))
								{
									if (g != 0)
										domList.append(", ");
									domList.append(domain->getName());
									g++;
								}
							}
							domList.append("}]");
							outString.append(getPostOperations(postValue, domList, ts));
							continue;
						}

						for (auto& domain : phys.getDomainsList())
						{
							std::string k = domain->getSolverKey();
							std::string m = domain->getTypeKey();
							if ((domain->getSolverKey() == postValue->getSolverKey()) and (postValue->isDefinedInDomain(getDomainGroup(domain))))
							{
								outString.append(getPostOperations(postValue, domain->getName(), ts));
							}
						}
						//write only one postop line for multiple port loadVoltage
						bool portWritten = false;

						for (auto& bc : phys.getBoundariesList())
						{
							//NEEDED FOR MULTIPLE PORt
						std:string multiPort_j, multiPort_i = "";
							

							if (!getBoundaryGroup(bc).find("Sur_UNIFPORT_N_")) {
								//bc.boundaryGroup = "Sur_UNIFPORT_N~{i}";
								multiPort_j = "Sur_UNIFPORT_N~{j}";
								multiPort_i = "Sur_UNIFPORT_N~{i}";

							}

							if ((bc->getSolverKey() == postValue->getSolverKey()) and (postValue->isDefinedInDomain(multiPort_j)))
							{

								if (!portWritten && postValue->needLoop()) {
									if (postValue->getStringId() == "acS~{j}~{k}") {
										outString.append("      For k In {1:#Sur_UNIFPORT_N_N_List()}\n");
									}
									outString.append("    For j In {1:#Sur_UNIFPORT_N_N_List()}\n");

									outString.append(getPostOperations(postValue, bc->getName(), ts));
									if (postValue->getStringId() == "acS~{j}~{k}") {
										outString.append("      EndFor\n");
									}
									outString.append("    EndFor\n");
									portWritten = true;
								}

							}
							if ((bc->getSolverKey() == postValue->getSolverKey()) and (postValue->isDefinedInDomain(multiPort_i)))
							{
								outString.append(getPostOperations(postValue, bc->getName(), ts));
							}

							if ((bc->getSolverKey() == postValue->getSolverKey()) and (postValue->isDefinedInDomain(getBoundaryGroup(bc))))
							{
								outString.append(getPostOperations(postValue, bc->getName(), ts));
							}

						}
					}
				}
			}
		}

		/*		for (auto prb : phys.probeList)
				{
					outString.append(prb.writePostOperationForGetdp());
				}*/
		
		outString.append("    } \n");
		
		outString.append("  } \n");
		if (slvKey == "getdpMicrowaves") {

			outString.append(" EndFor\n");
		}
		outString.append("} \n");
		

		// Adaptive part
		if (postType == PostType::SOLUTION)
		{
			outString.append("PostOperation{ ");
			std::string multPort = "";
			if (slvKey == "getdpMicrowaves") {
				multPort = "~{i}";
				outString.append("\nFor i In {1:#Sur_UNIFPORT_N_N_List()}\n");
			}
			outString.append("  {Name adaptivePostOp_" + slvKey + multPort + "; NameOfPostProcessing postProc_" + slvKey + multPort + "; \n");
			outString.append("    Operation{ \n");


			for (auto& postValue : phys.getEnabledPostValueList())
			{
				if (postValue->getSolverKey() == slvKey)
				{
					outString.append(getPostOperations(postValue, "", 0.0, PostValues::RESIDUAL_FOR_ADAPTIVE));
				}
			}

			outString.append("    } \n");
			outString.append("  } \n");
			if (slvKey == "getdpMicrowaves") {

				outString.append(" EndFor\n");
			}
			outString.append("} \n");

		}
	}
	postOpFile << outString;

	postOpFile.close();
	//runTimePostOpFile.close();
}


std::vector<double> GetdpSetup::getTimeList(PostType postType)
{
	double endTime = 0;
	int precision = 7;

	for (auto phys: casePhysics)
	{
		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			endTime =  std::stof(phys.getTimeParameters()["tend"].getValue());
		}
	}

	std::vector<double> retVec;

	std::ifstream infile(wDir + "/getdpresults/times.dat");
	double t;
	std::string stringRead;
	

	infile.seekg(0, std::ios::beg);

	//for transient
	std::vector<size_t> read_times;
	while (infile >> stringRead)
	{
		if (stringRead.front() == '(')
		{
			size_t found_comma =  stringRead.find(",");
			std::string complexVal = stringRead.substr(1, found_comma);
			t = std::stod(complexVal);
		}
		else
		{
			t = std::stod(stringRead);
		}

		// if t is not endtime and post type is not runtime post
		if (!((t == endTime) and (postType == PostType::RUNTIME))) {
			//if time is repeated then this time is skipped, used for multiport
			if (!(std::find(read_times.begin(), read_times.end(), t) != read_times.end())) {
				retVec.push_back(round(t * pow(10.0, precision)) / pow(10.0, precision));
			}
		}
		read_times.push_back(t);
		
			
	}

	infile.close();
	return retVec;
}


void GetdpSetup::addNetworkDomain()
{
	for (auto& network : networks)
	{
		std::map<std::string, std::vector<Domains_> > networkDomains; //group; vector of domains
		for (auto& phys : casePhysics)
		{
			for (auto& dom : phys.getDomainsList()) {
				for (auto& src : dom->getProperties())
				{
					if (src.first == network.getName())
					{
						std::string grName;

						if (dom->getGroupKey() == "none")
						{
							grName = dom->getName();
							dom->setGroupKey(grName);
						}
						else
						{
							grName = dom->getGroupKey();
						}
						networkDomains[grName].push_back(dom);
					}
				}
			}
		}

		int networkCount = 0;
		for (auto nd : networkDomains)
		{
			auto ents = mesh->getEntities();
			auto ent_maxid = std::max_element(ents.begin(), ents.end(),
				[](std::shared_ptr<MeshEntity> const lhs, std::shared_ptr<MeshEntity> const rhs) {
					return lhs->getId() < rhs->getId(); });
			int domid = (*ent_maxid)->getId() + 1;

			std::string domname = network.getName() + "_source" + std::to_string(networkCount);

			for (auto& phys : casePhysics)
			{
				if (phys.getSolverKey() == network.getSolverKey())
				{
					json domJson;
					domJson["id"] = domid;
					domJson["typeKey"] = network.getTypeKey();
					domJson["groupKey"] = nd.first;
					domJson["enabled"] = true;
					domJson["isPointDomain"] = true;
					Domains_ dom = std::make_shared<Domains>(domname, domJson, network.getSolverKey());
					dom->setProperties(nd.second[0]->getProperties());
					phys.addDomain(dom);
					// make entity with dummy shape (pnt)
					std::shared_ptr<TopoEntity> ent = std::make_shared<TopoEntity>(TopoDS_Shape(BRepBuilderAPI_MakeVertex(gp_Pnt(0,0, networkCount))), domname, domid);
					ent->setEnabled(false);
					dom->addEntity(ent);
					geom_data->addEntity(ent);


					std::shared_ptr<DomainGroup> dgroup = std::make_shared<DomainGroup>();
					dgroup->setName(domname);
					dgroup->setRole("_dummy_");
					dgroup->setId(domid);
					dgroup->addEntity(ent);
					dgroup->update();
					geom_data->addDomainGroup(dgroup);


					std::shared_ptr<MeshEntity> mesh_ent = std::make_shared<MeshEntity>(TopoDS_Shape(BRepBuilderAPI_MakeVertex(gp_Pnt(0,0, networkCount))), domname, domid);
					mesh_ent->setEnabled(false);
					mesh->addEntity(mesh_ent);

					networkCount++;

					json dtJson;
					dtJson["key"] = network.getTypeKey();
					Domains::addDomainTypes(dtJson);
				}

			}

		}
	}
}

//used for postprocessing to merge results with previous or write new
bool GetdpSetup::hasRestart()
{
	for (auto& phys: casePhysics)
	{
		for (auto& tp : phys.getTimeParameters())
		{
			if (tp.first == "continueFromLast")
			{
				if (tp.second.getValue() == "1")
					return true;

			}
			// For AD this needs to be false for all time steps
			if (phys.getSolverKey() == "getdpMicrowaves")
			{
				return false;
			}
		}
	}

	// if there is no continueFromLast, restart can also happen in scanning cases
	if (current_time > 0)
		return true;

	bool has_complex_motion = false;
	for (auto mt : caseMotions)
	{
		// skip if disabled
		if (!mt.isEnabled())
			continue;

		//skip if simple motion
		if (mt.getType() == Motion::SIMPLE)
			continue;

		has_complex_motion = true;
	}

	// first time step in scanning
	if ( !initial_step and has_complex_motion)
		return true;

	return false;
}

void GetdpSetup::setFlags()
{
	// set boundary flag values
	for (auto& phys : casePhysics)
	{
		for (auto& boundary : phys.getBoundariesList())
		{
			for (auto& boundaryType : Bcs::existingBoundaryTypes)
			{
				if (boundaryType.key == boundary->getTypeKey())
				{
					for (auto& prop : boundaryType.properties)
					{
						if (prop.front() == '_') 
						{
							if (boundary->getBcPropertyValue(prop) == "1")
							{
								phys.setFlagValue(prop, true);
							}
						}
					}
				}
			}
		}
	}


	//set material flag values
	for (auto& phys : casePhysics)
	{
		for (auto& material : phys.getMaterialsList())
		{
			for (auto& matData : material.getMaterialData())
			{
				if (matData.first == "BHModel")
				{
					phys.setFlagValue(matData.first, true);
				}
				else if (!matData.second.isBoolean())
				{
					if (matData.second.isTable())
					{
						phys.setFlagValue(matData.first, true);
					}
				}
			}
		}
	}				   				 
}

// Here we change boundary conditions.
// For solid conductors (PEC or normal) all the boundaries that are interface are changed to PEC or normal conductors
void GetdpSetup::updateBoundaryConditions()
{
	for (auto& phys : casePhysics)
	{
		//change type for non-set boundaries of conductor and pec domains
		if (phys.getPhysicsKey() == "physicsMicrowaves")
		{
			for (auto& dom : phys.getDomainsList())
			{
				//perfect electric conductor or (imperfect) condutor domain
				if (dom->getTypeKey() == "em_pec" or dom->getTypeKey() == "em_cc")
				{
					std::shared_ptr<DomainGroup> domain_group = geom_data->getDomainByName(dom->getName());

					//check all boundary groups belonging to domain_group
					for (auto bnd_group : geom_data->getBoundaries())
					{
						if (bnd_group->isSubShapeOf(*domain_group))
						{
							//find bnd corresponding to bnd_group and reassign boundary type
							auto bnd = phys.getBoundaryByName(bnd_group->getName());

							// skip if boundary type is assigned (interface is default)
							if (bnd->getTypeKey() != "interface")
								continue;

							// if domain is set as PEC, assign PEC to boundary.
							if (dom->getTypeKey() == "em_pec")
							{
								bnd->typeKey = "perfectConductor";
								bnd->addBcProperty("Efield", InputValue(0.0));
							}

							// if domain is set as conductor (imperfect), assign conductor to boundary.
							else if (dom->getTypeKey() == "em_cc")
							{
								bnd->typeKey = "conductingBoundary";
								bnd->has_material = true;
								bnd->material = phys.getMaterialById(dom->getMaterialId());
								bnd->material_id = dom->getMaterialId();
							}

						}
					}

				}
			}
		}
	}
}

void GetdpSetup::updateMeshEntities()
{
	// Assign shapes to MeshEntities (they have no shapes at read time for msh format)
	// Equalize ids for mesh and geometryData.
	for (auto ent : mesh->getEntities())
	{
		// CP-1263 Workaround for lastMesh.msh files that were written with 3.0.0 version (it contains CH entities)
		// in further versions this file will not contain these disabled entities
		std::string name = ent->getName();
		if (name.length() > 3 and name.substr(name.length() - 3) == "_CH")
			continue;




		try
		{
			auto gd_ent = geom_data->getEntityByName(ent->getName());
			ent->setShape(gd_ent->getShape());
			ent->setId(gd_ent->getId());
		}
		catch (std::exception& e)
		{
			// check for CP-1314
			// relevant only for conversions from 2.8.2. 
			// In 3.0.0 and 3.0.1 case conversion did not create new dump files, so mesh was ok,
			// but in scanning old dump files were used and that lead to inconsistency.
			throw(cenos_exception("Mesh inconsistend with geometry data. Try sending mesh again from Salome or contact support!"));
		}

	}
}

json GetdpSetup::getBlockGroupsJson() {
	json blockGroups = json::object();

	for (auto dom : casePhysics[0].getDomainsList()) {
		if (!dom->isPointDomain()) {
			std::string group = dom->getGroupKey();
			if (group != "none") {
				if (!blockGroups.contains(group)) {
					blockGroups[group]["blocks"] = json::array();
					blockGroups[group]["role"] = dom->getRole();
				}
				blockGroups[group]["blocks"].push_back(dom->getLabel());
				
			}
		}
	}
	return blockGroups;
}

json getBlocksJson(GeometryData_ geom_data) {
	json blocks = json::array();

	for (DomainGroup_ domain : geom_data->getDomains()) {
		json block = json::object();
		//block["name"] = domain->getName();
		block["label"] = domain->getLabel();
		block["role"] = domain->getRole();
		block["type"] = "domain";
		blocks.push_back(block);
	}
	for (BoundaryGroup_ boundary : geom_data->getBoundaries()) {
		json block = json::object();
		block["name"] = boundary->getName();
		block["label"] = boundary->getLabel();
		block["role"] = boundary->getRole();
		block["type"] = "boundary";
		block["group"] = nullptr;
		blocks.push_back(block);
	}

	return blocks;
}

json getFieldFormattingJson(std::vector<Physics> casePhysics) {
	json fieldFormatting = json::array();
	for (Physics phys : casePhysics) {

		for (std::shared_ptr<PostValues> post : phys.getEnabledPostValueList()) {
			if (!post->isIntermediate() and post->getType() == PostValues::NODE_FIELD) {
				json postValue = json::object();
				postValue["key"] = post->getStringId();
				postValue["name"] = post->getLabel();
				postValue["safeName"] = post->getName();
				postValue["symbol"] = post->getSymbol();
				postValue["units"] = post->getUnits();
				fieldFormatting.push_back(postValue);
			}

		}
	}
	return fieldFormatting;
}

json GetdpSetup::getGlobalPostValuesJson(std::vector<Physics> casePhysics, std::string wDir, bool isLegacyCase) {
	// This function a member of the GetdpSetup class only because of the getBoundaryGroup and getDomainGroup functions.

	json postValueArray = json::object();
	
	std::set<std::shared_ptr<PostValues>> allPostValues;
	std::set<Domains_> allDomains;
	std::set<Bcs_> allBoundaries;
	// Aggregate all postValues, domains and boundaries, we will check if they contain anything later
	for (Physics phys : casePhysics) {
		for (Domains_ domain : phys.getDomainsList()) {
			allDomains.insert(domain);
		}

		for (Bcs_ boundary : phys.getBoundariesList()) {
			allBoundaries.insert(boundary);
		}

		for (std::shared_ptr<PostValues> post : phys.getEnabledPostValueList()) {
			if (post->isIntermediate())
				continue;

			switch (post->getType()) {
			case PostValues::GLOBAL:
			case PostValues::INTEGRAL:
			case PostValues::GLOBAL_EXTERNAL:
			case PostValues::BOUNDARY_GLOBAL:
			case PostValues::BOUNDARY_INTEGRAL:
			case PostValues::CIRCULATION:
				allPostValues.insert(post);
			}
		}
	}

	for (std::shared_ptr<PostValues> post : allPostValues) {
		json postValueObj = json::object();

		postValueObj["name"] = post->getLabel();
		postValueObj["units"] = post->getUnits();
		postValueObj["domainValues"] = json::array();

		std::vector<double> timeValues;

		// Use this to check if a postvalue has any domains or boundaries where it is defined.
		// Do not add it to the output if there are no values in any domain or boundary.
		bool isDefinedAnywhere = false;

		std::vector<std::pair<double, double>> valuesSum;

		for (Domains_ domain : allDomains) {

			if (post->isDefinedInDomain(getDomainGroup(domain))) {
				isDefinedAnywhere = true;
				json domainObj = json::object();

				std::vector<std::pair<double, double>> values;
				if (isLegacyCase) {
					std::string domainName;

					if (domain->getName().rfind("ElectricalCircuit", 0) == 0) {
						// Instead of "ElectricalCircuit_source0", "winding_2..." is needed for inductor terminals.
						domainName = domain->getGroupKey();
					}
					else {
						domainName = domain->getName();
					}
					values = processGlobalValuesV2Legacy(post, domainName, wDir);


					if (values.size() == 0) {
						// The values of this domain are not in the globalResultsFile.json
						// Skip this domain. The processGlobalValuesV2Legacy logs this event as a warnings.
						isDefinedAnywhere = false;
						continue;
					}
				}
				else {
					values = processGlobalValuesV2(post, domain->getName(), wDir);
				}

				if (post->getSumAllDomains()) {
					// Loop over domain values and add them to the total
					if (valuesSum.empty()) {
						valuesSum = values;
					}
					else {
						for (int i = 0; i < values.size(); i++) {
							double timeVal = valuesSum[i].first;
							if (timeVal == values[i].first) {
								double sumVal = valuesSum[i].second + values[i].second;
								valuesSum[i] = std::pair <double, double>(timeVal, sumVal);
							}
							else {
								throw cenos_exception("Time values in " + post->getStringId() + " do not match previous ones");
							}
						}
					}
				}

				if (domain->isPointDomain()) {
					// This is the special group domain with the "..."
					domainObj["blockName"] = domain->getGroupKey();
					domainObj["blockGroup"] = "none";
				}
				else {
					// This domain is a normal one
					domainObj["blockName"] = domain->getName();
					domainObj["blockGroup"] = domain->getGroupKey();
				}

				domainObj["values"] = json::array();


				if (timeValues.empty()) {
					// Populate time values for the first time
					for (std::pair val : values) {
						timeValues.push_back(val.first);
						domainObj["values"].push_back(val.second);
					}
				}
				else {
					// Time values already populated, check that they match up
					int i = 0;
					for (std::pair val : values) {
						if (!timeValues[i] == val.first) {
							throw cenos_exception("Time values in " + post->getStringId() + " do not match previous ones");
						}
						domainObj["values"].push_back(val.second);
						i++;
					}
				}

				postValueObj["domainValues"].push_back(domainObj);

			}
		}

		if (post->getSumAllDomains()) {
			json domainObjSum = json::object();
			domainObjSum["values"] = json::array();

			domainObjSum["blockName"] = "All_Domains";
			domainObjSum["blockGroup"] = "none";

			int i = 0;
			for (std::pair val : valuesSum) {
				if (!timeValues[i] == val.first) {
					throw cenos_exception("Time values in " + post->getStringId() + " do not match previous ones");
				}
				domainObjSum["values"].push_back(val.second);
				i++;
			}

			postValueObj["domainValues"].push_back(domainObjSum);
		}
		
		for (Bcs_ boundary : allBoundaries) {
			if (post->isDefinedInDomain(getBoundaryGroup(boundary))) {
				isDefinedAnywhere = true;
				json domainObj = json::object();

				domainObj["blockName"] = boundary->getName();
				domainObj["blockGroup"] = "none";


				std::vector<std::pair<double, double>> values;
				if (isLegacyCase) {
					values = processGlobalValuesV2Legacy(post, boundary->getName(), wDir);

					if (values.size() == 0) {
						// The values of this domain are not in the globalResultsFile.json
						// Skip this domain. The processGlobalValuesV2Legacy logs this event as a warnings.
						isDefinedAnywhere = false;
						continue;
					}
				}
				else {
					values = processGlobalValuesV2(post, boundary->getName(), wDir);
				}
				domainObj["values"] = json::array();

				if (timeValues.empty()) {
					// Populate time values for the first time
					for (std::pair val : values) {
						timeValues.push_back(val.first);
						domainObj["values"].push_back(val.second);
					}
				}
				else {
					// Time values already populated, check that they match up
					int i = 0;
					for (std::pair val : values) {
						if (!timeValues[i] == val.first) {
							throw cenos_exception("Time values in " + post->getStringId() + " do not match previous ones");
						}
						domainObj["values"].push_back(val.second);
						i++;
					}
				}

				postValueObj["domainValues"].push_back(domainObj);
			}
		}

		json timeValueObj;
		if (misc::getCenosApp() == "ANTENNAS") {
			// Because the frequencies are in MHz, convert them to Hz
			int k = 1000000; // 10^6
			std::transform(timeValues.begin(), timeValues.end(), timeValues.begin(), [k](double& c) { return c * k; });

			timeValueObj = {
				{"name", "Frequency"},
				{"units", "Hz"},
				{"values", timeValues},
			};
		}
		else {
			timeValueObj = {
				{"name", "Time"},
				{"units", "s"},
				{"values", timeValues},
			};
		}

		if (isDefinedAnywhere) {
			postValueArray[post->getStringId()]["x"] = timeValueObj;
			postValueArray[post->getStringId()]["y"] = postValueObj;
		}
	}
	return postValueArray;
}

void GetdpSetup::outputResultsJson(bool isLegacyCase) {
	// Write the results.json file to the results folder.
	// isLegacyCase is used only for converting cases before this file format was implementd.
	std::string fname = wDir + "/results/results.json";
	json resultsJson;

	if (hasRestart()) {
		// Continuing from last time step, have to append new results to the existing file.
		std::ifstream input(fname);
		if (!input.is_open()){
			throw cenos_exception(fname + " not found");
		}

		std::stringstream filebuff;
		filebuff << input.rdbuf();
		resultsJson = json::parse(filebuff.str());

		// blocks dont change over time, so postValues need to be updated.
		json newpostValues = GetdpSetup::getGlobalPostValuesJson(casePhysics, wDir, isLegacyCase);

		for (auto it : resultsJson["globalPostValues"].items()) {
			// Update x values
			// Get the values under the same key in the new json
			json newXValues = newpostValues[it.key()]["x"]["values"];
			json post = it.value();
			post["x"]["values"].insert(post["x"]["values"].end(), newXValues.begin(), newXValues.end());

			// Update all y values
			json domains = post["y"]["domainValues"];
			int i = 0;
			for (json dom : domains) {
				json newValues = newpostValues[it.key()]["y"]["domainValues"][i]["values"];

				dom["values"].insert(dom["values"].end(), newValues.begin(), newValues.end());
				post["y"]["domainValues"][i] = dom;
				i++;
			}



			resultsJson["globalPostValues"][it.key()] = post;
		}
	}
	else {
		// A new calculation was run, create a new file.
		resultsJson = json::object();
		resultsJson["fileVersion"] = "1.0";
		resultsJson["general"] = json::object();
		resultsJson["general"]["app"] = misc::getCenosApp();
		resultsJson["general"]["isAxisymmetric"] = (casePhysics[0].getSymmetryType() == Physics::symmetryType::AXISYMMETRIC_2D);
		resultsJson["general"]["dimension"] = geom_data->getDimension();

		resultsJson["blocks"] = getBlocksJson(geom_data);
		resultsJson["blockGroups"] = getBlockGroupsJson();
		resultsJson["globalPostValues"] = GetdpSetup::getGlobalPostValuesJson(casePhysics, wDir, isLegacyCase);
		resultsJson["fieldFormatting"] = getFieldFormattingJson(casePhysics);
	}
	
	std::ofstream file;
	file.open(fname);
	file << resultsJson.dump(4);
	file.close();
}


void GetdpSetup::outputGlobals()
{
	std::cout << geom_data->isAxisymmetric() << std::endl;
	if (hasRestart())
		readGlobalsJson();

	std::string fname = wDir + "/results/globalResultsFile.json";
	std::ofstream globalsResFileJson;
	globalsResFileJson.open(fname);

	std::string outString;
	json resultJson;
	json domainsArray = json::array();
	json bcsArray = json::array();
	for (auto phys : casePhysics)
	{

		for (auto& postValue : phys.getEnabledPostValueList())
		{
			if (postValue->isIntermediate())
				continue;

			if ( (postValue->getType() == PostValues::TOTAL_VALUE))
			{
				json domainObj = json::object();
				bool domainIsInArray = false;

				//get globalValues for this postValue
				domainObj = processGlobalValues(postValue, "", wDir,"");
				domainObj["name"] = "All_Domains";
				domainObj["group"] = "none";

				//if such domain exists, update it. 
				//else add it to domainsArray
				for (auto& dm : domainsArray)
				{
					if (dm["name"] == "All_Domains")
					{
						if (dm.find("postValues") == dm.end())
							dm["postValues"] = json::object();

						
						if (domainObj.find("postValues") != domainObj.end())
							dm["postValues"].push_back(domainObj["postValues"][0]);
						if (domainObj.find("time") != domainObj.end())
							dm["time"] = domainObj["time"];
						domainIsInArray = true;
						break;
					}
				}

				if ((!domainIsInArray) and (domainObj.find("postValues") != domainObj.end()))
					domainsArray.push_back(domainObj);
			}

			for (auto& domain : phys.getDomainsList())
			{
				if ((postValue->getSolverKey() == phys.getSolverKey()) and (postValue->isDefinedInDomain(getDomainGroup(domain)))
					and (postValue->getType() != PostValues::TOTAL_VALUE))
				{
					json domainObj = json::object();
					bool domainIsInArray = false;
					std::string domName;
					std::string domGroup;

					//set name and group name for domain
					if (domain->isPointDomain())
					{
						domName = domain->getGroupKey();
						domGroup = "none";
					}
					else
					{
						domName = domain->getName();
						domGroup = domain->getGroupKey();
					}

					//get globalValues for this postValue
					domainObj = processGlobalValues(postValue, domain->getName(), wDir ,"");
					domainObj["name"] = domName;
					domainObj["group"] = domGroup;
					for (auto dom : geom_data->getDomains())
					{
						if (dom->getName() == domName)
						{
							if (dom->getRole() != "")
								domainObj["role"] = dom->getRole();
						}
					}

					//if such domain exists, update it. 
					//else add it to domainsArray
					for (auto& dm : domainsArray)
					{
						if (dm["name"] == domName)
						{
							if (dm.find("postValues") == dm.end())
								dm["postValues"] = json::object();


							if (domainObj.find("postValues") != domainObj.end())
								dm["postValues"].push_back(domainObj["postValues"][0]);
							if (domainObj.find("time") != domainObj.end())
								dm["time"] = domainObj["time"];
							domainIsInArray = true;
							break;
						}
					}

					if ( (!domainIsInArray) and (domainObj.find("postValues") != domainObj.end()))
						domainsArray.push_back(domainObj);
				}
			}
		}
	}
	resultJson["domains"] = domainsArray;

	for (auto phys : casePhysics)
	{
		for (auto& bc : phys.getBoundariesList())
		{
			std::string portNumber = phys.getLastPort();
			for (auto& postValue : phys.getEnabledPostValueList())
			{
				if (postValue->isIntermediate())
					continue;

				json bcObj = json::object();
				bool boundaryIsInArray = false;

				if ((postValue->getSolverKey() == phys.getSolverKey()) and (postValue->isDefinedInDomain(getBoundaryGroup(bc))))
				{
					
					if (phys.getSolverKey() == "getdpMicrowaves") {
						bcObj = processGlobalValues(postValue, bc->getName(), wDir, portNumber);
					}
					else {
						bcObj = processGlobalValues(postValue, bc->getName(), wDir, "");
					}
					bcObj["name"] = bc->getName();
					for (auto bnd : geom_data->getBoundaries())
					{
						if (bnd->getName() == bc->getName())
						{
							if (bnd->getRole() != "")
								bcObj["role"] = bnd->getRole();
						}
					}
				}
				if ((postValue->getSolverKey() == phys.getSolverKey()) and (postValue->isDefinedInDomainMultiPort(getBoundaryGroup(bc))))
				{
					bcObj = processGlobalValues(postValue, bc->getName(), wDir, portNumber);
					bcObj["name"] = bc->getName();
					for (auto bnd : geom_data->getBoundaries())
					{
						if (bnd->getName() == bc->getName())
						{
							if (bnd->getRole() != "")
								bcObj["role"] = bnd->getRole();
						}
					}
				}



				for (auto& bnd : bcsArray)
				{

					if (bnd["name"] == bc->getName())
					{

						if (bcObj.find("postValues") != bcObj.end()){
							
							for (int i = 0; i < bcObj["postValues"].size(); i++) {
								if (postValue->getStringId() == "acS~{j}~{k}") {
									break;
								}
								bnd["postValues"].push_back(bcObj["postValues"][i]);
							}
						}
						if (bcObj.find("time") != bcObj.end())
							bnd["time"]= bcObj["time"];

						if (postValue->getStringId() == "acS~{j}~{k}") {
							size_t portSize = std::stod(phys.getLastPort());
							size_t scatSize = portSize * portSize;
							size_t timeSize = bcObj["time"]["timeValues"].size()/scatSize;
							size_t Sj = 1;
							size_t Si = 1;
							for (int i = 0; i < scatSize; i++) {
								json bcObjScat = bcObj;
								std::vector<double> scatVal;
								bcObjScat["postValues"][0]["postName"] = "S" + std::to_string(Sj) + std::to_string(Si);
								bcObjScat["postValues"][0]["postStringId"] = "S" + std::to_string(Sj) + std::to_string(Si);
								if (Si == portSize) {
									Si = 1; Sj++; 
								}
								else { 
									Si++;
								}
								for (int time = 0; time < timeSize; time++)
								{
									scatVal.push_back(bcObj["postValues"][0]["values"][i+ scatSize*time]);
								}
								bcObjScat["postValues"][0]["values"] = scatVal;
								bnd["postValues"].push_back(bcObjScat["postValues"][0]);
							}
						}

						boundaryIsInArray = true;
						break;
					}
				}

				if ( (!boundaryIsInArray) and (bcObj.find("postValues") != bcObj.end()))
					bcsArray.push_back(bcObj);

			}
		}
	}
	resultJson["bcs"] = bcsArray;

	if (!globalsJson.empty())
	{
		appendTimeToGlobalsJson(resultJson, globalsJson);
	}
	else
	{
		globalsJson = resultJson;
	}


	std::string s = globalsJson.dump();
	globalsResFileJson << globalsJson.dump();
	globalsResFileJson.close();
	std::string csvName = wDir + "/results/globalResultsFile.csv";
	std::ofstream globalsResFileCsv;
	globalsResFileCsv.open(csvName);
	json mergedJson = mergeGroupedDomainsInJson(globalsJson);
	std::string sa = mergedJson.dump();

	outString.append(prepareCsvFromJson(mergedJson));
	globalsResFileCsv << outString;
	globalsResFileCsv.close();
}

void GetdpSetup::initialize()
{
}

//TODO: return value is not used, so can be removed!
bool GetdpSetup::remesh()
{
	if (caseMotions.empty())
		return true;

	// do not remesh (use original) if this is initial step and current time is zero
	if (current_time <= 0 and initial_step)
		return true;


	// ******************************
	//        Generate new mesh
	// ******************************
	
	messageQueue->push(messagedata(messagedata::info_lvl, messagedata::log, "Remeshing."));

	GeometryMover gm;

	bool old_remesh = false;
	bool new_remesh = false;
	std::vector<std::string> moved_domains;
	for (auto mt : caseMotions)
	{
		// skip if disabled
		if (!mt.isEnabled())
			continue;

		//skip if simple motion
		if (mt.getType() == Motion::SIMPLE)
			continue;

		//if motion type is old motion, move mesh immediately.
		// UI is not supposed to allow mixed type of motion (old+ new)
		// this is checked in UI
		if (mt.getType() == Motion::COMPLEX_WITH_PARAMETERS)
		{
			old_remesh = true;
			break;
		}

		//only new motion reaches this point
		new_remesh = true;

		//assing transformation for each domain in motion object
		auto domain_names = mt.getDomainNames();
		for (auto domn : domain_names)
		{
			moved_domains.push_back(domn);
			auto dom = geom_data->getDomainByName(domn);
			if (dom->getRole() == "air")
				throw(cenos_exception("Motion error: domains with role AIR cannot be moved! Please check your input."));
			for (auto ent : dom->getEntities())
				if (ent->isEnabled())
					gm.setTransform(ent, mt.getTransform(current_timestep));
		}
	}

	if (new_remesh)
	{
		//assign zero transform for non-moved entities
		for (auto dom : geom_data->getDomains())
		{
			if (std::find(moved_domains.begin(), moved_domains.end(), dom->getName()) != moved_domains.end())
				continue;
			if (dom->getRole() == "air")
				continue;

			gp_Trsf transform;
			transform.SetTranslation(gp_Vec(0, 0, 0));
			for (auto ent : dom->getEntities())
				if (ent->isEnabled())
					gm.setTransform(ent, transform);
		}

		gm.wDir = wDir;
		// modifies geom_data object!!!
		gm.moveGeometry(geom_data);
	}

	// everything fine if no remesh is necessary ( SIMPLE motion)
	if (!new_remesh and !old_remesh)
		return true;

	auto mesh_scale = mesh->getMeshScale();
	if (new_remesh)
	{
		MeshGenerator mg(geom_data);
		mg.setWdir(wDir);
		mg.setMesh(mesh);
		for (auto ft : gm.getFullTransform())
			mg.setTransform(ft.first, ft.second);
		mesh = mg.generateMesh();
	}

	if (old_remesh)
	{
		Motion* old_motion;
		for (auto& mt : caseMotions)
		{
			// only one old motion is allowed
			if (mt.getType() == Motion::COMPLEX_WITH_PARAMETERS)
			{
				old_motion = &mt;
				break;
			}
		}
		SalomeModule salome;
		salome.setWDir(wDir);
		salome.setRunMode(SalomeModule::BATCH_NO_RESPONSE);
		old_motion->readOldParameters_OM(wDir + "/geometry/SalomeMotionParameters.json");
		old_motion->calculateParameters_OM(current_time, current_timestep);
		salome.writeDumpScript(old_motion->getParameters_OM());
		messagedata msg = salome.runSalome();

		if (msg.msgStatus == messagedata::error)
			throw(cenos_exception(msg.text_data));

		old_motion->saveUpdatedParameters_OM(wDir + "/geometry/SalomeMotionParameters.json");

		try
		{
			mesh->clear();
			mesh->readMesh(wDir + "/geometry/meshFile.msh", std::string("GMSH"));
		}
		catch (cenos_exception& e)
		{
			throw(cenos_exception("Could not open last mesh file " + wDir + "/geometry/lastMesh.msh. Cannot continue calculation!"));
		}


		updateMeshEntities();

	}

	// ******************************
	// Generate cohomology cuts for 3D current supply!
	// ******************************

	for (auto& phys : casePhysics)
	{
		for (auto& boundary : phys.getBoundariesList())
		{
			for (auto& boundaryType : Bcs::existingBoundaryTypes)
			{
				if ((boundaryType.needs_cohomology) and (boundaryType.key == boundary->getTypeKey()))
				{
					for (auto ent : boundary->getEntities())
					{
						auto gd_ent = geom_data->getEntityByName(ent->getName() + "_CH");
						mesh->addCut(mesh->getEntityByName(ent->getName()));
						mesh->getEntityByName(ent->getName() + "_CH")->setId(gd_ent->getId());  // set correct id to mesh group for cut
					}
				}
			}
		}
	}

	// ******************************
	//        Write out mesh
	// ******************************
	//TO DO fix inconsistency in mesh scale (mesh_scale and 1/mesh_scale)
	mesh->setMeshScale(1/mesh_scale);
	mesh->writeMesh(wDir + "/calculation/calcMesh.msh", std::string("GMSH"));
	return true;
}

bool GetdpSetup::update()
{
	if (caseMotions.empty())
		return false;

	initial_step = false;
	bool old_remesh = false;
	bool new_remesh = false;
	for (auto mt : caseMotions)
	{
		// skip if disabled
		if (!mt.isEnabled())
			continue;

		//skip if simple motion
		if (mt.getType() == Motion::SIMPLE)
			continue;

		//if motion type is old motion, move mesh immediately.
		// UI is not supposed to allow mixed type of motion (old+ new)
		// this is checked in UI
		if (mt.getType() == Motion::COMPLEX_WITH_PARAMETERS)
		{
			old_remesh = true;
			break;
		}

		//only new motion reaches this point
		new_remesh = true;
	}


	if (!new_remesh and !old_remesh)
		return false;

	if (current_time >= end_time)
		return false;

	for (auto& phys : casePhysics)
	{ 
		if (phys.getTimeType() == Physics::TRANSIENT)
		{
			phys.setTimeParameter("tstart", current_time);
			phys.setTimeParameter("tend", current_time + current_timestep);
			if (current_time > 0)
				phys.changeInitializationForRestart();
			std::cout << "Current time : " << current_time << std::endl;
			std::cout << "Timestep :" << current_timestep << std::endl;
		}
	}



	writeMain(PostType::SOLUTION);
	writeResolution();
	WriteConstraint();
	return true;
}

void GetdpSetup::finalize()
{
	std::vector<double> tvec = getTimeList(PostType::ENDTIME);
	current_time = tvec.back();


	mesh->writeMesh(wDir + "/geometry/lastMesh.msh", std::string("GMSH"), Mesh_IO::WRITE_NON_SCALED);
}

void GetdpSetup::stop()
{
	if (getTimeList(PostType::ENDTIME).size() != 0)
		copyResFiles();
	else
		no_results = true;
	stopped_during_calculation = true;
}

void GetdpSetup::validateInput()
{
	// *************************************
	// 	   This part to be removed when power control with accurate is working
	  //check if power control is used

	std::map<std::string, double> entity_power;
	bool has_current = false;
	bool has_voltage = false;
	for (auto phys : casePhysics)
	{
		for (auto dom : phys.getDomainsList())
		{
			for (auto prop : dom->getProperties())
			{
				if (prop.first == "target_power")
				{
					std::string nm = dom->getName();
					if (dom->getGroupKey() != "")
						nm = dom->getGroupKey();
					entity_power[nm] = prop.second.getDoubleValue();
				}
				if (prop.first == "_current" and prop.second.getDoubleValue() == 1.0)
					has_current = true;
				if (prop.first == "_voltage" and prop.second.getDoubleValue() == 1.0)
					has_voltage = true;
			}
		}

		for (auto bnd : phys.getBoundariesList())
		{
			for (auto prop : bnd->getProperties())
			{
				if (prop.key == "target_power")
				{
					entity_power[bnd->getName()] = prop.value.getDoubleValue();
				}
			}

			if (bnd->getTypeKey() == "current")
				has_current = true;
			if (bnd->getTypeKey() == "voltage")
				has_voltage = true;
		}
	}

	 // check if power control is selected with accurate
	if (caseNumerics.getAlgorithm() == Numerics::Accurate and entity_power.size() > 0)
	{
		throw(cenos_exception("Power control currently available only with Fast algorithm!"));
	}


	if (entity_power.size() > 1 and casePhysics[0].getSymmetryType() != Physics::SECTOR_3D)
	{
		std::string msg = "Only one power controlled source is allowed in simulation! \n";
		throw(cenos_exception(msg));
	}

	if (entity_power.size() > 0)
	{
		double first_power = entity_power.begin()->second;
		for (auto ep : entity_power)
		{
			if (ep.second != first_power)
			{
				std::string msg = "Only one power controlled source is allowed in simulation! \n";
				msg.append(" Power in " + entity_power.begin()->first + " is " + std::to_string(first_power) + ", and Power in " + ep.first + " is " + std::to_string(ep.second));
				throw(cenos_exception(msg));
			}
		}
	}


	if (entity_power.size() > 0 and has_current)
	{
		std::string msg = "Mixed sources (power + current) are not allowed in one model! \n";
		throw(cenos_exception(msg));
	}

	if (entity_power.size() > 0 and has_voltage)
	{
		std::string msg = "Mixed sources (power + voltage) are not allowed in one model! \n";
		throw(cenos_exception(msg));
	}

	// *************************************

}

void GetdpSetup::readGlobalsJson()
{
	std::string globalResultsFile = wDir + "\\results\\globalResultsFile.json";
	std::ifstream input(globalResultsFile);
	if (!input.is_open())
	{
		LOG(ERROR) << "globalResultsFile.json not found";
		return;
	}
	std::stringstream filebuff;
	filebuff << input.rdbuf();
	std::string inputJson = filebuff.str();
	std::istringstream iss(inputJson);
	iss >> globalsJson;
}

bool GetdpSetup::copyResFiles()
{
	std::string line;
	std::string forCopy;
	copy_file(wDir + "/calculation/main.pre", wDir + "/calculation/mainRunTime.pre", boost::filesystem::copy_option::overwrite_if_exists); // copy pre file
	copy_file(wDir + "/calculation/main.res", wDir + "/calculation/temp.res", boost::filesystem::copy_option::overwrite_if_exists); // copy pre file

	std::ifstream tempResFile;
	std::ofstream runTimeResFile;
	tempResFile.open(wDir + "/calculation/temp.res");
	runTimeResFile.open(wDir + "/mainRunTime.res");
	if (tempResFile.is_open())
	{
		// copies only those solution blocks which are closed with $EndSolution
		while (getline(tempResFile, line))
		{
			forCopy.append(line + "\n");
			if (line.find("$EndSolution") != std::string::npos)
			{
				runTimeResFile << forCopy;
				forCopy = "";
			}
		}
	}
	else
	{
		return false;
	}
	tempResFile.close();
	runTimeResFile.close();
	return true;
}


messagedata GetdpSetup::calculate()
{
	messagedata returnData;
	std::string getdpArgs = "\"" + wDir + "/calculation/main.pro\" -msh \"" + wDir + "/calculation/calcMesh.msh \" -solve analysis";

	messageQueue->push(messagedata(messagedata::info_lvl, messagedata::log, "Calculating."));
	batchRunner runner(wDir, "getdp", getdpArgs);

	try {
		auto runFuture = std::async(std::launch::async, &batchRunner::run, std::ref(runner), outputAnalyzer);

		do
		{
			Sleep(100);
			// Make sure stopper is not nullptr
			if ( *stopper)
			{
				if (!this->hasMultipleCalculations())
				{
					runner.stop();
				}
				runFuture.wait();
				this->stop();
			}

		} while (runFuture.wait_for(std::chrono::milliseconds(1)) != std::future_status::ready);

		returnData = runFuture.get();
	}
	catch (std::exception& ex)
	{
		return messagedata(messagedata::error, messagedata::response, "Error during calculation : " + std::string(ex.what()));
	}

	return returnData;
}

messagedata GetdpSetup::convertCaseNewResultJson() {
	if (this->hasNoResults())
		return messagedata(messagedata::error, messagedata::response, "Case contains no results and does not need to be converted.");

	try
	{
		this->outputResultsJson(true);
	}
	catch (std::exception& ex)
	{
		LOG(ERROR) << ex.what();
		return messagedata(messagedata::error, messagedata::response, "Exception caught: " + std::string(ex.what()));
	}
	return messagedata(messagedata::success, messagedata::response, "null");
	

}

messagedata GetdpSetup::postProcess()
{
	if (this->hasNoResults())
		return messagedata(messagedata::error, messagedata::response, "Calculation stopped by user. No results created.");

	messageQueue->push(messagedata(messagedata::info_lvl, messagedata::log, "Post processing results."));
	this->writePostSetup();
	messagedata postReturnData;

	// check if stop is requested before post processing
	bool stop_requested_before = *stopper;
	std::string lastPort = "";
	for (auto& phys : this->getCasePhysics())
	{
		if (phys.getSolverKey() == "getdpMicrowaves") {
			lastPort = "_" + phys.getLastPort();
		}
		
		std::string getdp_post_args;
		if (this->stoppedDuringCalculation())
			getdp_post_args = "\"" + wDir + "/calculation/main.pro\" -res " + wDir + "/calculation/temp.res -msh \"" + wDir + "/calculation/calcMesh.msh \" -pos postOp_" + phys.getSolverKey() + lastPort;
		else
			getdp_post_args = "\"" + wDir + "/calculation/main.pro\" -msh \"" + wDir + "/calculation/calcMesh.msh \" -pos postOp_" + phys.getSolverKey() + lastPort;


		try
		{
			batchRunner post_runner(wDir, "getdp", getdp_post_args);
			auto runFuture = std::async(std::launch::async, &batchRunner::run, std::ref(post_runner), outputAnalyzer);

			do
			{
				Sleep(100);

				// if stop is requested now (not during calculation) print message
				if (*stopper and !stop_requested_before)
				{
					std::cout << "Cannot stop during post processing " << std::endl;
					runFuture.wait();
				}

			} while (runFuture.wait_for(std::chrono::milliseconds(1)) != std::future_status::ready);

			postReturnData = runFuture.get();
			if (postReturnData.msgStatus != messagedata::success)
				return postReturnData;

		}
		catch (std::exception& ex)
		{
			LOG(ERROR) << ex.what();
			return messagedata(messagedata::error, messagedata::response, "Exception caught: " + std::string(ex.what()));
		}
	}

	this->convertResults();
	messageQueue->push(messagedata(messagedata::info_lvl, messagedata::log, "Calculating global variables."));
	this->outputGlobals();
	this->outputResultsJson(false);

	return postReturnData;
}