diff --git a/beam.hpp b/beam.hpp index 79c4007..c50f12c 100644 --- a/beam.hpp +++ b/beam.hpp @@ -12,7 +12,7 @@ struct RectangularBeamDimensions { : b(width), h(height) { assert(width > 0 && height > 0); } - RectangularBeamDimensions() : b(1), h(1) {} + RectangularBeamDimensions() : b(0.002), h(0.002) {} }; struct CylindricalBeamDimensions { diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 6315f42..a485a96 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -8,18 +8,19 @@ #include #include +const bool debug = true; + void FormFinder::runUnitTests() { const std::filesystem::path groundOfTruthFolder{ "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"}; FormFinder formFinder; - formFinder.setTotalResidualForcesNormThreshold(1); // First example of the paper VCGEdgeMesh beam; // const size_t spanGridSize = 11; // mesh.createSpanGrid(spanGridSize); - beam.loadFromPly( + beam.loadPly( std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply")); std::unordered_map> fixedVertices; fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; @@ -39,8 +40,13 @@ void FormFinder::runUnitTests() { assert(CrossSectionType::name == CylindricalBeamDimensions::name); beamSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + Settings settings; + settings.Dtini = 0.1; + settings.xi = 0.9969; + settings.maxDRMIterations = 0; + settings.totalResidualForcesNormThreshold = 1; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation( - std::make_shared(beamSimulationJob), false, true); + std::make_shared(beamSimulationJob), settings); simpleBeam_simulationResults.save(); const std::string simpleBeamGroundOfTruthBinaryFilename = std::filesystem::path(groundOfTruthFolder) @@ -60,9 +66,9 @@ void FormFinder::runUnitTests() { VCGEdgeMesh shortSpanGrid; // const size_t spanGridSize = 11; // mesh.createSpanGrid(spanGridSize); - shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) - .append("shortSpanGridshell.ply") - .string()); + shortSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder) + .append("shortSpanGridshell.ply") + .string()); fixedVertices.clear(); //// Corner nodes @@ -92,7 +98,7 @@ void FormFinder::runUnitTests() { SimulationResults shortSpanGridshellSimulationResults = formFinder.executeSimulation( std::make_shared(shortSpanGridshellSimulationJob), - false, false); + settings); shortSpanGridshellSimulationResults.save(); const std::string groundOfTruthBinaryFilename = @@ -118,9 +124,9 @@ void FormFinder::runUnitTests() { // Third example of the paper VCGEdgeMesh longSpanGrid; - longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) - .append("longSpanGridshell.ply") - .string()); + longSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder) + .append("longSpanGridshell.ply") + .string()); const size_t spanGridSize = 11; fixedVertices.clear(); @@ -174,7 +180,7 @@ void FormFinder::runUnitTests() { SimulationResults longSpanGridshell_simulationResults = formFinder.executeSimulation( std::make_shared(longSpanGridshell_simulationJob), - false, false); + settings); longSpanGridshell_simulationResults.save(); const std::string longSpan_groundOfTruthBinaryFilename = @@ -200,12 +206,7 @@ void FormFinder::runUnitTests() { // polyscope::show(); } -void FormFinder::setTotalResidualForcesNormThreshold(double value) { - settings.totalResidualForcesNormThreshold = value; -} - void FormFinder::reset() { - settings.Dt = settings.Dtini; mCurrentSimulationStep = 0; history.clear(); constrainedVertices.clear(); @@ -214,9 +215,10 @@ void FormFinder::reset() { plotYValues.clear(); plotHandle.reset(); checkedForMaximumMoment = false; - shouldUseKineticEnergyThreshold = false; + mSettings.shouldUseTranslationalKineticEnergyThreshold = false; externalMomentsNorm = 0; - settings.drawingStep = 1; + mSettings.drawingStep = 1; + Dt = mSettings.Dtini; } VectorType FormFinder::computeDisplacementDifferenceDerivative( @@ -841,7 +843,7 @@ void FormFinder::updateResidualForcesOnTheFly( const double e_k = element.length - element.initialLength; const double e_kDeriv = computeDerivativeElementLength(e, dui); const double axialForce_dofi = - e_kDeriv * e_k * element.axialConstFactor; + e_kDeriv * e_k * element.rigidity.axial; // Torsional force computation const double theta1_j_deriv = @@ -852,7 +854,7 @@ void FormFinder::updateResidualForcesOnTheFly( const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; const double torsionalForce_dofi = - theta1DiffDerivative * theta1Diff * element.torsionConstFactor; + theta1DiffDerivative * theta1Diff * element.rigidity.torsional; // First bending force computation ////theta2_j derivative @@ -880,8 +882,7 @@ void FormFinder::updateResidualForcesOnTheFly( firstBendingForce_inBracketsTerm_2 + firstBendingForce_inBracketsTerm_3; const double firstBendingForce_dofi = - firstBendingForce_inBracketsTerm * - element.firstBendingConstFactor; + firstBendingForce_inBracketsTerm * element.rigidity.firstBending; // Second bending force computation ////theta2_j derivative @@ -908,7 +909,7 @@ void FormFinder::updateResidualForcesOnTheFly( secondBendingForce_inBracketsTerm_2 + secondBendingForce_inBracketsTerm_3; double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * - element.secondBendingConstFactor; + element.rigidity.secondBending; const bool shouldBreak = mCurrentSimulationStep == 118 && edgeNode.vi == 1 && dofi == 3; internalForcesContributionFromThisEdge[evi].second[dofi] = @@ -947,7 +948,7 @@ void FormFinder::updateResidualForcesOnTheFly( secondBendingForce_inBracketsTerm_3; const double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * - element.secondBendingConstFactor; + element.rigidity.secondBending; internalForcesContributionFromThisEdge[evi + 2].second[dofi] = secondBendingForce_dofi; } @@ -1162,24 +1163,24 @@ void FormFinder::updateNodalMasses() { assert(rotationalSumSk_J != 0); } mesh->nodes[v].translationalMass = - gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk; + gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk; mesh->nodes[v].rotationalMass_I2 = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2; mesh->nodes[v].rotationalMass_I3 = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3; mesh->nodes[v].rotationalMass_J = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J; - assert(std::pow(settings.Dtini, 2.0) * translationalSumSk / + assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk / mesh->nodes[v].translationalMass < 2); - assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I2 / + assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I2 / mesh->nodes[v].rotationalMass_I2 < 2); - assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I3 / + assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I3 / mesh->nodes[v].rotationalMass_I3 < 2); - assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_J / + assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_J / mesh->nodes[v].rotationalMass_J < 2); } @@ -1211,7 +1212,7 @@ void FormFinder::updateNodalVelocities() { for (VertexType &v : mesh->vert) { const VertexIndex vi = mesh->getIndex(v); Node &node = mesh->nodes[v]; - node.velocity = node.velocity + node.acceleration * settings.Dt; + node.velocity = node.velocity + node.acceleration * Dt; } updateKineticEnergy(); } @@ -1219,14 +1220,17 @@ void FormFinder::updateNodalVelocities() { void FormFinder::updateNodalDisplacements() { for (VertexType &v : mesh->vert) { Node &node = mesh->nodes[v]; - node.displacements = node.displacements + node.velocity * settings.Dt; - if (settings.beVerbose) { - std::cout << "Node " << node.vi << ":" << endl; - std::cout << node.displacements[0] << " " << node.displacements[0] << " " - << node.displacements[1] << " " << node.displacements[2] << " " - << node.displacements[3] << " " << node.displacements[4] << " " - << node.displacements[5] << std::endl; - } + node.displacements = node.displacements + node.velocity * Dt; + // if (mSettings.beVerbose) { + // std::cout << "Node " << node.vi << ":" << endl; + // std::cout << node.displacements[0] << " " << node.displacements[0] + // << " " + // << node.displacements[1] << " " << node.displacements[2] + // << " " + // << node.displacements[3] << " " << node.displacements[4] + // << " " + // << node.displacements[5] << std::endl; + // } } } @@ -1299,7 +1303,7 @@ void FormFinder::updateNodeNormal( const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0; if (!checkedForMaximumMoment && viHasMoments) { - shouldUseKineticEnergyThreshold = true; + mSettings.shouldUseTranslationalKineticEnergyThreshold = true; std::cout << "Maximum moment reached.The Kinetic energy of the system will " "be used as a convergence criterion" @@ -1340,7 +1344,7 @@ void FormFinder::applyDisplacements( updateNodeNormal(v, fixedVertices); } updateElementalFrames(); - if (settings.shouldDraw || true) { + if (mSettings.shouldDraw || true) { mesh->updateEigenEdgeAndVertices(); } } @@ -1392,6 +1396,7 @@ void FormFinder::applyForcedNormals( void FormFinder::updateKineticEnergy() { mesh->previousTotalKineticEnergy = mesh->currentTotalKineticEnergy; mesh->currentTotalKineticEnergy = 0; + mesh->currentTotalTranslationalKineticEnergy = 0; for (const VertexType &v : mesh->vert) { Node &node = mesh->nodes[v]; node.kineticEnergy = 0; @@ -1399,19 +1404,21 @@ void FormFinder::updateKineticEnergy() { const double translationalVelocityNorm = std::sqrt( std::pow(node.velocity[0], 2) + std::pow(node.velocity[1], 2) + std::pow(node.velocity[2], 2)); - node.kineticEnergy += + const double nodeTranslationalKineticEnergy = 0.5 * node.translationalMass * pow(translationalVelocityNorm, 2); + const double nodeRotationalKineticEnergy = + 0.5 * (node.rotationalMass_J * pow(node.velocity[3], 2) + + +node.rotationalMass_I3 * pow(node.velocity[4], 2) + + +node.rotationalMass_I2 * pow(node.velocity[5], 2)); + + node.kineticEnergy += + nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/; assert(node.kineticEnergy < 10000000000000000000); - // const double rotationalVelocityNorm = std::sqrt( - // std::pow(node.velocity[3], 2) + std::pow(node.velocity[4], 2) + - // std::pow(node.velocity[5], 2)); - // node.kineticEnergy += - // 0.5 * node.rotationalMass_J * pow(node.velocity[3], 2) + - // 0.5 * node.rotationalMass_I3 * pow(node.velocity[4], 2) + - // 0.5 * node.rotationalMass_I2 * pow(node.velocity[5], 2); mesh->currentTotalKineticEnergy += node.kineticEnergy; + mesh->currentTotalTranslationalKineticEnergy += + nodeTranslationalKineticEnergy; } // assert(mesh->currentTotalKineticEnergy < 100000000000000); } @@ -1460,24 +1467,24 @@ void FormFinder::updatePositionsOnTheFly( // assert(rotationalSumSk_J != 0); } mesh->nodes[v].translationalMass = - gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk; + gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk; mesh->nodes[v].rotationalMass_I2 = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2; mesh->nodes[v].rotationalMass_I3 = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3; mesh->nodes[v].rotationalMass_J = - gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J; + gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J; - // assert(std::pow(Dtini, 2.0) * translationalSumSk / + // assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk / // mesh->nodes[v].translationalMass < // 2); - // assert(std::pow(Dtini, 2.0) * rotationalSumSk_I2 / + // assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I2 / // mesh->nodes[v].rotationalMass_I2 < // 2); - // assert(std::pow(Dtini, 2.0) * rotationalSumSk_I3 / + // assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I3 / // mesh->nodes[v].rotationalMass_I3 < // 2); - // assert(std::pow(Dtini, 2.0) * rotationalSumSk_J / + // assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_J / // mesh->nodes[v].rotationalMass_J < // 2); } @@ -1503,13 +1510,13 @@ void FormFinder::updatePositionsOnTheFly( for (VertexType &v : mesh->vert) { Node &node = mesh->nodes[v]; - node.velocity = node.velocity + node.acceleration * settings.Dt; + node.velocity = node.velocity + node.acceleration * Dt; } updateKineticEnergy(); for (VertexType &v : mesh->vert) { Node &node = mesh->nodes[v]; - node.displacements = node.displacements + node.velocity * settings.Dt; + node.displacements = node.displacements + node.velocity * Dt; } for (VertexType &v : mesh->vert) { @@ -1636,45 +1643,43 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { polyscope::getCurveNetwork(meshPolyscopeLabel) ->addNodeScalarQuantity("Internal force norm", internalForcesNorm) ->setEnabled(true); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("First bending force-Translation", - internalFirstBendingTranslationForces) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("First bending force-Rotation", - internalFirstBendingRotationForces) - ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeVectorQuantity("First bending force-Translation", + // internalFirstBendingTranslationForces) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeVectorQuantity("First bending force-Rotation", + // internalFirstBendingRotationForces) + // ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("Second bending force-Translation", - internalSecondBendingTranslationForces) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("Second bending force-Rotation", - internalSecondBendingRotationForces); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->getQuantity("Second bending force-Rotation") - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeScalarQuantity("nR", nRs) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeScalarQuantity("theta3", theta3) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeScalarQuantity("theta2", theta2) - ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeVectorQuantity("Second bending force-Translation", + // internalSecondBendingTranslationForces) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeVectorQuantity("Second bending force-Rotation", + // internalSecondBendingRotationForces) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeScalarQuantity("nR", nRs) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeScalarQuantity("theta3", theta3) + // ->setEnabled(false); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeScalarQuantity("theta2", theta2) + // ->setEnabled(false); polyscope::getCurveNetwork(meshPolyscopeLabel) ->addNodeVectorQuantity("Residual force", residualForces) ->setEnabled(false); polyscope::getCurveNetwork(meshPolyscopeLabel) ->addNodeScalarQuantity("Residual force norm", residualForcesNorm) ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeScalarQuantity("Node acceleration x", accelerationX); + // polyscope::getCurveNetwork(meshPolyscopeLabel) + // ->addNodeScalarQuantity("Node acceleration x", accelerationX); // Edge quantities std::vector A(mesh->EN()); @@ -1706,9 +1711,9 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { // matching PopItemWidth() below. ImGui::InputInt("Simulation drawing step", - &settings.drawingStep); // set a int variable + &mSettings.drawingStep); // set a int variable ImGui::Checkbox("Enable drawing", - &settings.shouldDraw); // set a int variable + &mSettings.shouldDraw); // set a int variable ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); ImGui::PopItemWidth(); @@ -1733,25 +1738,35 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { } } +void FormFinder::printCurrentState() { + std::cout << "Simulation steps executed:" << mCurrentSimulationStep + << std::endl; + std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm + << std::endl; + std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy + << std::endl; +} + +void FormFinder::printDebugInfo() const { + std::cout << mesh->elements[0].rigidity.toString() << std::endl; + std::cout << "Number of dampings:" << numOfDampings << endl; +} + SimulationResults FormFinder::executeSimulation(const std::shared_ptr &pJob, - const bool &shouldDraw, const bool &beVerbose, - const bool &shouldCreatePlots, - const size_t &maxDRMIterations) { + const Settings &settings) { assert(pJob->pMesh.operator bool()); auto t1 = std::chrono::high_resolution_clock::now(); reset(); - settings.shouldDraw = shouldDraw; - settings.beVerbose = beVerbose; + mSettings = settings; - // if(!job.nodalExternalForces.empty()){ - // double externalForcesNorm=0; - // for(const auto& externalForce:job.nodalExternalForces) - // { - // externalForcesNorm+=externalForce.norm(); - // } - // setTotalResidualForcesNormThreshold(externalForcesNorm); - // } + if (!pJob->nodalExternalForces.empty()) { + double externalForcesNorm = 0; + for (const auto &externalForce : pJob->nodalExternalForces) { + externalForcesNorm += externalForce.second.norm(); + } + mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-3; + } constrainedVertices = pJob->constrainedVertices; if (!pJob->nodalForcedDisplacements.empty()) { @@ -1764,24 +1779,24 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, // if (!pJob->nodalForcedNormals.empty()) { // for (std::pair viNormalPair : // pJob->nodalForcedDisplacements) { - // assert(viNormalPair.second[2] >= 0); + // assert(viNormalPair3second[2] >= 0); // } // } mesh = std::make_unique(*pJob->pMesh); - if (beVerbose) { + if (mSettings.beVerbose) { std::cout << "Executing simulation for mesh with:" << mesh->VN() << " nodes and " << mesh->EN() << " elements." << std::endl; } computeRigidSupports(); - if (settings.shouldDraw || true) { + if (mSettings.shouldDraw || true) { if (!polyscope::state::initialized) { initPolyscope(); } polyscope::registerCurveNetwork( meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges()); - // registerWorldAxes(); + registerWorldAxes(); } for (auto fixedVertex : pJob->constrainedVertices) { assert(fixedVertex.first < mesh->VN()); @@ -1821,35 +1836,38 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, // applyForcedNormals(pJob->nodalForcedNormals); // } updateElementalLengths(); + mesh->previousTotalPotentialEnergykN = mesh->currentTotalPotentialEnergykN; + mesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; if (mCurrentSimulationStep != 0) { history.stepPulse(*mesh); } - if (mCurrentSimulationStep > 100000) { - shouldUseKineticEnergyThreshold = true; - } - if (mCurrentSimulationStep > 500000) { + if (std::isnan(mesh->currentTotalKineticEnergy)) { std::time_t now = time(0); std::tm *ltm = std::localtime(&now); - std::string dir = - "../ProblematicSimulationJobs/" + std::to_string(ltm->tm_hour) + ":" + - std::to_string(ltm->tm_min) + "_" + std::to_string(ltm->tm_mday) + - "_" + std::to_string(1 + ltm->tm_mon) + "_" + - std::to_string(1900 + ltm->tm_year); - std::filesystem::create_directories(dir); - pJob->save(dir); + std::string dir = "../ProblematicSimulationJobs/" + + std::to_string(ltm->tm_mday) + "_" + + std::to_string(1 + ltm->tm_mon) + "_" + + std::to_string(1900 + ltm->tm_year); + const std::string subDir = std::filesystem::path(dir) + .append(std::to_string(ltm->tm_hour) + + ":" + std::to_string(ltm->tm_min)) + .string(); + std::filesystem::create_directories(subDir); + pJob->save(subDir); std::cout << "Non terminating simulation found. Saved simulation job to:" << dir << std::endl; std::cout << "Exiting.." << std::endl; - FormFinder debug; - debug.executeSimulation(pJob, true, true, true); - std::terminate(); + // FormFinder debug; + // debug.executeSimulation(pJob, true, true, true); + // std::terminate(); + break; } - if (settings.shouldDraw && - mCurrentSimulationStep % settings.drawingStep == 0) /* && + if (mSettings.shouldDraw && + mCurrentSimulationStep % mSettings.drawingStep == 0) /* && currentSimulationStep > maxDRMIterations*/ { // std::string saveTo = std::filesystem::current_path() @@ -1866,63 +1884,75 @@ currentSimulationStep > maxDRMIterations*/ // shouldUseKineticEnergyThreshold = true; } - // if (mCurrentSimulationStep % 100 == 0 && mCurrentSimulationStep > - // 100000) { - // std::cout << "Current simulation step:" << mCurrentSimulationStep - // << std::endl; - // std::cout << "Residual forces norm: " << - // mesh->totalResidualForcesNorm - // << std::endl; - // std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy - // << std::endl; - // // auto yValues = history.residualForces; - // // auto xPlot = matplot::linspace(0, yValues.size(), - // yValues.size()); - // // plotHandle = matplot::scatter(xPlot, yValues); - // // std::string label = "Residual forces"; - // // plotHandle->legend_string(label); - // } + 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; - if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) { - if (/*mesh.totalPotentialEnergykN < 10 ||*/ - mesh->totalResidualForcesNorm < - settings.totalResidualForcesNormThreshold || - (shouldUseKineticEnergyThreshold && - mesh->currentTotalKineticEnergy < - settings.totalKineticEnergyThreshold)) { - if (beVerbose) { - std::cout << "Convergence after " << mCurrentSimulationStep - << " steps" << std::endl; - std::cout << "Residual force of magnitude:" - << mesh->previousTotalResidualForcesNorm << std::endl; - std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy - << std::endl; - mesh->totalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; - std::cout << "Total potential:" << mesh->totalPotentialEnergykN + // Kinetic energy convergence + if ((mSettings.shouldUseTranslationalKineticEnergyThreshold || + mCurrentSimulationStep > 100000) && + mesh->currentTotalTranslationalKineticEnergy < + mSettings.totalTranslationalKineticEnergyThreshold) { + if (mSettings.beVerbose) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + std::cout << "Total potential:" << mesh->currentTotalPotentialEnergykN + << " kNm" << std::endl; + } + std::cout << "Warning: The kinetic energy of the system was " + " used as a convergence criterion" + << std::endl; + break; + } + // Residual forces norm convergence + if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy + /*|| +mesh->previousTotalPotentialEnergykN > +mesh->currentTotalPotentialEnergykN*/ + /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { + if (mesh->totalResidualForcesNorm < + mSettings.totalResidualForcesNormThreshold) { + if (mSettings.beVerbose) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + std::cout << "Total potential:" << mesh->currentTotalPotentialEnergykN << " kNm" << std::endl; } - if (shouldUseKineticEnergyThreshold) { - std::cout << "Warning: The kinetic energy of the system was " - " used as a convergence criterion" - << std::endl; - } break; // } } + // 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(); - settings.Dt = settings.Dt * settings.xi; + Dt = Dt * mSettings.xi; + ++numOfDampings; + } + + if (debug) { + printDebugInfo(); } } if (mCurrentSimulationStep == maxDRMIterations) { std::cout << "Did not reach equilibrium before reaching the maximum number " "of DRM steps (" << maxDRMIterations << "). Breaking simulation" << std::endl; - } else if (beVerbose) { + } else if (mSettings.beVerbose) { auto t2 = std::chrono::high_resolution_clock::now(); double simulationDuration = std::chrono::duration_cast(t2 - t1).count(); @@ -1933,18 +1963,19 @@ currentSimulationStep > maxDRMIterations*/ << simulationDuration / (static_cast(mCurrentSimulationStep) * mesh->VN()) << "s" << std::endl; + std::cout << "Number of dampings:" << numOfDampings << endl; } // mesh.printVertexCoordinates(mesh.VN() / 2); - if (settings.shouldDraw) { + SimulationResults results = computeResults(pJob); + if (mSettings.shouldCreatePlots) { + SimulationResultsReporter reporter; + reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); + } + if (mSettings.shouldDraw) { draw(); } if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) { polyscope::removeCurveNetwork(meshPolyscopeLabel); } - SimulationResults results = computeResults(pJob); - if (shouldCreatePlots) { - SimulationResultsReporter reporter; - reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); - } return results; } diff --git a/beamformfinder.hpp b/beamformfinder.hpp index cc8aeb3..70ae184 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -21,27 +21,32 @@ struct DifferentiateWithRespectTo { const VertexType &v; const DoFType &dofi; }; + class FormFinder { +public: struct Settings { bool shouldDraw{false}; + bool beVerbose{false}; + bool shouldCreatePlots{false}; int drawingStep{1}; - const double totalKineticEnergyThreshold{1e-9}; + double totalTranslationalKineticEnergyThreshold{1e-7}; double totalResidualForcesNormThreshold{1e-3}; - const double Dtini{0.000001}; - double Dt{Dtini}; - const double xi{0.9969}; - bool beVerbose; + double Dtini{0.1}; + double xi{0.9969}; + int maxDRMIterations{0}; + bool shouldUseTranslationalKineticEnergyThreshold{false}; + Settings() {} }; - Settings settings; - private: - bool shouldUseKineticEnergyThreshold{false}; + Settings mSettings; + double Dt{mSettings.Dtini}; bool checkedForMaximumMoment{false}; double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; matplot::line_handle plotHandle; std::vector plotYValues; + size_t numOfDampings{0}; const std::string meshPolyscopeLabel{"Simulation mesh"}; std::unique_ptr mesh; @@ -189,13 +194,15 @@ private: void applyForcedNormals( const std::unordered_map nodalForcedRotations); + void printCurrentState(); + + void printDebugInfo() const; + public: FormFinder(); - SimulationResults executeSimulation( - const std::shared_ptr &pJob, - const bool &shouldDraw = false, const bool &beVerbose = false, - const bool &createPlots = false, - const size_t &maxDRMIterations = FormFinder::maxDRMIterations); + SimulationResults + executeSimulation(const std::shared_ptr &pJob, + const Settings &settings = Settings()); inline static const size_t maxDRMIterations{0}; static void runUnitTests(); diff --git a/edgemesh.cpp b/edgemesh.cpp index 2a51a68..c925f70 100644 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -174,7 +174,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth, return true; } -bool VCGEdgeMesh::loadFromPly(const std::string plyFilename) { +bool VCGEdgeMesh::loadPly(const std::string plyFilename) { std::string usedPath = plyFilename; if (std::filesystem::path(plyFilename).is_relative()) { @@ -352,5 +352,4 @@ void VCGEdgeMesh::registerForDrawing() const { const double drawingRadius = 0.0007; polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges()) ->setRadius(drawingRadius, false); - std::cout << "Registered:" << label << std::endl; } diff --git a/edgemesh.hpp b/edgemesh.hpp index 7f176a7..613deed 100644 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -66,7 +66,7 @@ public: bool loadUsingNanoply(const std::string &plyFilename); - bool loadFromPly(const std::string plyFilename); + virtual bool loadPly(const std::string plyFilename); bool savePly(const std::string plyFilename); diff --git a/elementalmesh.cpp b/elementalmesh.cpp index 16bc99a..beacd15 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -141,7 +141,7 @@ void SimulationMesh::reset() { const vcg::Segment3 s(p0, p1); element.initialLength = s.Length(); element.length = element.initialLength; - element.updateConstFactors(); + element.updateRigidity(); } for (const VertexType &v : vert) { @@ -194,7 +194,7 @@ void SimulationMesh::initializeElements() { element.initialLength = s.Length(); element.length = element.initialLength; // Initialize const factors - element.updateConstFactors(); + element.updateRigidity(); element.derivativeT1.resize( 2, std::vector(6, VectorType(0, 0, 0))); element.derivativeT2.resize( @@ -233,14 +233,14 @@ void SimulationMesh::setBeamCrossSection( for (size_t ei = 0; ei < EN(); ei++) { elements[ei].properties.dimensions = beamDimensions; elements[ei].properties.computeDimensionsProperties(beamDimensions); - elements[ei].updateConstFactors(); + elements[ei].updateRigidity(); } } void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { for (size_t ei = 0; ei < EN(); ei++) { elements[ei].properties.setMaterial(ElementMaterial{pr, ym}); - elements[ei].updateConstFactors(); + elements[ei].updateRigidity(); } } @@ -316,7 +316,7 @@ bool SimulationMesh::loadPly(const string &plyFilename) { if (!handleBeamProperties._handle->data.empty()) { for (size_t ei = 0; ei < EN(); ei++) { elements[ei].properties = Element::Properties(handleBeamProperties[ei]); - elements[ei].updateConstFactors(); + elements[ei].updateRigidity(); } } @@ -469,9 +469,9 @@ std::array Element::Properties::toArray() const { return std::array({E, G, A, I2, I3, J}); } -void Element::updateConstFactors() { - axialConstFactor = properties.E * properties.A / initialLength; - torsionConstFactor = properties.G * properties.J / initialLength; - firstBendingConstFactor = 2 * properties.E * properties.I2 / initialLength; - secondBendingConstFactor = 2 * properties.E * properties.I3 / initialLength; +void Element::updateRigidity() { + rigidity.axial = properties.E * properties.A / initialLength; + rigidity.torsional = properties.G * properties.J / initialLength; + rigidity.firstBending = 2 * properties.E * properties.I2 / initialLength; + rigidity.secondBending = 2 * properties.E * properties.I3 / initialLength; } diff --git a/elementalmesh.hpp b/elementalmesh.hpp index 7b66a16..b125642 100644 --- a/elementalmesh.hpp +++ b/elementalmesh.hpp @@ -33,15 +33,17 @@ public: std::vector getIncidentElements(const VCGEdgeMesh::VertexType &v); - bool loadPly(const string &plyFilename); + virtual bool loadPly(const string &plyFilename); std::vector getBeamDimensions(); std::vector getBeamMaterial(); double previousTotalKineticEnergy{0}; double previousTotalResidualForcesNorm{0}; double currentTotalKineticEnergy{0}; + double currentTotalTranslationalKineticEnergy{0}; double totalResidualForcesNorm{0}; - double totalPotentialEnergykN{0}; + double currentTotalPotentialEnergykN{0}; + double previousTotalPotentialEnergykN{0}; bool savePly(const std::string &plyFilename); void setBeamCrossSection(const CrossSectionType &beamDimensions); void setBeamMaterial(const double &pr, const double &ym); @@ -50,6 +52,7 @@ public: }; struct Element { + struct Properties { CrossSectionType dimensions; ElementMaterial material; @@ -81,17 +84,28 @@ struct Element { VectorType t3; }; - void updateConstFactors(); - EdgeIndex ei; double length{0}; Properties properties; double initialLength; LocalFrame frame; - double axialConstFactor; - double torsionConstFactor; - double firstBendingConstFactor; - double secondBendingConstFactor; + + struct Rigidity { + double axial; + double torsional; + double firstBending; + double secondBending; + std::string toString() const { + return std::string("Rigidity:") + std::string("\nAxial=") + + std::to_string(axial) + std::string("\nTorsional=") + + std::to_string(torsional) + std::string("\nFirstBending=") + + std::to_string(firstBending) + std::string("\nSecondBending=") + + std::to_string(secondBending); + } + }; + Rigidity rigidity; + void updateRigidity(); + VectorType f1_j; VectorType f1_jplus1; VectorType f2_j; diff --git a/flatpattern.cpp b/flatpattern.cpp index f237c0d..a942cf3 100644 --- a/flatpattern.cpp +++ b/flatpattern.cpp @@ -5,8 +5,16 @@ FlatPattern::FlatPattern() {} FlatPattern::FlatPattern(const std::string &filename, bool addNormalsIfAbsent) { - assert(std::filesystem::exists(std::filesystem::path(filename))); - loadFromPly(filename); + if (!std::filesystem::exists(std::filesystem::path(filename))) { + assert(false); + std::cerr << "No flat pattern with name " << filename << std::endl; + return; + } + if (!loadPly(filename)) { + assert(false); + std::cerr << "File could not be loaded " << filename << std::endl; + return; + } if (addNormalsIfAbsent) { bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; if (normalsAreAbsent) { diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index b7be5bc..a807f34 100644 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -85,15 +85,17 @@ struct SimulationResultsReporter { const std::string &saveTo = {}) { matplot::xlabel(xLabel); matplot::ylabel(yLabel); - matplot::figure(true); + // matplot::figure(true); // matplot::hold(matplot::on); matplot::grid(matplot::on); auto x = matplot::linspace(0, YvaluesToPlot.size(), YvaluesToPlot.size()); + matplot::scatter(x, YvaluesToPlot) // ->marker_indices(history.redMarks) // ->marker_indices(truncatedRedMarks) // .marker_color("r") - ->marker_size(1); + // ->marker_size(1) + ; // auto greenMarksY = matplot::transform( // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; // });