Refactoring. Changed reduced pattern to the single reduced for which in order to create the fan I connect to neighboring triangles. Changed Dome scenario.

This commit is contained in:
iasonmanolas 2021-03-26 11:58:13 +02:00
commit 4893f03667
3 changed files with 244 additions and 253 deletions

View File

@ -36,26 +36,25 @@ int main(int argc, char *argv[]) {
////Reduced pattern
const std::string filepath_reducedPattern = argv[2];
PatternGeometry reducedPattern(filepath_reducedPattern);
reducedPattern.setLabel(
std::filesystem::path(filepath_reducedPattern).stem().string());
reducedPattern.scale(0.03,interfaceNodeIndex);
reducedPattern.setLabel(std::filesystem::path(filepath_reducedPattern).stem().string());
reducedPattern.scale(0.03, interfaceNodeIndex);
// Set the optization settings
ReducedModelOptimization::xRange beamE{"E", 0.001, 10000};
ReducedModelOptimization::xRange beamA{"A", 0.001, 10000};
ReducedModelOptimization::xRange beamI2{"I2", 0.001,10000};
ReducedModelOptimization::xRange beamI3{"I3", 0.001,10000};
ReducedModelOptimization::xRange beamJ{"J", 0.001,10000};
ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.1, 0.8};
ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5};
ReducedModelOptimization::xRange beamE{"E", 0.001, 1000};
ReducedModelOptimization::xRange beamA{"A", 0.001, 1000};
ReducedModelOptimization::xRange beamI2{"I2", 0.001, 1000};
ReducedModelOptimization::xRange beamI3{"I3", 0.001, 1000};
ReducedModelOptimization::xRange beamJ{"J", 0.001, 1000};
ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
ReducedModelOptimization::Settings settings_optimization;
settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3,
innerHexagonSize, innerHexagonAngle};
const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls =
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.normalizationStrategy =
ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationStrategy
= ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationParameter = 0.0003;
settings_optimization.solutionAccuracy = 0.01;
@ -73,7 +72,7 @@ int main(int argc, char *argv[]) {
// Export results
const bool input_resultDirectoryDefined = argc >= 5;
std::string optimizationResultsDirectory =
input_resultDirectoryDefined ? argv[4] : "OptimizationResults";
input_resultDirectoryDefined ? argv[4] : std::filesystem::current_path().append("OptimizationResults");
std::string resultsOutputDir;
if (optimizationResults.numberOfSimulationCrashes != 0) {
const auto crashedJobsDirPath =

View File

@ -11,7 +11,8 @@ using namespace ReducedModelOptimization;
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<double> objectiveNormalizationValues;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
@ -31,19 +32,21 @@ struct GlobalOptimizationVariables {
int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{5};
ReducedModelOptimization::Settings optimizationSettings;
vcg::Triangle3<double> baseTriangle;
} global;
std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
double ReducedModelOptimizer::computeError(
double ReducedModelOptimizer::computeDisplacementError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor) {
const double rawError =
computeRawError(reducedPatternDisplacements, fullPatternDisplacements,
reducedToFullInterfaceViMap);
const double &normalizationFactor)
{
const double rawError = computeRawDisplacementError(reducedPatternDisplacements,
fullPatternDisplacements,
reducedToFullInterfaceViMap);
if (global.optimizationSettings.normalizationStrategy !=
Settings::NormalizationStrategy::NonNormalized) {
return rawError / normalizationFactor;
@ -51,11 +54,12 @@ double ReducedModelOptimizer::computeError(
return rawError;
}
double ReducedModelOptimizer::computeRawError(
double ReducedModelOptimizer::computeRawDisplacementError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) {
&reducedToFullInterfaceViMap)
{
double error = 0;
for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
const VertexIndex reducedModelVi = reducedFullViPair.first;
@ -82,12 +86,6 @@ double ReducedModelOptimizer::computeRawError(
return error;
}
double ReducedModelOptimizer::objective(double b, double r, double E) {
std::vector<double> x{b, r, E};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
double innerHexagonSize,
double innerHexagonRotationAngle) {
@ -95,19 +93,12 @@ double ReducedModelOptimizer::objective(double E,double A,double J,double I2,dou
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(double E,double A,
double innerHexagonSize,
double innerHexagonRotationAngle) {
std::vector<double> x{E,A, innerHexagonSize, innerHexagonRotationAngle};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(long n, const double *x) {
// std::cout.precision(17);
// for(int i=0;i<n;i++){
// std::cout<<x[i]<<" ";
// }
// std::cout<<std::endl;
// for (int i = 0; i < n; i++) {
// std::cout << x[i] << " ";
// }
// std::cout << std::endl;
// const Element &e =
// global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
@ -124,16 +115,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
// run simulations
double totalError = 0;
LinearSimulationModel simulator;
// FormFinder simulator;
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
//#ifdef POLYSCOPE_DEFINED
// global.reducedPatternSimulationJobs[simulationScenarioIndex]
// ->pMesh->registerForDrawing(colors.reducedInitial);
// reducedModelResults.registerForDrawing(colors.reducedDeformed);
// polyscope::show();
// reducedModelResults.unregister();
//#endif
// std::string filename;
if (!reducedModelResults.converged) {
totalError += std::numeric_limits<double>::max();
@ -142,12 +127,28 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
std::cout << "Failed simulation" << std::endl;
#endif
} else {
double simulationScenarioError = computeError(
reducedModelResults.displacements,
global.fullPatternDisplacements[simulationScenarioIndex],
global.reducedToFullInterfaceViMap,
global.objectiveNormalizationValues[simulationScenarioIndex]);
const bool usePotentialEnergy = false;
double simulationScenarioError;
if (usePotentialEnergy) {
simulationScenarioError = std::abs(
reducedModelResults.internalPotentialEnergy
- global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
} else {
simulationScenarioError = computeDisplacementError(
reducedModelResults.displacements,
global.fullPatternResults[simulationScenarioIndex].displacements,
global.reducedToFullInterfaceViMap,
global.objectiveNormalizationValues[simulationScenarioIndex]);
}
//#ifdef POLYSCOPE_DEFINED
// std::cout << "sim error:" << simulationScenarioError << std::endl;
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::reducedInitial);
// reducedModelResults.registerForDrawing(
// ReducedModelOptimization::Colors::reducedDeformed);
// polyscope::show();
// reducedModelResults.unregister();
//#endif
// if (global.optimizationSettings.normalizationStrategy !=
// NormalizationStrategy::Epsilon &&
// simulationScenarioError > 1) {
@ -186,10 +187,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
global.minX.assign(x, x + n);
}
#ifdef POLYSCOPE_DEFINED
if (++global.numberOfFunctionCalls % 100 == 0) {
std::cout << "Number of function calls:" << global.numberOfFunctionCalls
<< std::endl;
}
// if (++global.numberOfFunctionCalls % 100 == 0) {
// std::cout << "Number of function calls:" << global.numberOfFunctionCalls
// << std::endl;
// }
#endif
// compute error and return it
// global.objectiveValueHistory.push_back(totalError);
@ -256,9 +257,7 @@ void ReducedModelOptimizer::computeMaps(
const size_t baseTriangleInterfaceVi =
*(slotToNode.find(interfaceSlotIndex)->second.begin());
vcg::tri::Allocator<PatternGeometry>::PointerUpdater<
PatternGeometry::VertexPointer>
pu_fullModel;
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_fullModel;
fullPattern.deleteDanglingVertices(pu_fullModel);
const size_t fullModelBaseTriangleInterfaceVi =
pu_fullModel.remap.empty() ? baseTriangleInterfaceVi
@ -285,15 +284,13 @@ void ReducedModelOptimizer::computeMaps(
// Construct reduced->full and full->reduced interface vi map
reducedToFullInterfaceViMap.clear();
vcg::tri::Allocator<PatternGeometry>::PointerUpdater<
PatternGeometry::VertexPointer>
pu_reducedModel;
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_reducedModel;
reducedPattern.deleteDanglingVertices(pu_reducedModel);
const size_t reducedModelBaseTriangleInterfaceVi =
pu_reducedModel.remap[baseTriangleInterfaceVi];
const size_t reducedModelInterfaceVertexOffset =
reducedPattern.VN() - 1 /*- reducedModelBaseTriangleInterfaceVi*/;
reducedPattern.createFan();
const size_t reducedModelInterfaceVertexOffset
= reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/;
reducedPattern.createFan({1}); //TODO: should be an input parameter from main
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
reducedModelBaseTriangleInterfaceVi] =
@ -400,30 +397,31 @@ void ReducedModelOptimizer::initializePatterns(
PatternGeometry copyReducedPattern;
copyFullPattern.copy(fullPattern);
copyReducedPattern.copy(reducedPattern);
global.baseTriangle = copyReducedPattern.getBaseTriangle();
/*
* Here we create the vector that connects the central node with the bottom
* left node of the base triangle. During the optimization the vi%2==0 nodes
* move on these vectors.
* */
const double h = copyReducedPattern.getBaseTriangleHeight();
double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
0);
for (int rotationCounter = 0; rotationCounter < fanSize; rotationCounter++) {
VectorType rotatedVector =
vcg::RotationMatrix(patternPlaneNormal,
vcg::math::ToRad(rotationCounter * 60.0)) *
baseTriangle_leftBottomNode;
global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
}
const double innerHexagonInitialPos_x =
copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
// const double innerHexagonInitialPos_y =
// copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
global.innerHexagonInitialPos = innerHexagonInitialPos_x;
global.innerHexagonInitialRotationAngle =
30; /* NOTE: Here I assume that the CW reduced pattern is given as input.
This is not very generic */
// const double h = copyReducedPattern.getBaseTriangleHeight();
// double baseTriangle_bottomEdgeSize = 2 * h / tan(vcg::math::ToRad(60.0));
// VectorType baseTriangle_leftBottomNode(-baseTriangle_bottomEdgeSize / 2, -h,
// 0);
// for (int rotationCounter = 0; rotationCounter < fanSize; rotationCounter++) {
// VectorType rotatedVector =
// vcg::RotationMatrix(patternPlaneNormal,
// vcg::math::ToRad(rotationCounter * 60.0)) *
// baseTriangle_leftBottomNode;
// global.g_innerHexagonVectors[rotationCounter] = rotatedVector;
// }
// const double innerHexagonInitialPos_x =
// copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
// const double innerHexagonInitialPos_y =
// copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
// global.innerHexagonInitialPos = innerHexagonInitialPos_x;
// global.innerHexagonInitialRotationAngle =
// 30; /* NOTE: Here I assume that the CW reduced pattern is given as input.
// This is not very generic */
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
createSimulationMeshes(copyFullPattern, copyReducedPattern);
initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
@ -455,19 +453,38 @@ const double I3=global.initialParameters(4) * x[4];
}
assert(pReducedPatternSimulationMesh->EN() == 12);
assert(n>=2);
auto R = vcg::RotationMatrix(
ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x[n-1] - global.innerHexagonInitialRotationAngle));
CoordType center_barycentric(1, 0, 0);
CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
* x[n - 2]);
CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
* movableVertex_barycentric[0]
+ global.baseTriangle.cP(1)
* movableVertex_barycentric[1]
+ global.baseTriangle.cP(2)
* movableVertex_barycentric[2];
baseTriangleMovableVertexPosition
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, vcg::math::ToRad(x[n - 1]))
* baseTriangleMovableVertexPosition;
for (int rotationCounter = 0;
rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() =
R * global.g_innerHexagonVectors[rotationCounter] * x[n-2];
pReducedPatternSimulationMesh->vert[2 * rotationCounter + 1].P()
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* baseTriangleMovableVertexPosition;
}
pReducedPatternSimulationMesh->reset();
#ifdef POLYSCOPE_DEFINED
pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
#endif
//#ifdef POLYSCOPE_DEFINED
// pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
// pReducedPatternSimulationMesh->registerForDrawing();
// std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
// << std::endl;
// std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
// polyscope::show();
//#endif
}
@ -488,36 +505,35 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
void ReducedModelOptimizer::computeReducedModelSimulationJob(
const SimulationJob &simulationJobOfFullModel,
const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
SimulationJob &simulationJobOfReducedModel) {
const std::unordered_map<size_t, size_t> &fullToReducedMap,
SimulationJob &simulationJobOfReducedModel)
{
assert(simulationJobOfReducedModel.pMesh->VN() != 0);
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
reducedModelFixedVertices;
for (auto fullModelFixedVertex :
simulationJobOfFullModel.constrainedVertices) {
reducedModelFixedVertices[simulationJobFullToReducedMap.at(
fullModelFixedVertex.first)] = fullModelFixedVertex.second;
reducedModelFixedVertices[fullToReducedMap.at(fullModelFixedVertex.first)]
= fullModelFixedVertex.second;
}
std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
for (auto fullModelNodalForce :
simulationJobOfFullModel.nodalExternalForces) {
reducedModelNodalForces[simulationJobFullToReducedMap.at(
fullModelNodalForce.first)] = fullModelNodalForce.second;
reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)]
= fullModelNodalForce.second;
}
std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
for (auto fullForcedDisplacement : simulationJobOfFullModel.nodalForcedDisplacements) {
reducedNodalForcedDisplacements[fullToReducedMap.at(fullForcedDisplacement.first)]
= fullForcedDisplacement.second;
}
// std::unordered_map<VertexIndex, VectorType>
// reducedModelNodalForcedNormals; for (auto fullModelNodalForcedRotation :
// simulationJobOfFullModel.nodalForcedNormals) {
// reducedModelNodalForcedNormals[simulationJobFullToReducedMap.at(
// fullModelNodalForcedRotation.first)] =
// fullModelNodalForcedRotation.second;
// }
simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices;
simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces;
simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel();
// simulationJobOfReducedModel.nodalForcedNormals =
// reducedModelNodalForcedNormals;
simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements;
}
#if POLYSCOPE_DEFINED
@ -529,7 +545,7 @@ void ReducedModelOptimizer::visualizeResults(
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) {
FormFinder simulator;
DRMSimulationModel simulator;
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
fullPatternSimulationJobs[0]->pMesh;
pFullPatternSimulationMesh->registerForDrawing();
@ -559,9 +575,10 @@ void ReducedModelOptimizer::visualizeResults(
}
reducedModelResults.saveDeformedModel();
fullModelResults.saveDeformedModel();
double error = computeError(
reducedModelResults.displacements, fullModelResults.displacements,
reducedToFullInterfaceViMap, normalizationFactor);
double error = computeDisplacementError(reducedModelResults.displacements,
fullModelResults.displacements,
reducedToFullInterfaceViMap,
normalizationFactor);
std::cout << "Error of simulation scenario "
<< ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex] << " is "
<< error << std::endl;
@ -589,7 +606,7 @@ void ReducedModelOptimizer::registerResultsForDrawing(
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) {
FormFinder simulator;
DRMSimulationModel simulator;
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
pFullPatternSimulationJob->pMesh;
pFullPatternSimulationMesh->registerForDrawing();
@ -609,7 +626,7 @@ void ReducedModelOptimizer::registerResultsForDrawing(
simulator.executeSimulation(pReducedPatternSimulationJob);
// reducedModelResults.saveDeformedModel();
// fullModelResults.saveDeformedModel();
double error = computeRawError(
double error = computeRawDisplacementError(
reducedModelResults.displacements, fullModelResults.displacements,
reducedToFullInterfaceViMap/*,
global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
@ -692,12 +709,15 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
for (int simulationScenarioIndex:global.simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
const double error = computeError(
const double error = computeDisplacementError(
reducedModelResults.displacements,
global.fullPatternDisplacements[simulationScenarioIndex],
global.fullPatternResults[simulationScenarioIndex].displacements,
global.reducedToFullInterfaceViMap,
global.objectiveNormalizationValues[simulationScenarioIndex]);
results.rawObjectiveValue+=computeRawError(reducedModelResults.displacements,global.fullPatternDisplacements[simulationScenarioIndex],global.reducedToFullInterfaceViMap);
results.rawObjectiveValue += computeRawDisplacementError(
reducedModelResults.displacements,
global.fullPatternResults[simulationScenarioIndex].displacements,
global.reducedToFullInterfaceViMap);
results.objectiveValuePerSimulationScenario[simulationScenarioIndex] =
error;
}
@ -738,18 +758,6 @@ ReducedModelOptimizer::createScenarios(
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
// OldMethod
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType forceDirection =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// forceMagnitude * 10;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
scenarios[scenarioName] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
fixedVertices, nodalForces, {}));
@ -772,52 +780,16 @@ ReducedModelOptimizer::createScenarios(
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
// OldMethod
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType forceDirection = (v ^ patternPlaneNormal).Normalize();
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// 0.40 * forceMagnitude;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
scenarios[scenarioName] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
fixedVertices, nodalForces, {}));
// //// Torsion
// fixedVertices.clear();
// nodalForces.clear();
// for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
// viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
// const auto &viPair = *viPairIt;
// if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize();
// nodalForces[viPair.first] = Vector6d{
// 0, 0, 0, normalVDerivativeDir[0], normalVDerivativeDir[1], 0};
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// } else {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
// }
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces});
//// Bending
scenarioName = SimulationScenario::Bending;
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.5;
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.0005;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
@ -825,31 +797,48 @@ ReducedModelOptimizer::createScenarios(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
fixedVertices, nodalForces, {}));
//// Dome using moments
//// Dome
scenarioName = SimulationScenario::Dome;
fixedVertices.clear();
nodalForces.clear();
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
const auto viPair = *viPairIt;
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 2};
CoordType interfaceV = (pMesh->vert[viPair.first].cP()
- pMesh->vert[viPair.second].cP());
nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0],
-interfaceV[1],
0)
* 0.025;
nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceV[0],
interfaceV[1],
0)
* 0.025;
// CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// ^ CoordType(0, 0, -1).Normalize();
// nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
// * 0.0001;
// nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
// * 0.0001;
} else {
fixedVertices[viPair.first] = 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(0, 0, -1).Normalize();
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
* 0.0005;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
* 0.0005;
}
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
CoordType(0, 0, -1).Normalize();
nodalForces[viPair.first] =
Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 1;
nodalForces[viPair.second] =
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 1;
}
scenarios[scenarioName] = std::make_shared<SimulationJob>(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
fixedVertices, nodalForces, {}));
SimulationJob(pMesh,
simulationScenarioStrings[scenarioName],
fixedVertices,
nodalForces,
nodalForcedDisplacements));
//// Saddle
scenarioName = SimulationScenario::Saddle;
@ -862,18 +851,18 @@ ReducedModelOptimizer::createScenarios(
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
CoordType(0, 0, -1).Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
nodalForces[viPair.first] =
Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.2 * forceMagnitude;
nodalForces[viPair.second] =
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.2 * forceMagnitude;
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0002
* forceMagnitude;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0002
* forceMagnitude;
} else {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
nodalForces[viPair.first] =
Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.1 * forceMagnitude;
nodalForces[viPair.second] =
Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.1 * forceMagnitude;
nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0001
* forceMagnitude;
nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0001
* forceMagnitude;
}
}
scenarios[scenarioName] = std::make_shared<SimulationJob>(
@ -893,9 +882,9 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
double displacementNormSum = 0;
for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
const Vector6d &vertexDisplacement =
global.fullPatternDisplacements[simulationScenarioIndex]
[interfaceViPair.second];
const Vector6d &vertexDisplacement = global
.fullPatternResults[simulationScenarioIndex]
.displacements[interfaceViPair.second];
displacementNormSum += vertexDisplacement.getTranslation().norm();
}
fullPatternDisplacementNormSum[simulationScenarioIndex] =
@ -920,15 +909,16 @@ Results ReducedModelOptimizer::optimize(
global.simulationScenarioIndices = simulationScenarios;
if (global.simulationScenarioIndices.empty()) {
global.simulationScenarioIndices = {
SimulationScenario::Axial, SimulationScenario::Shear,
SimulationScenario::Bending, SimulationScenario::Dome,
SimulationScenario::Saddle};
global.simulationScenarioIndices = {SimulationScenario::Axial,
SimulationScenario::Shear,
SimulationScenario::Bending,
SimulationScenario::Dome,
SimulationScenario::Saddle};
}
global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios);
global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios);
global.fullPatternDisplacements.resize(NumberOfSimulationScenarios);
global.fullPatternResults.resize(NumberOfSimulationScenarios);
global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios);
global.minY = std::numeric_limits<double>::max();
global.numOfSimulationCrashes = 0;
@ -940,15 +930,19 @@ Results ReducedModelOptimizer::optimize(
NumberOfSimulationScenarios);
// polyscope::removeAllStructures();
FormFinder::Settings simulationSettings;
// settings.shouldDraw = true;
DRMSimulationModel::Settings simulationSettings;
simulationSettings.shouldDraw = false;
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::fullInitial);
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullModelResults = simulator.executeSimulation(
pFullPatternSimulationJob, simulationSettings);
global.fullPatternDisplacements[simulationScenarioIndex] =
fullModelResults.displacements;
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
// fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed);
// polyscope::show();
// fullPatternResults.unregister();
global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob,
@ -957,6 +951,7 @@ Results ReducedModelOptimizer::optimize(
global.reducedPatternSimulationJobs[simulationScenarioIndex] =
std::make_shared<SimulationJob>(reducedPatternSimulationJob);
}
// global.fullPatternSimulationJobs[0]->pMesh->unregister();
if (global.optimizationSettings.normalizationStrategy !=
Settings::NormalizationStrategy::NonNormalized) {

View File

@ -1,14 +1,14 @@
#ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP
#include "drmsimulationmodel.hpp"
#include "csvfile.hpp"
#include "drmsimulationmodel.hpp"
#include "edgemesh.hpp"
#include "simulationmesh.hpp"
#include "linearsimulationmodel.hpp"
#include "matplot/matplot.h"
#include <Eigen/Dense>
#include "reducedmodeloptimizer_structs.hpp"
#include "simulationmesh.hpp"
#include <Eigen/Dense>
#ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h"
@ -16,31 +16,32 @@
using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex;
class ReducedModelOptimizer {
std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
m_fullToReducedInterfaceViMap;
std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
m_fullPatternOppositeInterfaceViMap;
std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
class ReducedModelOptimizer
{
std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
m_fullToReducedInterfaceViMap;
std::unordered_map<FullPatternVertexIndex, FullPatternVertexIndex>
m_fullPatternOppositeInterfaceViMap;
std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
#ifdef POLYSCOPE_DEFINED
#endif // POLYSCOPE_DEFINED
public:
inline static int fanSize{6};
inline static VectorType patternPlaneNormal{0, 0, 1};
ReducedModelOptimization::Results optimize(const ReducedModelOptimization::Settings &xRanges,
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios =
std::vector<ReducedModelOptimization::SimulationScenario>());
ReducedModelOptimization::Results optimize(
const ReducedModelOptimization::Settings &xRanges,
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios
= std::vector<ReducedModelOptimization::SimulationScenario>());
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
static void computeReducedModelSimulationJob(
const SimulationJob &simulationJobOfFullModel,
const std::unordered_map<size_t, size_t> &simulationJobFullToReducedMap,
const std::unordered_map<size_t, size_t> &fullToReducedMap,
SimulationJob &simulationJobOfReducedModel);
SimulationJob
@ -75,15 +76,12 @@ public:
&fullToReducedInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
&fullPatternOppositeInterfaceViMap);
static void
visualizeResults(const std::vector<std::shared_ptr<SimulationJob>>
&fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>>
&reducedPatternSimulationJobs,
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static void visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static void registerResultsForDrawing(
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
@ -91,40 +89,39 @@ public:
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double
computeRawError(const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double computeRawDisplacementError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double
computeError(const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
static double computeDisplacementError(
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::vector<Vector6d> &fullPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor);
static double objective(double E, double A, double innerHexagonSize, double innerHexagonRotationAngle);
private:
static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static ReducedModelOptimization::Results runOptimization(const ReducedModelOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>>
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern,
const std::unordered_set<size_t> &reducedModelExcludedEges);
void createSimulationMeshes(PatternGeometry &fullModel,
PatternGeometry &reducedModel);
static void
initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh, const int &optimizationParamters);
static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static ReducedModelOptimization::Results runOptimization(
const ReducedModelOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>> createScenarios(
const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(PatternGeometry &fullModel,
PatternGeometry &reducedPattern,
const std::unordered_set<size_t> &reducedModelExcludedEges);
void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
static void initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh,
const int &optimizationParamters);
static double objective(long n, const double *x);
FormFinder simulator;
void computeObjectiveValueNormalizationFactors();
static double objective(long n, const double *x);
DRMSimulationModel simulator;
void computeObjectiveValueNormalizationFactors();
};
void updateMesh(long n, const double *x);
#endif // REDUCEDMODELOPTIMIZER_HPP