diff --git a/src/main.cpp b/src/main.cpp index 6812189..99e2949 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -43,6 +43,7 @@ int main(int argc, char *argv[]) { // pattern.savePly("fannedValid.ply"); // Create reduced models + // FormFinder::runUnitTests(); const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; std::vector singleBarReducedModelEdges{vcg::Point2i(0, 3)}; FlatPattern singleBarReducedModel(numberOfNodesPerSlot, @@ -63,25 +64,34 @@ int main(int argc, char *argv[]) { &CWReducedModel, &CCWReducedModel}; ReducedModelOptimizer optimizer(numberOfNodesPerSlot); + + // FlatPattern pattern("/home/iason/Models/TestSet_validPatterns/59.ply"); + // optimizer.initialize(pattern, *reducedModels[1], {}); + // optimizer.setInitialGuess({0.10073746137216238, 1.9611251949927333, + // 7.5940127902102867, 8.0260424171674654}); + // Eigen::VectorXd optimalParameters = optimizer.optimize(); + // optimizer.setInitialGuess({2, 2, 2, 2, 2}); + std::string fullPatternsTestSetDirectory = "/home/iason/Models/TestSet_validPatterns"; for (const auto &entry : filesystem::directory_iterator(fullPatternsTestSetDirectory)) { const auto filepath = // std::filesystem::path("/home/iason/Models/valid_6777.ply"); - // std::filesystem::path( - // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" + // std::filesystem::path( + // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" // "1v_0v_2e_1e_1c_6fan/3/Valid/222.ply"); // std::filesystem::path( // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" // "1v_0v_2e_1e_1c_6fan/2/Valid/33.ply"); - std::filesystem::path( - "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" - "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply"); - // std::filesystem::path( - // "/home/iason/Models/TestSet_validPatterns/865.ply"); - entry.path(); + // std::filesystem::path( + // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/" + // "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply"); + // std::filesystem::path( + // "/home/iason/Models/TestSet_validPatterns/59.ply"); + entry.path(); const auto filepathString = filepath.string(); + std::cout << "Full pattern:" << filepathString << std::endl; // Use only the base triangle version const std::string tiledSuffix = "_tiled.ply"; if (filepathString.compare(filepathString.size() - tiledSuffix.size(), @@ -91,18 +101,23 @@ int main(int argc, char *argv[]) { FlatPattern pattern(filepathString); pattern.setLabel(filepath.stem().string()); std::cout << "Testing Pattern:" << filepathString << std::endl; - for (FlatPattern *pReducedModel : reducedModels) { - // pReducedModel = reducedModels[1]; - 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); - Eigen::VectorXd optimalParameters = optimizer.optimize(); - } + // for (int reducedPatternIndex = 0; + // reducedPatternIndex < reducedModels.size(); + // reducedPatternIndex++) { + // FlatPattern *pReducedModel = reducedModels[reducedPatternIndex]; + std::unordered_set optimizationExcludedEi; + // if (pReducedModel != + // reducedModels[0]) { // assumes that the singleBar reduced model + // // is + // // the + // // first in the reducedModels vector + // optimizationExcludedEi.insert(0); + // } + std::cout << "Full pattern:" << filepathString << std::endl; + std::cout << "Reduced pattern:" << 1 << std::endl; + optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi); + Eigen::VectorXd optimalParameters = optimizer.optimize(); + // } } // // Full model simulation @@ -170,62 +185,65 @@ int main(int argc, char *argv[]) { // stiffnessFactor *= 1.5; // } - // Beam + // // Beam + // // registerWorldAxes(); + // VCGEdgeMesh mesh; + // // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); + // mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with + // flat " + // "patterns/RodModelOptimizationForPatterns/build/Debug/" + // "Fanned_Single bar reduced model.ply"); + // // vcg::Matrix44d R; + // // R.SetRotateDeg(45, {0, 0, 1}); + // // vcg::tri::UpdatePosition::Matrix(mesh, R); + // // mesh.updateEigenEdgeAndVertices(); + // FormFinder formFinder; + // formFinder.runUnitTests(); + // std::unordered_map> + // fixedVertices; fixedVertices[1] = std::unordered_set{2}; + // fixedVertices[2] = std::unordered_set{2}; + // fixedVertices[3] = std::unordered_set{1, 2}; + // fixedVertices[4] = std::unordered_set{2}; + // fixedVertices[5] = std::unordered_set{2}; + // fixedVertices[6] = std::unordered_set{0, 1, 2}; + // // Forced displacements + // std::unordered_map nodalForcedDisplacements; + // // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1); + // std::unordered_map nodalForcedNormal; + + // CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize(); + // CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize(); + // CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize(); + // const double forceMagnitude = 0.4; + // std::unordered_map nodalForces{ + // // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude}, + // // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) * + // forceMagnitude}, + // // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude}, + // // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) * + // forceMagnitude}, + // // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude}, + // {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}}; + // // nodalForcedNormal[10] = v; + // // nodalForcedNormal[0] = -v; + // // fixedVertices[10] = {0, 1}; + + // SimulationJob beamSimulationJob{std::make_shared(mesh), + // fixedVertices, + // nodalForces, + // {}, + // {}}; + // beamSimulationJob.mesh->setBeamMaterial(0.3, 1); + // beamSimulationJob.mesh->setBeamCrossSection( + // RectangularBeamDimensions{0.002, 0.002}); + // beamSimulationJob.registerForDrawing(); // registerWorldAxes(); - VCGEdgeMesh mesh; - // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply"); - mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/build/Debug/" - "Fanned_Single bar reduced model.ply"); - // vcg::Matrix44d R; - // R.SetRotateDeg(45, {0, 0, 1}); - // vcg::tri::UpdatePosition::Matrix(mesh, R); - // mesh.updateEigenEdgeAndVertices(); - FormFinder formFinder; - formFinder.runUnitTests(); - std::unordered_map> fixedVertices; - fixedVertices[1] = std::unordered_set{2}; - fixedVertices[2] = std::unordered_set{2}; - fixedVertices[3] = std::unordered_set{1, 2}; - fixedVertices[4] = std::unordered_set{2}; - fixedVertices[5] = std::unordered_set{2}; - fixedVertices[6] = std::unordered_set{0, 1, 2}; - // Forced displacements - std::unordered_map nodalForcedDisplacements; - // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1); - std::unordered_map nodalForcedNormal; - - CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize(); - CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize(); - CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize(); - const double forceMagnitude = 0.4; - std::unordered_map nodalForces{ - // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude}, - // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) * forceMagnitude}, - // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude}, - // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) * forceMagnitude}, - // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude}, - {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}}; - // nodalForcedNormal[10] = v; - // nodalForcedNormal[0] = -v; - // fixedVertices[10] = {0, 1}; - - SimulationJob beamSimulationJob{std::make_shared(mesh), - fixedVertices, - nodalForces, - {}, - {}}; - beamSimulationJob.mesh->setBeamMaterial(0.3, 1); - beamSimulationJob.mesh->setBeamCrossSection( - RectangularBeamDimensions{0.002, 0.002}); - beamSimulationJob.registerForDrawing(); - registerWorldAxes(); - polyscope::show(); - SimulationResults beamSimulationResults = - formFinder.executeSimulation(beamSimulationJob, true, true); - beamSimulationResults.simulationLabel = "Beam"; - beamSimulationResults.registerForDrawing(beamSimulationJob); - polyscope::show(); + // polyscope::show(); + // SimulationResults beamSimulationResults = + // formFinder.executeSimulation(beamSimulationJob, true, true); + // beamSimulationResults.label = "Beam"; + // beamSimulationResults.registerForDrawing(beamSimulationJob); + // polyscope::show(); return 0; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 88e35ae..63c36bd 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -16,7 +16,8 @@ matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; Eigen::Vector4d g_initialX; std::unordered_set g_reducedPatternExludedEdges; -double g_initialParameters; +// double g_initialParameters; +Eigen::VectorXd g_initialParameters; // struct OptimizationCallback { // double operator()(const size_t &iterations, const Eigen::VectorXd &x, @@ -96,6 +97,7 @@ double g_initialParameters; //}; double ReducedModelOptimizer::objective(long n, const double *x) { + std::cout.precision(17); for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) { std::cout << "x[" + std::to_string(parameterIndex) + "]=" @@ -107,13 +109,14 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // if (g_reducedPatternExludedEdges.contains(ei)) { // continue; // } - e.properties.E = g_initialParameters * x[0]; - // e.properties.G = g_initialParameters(1) * x[1]; - // e.properties.A = g_initialParameters(0, 0) * x[0]; - // e.properties.J = g_initialParameters(0, 1) * x[1]; - // e.properties.I2 = g_initialParameters(0, 2) * x[2]; - // e.properties.I3 = g_initialParameters(0, 3) * x[3]; - e.properties.G = e.properties.E / (2 * (1 + 0.3)); + // e.properties.E = g_initialParameters * x[ei]; + // e.properties.E = g_initialParameters(ei, 0) * x[ei]; + // e.properties.G = g_initialParameters(1) * x[1]; + e.properties.A = g_initialParameters(0) * x[0]; + e.properties.J = g_initialParameters(1) * x[1]; + e.properties.I2 = g_initialParameters(2) * x[2]; + e.properties.I3 = g_initialParameters(3) * x[3]; + // e.properties.G = e.properties.E / (2 * (1 + 0.3)); e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength; e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength; e.firstBendingConstFactor = @@ -145,10 +148,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) { error += errorVector.norm(); } - // gObjectiveValueHistory.push_back(error); - // auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), - // gObjectiveValueHistory.size()); - // gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); + gObjectiveValueHistory.push_back(error); + auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), + gObjectiveValueHistory.size()); + gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory); return error; } @@ -215,7 +218,7 @@ void ReducedModelOptimizer::computeMaps( // reducedPattern.setLabel("ReducedPattern"); // Create opposite vertex map m_fullPatternOppositeInterfaceViMap.clear(); - for (size_t fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) { + for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { const size_t vi0 = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; @@ -278,15 +281,19 @@ void ReducedModelOptimizer::computeMaps( void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel, FlatPattern &reducedModel) { - m_pReducedPatternElementalMesh = + if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { + std::cerr << "A rectangular cross section is expected." << std::endl; + return; + } + m_pReducedPatternSimulationMesh = std::make_shared(reducedModel); - m_pReducedPatternElementalMesh->setBeamCrossSection( + m_pReducedPatternSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - m_pReducedPatternElementalMesh->setBeamMaterial(0.3, 1); - m_pFullModelElementalMesh = std::make_shared(fullModel); - m_pFullModelElementalMesh->setBeamCrossSection( + m_pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1); + m_pFullModelSimulationMesh = std::make_shared(fullModel); + m_pFullModelSimulationMesh->setBeamCrossSection( CrossSectionType{0.002, 0.002}); - m_pFullModelElementalMesh->setBeamMaterial(0.3, 1); + m_pFullModelSimulationMesh->setBeamMaterial(0.3, 1); } ReducedModelOptimizer::ReducedModelOptimizer( @@ -312,7 +319,7 @@ void ReducedModelOptimizer::initialize( } void ReducedModelOptimizer::initializeStiffnesses() { - // g_initialParameters.resize(1, 2); + g_initialParameters.resize(4); // Save save the beam stiffnesses // for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) { // Element &e = pReducedModelElementalMesh->elements[ei]; @@ -323,18 +330,22 @@ void ReducedModelOptimizer::initializeStiffnesses() { // e.firstBendingConstFactor *= stiffnessFactor; // e.secondBendingConstFactor *= stiffnessFactor; // } - g_initialParameters = - m_pReducedPatternElementalMesh->elements[0].properties.E; - // g_initialParameters(1) = - // pReducedModelElementalMesh->elements[0].properties.G; - // g_initialParameters(0, 0) = - // pReducedModelElementalMesh->elements[0].properties.A; - // g_initialParameters(0, 1) = - // pReducedModelElementalMesh->elements[0].properties.J; - // g_initialParameters(0, 2) = - // pReducedModelElementalMesh->elements[0].properties.I2; - // g_initialParameters(0, 3) = - // pReducedModelElementalMesh->elements[0].properties.I3; + // g_initialParameters = + // m_pReducedPatternSimulationMesh->elements[0].properties.E; + // for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) { + // g_initialParameters(ei, 0) = + // m_pReducedPatternSimulationMesh->elements[ei].properties.E; + // } + // g_initialParameters(1) = + // pReducedModelElementalMesh->elements[0].properties.G; + g_initialParameters(0) = + m_pReducedPatternSimulationMesh->elements[0].properties.A; + g_initialParameters(1) = + m_pReducedPatternSimulationMesh->elements[0].properties.J; + g_initialParameters(2) = + m_pReducedPatternSimulationMesh->elements[0].properties.I2; + g_initialParameters(3) = + m_pReducedPatternSimulationMesh->elements[0].properties.I3; // } } @@ -363,7 +374,7 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( fullModelNodalForcedRotation.second; } - simulationJobOfReducedModel = SimulationJob{m_pReducedPatternElementalMesh, + simulationJobOfReducedModel = SimulationJob{m_pReducedPatternSimulationMesh, reducedModelFixedVertices, reducedModelNodalForces, {}, @@ -384,7 +395,7 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( optimalDisplacementsOfReducedModel.resize(0, 3); optimalDisplacementsOfReducedModel.resize( - m_pReducedPatternElementalMesh->VN(), 3); + m_pReducedPatternSimulationMesh->VN(), 3); optimalDisplacementsOfReducedModel.setZero( optimalDisplacementsOfReducedModel.rows(), optimalDisplacementsOfReducedModel.cols()); @@ -407,10 +418,10 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( // polyscope::show(); gObjectiveValueHistory.clear(); fullPatternSimulationJob.mesh->savePly( - "Fanned_" + m_pFullModelElementalMesh->getLabel() + ".ply"); + "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply"); SimulationResults fullModelResults = - simulator.executeSimulation(fullPatternSimulationJob, false, true); - fullModelResults.simulationLabel = "fullModel"; + simulator.executeSimulation(fullPatternSimulationJob, false); + fullModelResults.label = "fullModel"; fullModelResults.registerForDrawing(fullPatternSimulationJob); registerWorldAxes(); polyscope::show(); @@ -420,24 +431,27 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( g_reducedPatternSimulationJob); // gReducedPatternSimulationJob.registerForDrawing(); // polyscope::show(); - fullModelResults.registerForDrawing(fullPatternSimulationJob); - polyscope::show(); double (*pObjectiveFunction)(long, const double *) = &objective; - const size_t n = 1; - // g_initialParameters.rows(); - const size_t npt = ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + // const size_t n = m_pReducedPatternSimulationMesh->EN(); + const size_t n = 4; + const size_t npt = 2 * n; + // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; 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."); // Set initial guess of solution - const double initialGuess = 1; - std::vector x(n, initialGuess); + + std::vector x(n, 2); + if (!initialGuess.empty()) { + x = initialGuess; + } // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764}; // {1, 5.9277}; + // {0.0001, 2, 2.000000005047502, 1.3055270196964464}; // {initialGuess(0), initialGuess(1), initialGuess(2), // initialGuess(3)}; - const double xMin = 1; - const double xMax = 100; + const double xMin = 0.05; + const double xMax = 10; // assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) { // return d >= xMax || d <= xMin; // })); @@ -474,7 +488,7 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( std::cout << "Final objective value:" << error << std::endl; - reducedModelOptimizedResults.simulationLabel = "reducedModel"; + reducedModelOptimizedResults.label = "reducedModel"; reducedModelOptimizedResults.registerForDrawing( g_reducedPatternSimulationJob); polyscope::show(); @@ -486,65 +500,79 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob( return eigenX; } +void ReducedModelOptimizer::setInitialGuess(std::vector v) { + initialGuess = v; +} + std::vector ReducedModelOptimizer::createScenarios( const std::shared_ptr &pMesh) { std::vector scenarios; std::unordered_map> fixedVertices; std::unordered_map nodalForces; + // NOTE: Assuming that the first interface node lays on the y axis const double forceMagnitude = 1; // Assuming the patterns lays on the x-y plane const CoordType patternPlaneNormal(0, 0, 1); // Make the first interface node lay on the x axis - const size_t fullPatternFirstInterfaceNodeIndex = - m_fullPatternOppositeInterfaceViMap.begin()->second; - CoordType fullPatternFirstInterfaceNodePosition = - m_pFullModelElementalMesh->vert[fullPatternFirstInterfaceNodeIndex].cP(); - const vcg::Matrix33d R = vcg::RotationMatrix( - fullPatternFirstInterfaceNodePosition, - CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false); - std::for_each(m_pFullModelElementalMesh->vert.begin(), - m_pFullModelElementalMesh->vert.end(), [&](auto &v) { - v.P() = R * v.P(); - v.N() = R * v.N(); - }); - std::for_each(m_pReducedPatternElementalMesh->vert.begin(), - m_pReducedPatternElementalMesh->vert.end(), [&](auto &v) { - v.P() = R * v.P(); - v.N() = R * v.N(); - }); - m_pFullModelElementalMesh->updateEigenEdgeAndVertices(); - m_pReducedPatternElementalMesh->updateEigenEdgeAndVertices(); - - // //// 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 * 10; - // fixedVertices[viPair.second] = - // std::unordered_set{0, 1, 2, 3, 4, 5}; + // const size_t fullPatternFirstInterfaceNodeIndex = + // m_fullPatternOppositeInterfaceViMap.begin()->second; + // CoordType fullPatternFirstInterfaceNodePosition = + // m_pFullModelSimulationMesh->vert[fullPatternFirstInterfaceNodeIndex].cP(); + // CoordType centerOfMass(0, 0, 0); + // for (size_t vi = 0; vi < pMesh->VN(); vi++) { + // centerOfMass = centerOfMass + pMesh->vert[vi].P(); // } - // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + // centerOfMass /= pMesh->VN(); + // vcg::tri::UpdatePosition::Translate( + // *m_pFullModelSimulationMesh, -centerOfMass); + // vcg::tri::UpdatePosition::Translate( + // *m_pReducedPatternSimulationMesh, centerOfMass); - // //// In-plane Bending - // 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 = (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.push_back({pMesh, fixedVertices, nodalForces, {}}); + // const vcg::Matrix33d R = vcg::RotationMatrix( + // fullPatternFirstInterfaceNodePosition, + // CoordType(fullPatternFirstInterfaceNodePosition.Norm(), 0, 0), false); + // std::for_each(m_pFullModelSimulationMesh->vert.begin(), + // m_pFullModelSimulationMesh->vert.end(), [&](auto &v) { + // v.P() = R * v.P(); + // v.N() = R * v.N(); + // }); + // std::for_each(m_pReducedPatternSimulationMesh->vert.begin(), + // m_pReducedPatternSimulationMesh->vert.end(), [&](auto &v) { + // v.P() = R * v.P(); + // v.N() = R * v.N(); + // }); + + // m_pFullModelSimulationMesh->updateEigenEdgeAndVertices(); + // m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); + + //// 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 * 10; + fixedVertices[viPair.second] = + std::unordered_set{0, 1, 2, 3, 4, 5}; + } + scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + + //// In-plane Bending + 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 = (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.push_back({pMesh, fixedVertices, nodalForces, {}}); // //// Torsion // fixedVertices.clear(); @@ -569,38 +597,38 @@ std::vector ReducedModelOptimizer::createScenarios( // } // 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}) * - // 1; 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}) * 1; + fixedVertices[viPair.second] = + std::unordered_set{0, 1, 2, 3, 4, 5}; + } + scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); - // //// Double using moments - // fixedVertices.clear(); - // nodalForces.clear(); - // 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{1, 2}; - // } 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()) - // .Normalize(); - // nodalForces[viPair.first] = - // Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1; - // nodalForces[viPair.second] = - // Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1; - // } - // scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); + //// Double using moments + fixedVertices.clear(); + nodalForces.clear(); + 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}; + } 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()) + .Normalize(); + nodalForces[viPair.first] = + Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude * 0.1; + nodalForces[viPair.second] = + Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude * 0.1; + } + scenarios.push_back({pMesh, fixedVertices, nodalForces, {}}); //// Saddle fixedVertices.clear(); @@ -633,9 +661,11 @@ std::vector ReducedModelOptimizer::createScenarios( Eigen::VectorXd ReducedModelOptimizer::optimize() { std::vector simulationJobs = - createScenarios(m_pFullModelElementalMesh); + createScenarios(m_pFullModelSimulationMesh); std::vector results; - for (const SimulationJob &job : simulationJobs) { + for (int i = 0; i < simulationJobs.size(); i++) { + std::cout << "Running simulation job:" << i << std::endl; + const SimulationJob &job = simulationJobs[i]; polyscope::removeAllStructures(); auto result = optimizeForSimulationJob(job); results.push_back(result); diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index ceb1b81..8acf9f5 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -11,14 +11,15 @@ using FullPatternVertexIndex = VertexIndex; using ReducedPatternVertexIndex = VertexIndex; class ReducedModelOptimizer { - std::shared_ptr m_pReducedPatternElementalMesh; - std::shared_ptr m_pFullModelElementalMesh; + std::shared_ptr m_pReducedPatternSimulationMesh; + std::shared_ptr m_pFullModelSimulationMesh; std::unordered_map m_fullToReducedInterfaceViMap; std::unordered_map m_fullPatternOppositeInterfaceViMap; std::unordered_map nodeToSlot; std::unordered_map> slotToNode; + std::vector initialGuess; public: Eigen::VectorXd optimize(); @@ -35,6 +36,8 @@ public: void initialize(FlatPattern &fullPattern, FlatPattern &reducedPatterm, const std::unordered_set &reducedModelExcludedEges); + void setInitialGuess(std::vector v); + private: void computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults,