diff --git a/CMakeLists.txt b/CMakeLists.txt index fec2faa..7ce4615 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,7 +56,11 @@ find_package(OpenMP) set(MYSOURCESDIR "/home/iason/Coding/Libraries/MySources") file(GLOB MYSOURCES ${MYSOURCESDIR}/*.hpp ${MYSOURCESDIR}/*.cpp) -add_executable(${PROJECT_NAME} ${SOURCES} ${MYSOURCES} ) +file(GLOB EXT_SOURCES ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp + ) + + +add_executable(${PROJECT_NAME} ${SOURCES} ${MYSOURCES} ${EXT_SOURCES} ) target_include_directories(${PROJECT_NAME} PRIVATE ${vcglib_devel_SOURCE_DIR} diff --git a/src/main.cpp b/src/main.cpp index 1d03182..7e38d96 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,37 +13,79 @@ #include int main(int argc, char *argv[]) { - FlatPattern pattern("/home/iason/Models/valid_6777.ply"); + // FlatPattern pattern("/home/iason/Models/valid_6777.ply"); // FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply"); - pattern.createFan(); // pattern.savePly("fannedValid.ply"); - FlatPattern reducedModel("/home/iason/Models/reducedModel.ply"); - reducedModel.createFan(); - // reducedModel.savePly("fannedReduced.ply"); + registerWorldAxes(); + // Create reduced models + const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; + std::vector singleBarReducedModelEdges{vcg::Point2i(0, 3)}; + FlatPattern singleBarReducedModel(numberOfNodesPerSlot, + singleBarReducedModelEdges); + singleBarReducedModel.setLabel("Single bar reduced model"); - std::unordered_map reducedToFullViMap{ - {1, 3}, {2, 7}, {3, 11}, {4, 15}, {5, 19}, {6, 23}}; - ReducedModelOptimizer optimizer(pattern, reducedModel, reducedToFullViMap, - {}); // 6 - // fan reduced model - // {{1, 3}, {2, 15}}, {}); // 2 fan reduced model - Eigen::VectorXd optimalParameters = optimizer.optimize(); + std::vector CCWreducedModelEdges{vcg::Point2i(1, 5), + vcg::Point2i(3, 5)}; + FlatPattern CWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges); + CWReducedModel.setLabel("CW reduced model"); - // Full model simulation - std::unordered_map> fixedVertices; - // fixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; - fixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; - // fixedVertices[7] = std::unordered_set{0, 1, 2}; - std::unordered_map nodalForces{ - {15, {0, 0, 2000, 0, 0, 0}}}; - SimulationJob fullModelSimulationJob{ - std::make_shared(pattern), fixedVertices, nodalForces, {}}; - Simulator formFinder; - SimulationResults fullModelResults = - formFinder.executeSimulation(fullModelSimulationJob); - fullModelResults.simulationLabel = "Full Model"; - fullModelResults.draw(fullModelSimulationJob); + std::vector CWreducedModelEdges{vcg::Point2i(1, 5), + vcg::Point2i(3, 1)}; + FlatPattern CCWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges); + CCWReducedModel.setLabel("CCW reduced model"); + + std::vector reducedModels{&singleBarReducedModel, + &CWReducedModel, &CCWReducedModel}; + + ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + std::string fullPatternsTestSetDirectory = + "/home/iason/Models/TestSet_validPatterns"; + for (const auto &entry : + filesystem::directory_iterator(fullPatternsTestSetDirectory)) { + const auto filepath = + // entry.path(); + std::filesystem::path("/home/iason/Models/valid_6777.ply"); + const auto filepathString = filepath.string(); + // Use only the base triangle version + const std::string tiledSuffix = "_tiled.ply"; + if (filepathString.compare(filepathString.size() - tiledSuffix.size(), + tiledSuffix.size(), tiledSuffix) == 0) { + continue; + } + FlatPattern pattern(filepathString); + pattern.setLabel(filepath.stem().string()); + std::cout << "Testing Pattern:" << filepathString << std::endl; + for (FlatPattern *pReducedModel : reducedModels) { + pReducedModel = reducedModels[0]; + std::unordered_set optimizationExcludedEi; + if (pReducedModel != + reducedModels[0]) { // assumes that the singleBar reduced model is the + // first in the reducedModels vector + optimizationExcludedEi.insert(0); + } + optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi); + optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi); + Eigen::VectorXd optimalParameters = optimizer.optimize(); + } + } + + // // Full model simulation + // std::unordered_map> + // fixedVertices; + // // fixedVertices[0] = std::unordered_set{0, 1, 2, 3, 4, 5}; + // fixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; + // // fixedVertices[7] = std::unordered_set{0, 1, 2}; + // std::unordered_map nodalForces{ + // {15, {0, 0, 2000, 0, 0, 0}}}; + // SimulationJob fullModelSimulationJob{ + // std::make_shared(pattern), fixedVertices, nodalForces, + // {}}; + // Simulator formFinder; + // SimulationResults fullModelResults = + // formFinder.executeSimulation(fullModelSimulationJob); + // fullModelResults.simulationLabel = "Full Model"; + // fullModelResults.draw(fullModelSimulationJob); // fullModelSimulationJob.draw(); // double stiffnessFactor = 1; diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 8f9fc4d..c026a9f 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -1,19 +1,24 @@ -#include "reducedmodeloptimizer.hpp" +#include "reducedmodeloptimizer.hpp" #include "bobyqa.h" +#include "flatpattern.hpp" #include "gradientDescent.h" #include "matplot/matplot.h" #include "simulationhistoryplotter.hpp" +#include "trianglepattterntopology.hpp" -const bool draw = true; -size_t numberOfOptimizationRounds{0}; -Simulator simulator; -Eigen::MatrixX3d optimalReducedModelDisplacements; -SimulationJob reducedModelSimulationJob; +const bool gShouldDraw = true; +size_t g_numberOfOptimizationRounds{0}; +FormFinder simulator; +Eigen::MatrixX3d g_optimalReducedModelDisplacements; +SimulationJob gReducedPatternSimulationJob; std::unordered_map - reducedToFullVertexIndexMap; -matplot::line_handle plotHandle; -std::vector fValueHistory; -Eigen::Vector4d initialParameters; + g_reducedToFullInterfaceViMap; +matplot::line_handle gPlotHandle; +std::vector gObjectiveValueHistory; +Eigen::Vector4d g_initialX; +std::unordered_set g_reducedPatternExludedEdges; +Eigen::MatrixX4d g_initialStiffnessFactors; + struct OptimizationCallback { double operator()(const size_t &iterations, const Eigen::VectorXd &x, const double &fval, Eigen::VectorXd &gradient) const { @@ -21,10 +26,10 @@ struct OptimizationCallback { // SimulationResults reducedModelResults = // simulator.executeSimulation(reducedModelSimulationJob); // reducedModelResults.draw(reducedModelSimulationJob); - fValueHistory.push_back(fval); - auto xPlot = - matplot::linspace(0, fValueHistory.size(), fValueHistory.size()); - plotHandle = matplot::scatter(xPlot, fValueHistory); + gObjectiveValueHistory.push_back(fval); + auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + gObjectiveValueHistory.size()); + gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); // const std::string plotImageFilename = "objectivePlot.png"; // matplot::save(plotImageFilename); // if (numberOfOptimizationRounds % 30 == 0) { @@ -50,24 +55,21 @@ struct Objective { // drawSimulationJob(simulationJob); // Set mesh from x - std::shared_ptr reducedModel = - reducedModelSimulationJob.mesh; + std::shared_ptr reducedModel = + gReducedPatternSimulationJob.mesh; for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { - BeamFormFinder::Element &e = reducedModel->elements[ei]; - e.properties.A = initialParameters(0) * x(0); - e.properties.J = initialParameters(1) * x(1); - e.properties.I2 = initialParameters(2) * x(2); - e.properties.I3 = initialParameters(3) * x(3); - e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; - e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; - e.firstBendingConstFactor = - 2 * e.properties.E * e.properties.I2 / e.initialLength; - e.secondBendingConstFactor = - 2 * e.properties.E * e.properties.I3 / e.initialLength; + if (g_reducedPatternExludedEdges.contains(ei)) { + continue; + } + Element &e = reducedModel->elements[ei]; + e.axialConstFactor *= x(0); + e.torsionConstFactor *= x(1); + e.firstBendingConstFactor *= x(2); + e.secondBendingConstFactor *= x(3); } // run simulation SimulationResults reducedModelResults = - simulator.executeSimulation(reducedModelSimulationJob); + simulator.executeSimulation(gReducedPatternSimulationJob); // std::stringstream ss; // ss << x; // reducedModelResults.simulationLabel = ss.str(); @@ -77,7 +79,7 @@ struct Objective { // std::filesystem::current_path().append("Results")); // compute error and return it double error = 0; - for (const auto reducedFullViPair : reducedToFullVertexIndexMap) { + for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { VertexIndex reducedModelVi = reducedFullViPair.first; Eigen::Vector3d vertexDisplacement( reducedModelResults.displacements[reducedModelVi][0], @@ -85,7 +87,7 @@ struct Objective { reducedModelResults.displacements[reducedModelVi][2]); Eigen::Vector3d errorVector = Eigen::Vector3d( - optimalReducedModelDisplacements.row(reducedModelVi)) - + g_optimalReducedModelDisplacements.row(reducedModelVi)) - vertexDisplacement; error += errorVector.norm(); } @@ -99,76 +101,217 @@ double objective(long n, const double *x) { for (size_t xi = 0; xi < n; xi++) { eigenX(xi) = x[xi]; } - // Eigen::VectorXd emptyGradient; - // return operator()(eigenX, emptyGradient); - // drawSimulationJob(simulationJob); - // Set mesh from x - std::shared_ptr reducedModel = reducedModelSimulationJob.mesh; - for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) { - BeamFormFinder::Element &e = reducedModel->elements[ei]; - e.properties.A = initialParameters(0) * eigenX(0); - e.properties.J = initialParameters(1) * eigenX(1); - e.properties.I2 = initialParameters(2) * eigenX(2); - e.properties.I3 = initialParameters(3) * eigenX(3); - e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; - e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; - e.firstBendingConstFactor = - 2 * e.properties.E * e.properties.I2 / e.initialLength; - e.secondBendingConstFactor = - 2 * e.properties.E * e.properties.I3 / e.initialLength; + std::shared_ptr reducedPattern = + gReducedPatternSimulationJob.mesh; + for (EdgeIndex ei = 0; ei < reducedPattern->EN(); ei++) { + Element &e = reducedPattern->elements[ei]; + if (g_reducedPatternExludedEdges.contains(ei)) { + continue; + } + e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * eigenX(0); + e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * eigenX(1); + e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * eigenX(2); + e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * eigenX(3); } // run simulation SimulationResults reducedModelResults = - simulator.executeSimulation(reducedModelSimulationJob); - // std::stringstream ss; - // ss << x; - // reducedModelResults.simulationLabel = ss.str(); - // SimulationResultsReporter resultsReporter; - // resultsReporter.reportResults( - // {reducedModelResults}, - // std::filesystem::current_path().append("Results")); + simulator.executeSimulation(gReducedPatternSimulationJob, false, false); + // compute error and return it double error = 0; - for (const auto reducedFullViPair : reducedToFullVertexIndexMap) { + for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) { VertexIndex reducedModelVi = reducedFullViPair.first; Eigen::Vector3d vertexDisplacement( reducedModelResults.displacements[reducedModelVi][0], reducedModelResults.displacements[reducedModelVi][1], reducedModelResults.displacements[reducedModelVi][2]); Eigen::Vector3d errorVector = - Eigen::Vector3d(optimalReducedModelDisplacements.row(reducedModelVi)) - + Eigen::Vector3d( + g_optimalReducedModelDisplacements.row(reducedModelVi)) - vertexDisplacement; error += errorVector.norm(); } + gObjectiveValueHistory.push_back(error); + auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + gObjectiveValueHistory.size()); + gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); + return error; } +void ReducedModelOptimizer::computeMaps( + FlatPattern &fullPattern, FlatPattern &reducedPattern, + const std::unordered_set &reducedModelExcludedEges) { + // Compute the offset between the interface nodes + const size_t interfaceSlotIndex = 4; // bottom edge + assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && + slotToNode.find(interfaceSlotIndex)->second.size() == 1); + // Assuming that in the bottom edge there is only one vertex which is also the + // interface + const size_t baseTriangleInterfaceVi = + *(slotToNode.find(interfaceSlotIndex)->second.begin()); + + vcg::tri::Allocator::PointerUpdater + pu_fullModel; + fullPattern.deleteDanglingVertices(pu_fullModel); + const size_t fullModelBaseTriangleInterfaceVi = + pu_fullModel.remap.empty() ? baseTriangleInterfaceVi + : pu_fullModel.remap[baseTriangleInterfaceVi]; + const size_t fullModelBaseTriangleVN = fullPattern.VN(); + fullPattern.createFan(); + const size_t duplicateVerticesPerFanPair = + fullModelBaseTriangleVN - fullPattern.VN() / 6; + const size_t fullPatternInterfaceVertexOffset = + fullModelBaseTriangleVN - duplicateVerticesPerFanPair; + // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << + // std::endl; + + // Save excluded edges + g_reducedPatternExludedEdges.clear(); + const size_t fanSize = 6; + const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); + for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { + for (const size_t ei : reducedModelExcludedEges) { + g_reducedPatternExludedEdges.insert( + fanIndex * reducedBaseTriangleNumberOfEdges + ei); + } + } + + // Construct reduced->full and full->reduced interface vi map + g_reducedToFullInterfaceViMap.clear(); + 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(); + for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { + g_reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + + reducedModelBaseTriangleInterfaceVi] = + fullModelBaseTriangleInterfaceVi + + fanIndex * fullPatternInterfaceVertexOffset; + } + m_fullToReducedInterfaceViMap.clear(); + constructInverseMap(g_reducedToFullInterfaceViMap, + m_fullToReducedInterfaceViMap); + + // fullPattern.setLabel("FullPattern"); + // reducedPattern.setLabel("ReducedPattern"); + // Create opposite vertex map + m_fullPatternOppositeInterfaceViMap.clear(); + for (size_t fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) { + const size_t vi0 = fullModelBaseTriangleInterfaceVi + + fanIndex * fullPatternInterfaceVertexOffset; + const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; + assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); + m_fullPatternOppositeInterfaceViMap[vi0] = vi1; + } + + const bool debugMapping = false; + if (debugMapping) { + reducedPattern.registerForDrawing(); + std::vector colors_reducedPatternExcludedEdges( + reducedPattern.EN(), glm::vec3(0, 0, 0)); + for (const size_t ei : g_reducedPatternExludedEdges) { + colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); + } + const std::string label = reducedPattern.getLabel(); + polyscope::getCurveNetwork(label) + ->addEdgeColorQuantity("Excluded edges", + colors_reducedPatternExcludedEdges) + ->setEnabled(true); + polyscope::show(); + + std::vector nodeColorsOpposite(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (const std::pair oppositeVerts : + m_fullPatternOppositeInterfaceViMap) { + auto color = polyscope::getNextUniqueColor(); + nodeColorsOpposite[oppositeVerts.first] = color; + nodeColorsOpposite[oppositeVerts.second] = color; + } + fullPattern.registerForDrawing(); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite) + ->setEnabled(true); + polyscope::show(); + + std::vector nodeColorsReducedToFull_reduced(reducedPattern.VN(), + glm::vec3(0, 0, 0)); + std::vector nodeColorsReducedToFull_full(fullPattern.VN(), + glm::vec3(0, 0, 0)); + for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { + if (g_reducedToFullInterfaceViMap.contains(vi)) { + + auto color = polyscope::getNextUniqueColor(); + nodeColorsReducedToFull_reduced[vi] = color; + nodeColorsReducedToFull_full[g_reducedToFullInterfaceViMap[vi]] = color; + } + } + polyscope::getCurveNetwork(reducedPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_reduced", + nodeColorsReducedToFull_reduced) + ->setEnabled(true); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_full", + nodeColorsReducedToFull_full) + ->setEnabled(true); + polyscope::show(); + } +} + +void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, + FlatPattern &reducedModel) { + pReducedModelElementalMesh = std::make_shared(reducedModel); + pFullModelElementalMesh = std::make_shared(fullModel); +} + ReducedModelOptimizer::ReducedModelOptimizer( - ConstVCGEdgeMesh &fullModel, ConstVCGEdgeMesh &reducedModel, - const std::unordered_map - &viMap, - const std::vector> &oppositeVertices) - : reducedToFullViMap(viMap) { - reducedToFullVertexIndexMap = viMap; + const std::vector &numberOfNodesPerSlot) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); + FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); +} - // construct fullToReducedViMap - for (const std::pair &viPair : - viMap) { - fullToReducedViMap[viPair.second] = viPair.first; +void ReducedModelOptimizer::initialize( + FlatPattern &fullPattern, FlatPattern &reducedPattern, + const std::unordered_set &reducedModelExcludedEdges) { + assert(fullPattern.VN() == reducedPattern.VN() && + fullPattern.EN() > reducedPattern.EN()); + polyscope::removeAllStructures(); + // Create copies of the input models + FlatPattern copyFullPattern; + FlatPattern copyReducedPattern; + copyFullPattern.copy(fullPattern); + copyReducedPattern.copy(reducedPattern); + computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); + createSimulationMeshes(copyFullPattern, copyReducedPattern); + initializeStiffnesses(); + + int i = 0; + i++; +} + +void ReducedModelOptimizer::initializeStiffnesses() { + g_initialStiffnessFactors.resize(pReducedModelElementalMesh->EN(), 4); + // Save save the beam stiffnesses + for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { + Element &e = pReducedModelElementalMesh->elements[ei]; + // if (g_reducedPatternExludedEdges.contains(ei)) { + // const double stiffnessFactor = 1; + // e.axialConstFactor *= stiffnessFactor; + // e.torsionConstFactor *= stiffnessFactor; + // e.firstBendingConstFactor *= stiffnessFactor; + // e.secondBendingConstFactor *= stiffnessFactor; + // } + g_initialStiffnessFactors(ei, 0) = e.axialConstFactor; + g_initialStiffnessFactors(ei, 1) = e.torsionConstFactor; + g_initialStiffnessFactors(ei, 2) = e.firstBendingConstFactor; + g_initialStiffnessFactors(ei, 3) = e.secondBendingConstFactor; } - - // construct opposite vertex map - for (const std::pair &viPair : oppositeVertices) { - oppositeVertexMap[viPair.first] = viPair.second; - oppositeVertexMap[viPair.second] = viPair.first; - } - - reducedModelElementalMesh = std::make_shared(reducedModel); - pFullModelElementalMesh = std::make_shared(fullModel); - optimalReducedModelDisplacements.resize(reducedModelElementalMesh->VN(), 3); } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -177,18 +320,18 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( std::unordered_map> reducedModelFixedVertices; for (auto fullModelFixedVertex : simulationJobOfFullModel.fixedVertices) { - reducedModelFixedVertices[fullToReducedViMap.at( + reducedModelFixedVertices[m_fullToReducedInterfaceViMap.at( fullModelFixedVertex.first)] = fullModelFixedVertex.second; } std::unordered_map reducedModelNodalForces; for (auto fullModelNodalForce : simulationJobOfFullModel.nodalExternalForces) { - reducedModelNodalForces[fullToReducedViMap.at(fullModelNodalForce.first)] = - fullModelNodalForce.second; + reducedModelNodalForces[m_fullToReducedInterfaceViMap.at( + fullModelNodalForce.first)] = fullModelNodalForce.second; } - simulationJobOfReducedModel = SimulationJob{reducedModelElementalMesh, + simulationJobOfReducedModel = SimulationJob{pReducedModelElementalMesh, reducedModelFixedVertices, reducedModelNodalForces, {}}; @@ -204,12 +347,19 @@ SimulationJob ReducedModelOptimizer::getReducedSimulationJob( void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, - Eigen::MatrixX3d &optimaldisplacementsOfReducedModel) { - for (auto reducedFullViPair : reducedToFullVertexIndexMap) { + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { + + optimalDisplacementsOfReducedModel.resize(pReducedModelElementalMesh->VN(), + 3); + optimalDisplacementsOfReducedModel.setZero( + optimalDisplacementsOfReducedModel.rows(), + optimalDisplacementsOfReducedModel.cols()); + + for (auto reducedFullViPair : g_reducedToFullInterfaceViMap) { const VertexIndex fullModelVi = reducedFullViPair.second; - Vector6d fullModelViDisplacements = + const Vector6d fullModelViDisplacements = fullModelResults.displacements[fullModelVi]; - optimaldisplacementsOfReducedModel.row(reducedFullViPair.first) = + optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) = Eigen::Vector3d(fullModelViDisplacements[0], fullModelViDisplacements[1], fullModelViDisplacements[2]); @@ -218,27 +368,25 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( const SimulationJob &fullModelSimulationJob) { + + gObjectiveValueHistory.clear(); SimulationResults fullModelResults = simulator.executeSimulation(fullModelSimulationJob, false, false); fullModelResults.simulationLabel = "fullModel"; - if (draw) { - drawWorldAxes(); - fullModelResults.draw(fullModelSimulationJob); - } computeDesiredReducedModelDisplacements(fullModelResults, - optimalReducedModelDisplacements); + g_optimalReducedModelDisplacements); computeReducedModelSimulationJob(fullModelSimulationJob, - reducedModelSimulationJob); + gReducedPatternSimulationJob); + // fullModelSimulationJob.registerForDrawing(); + // polyscope::show(); + // gReducedPatternSimulationJob.registerForDrawing(); + // polyscope::show(); + fullModelResults.registerForDrawing(fullModelSimulationJob); + polyscope::show(); // Set initial guess of solution Eigen::VectorXd initialGuess(4); - auto propertiesOfFirstBeamOfFullModel = - pFullModelElementalMesh->elements[0].properties; - initialParameters(0) = propertiesOfFirstBeamOfFullModel.A; - initialParameters(1) = propertiesOfFirstBeamOfFullModel.J; - initialParameters(2) = propertiesOfFirstBeamOfFullModel.I2; - initialParameters(3) = propertiesOfFirstBeamOfFullModel.I3; - const double stifnessFactor = 2; + const double stifnessFactor = 1; initialGuess(0) = stifnessFactor; initialGuess(1) = stifnessFactor; initialGuess(2) = stifnessFactor; @@ -282,9 +430,10 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( std::cout << "Final xval: " << result.xval.transpose() << std::endl; SimulationResults reducedModelOptimizedResults = - simulator.executeSimulation(reducedModelSimulationJob); + simulator.executeSimulation(gReducedPatternSimulationJob); reducedModelOptimizedResults.simulationLabel = "reducedModel"; - reducedModelOptimizedResults.draw(reducedModelSimulationJob); + reducedModelOptimizedResults.registerForDrawing( + gReducedPatternSimulationJob); return result.xval; } else { // use bobyqa @@ -294,8 +443,9 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2); assert(npt < 2 * n + 1 && "The choice of the number of interpolation " "conditions is not recommended."); - std::vector x{initialGuess(0), initialGuess(1), initialGuess(2), - initialGuess(3)}; + std::vector x + // {1.03424, 0.998456, 0.619916, -0.202997}; + {initialGuess(0), initialGuess(1), initialGuess(2), initialGuess(3)}; std::vector xLow(x.size(), -100); std::vector xUpper(x.size(), 100); const double maxX = *std::max_element( @@ -318,81 +468,147 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( } SimulationResults reducedModelOptimizedResults = - simulator.executeSimulation(reducedModelSimulationJob); + simulator.executeSimulation(gReducedPatternSimulationJob); reducedModelOptimizedResults.simulationLabel = "reducedModel"; - reducedModelOptimizedResults.draw(reducedModelSimulationJob); + reducedModelOptimizedResults.registerForDrawing( + gReducedPatternSimulationJob); + polyscope::show(); return eigenX; } } std::vector ReducedModelOptimizer::createScenarios( - const std::shared_ptr &pMesh) { + const std::shared_ptr &pMesh) { std::vector scenarios; - // Pull up - ////1Vertex - std::unordered_map> - pullUp1_FixedVertices; - pullUp1_FixedVertices[3] = std::unordered_set{0, 1, 2, 3, 4, 5}; - std::unordered_map pullUp1_NodalForces{ - {15, {0, 0, 100, 0, 0, 0}}}; - scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, {}}); + std::unordered_map> fixedVertices; + std::unordered_map nodalForces; + const double forceMagnitude = 250; - // /// 2Vertices - // std::unordered_map> - // pullUp1_FixedVertices; - // saddleFixedVertices[3] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; - // std::unordered_map pullUp1_NodalForces{ - // {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}}, - // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, - // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; - // scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, - // {}}); - - /// 3Vertices - // std::unordered_map> - // pullUp1_FixedVertices; - // saddleFixedVertices[3] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[11] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[15] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[19] = std::unordered_set{0, 1, 2}; - // saddleFixedVertices[23] = std::unordered_set{0, 1, 2}; - // std::unordered_map pullUp1_NodalForces{ - // {15, {0, 0, 0, 0, -900, 0}}, {3, {0, 0, 0, 0, 900, 0}}, - // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, - // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; - // scenarios.push_back({pMesh, pullUp1_FixedVertices, pullUp1_NodalForces, - // {}}); - - // Single moment - // std::unordered_map> - // fixedVertices; fixedVertices[3] = std::unordered_set{0, 1, 2, 3, - // 4, 5}; fixedVertices[15] = std::unordered_set{0, 1, 2}; - // std::unordered_map nodalForces{ - // {15, {0, 0, 0, 0, 700, 0}}}; - // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}}); + // // Axial + // 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; + // fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + // // In-plane Bending + // // Assuming the patterns lay on the x-y plane + // const CoordType patternPlane(0, 0, 1); + // fixedVertices.clear(); + // nodalForces.clear(); + // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { + // CoordType v = + // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // .Normalize(); + // CoordType forceDirection = patternPlane ^ v; + // nodalForces[viPair.first] = Vector6d({forceDirection[0], + // forceDirection[1], + // forceDirection[2], 0, 0, 0}) * + // forceMagnitude; + // fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // scenarios.push_back({pMesh, 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(); + // nodalForces[viPair.first] = + // Vector6d({0, 0, 0, v[1], v[0], 0}) * forceMagnitude; + // } else { + // fixedVertices[viPair.first] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // fixedVertices[viPair.second] = + // std::unordered_set{0, 1, 2, 3, 4, 5}; + // } + // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + // } + // Out-of-plane bending . Pull towards Z + fixedVertices.clear(); + nodalForces.clear(); + for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { + nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}); + fixedVertices[viPair.second] = + std::unordered_set{0, 1, 2, 3, 4, 5}; + } + scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + // // Dou?? + // fixedVertices.clear(); + // nodalForces.clear(); + // for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { + // CoordType v = + // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // .Normalize(); + // CoordType momentDirection = patternPlane ^ v; + // nodalForces[viPair.first] = + // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * + // forceMagnitude; + // fixedVertices[viPair.first] = std::unordered_set{2}; + // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + // } + // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + // // Saddle + // fixedVertices.clear(); + // nodalForces.clear(); + // for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); + // viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + // const auto &viPair = *viPairIt; + // CoordType v = + // (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // .Normalize(); + // CoordType momentDirection = patternPlane ^ v; + // if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + // nodalForces[viPair.first] = + // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3 + // * forceMagnitude; + // nodalForces[viPair.second] = + // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3 + // * + // (-forceMagnitude); + // } else { + // fixedVertices[viPair.first] = std::unordered_set{2}; + // nodalForces[viPair.first] = + // Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * + // (-forceMagnitude); + // fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + // } + // } // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); - // Saddle // std::unordered_map> - // saddle_fixedVertices; + // saddleFixedVertices; // // saddle_fixedVertices[3] = std::unordered_set{0, 1, 2}; - // saddle_fixedVertices[7] = std::unordered_set{0, 1, 2}; - // saddle_fixedVertices[11] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[7] = std::unordered_set{0, 1, 2}; + // saddleFixedVertices[11] = std::unordered_set{0, 1, 2}; // // saddle_fixedVertices[15] = std::unordered_set{0, 1, 2}; // // saddle_fixedVertices[19] = std::unordered_set{0, 1, 2}; // // saddle_fixedVertices[23] = std::unordered_set{0, 1, 2}; - // std::unordered_map saddle_nodalForces{ - // {15, {0, 0, 0, 0, -1400, 0}}, {3, {0, 0, 0, 0, 1400, 0}}, - // {7, {0, 0, 0, 700, 0, 0}}, {11, {0, 0, 0, 700, 0, 0}}, - // {19, {0, 0, 0, -700, 0, 0}}, {23, {0, 0, 0, -700, 0, 0}}}; - // // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, - // {}}); std::unordered_map - // saddle_forcedDisplacements; scenarios.push_back({pMesh, - // saddle_fixedVertices, saddle_nodalForces, - // saddle_forcedDisplacements}); + // std::unordered_map saddleNodalForces{ + // {15, {0, 0, 0, 0, -4 * 90, 0}}, {3, {0, 0, 0, 0, 4 * 90, 0}}, + // {7, {0, 0, 0, 4 * 70, 0, 0}}, {11, {0, 0, 0, 4 * 70, 0, 0}}, + // {19, {0, 0, 0, -4 * 70, 0, 0}}, {23, {0, 0, 0, -4 * 70, 0, 0}}}; + // scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}}); return scenarios; } @@ -402,7 +618,7 @@ Eigen::VectorXd ReducedModelOptimizer::optimize() { createScenarios(pFullModelElementalMesh); std::vector results; for (const SimulationJob &job : simulationJobs) { - // fullModelSimulationJob.draw(); + polyscope::removeAllStructures(); auto result = optimizeForSimulationJob(job); results.push_back(result); } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 675d87f..c29d5e3 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -9,34 +9,21 @@ using FullModelVertexIndex = VertexIndex; using ReducedModelVertexIndex = VertexIndex; -struct SimulationScenarios { - std::vector scenarios; - std::shared_ptr pMesh; - SimulationScenarios(std::shared_ptr pMesh) : pMesh(pMesh) {} - - std::vector getSimulationJobs() const { return scenarios; } -}; - class ReducedModelOptimizer { - std::shared_ptr reducedModelElementalMesh; - std::shared_ptr pFullModelElementalMesh; - const std::unordered_map - &reducedToFullViMap; + std::shared_ptr pReducedModelElementalMesh; + std::shared_ptr pFullModelElementalMesh; std::unordered_map - fullToReducedViMap; + m_fullToReducedInterfaceViMap; std::unordered_map - oppositeVertexMap; + m_fullPatternOppositeInterfaceViMap; + std::unordered_map nodeToSlot; + std::unordered_map> slotToNode; public: Eigen::VectorXd optimize(); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; - ReducedModelOptimizer( - ConstVCGEdgeMesh &fullModel, ConstVCGEdgeMesh &reducedModel, - const std::unordered_map - &viMap, - const std::vector> &oppositeVertices); + ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); void computeReducedModelSimulationJob( const SimulationJob &simulationJobOfFullModel, SimulationJob &simulationJobOfReducedModel); @@ -44,14 +31,22 @@ public: SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); + void initialize(FlatPattern &fullPattern, FlatPattern &reducedPatterm, + const std::unordered_set &reducedModelExcludedEges); + private: void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, - Eigen::MatrixX3d &optimaldisplacementsOfReducedModel); + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); Eigen::VectorXd optimizeForSimulationJob(const SimulationJob &fullModelSimulationJob); std::vector - createScenarios(const std::shared_ptr &pMesh); + createScenarios(const std::shared_ptr &pMesh); + void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern, + const std::unordered_set &reducedModelExcludedEges); + void createSimulationMeshes(FlatPattern &fullModel, + FlatPattern &reducedModel); + void initializeStiffnesses(); }; #endif // REDUCEDMODELOPTIMIZER_HPP