manolas
/
ReducedModelOptimization
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Replaced NumberOfSimulationScenarios with a variable calld totalNumberOfSimulationScenarios.This var will hold the number of base sim scenarios times the sub-sim scenarios for each case.

This commit is contained in:
iasonmanolas 2021-04-05 12:41:05 +03:00
parent d7d1951be6
commit 49494ccef8
4 changed files with 401 additions and 334 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists.txt
View File

@ -1,5 +1,8 @@
cmake_minimum_required(VERSION 2.8) cmake_minimum_required(VERSION 2.8)
project(ReducedModelOptimization) project(ReducedModelOptimization)
set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
message(STATUS "The compiler ${CMAKE_CXX_COMPILER}")
#Add the project cmake scripts to the module path #Add the project cmake scripts to the module path
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)

3
src/main.cpp
View File

@ -66,8 +66,7 @@ int main(int argc, char *argv[]) {
assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]); assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]);
ReducedModelOptimizer optimizer(numberOfNodesPerSlot); ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size()); optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size());
ReducedModelOptimization::Results optimizationResults = ReducedModelOptimization::Results optimizationResults = optimizer.optimize(settings_optimization);
optimizer.optimize(settings_optimization);
// Export results // Export results
const bool input_resultDirectoryDefined = argc >= 5; const bool input_resultDirectoryDefined = argc >= 5;

516
src/reducedmodeloptimizer.cpp
View File

