Uses chronos simulation model. Refactoring

This commit is contained in:
iasonmanolas 2022-07-19 14:59:05 +03:00
parent ef18646cfd
commit 7a0124155c
4 changed files with 90 additions and 74 deletions

View File

@ -331,19 +331,19 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
pTiledFullPattern->setLabel("Tiled_full_patterns");
// pTiledFullPattern->registerForDrawing();
//Tile reduced pattern into surface
// PatternGeometry reducedPattern;
ReducedModel reducedModel;
// reducedModel.registerForDrawing();
// polyscope::show();
reducedModel.deleteDanglingVertices();
// reducedPattern.interfaceNodeIndex = 1;
// assert(reducedPattern.interfaceNodeIndex == 1);
std::unordered_map<std::string, double>
optimalXVariables_set(optimizationResult.optimalXNameValuePairs.begin(),
optimizationResult.optimalXNameValuePairs.end());
reducedModel.updateBaseTriangleGeometry_R(optimalXVariables_set.at("R"));
reducedModel.updateBaseTriangleGeometry_theta(optimalXVariables_set.at("Theta"));
// const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle();
// ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned(
// optimizationResult, reducedPatternBaseTriangle, reducedModel);
std::vector<PatternGeometry> reducedPatterns(1);
reducedPatterns[0].copy(reducedModel);
const auto reducedPatternBaseTriangle = reducedModel.computeBaseTriangle();
ReducedModelOptimization::Results::applyOptimizationResults_reducedModel_nonFanned(
optimizationResult, reducedPatternBaseTriangle, reducedPatterns[0]);
std::vector<std::vector<size_t>> perPatternIndexTiledReducedPatternEdgeIndices;
std::vector<size_t> tileIntoEdgeToTiledReducedVi;
std::shared_ptr<PatternGeometry> pTiledReducedPattern
@ -354,8 +354,9 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
tileIntoEdgeToTiledReducedVi,
perPatternIndexTiledReducedPatternEdgeIndices);
pTiledReducedPattern->setLabel("Tiled_reduced_patterns");
// pTiledReducedPattern->registerForDrawing();
// polyscope::show();
#ifdef POLYSCOPE_DEFINED
pTiledReducedPattern->registerForDrawing();
#endif
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
fullToReducedViMap; //of only the common vertices
@ -379,9 +380,10 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
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
constexpr double beamWidth = 0.001;
constexpr double beamHeight = 0.001;
pTiledFullPattern_simulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight});
// const double beamWidth = optimizationResult.settings.patternBeamDimensions.getWidth();
// const double beamHeight = optimizationResult.settings.patternBeamDimensions.getHeight();
pTiledFullPattern_simulationMesh->setBeamCrossSection(
optimizationResult.settings.patternBeamDimensions);
if (optimizationResult.fullPatternYoungsModulus == 0) {
std::cerr << "Full pattern's young modulus not found." << std::endl;
std::terminate();
@ -419,7 +421,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
#endif
constexpr bool shouldRerunFullPatternSimulation = false;
enum PatternSimulationModelTag { DRM, Chronos };
const PatternSimulationModelTag simulationModelTag{DRM};
const PatternSimulationModelTag simulationModelTag{Chronos};
const std::string simulationModelLabel = [&]() {
switch (simulationModelTag) {
case DRM:
@ -431,7 +433,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
// for (int jobIndex = 0; jobIndex < scenariosTestSetLabels.size(); jobIndex++) {
std::for_each(
//#ifndef POLYSCOPE_DEFINED
std::execution::par_unseq,
// std::execution::par_unseq,
//#endif
scenariosTestSetLabels.begin(),
scenariosTestSetLabels.end(),
@ -460,9 +462,9 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
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();
// pJob_tiledReducedPattern->registerForDrawing(pTiledReducedPattern->getLabel());
// pJob_tiledFullPattern->registerForDrawing(pTiledFullPattern->getLabel());
// polyscope::show();
const std::filesystem::path surfaceFolderPath
= std::filesystem::path(fullPatternTessellatedResultsDirectoryPath)
.append(simulationModelLabel + "_" + pTileIntoSurface->getLabel());
@ -501,12 +503,12 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
scenarioDirectoryPath)
.append(patternLabel)
.append("SimulationJob");
SimulationResults simulationResults_tiledPattern;
SimulationResults simulationResults_tilledPattern;
if (std::filesystem::exists(tilledPatternResultsFolderPath)) {
//Load full pattern results
assert(std::filesystem::exists(fullPatternJobFolderPath));
simulationResults_tiledPattern.load(tilledPatternResultsFolderPath,
fullPatternJobFolderPath);
simulationResults_tilledPattern.load(tilledPatternResultsFolderPath,
fullPatternJobFolderPath);
//#ifdef POLYSCOPE_DEFINED
// std::array<double, 3> resultsColor({28.0, 99.0, 227.0});
// simulationResults_tiledFullPattern_drm.registerForDrawing(resultsColor);
@ -532,7 +534,7 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
// simulationResults_tiledFullPattern_drm.unregister();
//#endif
simulationResults_tiledPattern.converged = true;
simulationResults_tilledPattern.converged = true;
} else {
std::cout << "Tilled pattern simulation results not found in:"
<< tilledPatternResultsFolderPath << std::endl;
@ -541,20 +543,20 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
switch (simulationModelTag) {
case DRM: {
DRMSimulationModel drmSimulationModel;
simulationResults_tiledPattern = drmSimulationModel
.executeSimulation(pJob_tiledFullPattern,
drmSimulationSettings);
simulationResults_tilledPattern = drmSimulationModel
.executeSimulation(pJob_tiledFullPattern,
drmSimulationSettings);
break;
}
case Chronos: {
ChronosEulerSimulationModel chronosSimulationModel;
simulationResults_tiledPattern = chronosSimulationModel.executeSimulation(
simulationResults_tilledPattern = chronosSimulationModel.executeSimulation(
pJob_tiledFullPattern);
break;
}
}
}
if (!simulationResults_tiledPattern.converged) {
if (!simulationResults_tilledPattern.converged) {
std::cerr << "Full pattern simulation failed." << std::endl;
std::cerr << "Not saving results" << std::endl;
// continue;
@ -563,22 +565,25 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
std::filesystem::create_directories(tilledPatternResultsFolderPath);
const std::filesystem::path drmResultsOutputPath
= std::filesystem::path(scenarioDirectoryPath).append(patternLabel);
simulationResults_tiledPattern.save(drmResultsOutputPath);
simulationResults_tilledPattern.save(drmResultsOutputPath);
LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_tiledReducedPattern
SimulationResults simulationResults_tiledReducedModel
= linearSimulationModel.executeSimulation(pJob_tiledReducedPattern);
// simulationResults_tiledReducedPattern.registerForDrawing();
// simulationResults_tiledFullPattern_drm.registerForDrawing();
// polyscope::show();
#ifdef POLYSCOPE_DEFINED
simulationResults_tiledReducedModel.registerForDrawing();
simulationResults_tilledPattern.registerForDrawing();
polyscope::show();
#endif
//compute the full2reduced distance
const double distance_patternToReduced = simulationResults_tiledPattern.computeDistance(
simulationResults_tiledReducedPattern, fullToReducedViMap);
const double distance_patternToReduced
= simulationResults_tilledPattern
.computeDistance(simulationResults_tiledReducedModel, fullToReducedViMap);
double distance_patternSumOfAllVerts = 0;
for (std::pair<size_t, size_t> fullToReducedPair : fullToReducedViMap) {
distance_patternSumOfAllVerts += simulationResults_tiledPattern
distance_patternSumOfAllVerts += simulationResults_tilledPattern
.displacements[fullToReducedPair.first]
.getTranslation()
.norm();

View File

@ -362,8 +362,9 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x,
}
void ReducedModelOptimizer::createSimulationMeshes(
PatternGeometry &fullModel,
PatternGeometry &pattern,
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions,
std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)
{
@ -373,23 +374,25 @@ void ReducedModelOptimizer::createSimulationMeshes(
}
// Full pattern
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight});
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(pattern);
pFullPatternSimulationMesh->setBeamCrossSection(beamDimensions);
pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus);
pFullPatternSimulationMesh->reset();
// Reduced pattern
pReducedPatternSimulationMesh = std::make_shared<SimulationMesh>(reducedModel);
pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{beamWidth, beamHeight});
pReducedPatternSimulationMesh->setBeamCrossSection(beamDimensions);
pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus);
pReducedPatternSimulationMesh->reset();
}
void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel,
PatternGeometry &reducedModel)
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions)
{
ReducedModelOptimizer::createSimulationMeshes(fullModel,
reducedModel,
beamDimensions,
m_pFullPatternSimulationMesh,
m_pReducedModelSimulationMesh);
}
@ -558,9 +561,11 @@ ReducedModelOptimizer::ReducedModelOptimizer()
void ReducedModelOptimizer::initializePatterns(
PatternGeometry &fullPattern,
ReducedModel &reducedModel,
const std::array<xRange, NumberOfOptimizationVariables> &optimizationParameters)
const ReducedModelOptimization::Settings &optimizationSettings)
{
assert(fullPattern.VN() == reducedModel.VN() && fullPattern.EN() >= reducedModel.EN());
const std::array<xRange, NumberOfOptimizationVariables> &optimizationParameters
= optimizationSettings.variablesRanges;
#ifdef POLYSCOPE_DEFINED
polyscope::removeAllStructures();
@ -580,7 +585,9 @@ void ReducedModelOptimizer::initializePatterns(
optimizationState.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs;
g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs
= m_fullPatternOppositeInterfaceViPairs;
createSimulationMeshes(copyFullPattern, copyReducedModel);
createSimulationMeshes(copyFullPattern,
copyReducedModel,
optimizationSettings.patternBeamDimensions);
initializeUpdateReducedPatternFunctions();
initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
}
@ -596,7 +603,7 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern,
auto start = std::chrono::system_clock::now();
optimizationResults.label = fullPattern.getLabel();
optimizationResults.baseTriangleFullPattern.copy(fullPattern);
initializePatterns(fullPattern, reducedModel, optimizationSettings.variablesRanges);
initializePatterns(fullPattern, reducedModel, optimizationSettings);
optimize(optimizationSettings, optimizationResults);
optimizationResults.settings
= optimizationSettings; //NOTE: currently I set the max force base magn in optimize thus this must be called after the optimize function
@ -668,7 +675,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
const double beamHeight = beamWidth;
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
e.setDimensions(CrossSectionType(beamWidth, beamHeight));
}
};

