From b4d80ad6b88ea8f8a7f330aa3e7b253382f12260 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 13 Jul 2021 09:02:55 +0300 Subject: [PATCH 1/5] Windows refactoring --- simulationhistoryplotter.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index 371caec..d43bf67 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -64,7 +64,7 @@ struct SimulationResultsReporter { { const auto simulationResultPath = std::filesystem::path(reportFolderPath).append(history.label); std::filesystem::create_directories(simulationResultPath); - createPlots(history, simulationResultPath, graphSuffix); + createPlots(history, simulationResultPath.string(), graphSuffix); } void reportResults(const std::vector &results, From 8391646d564a6148495dda6e0330dc3af06f4476 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Wed, 14 Jul 2021 19:17:02 +0300 Subject: [PATCH 2/5] Average residual forces norm. Alpha angles stored in vectors. Residual forces criterion applicable not only at local kinetic energy maxima --- drmsimulationmodel.cpp | 161 +++++++++++++++++++++++------------------ drmsimulationmodel.hpp | 11 +-- simulationmesh.cpp | 8 +- simulationmesh.hpp | 9 ++- 4 files changed, 107 insertions(+), 82 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 4dd80cb..0c41318 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -222,7 +222,6 @@ void DRMSimulationModel::reset() mCurrentSimulationStep = 0; history.clear(); constrainedVertices.clear(); - rigidSupports.clear(); pMesh.reset(); plotYValues.clear(); plotHandle.reset(); @@ -232,6 +231,8 @@ void DRMSimulationModel::reset() numOfDampings = 0; shouldTemporarilyDampForces = false; externalLoadStep = 1; + isVertexConstrained.clear(); + minTotalResidualForcesNorm = std::numeric_limits::max(); } VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative( @@ -544,13 +545,18 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v) // Use nR as theta3 only for the first star edge return nR; } - std::vector incidentElementsIndices(node.incidentElements.size()); - for (int iei = 0; iei < incidentElementsIndices.size(); iei++) { - incidentElementsIndices[iei] = pMesh->getIndex(node.incidentElements[iei]); - } + // std::vector incidentElementsIndices(node.incidentElements.size()); + // for (int iei = 0; iei < incidentElementsIndices.size(); iei++) { + // incidentElementsIndices[iei] = pMesh->getIndex(node.incidentElements[iei]); + // } assert(pMesh->getIndex(e) == ei); - assert(node.alphaAngles.contains(ei)); - const double alphaAngle = node.alphaAngles.at(elem.ei); + // assert(node.alphaAngles.contains(ei)); + const double alphaAngle = std::find_if(node.alphaAngles.begin(), + node.alphaAngles.end(), + [&](const std::pair &p) { + return elem.ei == p.first; + }) + ->second; const EdgeType &refElem = *node.referenceElement; const double rotationAngle = nR + alphaAngle; @@ -591,8 +597,7 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, { const Node &node = pMesh->nodes[v]; const VertexIndex &vi = pMesh->nodes[v].vi; - const bool isRigidSupport = rigidSupports.contains(vi); - if (&e == node.referenceElement && !isRigidSupport) { + if (&e == node.referenceElement && !isRigidSupport[vi]) { if (dui.dofi == DoF::Nr && &dui.v == &v) { return 1; } else { @@ -610,7 +615,7 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, const VertexPointer &vp_jplus1 = e.cV(1); double derivativeTheta3_dofi = 0; - if (isRigidSupport) { + if (isRigidSupport[vi]) { const VectorType &t1Initial = computeT1Vector(pMesh->nodes[vp_j].initialLocation, pMesh->nodes[vp_jplus1].initialLocation); VectorType g1 = Cross(t1, t1Initial); @@ -800,9 +805,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(8) - //#endif +#ifdef ENABLE_OPENMP +#pragma omp parallel for schedule(static) num_threads(5) +#endif for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; const SimulationMesh::VertexType &ev_j = *e.cV(0); @@ -843,7 +848,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const size_t vi = edgeNode.vi; // #pragma omp parallel for schedule(static) num_threads(6) for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { - const bool isDofConstrainedFor_ev = fixedVertices.contains(edgeNode.vi) + const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi] && fixedVertices.at(edgeNode.vi).contains(dofi); if (!isDofConstrainedFor_ev) { DifferentiateWithRespectTo dui{ev, dofi}; @@ -924,7 +929,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( } if (edgeNode.referenceElement != &e) { const bool isDofConstrainedFor_refElemOtherVertex - = fixedVertices.contains(refElemOtherVertexNode.vi) + = isVertexConstrained[refElemOtherVertexNode.vi] && fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi); if (!isDofConstrainedFor_refElemOtherVertex) { DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi}; @@ -1012,6 +1017,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( } pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; pMesh->totalResidualForcesNorm = totalResidualForcesNorm; + if (mSettings.beVerbose) { + if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) { + minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; + } + } pMesh->averageResidualForcesNorm = totalResidualForcesNorm / pMesh->VN(); // pMesh->averageResidualForcesNorm = sumOfResidualForces.norm() / pMesh->VN(); @@ -1078,6 +1088,8 @@ void DRMSimulationModel::updateResidualForces() void DRMSimulationModel::computeRigidSupports() { + isRigidSupport.clear(); + isRigidSupport.resize(pMesh->VN(), false); for (const VertexType &v : pMesh->vert) { const VertexIndex vi = pMesh->nodes[v].vi; const bool isVertexConstrained = constrainedVertices.contains(vi); @@ -1090,7 +1102,7 @@ void DRMSimulationModel::computeRigidSupports() && constrainedDoFType.contains(DoF::Ny) && constrainedDoFType.contains(DoF::Nr); if (hasAllDoFTypeConstrained) { - rigidSupports.insert(vi); + isRigidSupport[vi] = true; } } } @@ -1180,8 +1192,9 @@ void DRMSimulationModel::updateNodalMasses() if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) { gamma *= 1e6 * (1 - static_cast(mCurrentSimulationStep) / untilStep); } - if (mCurrentSimulationStep == untilStep && shouldTemporarilyDampForces) { - Dt = mSettings.Dtini * 0.95; + if (mCurrentSimulationStep == static_cast(1.2 * untilStep) + && shouldTemporarilyDampForces) { + Dt = mSettings.Dtini; } for (VertexType &v : pMesh->vert) { const size_t vi = pMesh->getIndex(v); @@ -1328,7 +1341,7 @@ void DRMSimulationModel::updateNodeNr(VertexType &v) { const VertexIndex &vi = pMesh->nodes[v].vi; Node &node = pMesh->nodes[v]; - if (!rigidSupports.contains(vi)) { + if (!isRigidSupport[vi]) { node.nR = node.displacements[5]; } else { const EdgePointer &refElem = node.referenceElement; @@ -1509,11 +1522,11 @@ void DRMSimulationModel::resetVelocities() { for (const VertexType &v : pMesh->vert) { pMesh->nodes[v].velocity = - // pMesh->nodes[v].acceleration - // * Dt; // NOTE: Do I reset the previous - // velocity? - // reset - // current to 0 or this? + // pMesh->nodes[v].acceleration * Dt + // * 0.5; // NOTE: Do I reset the previous + // // velocity? + // // reset + // // current to 0 or this? 0; } updateKineticEnergy(); @@ -1616,7 +1629,7 @@ void DRMSimulationModel::updatePositionsOnTheFly( VectorType newNormal(nx, ny, nz); v.N() = newNormal; } - if (!rigidSupports.contains(vi)) { + if (!isRigidSupport[vi]) { node.nR = node.displacements[5]; } else { } @@ -1650,15 +1663,18 @@ void DRMSimulationModel::printCurrentState() const << std::endl; std::cout << "Kinetic energy:" << pMesh->currentTotalKineticEnergy << std::endl; static std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); - const auto timePerNodePerIteration = std::chrono::duration_cast( - std::chrono::steady_clock::now() - begin) - .count() - * 1e-6 / (mCurrentSimulationStep * pMesh->VN()); + const double timePerNodePerIteration = std::chrono::duration_cast( + std::chrono::steady_clock::now() - begin) + .count() + * 1e-6 + / (static_cast(mCurrentSimulationStep) + * pMesh->VN()); std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl; std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << std::endl; std::cout << "Mov aver deriv norm:" << pMesh->residualForcesMovingAverageDerivativeNorm << std::endl; std::cout << "xi:" << mSettings.xi << std::endl; + std::cout << "Dt:" << Dt << std::endl; } void DRMSimulationModel::printDebugInfo() const @@ -1936,7 +1952,7 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue node.displacements[5] = nr; } // const double nr_asin = q_nr.x() - if (rigidSupports.contains(vi)) { + if (isRigidSupport[vi]) { const EdgePointer &refElem = node.referenceElement; const VectorType &refT1 = computeT1Vector(refElem->cP(0), refElem->cP(1)); @@ -2046,8 +2062,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr::Box(*pMesh); computeRigidSupports(); + isVertexConstrained.resize(pMesh->VN(), false); for (auto fixedVertex : pJob->constrainedVertices) { - assert(fixedVertex.first < pMesh->VN()); + isVertexConstrained[fixedVertex.first] = true; } #ifdef POLYSCOPE_DEFINED @@ -2058,7 +2075,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrgetEigenEdges()); polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, pMesh->getEigenVertices(), - pMesh->getEigenEdges()); + pMesh->getEigenEdges()) + ->setRadius(0.002); // registerWorldAxes(); } #endif @@ -2120,9 +2138,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrsave("./PatternOptimizationNonConv"); // Dt = mSettings.Dtini; } - if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { - Dt = mSettings.Dtini; - } + // if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { + // Dt = mSettings.Dtini; + // } // while (true) { updateNormalDerivatives(); updateT1Derivatives(); @@ -2208,19 +2226,19 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr(t2 - beginTime).count(); // std::cout << "Execution time(min):" << executionTime_min << std::endl; if (mSettings.useAverage) { - std::cout << "Percentage of target (average):" + std::cout << "Best percentage of target (average):" << 100 * mSettings.averageResidualForcesCriterionThreshold * totalExternalForcesNorm - / (pMesh->totalResidualForcesNorm / pMesh->VN()) + / (minTotalResidualForcesNorm / pMesh->VN()) << "%" << std::endl; } - std::cout << "Percentage of target:" + std::cout << "Best percentage of target:" << 100 * mSettings.totalExternalForcesNormPercentageTermination - * totalExternalForcesNorm / pMesh->totalResidualForcesNorm + * totalExternalForcesNorm / minTotalResidualForcesNorm << "%" << std::endl; - SimulationResultsReporter::createPlot("Number of Steps", - "Residual Forces mov aver", - history.residualForcesMovingAverage); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver", + // history.residualForcesMovingAverage); // SimulationResultsReporter::createPlot("Number of Steps", // "Residual Forces mov aver deriv", // movingAveragesDerivatives); @@ -2297,8 +2315,18 @@ currentSimulationStep > maxDRMIterations*/ // << std::endl; // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // << std::endl; + const bool fullfillsResidualForcesNormTerminationCriterion + = !mSettings.useAverage + && pMesh->totalResidualForcesNorm / totalExternalForcesNorm + < mSettings.totalExternalForcesNormPercentageTermination; + const bool fullfillsAverageResidualForcesNormTerminationCriterion + = mSettings.useAverage + && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm + < mSettings.averageResidualForcesCriterionThreshold; // Residual forces norm convergence - if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + || fullfillsAverageResidualForcesNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion) // && mCurrentSimulationStep > movingAverageSampleSize && (pJob->nodalForcedDisplacements.empty() || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) @@ -2314,13 +2342,6 @@ mesh->currentTotalPotentialEnergykN*/ && pMesh->currentTotalTranslationalKineticEnergy < mSettings.totalTranslationalKineticEnergyThreshold && mCurrentSimulationStep > 20 && numOfDampings > 0; - const bool fullfillsResidualForcesNormTerminationCriterion - = pMesh->totalResidualForcesNorm / totalExternalForcesNorm - < mSettings.totalExternalForcesNormPercentageTermination; - const bool fullfillsAverageResidualForcesNormTerminationCriterion - = mSettings.useAverage - && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm - < mSettings.averageResidualForcesCriterionThreshold; const bool fullfillsMovingAverageNormTerminationCriterion = pMesh->residualForcesMovingAverage < mSettings.residualForcesMovingAverageNormThreshold; @@ -2341,7 +2362,7 @@ mesh->currentTotalPotentialEnergykN*/ std::cout << "Converged using residual forces norm threshold criterion" << std::endl; } else if (fullfillsKineticEnergyTerminationCriterion) { - std::cout << "Warning: The kinetic energy of the system was " + std::cout << "The kinetic energy of the system was " " used as a convergence criterion" << std::endl; } else if (fullfillsMovingAverageNormTerminationCriterion) { @@ -2377,25 +2398,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 * 0.5 * 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; @@ -2408,11 +2429,11 @@ mesh->currentTotalPotentialEnergykN*/ // + xi_init - triggerPercentage * xi_min) // / (1 - triggerPercentage); // } + resetVelocities(); Dt *= mSettings.xi; // if (mSettings.isDebugMode) { // std::cout << Dt << std::endl; // } - resetVelocities(); ++numOfDampings; if (mSettings.shouldCreatePlots) { history.redMarks.push_back(mCurrentSimulationStep); diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 6b98674..f415b80 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -45,7 +45,7 @@ public: double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; bool useAverage{false}; - double averageResidualForcesCriterionThreshold{1e-3}; + double averageResidualForcesCriterionThreshold{1e-5}; Settings() {} }; @@ -61,6 +61,9 @@ private: std::vector plotYValues; size_t numOfDampings{0}; int externalLoadStep{1}; + std::vector isVertexConstrained; + std::vector isRigidSupport; + double minTotalResidualForcesNorm{std::numeric_limits::max()}; const std::string meshPolyscopeLabel{"Simulation mesh"}; std::unique_ptr pMesh; @@ -69,7 +72,6 @@ private: // Eigen::Tensor theta3Derivatives; // std::unordered_map theta3Derivatives; bool shouldApplyInitialDistortion = false; - std::unordered_set rigidSupports; void reset(); void updateNodalInternalForces( @@ -156,10 +158,9 @@ private: void updateElementalFrames(); - VectorType computeDerivativeOfR(const EdgeType &e, - const DifferentiateWithRespectTo &dui) const; + VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const; - bool isRigidSupport(const VertexType &v) const; + // bool isRigidSupport(const VertexType &v) const; static double computeDerivativeOfNorm(const VectorType &x, const VectorType &derivativeOfX); diff --git a/simulationmesh.cpp b/simulationmesh.cpp index 997364d..1bd4076 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -140,6 +140,7 @@ void SimulationMesh::initializeNodes() { const EdgeType &referenceElement = *node.referenceElement; const VectorType t01 = computeT1Vector(referenceElement.cP(0), referenceElement.cP(1)); const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize(); + node.alphaAngles.reserve(incidentElements.size()); for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) { assert(referenceElement.cV(0) == ep->cV(0) || referenceElement.cV(0) == ep->cV(1) @@ -148,7 +149,7 @@ void SimulationMesh::initializeNodes() { const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize(); const EdgeIndex ei = getIndex(ep); const double alphaAngle = computeAngle(f01, f1, v.cN()); - node.alphaAngles[ei] = alphaAngle; + node.alphaAngles.emplace_back(std::make_pair(ei, alphaAngle)); } } } @@ -182,7 +183,8 @@ void SimulationMesh::reset() { const EdgeType &referenceElement = *getReferenceElement(v); const VectorType t01 = computeT1Vector(referenceElement.cP(0), referenceElement.cP(1)); const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize(); - + node.alphaAngles.clear(); + node.alphaAngles.reserve(node.incidentElements.size()); for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) { assert(referenceElement.cV(0) == ep->cV(0) || referenceElement.cV(0) == ep->cV(1) || referenceElement.cV(1) == ep->cV(0) || referenceElement.cV(1) == ep->cV(1)); @@ -190,7 +192,7 @@ void SimulationMesh::reset() { const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize(); const EdgeIndex ei = getIndex(ep); const double alphaAngle = computeAngle(f01, f1, v.cN()); - node.alphaAngles[ei] = alphaAngle; + node.alphaAngles.emplace_back(std::make_pair(ei, alphaAngle)); } } } diff --git a/simulationmesh.hpp b/simulationmesh.hpp index af36c7c..cecadf8 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -157,10 +157,11 @@ struct Node { double kineticEnergy{0}; Vector6d displacements{0}; double nR{0}; - std::unordered_map - alphaAngles; // contains the initial angles between the first star element - // incident to this node and the other elements of the star - // has size equal to the valence of the vertex + // std::unordered_map + // alphaAngles; // contains the initial angles between the first star element + // // incident to this node and the other elements of the star + // // has size equal to the valence of the vertex + std::vector> alphaAngles; std::vector incidentElements; std::vector derivativeOfNormal; From 28db0fa63e38f5cb47ed777ad8b7e2247b70a8c2 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Fri, 16 Jul 2021 14:15:47 +0300 Subject: [PATCH 3/5] Removed resetting of dt when guessing the solution with the linear model --- drmsimulationmodel.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 0c41318..4e682d1 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -805,9 +805,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(5) -#endif + //#ifdef ENABLE_OPENMP + //#pragma omp parallel for schedule(static) num_threads(5) + //#endif for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; const SimulationMesh::VertexType &ev_j = *e.cV(0); @@ -1005,8 +1005,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( Node::Forces &force = pMesh->nodes[vi].force; const Vector6d &nodeResidualForce = force.residual; sumOfResidualForces = sumOfResidualForces + nodeResidualForce; - const double residualForceNorm = nodeResidualForce.norm(); - // const double residualForceNorm = nodeResidualForce.getTranslation().norm(); + // const double residualForceNorm = nodeResidualForce.norm(); + const double residualForceNorm = nodeResidualForce.getTranslation().norm(); const bool shouldTrigger = std::isnan(residualForceNorm); totalResidualForcesNorm += residualForceNorm; #ifdef POLYSCOPE_DEFINED @@ -1192,10 +1192,10 @@ void DRMSimulationModel::updateNodalMasses() if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) { gamma *= 1e6 * (1 - static_cast(mCurrentSimulationStep) / untilStep); } - if (mCurrentSimulationStep == static_cast(1.2 * untilStep) - && shouldTemporarilyDampForces) { - Dt = mSettings.Dtini; - } + // if (mCurrentSimulationStep == static_cast(1.4 * untilStep) + // && shouldTemporarilyDampForces) { + // Dt = mSettings.Dtini; + // } for (VertexType &v : pMesh->vert) { const size_t vi = pMesh->getIndex(v); double translationalSumSk = 0; @@ -2354,7 +2354,7 @@ mesh->currentTotalPotentialEnergykN*/ // || fullfillsMovingAverageDerivativeNormTerminationCriterion || fullfillsAverageResidualForcesNormTerminationCriterion || fullfillsResidualForcesNormTerminationCriterion; - if (shouldTerminate) { + if (shouldTerminate && mCurrentSimulationStep > 100) { if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) { std::cout << "Simulation converged." << std::endl; printCurrentState(); From 613bb85a4dfb2e0d48da6ed27bce45ecc69259c2 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Fri, 16 Jul 2021 14:16:53 +0300 Subject: [PATCH 4/5] Added kinetic damping markers to the plots --- simulationhistoryplotter.hpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index 371caec..7304715 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -152,7 +152,8 @@ struct SimulationResultsReporter { history.kineticEnergy, std::filesystem::path(graphsFolder) .append("KineticEnergyLog_" + graphSuffix + ".png") - .string()); + .string(), + history.redMarks); } if (!history.logResidualForces.empty()) { @@ -161,7 +162,8 @@ struct SimulationResultsReporter { history.logResidualForces, std::filesystem::path(graphsFolder) .append("ResidualForcesLog_" + graphSuffix + ".png") - .string()); + .string(), + history.redMarks); } if (!history.potentialEnergies.empty()) { @@ -170,7 +172,8 @@ struct SimulationResultsReporter { history.potentialEnergies, std::filesystem::path(graphsFolder) .append("PotentialEnergy_" + graphSuffix + ".png") - .string()); + .string(), + history.redMarks); } if (!history.residualForcesMovingAverage.empty()) { createPlot("Number of Iterations", @@ -178,7 +181,8 @@ struct SimulationResultsReporter { history.residualForcesMovingAverage, std::filesystem::path(graphsFolder) .append("ResidualForcesMovingAverage_" + graphSuffix + ".png") - .string()); + .string(), + history.redMarks); } // if (!history.residualForcesMovingAverageDerivativesLog.empty()) { // createPlot("Number of Iterations", @@ -194,7 +198,8 @@ struct SimulationResultsReporter { history.sumOfNormalizedDisplacementNorms, std::filesystem::path(graphsFolder) .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") - .string()); + .string(), + history.redMarks); } } }; From 846289448bb1c5eb0f77843795b6db31367db649 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 20 Jul 2021 16:32:25 +0300 Subject: [PATCH 5/5] Refactoring. Fixed bug --- drmsimulationmodel.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 4e682d1..e33b7b4 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -692,8 +692,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, derivativeF1 = derivativeT1_k - (n * (derivativeT1_k * n)); const double f1Norm = f1.Norm(); - const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1) - / (f1Norm * f1Norm * f1Norm) + const double derivativeF1Norm = f1 * derivativeF1 / f1Norm; + const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm) + derivativeF1 / f1Norm; derivativeF2 = derivativeF1Normalized * cosRotationAngle @@ -713,8 +713,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e, const VectorType &derivativeOfNormal = node.derivativeOfNormal[dofi]; derivativeF1 = -(n * (t1_k * derivativeOfNormal) + derivativeOfNormal * (t1_k * n)); const double f1Norm = f1.Norm(); - const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1) - / (f1Norm * f1Norm * f1Norm) + const double derivativeF1Norm = f1 * derivativeF1 / f1Norm; + const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm) + derivativeF1 / f1Norm; derivativeF2 = derivativeF1Normalized * cosRotationAngle + @@ -1005,9 +1005,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( Node::Forces &force = pMesh->nodes[vi].force; const Vector6d &nodeResidualForce = force.residual; sumOfResidualForces = sumOfResidualForces + nodeResidualForce; - // const double residualForceNorm = nodeResidualForce.norm(); - const double residualForceNorm = nodeResidualForce.getTranslation().norm(); - const bool shouldTrigger = std::isnan(residualForceNorm); + const double residualForceNorm = nodeResidualForce.norm(); + // const double residualForceNorm = nodeResidualForce.getTranslation().norm(); + const bool shouldTrigger = mCurrentSimulationStep == 500; totalResidualForcesNorm += residualForceNorm; #ifdef POLYSCOPE_DEFINED if (std::isnan(force.residual.norm())) { @@ -1751,7 +1751,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0]; } meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false); - meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->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); @@ -2063,7 +2063,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr::Box(*pMesh); computeRigidSupports(); isVertexConstrained.resize(pMesh->VN(), false); - for (auto fixedVertex : pJob->constrainedVertices) { + for (auto fixedVertex : constrainedVertices) { isVertexConstrained[fixedVertex.first] = true; }