@ -10,7 +10,6 @@
using namespace ReducedModelOptimization; using namespace ReducedModelOptimization;
struct GlobalOptimizationVariables { struct GlobalOptimizationVariables {
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
// std::vector<std::vector<Vector6d>> fullPatternDisplacements; // std::vector<std::vector<Vector6d>> fullPatternDisplacements;
std::vector<SimulationResults> fullPatternResults; std::vector<SimulationResults> fullPatternResults;
std::vector<double> translationalDisplacementNormalizationValues; std::vector<double> translationalDisplacementNormalizationValues;
@ -22,8 +21,7 @@ struct GlobalOptimizationVariables {
matplot::line_handle gPlotHandle; matplot::line_handle gPlotHandle;
std::vector<double> objectiveValueHistory; std::vector<double> objectiveValueHistory;
Eigen::VectorXd initialParameters; Eigen::VectorXd initialParameters;
std::vector<SimulationScenario> std::vector<int> simulationScenarioIndices;
simulationScenarioIndices;
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)}; std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
double innerHexagonInitialRotationAngle{30}; double innerHexagonInitialRotationAngle{30};
double innerHexagonInitialPos{0}; double innerHexagonInitialPos{0};
@ -36,8 +34,6 @@ struct GlobalOptimizationVariables {
vcg::Triangle3<double> baseTriangle; vcg::Triangle3<double> baseTriangle;
} global; } global;
std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
double ReducedModelOptimizer::computeDisplacementError( double ReducedModelOptimizer::computeDisplacementError(
const std::vector<Vector6d> &fullPatternDisplacements, const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements, const std::vector<Vector6d> &reducedPatternDisplacements,
@ -48,11 +44,7 @@ double ReducedModelOptimizer::computeDisplacementError(
const double rawError = computeRawDisplacementError(fullPatternDisplacements, const double rawError = computeRawDisplacementError(fullPatternDisplacements,
reducedPatternDisplacements, reducedPatternDisplacements,
reducedToFullInterfaceViMap); reducedToFullInterfaceViMap);
if (global.optimizationSettings.normalizationStrategy != return rawError / normalizationFactor;
Settings::NormalizationStrategy::NonNormalized) {
return rawError / normalizationFactor;
}
return rawError;
} }
double ReducedModelOptimizer::computeRawDisplacementError( double ReducedModelOptimizer::computeRawDisplacementError(
@ -75,7 +67,6 @@ double ReducedModelOptimizer::computeRawDisplacementError(
fullPatternDisplacements[fullModelVi][2]); fullPatternDisplacements[fullModelVi][2]);
Eigen::Vector3d errorVector = Eigen::Vector3d errorVector =
fullVertexDisplacement - reducedVertexDisplacement; fullVertexDisplacement - reducedVertexDisplacement;
// error += errorVector.squaredNorm();
error += errorVector.norm(); error += errorVector.norm();
} }
@ -181,7 +172,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
reducedModelResults.internalPotentialEnergy reducedModelResults.internalPotentialEnergy
- global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
} else { } else {
simulationScenarioError += computeError( simulationScenarioError = computeError(
global.fullPatternResults[simulationScenarioIndex], global.fullPatternResults[simulationScenarioIndex],
reducedModelResults, reducedModelResults,
global.reducedToFullInterfaceViMap, global.reducedToFullInterfaceViMap,
@ -287,14 +278,15 @@ void ReducedModelOptimizer::createSimulationMeshes(
void ReducedModelOptimizer::computeMaps( void ReducedModelOptimizer::computeMaps(
const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode, const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &fullPattern, PatternGeometry &reducedPattern, PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap, &reducedToFullInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
&fullToReducedInterfaceViMap, &fullToReducedInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
&fullPatternOppositeInterfaceViMap) { &fullPatternOppositeInterfaceViMap)
{
// Compute the offset between the interface nodes // Compute the offset between the interface nodes
const size_t interfaceSlotIndex = 4; // bottom edge const size_t interfaceSlotIndex = 4; // bottom edge
assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() &&
@ -349,12 +341,14 @@ void ReducedModelOptimizer::computeMaps(
// Create opposite vertex map // Create opposite vertex map
fullPatternOppositeInterfaceViMap.clear(); fullPatternOppositeInterfaceViMap.clear();
for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { fullPatternOppositeInterfaceViMap.reserve(fanSize / 2);
const size_t vi0 = fullModelBaseTriangleInterfaceVi + // for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
fanIndex * fullPatternInterfaceVertexOffset; for (int fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) {
const size_t vi0 = fullModelBaseTriangleInterfaceVi
+ fanIndex * fullPatternInterfaceVertexOffset;
const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
fullPatternOppositeInterfaceViMap[vi0] = vi1; fullPatternOppositeInterfaceViMap.emplace_back(std::make_pair(vi0, vi1));
} }
const bool debugMapping = false; const bool debugMapping = false;
@ -422,11 +416,11 @@ void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
reducedPattern, reducedPattern,
global.reducedToFullInterfaceViMap, global.reducedToFullInterfaceViMap,
m_fullToReducedInterfaceViMap, m_fullToReducedInterfaceViMap,
m_fullPatternOppositeInterfaceViMap); m_fullPatternOppositeInterfaceViPairs);
} }
ReducedModelOptimizer::ReducedModelOptimizer( ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot)
const std::vector<size_t> &numberOfNodesPerSlot) { {
FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot);
FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
} }
@ -561,103 +555,102 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements; simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements;
} }
#if POLYSCOPE_DEFINED //#if POLYSCOPE_DEFINED
void ReducedModelOptimizer::visualizeResults( //void ReducedModelOptimizer::visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> // const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
&fullPatternSimulationJobs, // const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> // const std::vector<BaseSimulationScenario> &simulationScenarios,
&reducedPatternSimulationJobs, // const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios, // &reducedToFullInterfaceViMap)
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> //{
&reducedToFullInterfaceViMap) { // DRMSimulationModel simulator;
DRMSimulationModel simulator; // std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh = // fullPatternSimulationJobs[0]->pMesh;
fullPatternSimulationJobs[0]->pMesh; // pFullPatternSimulationMesh->registerForDrawing();
pFullPatternSimulationMesh->registerForDrawing(); // pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply");
pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply"); // reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel()
reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel() // + "_undeformed.ply");
+ "_undeformed.ply"); // double totalError = 0;
double totalError = 0; // for (const int simulationScenarioIndex : simulationScenarios) {
for (const int simulationScenarioIndex : simulationScenarios) { // const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob = // fullPatternSimulationJobs[simulationScenarioIndex];
fullPatternSimulationJobs[simulationScenarioIndex]; // pFullPatternSimulationJob->registerForDrawing(
pFullPatternSimulationJob->registerForDrawing( // pFullPatternSimulationMesh->getLabel());
pFullPatternSimulationMesh->getLabel()); // SimulationResults fullModelResults =
SimulationResults fullModelResults = // simulator.executeSimulation(pFullPatternSimulationJob);
simulator.executeSimulation(pFullPatternSimulationJob); // fullModelResults.registerForDrawing();
fullModelResults.registerForDrawing(); // // fullModelResults.saveDeformedModel();
// fullModelResults.saveDeformedModel(); // const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob = // reducedPatternSimulationJobs[simulationScenarioIndex];
reducedPatternSimulationJobs[simulationScenarioIndex]; // SimulationResults reducedModelResults =
SimulationResults reducedModelResults = // simulator.executeSimulation(pReducedPatternSimulationJob);
simulator.executeSimulation(pReducedPatternSimulationJob); // double normalizationFactor = 1;
double normalizationFactor = 1; // if (global.optimizationSettings.normalizationStrategy !=
if (global.optimizationSettings.normalizationStrategy != // ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { // normalizationFactor
normalizationFactor // = global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
= global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; // }
} // reducedModelResults.saveDeformedModel();
reducedModelResults.saveDeformedModel(); // fullModelResults.saveDeformedModel();
fullModelResults.saveDeformedModel(); // double error = computeDisplacementError(reducedModelResults.displacements,
double error = computeDisplacementError(reducedModelResults.displacements, // fullModelResults.displacements,
fullModelResults.displacements, // reducedToFullInterfaceViMap,
reducedToFullInterfaceViMap, // normalizationFactor);
normalizationFactor); // std::cout << "Error of simulation scenario "
std::cout << "Error of simulation scenario " // <baseSimulationScenarioNames[simulationScenarioIndex]
<< ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex] << " is " // << " is " << error << std::endl;
<< error << std::endl; // totalError += error;
totalError += error; // reducedModelResults.registerForDrawing();
reducedModelResults.registerForDrawing(); // // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
// firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); // // registerWorldAxes();
// registerWorldAxes(); // const std::string screenshotFilename
const std::string screenshotFilename = // = "/home/iason/Coding/Projects/Approximating shapes with flat "
"/home/iason/Coding/Projects/Approximating shapes with flat " // "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
"patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" // "Images/"
"Images/" + // + pFullPatternSimulationMesh->getLabel() + "_"
pFullPatternSimulationMesh->getLabel() + "_" + // + baseSimulationScenarioNames[simulationScenarioIndex];
ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex]; // polyscope::show();
polyscope::show(); // polyscope::screenshot(screenshotFilename, false);
polyscope::screenshot(screenshotFilename, false); // fullModelResults.unregister();
fullModelResults.unregister(); // reducedModelResults.unregister();
reducedModelResults.unregister(); // // firstOptimizationRoundResults[simulationScenarioIndex].unregister();
// firstOptimizationRoundResults[simulationScenarioIndex].unregister(); // }
} // std::cout << "Total error:" << totalError << std::endl;
std::cout << "Total error:" << totalError << std::endl; //}
}
void ReducedModelOptimizer::registerResultsForDrawing( //void ReducedModelOptimizer::registerResultsForDrawing(
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob, // const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob, // const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> // const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) { // &reducedToFullInterfaceViMap) {
DRMSimulationModel simulator; // DRMSimulationModel simulator;
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh = // std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
pFullPatternSimulationJob->pMesh; // pFullPatternSimulationJob->pMesh;
pFullPatternSimulationMesh->registerForDrawing(); // pFullPatternSimulationMesh->registerForDrawing();
// pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() // // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel()
// + // // +
// "_undeformed.ply"); // // "_undeformed.ply");
// reducedPatternSimulationJobs[0]->pMesh->savePly( // // reducedPatternSimulationJobs[0]->pMesh->savePly(
// reducedPatternSimulationJobs[0]->pMesh->getLabel() + // // reducedPatternSimulationJobs[0]->pMesh->getLabel() +
// "_undeformed.ply"); // // "_undeformed.ply");
pFullPatternSimulationJob->registerForDrawing( // pFullPatternSimulationJob->registerForDrawing(
pFullPatternSimulationMesh->getLabel()); // pFullPatternSimulationMesh->getLabel());
SimulationResults fullModelResults = // SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob); // simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing(); // fullModelResults.registerForDrawing();
// fullModelResults.saveDeformedModel(); // // fullModelResults.saveDeformedModel();
SimulationResults reducedModelResults = // SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob); // simulator.executeSimulation(pReducedPatternSimulationJob);
// reducedModelResults.saveDeformedModel(); // // reducedModelResults.saveDeformedModel();
// fullModelResults.saveDeformedModel(); // // fullModelResults.saveDeformedModel();
double error = computeRawDisplacementError( // double error = computeRawDisplacementError(
reducedModelResults.displacements, fullModelResults.displacements, // reducedModelResults.displacements, fullModelResults.displacements,
reducedToFullInterfaceViMap/*, // reducedToFullInterfaceViMap/*,
global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); // global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
std::cout << "Error is " << error << std::endl; // std::cout << "Error is " << error << std::endl;
reducedModelResults.registerForDrawing(); // reducedModelResults.registerForDrawing();
} //}
#endif // POLYSCOPE_DEFINED //#endif // POLYSCOPE_DEFINED
void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults, const SimulationResults &fullModelResults,
@ -709,31 +702,31 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
results.optimalXNameValuePairs[xVariableIndex] results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex]); = std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex]);
} }
results.objectiveValue = global.minY;
#ifdef POLYSCOPE_DEFINED
std::cout<<"Total optimal objective value:"<<global.minY<<std::endl;
if (global.minY != result.y) {
std::cerr << "Global min objective is not equal to result objective"
<< std::endl;
}
#endif
// Compute obj value per simulation scenario and the raw objective value // Compute obj value per simulation scenario and the raw objective value
results.rawObjectiveValue=0; results.rawObjectiveValue=0;
std::vector<double> optimalX(results.optimalXNameValuePairs.size()); std::vector<double> optimalX(results.optimalXNameValuePairs.size());
for(int xVariableIndex=0;xVariableIndex<global.numberOfOptimizationParameters;xVariableIndex++){ for (int xVariableIndex = 0; xVariableIndex < global.numberOfOptimizationParameters;
optimalX[xVariableIndex]= xVariableIndex++) {
global.minX[xVariableIndex]; optimalX[xVariableIndex] = global.minX[xVariableIndex];
assert(global.minX[xVariableIndex] == result.x(xVariableIndex));
} }
updateMesh(optimalX.size(), optimalX.data()); updateMesh(optimalX.size(), optimalX.data());
results.objectiveValuePerSimulationScenario.resize( results.objectiveValuePerSimulationScenario.resize(totalNumberOfSimulationScenarios);
ReducedModelOptimization::NumberOfSimulationScenarios);
LinearSimulationModel simulator; LinearSimulationModel simulator;
for (int simulationScenarioIndex:global.simulationScenarioIndices) { double totalObjectiveValue = 0;
double totalRawObjectiveValue = 0;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation( SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex]); global.reducedPatternSimulationJobs[simulationScenarioIndex]);
totalObjectiveValue += computeError(
global.fullPatternResults[simulationScenarioIndex],
reducedModelResults,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
const double rawTranslationalError = computeRawDisplacementError( const double rawTranslationalError = computeRawDisplacementError(
global.fullPatternResults[simulationScenarioIndex].displacements, global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements, reducedModelResults.displacements,
@ -748,17 +741,55 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
const double normalizedTranslationalError const double normalizedTranslationalError
= rawTranslationalError = rawTranslationalError
/ global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; / global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
const double test_normalizedTranslationError = computeDisplacementError(
global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
const double normalizedRotationalError const double normalizedRotationalError
= rawRotationalError = rawRotationalError
/ global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]; / global.rotationalDisplacementNormalizationValues[simulationScenarioIndex];
std::cout << "Simulation scenario:" << simulationScenarioStrings[simulationScenarioIndex] const double test_normalizedRotationalError = computeRotationalError(
global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
reducedModelResults.rotationalDisplacementQuaternion,
global.reducedToFullInterfaceViMap,
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
assert(test_normalizedTranslationError == normalizedTranslationalError);
assert(test_normalizedRotationalError == normalizedRotationalError);
std::cout << "Simulation scenario:" << baseSimulationScenarioNames[simulationScenarioIndex]
<< std::endl;
std::cout << "raw translational error:" << rawTranslationalError << std::endl;
std::cout << "translation normalization value:"
<< global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
<< std::endl; << std::endl;
std::cout << "Translational error:" << normalizedTranslationalError << std::endl; std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
std::cout << "raw Rotational error:" << rawRotationalError << std::endl;
std::cout << "rotational normalization value:"
<< global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
<< std::endl;
std::cout << "Rotational error:" << normalizedRotationalError << std::endl; std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
std::cout << std::endl;
results.objectiveValuePerSimulationScenario[simulationScenarioIndex] results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
= normalizedTranslationalError + normalizedRotationalError; = normalizedTranslationalError + normalizedRotationalError;
std::cout << "Objective value:"
<< results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
<< std::endl;
// totalObjectiveValue += results.objectiveValuePerSimulationScenario[simulationScenarioIndex];
totalRawObjectiveValue += rawTranslationalError + rawRotationalError;
std::cout << std::endl;
} }
assert(result.y == global.minY);
results.objectiveValue = result.y;
auto functionReturnValue = objective(optimalX.size(), optimalX.data());
std::cout << "Total (functionReturn)" << functionReturnValue << std::endl;
std::cout << "Total:" << totalObjectiveValue << std::endl;
#ifdef POLYSCOPE_DEFINED
std::cout << "Total optimal objective value:" << results.objectiveValue << std::endl;
// std::cout << "Total raw objective value:" << global.minY << std::endl;
if (global.minY != result.y) {
std::cerr << "Global min objective is not equal to result objective" << std::endl;
}
#endif
results.time = elapsed.count() / 1000.0; results.time = elapsed.count() / 1000.0;
const bool printDebugInfo = false; const bool printDebugInfo = false;
if (printDebugInfo) { if (printDebugInfo) {
@ -775,75 +806,82 @@ std::vector<std::shared_ptr<SimulationJob>>
ReducedModelOptimizer::createScenarios( ReducedModelOptimizer::createScenarios(
const std::shared_ptr<SimulationMesh> &pMesh) { const std::shared_ptr<SimulationMesh> &pMesh) {
std::vector<std::shared_ptr<SimulationJob>> scenarios; std::vector<std::shared_ptr<SimulationJob>> scenarios;
scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); scenarios.resize(totalNumberOfSimulationScenarios);
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices; std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalForces; std::unordered_map<VertexIndex, Vector6d> nodalForces;
const double forceMagnitude = 10;
//// Axial //// Axial
SimulationScenario scenarioName = SimulationScenario::Axial; int simulationScenarioIndex = BaseSimulationScenario::Axial;
// NewMethod for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(1, 0, 0); CoordType forceDirection(1, 0, 0);
const auto viPair = *viPairIt; const auto viPair = *viPairIt;
nodalForces[viPair.first] = nodalForces[viPair.first]
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
0, 0}) * * 500;
forceMagnitude * 10;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
} }
scenarios[scenarioName] = std::make_shared<SimulationJob>( scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName], SimulationJob(pMesh,
fixedVertices, nodalForces, {})); "Axial_" + std::to_string(simulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
//// Shear //// Shear
scenarioName = SimulationScenario::Shear; simulationScenarioIndex = BaseSimulationScenario::Shear;
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
// NewMethod // NewMethod
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(0, 1, 0); CoordType forceDirection(0, 1, 0);
const auto viPair = *viPairIt; const auto viPair = *viPairIt;
nodalForces[viPair.first] = nodalForces[viPair.first]
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * 40;
0, 0}) *
forceMagnitude * 4;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
} }
scenarios[scenarioName] = std::make_shared<SimulationJob>( scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName], SimulationJob(pMesh,
fixedVertices, nodalForces, {})); baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
//// Bending //// Bending
scenarioName = SimulationScenario::Bending; simulationScenarioIndex = BaseSimulationScenario::Bending;
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) {
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.0005; nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * 0.005;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
scenarios[scenarioName] = std::make_shared<SimulationJob>( scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName], SimulationJob(pMesh,
fixedVertices, nodalForces, {})); baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
//// Dome //// Dome
scenarioName = SimulationScenario::Dome; simulationScenarioIndex = BaseSimulationScenario::Dome;
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements; std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto viPair = *viPairIt; const auto viPair = *viPairIt;
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType interfaceV = (pMesh->vert[viPair.first].cP() CoordType interfaceV = (pMesh->vert[viPair.first].cP()
- pMesh->vert[viPair.second].cP()); - pMesh->vert[viPair.second].cP());
nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0], nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0],
@ -865,47 +903,45 @@ ReducedModelOptimizer::createScenarios(
fixedVertices[viPair.second] = std::unordered_set<DoFType>{2}; fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
^ CoordType(0, 0, -1).Normalize(); ^ CoordType(0, 0, -1).Normalize();
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.005;
* 0.0005; nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.005;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
* 0.0005;
} }
} }
scenarios[scenarioName] = std::make_shared<SimulationJob>( scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, SimulationJob(pMesh,
simulationScenarioStrings[scenarioName], baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices, fixedVertices,
nodalForces, nodalForces,
nodalForcedDisplacements)); nodalForcedDisplacements));
//// Saddle //// Saddle
scenarioName = SimulationScenario::Saddle; simulationScenarioIndex = BaseSimulationScenario::Saddle;
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto &viPair = *viPairIt; const auto &viPair = *viPairIt;
CoordType v = CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^ (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
CoordType(0, 0, -1).Normalize(); CoordType(0, 0, -1).Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0002 nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.002;
* forceMagnitude; nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.002;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0002
* forceMagnitude;
} else { } else {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2}; fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2}; fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0001 nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.001;
* forceMagnitude; nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.001;
nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0001
* forceMagnitude;
} }
} }
scenarios[scenarioName] = std::make_shared<SimulationJob>( scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName], SimulationJob(pMesh,
fixedVertices, nodalForces, {})); baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
return scenarios; return scenarios;
} }
@ -913,11 +949,11 @@ ReducedModelOptimizer::createScenarios(
void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
// Compute the sum of the displacement norms // Compute the sum of the displacement norms
std::vector<double> fullPatternTranslationalDisplacementNormSum(NumberOfSimulationScenarios); std::vector<double> fullPatternTranslationalDisplacementNormSum(
std::vector<double> fullPatternAngularDistance(NumberOfSimulationScenarios); totalNumberOfSimulationScenarios);
std::vector<double> fullPatternAngularDistance(totalNumberOfSimulationScenarios);
for (int simulationScenarioIndex : global.simulationScenarioIndices) { for (int simulationScenarioIndex : global.simulationScenarioIndices) {
double translationalDisplacementNormSum = 0; double translationalDisplacementNormSum = 0;
double angularDistanceSum = 0;
for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
const int fullPatternVi = interfaceViPair.second; const int fullPatternVi = interfaceViPair.second;
//If the full pattern vertex is translationally constrained dont take it into account //If the full pattern vertex is translationally constrained dont take it into account
@ -936,6 +972,10 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
.fullPatternResults[simulationScenarioIndex] .fullPatternResults[simulationScenarioIndex]
.displacements[fullPatternVi]; .displacements[fullPatternVi];
translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm();
}
double angularDistanceSum = 0;
for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
const int fullPatternVi = interfaceViPair.second;
//If the full pattern vertex is rotationally constrained dont take it into account //If the full pattern vertex is rotationally constrained dont take it into account
if (global.fullPatternSimulationJobs[simulationScenarioIndex] if (global.fullPatternSimulationJobs[simulationScenarioIndex]
->constrainedVertices.contains(fullPatternVi)) { ->constrainedVertices.contains(fullPatternVi)) {
@ -951,6 +991,7 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
.rotationalDisplacementQuaternion[fullPatternVi] .rotationalDisplacementQuaternion[fullPatternVi]
.angularDistance(Eigen::Quaterniond::Identity()); .angularDistance(Eigen::Quaterniond::Identity());
} }
fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex]
= translationalDisplacementNormSum; = translationalDisplacementNormSum;
fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum;
@ -964,72 +1005,91 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex] global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
epsilon_translationalDisplacement); epsilon_translationalDisplacement);
const double epsilon_rotationalDisplacement = vcg::math::ToRad(2.0); // const double epsilon_rotationalDisplacement = vcg::math::ToRad(10.0);
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternAngularDistance[simulationScenarioIndex], = /*std::max(*/ fullPatternAngularDistance[simulationScenarioIndex] /*,
epsilon_rotationalDisplacement); epsilon_rotationalDisplacement)*/
;
} else { } else {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
} }
} }
} }
Results ReducedModelOptimizer::optimize( Results ReducedModelOptimizer::optimize(
const Settings &optimizationSettings, const Settings &optimizationSettings,
const std::vector<SimulationScenario> &simulationScenarios) { const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
{
global.simulationScenarioIndices = simulationScenarios; for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) {
if (global.simulationScenarioIndices.empty()) { //Increase the size of the vector holding the simulation scenario indices
global.simulationScenarioIndices = {SimulationScenario::Axial, global.simulationScenarioIndices.resize(
SimulationScenario::Shear, global.simulationScenarioIndices.size()
SimulationScenario::Bending, + simulationScenariosResolution[baseSimulationScenarioIndex]);
SimulationScenario::Dome, //Add the simulation scenarios indices that correspond to this base simulation scenario
SimulationScenario::Saddle}; std::iota(global.simulationScenarioIndices.end()
- simulationScenariosResolution[baseSimulationScenarioIndex],
global.simulationScenarioIndices.end(),
std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin()
+ baseSimulationScenarioIndex,
0));
}
if (desiredBaseSimulationScenarioIndices.empty()) {
global.simulationScenarioIndices.resize(totalNumberOfSimulationScenarios);
std::iota(global.simulationScenarioIndices.begin(),
global.simulationScenarioIndices.end(),
0);
} }
global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios);
global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); global.fullPatternResults.resize(totalNumberOfSimulationScenarios);
global.fullPatternResults.resize(NumberOfSimulationScenarios); global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
global.translationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios); global.rotationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
global.rotationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios);
global.minY = std::numeric_limits<double>::max(); global.minY = std::numeric_limits<double>::max();
global.numOfSimulationCrashes = 0; global.numOfSimulationCrashes = 0;
global.numberOfFunctionCalls = 0; global.numberOfFunctionCalls = 0;
global.optimizationSettings = optimizationSettings; global.optimizationSettings = optimizationSettings;
global.fullPatternSimulationJobs = global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh);
createScenarios(m_pFullPatternSimulationMesh);
reducedPatternMaximumDisplacementSimulationJobs.resize(
NumberOfSimulationScenarios);
// polyscope::removeAllStructures(); // polyscope::removeAllStructures();
DRMSimulationModel::Settings simulationSettings; DRMSimulationModel::Settings simulationSettings;
simulationSettings.shouldDraw = false; simulationSettings.shouldDraw = false;
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::fullInitial); // ReducedModelOptimization::Colors::fullInitial);
// LinearSimulationModel linearSimulator;
for (int simulationScenarioIndex : global.simulationScenarioIndices) { for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob = const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
global.fullPatternSimulationJobs[simulationScenarioIndex]; = global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings); simulationSettings);
// fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed); // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
// fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
// polyscope::show(); // polyscope::show();
// fullPatternResults.unregister(); // fullPatternResults.unregister();
// fullPatternResults_linear.unregister();
global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
SimulationJob reducedPatternSimulationJob; SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob, computeReducedModelSimulationJob(*pFullPatternSimulationJob,
m_fullToReducedInterfaceViMap, m_fullToReducedInterfaceViMap,
reducedPatternSimulationJob); reducedPatternSimulationJob);
global.reducedPatternSimulationJobs[simulationScenarioIndex] = global.reducedPatternSimulationJobs[simulationScenarioIndex]
std::make_shared<SimulationJob>(reducedPatternSimulationJob); = std::make_shared<SimulationJob>(reducedPatternSimulationJob);
} }
// global.fullPatternSimulationJobs[0]->pMesh->unregister(); // global.fullPatternSimulationJobs[0]->pMesh->unregister();
if (global.optimizationSettings.normalizationStrategy != // if (global.optimizationSettings.normalizationStrategy
Settings::NormalizationStrategy::NonNormalized) { // != Settings::NormalizationStrategy::NonNormalized) {
computeObjectiveValueNormalizationFactors(); computeObjectiveValueNormalizationFactors();
} // }
Results optResults = runOptimization(optimizationSettings); Results optResults = runOptimization(optimizationSettings);
for (int simulationScenarioIndex : global.simulationScenarioIndices) { for (int simulationScenarioIndex : global.simulationScenarioIndices) {
optResults.fullPatternSimulationJobs.push_back( optResults.fullPatternSimulationJobs.push_back(

213
src/reducedmodeloptimizer.hpp
View File

@ -22,127 +22,132 @@ class ReducedModelOptimizer
std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh; std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
m_fullToReducedInterfaceViMap; m_fullToReducedInterfaceViMap;
std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex> std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
m_fullPatternOppositeInterfaceViMap; m_fullPatternOppositeInterfaceViPairs;
std::unordered_map<size_t, size_t> nodeToSlot; std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode; std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
#ifdef POLYSCOPE_DEFINED
#endif // POLYSCOPE_DEFINED
public: public:
inline static int fanSize{6}; constexpr static std::array<int, 5> simulationScenariosResolution = {1, 1, 1, 1, 1};
inline static VectorType patternPlaneNormal{0, 0, 1}; inline static int totalNumberOfSimulationScenarios
ReducedModelOptimization::Results optimize( = std::accumulate(simulationScenariosResolution.begin(),
const ReducedModelOptimization::Settings &xRanges, simulationScenariosResolution.end(),
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios 0);
= std::vector<ReducedModelOptimization::SimulationScenario>()); inline static int fanSize{6};
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; inline static VectorType patternPlaneNormal{0, 0, 1};
ReducedModelOptimization::Results optimize(
const ReducedModelOptimization::Settings &xRanges,
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios
= std::vector<ReducedModelOptimization::BaseSimulationScenario>());
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot); ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
static void computeReducedModelSimulationJob( static void computeReducedModelSimulationJob(
const SimulationJob &simulationJobOfFullModel, const SimulationJob &simulationJobOfFullModel,
const std::unordered_map<size_t, size_t> &fullToReducedMap, const std::unordered_map<size_t, size_t> &fullToReducedMap,
SimulationJob &simulationJobOfReducedModel); SimulationJob &simulationJobOfReducedModel);
SimulationJob SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
void initializePatterns(PatternGeometry &fullPattern, void initializePatterns(PatternGeometry &fullPattern,
PatternGeometry &reducedPatterm, PatternGeometry &reducedPatterm,
const int &optimizationParameters); const int &optimizationParameters);
static void runSimulation(const std::string &filename, static void runSimulation(const std::string &filename, std::vector<double> &x);
std::vector<double> &x);
static double objective(double x2, double A, double J, double I2, double I3, double x3, static double objective(double x2,
double innerHexagonRotationAngle); double A,
static double objective(double b, double r, double E); double J,
double I2,
double I3,
double x3,
double innerHexagonRotationAngle);
static double objective(double b, double r, double E);
static std::vector<std::shared_ptr<SimulationJob>> static std::vector<std::shared_ptr<SimulationJob>> createScenarios(
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh, const std::shared_ptr<SimulationMesh> &pMesh,
const std::unordered_map<size_t, size_t> const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);
&fullPatternOppositeInterfaceViMap);
static void createSimulationMeshes( static void createSimulationMeshes(
PatternGeometry &fullModel, PatternGeometry &reducedModel, PatternGeometry &fullModel,
std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh, PatternGeometry &reducedModel,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh); std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
static void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
PatternGeometry &fullPattern, static void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &reducedPattern, PatternGeometry &fullPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> PatternGeometry &reducedPattern,
&reducedToFullInterfaceViMap, std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> &reducedToFullInterfaceViMap,
&fullToReducedInterfaceViMap, std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> &fullToReducedInterfaceViMap,
&fullPatternOppositeInterfaceViMap); std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
static void visualizeResults( &fullPatternOppositeInterfaceViMap);
const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs, static void visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs, const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios, const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios,
&reducedToFullInterfaceViMap); const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
static void registerResultsForDrawing( &reducedToFullInterfaceViMap);
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob, static void registerResultsForDrawing(
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob, const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex, const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap); &reducedToFullInterfaceViMap);
static double computeRawDisplacementError( static double computeRawDisplacementError(
const std::vector<Vector6d> &fullPatternDisplacements, const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements, const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap); &reducedToFullInterfaceViMap);
static double computeDisplacementError( static double computeDisplacementError(
const std::vector<Vector6d> &fullPatternDisplacements, const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements, const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap, &reducedToFullInterfaceViMap,
const double &normalizationFactor); const double &normalizationFactor);
static double objective(double E, static double objective(double E,
double A, double A,
double innerHexagonSize, double innerHexagonSize,
double innerHexagonRotationAngle); double innerHexagonRotationAngle);
static double computeRawRotationalError( static double computeRawRotationalError(
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern, const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern, const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap); &reducedToFullInterfaceViMap);
static double computeRotationalError( static double computeRotationalError(
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern, const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern, const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap, &reducedToFullInterfaceViMap,
const double &normalizationFactor); const double &normalizationFactor);
static double computeError( static double computeError(
const SimulationResults &simulationResults_fullPattern, const SimulationResults &simulationResults_fullPattern,
const SimulationResults &simulationResults_reducedPattern, const SimulationResults &simulationResults_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap, &reducedToFullInterfaceViMap,
const double &normalizationFactor_translationalDisplacement, const double &normalizationFactor_translationalDisplacement,
const double &normalizationFactor_rotationalDisplacement); const double &normalizationFactor_rotationalDisplacement);
private: private:
static void computeDesiredReducedModelDisplacements( static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults, const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap, const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static ReducedModelOptimization::Results runOptimization( static ReducedModelOptimization::Results runOptimization(
const ReducedModelOptimization::Settings &settings); const ReducedModelOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>> createScenarios( std::vector<std::shared_ptr<SimulationJob>> createScenarios(
const std::shared_ptr<SimulationMesh> &pMesh); const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh, static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
const int &optimizationParamters); const int &optimizationParamters);
static double objective(long n, const double *x); static double objective(long n, const double *x);
DRMSimulationModel simulator; DRMSimulationModel simulator;
void computeObjectiveValueNormalizationFactors(); void computeObjectiveValueNormalizationFactors();
}; };
void updateMesh(long n, const double *x); void updateMesh(long n, const double *x);
#endif // REDUCEDMODELOPTIMIZER_HPP #endif // REDUCEDMODELOPTIMIZER_HPP