diff --git a/src/main.cpp b/src/main.cpp index 3179038..fe4cdbb 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -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 = diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index e0cfe5d..fb0f687 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -11,7 +11,8 @@ using namespace ReducedModelOptimization; struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; - std::vector> fullPatternDisplacements; + // std::vector> fullPatternDisplacements; + std::vector fullPatternResults; std::vector objectiveNormalizationValues; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; @@ -31,19 +32,21 @@ struct GlobalOptimizationVariables { int numberOfFunctionCalls{0}; int numberOfOptimizationParameters{5}; ReducedModelOptimization::Settings optimizationSettings; + vcg::Triangle3 baseTriangle; } global; std::vector reducedPatternMaximumDisplacementSimulationJobs; -double ReducedModelOptimizer::computeError( +double ReducedModelOptimizer::computeDisplacementError( const std::vector &reducedPatternDisplacements, const std::vector &fullPatternDisplacements, const std::unordered_map &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 &reducedPatternDisplacements, const std::vector &fullPatternDisplacements, const std::unordered_map - &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 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 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;ipMesh->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::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::PointerUpdater< - PatternGeometry::VertexPointer> - pu_fullModel; + vcg::tri::Allocator::PointerUpdater 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::PointerUpdater< - PatternGeometry::VertexPointer> - pu_reducedModel; + vcg::tri::Allocator::PointerUpdater 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 &simulationJobFullToReducedMap, - SimulationJob &simulationJobOfReducedModel) { + const std::unordered_map &fullToReducedMap, + SimulationJob &simulationJobOfReducedModel) +{ assert(simulationJobOfReducedModel.pMesh->VN() != 0); std::unordered_map> reducedModelFixedVertices; for (auto fullModelFixedVertex : simulationJobOfFullModel.constrainedVertices) { - reducedModelFixedVertices[simulationJobFullToReducedMap.at( - fullModelFixedVertex.first)] = fullModelFixedVertex.second; + reducedModelFixedVertices[fullToReducedMap.at(fullModelFixedVertex.first)] + = fullModelFixedVertex.second; } std::unordered_map reducedModelNodalForces; for (auto fullModelNodalForce : simulationJobOfFullModel.nodalExternalForces) { - reducedModelNodalForces[simulationJobFullToReducedMap.at( - fullModelNodalForce.first)] = fullModelNodalForce.second; + reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)] + = fullModelNodalForce.second; + } + + std::unordered_map reducedNodalForcedDisplacements; + for (auto fullForcedDisplacement : simulationJobOfFullModel.nodalForcedDisplacements) { + reducedNodalForcedDisplacements[fullToReducedMap.at(fullForcedDisplacement.first)] + = fullForcedDisplacement.second; } - // std::unordered_map - // 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 &simulationScenarios, const std::unordered_map &reducedToFullInterfaceViMap) { - FormFinder simulator; + DRMSimulationModel simulator; std::shared_ptr 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 &pReducedPatternSimulationJob, const std::unordered_map &reducedToFullInterfaceViMap) { - FormFinder simulator; + DRMSimulationModel simulator; std::shared_ptr 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{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{0, 1, 2, 3, 4, 5}; - // } scenarios[scenarioName] = std::make_shared( SimulationJob(pMesh, simulationScenarioStrings[scenarioName], fixedVertices, nodalForces, {})); @@ -772,52 +780,16 @@ ReducedModelOptimizer::createScenarios( std::unordered_set{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{0, 1, 2, 3, 4, 5}; - // } scenarios[scenarioName] = std::make_shared( 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{0, 1, 2, 3, 4, 5}; - // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; - // } else { - // fixedVertices[viPair.first] = std::unordered_set{0, 1, 2}; - // fixedVertices[viPair.second] = std::unordered_set{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{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 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{0, 1, 2}; - fixedVertices[viPair.second] = std::unordered_set{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{2}; fixedVertices[viPair.second] = std::unordered_set{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(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{2}; fixedVertices[viPair.second] = std::unordered_set{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( @@ -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::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 &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(reducedPatternSimulationJob); } + // global.fullPatternSimulationJobs[0]->pMesh->unregister(); if (global.optimizationSettings.normalizationStrategy != Settings::NormalizationStrategy::NonNormalized) { diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 4b3e5e7..4ea8f4e 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -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 #include "reducedmodeloptimizer_structs.hpp" +#include "simulationmesh.hpp" +#include #ifdef POLYSCOPE_DEFINED #include "polyscope/color_management.h" @@ -16,31 +16,32 @@ using FullPatternVertexIndex = VertexIndex; using ReducedPatternVertexIndex = VertexIndex; -class ReducedModelOptimizer { - - std::shared_ptr m_pReducedPatternSimulationMesh; - std::shared_ptr m_pFullPatternSimulationMesh; - std::unordered_map - m_fullToReducedInterfaceViMap; - std::unordered_map - m_fullPatternOppositeInterfaceViMap; - std::unordered_map nodeToSlot; - std::unordered_map> slotToNode; +class ReducedModelOptimizer +{ + std::shared_ptr m_pReducedPatternSimulationMesh; + std::shared_ptr m_pFullPatternSimulationMesh; + std::unordered_map + m_fullToReducedInterfaceViMap; + std::unordered_map + m_fullPatternOppositeInterfaceViMap; + std::unordered_map nodeToSlot; + std::unordered_map> 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 &simulationScenarios = - std::vector()); + ReducedModelOptimization::Results optimize( + const ReducedModelOptimization::Settings &xRanges, + const std::vector &simulationScenarios + = std::vector()); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); static void computeReducedModelSimulationJob( const SimulationJob &simulationJobOfFullModel, - const std::unordered_map &simulationJobFullToReducedMap, + const std::unordered_map &fullToReducedMap, SimulationJob &simulationJobOfReducedModel); SimulationJob @@ -75,15 +76,12 @@ public: &fullToReducedInterfaceViMap, std::unordered_map &fullPatternOppositeInterfaceViMap); - static void - visualizeResults(const std::vector> - &fullPatternSimulationJobs, - const std::vector> - &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, - const std::unordered_map - &reducedToFullInterfaceViMap); + static void visualizeResults( + const std::vector> &fullPatternSimulationJobs, + const std::vector> &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap); static void registerResultsForDrawing( const std::shared_ptr &pFullPatternSimulationJob, const std::shared_ptr &pReducedPatternSimulationJob, @@ -91,40 +89,39 @@ public: FullPatternVertexIndex> &reducedToFullInterfaceViMap); - static double - computeRawError(const std::vector &reducedPatternDisplacements, - const std::vector &fullPatternDisplacements, - const std::unordered_map - &reducedToFullInterfaceViMap); + static double computeRawDisplacementError( + const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap); - static double - computeError(const std::vector &reducedPatternDisplacements, - const std::vector &fullPatternDisplacements, - const std::unordered_map - &reducedToFullInterfaceViMap, - const double &normalizationFactor); + static double computeDisplacementError( + const std::vector &reducedPatternDisplacements, + const std::vector &fullPatternDisplacements, + const std::unordered_map + &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 &displacementsReducedToFullMap, - Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static ReducedModelOptimization::Results runOptimization(const ReducedModelOptimization::Settings &settings); - std::vector> - createScenarios(const std::shared_ptr &pMesh); - void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern, - const std::unordered_set &reducedModelExcludedEges); - void createSimulationMeshes(PatternGeometry &fullModel, - PatternGeometry &reducedModel); - static void - initializeOptimizationParameters(const std::shared_ptr &mesh, const int &optimizationParamters); + static void computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); + static ReducedModelOptimization::Results runOptimization( + const ReducedModelOptimization::Settings &settings); + std::vector> createScenarios( + const std::shared_ptr &pMesh); + void computeMaps(PatternGeometry &fullModel, + PatternGeometry &reducedPattern, + const std::unordered_set &reducedModelExcludedEges); + void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); + static void initializeOptimizationParameters(const std::shared_ptr &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