View File

@ -90,8 +90,8 @@ private:
std::vector<bool> scenarioIsSymmetrical;
int fullPatternNumberOfEdges;
constexpr static double youngsModulus{1 * 1e9};
constexpr static double beamWidth{0.001};
constexpr static double beamHeight{0.001};
constexpr static double defaultBeamWidth{0.002};
constexpr static double defaultBeamHeight{0.002};
std::string fullPatternLabel;
// ReducedModelOptimization::Settings optimizationSettings;
@ -147,8 +147,9 @@ public:
const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);
static void createSimulationMeshes(
PatternGeometry &fullModel,
PatternGeometry &pattern,
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions,
std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
@ -288,10 +289,8 @@ private:
ReducedModelOptimization::S}));
void initializePatterns(PatternGeometry &fullPattern,
ReducedModel &reducedPattern,
const std::array<ReducedModelOptimization::xRange,
ReducedModelOptimization::NumberOfOptimizationVariables>
&optimizationParameters);
ReducedModel &reducedModel,
const ReducedModelOptimization::Settings &optimizationSettings);
static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
@ -299,7 +298,9 @@ private:
void runOptimization(const ReducedModelOptimization::Settings &settings,
ReducedModelOptimization::Results &results);
void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel);
void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
void createSimulationMeshes(PatternGeometry &fullModel,
PatternGeometry &reducedModel,
const RectangularBeamDimensions &beamDimensions);
void initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh,
const std::array<ReducedModelOptimization::xRange,

