PatternTillingReducedModel/src/reducedpatternsimulator.cpp

1616 lines
79 KiB
C++

#include "reducedpatternsimulator.hpp"
#include "imgui.h"
#include "misc/cpp/imgui_stdlib.h"
#include "simulationhistoryplotter.hpp"
#include <wrap/io_trimesh/export_obj.h>
#include <algorithm>
void ReducedPatternSimulator::applyOptimizationResults_elements(
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh,
const std::vector<size_t> &elementIndices)
{
assert(CrossSectionType::name == RectangularBeamDimensions::name);
// Set reduced parameters fron the optimization results
std::unordered_map<std::string, double>
optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
reducedPattern_optimizationResults.optimalXNameValuePairs.end());
const std::string ALabel = "A";
assert(optimalXVariables.contains(ALabel));
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
CrossSectionType elementDimensions(beamWidth, beamHeight);
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";
assert(optimalXVariables.contains(ymLabel));
const double E = optimalXVariables.at(ymLabel);
const ElementMaterial elementMaterial(poissonsRatio, E);
const std::string JLabel = "J";
assert(optimalXVariables.contains(JLabel));
const double J = optimalXVariables.at(JLabel);
const std::string I2Label = "I2";
assert(optimalXVariables.contains(I2Label));
const double I2 = optimalXVariables.at(I2Label);
const std::string I3Label = "I3";
assert(optimalXVariables.contains(I3Label));
const double I3 = optimalXVariables.at(I3Label);
for (size_t ei : elementIndices) {
Element &e = pTiledReducedPattern_simulationMesh->elements[ei];
e.setDimensions(elementDimensions);
e.setMaterial(elementMaterial);
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
}
std::pair<DRMFullSimulationResults, ReducedSimulationResults> ReducedPatternSimulator::runSimulations(
const std::string &scenarioLabel,
const bool &shouldRerunFullPatternSimulation,
const bool &showReducedPatternResultsBeforeRunningDRM)
{
if (pJob_tiledFullPattern->isEmpty() || pJob_tiledFullPattern->isEmpty()) {
std::cout << "Empty simulation job. Load or create a simulation job to run a simulation."
<< std::endl;
return std::make_pair(SimulationResults(), SimulationResults());
}
// std::cout << "Running " << scenarioLabel << std::endl;
pJob_tiledFullPattern->label = scenarioLabel;
pJob_tiledReducedPattern->label = scenarioLabel;
const std::filesystem::path surfaceOutputPath
= std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/Scenarios/")
.append(surfaceLabel);
//Tesselation
const std::filesystem::path tesselationFolderPath = std::filesystem::path(surfaceOutputPath)
.append("Tessellations");
saveTesselation(tesselationFolderPath);
const std::filesystem::path scenarioOutputPath(
std::filesystem::path(surfaceOutputPath).append(scenarioLabel));
const std::string simJobsOutputFolderPath(
std::filesystem::path(scenarioOutputPath).append("SimulationJobs").string());
std::filesystem::create_directories(simJobsOutputFolderPath);
// if (!std::filesystem::exists(simJobsOutputFolderPath)) {
saveJobs(simJobsOutputFolderPath);
// }
//Reduced
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_reducedPattern = linearSimulationModel.executeSimulation(
pJob_tiledReducedPattern);
simulationResults_reducedPattern.setLabelPrefix("Reduced");
const bool wasSuccessful = simulationResults_reducedPattern.saveDeformedModel(
scenarioOutputPath.string());
SimulationResults simulationResults_fullPattern_linear = linearSimulationModel.executeSimulation(
pJob_tiledFullPattern);
simulationResults_fullPattern_linear.setLabelPrefix("FullLinear");
simulationResults_fullPattern_linear.saveDeformedModel(scenarioOutputPath.string());
if (showReducedPatternResultsBeforeRunningDRM) {
simulationResults_reducedPattern.registerForDrawing();
simulationResults_fullPattern_linear.registerForDrawing();
polyscope::show();
// polyscope::requestRedraw();
simulationResults_reducedPattern.unregister();
simulationResults_fullPattern_linear.unregister();
}
// const auto reducedResultsFolderPath = std::filesystem::path(outputPath)
// .append("ReducedResults")
// .append(fullPatternsLabel);
// std::filesystem::create_directories(reducedResultsFolderPath);
// simulationResults_reducedPattern.save(reducedResultsFolderPath);
//Full
const auto fullResultsFolderPath
= std::filesystem::path(scenarioOutputPath).append("FullResults").append(fullPatternsLabel);
if (shouldRerunFullPatternSimulation && std::filesystem::exists(fullResultsFolderPath)) {
std::filesystem::remove_all(fullResultsFolderPath);
}
const std::filesystem::path fullPatternJobFolderPath
= std::filesystem::path(simJobsOutputFolderPath).append(fullPatternsLabel);
SimulationResults simulationResults_fullPattern_drm;
if (std::filesystem::exists(fullResultsFolderPath) && !randomTesselationIsEnabled) {
//Load full pattern results
assert(std::filesystem::exists(fullPatternJobFolderPath));
simulationResults_fullPattern_drm.load(fullResultsFolderPath, fullPatternJobFolderPath);
simulationResults_fullPattern_drm.converged = true;
} else {
//Full
DRMSimulationModel drmSimulationModel;
DRMSimulationModel::Settings drmSimulationSettings;
drmSimulationSettings.isDebugMode = true;
drmSimulationSettings.debugModeStep = 50000;
// drmSimulationSettings.shouldCreatePlots = true;
// drmSimulationSettings.drawingStep = 50000;
// drmSimulationSettings.shouldDraw = true;
// drmSimulationSettings.drawingStep = 100 * drmSimulationSettings.debugModeStep;
drmSimulationSettings.beVerbose = true;
// drmSimulationSettings.desiredGradualExternalLoadsSteps = 10;
// drmSimulationSettings.shouldCreatePlots = true;
drmSimulationSettings.useAverage = true;
// drmSimulationSettings.displacementCap = 0.1 * surfaceBaseTriangleHeight;
// drmSimulationSettings.xi = 0.97;
// drmSimulationSettings.Dtini = 0.06;
// drmSimulationSettings.gamma = 0.3;
drmSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-1;
// drmSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// drmSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-15;
// PolyscopeInterface::deinitPolyscope();
const double forceScaleFactor = 0.5;
for (auto &externalForce : pJob_tiledFullPattern->nodalExternalForces) {
externalForce.second = externalForce.second * forceScaleFactor;
}
SimulationResults simulationResults_fullPatternLinearModel
= linearSimulationModel.executeSimulation(pJob_tiledFullPattern);
// simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob_tiledFullPattern->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
simulationResults_fullPattern_drm
= drmSimulationModel.executeSimulation(pJob_tiledFullPattern,
drmSimulationSettings,
simulationResults_fullPatternLinearModel);
// simulationResults_fullPattern_drm
// = drmSimulationModel.executeSimulation(pJob_tiledFullPattern, drmSimulationSettings);
if (!simulationResults_fullPattern_drm.converged) {
std::cerr << "Full pattern simulation failed." << std::endl;
} else if (!randomTesselationIsEnabled) {
std::filesystem::create_directories(fullResultsFolderPath);
simulationResults_fullPattern_drm.save(fullResultsFolderPath.string());
pJob_tiledFullPattern->save(fullPatternJobFolderPath.string());
}
}
return std::make_pair(simulationResults_fullPattern_drm, simulationResults_reducedPattern);
}
void ReducedPatternSimulator::tileReducedPatterns(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<int> &perSurfaceFacePatternIndices,
std::shared_ptr<SimulationMesh> &pSimulationMesh_tiledReduced,
std::vector<size_t> &tileIntoEdgeToTiledReducedVi)
{
std::vector<PatternGeometry> reducedPatterns(optimizationResults.size());
for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size();
reducedPatternIndex++) {
reducedPatterns[reducedPatternIndex].copy(reducedPattern);
ReducedPatternOptimization::Results::applyOptimizationResults_innerHexagon(
optimizationResults[reducedPatternIndex],
reducedPatternBaseTriangle,
reducedPatterns[reducedPatternIndex]);
}
std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
std::shared_ptr<PatternGeometry> pTiledReducedPattern
= PatternGeometry::tilePattern(reducedPatterns,
{0},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledReducedVi,
perPatternIndexTiledReducedPatternEdgeIndices);
pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
std::size_t shuffleHash = computeHashOrdered(perSurfaceFacePatternIndices);
if (!shuffleToNumOfOccur.contains(shuffleHash)) {
shuffleToNumOfOccur[shuffleHash] = 0;
} else {
shuffleToNumOfOccur[shuffleHash]++;
}
if (shuffleToNumOfOccur[shuffleHash] != 0) {
std::cout << shuffleToNumOfOccur[shuffleHash] << std::endl;
}
// polyscope::show();
// pTiledReducedPattern->unregister();
pSimulationMesh_tiledReduced.reset();
pSimulationMesh_tiledReduced = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size();
reducedPatternIndex++) {
const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
= perPatternIndexTiledReducedPatternEdgeIndices[reducedPatternIndex];
applyOptimizationResults_elements(optimizationResults[reducedPatternIndex],
pTiledReducedPattern_simulationMesh,
tiledPatternElementIndicesForReducedPattern);
}
pTiledReducedPattern_simulationMesh->reset();
}
void ReducedPatternSimulator::shuffleReducedPatterns(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<size_t> &tileIntoEdgeToTiledFullPattern,
const std::vector<int> &perSurfaceFacePatternIndices)
{
std::vector<int> randomShufflePerSurfaceFacePatternIndices = perSurfaceFacePatternIndices;
std::random_device rd;
std::mt19937 g(rd());
std::shuffle(randomShufflePerSurfaceFacePatternIndices.begin(),
randomShufflePerSurfaceFacePatternIndices.end(),g);
std::vector<size_t> tileIntoEdgesToTiledReducedVi;
tileReducedPatterns(optimizationResults,
randomShufflePerSurfaceFacePatternIndices,
pTiledReducedPattern_simulationMesh,
tileIntoEdgesToTiledReducedVi);
constructViMaps(tileIntoEdgeToTiledFullPattern, tileIntoEdgesToTiledReducedVi);
}
double ReducedPatternSimulator::computeDistance(const std::string &scenarioLabel,
const bool &shouldRerunDRMSimulation,
const bool &shouldDraw)
{
loadScenario(scenarioLabel);
const auto results = runSimulations(scenarioLabel, shouldRerunDRMSimulation);
if (shouldDraw) {
results.first.registerForDrawing();
results.second.registerForDrawing();
polyscope::show();
results.first.unregister();
results.second.unregister();
}
const double distance = computeDistance(results.first, results.second, fullToReducedViMap)
;
return distance;
}
std::vector<double> ReducedPatternSimulator::computeDistancesPerScenario(
const std::vector<std::string> &scenarioLabels,
const bool &shouldRerunDRMSimulation,
const bool &shouldDraw)
{
std::vector<double> fullReducedDistancesPerScenario(scenarioLabels.size());
for (size_t scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
const std::string &scenarioLabel = scenarioLabels[scenarioIndex];
if (shouldRerunDRMSimulation) {
std::cout << "Running:" << scenarioLabel << std::endl;
}
fullReducedDistancesPerScenario[scenarioIndex] = computeDistance(scenarioLabel,
shouldRerunDRMSimulation,
shouldDraw);
// std::cout << "Distance:" << fullReducedDistancesPerScenario[scenarioIndex] << std::endl;
reset();
}
return fullReducedDistancesPerScenario;
}
std::vector<int> ReducedPatternSimulator::computePerSurfaceFacePatternsIndices(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults) const
{
std::vector<int> perSurfaceFacePatternIndices(pTileIntoSurface->FN());
if (randomTesselationIsEnabled) {
for (size_t tileIntoFi = 0; tileIntoFi < pTileIntoSurface->FN(); tileIntoFi++) {
const size_t randomPatternIndex = (rand()
% static_cast<size_t>(optimizationResults.size()));
perSurfaceFacePatternIndices[tileIntoFi] = randomPatternIndex;
}
} else {
const size_t faceIndexStep = pTileIntoSurface->FN() / optimizationResults.size();
for (OptimizationResultsIndex optResIndex = 0; optResIndex < optimizationResults.size();
optResIndex++) {
std::fill_n(perSurfaceFacePatternIndices.begin() + optResIndex * faceIndexStep,
faceIndexStep,
optResIndex);
}
}
return perSurfaceFacePatternIndices;
}
std::vector<ReducedPatternOptimization::Results> ReducedPatternSimulator::loadOptimizationResults(
const std::filesystem::path &optimizationResultsFolderPath)
{
std::set<std::filesystem::path> optimizationResultPaths;
for (const auto &dirEntry : std::filesystem::directory_iterator(optimizationResultsFolderPath)) {
// std::cout << "Loading:" << dirEntry.path() << std::endl;
optimizationResultPaths.insert(dirEntry.path());
}
std::vector<ReducedPatternOptimization::Results> optimizationResults(
optimizationResultPaths.size());
size_t optimizationResultIndex = 0;
for (const auto &p : optimizationResultPaths) {
optimizationResults[optimizationResultIndex++].load(p.string());
}
return optimizationResults;
}
void ReducedPatternSimulator::tileFullPatterns(
std::vector<ConstPatternGeometry> &fullPatterns,
const std::vector<int> &perSurfaceFacePatternIndices,
std::shared_ptr<SimulationMesh> &pTiledFullPattern_simulationMesh,
std::vector<size_t> &tileIntoEdgeToTiledFullVi)
{
///Create tiled configurations
std::vector<std::vector<size_t>> perPatternIndexTiledFullPatternEdgeIndices;
std::shared_ptr<PatternGeometry> pTiledFullPattern
= PatternGeometry::tilePattern(fullPatterns,
{},
*pTileIntoSurface,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledFullVi,
perPatternIndexTiledFullPatternEdgeIndices);
pTiledFullPattern->setLabel("Tiled_full_patterns");
//#ifdef POLYSCOPE_DEFINED
// pTiledFullPattern->registerForDrawing();
// polyscope::show();
// pTiledFullPattern->unregister();
//#endif
pTiledFullPattern_simulationMesh.reset();
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});
pTiledFullPattern_simulationMesh->setBeamMaterial(0.3, 1 * 1e9);
}
void ReducedPatternSimulator::constructViMaps(
const std::vector<size_t> &tileIntoEdgeToTiledFullPattern,
const std::vector<size_t> &tileIntoEdgeToTiledReducedPattern)
{
fullToReducedViMap.clear();
for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
fullToReducedViMap[tileIntoEdgeToTiledFullPattern[ei]]
= tileIntoEdgeToTiledReducedPattern[ei];
}
constructInverseMap(fullToReducedViMap, reducedToFullViMap);
tilledFullPatternInterfaceVi.clear();
tilledFullPatternInterfaceVi.resize(fullToReducedViMap.size());
size_t viIndex = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
tilledFullPatternInterfaceVi[viIndex++] = fullToReducedPair.first;
}
}
void ReducedPatternSimulator::createShufflings(
std::vector<ReducedPatternOptimization::Results> &optimizationResults)
{
std::vector<int> shuffledPerSurfaceFacePatternIndices = computePerSurfaceFacePatternsIndices(
optimizationResults);
std::random_device rd;
std::mt19937 g(rd());
std::shuffle(shuffledPerSurfaceFacePatternIndices.begin(),
shuffledPerSurfaceFacePatternIndices.end(),g);
createTiledSimulationMeshes(pTileIntoSurface,
optimizationResults,
shuffledPerSurfaceFacePatternIndices);
}
void ReducedPatternSimulator::runWeightEvaluation(
const std::filesystem::path &optimizationResultsFolderPath)
{
assert(pTileIntoSurface->VN() == 62 && pTileIntoSurface->EN() == 83
&& pTileIntoSurface->FN() == 22);
std::filesystem::path bestPerformingFolder;
double minDist = std::numeric_limits<double>::max();
std::vector<std::vector<int>> shufflings = csvFile::parse<int>(
std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo")
.append("shufflings.csv"));
const size_t numberOfShufflings = shufflings.size();
for (auto &p : std::filesystem::directory_iterator(optimizationResultsFolderPath)) {
if (!std::filesystem::is_directory(std::filesystem::path(p))) {
continue;
}
std::vector<ReducedPatternOptimization::Results> optimizationResults
= loadOptimizationResults(std::filesystem::path(p).append("ConvergedJobs"));
double averageDistanceOverAllShufflings = 0;
// computePerSurfaceFacePatternsIndices(optimizationResults);
for (std::vector<int> shuffledPerSurfaceFacePatternIndices : shufflings) {
// for (const auto &el : shuffledPerSurfaceFacePatternIndices) {
// std::cout << el << " ";
// }
// std::cout << std::endl;
createTiledSimulationMeshes(pTileIntoSurface,
optimizationResults,
shuffledPerSurfaceFacePatternIndices);
std::vector<double> perScenarioDistances = computeDistancesPerScenario(
scenariosTestSetLabels);
averageDistanceOverAllShufflings += std::accumulate(perScenarioDistances.begin(),
perScenarioDistances.end(),
0.0)
/ numberOfShufflings;
}
if (averageDistanceOverAllShufflings < minDist) {
bestPerformingFolder = p;
minDist = averageDistanceOverAllShufflings;
}
std::cout << "Evaluated:" << p.path().string() << std::endl;
std::cout << "Dist:" << averageDistanceOverAllShufflings << std::endl;
}
std::cout << "The best performing set was at:" << bestPerformingFolder << std::endl;
}
void ReducedPatternSimulator::runShufflingEvaluation(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<size_t> &tileIntoEiToTiledFullVi)
{
assert(pTileIntoSurface->VN() == 62 && pTileIntoSurface->EN() == 83
&& pTileIntoSurface->FN() == 22);
std::vector<int> perSurfaceFacePatternIndices = computePerSurfaceFacePatternsIndices(
optimizationResults);
std::vector<double> nonShuffledDistances = computeDistancesPerScenario(scenariosTestSetLabels);
const size_t numberOfShufflings = 500;
std::vector<std::vector<double>> shuffledDistances(numberOfShufflings,
std::vector<double>(
scenariosTestSetLabels.size()));
for (size_t shufflingIndex = 0; shufflingIndex < numberOfShufflings; shufflingIndex++) {
shuffleReducedPatterns(optimizationResults,
tileIntoEiToTiledFullVi,
perSurfaceFacePatternIndices);
shuffledDistances[shufflingIndex] = computeDistancesPerScenario(scenariosTestSetLabels);
if (shufflingIndex % (numberOfShufflings / 10) == 0) {
std::cout << "Shuffling index:" << shufflingIndex << std::endl;
}
// for (size_t scenarioIndex = 0; scenarioIndex < scenariosTestSetLabels.size();
// scenarioIndex++) {
// const bool shuffledPerformsBetter = shuffledDistances[shufflingIndex][scenarioIndex]
// < nonShuffledDistances[scenarioIndex];
// if (shuffledPerformsBetter) {
// computeDistance(scenariosTestSetLabels[scenarioIndex], false, true);
// }
// }
}
//Compute results
size_t numberOfShuffledTesselationsPerformingBetter = 0;
size_t totalNumberOfScenariosPerformingBetter = 0;
std::vector<double> numberOfBetterPerformingShufflinsPerScenario(scenariosTestSetLabels.size(),
0);
for (size_t shufflingIndex = 0; shufflingIndex < numberOfShufflings; shufflingIndex++) {
size_t scenariosPerformingBetter = 0;
for (size_t scenarioIndex = 0; scenarioIndex < scenariosTestSetLabels.size();
scenarioIndex++) {
const bool shuffledPerformsBetter = shuffledDistances[shufflingIndex][scenarioIndex]
< nonShuffledDistances[scenarioIndex];
// if (shuffledPerformsBetter) {
// int i = 0;
// i++;
// }
scenariosPerformingBetter += shuffledPerformsBetter;
numberOfBetterPerformingShufflinsPerScenario[scenarioIndex] += shuffledPerformsBetter;
}
// const float performingBetterThreshold = 0.8;
// if (static_cast<double>(scenariosPerformingBetter) / scenariosTestSetLabels.size()
// >= performingBetterThreshold) {
// numberOfShuffledTesselationsPerformingBetter++;
// }
totalNumberOfScenariosPerformingBetter += scenariosPerformingBetter;
}
matplot::bar(numberOfBetterPerformingShufflinsPerScenario);
std::vector<std::string> scenarioTestSetSimplifiedLabels(scenariosTestSetLabels.size());
std::transform(
scenariosTestSetLabels.begin(),
scenariosTestSetLabels.end(),
scenarioTestSetSimplifiedLabels.begin(),
[](const std::string &label) {
std::string simplifiedLabel(label.begin() + 6, label.end());
simplifiedLabel.erase(std::remove(simplifiedLabel.begin(), simplifiedLabel.end(), '_'),
simplifiedLabel.end());
const size_t beginPos_pov = simplifiedLabel.find("pullOppositeVerts");
if (beginPos_pov != simplifiedLabel.npos) {
simplifiedLabel.replace(beginPos_pov,
std::string("pullOppositeVerts").size(),
"POV");
}
const size_t beginPos_random = simplifiedLabel.find("randomBending");
if (beginPos_random != simplifiedLabel.npos) {
simplifiedLabel.replace(beginPos_random,
std::string("randomBending").size(),
"random");
}
return simplifiedLabel;
});
matplot::gca()->x_axis().ticklabels(scenarioTestSetSimplifiedLabels);
matplot::gca()->y_axis().limits(std::array<double, 2>{0, numberOfShufflings});
matplot::show();
std::cout << "Number of shufflings:" << numberOfShufflings << std::endl;
std::cout << "Average scenarios performing better/shuffling:"
<< static_cast<double>(totalNumberOfScenariosPerformingBetter) / numberOfShufflings
<< std::endl;
// std::cout << "Number of shufflings performing better:"
// << numberOfShuffledTesselationsPerformingBetter << std::endl;
// std::cout << "Percentage:"
// << static_cast<double>(numberOfShuffledTesselationsPerformingBetter)
// / numberOfShufflings
// << std::endl;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>
ReducedPatternSimulator::getPickedInterfaceVertices(
const std::pair<std::string, size_t> &selection) const
{
auto tiledFullPattern_polyscopeHandle = polyscope::getCurveNetwork(
pTiledFullPattern_simulationMesh->getLabel());
auto tiledReducedPattern_polyscopeHandle = polyscope::getCurveNetwork(
pTiledReducedPattern_simulationMesh->getLabel());
const std::string &pickedStructureName = selection.first;
const bool pickedFull = pickedStructureName == tiledFullPattern_polyscopeHandle->name;
const bool pickedReduced = pickedStructureName == tiledReducedPattern_polyscopeHandle->name;
if (!(pickedFull || pickedReduced)) {
return std::make_pair(-1, -1);
}
const bool vertexWasPicked = pickedFull
? selection.second < tiledFullPattern_polyscopeHandle->nNodes()
: selection.second
< tiledReducedPattern_polyscopeHandle->nNodes();
if (!vertexWasPicked) {
return std::make_pair(-1, -1);
}
const int &vi = selection.second;
assert(!pickedStructureName.empty());
int fullVi = -1;
int reducedVi = -1;
if (pickedFull) {
fullVi = vi;
if (!fullToReducedViMap.contains(fullVi)) {
return std::make_pair(-1, -1);
}
reducedVi = fullToReducedViMap.at(vi);
} else if (pickedReduced) {
reducedVi = vi;
if (!reducedToFullViMap.contains(reducedVi)) {
return std::make_pair(-1, -1);
}
fullVi = reducedToFullViMap.at(vi);
}
return std::make_pair(fullVi, reducedVi);
}
void ReducedPatternSimulator::removeDrawnSimulationJobs()
{
//Keep only initial patterns
for (std::pair<std::string, std::map<std::string, polyscope::Structure *>>
polyscopeStructureCategory : polyscope::state::structures) {
for (std::pair<std::string, polyscope::Structure *> polyscopeStructure :
polyscopeStructureCategory.second) {
if ((pTiledReducedPattern_simulationMesh
&& polyscopeStructure.first != pTiledReducedPattern_simulationMesh->getLabel())
&& polyscopeStructure.first != pTiledFullPattern_simulationMesh->getLabel()
&& polyscopeStructure.first != pTileIntoSurface->getLabel()) {
polyscope::removeStructure(polyscopeStructure.second);
} else {
if (dynamic_cast<polyscope::SurfaceMesh *>(polyscopeStructure.second)) {
dynamic_cast<polyscope::SurfaceMesh *>(polyscopeStructure.second)
->removeAllQuantities();
} else if (dynamic_cast<polyscope::CurveNetwork *>(polyscopeStructure.second))
dynamic_cast<polyscope::CurveNetwork *>(polyscopeStructure.second)
->removeAllQuantities();
}
}
}
//Remove simulation job quantities from the initial meshes
if (pTiledReducedPattern_simulationMesh && pJob_tiledReducedPattern) {
pJob_tiledReducedPattern->unregister(pTiledReducedPattern_simulationMesh->getLabel());
pJob_tiledReducedPattern->clear();
}
if (pTiledFullPattern_simulationMesh && pJob_tiledFullPattern) {
pJob_tiledFullPattern->unregister(pTiledFullPattern_simulationMesh->getLabel());
pJob_tiledFullPattern->clear();
}
// if (pTiledReducedPattern_simulationMesh) {
// pTiledReducedPattern_simulationMesh->unregister();
// }
// if (pTiledFullPattern_simulationMesh) {
// pTiledFullPattern_simulationMesh->unregister();
// }
}
void ReducedPatternSimulator::resetUserSelectedVertices()
{
gui_fullPatternSelectedVertices.clear();
gui_fullVerticesColors.clear();
if (pTiledFullPattern_simulationMesh) {
gui_fullVerticesColors.resize(pTiledFullPattern_simulationMesh->VN());
}
gui_reducedVerticesColors.clear();
if (pTiledReducedPattern_simulationMesh) {
gui_reducedVerticesColors.resize(pTiledReducedPattern_simulationMesh->VN());
}
polyscope::pick::resetSelection();
}
void ReducedPatternSimulator::resetUserSelectedFaces()
{
gui_faceToPatternIndex.clear();
gui_faceToPatternIndex.resize(pTileIntoSurface->FN(), -1);
gui_currentColorPatternIndexPair = {glm::vec3(1, 0, 0), 0};
gui_colorsPerFace.clear();
gui_colorsPerFace.resize(pTileIntoSurface->FN(), gui_currentColorPatternIndexPair.first);
}
void ReducedPatternSimulator::removeTesselatedPatterns()
{
if (pTiledReducedPattern_simulationMesh) {
pTiledReducedPattern_simulationMesh->unregister();
pTiledReducedPattern_simulationMesh.reset();
}
if (pTiledFullPattern_simulationMesh) {
pTiledFullPattern_simulationMesh->unregister();
pTiledFullPattern_simulationMesh.reset();
}
}
void ReducedPatternSimulator::reset()
{
removeDrawnSimulationJobs();
removeTesselatedPatterns();
resetUserSelectedFaces();
resetUserSelectedVertices();
polyscope::requestRedraw();
}
double ReducedPatternSimulator::computeDistance(
const SimulationResults &resultsA,
const SimulationResults &resultsB,
const std::unordered_map<VertexIndex, VertexIndex> &resultsAToResultsBViMap) const
{
double distance = 0;
for (std::pair<int, int> resultsAToResultsBViPair : resultsAToResultsBViMap) {
const double vertexToVertexDistance
= (resultsA.displacements[resultsAToResultsBViPair.first].getTranslation()
- resultsB.displacements[resultsAToResultsBViPair.second].getTranslation())
.norm();
distance += vertexToVertexDistance;
}
return distance;
}
void ReducedPatternSimulator::saveJobs(const std::filesystem::path &saveToPath)
{
std::filesystem::create_directories(saveToPath);
const auto fullJobFolderPath = std::filesystem::path(saveToPath).append(fullPatternsLabel);
std::filesystem::create_directories(fullJobFolderPath);
pJob_tiledFullPattern->save(fullJobFolderPath.string());
const auto reducedJobFolderPath = std::filesystem::path(saveToPath).append("Reduced");
std::filesystem::create_directories(reducedJobFolderPath);
pJob_tiledReducedPattern->save(reducedJobFolderPath.string());
}
void ReducedPatternSimulator::saveResults(const std::string &outputFolderPath,
SimulationResults &fullResults,
SimulationResults &reducedResults)
{
std::filesystem::create_directories(outputFolderPath);
const auto fullResultsFolderPath = std::filesystem::path(outputFolderPath).append("Full");
std::filesystem::create_directories(fullResultsFolderPath);
const auto reducedResultsFolderPath = std::filesystem::path(outputFolderPath).append("Reduced");
std::filesystem::create_directories(reducedResultsFolderPath);
fullResults.save(fullResultsFolderPath.string());
reducedResults.save(reducedResultsFolderPath.string());
}
void ReducedPatternSimulator::generateRandomSimulationScenario(
const std::array<float, 4> &randomScenarioParameters)
{
const float percentageOfRigidVertices = gui_randomnessParams[0];
const float percentageOfFixedVertices = gui_randomnessParams[1];
const float percentageOfNodalForcesVertices = gui_randomnessParams[2];
assert(percentageOfFixedVertices <= 1 && percentageOfNodalForcesVertices <= 1
&& percentageOfRigidVertices <= 1);
const float maxAbsForceMagnitude = gui_randomnessParams[3];
assert(maxAbsForceMagnitude >= 0);
std::random_device dev;
std::mt19937 rng(dev());
std::uniform_int_distribution<std::mt19937::result_type>
distInt(0, tilledFullPatternInterfaceVi.size() - 1);
std::unordered_set<size_t> usedVertices;
size_t numberOfFixedVertices = 0;
while (numberOfFixedVertices / static_cast<double>(tilledFullPatternInterfaceVi.size())
< percentageOfFixedVertices) {
const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)];
if (!usedVertices.contains(fullTilledVi)) {
pJob_tiledFullPattern->constrainedVertices[fullTilledVi] = {0, 1, 2};
numberOfFixedVertices++;
usedVertices.insert(fullTilledVi);
}
}
size_t numberOfRigidVertices = 0;
while (numberOfRigidVertices / static_cast<double>(tilledFullPatternInterfaceVi.size())
< percentageOfRigidVertices) {
const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)];
if (!usedVertices.contains(fullTilledVi)) {
pJob_tiledFullPattern->constrainedVertices[fullTilledVi] = {0, 1, 2, 3, 4, 5};
numberOfRigidVertices++;
usedVertices.insert(fullTilledVi);
}
}
std::uniform_real_distribution<> distFloat(-maxAbsForceMagnitude, maxAbsForceMagnitude);
while (pJob_tiledFullPattern->nodalExternalForces.size()
/ static_cast<double>(tilledFullPatternInterfaceVi.size())
< percentageOfNodalForcesVertices) {
const size_t fullTilledVi = tilledFullPatternInterfaceVi[distInt(rng)];
if (usedVertices.contains(fullTilledVi)) {
continue;
}
// const double inPlaneForceMagnitude = 100 * std::abs(distFloat(rng));
// const CoordType vertexPosition
// = pTiledFullPattern_simulationMesh->vert[fullTilledVi].cP().Normalize();
// Vector6d externalForce({inPlaneForceMagnitude * vertexPosition[0],
// inPlaneForceMagnitude * vertexPosition[1],
// distFloat(rng),
// 0,
// 0,
// 0});
Vector6d externalForce;
if (distFloat(rng) > 0) {
externalForce = Vector6d({0, 0, distFloat(rng), 0, 0, 0});
} else {
externalForce = Vector6d({distFloat(rng), distFloat(rng), 0, 0, 0, 0});
}
minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm());
pJob_tiledFullPattern->nodalExternalForces[fullTilledVi] = externalForce;
usedVertices.insert(fullTilledVi);
}
pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
pJob_tiledFullPattern->remap(fullToReducedViMap, *pJob_tiledReducedPattern);
pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
}
void ReducedPatternSimulator::reportDistances(
const std::tuple<DRMFullSimulationResults, ReducedSimulationResults, LinearFullSimulationResults>
&simulationResults,
const std::string &csvFileName)
{
DRMFullSimulationResults drmFullTilledResults;
ReducedSimulationResults reducedTilledResults;
LinearFullSimulationResults linearFullTilledResults;
std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults) = simulationResults;
const double distance_fullDrmToReduced = computeDistance(drmFullTilledResults,
reducedTilledResults,
fullToReducedViMap);
// std::cout << "Distance full to reduced" << distance_fullDrmToReduced<< std::endl;
double distance_fullSumOfAllVerts = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
distance_fullSumOfAllVerts
+= drmFullTilledResults.displacements[fullToReducedPair.first].getTranslation().norm();
}
// csvFile csv_results({}, false);
// csv_results << "Full Pattern Names"
// << "Simulation job label"
// << "Distance";
// csv_results.endrow();
// csv_results << fullPatternsSurfacelabel << gui_jobLabel << std::to_string(distance);
// csv_results.endrow();
std::unordered_map<VertexIndex, VertexIndex> fullToFullDummyViMap = fullToReducedViMap;
for (auto &viPair : fullToFullDummyViMap) {
viPair.second = viPair.first;
}
const double distance_fullLinearToDrm = computeDistance(drmFullTilledResults,
linearFullTilledResults,
fullToFullDummyViMap)
/*/ (2 * surfaceBaseTriangleHeight * fullToReducedViMap.size())*/;
// std::cout << "Distance tilled full linear vs drm:" << distance_fullLinearToDrm << std::endl;
// std::cout << "Distance full to reduced norm:"
// << distance_fullTilledToReducedTilled / distance_fullSumOfAllVerts << std::endl;
// std::cout << "Distance tilled full linear vs drm norm:"
// << distance_fullLinearToDrm / distance_fullSumOfAllVerts << std::endl;
// Write results in csv
csvFile csv_results(csvFileName, false);
// csvFile csv_results(std::filesystem::path(dirPath_thisOptimization)
// .append("results.csv")
// .string(),
// false);
csv_results << "Job Label"
<< "drm2Reduced"
<< "drm2Linear"
<< "norm_drm2Reduced"
<< "norm_drm2Linear";
csv_results << endrow;
csv_results << pJob_tiledFullPattern->getLabel() << distance_fullDrmToReduced
<< distance_fullLinearToDrm
<< distance_fullDrmToReduced / distance_fullSumOfAllVerts
<< distance_fullLinearToDrm / distance_fullSumOfAllVerts;
csv_results << endrow;
}
void ReducedPatternSimulator::reportDistances(const std::vector<std::string> &scenarioLabels)
{
for (const std::string &scenario : scenarioLabels) {
reset();
gui_jobLabel = scenario;
// generateRandomSimulationScenario(gui_randomnessParams);
loadScenario(gui_jobLabel);
DRMFullSimulationResults drmFullTilledResults;
ReducedSimulationResults reducedTilledResults;
LinearFullSimulationResults linearFullTilledResults;
std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults)
= runAllSimulations(gui_jobLabel, gui_shouldRerunFullPatternSimulation, false);
reportDistances({drmFullTilledResults, reducedTilledResults, linearFullTilledResults},
"distances.csv");
}
}
std::tuple<DRMFullSimulationResults, ReducedSimulationResults, LinearFullSimulationResults>
ReducedPatternSimulator::runAllSimulations(const std::string &jobLabel,
const bool &shouldRerunDRMFullPatternSimulation,
const bool &showReducedPatternResultsBeforeDRM)
{
std::cout << "Executing scenario:" << jobLabel << std::endl;
std::pair<DRMFullSimulationResults, ReducedSimulationResults> results
= runSimulations(jobLabel,
shouldRerunDRMFullPatternSimulation,
showReducedPatternResultsBeforeDRM);
if (results.first.converged == false || results.second.converged == false) {
std::cerr << "Simulation did not converge" << std::endl;
return std::tuple<DRMFullSimulationResults,
ReducedSimulationResults,
LinearFullSimulationResults>();
}
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullTilledLinearModel
= linearSimulationModel.executeSimulation(pJob_tiledFullPattern);
simulationResults_fullTilledLinearModel.setLabelPrefix("LinearSimModel");
return std::make_tuple(results.first, results.second, simulationResults_fullTilledLinearModel);
}
void ReducedPatternSimulator::saveTesselation(const std::filesystem::path &saveTo)
{
std::filesystem::create_directories(saveTo);
nlohmann::json json;
const std::string jsonLabel_tesselation{"faceToPatternIndices"};
json[jsonLabel_tesselation] = gui_faceToPatternIndex;
std::filesystem::path jsonFilePath(
std::filesystem::path(saveTo).append(fullPatternsSurfacelabel + ".json"));
std::ofstream jsonFile(jsonFilePath.string());
jsonFile << json;
}
void ReducedPatternSimulator::computeSurfaceColorsFromPerFacePatterns(
const std::vector<int> &faceToPatternIndex)
{
std::set uniquePatternSet(faceToPatternIndex.begin(), faceToPatternIndex.end());
const int numberOfDistinctPatterns = uniquePatternSet.size();
std::vector<glm::vec3> surfaceColors(numberOfDistinctPatterns);
for (int colorIndex = 0; colorIndex < numberOfDistinctPatterns; colorIndex++) {
surfaceColors[colorIndex] = polyscope::getNextUniqueColor();
}
// gui_currentColorPatternIndexPair{glm::vec3(1, 0, 0), 0};
for (int fi = 0; fi < pTileIntoSurface->FN(); fi++) {
const int patternIndex = gui_faceToPatternIndex[fi];
// if (patternIndex != gui_currentColorPatternIndexPair.second) {
// gui_currentColorPatternIndexPair.first = surfaceColors[patternIndex];
// gui_current
// }
gui_colorsPerFace[fi] = surfaceColors[patternIndex];
}
}
void ReducedPatternSimulator::loadTessellation(const std::filesystem::path &jsonFilePath)
{
nlohmann::json json;
std::ifstream ifs(jsonFilePath);
ifs >> json;
const std::string jsonLabel_tesselation{"faceToPatternIndices"};
gui_faceToPatternIndex = json[jsonLabel_tesselation].get<std::vector<int>>();
computeSurfaceColorsFromPerFacePatterns(gui_faceToPatternIndex);
}
void ReducedPatternSimulator::updateTesselationColors()
{
computeSurfaceColorsFromPerFacePatterns(gui_faceToPatternIndex);
auto tileIntoSurface_polyscopeHandle = polyscope::getSurfaceMesh(pTileIntoSurface->getLabel());
tileIntoSurface_polyscopeHandle->addFaceColorQuantity("Selected faces", gui_colorsPerFace)
->setEnabled(true);
}
void ReducedPatternSimulator::createGuiMenu()
{
// std::shared_ptr<SimulationMesh> &pTiledFullPattern_simulationMesh = pJob_tiledFullPattern->pMesh;
// assert(pTiledFullPattern_simulationMesh != nullptr);
// std::shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh = pJob_tiledReducedPattern
// ->pMesh;
// assert(pTiledReducedPattern_simulationMesh != nullptr);
PolyscopeInterface::addUserCallback([&]() {
// ImGuiIO &io = ImGui::GetIO();
if (ImGui::Button("Fix vertex")) {
if (!gui_fullPatternSelectedVertices.empty()) {
for (const int fullVi : gui_fullPatternSelectedVertices) {
pJob_tiledFullPattern->constrainedVertices[fullVi] = {0, 1, 2, 3, 4, 5};
pJob_tiledReducedPattern->constrainedVertices[fullToReducedViMap.at(fullVi)]
= {0, 1, 2, 3, 4, 5};
}
resetUserSelectedVertices();
} else {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string &pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {0, 1, 2, 3, 4, 5};
pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second]
= {0, 1, 2, 3, 4, 5};
}
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
static int gui_dof = 0;
ImGui::InputInt("DoF", &gui_dof);
if (ImGui::Button("Constrain DoF of vertex")) {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string &pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {gui_dof};
pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second] = {gui_dof};
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
if (ImGui::Button("Constrain translation of vertex")) {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string &pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->constrainedVertices[pickedVertices.first] = {0, 1, 2};
pJob_tiledReducedPattern->constrainedVertices[pickedVertices.second] = {0, 1, 2};
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
if (ImGui::Button("Remove constrains on vertex")) {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string &pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->constrainedVertices.erase(pickedVertices.first);
pJob_tiledReducedPattern->constrainedVertices.erase(pickedVertices.second);
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
ImGui::InputFloat("Scale factor", &gui_scaleFactor);
if (ImGui::Button("Scale forces/moments")) {
if (!gui_fullPatternSelectedVertices.empty()) {
for (const int fullVi : gui_fullPatternSelectedVertices) {
pJob_tiledFullPattern->nodalExternalForces[fullVi]
= pJob_tiledFullPattern->nodalExternalForces[fullVi] * gui_scaleFactor;
const int reducedVi = fullToReducedViMap.at(fullVi);
pJob_tiledReducedPattern->nodalExternalForces[reducedVi]
= pJob_tiledReducedPattern->nodalExternalForces[reducedVi]
* gui_scaleFactor;
}
resetUserSelectedVertices();
} else {
for (auto &externalLoad : pJob_tiledFullPattern->nodalExternalForces) {
externalLoad.second = externalLoad.second * gui_scaleFactor;
}
for (auto &externalLoad : pJob_tiledReducedPattern->nodalExternalForces) {
externalLoad.second = externalLoad.second * gui_scaleFactor;
}
}
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
if (ImGui::Button("Set all forces/moments")) {
const Vector6d desiredExternalLoad({gui_externalForce[0],
gui_externalForce[1],
gui_externalForce[2],
gui_externalMoment[0],
gui_externalMoment[1],
gui_externalMoment[2]});
for (auto &externalLoad : pJob_tiledFullPattern->nodalExternalForces) {
externalLoad.second = desiredExternalLoad;
}
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
for (auto &externalLoad : pJob_tiledReducedPattern->nodalExternalForces) {
externalLoad.second = desiredExternalLoad;
}
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
minInputForceMagnitude = std::min(minInputForceMagnitude, externalForce.norm());
polyscope::requestRedraw();
}
if (ImGui::Button("Add force/moment")) {
Vector6d externalForce({gui_externalForce[0],
gui_externalForce[1],
gui_externalForce[2],
gui_externalMoment[0],
gui_externalMoment[1],
gui_externalMoment[2]});
if (!gui_fullPatternSelectedVertices.empty()) {
for (const int fullVi : gui_fullPatternSelectedVertices) {
pJob_tiledFullPattern->nodalExternalForces[fullVi] = externalForce;
pJob_tiledReducedPattern->nodalExternalForces[fullToReducedViMap.at(fullVi)]
= externalForce;
}
resetUserSelectedVertices();
} else {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string &pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->nodalExternalForces[pickedVertices.first] = externalForce;
pJob_tiledReducedPattern->nodalExternalForces[pickedVertices.second] = externalForce;
}
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
ImGui::InputFloat3("Force", gui_externalForce.data());
ImGui::InputFloat3("Moment", gui_externalMoment.data());
if (ImGui::Button("Remove Force/Moment")) {
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
const std::string pickedStructureName = selection.first;
if (pickedStructureName.empty()) {
return;
}
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
pJob_tiledFullPattern->nodalExternalForces.erase(pickedVertices.first);
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern->nodalExternalForces.erase(pickedVertices.second);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
polyscope::requestRedraw();
}
ImGui::InputText("Job label", &gui_jobLabel);
if (ImGui::Button("Load job")) {
// reset();
removeDrawnSimulationJobs();
resetUserSelectedFaces();
resetUserSelectedVertices();
loadScenario(gui_jobLabel);
polyscope::requestRedraw();
}
ImGui::InputText("Tesselation label", &gui_tessellationLabel);
if (ImGui::Button("Load tesselation")) {
if (gui_tessellationLabel.empty()) {
return;
}
reset();
// const std::string fullPatternsOrderLabel = gui_tessellationLabel;
// const std::string tessellationFilename = computeFullPatternSetLabel(mOptimizationResults)
// + "_" + fullPatternsOrderLabel + "_"
// + pTileIntoSurface->getLabel();
const std::string tessellationFilename = gui_tessellationLabel;
std::filesystem::path tessellationJsonPath(
std::filesystem::path(
"/home/iason/Coding/build/PatternTillingReducedModel/Scenarios")
.append(pTileIntoSurface->getLabel())
.append("Tessellations")
.append(tessellationFilename + ".json"));
loadTessellation(tessellationJsonPath);
updateTesselationColors();
resetTilledMeshes();
polyscope::requestRedraw();
}
if (ImGui::Button("Random tesselation")) {
reset();
gui_faceToPatternIndex = computePerSurfaceFacePatternsIndices(mOptimizationResults);
updateTesselationColors();
resetTilledMeshes();
}
ImGui::Checkbox("Re-run full pattern simulation", &gui_shouldRerunFullPatternSimulation);
if (ImGui::Button("Run Simulations")) {
if (pJob_tiledFullPattern->isEmpty() || pJob_tiledFullPattern->isEmpty()) {
std::cerr
<< "Empty simulation job. Load or create a simulation job to run a simulation."
<< std::endl;
return;
}
DRMFullSimulationResults drmFullTilledResults;
ReducedSimulationResults reducedTilledResults;
LinearFullSimulationResults linearFullTilledResults;
std::tie(drmFullTilledResults, reducedTilledResults, linearFullTilledResults)
= runAllSimulations(gui_jobLabel, gui_shouldRerunFullPatternSimulation, true);
if (!drmFullTilledResults.converged || !reducedTilledResults.converged
|| !linearFullTilledResults.converged) {
std::cerr << "A simulation did not converge." << std::endl;
return;
}
reportDistances({drmFullTilledResults, reducedTilledResults, linearFullTilledResults});
drmFullTilledResults.registerForDrawing();
reducedTilledResults.registerForDrawing();
linearFullTilledResults.registerForDrawing();
polyscope::requestRedraw();
}
ImGui::InputFloat4("Randomness params", gui_randomnessParams.data());
if (ImGui::Button("Generate random scenario")) {
gui_jobLabel = "22Hex_random";
gui_shouldRerunFullPatternSimulation = true;
reset();
generateRandomSimulationScenario(gui_randomnessParams);
pJob_tiledFullPattern->label = gui_jobLabel;
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
true);
pJob_tiledReducedPattern
->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(), true);
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullTilledLinearModel
= linearSimulationModel.executeSimulation(pJob_tiledFullPattern);
simulationResults_fullTilledLinearModel.setLabelPrefix("LinearSimModel");
simulationResults_fullTilledLinearModel.registerForDrawing();
polyscope::requestRedraw();
}
if (ImGui::Button("Reset")) {
reset();
}
// Brush
ImGuiIO &io = ImGui::GetIO();
if (io.KeyShift && polyscope::render::engine->isKeyPressed('b')) {
ImVec2 p = ImGui::GetMousePos();
std::pair<polyscope::Structure *, size_t> pickResult
= polyscope::pick::evaluatePickQuery(io.DisplayFramebufferScale.x * p.x,
io.DisplayFramebufferScale.y * p.y);
if (pickResult.first == nullptr) {
return;
}
polyscope::pick::setSelection(pickResult);
auto tileIntoSurface_polyscopeHandle = polyscope::getSurfaceMesh(
pTileIntoSurface->getLabel());
const bool pickedFace = pickResult.first->name == tileIntoSurface_polyscopeHandle->name
&& pickResult.second
>= tileIntoSurface_polyscopeHandle->facePickIndStart
&& pickResult.second
< tileIntoSurface_polyscopeHandle->edgePickIndStart;
if (pickedFace) {
const size_t faceIndex = pickResult.second
- tileIntoSurface_polyscopeHandle->facePickIndStart;
// std::cout << "Picked face:" << faceIndex << std::endl;
gui_colorsPerFace[faceIndex] = gui_currentColorPatternIndexPair.first;
gui_faceToPatternIndex[faceIndex] = gui_currentColorPatternIndexPair.second;
std::vector<size_t> tileIntoEiToTiledFullPatternVi;
fillFullPatterns(mOptimizationResults);
tileFullPatterns(fullPatterns,
gui_faceToPatternIndex,
pTiledFullPattern_simulationMesh,
tileIntoEiToTiledFullPatternVi);
///Create simulation jobs
// pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
// pJob_tiledFullPattern->label="TiledFull";
// pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
// pJob_tiledReducedPattern->setLabel("TiledReduced");
// createSimulationJobs(externalForce, *pJob_tiledFullPattern, *pJob_tiledReducedPattern);
// pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
if (pTiledFullPattern_simulationMesh
&& polyscope::hasStructure(polyscope::CurveNetwork::structureTypeName,
pTiledFullPattern_simulationMesh->getLabel())) {
bool shouldEnableFullPattern = polyscope::getCurveNetwork(
pTiledFullPattern_simulationMesh->getLabel())
->isEnabled();
pTiledFullPattern_simulationMesh
->registerForDrawing(color_tesselatedFullPatterns)
->setEnabled(shouldEnableFullPattern);
} else {
pTiledFullPattern_simulationMesh
->registerForDrawing(color_tesselatedFullPatterns)
->setEnabled(true);
}
// pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
// pTiledReducedPattern_simulationMesh->registerForDrawing();
}
tileIntoSurface_polyscopeHandle
->addFaceColorQuantity("Selected faces", gui_colorsPerFace)
->setEnabled(true);
polyscope::requestRedraw();
}
if (io.KeyShift && polyscope::render::engine->isKeyPressed('n')) {
gui_currentColorPatternIndexPair.first = polyscope::getNextUniqueColor();
gui_currentColorPatternIndexPair.second++;
}
if (io.KeyShift && polyscope::render::engine->isKeyPressed('t')) {
resetTilledMeshes();
}
if (io.KeyCtrl && io.KeyShift) {
if (polyscope::pick::haveSelection()
&& (!pTiledFullPattern_simulationMesh || !pTiledReducedPattern_simulationMesh)) {
std::cerr << "Tiled patterns not registered. Tessellated patterns are required for "
"selecting a vertex."
<< std::endl;
polyscope::pick::resetSelection();
return;
}
// ImVec2 p = ImGui::GetMousePos();
// std::pair<polyscope::Structure *, size_t> pickResult
// = polyscope::pick::evaluatePickQuery(io.DisplayFramebufferScale.x * p.x,
// io.DisplayFramebufferScale.y * p.y);
// if (pickResult.first == nullptr) {
// return;
// }
const std::pair<std::string, size_t> selection = PolyscopeInterface::getSelection();
std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex> pickedVertices
= getPickedInterfaceVertices(selection);
if (pickedVertices.first == -1 && pickedVertices.second == -1) {
return;
}
gui_fullPatternSelectedVertices.push_back(pickedVertices.first);
gui_fullVerticesColors[pickedVertices.first] = glm::vec3(1, 0, 0);
auto tiledFullPattern_polyscopeHandle = polyscope::getCurveNetwork(
pTiledFullPattern_simulationMesh->getLabel());
tiledFullPattern_polyscopeHandle
->addNodeColorQuantity("Selected Vertices", gui_fullVerticesColors)
->setEnabled(true);
gui_reducedVerticesColors[pickedVertices.second] = glm::vec3(1, 0, 0);
auto tiledReducedPattern_polyscopeHandle = polyscope::getCurveNetwork(
pTiledReducedPattern_simulationMesh->getLabel());
tiledReducedPattern_polyscopeHandle
->addNodeColorQuantity("Selected Vertices", gui_reducedVerticesColors)
->setEnabled(true);
polyscope::requestRedraw();
}
});
}
void ReducedPatternSimulator::fillFullPatterns(
std::vector<ReducedPatternOptimization::Results> &optimizationResults)
{
for (size_t fullPatternIndex = 0; fullPatternIndex < optimizationResults.size();
fullPatternIndex++) {
// const vcg::Triangle3<double> fullPatternBaseTriangle = fullPattern.computeBaseTriangle();
PatternGeometry &baseTriangleFullPattern = optimizationResults[fullPatternIndex]
.baseTriangleFullPattern;
fullPatterns[fullPatternIndex].copy(baseTriangleFullPattern);
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<VCGEdgeMesh::VertexPointer> pu;
fullPatterns[fullPatternIndex].deleteDanglingVertices(pu);
if (!pu.remap.empty()) {
fullPatterns[fullPatternIndex].interfaceNodeIndex = pu.remap[3];
} else {
fullPatterns[fullPatternIndex].interfaceNodeIndex = 3;
}
}
}
void ReducedPatternSimulator::createTiledSimulationMeshes(
const std::shared_ptr<ConstVCGPolyMesh> &pTileIntoSurface,
std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<OptimizationResultsIndex> &perSurfaceFacePatternIndices)
{
std::vector<size_t> tileIntoEdgeToTiledFullPattern;
createTiledSimulationMeshes(pTileIntoSurface,
optimizationResults,
perSurfaceFacePatternIndices,
tileIntoEdgeToTiledFullPattern);
}
void ReducedPatternSimulator::createTiledSimulationMeshes(
const std::shared_ptr<ConstVCGPolyMesh> &pTileIntoSurface,
std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<OptimizationResultsIndex> &perSurfaceFacePatternIndices,
std::vector<size_t> &tileIntoEdgeToTiledFullPattern)
{
assert(perSurfaceFacePatternIndices.size() == pTileIntoSurface->FN());
fillFullPatterns(optimizationResults);
tileFullPatterns(fullPatterns,
perSurfaceFacePatternIndices,
pTiledFullPattern_simulationMesh,
tileIntoEdgeToTiledFullPattern);
std::vector<size_t> tileIntoEdgeToTiledReducedPattern;
tileReducedPatterns(optimizationResults,
perSurfaceFacePatternIndices,
pTiledReducedPattern_simulationMesh,
tileIntoEdgeToTiledReducedPattern);
constructViMaps(tileIntoEdgeToTiledFullPattern, tileIntoEdgeToTiledReducedPattern);
computeLabels(optimizationResults, perSurfaceFacePatternIndices);
// std::vector<PatternGeometry> reducedPatterns(optimizationResults.size());
// for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size();
// reducedPatternIndex++) {
// reducedPatterns[reducedPatternIndex].copy(reducedPattern);
// ReducedPatternOptimization::Results::applyOptimizationResults_innerHexagon(
// optimizationResults[reducedPatternIndex],
// reducedPatternBaseTriangle,
// reducedPatterns[reducedPatternIndex]);
// }
// std::vector<OptimizationResultsIndex> randomShufflePerSurfaceFacePatternIndices
// = perSurfaceFacePatternIndices;
// std::random_shuffle(randomShufflePerSurfaceFacePatternIndices.begin(),
// randomShufflePerSurfaceFacePatternIndices.end());
// std::vector<int> tileIntoEdgeToTiledReducedPattern;
// std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
// std::shared_ptr<PatternGeometry> pTiledReducedPattern = PatternGeometry::tilePattern(
// reducedPatterns,
// {0},
// *pTileIntoSurface,
// // randomShufflePerSurfaceFacePatternIndices,
// perSurfaceFacePatternIndices,
// tileIntoEdgeToTiledReducedPattern,
// perPatternIndexTiledReducedPatternEdgeIndices);
// pTiledReducedPattern->setLabel("Tiled_" + reducedPatterns[0].getLabel());
//#ifdef POLYSCOPE_DEFINED
// pTiledReducedPattern->registerForDrawing();
//#endif
// pTiledReducedPattern_simulationMesh = std::make_shared<SimulationMesh>(*pTiledReducedPattern);
// for (size_t reducedPatternIndex = 0; reducedPatternIndex < optimizationResults.size();
// reducedPatternIndex++) {
// const std::vector<size_t> &tiledPatternElementIndicesForReducedPattern
// = perPatternIndexTiledReducedPatternEdgeIndices[reducedPatternIndex];
// applyOptimizationResults_elements(optimizationResults[reducedPatternIndex],
// pTiledReducedPattern_simulationMesh,
// tiledPatternElementIndicesForReducedPattern);
// }
// pTiledReducedPattern_simulationMesh->reset();
// for (int ei = 0; ei < pTileIntoSurface->EN(); ei++) {
// fullToReducedViMap[tileIntoEdgeToTiledFullPattern[ei]]
// = tileIntoEdgeToTiledReducedPattern[ei];
// }
// constructInverseMap(fullToReducedViMap, reducedToFullViMap);
}
void ReducedPatternSimulator::loadScenario(const string &scenarioLabel)
{
std::filesystem::path simulationJobsInputFolderPath(
std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/Scenarios")
.append(scenarioLabel)
.append("SimulationJobs"));
if (!std::filesystem::exists(simulationJobsInputFolderPath)
|| !std::filesystem::is_directory(simulationJobsInputFolderPath)) {
simulationJobsInputFolderPath
= std::filesystem::path("/home/iason/Coding/build/PatternTillingReducedModel/Scenarios")
.append(surfaceLabel)
.append(scenarioLabel)
.append("SimulationJobs");
}
if (!std::filesystem::is_directory(simulationJobsInputFolderPath)) {
std::cerr << "The given path is not a directory:" << simulationJobsInputFolderPath
<< std::endl;
return;
}
// //Full
// std::filesystem::path fullPatternSimulationJobFilePath;
// for (auto &p : std::filesystem::recursive_directory_iterator(
// std::filesystem::path(inputFolderPath).append("Full"))) {
// if (p.path().extension() == ".json")
// fullPatternSimulationJobFilePath = p.path();
// }
// assert(pTiledFullPattern_simulationMesh != nullptr);
// pJob_tiledFullPattern->load(fullPatternSimulationJobFilePath.string());
// // pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
// assert(pTiledFullPattern_simulationMesh != nullptr);
// pTiledFullPattern_simulationMesh->registerForDrawing();
// pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(),
// true);
// for (const auto &externalLoad : pJob_tiledFullPattern->nodalExternalForces) {
// minInputForceMagnitude = std::min(minInputForceMagnitude,
// externalLoad.second.norm());
// }
//Reduced
std::filesystem::path reducedPatternSimulationJobFilePath;
for (auto &p : std::filesystem::directory_iterator(std::filesystem::path(
std::filesystem::path(simulationJobsInputFolderPath).append("Reduced")))) {
if (p.path().extension() == ".json")
reducedPatternSimulationJobFilePath = p.path();
}
pJob_tiledReducedPattern->load(reducedPatternSimulationJobFilePath.string(), false);
//NOTE: Is there a way to not apply the optimization results again? I only need the Job and not the mesh
// pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
// applyOptimizationResults_elements(optimizationResults,
// pTiledReducedPattern_simulationMesh);
pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
assert(pTiledReducedPattern_simulationMesh != nullptr);
// pTiledReducedPattern_simulationMesh->registerForDrawing();
pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern_simulationMesh->getLabel(),
true);
pJob_tiledFullPattern->label = scenarioLabel;
pJob_tiledReducedPattern->label = scenarioLabel;
//Apply the reduced job to the full pattern
pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
pJob_tiledReducedPattern->remap(reducedToFullViMap, *pJob_tiledFullPattern);
pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern_simulationMesh->getLabel(), true);
}
std::string ReducedPatternSimulator::computeFullPatternSetLabel(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults)
{
std::string fullPatternSetLabel;
for (OptimizationResultsIndex optResIndex = 0; optResIndex < optimizationResults.size();
optResIndex++) {
fullPatternSetLabel += optimizationResults[optResIndex].baseTriangleFullPattern.getLabel()
+ "_";
}
fullPatternSetLabel.pop_back();
return fullPatternSetLabel;
}
void ReducedPatternSimulator::computeLabels(
const std::vector<ReducedPatternOptimization::Results> &optimizationResults,
const std::vector<OptimizationResultsIndex> &perSurfaceFacePatternIndex)
{
fullPatternsLabel.clear();
const std::string fullPatternsSetLabel = computeFullPatternSetLabel(optimizationResults);
const std::string fullPatternsOrderLabel = std::to_string(
computeHashOrdered(perSurfaceFacePatternIndex));
fullPatternsLabel = fullPatternsSetLabel + "_" + fullPatternsOrderLabel;
// fullPatternsLabel += "_" + std::to_string(optimizationResults[0].settings.objectiveWeights.translational);
surfaceLabel = pTileIntoSurface->getLabel();
fullPatternsSurfacelabel = fullPatternsLabel + "_" + surfaceLabel;
}
void ReducedPatternSimulator::resetTilledMeshes()
{
std::vector<size_t> tileIntoEiToTiledFullPatternVi;
createTiledSimulationMeshes(
pTileIntoSurface,
mOptimizationResults,
gui_faceToPatternIndex,
// std::vector<size_t>(tileIntoSurface->FN(), 0),
tileIntoEiToTiledFullPatternVi);
///Create simulation jobs
pJob_tiledFullPattern = std::make_shared<SimulationJob>(SimulationJob());
// pJob_tiledFullPattern->label="TiledFull";
pJob_tiledReducedPattern = std::make_shared<SimulationJob>(SimulationJob());
// pJob_tiledReducedPattern->setLabel("TiledReduced");
// createSimulationJobs(externalForce, *pJob_tiledFullPattern, *pJob_tiledReducedPattern);
pJob_tiledFullPattern->pMesh = pTiledFullPattern_simulationMesh;
pTiledFullPattern_simulationMesh->registerForDrawing(color_tesselatedFullPatterns);
pJob_tiledReducedPattern->pMesh = pTiledReducedPattern_simulationMesh;
pTiledReducedPattern_simulationMesh->registerForDrawing(color_tesselatedReducedPatterns);
gui_fullVerticesColors.resize(pTiledFullPattern_simulationMesh->VN());
gui_reducedVerticesColors.resize(pTiledReducedPattern_simulationMesh->VN());
}
ReducedPatternSimulator::ReducedPatternSimulator(
std::vector<ReducedPatternOptimization::Results> &optimizationResults)
: fullPatterns(optimizationResults.size()), mOptimizationResults(optimizationResults)
{}
void ReducedPatternSimulator::simulate(
std::shared_ptr<VCGPolyMesh> &tileIntoSurface,
std::vector<ReducedPatternOptimization::Results> &optimizationResults,
PatternGeometry &reducedPattern,
const std::vector<OptimizationResultsIndex> &inputPerSurfaceFacePatternIndices)
{
reducedPatternBaseTriangle = reducedPattern.computeBaseTriangle();
this->reducedPattern.copy(reducedPattern);
this->reducedPattern.deleteDanglingVertices();
this->reducedPattern.interfaceNodeIndex = 1;
this->pTileIntoSurface = tileIntoSurface;
///Preprocess tile into surface
/// NOTE: I am assuming here that the base triangle is the same for all optimization results
const double optimizedBaseTriangleHeight
= vcg::Distance(optimizationResults[0].baseTriangle.cP(0),
(optimizationResults[0].baseTriangle.cP(1)
+ optimizationResults[0].baseTriangle.cP(2))
/ 2);
preprocessSurface(optimizedBaseTriangleHeight, *pTileIntoSurface);
// vcg::tri::io::ExporterOBJ<VCGPolyMesh>::Save(*pTileIntoSurface,
// "tileIntoSurface.obj",
// vcg::tri::io::Mask::IOM_BITPOLYGONAL);
pTileIntoSurface->registerForDrawing()->setEdgeWidth(1);
// TODO:REMOVE
reset();
gui_faceToPatternIndex = computePerSurfaceFacePatternsIndices(mOptimizationResults);
resetTilledMeshes();
removeDrawnSimulationJobs();
loadScenario("ex1.1");
gui_colorsPerFace.resize(pTileIntoSurface->FN(), glm::vec3(0, 1, 0));
gui_faceToPatternIndex.resize(pTileIntoSurface->FN(), -1);
createGuiMenu();
polyscope::show();
}
void ReducedPatternSimulator::preprocessSurface(const double &desiredBaseTriangleHeight,
VCGPolyMesh &surface)
{
surface.moveToCenter();
// const double desiredAverageBaseTriangleSize = des; //TODO: obtain from optimization results
vcg::tri::UpdatePosition<VCGPolyMesh>::Scale(surface,
desiredBaseTriangleHeight
/ surface.getAverageFaceRadius());
surfaceBaseTriangleHeight = desiredBaseTriangleHeight;
}