From e9707e2cfb60d3c52e2de0f2449a1e735dd78965 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Fri, 6 May 2022 16:16:43 +0300 Subject: [PATCH] Refactoring --- drmsimulationmodel.cpp | 690 ++++++++++++++++++++++++----------------- drmsimulationmodel.hpp | 22 +- 2 files changed, 417 insertions(+), 295 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index ae36bd8..4a39a3b 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -263,10 +263,10 @@ void DRMSimulationModel::reset(const std::shared_ptr &pJob, const #ifdef POLYSCOPE_DEFINED if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { PolyscopeInterface::init(); - polyscope::registerCurveNetwork(meshPolyscopeLabel, + polyscope::registerCurveNetwork(meshPolyscopeLabel + "_" + pMesh->getLabel(), pMesh->getEigenVertices(), pMesh->getEigenEdges()); - polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, + polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(), pMesh->getEigenVertices(), pMesh->getEigenEdges()) ->setRadius(0.002); @@ -274,9 +274,9 @@ void DRMSimulationModel::reset(const std::shared_ptr &pJob, const } #endif - if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { - shouldApplyInitialDistortion = true; - } + // if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { + // shouldApplyInitialDistortion = true; + // } updateElementalFrames(); } @@ -584,8 +584,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v) const double &nR = node.nR; // TODO: This makes the function non-const. // Should be refactored in the future - double theta3; - const bool shouldBreak = mCurrentSimulationStep == 12970; + // const bool shouldBreak = mCurrentSimulationStep == 12970; if (&e == node.referenceElement) { // Use nR as theta3 only for the first star edge return nR; @@ -631,7 +630,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v) element.cosRotationAngle_jplus1 = cosRotationAngle; element.sinRotationAngle_jplus1 = sinRotationAngle; } - theta3 = f3.dot(n); + const double theta3 = f3.dot(n); return theta3; } @@ -845,9 +844,18 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const std::unordered_map> &fixedVertices) { // std::vector internalForcesParallel(mesh->vert.size()); + std::for_each(pMesh->nodes._handle->data.begin(), + pMesh->nodes._handle->data.end(), + [](Node &node) { + node.force.internalAxial = 0; + node.force.internalTorsion = 0; + node.force.internalFirstBending = 0; + node.force.internalSecondBending = 0; + }); std::vector>> internalForcesContributionsFromEachEdge( pMesh->EN(), std::vector>(4, {-1, Vector6d()})); + // omp_lock_t writelock; // omp_init_lock(&writelock); //#ifdef ENABLE_OPENMP @@ -909,6 +917,7 @@ void DRMSimulationModel::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.rigidity.axial; + pMesh->nodes[edgeNode.vi].force.internalAxial[dofi] += axialForce_dofi; #ifdef POLYSCOPE_DEFINED if (std::isnan(axialForce_dofi)) { @@ -923,9 +932,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff * element.rigidity.torsional; + pMesh->nodes[edgeNode.vi].force.internalTorsion[dofi] += torsionalForce_dofi; #ifdef POLYSCOPE_DEFINED if (std::isnan(torsionalForce_dofi)) { std::cerr << "nan in torsional" << evi << std::endl; + std::terminate(); } #endif @@ -954,6 +965,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( + firstBendingForce_inBracketsTerm_3; const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm * element.rigidity.firstBending; + pMesh->nodes[edgeNode.vi].force.internalFirstBending[dofi] + += firstBendingForce_dofi; // Second bending force computation ////theta2_j derivative @@ -982,6 +995,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( + secondBendingForce_inBracketsTerm_3; double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * element.rigidity.secondBending; + pMesh->nodes[edgeNode.vi].force.internalSecondBending[dofi] + += secondBendingForce_dofi; internalForcesContributionFromThisEdge[evi].second[dofi] = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi + torsionalForce_dofi; @@ -1021,6 +1036,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( * element.rigidity.secondBending; internalForcesContributionFromThisEdge[evi + 2].second[dofi] = secondBendingForce_dofi; + pMesh->nodes[refElemOtherVertexNode.vi].force.internalSecondBending[dofi] + += secondBendingForce_dofi; } } } @@ -1028,18 +1045,14 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge; }); - //#pragma omp parallel for schedule(static) num_threads(8) + std::for_each(pMesh->nodes._handle->data.begin(), + pMesh->nodes._handle->data.end(), + [](Node &node) { + Node::Forces &force = node.force; + force.residual = force.external; + force.internal = 0; + }); - for (size_t vi = 0; vi < pMesh->VN(); vi++) { - Node::Forces &force = pMesh->nodes[vi].force; - // Vector6d ext = force.external; - // if (mCurrentSimulationStep <= 100) { - // ext[3] = 0; - // ext[4] = 0; - // } - force.residual = force.external; - force.internal = 0; - } for (size_t ei = 0; ei < pMesh->EN(); ei++) { for (int i = 0; i < 4; i++) { std::pair internalForcePair @@ -1054,10 +1067,23 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( #ifdef POLYSCOPE_DEFINED if (std::isnan(internalForcePair.second.norm())) { std::cerr << "nan on edge" << ei << std::endl; + std::terminate(); } #endif } } +#ifdef POLYSCOPE_DEFINED + std::for_each(pMesh->nodes._handle->data.begin(), + pMesh->nodes._handle->data.end(), + [](Node &node) { + Node::Forces &force = node.force; + const Vector6d debug_residual = force.external - force.internalAxial + - force.internalTorsion + - force.internalFirstBending + - force.internalSecondBending; + assert((force.residual - debug_residual).norm() < 1e-5); + }); +#endif double totalResidualForcesNorm = 0; Vector6d sumOfResidualForces(0); for (size_t vi = 0; vi < pMesh->VN(); vi++) { @@ -1070,7 +1096,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( #ifdef POLYSCOPE_DEFINED if (std::isnan(force.residual.norm())) { std::cout << "residual " << vi << ":" << force.residual.toString() << std::endl; + std::terminate(); } +// std::cout << "residual " << vi << ":" << force.residual.norm() << std::endl; #endif } pMesh->totalResidualForcesNorm = totalResidualForcesNorm; @@ -1448,6 +1476,7 @@ void DRMSimulationModel::updateT2Derivatives() { for (const EdgeType &e : pMesh->edge) { for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + const auto ei = pMesh->getIndex(e); DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; pMesh->elements[e].derivativeT2_j[dofi] = computeDerivativeT2(e, dui_v0); DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; @@ -1830,8 +1859,11 @@ void DRMSimulationModel::applyForcedNormals( void DRMSimulationModel::updateKineticEnergy() { pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy; + pMesh->previousTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy; + pMesh->previousTotalRotationalKineticEnergy = pMesh->currentTotalRotationalKineticEnergy; pMesh->currentTotalKineticEnergy = 0; pMesh->currentTotalTranslationalKineticEnergy = 0; + pMesh->currentTotalRotationalKineticEnergy = 0; for (const VertexType &v : pMesh->vert) { Node &node = pMesh->nodes[v]; node.kineticEnergy = 0; @@ -1853,6 +1885,7 @@ void DRMSimulationModel::updateKineticEnergy() pMesh->currentTotalKineticEnergy += node.kineticEnergy; pMesh->currentTotalTranslationalKineticEnergy += nodeTranslationalKineticEnergy; + pMesh->currentTotalRotationalKineticEnergy += nodeRotationalKineticEnergy; } // assert(mesh->currentTotalKineticEnergy < 100000000000000); } @@ -2009,7 +2042,7 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptrcurrentTotalKineticEnergy)) { @@ -2018,46 +2051,46 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptrconstrainedVertices); - std::vector perVertexInternalForces_axial = [=]() { - std::vector axialInternalForces; - axialInternalForces.reserve(results.perVertexInternalForces.size()); - for (const std::array &internalForces : results.perVertexInternalForces) { - axialInternalForces.push_back(internalForces[0].norm()); - } - return axialInternalForces; - }(); - const auto pCurvNet = pMesh->registerForDrawing(); - pCurvNet->addNodeScalarQuantity("axial forces", perVertexInternalForces_axial); - std::vector perVertexInternalForces_torsion = [=]() { - std::vector axialInternalForces; - axialInternalForces.reserve(results.perVertexInternalForces.size()); - for (const std::array &internalForces : results.perVertexInternalForces) { - axialInternalForces.push_back(internalForces[1].norm()); - } - return axialInternalForces; - }(); - pCurvNet->addNodeScalarQuantity("torsional forces", perVertexInternalForces_torsion); - std::vector perVertexInternalForces_firstBending = [=]() { - std::vector axialInternalForces; - axialInternalForces.reserve(results.perVertexInternalForces.size()); - for (const std::array &internalForces : results.perVertexInternalForces) { - axialInternalForces.push_back(internalForces[2].norm()); - } - return axialInternalForces; - }(); - pCurvNet->addNodeScalarQuantity("first bending forces", - perVertexInternalForces_firstBending); - std::vector perVertexInternalForces_secondBending = [=]() { - std::vector axialInternalForces; - axialInternalForces.reserve(results.perVertexInternalForces.size()); - for (const std::array &internalForces : results.perVertexInternalForces) { - axialInternalForces.push_back(internalForces[3].norm()); - } - return axialInternalForces; - }(); - pCurvNet->addNodeScalarQuantity("second bending forces", - perVertexInternalForces_secondBending); - polyscope::show(); +// std::vector perVertexInternalForces_axial = [=]() { +// std::vector axialInternalForces; +// axialInternalForces.reserve(results.perVertexInternalForces.size()); +// for (const std::array &internalForces : results.perVertexInternalForces) { +// axialInternalForces.push_back(internalForces[0].norm()); +// } +// return axialInternalForces; +// }(); +// const auto pCurvNet = pMesh->registerForDrawing(); +// pCurvNet->addNodeScalarQuantity("axial forces", perVertexInternalForces_axial); +// std::vector perVertexInternalForces_torsion = [=]() { +// std::vector axialInternalForces; +// axialInternalForces.reserve(results.perVertexInternalForces.size()); +// for (const std::array &internalForces : results.perVertexInternalForces) { +// axialInternalForces.push_back(internalForces[1].norm()); +// } +// return axialInternalForces; +// }(); +// pCurvNet->addNodeScalarQuantity("torsional forces", perVertexInternalForces_torsion); +// std::vector perVertexInternalForces_firstBending = [=]() { +// std::vector axialInternalForces; +// axialInternalForces.reserve(results.perVertexInternalForces.size()); +// for (const std::array &internalForces : results.perVertexInternalForces) { +// axialInternalForces.push_back(internalForces[2].norm()); +// } +// return axialInternalForces; +// }(); +// pCurvNet->addNodeScalarQuantity("first bending forces", +// perVertexInternalForces_firstBending); +// std::vector perVertexInternalForces_secondBending = [=]() { +// std::vector axialInternalForces; +// axialInternalForces.reserve(results.perVertexInternalForces.size()); +// for (const std::array &internalForces : results.perVertexInternalForces) { +// axialInternalForces.push_back(internalForces[3].norm()); +// } +// return axialInternalForces; +// }(); +// pCurvNet->addNodeScalarQuantity("second bending forces", +// perVertexInternalForces_secondBending); +// polyscope::show(); #endif results.rotationalDisplacementQuaternion.resize(pMesh->VN()); results.debug_q_f1.resize(pMesh->VN()); @@ -2191,11 +2224,14 @@ void DRMSimulationModel::printDebugInfo() const } #ifdef POLYSCOPE_DEFINED -void DRMSimulationModel::draw(const std::string &screenshotsFolder) +void DRMSimulationModel::draw(const std::shared_ptr &pJob, + const std::string &screenshotsFolder) { // update positions // polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false); - polyscope::CurveNetwork *meshPolyscopeHandle = polyscope::getCurveNetwork(meshPolyscopeLabel); + auto deformedMeshPolyscopeLabel = meshPolyscopeLabel + "_" + pMesh->getLabel(); + polyscope::CurveNetwork *meshPolyscopeHandle = polyscope::getCurveNetwork( + deformedMeshPolyscopeLabel); meshPolyscopeHandle->updateNodePositions(pMesh->getEigenVertices()); // Vertex quantities @@ -2259,11 +2295,12 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false); meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(true); meshPolyscopeHandle->addNodeVectorQuantity("Internal force", internalForces)->setEnabled(false); - meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false); - meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true); meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm) ->setEnabled(true); meshPolyscopeHandle->setRadius(0.002); + // meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(true); + // meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true); + pJob->registerForDrawing(deformedMeshPolyscopeLabel, true); // polyscope::getCurveNetwork(meshPolyscopeLabel) // ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) // ->setEnabled(false); @@ -2284,7 +2321,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) // ->addNodeVectorQuantity("Second bending force-Rotation", // internalSecondBendingRotationForces) // ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel) ->addNodeScalarQuantity("nR", nRs) ->setEnabled(false); // polyscope::getCurveNetwork(meshPolyscopeLabel) @@ -2293,10 +2330,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) // polyscope::getCurveNetwork(meshPolyscopeLabel) // ->addNodeScalarQuantity("theta2", theta2) // ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel) ->addNodeVectorQuantity("Residual force", residualForces) ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel) ->addNodeScalarQuantity("Residual force norm", residualForcesNorm) ->setEnabled(false); // polyscope::getCurveNetwork(meshPolyscopeLabel) @@ -2315,11 +2352,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) I3[ei] = pMesh->elements[e].dimensions.inertia.I3; } - polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); - polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addEdgeScalarQuantity("I2", I2); - polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3); + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("A", A); + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("J", J); + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("I2", I2); + polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3); // Specify the callback static bool calledOnce = false; @@ -2574,7 +2610,19 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr &pJob) { auto beginTime = std::chrono::high_resolution_clock::now(); + updateNodalMasses(); + // constexpr bool useDRM = true; + //#ifdef USE_ENSMALLEN + // if (!useDRM) { + // setJob(pJob); + // // ens::L_BFGS optimizer(20); + // ens::SA optimizer; + // arma::mat x(pJob->pMesh->VN() * NumDoF, 1); + // optimizer.Optimize(*this, x); + // // getD + // } else { + //#endif // std::unordered_map nodalExternalForces = pJob->nodalExternalForces; // double totalExternalForcesNorm = 0; // Vector6d sumOfExternalForces(0); @@ -2594,197 +2642,224 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN() << " nodes and " - << pMesh->EN() << " elements." << std::endl; - } + if (mSettings.beVerbose) { + std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " + << pMesh->EN() << " elements." << std::endl; + } - const size_t movingAverageSampleSize = 200; - std::queue residualForcesMovingAverageHistorySample; - std::vector percentageOfConvergence; - // double residualForcesMovingAverageDerivativeMax = 0; - while (!mSettings.maxDRMIterations.has_value() - || mCurrentSimulationStep < mSettings.maxDRMIterations.value()) { - if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { + const size_t movingAverageSampleSize = 200; + std::queue residualForcesMovingAverageHistorySample; + std::vector percentageOfConvergence; + // double residualForcesMovingAverageDerivativeMax = 0; + while (!mSettings.maxDRMIterations.has_value() + || mCurrentSimulationStep < mSettings.maxDRMIterations.value()) { + // if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { // std::filesystem::create_directory("./PatternOptimizationNonConv"); // pJob->save("./PatternOptimizationNonConv"); // Dt = mSettings.Dtini; - } - // if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { - // Dt = mSettings.Dtini; - // } - // while (true) { - updateNormalDerivatives(); - updateT1Derivatives(); - updateRDerivatives(); - updateT2Derivatives(); - updateT3Derivatives(); - const bool shouldBreak = mCurrentSimulationStep == 12970; - updateResidualForcesOnTheFly(constrainedVertices); - - // 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( - - pJob->nodalForcedDisplacements); - } - // if (!pJob->nodalForcedNormals.empty()) { - // applyForcedNormals(pJob->nodalForcedNormals); - // } - updateElementalLengths(); - mCurrentSimulationStep++; - if (std::isnan(pMesh->currentTotalKineticEnergy)) { - std::cout << pMesh->currentTotalKineticEnergy << std::endl; - if (mSettings.beVerbose) { - std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep - << ". Exiting.." << std::endl; + // } + // if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { + // Dt = mSettings.Dtini; + // } + // while (true) { + updateNormalDerivatives(); + updateT1Derivatives(); + updateRDerivatives(); + updateT2Derivatives(); + updateT3Derivatives(); + const bool shouldBreak = mCurrentSimulationStep == 3935; + if (shouldBreak) { + int i = 0; + i++; } - break; - } + updateResidualForcesOnTheFly(constrainedVertices); - if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) { - minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; - static int lastSavedStep = mCurrentSimulationStep; - if (mSettings.saveIntermediateBestStates.has_value() - && mSettings.saveIntermediateBestStates.value() - && mCurrentSimulationStep - lastSavedStep > 20000) { - lastSavedStep = mCurrentSimulationStep; - computeResults(pJob).save("./IntermediateResults_" + pJob->getLabel()); + // 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( + + pJob->nodalForcedDisplacements); } - } - - //normalized sum of displacements - // 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; - - //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->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.debugModeStep.has_value() - && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { - printCurrentState(); - - // auto t2 = std::chrono::high_resolution_clock::now(); - // const double executionTime_min - // = std::chrono::duration_cast(t2 - beginTime).count(); - // std::cout << "Execution time(min):" << executionTime_min << std::endl; - if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { - std::cout << "Best percentage of target (average):" - << 100 * mSettings.averageResidualForcesCriterionThreshold.value() - * pMesh->totalExternalForcesNorm * pMesh->VN() - / minTotalResidualForcesNorm - << "%" << std::endl; - } else { - std::cout << "Best percentage of target:" - << 100 * mSettings.totalExternalForcesNormPercentageTermination - * pMesh->totalExternalForcesNorm / minTotalResidualForcesNorm - << "%" << std::endl; + // if (!pJob->nodalForcedNormals.empty()) { + // applyForcedNormals(pJob->nodalForcedNormals); + // } + updateElementalLengths(); + mCurrentSimulationStep++; + if (std::isnan(pMesh->currentTotalKineticEnergy)) { + std::cout << pMesh->currentTotalKineticEnergy << std::endl; + if (mSettings.beVerbose) { + std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep + << ". Exiting.." << std::endl; + } + break; } - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver", - // history.residualForcesMovingAverage); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver deriv", - // movingAveragesDerivatives); - // draw(); - // SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps", - // "Time/(#nodes*#iterations)", - // timePerNodePerIterationHistory); - } - if ((mSettings.shouldCreatePlots || mSettings.debugModeStep.has_value()) - && mCurrentSimulationStep != 0) { - history.stepPulse(*pMesh); - percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold - / pMesh->totalResidualForcesNorm); - } + if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) { + minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; + static int lastSavedStep = mCurrentSimulationStep; + if (mSettings.saveIntermediateBestStates.has_value() + && mSettings.saveIntermediateBestStates.value() + && mCurrentSimulationStep - lastSavedStep > 20000) { + lastSavedStep = mCurrentSimulationStep; + computeResults(pJob).save("./IntermediateResults_" + pJob->getLabel()); + } + } - if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() - && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver deriv", - // movingAveragesDerivatives_norm); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver", - // history.residualForcesMovingAverage, - // {}, - // history.redMarks); - SimulationResultsReporter::createPlot("Number of Steps", - "Log of Residual Forces", - history.logResidualForces, - {}, - history.redMarks); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Log of Kinetic energy", - // history.kineticEnergy, - // {}, - // history.redMarks); - // SimulationResultsReporter reporter; - // reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel()); - // SimulationResultsReporter::createPlot("Number of Iterations", - // "Sum of normalized displacement norms", - // history.sumOfNormalizedDisplacementNorms /*, - // std::filesystem::path("./") - // .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") - // .string()*/ - // , - // {}, - // history.redMarks); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Percentage of convergence", - // percentageOfConvergence, - // {}, - // history.redMarks); - } + // pMesh->previousTotalPotentialEnergykN = + // pMesh->currentTotalPotentialEnergykN; + // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + // timePerNodePerIterationHistor.push_back(timePerNodePerIteration); + if (mSettings.beVerbose && mSettings.debugModeStep.has_value() + && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { + printCurrentState(); + + // auto t2 = std::chrono::high_resolution_clock::now(); + // const double executionTime_min + // = std::chrono::duration_cast(t2 - beginTime).count(); + // std::cout << "Execution time(min):" << executionTime_min << std::endl; + if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { + std::cout << "Best percentage of target (average):" + << 100 * mSettings.averageResidualForcesCriterionThreshold.value() + * pMesh->totalExternalForcesNorm * pMesh->VN() + / minTotalResidualForcesNorm + << "%" << std::endl; + } else { + std::cout << "Best percentage of target:" + << 100 * mSettings.totalExternalForcesNormPercentageTermination + * pMesh->totalExternalForcesNorm / minTotalResidualForcesNorm + << "%" << std::endl; + } + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver", + // history.residualForcesMovingAverage); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // movingAveragesDerivatives); + // draw(); + // SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps", + // "Time/(#nodes*#iterations)", + // timePerNodePerIterationHistory); + } + + if ((mSettings.shouldCreatePlots || mSettings.debugModeStep.has_value()) + && mCurrentSimulationStep != 0) { + //normalized sum of displacements + double sumOfDisplacementsNorm = 0; + for (size_t vi = 0; vi < pMesh->VN(); vi++) { + sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); + } + pMesh->perVertexAverageNormalizedDisplacementNorm = sumOfDisplacementsNorm + / (pMesh->bbox.Diag() + * pMesh->VN()); + + //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->residualForcesMovingAverageDerivativeNorm + = std::abs((residualForcesMovingAverageHistorySample.back() + - residualForcesMovingAverageHistorySample.front())) + / (movingAverageSampleSize); + // residualForcesMovingAverageDerivativeMax + // = std::max(pMesh->residualForcesMovingAverageDerivativeNorm, + // residualForcesMovingAverageDerivativeMax); + // pMesh->residualForcesMovingAverageDerivativeNorm + // /= residualForcesMovingAverageDerivativeMax; + // std::cout << "Normalized derivative:" + // << residualForcesMovingAverageDerivativeNorm + // << std::endl; + } + history.stepPulse(*pMesh); + percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold + / pMesh->totalResidualForcesNorm); + } + + if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() + && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // history + // .residualForcesMovingAverageDerivativesLog); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver", + // history.residualForcesMovingAverage, + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Sum of axial forces norms", + // history.logOfSumOfAxialForcesNorm, + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Sum of torsion forces norms", + // history.logOfSumOfTorsionForcesNorm, + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Sum of 1stBending forces norms", + // history.logOfSumOfFirstBendingForcesNorm, + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Sum of 2ndBending forces norms", + // history.logOfSumOfSecondBendingForcesNorm, + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Residual Forces", + // history.logResidualForces, + // {}, + // history.redMarks); + SimulationResultsReporter::createPlot("Number of Steps", + "Log of Kinetic energy", + history.kineticEnergy, + {}, + history.redMarks); + // SimulationResultsReporter reporter; + // reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel()); + // SimulationResultsReporter::createPlot( + // "Number of Iterations", + // "Per vertex average normalized displacement norms", + // history.perVertexAverageNormalizedDisplacementNorm /*, + // std::filesystem::path("./") + // .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") + // .string()*/ + // , + // {}, + // history.redMarks); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Percentage of convergence", + // percentageOfConvergence, + // {}, + // history.redMarks); + } #ifdef POLYSCOPE_DEFINED if (mSettings.shouldDraw && mSettings.debugModeStep.has_value() @@ -2796,7 +2871,7 @@ currentSimulationStep > maxDRMIterations*/ // .append("Debugging_files") // .append("Screenshots") // .string(); - draw(); + draw(pJob); // yValues = history.kineticEnergy; // plotHandle = matplot::scatter(xPlot, yValues); // label = "Log of Kinetic energy"; @@ -2837,17 +2912,14 @@ currentSimulationStep > maxDRMIterations*/ } break; } - // Residual forces norm convergence - if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - // && mCurrentSimulationStep > movingAverageSampleSize - ) - /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 - && mCurrentSimulationStep > 1)*/ - /*|| -mesh->previousTotalPotentialEnergykN > -mesh->currentTotalPotentialEnergykN*/ - /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { - // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + + const bool isKineticEnergyPeak = (mSettings.useTotalRotationalKineticEnergyForKineticDamping + && pMesh->previousTotalRotationalKineticEnergy + > pMesh->currentTotalRotationalKineticEnergy) + || (mSettings.useTranslationalKineticEnergyForKineticDamping + && pMesh->previousTranslationalKineticEnergy + > pMesh->currentTotalTranslationalKineticEnergy); + if (isKineticEnergyPeak) { const bool fullfillsKineticEnergyTerminationCriterion = mSettings.totalTranslationalKineticEnergyThreshold.has_value() && pMesh->currentTotalTranslationalKineticEnergy @@ -2899,7 +2971,7 @@ mesh->currentTotalPotentialEnergykN*/ // } } - if (mSettings.useKineticDamping) { + if (mSettings.useTranslationalKineticEnergyForKineticDamping) { applyKineticDamping(pJob); Dt *= mSettings.xi; if (mSettings.shouldCreatePlots) { @@ -2912,22 +2984,24 @@ mesh->currentTotalPotentialEnergykN*/ // std::cout << "Number of dampings:" << numOfDampings << endl; // draw(); } - } - auto endTime = std::chrono::high_resolution_clock::now(); + } + //#ifdef USE_ENSMALLEN + // } + //#endif + auto endTime = std::chrono::high_resolution_clock::now(); + SimulationResults results = computeResults(pJob); + results.executionTime + = std::chrono::duration_cast(endTime - beginTime).count(); - SimulationResults results = computeResults(pJob); - results.executionTime = std::chrono::duration_cast(endTime - beginTime) - .count(); - - if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) { - SimulationResultsReporter reporter; - reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); - } + if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) { + SimulationResultsReporter reporter; + reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); + } #ifdef POLYSCOPE_DEFINED if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { - polyscope::removeCurveNetwork(meshPolyscopeLabel); - polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); + polyscope::removeCurveNetwork(meshPolyscopeLabel + "_" + pMesh->getLabel()); + polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel()); } #endif return results; @@ -2936,6 +3010,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh != nullptr); @@ -2953,16 +3028,26 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces) { - externalForce.second = externalForce.second / forceScaleFactor; - } - applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); - shouldTemporarilyDampForces = true; + if (simulationResults_fullPatternLinearModel.converged) { + // simulationResults_fullPatternLinearModel.save(fullResultsFolderPath); + for (auto &externalForce : pJob->nodalExternalForces) { + externalForce.second = externalForce.second / forceScaleFactor; + } + + applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); + shouldTemporarilyDampForces = true; + } } SimulationResults results = executeSimulation(pJob); //calling the virtual function + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + if (settings.beVerbose) { + std::cout << "Simulation converged in " + << std::chrono::duration_cast(end - begin).count() + / (60 * 1000.0) + << " min" << std::endl; + } return results; } @@ -2983,7 +3068,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob) //{ // for (int vi = 0; vi < pMesh->VN(); vi++) { // Node &node = pMesh->nodes[vi]; @@ -3000,6 +3086,32 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN(); vi++) { +// Node &node = pMesh->nodes[vi]; +// for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { +// node.displacements[dofi] = x(vi * DoF::NumDoF + dofi); +// } +// } + +// applyDisplacements(constrainedVertices); +// if (!pJob->nodalForcedDisplacements.empty()) { +// applyForcedDisplacements(pJob->nodalForcedDisplacements); +// } +// updateElementalLengths(); +// updateNormalDerivatives(); +// updateT1Derivatives(); +// updateRDerivatives(); +// updateT2Derivatives(); +// updateT3Derivatives(); +// updateResidualForcesOnTheFly(constrainedVertices); +// const double PE = computeTotalPotentialEnergy(); +// // const double PE = computeTotalInternalPotentialEnergy(); +// // std::cout << "PE:" << PE << std::endl; +// return pMesh->totalResidualForcesNorm; +//} + //double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g) //{ // for (int vi = 0; vi < pMesh->VN(); vi++) { @@ -3028,7 +3140,13 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh->VN() * NumDoF; xi++) { +// sumOfDisplacements += x(xi) * x(xi); +// } +// std::cout << "Sum of disp:" << std::sqrt(sumOfDisplacements) << std::endl; // return PE; //} @@ -3043,7 +3161,8 @@ void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePat json[GET_VARIABLE_NAME(Dtini)] = Dtini; json[GET_VARIABLE_NAME(xi)] = xi; json[GET_VARIABLE_NAME(gamma)] = gamma; - json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold; + json[GET_VARIABLE_NAME(useTranslationalKineticEnergyForKineticDamping)] + = useTranslationalKineticEnergyForKineticDamping; if (maxDRMIterations.has_value()) { json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value(); } @@ -3064,7 +3183,8 @@ void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePat if (viscousDampingFactor.has_value()) { json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value(); } - json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping; + json[GET_VARIABLE_NAME(useTranslationalKineticEnergyForKineticDamping)] + = useTranslationalKineticEnergyForKineticDamping; // if (intermediateResultsSaveStep.has_value()) { // json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value(); // } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index c1ed842..a49391e 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -16,7 +16,7 @@ #include #endif -struct SimulationJob; +class SimulationJob; enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF }; using DoFType = int; @@ -33,6 +33,8 @@ class DRMSimulationModel : public SimulationModel public: struct Settings { + bool useTranslationalKineticEnergyForKineticDamping{true}; + bool useTotalRotationalKineticEnergyForKineticDamping{false}; bool shouldDraw{false}; bool beVerbose{false}; bool shouldCreatePlots{false}; @@ -45,7 +47,7 @@ public: int gradualForcedDisplacementSteps{50}; // int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; - double totalResidualForcesNormThreshold{1e-8}; + double totalResidualForcesNormThreshold{1e-20}; double totalExternalForcesNormPercentageTermination{1e-5}; std::optional maxDRMIterations; std::optional debugModeStep; @@ -55,7 +57,6 @@ public: // std::optional intermediateResultsSaveStep; std::optional saveIntermediateBestStates; std::optional viscousDampingFactor; - bool useKineticDamping{true}; Settings() {} void save(const std::filesystem::path &jsonFilePath) const; bool load(const std::filesystem::path &filePath); @@ -153,7 +154,7 @@ private: const std::unordered_map> &fixedVertices); #ifdef POLYSCOPE_DEFINED - void draw(const std::string &screenshotsFolder = {}); + void draw(const std::shared_ptr &pJob, const std::string &screenshotsFolder = {}); #endif void updateNodalInternalForce(Vector6d &nodalInternalForce, @@ -262,12 +263,13 @@ private: SimulationResults executeSimulation(const std::shared_ptr &pJob, const Settings &settings, const SimulationResults &solutionGuess = SimulationResults()); -#ifdef USE_ENSMALLEN - double EvaluateWithGradient(const arma::mat &x, arma::mat &g); - void setJob(const std::shared_ptr &pJob); - SimulationMesh *getDeformedMesh(const arma::mat &x); - -#endif + //#ifdef USE_ENSMALLEN + // std::shared_ptr pJob; + // double EvaluateWithGradient(const arma::mat &x, arma::mat &g); + // void setJob(const std::shared_ptr &pJob); + // SimulationMesh *getDeformedMesh(const arma::mat &x, const std::shared_ptr &pJob); + // double Evaluate(const arma::mat &x); + //#endif static void runUnitTests(); };