View File

@ -116,15 +116,11 @@ struct Settings
// 0.0127508,
// 1.18079,
// 0}; //final
std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
0.590241 / 6, 0.588372 / 6, 0.368304 / 2, 0.05, 1.18 / 4, 0}; //final b,h= 0.001
// std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
// 0.590241 / 6, 0.588372 / 6, 0.368304 / 2, 0.05, 1.18 / 4, 0}; //final b,h= 0.001
// std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0.590241 / 3,
// 0.588372 / 3,
// 0.368304,
// 0.1,
// 1.18 / 2,
// 0}; //final b,h= 0.002
std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{
0.590241 / 3, 0.588372 / 3, 0.368304, 0.1, 1.18 / 2, 0}; //final b,h= 0.002
// std::array<double, NumberOfBaseSimulationScenarios> baseScenarioMaxMagnitudes{0, 0, 0, 0.1, 0};
// std::array<double, NumberOfBaseSimulationScenarios>
@ -233,6 +229,7 @@ struct Settings
// double angularDistanceEpsilon{vcg::math::ToRad(2.0)};
double angularDistanceEpsilon{vcg::math::ToRad(0.0)};
double targetBaseTriangleSize{0.03};
RectangularBeamDimensions patternBeamDimensions{0.002, 0.002};
std::filesystem::path intermediateResultsDirectoryPath;
struct ObjectiveWeights
{
@ -310,6 +307,9 @@ struct Settings
json[GET_VARIABLE_NAME(angularDistanceEpsilon)] = vcg::math::ToDeg(angularDistanceEpsilon);
json[GET_VARIABLE_NAME(targetBaseTriangleSize)] = targetBaseTriangleSize;
json[GET_VARIABLE_NAME(baseScenarioMaxMagnitudes)] = baseScenarioMaxMagnitudes;
nlohmann::json json_dimensions;
patternBeamDimensions.to_json(json_dimensions, patternBeamDimensions);
json.update(json_dimensions);
#ifdef USE_ENSMALLEN
#ifdef USE_PSO
@ -412,6 +412,8 @@ struct Settings
pso.numberOfParticles = static_cast<int>(json[GET_VARIABLE_NAME(pso.numberOfParticles)]);
}
patternBeamDimensions.from_json(json, patternBeamDimensions);
// perBaseScenarioObjectiveWeights = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.translational = json.at(JsonKeys::ObjectiveWeights);
// objectiveWeights.rotational = 2 - objectiveWeights.translational;
@ -555,6 +557,7 @@ struct Settings
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
double fullPatternYoungsModulus{0};
//Full pattern
PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern
vcg::Triangle3<double> baseTriangle;
// std::string notes;
@ -808,9 +811,9 @@ struct Settings
JsonKeys::FullPatternLabel);
if (!baseTriangleFullPattern.load(
std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply").string())) {
!baseTriangleFullPattern.load(std::filesystem::path(loadFromPath)
.append(loadFromPath.stem().string() + ".ply")
.string());
baseTriangleFullPattern.load(std::filesystem::path(loadFromPath)
.append(loadFromPath.stem().string() + ".ply")
.string());
}
std::vector<double> baseTriangleVertices = json_optimizationResults.at(
@ -872,7 +875,7 @@ struct Settings
static void applyOptimizationResults_reducedModel_nonFanned(
const ReducedModelOptimization::Results &reducedPattern_optimizationResults,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
MeshType &reducedModel)
{
std::unordered_map<std::string, double>
optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
@ -884,13 +887,13 @@ struct Settings
applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("HexSize"),
optimalXVariables.at("HexAngle"),
patternBaseTriangle,
reducedPattern);
reducedModel);
return;
}
applyOptimizationResults_reducedModel_nonFanned(optimalXVariables.at("R"),
optimalXVariables.at("Theta"),
patternBaseTriangle,
reducedPattern);
reducedModel);
}
template<typename MeshType>
@ -898,7 +901,7 @@ struct Settings
const double &hexSize,
const double &hexAngle,
const vcg::Triangle3<double> &patternBaseTriangle,
MeshType &reducedPattern)
MeshType &reducedModel)
{
//Set optimal geometrical params of the reduced pattern
// CoordType center_barycentric(1, 0, 0);
@ -907,13 +910,13 @@ struct Settings
// (center_barycentric * (1 - hexSize) + interfaceEdgeMiddle_barycentric));
CoordType movableVertex_barycentric(1 - hexSize, hexSize / 2, hexSize / 2);
reducedPattern.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0]
+ patternBaseTriangle.cP(1) * movableVertex_barycentric[1]
+ patternBaseTriangle.cP(2) * movableVertex_barycentric[2];
reducedModel.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0]
+ patternBaseTriangle.cP(1) * movableVertex_barycentric[1]
+ patternBaseTriangle.cP(2) * movableVertex_barycentric[2];
if (hexAngle != 0) {
reducedPattern.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1},
vcg::math::ToRad(hexAngle))
* reducedPattern.vert[0].cP();
reducedModel.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1},
vcg::math::ToRad(hexAngle))
* reducedModel.vert[0].cP();
}
// for (int rotationCounter = 0;
// rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {