MySources/reducedmodelevaluator.cpp

#include "reducedmodelevaluator.hpp"
#include "hexagonremesher.hpp"
#include "trianglepatterngeometry.hpp"
#include <filesystem>

using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;

ReducedModelEvaluator::ReducedModelEvaluator()
{

}

std::vector<double> ReducedModelEvaluator::evaluateReducedModel(
		ReducedModelOptimization::Results &optimizationResults)
{
	//    std::shared_ptr<VCGPolyMesh> pTileIntoSurface = std::make_shared<VCGPolyMesh>();
	//    std::filesystem::path tileIntoSurfaceFilePath{
	//        "/home/iason/Coding/Projects/Approximating shapes with flat "
	//        "patterns/ReducedModelOptimization/TestSet/TileIntoSurface/plane_34Polygons.ply"};
	//    assert(std::filesystem::exists(tileIntoSurfaceFilePath));
	//    pTileIntoSurface->load(tileIntoSurfaceFilePath);

	//Set required file paths
	const std::filesystem::path tileInto_triMesh_filename
		= "/home/iason/Coding/build/PatternTillingReducedModel/Meshes/"
		"instantMeshes_plane_34.ply";
	const std::filesystem::path reducedPatternFilePath
		= "/home/iason/Coding/Projects/Approximating shapes with flat "
		"patterns/ReducedModelOptimization/TestSet/ReducedPatterns/single_reduced.ply";
	const std::filesystem::path intermediateResultsDirectoryPath
		= "/home/iason/Coding/build/ReducedModelOptimization/IntermediateResults";
	//    const std::filesystem::path drmSettingsFilePath
	//        = "/home/iason/Coding/Projects/Approximating shapes with flat "
	//          "patterns/ReducedModelOptimization/DefaultSettings/DRMSettings/"
	//          "defaultDRMSimulationSettings.json";
	DRMSimulationModel::Settings drmSimulationSettings;
	//    drmSimulationSettings.load(drmSettingsFilePath);
	drmSimulationSettings.linearGuessForceScaleFactor = 0.7;
	drmSimulationSettings.debugModeStep = 10000;
	drmSimulationSettings.beVerbose = true;
	constexpr bool shouldRerunFullPatternSimulation = false;
	const std::vector<std::string> scenariosTestSetLabels{"22Hex_randomBending0",
		"22Hex_randomBending1",
		"22Hex_randomBending2",
		"22Hex_randomBending3",
		"22Hex_randomBending4",
		"22Hex_randomBending5",
		"22Hex_randomBending6",
		"22Hex_randomBending7",
		"22Hex_randomBending8",
		"22Hex_randomBending9",
		"22Hex_randomBending10",
		"22Hex_randomBending11",
		"22Hex_randomBending12",
		"22Hex_randomBending13",
		"22Hex_randomBending14",
		"22Hex_randomBending15",
		"22Hex_randomBending16",
		"22Hex_randomBending17",
		"22Hex_randomBending18",
		"22Hex_randomBending19",
		"22Hex_randomBending20",
		"22Hex_bending_0.005N",
		"22Hex_bending_0.01N",
		"22Hex_bending_0.03N",
		"22Hex_bending_0.05N",
		"22Hex_pullOppositeVerts_0.05N",
		"22Hex_pullOppositeVerts_0.1N",
		"22Hex_pullOppositeVerts_0.3N",
		"22Hex_shear_2N",
		"22Hex_shear_5N",
		"22Hex_axial_10N",
		"22Hex_axial_20N",
		"22Hex_cylinder_0.05N",
		"22Hex_cylinder_0.1N",
		"22Hex_s_0.05N",
		"22Hex_s_0.1N"};

	//Load surface
	std::shared_ptr<VCGPolyMesh> pTileIntoSurface = [&]() {
		VCGTriMesh tileInto_triMesh;
		const bool surfaceLoadSuccessfull = tileInto_triMesh.load(tileInto_triMesh_filename);
		assert(surfaceLoadSuccessfull);
		return PolygonalRemeshing::remeshWithPolygons(tileInto_triMesh);
	}();
	const double optimizedBaseTriangleHeight = vcg::Distance(optimizationResults.baseTriangle.cP(0),
			(optimizationResults.baseTriangle.cP(1)
			 + optimizationResults.baseTriangle.cP(
				 2))
			/ 2);
	pTileIntoSurface->moveToCenter();
	const double scaleFactor = optimizedBaseTriangleHeight
		/ pTileIntoSurface->getAverageFaceRadius();
	vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(*pTileIntoSurface, scaleFactor);
	//    tileIntoSurface.registerForDrawing();
	//    polyscope::show();

	//Tile full pattern into surface
	std::vector<PatternGeometry> fullPatterns(1);
	fullPatterns[0].copy(optimizationResults.baseTriangleFullPattern);
	////        Base triangle pattern might contain dangling vertices.Remove those
	fullPatterns[0].interfaceNodeIndex = 3;
	fullPatterns[0].deleteDanglingVertices();
	std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN(), 0);
	std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
	std::vector<size_t> tileIntoEdgeToTiledFullVi;
	std::shared_ptr<PatternGeometry> pTiledFullPattern
		= PatternGeometry::tilePattern(fullPatterns,
				{},
				*pTileIntoSurface,
				perSurfaceFacePatternIndices,
				tileIntoEdgeToTiledFullVi,
				perPatternIndexTiledFullPatternEdgeIndices);
	pTiledFullPattern->setLabel("Tiled_full_patterns");
	//    pTiledFullPattern->registerForDrawing();
	//Tile reduced pattern into surface
	PatternGeometry reducedPattern;
	reducedPattern.load(reducedPatternFilePath);
	reducedPattern.deleteDanglingVertices();
	assert(reducedPattern.interfaceNodeIndex == 1);
	std::vector<PatternGeometry> reducedPatterns(1);
	reducedPatterns[0].copy(reducedPattern);
	const auto reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
	ReducedModelOptimization::Results::applyOptimizationResults_innerHexagon(
			optimizationResults, reducedPatternBaseTriangle, reducedPatterns[0]);

	std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
	std::vector<size_t> tileIntoEdgeToTiledReducedVi;
	std::shared_ptr<PatternGeometry> pTiledReducedPattern
		= PatternGeometry::tilePattern(reducedPatterns,
				{0},
				*pTileIntoSurface,
				perSurfaceFacePatternIndices,
				tileIntoEdgeToTiledReducedVi,
				perPatternIndexTiledReducedPatternEdgeIndices);
	pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
	//        pTiledReducedPattern->registerForDrawing();

	std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
		fullToReducedViMap; //of only the common vertices
	std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
		reducedToFullViMap; //of only the common vertices
	for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
		fullToReducedViMap[tileIntoEdgeToTiledFullVi[ei]] = tileIntoEdgeToTiledReducedVi[ei];
	}
	constructInverseMap(fullToReducedViMap, reducedToFullViMap);

	std::vector<size_t> tilledFullPatternInterfaceVi;
	tilledFullPatternInterfaceVi.clear();
	tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
	size_t viIndex = 0;
	for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
		tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
	}

	//Create simulation meshes
	////Tessellated full pattern simulation mesh
	std::shared_ptr<SimulationMesh> pTiledFullPattern_simulationMesh;
	pTiledFullPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledFullPattern);
	//NOTE: Those should be derived from the optimization results instead of hardcoding them
	pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
	if (optimizationResults.fullPatternYoungsModulus == 0) {
		std::cerr << "Full pattern's young modulus not found." << std::endl;
		std::terminate();
	}
	pTiledFullPattern_simulationMesh->setBeamMaterial(0.3,
			optimizationResults.fullPatternYoungsModulus);
	pTiledFullPattern_simulationMesh->reset();

	////Tessellated reduced pattern simulation mesh
	std::shared_ptr<SimulationMesh> pTiledReducedPattern_simulationMesh;
	pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
	const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
		= perPatternIndexTiledReducedPatternEdgeIndices[0];
	ReducedModelOptimization::Results::applyOptimizationResults_elements(
			optimizationResults, pTiledReducedPattern_simulationMesh);
	pTiledReducedPattern_simulationMesh->reset();
	std::vector<double> distances_drm2reduced(scenariosTestSetLabels.size(), 0);
	std::vector<double> distances_normalizedDrm2reduced(scenariosTestSetLabels.size(), 0);
	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
		const std::filesystem::path scenariosDirectoryPath
			= "/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/";
		const std::filesystem::path tiledReducedPatternJobFilePath
			= std::filesystem::path(scenariosDirectoryPath)
			.append(jobLabel)
			.append("SimulationJobs")
			.append("Reduced")
			.append(SimulationJob::jsonDefaultFileName);
		//set jobs
		std::shared_ptr<SimulationJob> pJob_tiledReducedPattern;
		pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
		pJob_tiledReducedPattern->load(tiledReducedPatternJobFilePath, false);
		pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
		std::shared_ptr<SimulationJob> pJob_tiledFullPattern;
		pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
		pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
		pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
		//        pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
		//        pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
		//    polyscope::show();
		//Save reduced job
		const std::filesystem::path tesellatedResultsFolderPath
			= std::filesystem::path(intermediateResultsDirectoryPath).append("TessellatedResults");
		const std::filesystem::path surfaceFolderPath = std::filesystem::path(
				tesellatedResultsFolderPath)
			.append(pTileIntoSurface->getLabel());
		const std::string scenarioLabel = pJob_tiledFullPattern->getLabel();
		const std::filesystem::path scenarioDirectoryPath = std::filesystem::path(surfaceFolderPath)
			.append(scenarioLabel);
		const std::filesystem::path reducedJobDirectoryPath
			= std::filesystem::path(scenarioDirectoryPath).append("ReducedJob");
		std::filesystem::create_directories(reducedJobDirectoryPath);
		pJob_tiledReducedPattern->save(reducedJobDirectoryPath);
		//Run scenario

		////Full

		const std::string patternLabel = [&]() {
			const std::string patternLabel = optimizationResults.baseTriangleFullPattern.getLabel();
			const int numberOfOccurences = std::count_if(patternLabel.begin(),
					patternLabel.end(),
					[](char c) { return c == '#'; });
			if (numberOfOccurences == 0) {
				return std::to_string(optimizationResults.baseTriangleFullPattern.EN()) + "#"
					+ optimizationResults.baseTriangleFullPattern.getLabel();
			} else if (numberOfOccurences == 1) {
				return optimizationResults.baseTriangleFullPattern.getLabel();
			}
		}();
		const auto fullResultsFolderPath
			= std::filesystem::path(scenarioDirectoryPath).append(patternLabel).append("Results");
		if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
			std::filesystem::remove_all(fullResultsFolderPath);
		}

		const std::filesystem::path fullPatternJobFolderPath = std::filesystem::path(
				scenarioDirectoryPath)
			.append(patternLabel)
			.append("SimulationJob");
		SimulationResults simulationResults_tiledFullPattern_drm;
		if (std::filesystem::exists(fullResultsFolderPath)) {
			//Load full pattern results
			assert(std::filesystem::exists(fullPatternJobFolderPath));
			simulationResults_tiledFullPattern_drm.load(fullResultsFolderPath,
					fullPatternJobFolderPath);
			simulationResults_tiledFullPattern_drm.converged = true;
		} else {
			//Full
			std::cout << "Executing:" << jobLabel << std::endl;
			DRMSimulationModel drmSimulationModel;
			simulationResults_tiledFullPattern_drm
				= drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
			simulationResults_tiledFullPattern_drm.setLabelPrefix("DRM");
		}
		std::filesystem::create_directories(fullResultsFolderPath);
		simulationResults_tiledFullPattern_drm.save(
				std::filesystem::path(scenarioDirectoryPath).append(patternLabel));
		if (!simulationResults_tiledFullPattern_drm.converged) {
			std::cerr << "Full pattern simulation failed." << std::endl;
		}

		LinearSimulationModel linearSimulationModel;
		SimulationResults simulationResults_tiledReducedPattern
			= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);

		//        simulationResults_tiledReducedPattern.registerForDrawing();
		//        simulationResults_tiledFullPattern_drm.registerForDrawing();
		//        polyscope::show();
		//measure distance
		const double distance_fullDrmToReduced
			= simulationResults_tiledFullPattern_drm
			.computeDistance(simulationResults_tiledReducedPattern, fullToReducedViMap);
		double distance_fullSumOfAllVerts = 0;
		for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
			distance_fullSumOfAllVerts += simulationResults_tiledFullPattern_drm
				.displacements[fullToReducedPair.first]
				.getTranslation()
				.norm();
		}
		const double distance_normalizedFullDrmToReduced = distance_fullDrmToReduced
			/ distance_fullSumOfAllVerts;
		distances_drm2reduced[jobIndex] = distance_fullDrmToReduced;
		distances_normalizedDrm2reduced[jobIndex] = distance_normalizedFullDrmToReduced;
	}
	//#ifndef POLYSCOPE_DEFINED
	//    return distances_drm2reduced;
	//#else

	//report distance
	//        csvFile csv_results("", flse);
	csvFile csv_results({}, false);
	//    csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
	//                            .append("results.csv")
	//                            .string(),
	//                        false);
	csv_results /*<< "Job Label"*/
		<< "drm2Reduced"
		<< "norm_drm2Reduced";
	csv_results << endrow;
	double sumOfNormalizedFull2Reduced = 0;
	for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
		const std::string &jobLabel = scenariosTestSetLabels[jobIndex];
		const double &distance_fullDrmToReduced = distances_drm2reduced[jobIndex];
		const double &distance_normalizedFullDrmToReduced = distances_normalizedDrm2reduced[jobIndex];
		csv_results /*<< jobLabel*/ << distance_fullDrmToReduced
			<< distance_normalizedFullDrmToReduced;
		csv_results << endrow;
		sumOfNormalizedFull2Reduced += distance_normalizedFullDrmToReduced;
	}

	std::cout << "Average normalized error per scenario:"
		<< sumOfNormalizedFull2Reduced / scenariosTestSetLabels.size() << std::endl;
	return distances_normalizedDrm2reduced;
	//#endif
}