diff --git a/csvfile.hpp b/csvfile.hpp index de28e8c..bf17814 100755 --- a/csvfile.hpp +++ b/csvfile.hpp @@ -1,10 +1,13 @@ #ifndef CSVFILE_HPP #define CSVFILE_HPP +#include +#include #include #include #include #include #include +#include class csvFile; @@ -60,7 +63,46 @@ public: csvFile &operator<<(const std::string &val) { return write(escape(val)); } - template csvFile &operator<<(const T &val) { return write(val); } + template + csvFile &operator<<(const T &val) + { + return write(val); + } + + template + static std::vector> parse(const std::filesystem::path &csvFilepath) + { + std::vector> resultCSV; + if (!std::filesystem::exists(csvFilepath)) { + std::cerr << "The file does not exist:" << csvFilepath.string() << std::endl; + return resultCSV; + } + + std::ifstream inputfile(csvFilepath.string().c_str()); + if (!inputfile.is_open()) { + std::cerr << "Can't open file:" << csvFilepath.string() << std::endl; + return resultCSV; + } + std::vector row; + std::string line; + using Tokenizer = boost::tokenizer>; + while (std::getline(inputfile, line)) { + Tokenizer tokenizer(line); + row.resize(std::distance(tokenizer.begin(), tokenizer.end())); + std::transform(tokenizer.begin(), tokenizer.end(), row.begin(), [](const std::string &el) { + return boost::lexical_cast(el); + }); + // std::cout << std::endl; + // row.assign(tokenizer.begin(), tokenizer.end()); + // for (const auto &el : row) { + // std::cout << el << " "; + // } + // std::cout << std::endl; + resultCSV.push_back(row); + } + + return resultCSV; + } private: template csvFile &write(const T &val) { diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 33b40c5..ddeadb6 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -22,7 +22,7 @@ void DRMSimulationModel::runUnitTests() // mesh.createSpanGrid(spanGridSize); beam.load(std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string()); std::unordered_map> fixedVertices; - fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; + fixedVertices[0] = std::unordered_set{0, 1, 2}; fixedVertices[beam.VN() - 1] = std::unordered_set{1, 2}; std::unordered_map nodalForces{ {beam.VN() / 2, Vector6d({0, 1000, 1000, 0, 0, 0})}}; @@ -45,9 +45,10 @@ void DRMSimulationModel::runUnitTests() settings.Dtini = 0.1; settings.xi = 0.9969; settings.maxDRMIterations = 0; - totalResidualForcesNormThreshold = 1; - // settings.shouldDraw = true; - settings.beVerbose = true; + formFinder.totalResidualForcesNormThreshold = 1; + settings.shouldDraw = false; + settings.beVerbose = false; + // settings.shouldCreatePlots = true; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(std::make_shared(beamSimulationJob), settings); simpleBeam_simulationResults.save(); @@ -58,10 +59,19 @@ void DRMSimulationModel::runUnitTests() assert(std::filesystem::exists(std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename))); Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements; Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename, simpleBeam_groundOfTruthDisplacements); - if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements)) { - std::cerr << "Error:Beam simulation produces wrong results." << std::endl; + + double error; + if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements, error)) { + std::cerr << "Error:Beam simulation produces wrong results. Error:" << error << std::endl; // return; } +#ifdef POLYSCOPE_DEFINED + beam.registerForDrawing(); + simpleBeam_simulationResults.registerForDrawing(); + polyscope::show(); + beam.unregister(); + simpleBeam_simulationResults.unregister(); +#endif // Second example of the paper VCGEdgeMesh shortSpanGrid; @@ -109,13 +119,21 @@ void DRMSimulationModel::runUnitTests() // shortSpanGridshellSimulationResults.registerForDrawing( // shortSpanGridshellSimulationJob); // polyscope::show(); - assert( - shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements)); - if (!shortSpanGridshellSimulationResults.isEqual( - shortSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Short span simulation produces wrong results." << std::endl; - return; + assert(shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements, + error)); + if (!shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements, + error)) { + std::cerr << "Error:Short span simulation produces wrong results. Error:" << error + << std::endl; + // return; } +#ifdef POLYSCOPE_DEFINED + shortSpanGrid.registerForDrawing(); + shortSpanGridshellSimulationResults.registerForDrawing(); + polyscope::show(); + shortSpanGrid.unregister(); + shortSpanGridshellSimulationResults.unregister(); +#endif // Third example of the paper VCGEdgeMesh longSpanGrid; @@ -163,6 +181,7 @@ void DRMSimulationModel::runUnitTests() << std::endl; } longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + SimulationResults longSpanGridshell_simulationResults = formFinder.executeSimulation(std::make_shared( longSpanGridshell_simulationJob), @@ -177,12 +196,21 @@ void DRMSimulationModel::runUnitTests() Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements; Eigen::read_binary(longSpan_groundOfTruthBinaryFilename, longSpanGridshell_groundOfTruthDisplacements); - assert( - longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements)); - if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Long span simulation produces wrong results." << std::endl; - return; + assert(longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements, + error)); + if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements, + error)) { + std::cerr << "Error:Long span simulation produces wrong results. Error:" << error + << std::endl; + // return; } +#ifdef POLYSCOPE_DEFINED + longSpanGrid.registerForDrawing(); + longSpanGridshell_simulationResults.registerForDrawing(); + polyscope::show(); + longSpanGrid.unregister(); + longSpanGridshell_simulationResults.unregister(); +#endif std::cout << "Form finder unit tests succesufully passed." << std::endl; @@ -1130,14 +1158,11 @@ void DRMSimulationModel::updateNodalMasses() const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const long double SkRotational_I2 - = elem.material.youngsModulus * elem.I2 - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const long double SkRotational_I3 - = elem.material.youngsModulus * elem.I3 - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const long double SkRotational_J = elem.material.youngsModulus * elem.J - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const long double SkRotational_I2 = elem.material.youngsModulus * elem.I2 + / lengthToThe3; + const long double SkRotational_I3 = elem.material.youngsModulus * elem.I3 + / lengthToThe3; + const long double SkRotational_J = elem.material.youngsModulus * elem.J / lengthToThe3; rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_J += SkRotational_J; @@ -1336,9 +1361,8 @@ void DRMSimulationModel::updateElementalFrames() } void DRMSimulationModel::applyForcedDisplacements( - const std::unordered_map nodalForcedDisplacements) + const std::unordered_map &nodalForcedDisplacements) { - const int gradualDisplacementSteps = 500; for (const std::pair vertexIndexDisplacementPair : nodalForcedDisplacements) { const VertexIndex vi = vertexIndexDisplacementPair.first; @@ -1347,9 +1371,9 @@ void DRMSimulationModel::applyForcedDisplacements( VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1), vertexDisplacement(2)); - if (mCurrentSimulationStep < gradualDisplacementSteps) { + if (mCurrentSimulationStep < mSettings.gradualForcedDisplacementSteps) { displacementVector *= mCurrentSimulationStep - / static_cast(gradualDisplacementSteps); + / static_cast(mSettings.gradualForcedDisplacementSteps); } pMesh->vert[vi].P() = node.initialLocation + displacementVector; node.displacements = Vector6d({displacementVector[0], @@ -1553,6 +1577,13 @@ void DRMSimulationModel::printCurrentState() const std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl; std::cout << "Residual forces norm: " << pMesh->totalResidualForcesNorm << std::endl; std::cout << "Kinetic energy:" << pMesh->currentTotalKineticEnergy << std::endl; + static std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); + const auto timePerNodePerIteration = std::chrono::duration_cast( + std::chrono::steady_clock::now() - begin) + .count() + * 1e-6 / (mCurrentSimulationStep * pMesh->VN()); + std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl; + std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << std::endl; } void DRMSimulationModel::printDebugInfo() const @@ -1694,22 +1725,26 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) ->addEdgeScalarQuantity("I3", I3); // Specify the callback - PolyscopeInterface::addUserCallback([&]() { - // Since options::openImGuiWindowForUserCallback == true by default, - // we can immediately start using ImGui commands to build a UI + static bool calledOnce = false; + if (!calledOnce) { + PolyscopeInterface::addUserCallback([&]() { + // Since options::openImGuiWindowForUserCallback == true by default, + // we can immediately start using ImGui commands to build a UI - ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, - // instead of full width. Must have - // matching PopItemWidth() below. + ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, + // instead of full width. Must have + // matching PopItemWidth() below. - ImGui::InputInt("Simulation drawing step", - &mSettings.drawingStep); // set a int variable - ImGui::Checkbox("Enable drawing", - &mSettings.shouldDraw); // set a int variable - ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); + ImGui::InputInt("Simulation drawing step", + &mSettings.drawingStep); // set a int variable + ImGui::Checkbox("Enable drawing", + &mSettings.shouldDraw); // set a int variable + ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); - ImGui::PopItemWidth(); - }); + ImGui::PopItemWidth(); + }); + calledOnce = true; + } if (!screenshotsFolder.empty()) { static bool firstDraw = true; @@ -1742,16 +1777,27 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &externalForce : pJob->nodalExternalForces) { totalExternalForcesNorm += externalForce.second.norm(); } - totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination - * totalExternalForcesNorm; - // if (!pJob->nodalExternalForces.empty()) { - // double externalForcesNorm = 0; - // for (const auto &externalForce : pJob->nodalExternalForces) { - // externalForcesNorm += externalForce.second.norm(); - // } - // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; - // } + if (!pJob->nodalExternalForces.empty()) { + totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination + * totalExternalForcesNorm; + } else { + totalResidualForcesNormThreshold = 1; + std::cout << "No forces setted default residual forces norm threshold" << std::endl; + } + + if (mSettings.beVerbose) { + std::cout << "totalResidualForcesNormThreshold:" << totalResidualForcesNormThreshold + << std::endl; + } + + // if (!pJob->nodalExternalForces.empty()) { + // double externalForcesNorm = 0; + // for (const auto &externalForce : pJob->nodalExternalForces) { + // externalForcesNorm += externalForce.second.norm(); + // } + // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; + // } constrainedVertices = pJob->constrainedVertices; if (!pJob->nodalForcedDisplacements.empty()) { @@ -1773,18 +1819,20 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN() << " nodes and " << pMesh->EN() << " elements." << std::endl; } + vcg::tri::UpdateBounding::Box(*pMesh); computeRigidSupports(); for (auto fixedVertex : pJob->constrainedVertices) { assert(fixedVertex.first < pMesh->VN()); } #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw ) { - PolyscopeInterface::init(); - polyscope::registerCurveNetwork( - meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges()); - // registerWorldAxes(); - } + if (mSettings.shouldDraw || mSettings.isDebugMode) { + PolyscopeInterface::init(); + polyscope::registerCurveNetwork(meshPolyscopeLabel, + pMesh->getEigenVertices(), + pMesh->getEigenEdges()); + // registerWorldAxes(); + } #endif updateElementalFrames(); @@ -1793,130 +1841,203 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces, constrainedVertices); - while (settings.maxDRMIterations == 0 || - mCurrentSimulationStep < settings.maxDRMIterations) { - // while (true) { - updateNormalDerivatives(); - updateT1Derivatives(); - updateRDerivatives(); - updateT2Derivatives(); - updateT3Derivatives(); - updateResidualForcesOnTheFly(constrainedVertices); + const size_t movingAverageSampleSize = 50; + std::queue residualForcesMovingAverageHistorySample; - // TODO: write parallel function for updating positions - // TODO: Make a single function out of updateResidualForcesOnTheFly - // updatePositionsOnTheFly - // updatePositionsOnTheFly(constrainedVertices); - updateNodalMasses(); - updateNodalAccelerations(); - updateNodalVelocities(); - updateNodalDisplacements(); - applyDisplacements(constrainedVertices); - if (!pJob->nodalForcedDisplacements.empty()) { - applyForcedDisplacements( + double residualForcesMovingAverageDerivativeMax = 0; + while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { + // while (true) { + updateNormalDerivatives(); + updateT1Derivatives(); + updateRDerivatives(); + updateT2Derivatives(); + updateT3Derivatives(); + updateResidualForcesOnTheFly(constrainedVertices); - pJob->nodalForcedDisplacements); - } - // if (!pJob->nodalForcedNormals.empty()) { - // applyForcedNormals(pJob->nodalForcedNormals); - // } - updateElementalLengths(); - // pMesh->previousTotalPotentialEnergykN = - // pMesh->currentTotalPotentialEnergykN; - // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + // TODO: write parallel function for updating positions + // TODO: Make a single function out of updateResidualForcesOnTheFly + // updatePositionsOnTheFly + // updatePositionsOnTheFly(constrainedVertices); + updateNodalMasses(); + updateNodalAccelerations(); + updateNodalVelocities(); + updateNodalDisplacements(); + applyDisplacements(constrainedVertices); + if (!pJob->nodalForcedDisplacements.empty()) { + applyForcedDisplacements( - // if (mCurrentSimulationStep != 0) { - // history.stepPulse(*pMesh); - // } - - if (std::isnan(pMesh->currentTotalKineticEnergy)) { - if (mSettings.beVerbose) { - std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; + pJob->nodalForcedDisplacements); + } + // if (!pJob->nodalForcedNormals.empty()) { + // applyForcedNormals(pJob->nodalForcedNormals); + // } + updateElementalLengths(); + mCurrentSimulationStep++; + + double sumOfDisplacementsNorm = 0; + for (size_t vi = 0; vi < pMesh->VN(); vi++) { + sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); + } + sumOfDisplacementsNorm /= pMesh->bbox.Diag(); + pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; + // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm + // + mCurrentSimulationStep + // * pMesh->residualForcesMovingAverage) + // / (1 + mCurrentSimulationStep); + pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm + / movingAverageSampleSize; + residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage); + if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) { + const double firstElementValue = residualForcesMovingAverageHistorySample.front(); + pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize; + residualForcesMovingAverageHistorySample.pop(); + + pMesh->residualForcesMovingAverageDerivativeNorm + = std::abs((residualForcesMovingAverageHistorySample.back() + - residualForcesMovingAverageHistorySample.front())) + / (movingAverageSampleSize); + residualForcesMovingAverageDerivativeMax + = std::max(pMesh->residualForcesMovingAverageDerivativeNorm, + residualForcesMovingAverageDerivativeMax); + pMesh->residualForcesMovingAverageDerivativeNorm + /= residualForcesMovingAverageDerivativeMax; + // std::cout << "Normalized derivative:" + // << residualForcesMovingAverageDerivativeNorm + // << std::endl; + } + + // pMesh->previousTotalPotentialEnergykN = + // pMesh->currentTotalPotentialEnergykN; + // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + // timePerNodePerIterationHistor.push_back(timePerNodePerIteration); + if (mSettings.beVerbose && mSettings.isDebugMode + && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + printCurrentState(); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver", + // movingAverages); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // movingAveragesDerivatives); + // draw(); + std::cout << "Mov aver deriv:" << pMesh->residualForcesMovingAverageDerivativeNorm + << std::endl; + // SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps", + // "Time/(#nodes*#iterations)", + // timePerNodePerIterationHistory); + } + + if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) { + history.stepPulse(*pMesh); + } + + if (std::isnan(pMesh->currentTotalKineticEnergy)) { + if (mSettings.beVerbose) { + std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; + } + break; } - break; - } #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) /* && + if ((mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) + || (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0)) /* && + currentSimulationStep > maxDRMIterations*/ - { - // std::string saveTo = std::filesystem::current_path() - // .append("Debugging_files") - // .append("Screenshots") - // .string(); - draw(); - // yValues = history.kineticEnergy; - // plotHandle = matplot::scatter(xPlot, yValues); - // label = "Log of Kinetic energy"; - // plotHandle->legend_string(label); + { + // std::string saveTo = std::filesystem::current_path() + // .append("Debugging_files") + // .append("Screenshots") + // .string(); + draw(); + // yValues = history.kineticEnergy; + // plotHandle = matplot::scatter(xPlot, yValues); + // label = "Log of Kinetic energy"; + // plotHandle->legend_string(label); - // shouldUseKineticEnergyThreshold = true; - } + // shouldUseKineticEnergyThreshold = true; + } #endif - if (mSettings.shouldCreatePlots && mCurrentSimulationStep % 10 == 0) { - printCurrentState(); - SimulationResultsReporter::createPlot( - "Number of Steps", "Log of Kinetic energy", history.kineticEnergy); - } - // t = t + Dt; - mCurrentSimulationStep++; - // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy - // << std::endl; - // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm - // << std::endl; - // Kinetic energy convergence - if ((mSettings.shouldUseTranslationalKineticEnergyThreshold || mCurrentSimulationStep > 500000) - && pMesh->currentTotalTranslationalKineticEnergy - < mSettings.totalTranslationalKineticEnergyThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" - << std::endl; - } -#ifdef POLYSCOPE_DEFINED - std::cout << "Warning: The kinetic energy of the system was " - " used as a convergence criterion" - << std::endl; -#endif - break; - } - // Residual forces norm convergence - if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - /*|| + if (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // movingAveragesDerivatives_norm); + SimulationResultsReporter::createPlot("Number of Steps", + "Residual Forces mov aver", + history.residualForcesMovingAverage); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces", + // history.residualForces); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Kinetic energy", + // history.kineticEnergy); + } + // t = t + Dt; + // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy + // << std::endl; + // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm + // << std::endl; + // Residual forces norm convergence + if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + && mCurrentSimulationStep > movingAverageSampleSize + && (pJob->nodalForcedDisplacements.empty() + || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps) + /*|| mesh->previousTotalPotentialEnergykN > mesh->currentTotalPotentialEnergykN*/ - /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { - if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" - << std::endl; - } - break; - // } - } - // printCurrentState(); - // for (VertexType &v : mesh->vert) { - // Node &node = mesh->nodes[v]; - // node.displacements = node.displacements - node.velocity * Dt; - // } - // applyDisplacements(constrainedVertices); - // if (!pJob->nodalForcedDisplacements.empty()) { - // applyForcedDisplacements( + /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { + // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + const bool fullfillsKineticEnergyTerminationCriterion + = mSettings.shouldUseTranslationalKineticEnergyThreshold + && pMesh->currentTotalTranslationalKineticEnergy + < mSettings.totalTranslationalKineticEnergyThreshold + && mCurrentSimulationStep > 20 && numOfDampings > 2; + const bool fullfillsResidualForcesNormTerminationCriterion + = pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold; + const bool fullfillsMovingAverageNormTerminationCriterion + = pMesh->residualForcesMovingAverage + < mSettings.residualForcesMovingAverageNormThreshold; + /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ + const bool shouldTerminate + = fullfillsKineticEnergyTerminationCriterion + // || fullfillsMovingAverageNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion; + if (shouldTerminate) { + if (mSettings.beVerbose) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + if (fullfillsResidualForcesNormTerminationCriterion) { + std::cout << "Converged using residual forces norm threshold criterion" + << std::endl; + } else if (fullfillsKineticEnergyTerminationCriterion) { + std::cout << "Warning: The kinetic energy of the system was " + " used as a convergence criterion" + << std::endl; + } else if (fullfillsMovingAverageNormTerminationCriterion) { + std::cout << "Converged using residual forces moving average derivative norm " + "threshold criterion" + << std::endl; + } + } + break; + // } + } + // printCurrentState(); + // for (VertexType &v : mesh->vert) { + // Node &node = mesh->nodes[v]; + // node.displacements = node.displacements - node.velocity * Dt; + // } + // applyDisplacements(constrainedVertices); + // if (!pJob->nodalForcedDisplacements.empty()) { + // applyForcedDisplacements( - // pJob->nodalForcedDisplacements); - // } - // updateElementalLengths(); - resetVelocities(); - Dt = Dt * mSettings.xi; - ++numOfDampings; - } - - if (mSettings.debugMessages) { - printDebugInfo(); - } + // pJob->nodalForcedDisplacements); + // } + // updateElementalLengths(); + resetVelocities(); + Dt = Dt * mSettings.xi; + ++numOfDampings; + // std::cout << "Number of dampings:" << numOfDampings << endl; + } } SimulationResults results = computeResults(pJob); @@ -1929,20 +2050,21 @@ mesh->currentTotalPotentialEnergykN*/ << std::endl; results.converged = false; } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { - results.converged = false; + std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; + results.converged = false; } else if (mSettings.beVerbose) { - auto t2 = std::chrono::high_resolution_clock::now(); - double simulationDuration = - std::chrono::duration_cast(t2 - t1).count(); - simulationDuration /= 1000; - std::cout << "Simulation converged after " << simulationDuration << "s" - << std::endl; - std::cout << "Time/(node*iteration) " - << simulationDuration / - (static_cast(mCurrentSimulationStep) * pMesh->VN()) - << "s" << std::endl; - std::cout << "Number of dampings:" << numOfDampings << endl; + auto t2 = std::chrono::high_resolution_clock::now(); + double simulationDuration = std::chrono::duration_cast(t2 - t1) + .count(); + simulationDuration /= 1000; + std::cout << "Simulation converged after " << simulationDuration << "s" << std::endl; + std::cout << "Time/(node*iteration) " + << simulationDuration / (static_cast(mCurrentSimulationStep) * pMesh->VN()) + << "s" << std::endl; + std::cout << "Number of dampings:" << numOfDampings << endl; + std::cout << "Number of steps:" << mCurrentSimulationStep << endl; + results.converged = true; } // mesh.printVertexCoordinates(mesh.VN() / 2); #ifdef POLYSCOPE_DEFINED @@ -1994,7 +2116,7 @@ mesh->currentTotalPotentialEnergykN*/ } } - if (mSettings.shouldCreatePlots) { + if (mSettings.shouldCreatePlots && results.converged) { SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 928ef08..4e790fc 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -23,28 +23,33 @@ struct DifferentiateWithRespectTo { class DRMSimulationModel { public: - struct Settings { - bool debugMessages{false}; - bool shouldDraw{false}; - bool beVerbose{false}; - bool shouldCreatePlots{false}; - int drawingStep{1}; - double totalTranslationalKineticEnergyThreshold{1e-10}; - double totalExternalForcesNormPercentageTermination{1e-3}; - double Dtini{0.1}; - double xi{0.9969}; - int maxDRMIterations{0}; - bool shouldUseTranslationalKineticEnergyThreshold{false}; - Settings() {} - }; + struct Settings + { + bool isDebugMode{false}; + int debugModeStep{10000}; + bool shouldDraw{false}; + bool beVerbose{false}; + bool shouldCreatePlots{false}; + int drawingStep{1}; + double totalTranslationalKineticEnergyThreshold{1e-10}; + double totalExternalForcesNormPercentageTermination{1e-3}; + double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; + double residualForcesMovingAverageNormThreshold{1e-8}; + double Dtini{0.1}; + double xi{0.9969}; + int maxDRMIterations{0}; + bool shouldUseTranslationalKineticEnergyThreshold{false}; + int gradualForcedDisplacementSteps{100}; + Settings() {} + }; + double totalResidualForcesNormThreshold{1}; -private: + private: Settings mSettings; double Dt{mSettings.Dtini}; bool checkedForMaximumMoment{false}; double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; - double totalResidualForcesNormThreshold{1e-3}; matplot::line_handle plotHandle; std::vector plotYValues; size_t numOfDampings{0}; @@ -80,8 +85,7 @@ private: void updateNodalDisplacements(); void applyForcedDisplacements( - const std::unordered_map - nodalForcedDisplacements); + const std::unordered_map &nodalForcedDisplacements); Vector6d computeElementTorsionalForce( const Element &element, const Vector6d &displacementDifference, @@ -208,8 +212,7 @@ private: executeSimulation(const std::shared_ptr &pJob, const Settings &settings = Settings()); - void runUnitTests(); - void setTotalResidualForcesNormThreshold(double value); + static void runUnitTests(); }; template PointType Cross(PointType p1, PointType p2) { diff --git a/edgemesh.cpp b/edgemesh.cpp index d1f4f2d..0fe03e1 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -200,7 +200,6 @@ bool VCGEdgeMesh::load(const string &plyFilename) { getEdges(eigenEdges); getVertices(eigenVertices); vcg::tri::UpdateTopology::VertexEdge(*this); - label = std::filesystem::path(plyFilename).stem().string(); const bool printDebugInfo = false; if (printDebugInfo) { @@ -208,7 +207,7 @@ bool VCGEdgeMesh::load(const string &plyFilename) { std::cout << "Mesh has " << EN() << " edges." << std::endl; } - label=std::filesystem::path(plyFilename).stem(); + label = std::filesystem::path(plyFilename).stem().string(); return true; } diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index fa881f5..7c3250a 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -126,12 +126,14 @@ struct Colors } xRange x; x.fromString(json.at(jsonXRangeKey)); + xRanges.push_back(x); xRangeIndex++; } numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); solverAccuracy = json.at(JsonKeys::SolverAccuracy); objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); + objectiveWeights.rotational = 2 - objectiveWeights.translational; return true; } @@ -197,6 +199,7 @@ struct Colors inline void updateMeshWithOptimalVariables(const std::vector &x, SimulationMesh &m) { + assert(CrossSectionType::name == RectangularBeamDimensions::name); const double E = x[0]; const double A = x[1]; const double beamWidth = std::sqrt(A); @@ -207,7 +210,7 @@ struct Colors const double I3 = x[4]; for (EdgeIndex ei = 0; ei < m.EN(); ei++) { Element &e = m.elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); + e.setDimensions(CrossSectionType(beamWidth, beamHeight)); e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); e.J = J; e.I2 = I2; @@ -494,6 +497,7 @@ struct Colors const ReducedPatternOptimization::Results &reducedPattern_optimizationResults, const shared_ptr &pTiledReducedPattern_simulationMesh) { + assert(CrossSectionType::name == RectangularBeamDimensions::name); // Set reduced parameters fron the optimization results std::unordered_map optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(), @@ -504,7 +508,7 @@ struct Colors const double A = optimalXVariables.at(ALabel); const double beamWidth = std::sqrt(A); const double beamHeight = beamWidth; - RectangularBeamDimensions elementDimensions(beamWidth, beamHeight); + CrossSectionType elementDimensions(beamWidth, beamHeight); const double poissonsRatio = 0.3; const std::string ymLabel = "E"; diff --git a/simulation_structs.hpp b/simulation_structs.hpp index fd75566..9d60766 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -9,11 +9,14 @@ struct SimulationHistory { size_t numberOfSteps{0}; std::string label; - std::vector residualForces; + std::vector residualForcesLog; std::vector kineticEnergy; std::vector potentialEnergies; std::vector redMarks; std::vector greenMarks; + std::vector residualForcesMovingAverage; + std::vector sumOfNormalizedDisplacementNorms; + // std::vector residualForcesMovingAverageDerivativesLog; void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); } @@ -22,13 +25,21 @@ struct SimulationHistory { void stepPulse(const SimulationMesh &mesh) { kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy)); // potentialEnergy.push_back(mesh.totalPotentialEnergykN); - residualForces.push_back(mesh.totalResidualForcesNorm); + residualForcesLog.push_back(std::log(mesh.totalResidualForcesNorm)); + residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage); + sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms); + // residualForcesMovingAverageDerivativesLog.push_back( + // std::log(mesh.residualForcesMovingAverageDerivativeNorm)); } - void clear() { - residualForces.clear(); - kineticEnergy.clear(); - potentialEnergies.clear(); + void clear() + { + residualForcesLog.clear(); + kineticEnergy.clear(); + potentialEnergies.clear(); + residualForcesMovingAverage.clear(); + sumOfNormalizedDisplacementNorms.clear(); + // residualForcesMovingAverageDerivativesLog.clear(); } }; @@ -399,6 +410,9 @@ public: } void unregister(const std::string &meshLabel) { + if (polyscope::getCurveNetwork(meshLabel) == nullptr) { + return; + } polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label); polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label); } @@ -478,23 +492,29 @@ struct SimulationResults void save(const std::string &outputFolder = std::string()) { - const std::filesystem::path outputPath(outputFolder); - job->save(std::filesystem::path(outputPath).append("SimulationJob").string()); + const std::filesystem::path outputFolderPath = outputFolder.empty() + ? std::filesystem::current_path() + : std::filesystem::path(outputFolder); + std::filesystem::path simulationJobOutputFolderPath + = std::filesystem::path(outputFolderPath).append("SimulationJob"); + std::filesystem::create_directories(simulationJobOutputFolderPath); + job->save(simulationJobOutputFolderPath); const std::string filename(getLabel() + "_displacements.eigenBin"); Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); - Eigen::write_binary(std::filesystem::path(outputPath).append(filename).string(), m); + Eigen::write_binary(std::filesystem::path(outputFolderPath).append(filename).string(), m); - saveDeformedModel(outputFolder); + saveDeformedModel(outputFolderPath); } // The comparison of the results happens comparing the 6-dof nodal // displacements - bool isEqual(const Eigen::MatrixXd &nodalDisplacements) + bool isEqual(const Eigen::MatrixXd &nodalDisplacements, double &error) { assert(nodalDisplacements.cols() == 6); Eigen::MatrixXd eigenDisplacements = Utilities::toEigenMatrix(this->displacements); const double errorNorm = (eigenDisplacements - nodalDisplacements).norm(); + error = errorNorm; return errorNorm < 1e-10; // return eigenDisplacements.isApprox(nodalDisplacements); } diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index 9beb5fd..422a4a8 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -38,11 +38,11 @@ struct SimulationResultsReporter { file << "\n"; } - if (!history.residualForces.empty()) { + if (!history.residualForcesLog.empty()) { file << "Residual forces" << "\n"; for (size_t step = 0; step < numberOfSteps; step++) { - file << history.residualForces[step] << "\n"; + file << history.residualForcesLog[step] << "\n"; } file << "\n"; } @@ -92,27 +92,31 @@ struct SimulationResultsReporter { } } - static void createPlot(const std::string &xLabel, const std::string &yLabel, + static void createPlot(const std::string &xLabel, + const std::string &yLabel, const std::vector &YvaluesToPlot, - const std::string &saveTo = {}) { - auto x = - matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); - createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo); - // matplot::figure(true); - // matplot::hold(matplot::on); + const std::string &saveTo = {}) + { + if (YvaluesToPlot.size() < 2) { + return; + } + auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); + createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo); + // matplot::figure(true); + // matplot::hold(matplot::on); - // ->marker_indices(history.redMarks) - // ->marker_indices(truncatedRedMarks) - // .marker_color("r") - // ->marker_size(1) - ; - // auto greenMarksY = matplot::transform( - // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; - // }); - // matplot::scatter(history.greenMarks, greenMarksY) - // ->color("green") - // .marker_size(10); - // matplot::hold(matplot::off); + // ->marker_indices(history.redMarks) + // ->marker_indices(truncatedRedMarks) + // .marker_color("r") + // ->marker_size(1) + ; + // auto greenMarksY = matplot::transform( + // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; + // }); + // matplot::scatter(history.greenMarks, greenMarksY) + // ->color("green") + // .marker_size(10); + // matplot::hold(matplot::off); } void createPlots(const SimulationHistory &history, @@ -124,27 +128,54 @@ struct SimulationResultsReporter { std::filesystem::create_directory(graphsFolder.string()); if (!history.kineticEnergy.empty()) { - createPlot("Number of Iterations", "Log of Kinetic Energy", - history.kineticEnergy, - std::filesystem::path(graphsFolder) - .append("KineticEnergy" + graphSuffix + ".png") - .string()); + createPlot("Number of Iterations", + "Log of Kinetic Energy log", + history.kineticEnergy, + std::filesystem::path(graphsFolder) + .append("KineticEnergyLog_" + graphSuffix + ".png") + .string()); } - if (!history.residualForces.empty()) { - createPlot("Number of Iterations", "Residual Forces norm", - history.residualForces, - std::filesystem::path(graphsFolder) - .append("ResidualForces" + graphSuffix + ".png") - .string()); + if (!history.residualForcesLog.empty()) { + createPlot("Number of Iterations", + "Residual Forces norm log", + history.residualForcesLog, + std::filesystem::path(graphsFolder) + .append("ResidualForcesLog_" + graphSuffix + ".png") + .string()); } if (!history.potentialEnergies.empty()) { - createPlot("Number of Iterations", "Potential energy", - history.potentialEnergies, - std::filesystem::path(graphsFolder) - .append("PotentialEnergy" + graphSuffix + ".png") - .string()); + createPlot("Number of Iterations", + "Potential energy", + history.potentialEnergies, + std::filesystem::path(graphsFolder) + .append("PotentialEnergy_" + graphSuffix + ".png") + .string()); + } + if (!history.residualForcesMovingAverage.empty()) { + createPlot("Number of Iterations", + "Residual forces moving average", + history.residualForcesMovingAverage, + std::filesystem::path(graphsFolder) + .append("ResidualForcesMovingAverage_" + graphSuffix + ".png") + .string()); + } + // if (!history.residualForcesMovingAverageDerivativesLog.empty()) { + // createPlot("Number of Iterations", + // "Residual forces moving average derivative log", + // history.residualForcesMovingAverageDerivativesLog, + // std::filesystem::path(graphsFolder) + // .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png") + // .string()); + // } + if (!history.sumOfNormalizedDisplacementNorms.empty()) { + createPlot("Number of Iterations", + "Sum of normalized displacement norms", + history.sumOfNormalizedDisplacementNorms, + std::filesystem::path(graphsFolder) + .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") + .string()); } } }; diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 95e02b6..5dea7af 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -9,7 +9,7 @@ struct Element; struct Node; using CrossSectionType = RectangularBeamDimensions; -// using CrossSectionType = CylindricalBeamDimensions; +//using CrossSectionType = CylindricalBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: @@ -44,6 +44,9 @@ public: double totalResidualForcesNorm{0}; double currentTotalPotentialEnergykN{0}; double previousTotalPotentialEnergykN{0}; + double residualForcesMovingAverageDerivativeNorm{0}; + double residualForcesMovingAverage{0}; + double sumOfNormalizedDisplacementNorms{0}; bool save(const std::string &plyFilename = std::string()); void setBeamCrossSection(const CrossSectionType &beamDimensions); void setBeamMaterial(const double &pr, const double &ym); diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index 42f242f..9aa76d9 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -464,7 +464,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 5); if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { - if (debugIsOn || savePlyFiles) { + if (debugIsOn /*|| savePlyFiles*/) { if (savePlyFiles) { auto danglingEdgesPath = std::filesystem::path(resultsPath) .append("ExceedingAngleThreshold"); diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index 9eb7e61..e1cbe47 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -435,7 +435,7 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector &nu thetaAnglesOfIncidentVectors[thetaAngleIndex] - thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]); if (absAngleDifference < angleThreshold_rad - && absAngleDifference > vcg::math::ToRad(0.01)) { + /*&& absAngleDifference > vcg::math::ToRad(0.01)*/) { // std::cout << "Found angDiff:" << absAngleDifference << std::endl; // vert[vi].C() = vcg::Color4b::Magenta; // hasAngleSmallerThan = true; @@ -817,7 +817,7 @@ std::shared_ptr PatternGeometry::tilePattern( const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, const std::vector &perSurfaceFacePatternIndices, - std::vector &tileIntoEdgesToTiledVi, + std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexToTiledPatternEdgeIndex) { perPatternIndexToTiledPatternEdgeIndex.resize(patterns.size()); diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 4b19f44..5bbc857 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -121,7 +121,7 @@ private: const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, const std::vector &perSurfaceFacePatternIndices, - std::vector &tileIntoEdgesToTiledVi, + std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexTiledPatternEdgeIndex); }; diff --git a/utilities.hpp b/utilities.hpp index 8a1be83..0289725 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -130,13 +130,13 @@ inline void parseIntegers(const std::string &str, std::vector &result) { } inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { - Eigen::MatrixXd m(v.size(), 6); + Eigen::MatrixXd m(v.size(), 6); - for (size_t vi = 0; vi < v.size(); vi++) { - const Vector6d &vec = v[vi]; - for (size_t i = 0; i < 6; i++) { - m(vi, i) = vec[i]; - } + for (size_t vi = 0; vi < v.size(); vi++) { + const Vector6d &vec = v[vi]; + for (size_t i = 0; i < 6; i++) { + m(vi, i) = vec[i]; + } } return m; @@ -220,7 +220,7 @@ inline void deinitPolyscope() return; } - polyscope::shutdown(); + polyscope::render::engine->shutdownImGui(); } inline void init() @@ -302,4 +302,24 @@ std::string to_string_with_precision(const T a_value, const int n = 2) return out.str(); } +template +size_t computeHashUnordered(const std::vector &v) +{ + size_t hash = 0; + for (const auto &el : v) { + hash += std::hash{}(el); + } + return hash; +} + +inline size_t computeHashOrdered(const std::vector &v) +{ + std::string elementsString; + for (const auto &el : v) { + elementsString += std::to_string(el); + } + + return std::hash{}(elementsString); +} + #endif // UTILITIES_H