diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index eb7d742..4dd80cb 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -640,6 +640,9 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, ? node.derivativeOfNormal[dui.dofi] : VectorType(0, 0, 0); derivativeTheta3_dofi = derivativeG1 * n + g1 * derivativeNormal; + if (std::isnan(derivativeTheta3_dofi)) { + std::cerr << "nan derivative theta3 rigid" << std::endl; + } return derivativeTheta3_dofi; } EdgeType &refElem = *node.referenceElement; @@ -718,6 +721,9 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, derivativeTheta3_dofi = derivativeF3 * n + f3 * derivativeOfNormal; } } + if (std::isnan(derivativeTheta3_dofi)) { + std::cerr << "nan derivative theta3" << std::endl; + } return derivativeTheta3_dofi; } @@ -794,9 +800,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( pMesh->EN(), std::vector>(4, {-1, Vector6d()})); // omp_lock_t writelock; // omp_init_lock(&writelock); -#ifdef ENABLE_OPENMP -#pragma omp parallel for schedule(static) num_threads(6) -#endif + //#ifdef ENABLE_OPENMP + //#pragma omp parallel for schedule(static) num_threads(8) + //#endif for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; const SimulationMesh::VertexType &ev_j = *e.cV(0); @@ -846,6 +852,12 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const double e_kDeriv = computeDerivativeElementLength(e, dui); const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial; +#ifdef POLYSCOPE_DEFINED + if (std::isnan(axialForce_dofi)) { + std::cerr << "nan in axial" << evi << std::endl; + } +#endif + // Torsional force computation const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi); const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi); @@ -853,6 +865,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff * element.rigidity.torsional; +#ifdef POLYSCOPE_DEFINED + if (std::isnan(torsionalForce_dofi)) { + std::cerr << "nan in torsional" << evi << std::endl; + } +#endif // First bending force computation ////theta2_j derivative @@ -970,6 +987,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( Node::Forces &force = pMesh->nodes[vi].force; force.internal = force.internal + internalForcePair.second; force.residual = force.residual + (internalForcePair.second * -1); +#ifdef POLYSCOPE_DEFINED + if (std::isnan(internalForcePair.second.norm())) { + std::cerr << "nan on edge" << ei << std::endl; + } +#endif } } double totalResidualForcesNorm = 0; @@ -982,6 +1004,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( // const double residualForceNorm = nodeResidualForce.getTranslation().norm(); const bool shouldTrigger = std::isnan(residualForceNorm); totalResidualForcesNorm += residualForceNorm; +#ifdef POLYSCOPE_DEFINED + if (std::isnan(force.residual.norm())) { + std::cout << "residual " << vi << ":" << force.residual.toString() << std::endl; + } +#endif } pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; pMesh->totalResidualForcesNorm = totalResidualForcesNorm; @@ -1218,6 +1245,11 @@ void DRMSimulationModel::updateNodalAccelerations() node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2; } } +#ifdef POLYSCOPE_DEFINED + if (std::isnan(node.acceleration.norm())) { + std::cout << "acceleration " << vi << ":" << node.acceleration.toString() << std::endl; + } +#endif // if (vi == 10) { // std::cout << "Acceleration:" << node.acceleration[0] << " " << node.acceleration[1] // << " " << node.acceleration[2] << std::endl; @@ -1234,10 +1266,12 @@ void DRMSimulationModel::updateNodalVelocities() for (VertexType &v : pMesh->vert) { const VertexIndex vi = pMesh->getIndex(v); Node &node = pMesh->nodes[v]; - node.velocity = node.velocity + node.acceleration * Dt; +#ifdef POLYSCOPE_DEFINED if (std::isnan(node.velocity.norm())) { - std::cout << "Velocity:" << node.velocity.toString() << std::endl; + std::cout << "Velocity " << vi << ":" << node.velocity.toString() << std::endl; } +#endif + node.velocity = node.velocity + node.acceleration * Dt; } updateKineticEnergy(); } @@ -1609,6 +1643,11 @@ void DRMSimulationModel::printCurrentState() const { std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl; std::cout << "Residual forces norm: " << pMesh->totalResidualForcesNorm << std::endl; + std::cout << "Average Residual forces norm/extForcesNorm: " + << pMesh->totalResidualForcesNorm / pMesh->VN() / pMesh->totalExternalForcesNorm + << std::endl; + std::cout << "Moving average residual forces:" << pMesh->residualForcesMovingAverage + << 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( @@ -1758,8 +1797,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J); polyscope::getCurveNetwork(meshPolyscopeLabel) ->addEdgeScalarQuantity("I2", I2); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addEdgeScalarQuantity("I3", I3); + polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3); // Specify the callback static bool calledOnce = false; @@ -1798,9 +1836,8 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) .append(std::to_string(mCurrentSimulationStep) + ".png") .string(), false); - } else { - polyscope::show(); } + polyscope::show(); } #endif @@ -1873,6 +1910,9 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue Eigen::Quaternion q_normal; q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); Eigen::Quaternion q_nr = q_f1.inverse() * q_normal.inverse() * q; + q_nr.w() = q_nr.w() >= 1 ? 1 : q_nr.w(); + q_nr.w() = q_nr.w() <= -1 ? -1 : q_nr.w(); + const double nr_0To2pi_pos = 2 * std::acos(q_nr.w()); // const double nr_0To2pi_groundOfTruth = 2 * std::acos(solutionGuess.debug_q_nr[vi].w()); const double nr_0To2pi = nr_0To2pi_pos <= M_PI ? nr_0To2pi_pos : (nr_0To2pi_pos - 2 * M_PI); @@ -1880,6 +1920,16 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue q_nr.y() * sin(nr_0To2pi_pos / 2), q_nr.z() * sin(nr_0To2pi_pos / 2)); /*deformedNormal_debug =*/deformedNormal_debug.normalize(); +#ifdef POLYSCOPE_DEFINED + if (std::isnan(deformedNormal_debug.norm())) { + std::cerr << "nr_0To2pi_pos:" << nr_0To2pi_pos << std::endl; + std::cerr << "q_nrx:" << q_nr.x() << std::endl; + std::cerr << "q_nry:" << q_nr.y() << std::endl; + std::cerr << "q_nrz:" << q_nr.z() << std::endl; + std::cerr << "q_nrw:" << q_nr.w() << std::endl; + std::cerr << "nan deformedNormal in guess" << std::endl; + } +#endif const double nr = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? nr_0To2pi : -nr_0To2pi; if (!pJob->constrainedVertices.contains(vi) || !pJob->constrainedVertices.at(vi).contains(5)) { @@ -2006,6 +2056,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrgetEigenVertices(), pMesh->getEigenEdges()); + polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, + pMesh->getEigenVertices(), + pMesh->getEigenEdges()); // registerWorldAxes(); } #endif @@ -2039,15 +2092,11 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrtotalExternalForcesNorm = totalExternalForcesNorm; updateNodalExternalForces(nodalExternalForces, constrainedVertices); - double averageExternalForcesNorm = - // sumOfExternalForces.norm() / pMesh->VN(); - totalExternalForcesNorm / pMesh->VN(); if (!nodalExternalForces.empty()) { mSettings.totalResidualForcesNormThreshold - = std::min(settings.totalExternalForcesNormPercentageTermination - * totalExternalForcesNorm, - 1e-3); + = settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm; } else { mSettings.totalResidualForcesNormThreshold = 1e-3; std::cout << "No forces setted default residual forces norm threshold" << std::endl; @@ -2055,9 +2104,13 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr residualForcesMovingAverageHistorySample; + const size_t movingAverageSampleSize = 200; + std::queue residualForcesMovingAverageHistorySample; std::vector percentageOfConvergence; // double residualForcesMovingAverageDerivativeMax = 0; @@ -2107,37 +2160,40 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN(); vi++) { // sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); // } // sumOfDisplacementsNorm /= pMesh->bbox.Diag(); // pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; + + //moving average // // 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->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->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; @@ -2151,19 +2207,20 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr(t2 - beginTime).count(); // std::cout << "Execution time(min):" << executionTime_min << std::endl; - if (mSettings.useAverageResidualForcesNorm) { - const double percOfConv = averageExternalForcesNorm - / pMesh->averageResidualForcesNorm; - std::cout << "Percentage of target:" << percOfConv << "%" << std::endl; - } else { - std::cout << "Percentage of target:" - << 100 * mSettings.totalResidualForcesNormThreshold - / pMesh->totalResidualForcesNorm + if (mSettings.useAverage) { + std::cout << "Percentage of target (average):" + << 100 * mSettings.averageResidualForcesCriterionThreshold + * totalExternalForcesNorm + / (pMesh->totalResidualForcesNorm / pMesh->VN()) << "%" << std::endl; } - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver", - // movingAverages); + std::cout << "Percentage of target:" + << 100 * mSettings.totalExternalForcesNormPercentageTermination + * totalExternalForcesNorm / pMesh->totalResidualForcesNorm + << "%" << std::endl; + SimulationResultsReporter::createPlot("Number of Steps", + "Residual Forces mov aver", + history.residualForcesMovingAverage); // SimulationResultsReporter::createPlot("Number of Steps", // "Residual Forces mov aver deriv", // movingAveragesDerivatives); @@ -2184,14 +2241,14 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrcurrentTotalPotentialEnergykN*/ && pMesh->currentTotalTranslationalKineticEnergy < mSettings.totalTranslationalKineticEnergyThreshold && mCurrentSimulationStep > 20 && numOfDampings > 0; - const bool fullfillsAverageResidualForcesNormCriterion - = mSettings.useAverageResidualForcesNorm - && pMesh->averageResidualForcesNorm / averageExternalForcesNorm < 1e-2; const bool fullfillsResidualForcesNormTerminationCriterion - = mSettings.useResidualForcesNorm - && pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; + = pMesh->totalResidualForcesNorm / totalExternalForcesNorm + < mSettings.totalExternalForcesNormPercentageTermination; + const bool fullfillsAverageResidualForcesNormTerminationCriterion + = mSettings.useAverage + && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm + < mSettings.averageResidualForcesCriterionThreshold; const bool fullfillsMovingAverageNormTerminationCriterion = pMesh->residualForcesMovingAverage < mSettings.residualForcesMovingAverageNormThreshold; @@ -2273,10 +2331,10 @@ mesh->currentTotalPotentialEnergykN*/ = fullfillsKineticEnergyTerminationCriterion // || fullfillsMovingAverageNormTerminationCriterion // || fullfillsMovingAverageDerivativeNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion - || fullfillsAverageResidualForcesNormCriterion; + || fullfillsAverageResidualForcesNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion; if (shouldTerminate) { - if (mSettings.beVerbose && !mSettings.isDebugMode) { + if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) { std::cout << "Simulation converged." << std::endl; printCurrentState(); if (fullfillsResidualForcesNormTerminationCriterion) { @@ -2308,7 +2366,6 @@ mesh->currentTotalPotentialEnergykN*/ } updateNodalExternalForces(nodalExternalForces, constrainedVertices); - averageExternalForcesNorm = totalExternalForcesNorm / pMesh->VN(); if (!nodalExternalForces.empty()) { mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm; } @@ -2320,25 +2377,25 @@ mesh->currentTotalPotentialEnergykN*/ // } } - const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); - for (VertexType &v : pMesh->vert) { - Node &node = pMesh->nodes[v]; - Vector6d stepDisplacement = node.velocity * Dt; - if (shouldCapDisplacements - && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { - stepDisplacement = stepDisplacement - * (*mSettings.displacementCap - / stepDisplacement.getTranslation().norm()); - } - node.displacements = node.displacements - stepDisplacement; - } - applyDisplacements(constrainedVertices); - if (!pJob->nodalForcedDisplacements.empty()) { - applyForcedDisplacements( + // const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); + // for (VertexType &v : pMesh->vert) { + // Node &node = pMesh->nodes[v]; + // Vector6d stepDisplacement = node.velocity * Dt; + // if (shouldCapDisplacements + // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { + // stepDisplacement = stepDisplacement + // * (*mSettings.displacementCap + // / stepDisplacement.getTranslation().norm()); + // } + // node.displacements = node.displacements - stepDisplacement; + // } + // applyDisplacements(constrainedVertices); + // if (!pJob->nodalForcedDisplacements.empty()) { + // applyForcedDisplacements( - pJob->nodalForcedDisplacements); - } - updateElementalLengths(); + // pJob->nodalForcedDisplacements); + // } + // updateElementalLengths(); // const double triggerPercentage = 0.01; // const double xi_min = 0.55; @@ -2491,5 +2548,12 @@ mesh->currentTotalPotentialEnergykN*/ reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); } - return results; +#ifdef POLYSCOPE_DEFINED + if (mSettings.shouldDraw || mSettings.isDebugMode) { + polyscope::removeCurveNetwork(meshPolyscopeLabel); + polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); + } +#endif + + return results; } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 96422bf..6b98674 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -38,14 +38,14 @@ public: double xi{0.9969}; int maxDRMIterations{0}; bool shouldUseTranslationalKineticEnergyThreshold{false}; - int gradualForcedDisplacementSteps{100}; + int gradualForcedDisplacementSteps{50}; int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; std::optional displacementCap; double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; - bool useResidualForcesNorm{true}; - bool useAverageResidualForcesNorm{false}; + bool useAverage{false}; + double averageResidualForcesCriterionThreshold{1e-3}; Settings() {} }; diff --git a/simulation_structs.hpp b/simulation_structs.hpp index f27e8be..9d3b1df 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -53,8 +53,8 @@ struct SimulationHistory { kineticEnergy.push_back(std::log10(mesh.currentTotalKineticEnergy)); // potentialEnergy.push_back(mesh.totalPotentialEnergykN); logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm)); - residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage); - sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms); + residualForcesMovingAverage.push_back(std::log10(mesh.residualForcesMovingAverage)); + sumOfNormalizedDisplacementNorms.push_back(std::log10(mesh.sumOfNormalizedDisplacementNorms)); // residualForcesMovingAverageDerivativesLog.push_back( // std::log(mesh.residualForcesMovingAverageDerivativeNorm)); } diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 150f8d3..af36c7c 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -43,6 +43,7 @@ public: double currentTotalKineticEnergy{0}; double currentTotalTranslationalKineticEnergy{0}; double totalResidualForcesNorm{0}; + double totalExternalForcesNorm{0}; double averageResidualForcesNorm{0}; double currentTotalPotentialEnergykN{0}; double previousTotalPotentialEnergykN{0};