diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 6f3cc15..0f868bc 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -794,7 +794,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( // omp_lock_t writelock; // omp_init_lock(&writelock); #ifdef ENABLE_OPENMP -#pragma omp parallel for schedule(static) num_threads(8) +#pragma omp parallel for schedule(static) num_threads(6) #endif for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; @@ -1151,6 +1151,9 @@ void DRMSimulationModel::updateNodalMasses() if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) { gamma *= 1e6 * (1 - static_cast(mCurrentSimulationStep) / untilStep); } + if (mCurrentSimulationStep == untilStep && shouldTemporarilyDampForces) { + Dt = mSettings.Dtini * 0.95; + } for (VertexType &v : pMesh->vert) { const size_t vi = pMesh->getIndex(v); double translationalSumSk = 0; @@ -2037,6 +2040,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrsave("./PatternOptimizationNonConv"); + // Dt = mSettings.Dtini; + } + if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { + Dt = mSettings.Dtini; + } // while (true) { updateNormalDerivatives(); updateT1Derivatives(); @@ -2144,6 +2156,42 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrpMesh->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); + } + #ifdef POLYSCOPE_DEFINED if (mSettings.shouldDraw && mSettings.isDebugMode && mCurrentSimulationStep % mSettings.drawingStep == 0) /* && @@ -2163,253 +2211,231 @@ currentSimulationStep > maxDRMIterations*/ // shouldUseKineticEnergyThreshold = true; } #endif - if (mSettings.shouldCreatePlots && mSettings.isDebugMode - && mCurrentSimulationStep % mSettings.debugModeStep == 0) { - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver deriv", - // movingAveragesDerivatives_norm); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Residual Forces mov aver", - // history.residualForcesMovingAverage); - SimulationResultsReporter::createPlot("Number of Steps", - "Log of Residual Forces", - history.logResidualForces, - {}, - history.redMarks); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Log of Kinetic energy", - // history.kineticEnergy); - // SimulationResultsReporter reporter; - // reporter.reportHistory(history, "Results"); - // SimulationResultsReporter::createPlot("Number of Steps", - // "Percentage of convergence", - // percentageOfConvergence, - // {}, - // history.redMarks); - } - // t = t + Dt; - // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy - // << std::endl; - // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm - // << std::endl; - // Residual forces norm convergence - if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - && mCurrentSimulationStep > movingAverageSampleSize - && (pJob->nodalForcedDisplacements.empty() - || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps) - /*|| + + // t = t + Dt; + // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy + // << std::endl; + // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm + // << std::endl; + // Residual forces norm convergence + if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + // && mCurrentSimulationStep > movingAverageSampleSize + && (pJob->nodalForcedDisplacements.empty() + || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) + /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 + && mCurrentSimulationStep > 1)*/ + /*|| mesh->previousTotalPotentialEnergykN > mesh->currentTotalPotentialEnergykN*/ - /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { - // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { - const bool fullfillsKineticEnergyTerminationCriterion - = mSettings.shouldUseTranslationalKineticEnergyThreshold - && pMesh->currentTotalTranslationalKineticEnergy - < mSettings.totalTranslationalKineticEnergyThreshold - && mCurrentSimulationStep > 20 && numOfDampings > 0; - const bool fullfillsResidualForcesNormTerminationCriterion - = pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; - const bool fullfillsMovingAverageNormTerminationCriterion - = pMesh->residualForcesMovingAverage - < mSettings.residualForcesMovingAverageNormThreshold; - /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ - const bool fullfillsMovingAverageDerivativeNormTerminationCriterion - = pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4; - const bool shouldTerminate - = fullfillsKineticEnergyTerminationCriterion - // || fullfillsMovingAverageNormTerminationCriterion - // || fullfillsMovingAverageDerivativeNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion; - if (shouldTerminate) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - if (fullfillsResidualForcesNormTerminationCriterion) { - std::cout << "Converged using residual forces norm threshold criterion" - << std::endl; - } else if (fullfillsKineticEnergyTerminationCriterion) { - std::cout << "Warning: The kinetic energy of the system was " - " used as a convergence criterion" - << std::endl; - } else if (fullfillsMovingAverageNormTerminationCriterion) { - std::cout << "Converged using residual forces moving average derivative norm " - "threshold criterion" - << std::endl; - } - } - break; - // } - } + /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { + // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + const bool fullfillsKineticEnergyTerminationCriterion + = mSettings.shouldUseTranslationalKineticEnergyThreshold + && pMesh->currentTotalTranslationalKineticEnergy + < mSettings.totalTranslationalKineticEnergyThreshold + && mCurrentSimulationStep > 20 && numOfDampings > 0; + const bool fullfillsResidualForcesNormTerminationCriterion + = pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold; + const bool fullfillsMovingAverageNormTerminationCriterion + = pMesh->residualForcesMovingAverage + < mSettings.residualForcesMovingAverageNormThreshold; + /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ + const bool fullfillsMovingAverageDerivativeNormTerminationCriterion + = pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4; + const bool shouldTerminate + = fullfillsKineticEnergyTerminationCriterion + // || fullfillsMovingAverageNormTerminationCriterion + // || fullfillsMovingAverageDerivativeNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion; + if (shouldTerminate) { + if (mSettings.beVerbose && !mSettings.isDebugMode) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + if (fullfillsResidualForcesNormTerminationCriterion) { + std::cout << "Converged using residual forces norm threshold criterion" + << std::endl; + } else if (fullfillsKineticEnergyTerminationCriterion) { + std::cout << "Warning: The kinetic energy of the system was " + " used as a convergence criterion" + << std::endl; + } else if (fullfillsMovingAverageNormTerminationCriterion) { + std::cout + << "Converged using residual forces moving average derivative norm " + "threshold criterion" + << std::endl; + } + } + break; + // } + } - // if (mSettings.beVerbose) { - // printCurrentState(); - // } - // 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.001; - // const double xi_min = 0.85; - // const double xi_init = 0.9969; - // if (mSettings.totalResidualForcesNormThreshold / pMesh->totalResidualForcesNorm - // > triggerPercentage) { - // mSettings.xi = ((xi_min - xi_init) - // * (mSettings.totalResidualForcesNormThreshold - // / pMesh->totalResidualForcesNorm) - // + xi_init - triggerPercentage * xi_min) - // / (1 - triggerPercentage); - // } - Dt *= mSettings.xi; - resetVelocities(); - ++numOfDampings; - if (mSettings.shouldCreatePlots) { - history.redMarks.push_back(mCurrentSimulationStep); - } - // std::cout << "Number of dampings:" << numOfDampings << endl; - // draw(); - } + // const double triggerPercentage = 0.01; + // const double xi_min = 0.55; + // const double xi_init = 0.9969; + // if (mSettings.totalResidualForcesNormThreshold / pMesh->totalResidualForcesNorm + // > triggerPercentage) { + // mSettings.xi = ((xi_min - xi_init) + // * (mSettings.totalResidualForcesNormThreshold + // / pMesh->totalResidualForcesNorm) + // + xi_init - triggerPercentage * xi_min) + // / (1 - triggerPercentage); + // } + Dt *= mSettings.xi; + // if (mSettings.isDebugMode) { + // std::cout << Dt << std::endl; + // } + resetVelocities(); + ++numOfDampings; + if (mSettings.shouldCreatePlots) { + history.redMarks.push_back(mCurrentSimulationStep); + } + // std::cout << "Number of dampings:" << numOfDampings << endl; + // draw(); + } } - SimulationResults results = computeResults(pJob); + SimulationResults results = computeResults(pJob); - if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) { - std::cout << "Did not reach equilibrium before reaching the maximum number " - "of DRM steps (" - << settings.maxDRMIterations << "). Breaking simulation" << std::endl; - results.converged = false; - } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { - std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; - results.converged = false; - - } else if (mSettings.beVerbose) { - auto t2 = std::chrono::high_resolution_clock::now(); - double simulationDuration = std::chrono::duration_cast(t2 - t1) - .count(); - simulationDuration /= 1000; - std::cout << "Simulation converged after " << simulationDuration << "s" << std::endl; - std::cout << "Time/(node*iteration) " - << simulationDuration / (static_cast(mCurrentSimulationStep) * pMesh->VN()) - << "s" << std::endl; - std::cout << "Number of dampings:" << numOfDampings << endl; - std::cout << "Number of steps:" << mCurrentSimulationStep << endl; - } - // mesh.printVertexCoordinates(mesh.VN() / 2); + if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) { + if (mSettings.beVerbose) { + std::cout << "Did not reach equilibrium before reaching the maximum number " + "of DRM steps (" + << settings.maxDRMIterations << "). Breaking simulation" << std::endl; + } + results.converged = false; + } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { + if (mSettings.beVerbose) { + std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; + } + results.converged = false; + } + // mesh.printVertexCoordinates(mesh.VN() / 2); #ifdef POLYSCOPE_DEFINED - if (settings.shouldDraw) { - draw(); - } + if (mSettings.shouldDraw && !mSettings.isDebugMode) { + draw(); + } #endif - if (results.converged) { - results.debug_drmDisplacements = results.displacements; - results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); + if (!std::isnan(pMesh->currentTotalKineticEnergy)) { + results.debug_drmDisplacements = results.displacements; + results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); - results.rotationalDisplacementQuaternion.resize(pMesh->VN()); - results.debug_q_f1.resize(pMesh->VN()); - results.debug_q_normal.resize(pMesh->VN()); - results.debug_q_nr.resize(pMesh->VN()); - for (int vi = 0; vi < pMesh->VN(); vi++) { - const Node &node = pMesh->nodes[vi]; - const Eigen::Vector3d nInitial_eigen = node.initialNormal.ToEigenVector(); - const Eigen::Vector3d nDeformed_eigen - = pMesh->vert[vi].cN().ToEigenVector(); + results.rotationalDisplacementQuaternion.resize(pMesh->VN()); + results.debug_q_f1.resize(pMesh->VN()); + results.debug_q_normal.resize(pMesh->VN()); + results.debug_q_nr.resize(pMesh->VN()); + for (int vi = 0; vi < pMesh->VN(); vi++) { + const Node &node = pMesh->nodes[vi]; + const Eigen::Vector3d nInitial_eigen = node.initialNormal + .ToEigenVector(); + const Eigen::Vector3d nDeformed_eigen + = pMesh->vert[vi].cN().ToEigenVector(); - Eigen::Quaternion q_normal; - q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); - Eigen::Quaternion q_nr_nDeformed; - q_nr_nDeformed = Eigen::AngleAxis(pMesh->nodes[vi].nR, nDeformed_eigen); - Eigen::Quaternion q_nr_nInit; - q_nr_nInit = Eigen::AngleAxis(pMesh->nodes[vi].nR, nInitial_eigen); - const auto w = q_nr_nDeformed.w(); - const auto theta = 2 * acos(q_nr_nDeformed.w()); - const auto nr = pMesh->nodes[vi].nR; - Eigen::Vector3d deformedNormal_debug(q_nr_nDeformed.x() * sin(theta / 2), - q_nr_nDeformed.y() * sin(theta / 2), - q_nr_nDeformed.z() * sin(theta / 2)); - deformedNormal_debug.normalize(); - // const double nr = nr_0To2pi <= M_PI ? nr_0To2pi : (nr_0To2pi - 2 * M_PI); - const double nr_debug = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? theta : -theta; - assert(pMesh->nodes[vi].nR - nr_debug < 1e-6); - VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1; - const VectorType &nDeformed = pMesh->vert[vi].cN(); - const VectorType referenceF1_deformed = (referenceT1_deformed - - (node.initialNormal - * (referenceT1_deformed * node.initialNormal))) - .Normalize(); + Eigen::Quaternion q_normal; + q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); + Eigen::Quaternion q_nr_nDeformed; + q_nr_nDeformed = Eigen::AngleAxis(pMesh->nodes[vi].nR, nDeformed_eigen); + Eigen::Quaternion q_nr_nInit; + q_nr_nInit = Eigen::AngleAxis(pMesh->nodes[vi].nR, nInitial_eigen); + const auto w = q_nr_nDeformed.w(); + const auto theta = 2 * acos(q_nr_nDeformed.w()); + const auto nr = pMesh->nodes[vi].nR; + Eigen::Vector3d deformedNormal_debug(q_nr_nDeformed.x() * sin(theta / 2), + q_nr_nDeformed.y() * sin(theta / 2), + q_nr_nDeformed.z() * sin(theta / 2)); + deformedNormal_debug.normalize(); + // const double nr = nr_0To2pi <= M_PI ? nr_0To2pi : (nr_0To2pi - 2 * M_PI); + const double nr_debug = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? theta : -theta; + assert(pMesh->nodes[vi].nR - nr_debug < 1e-6); + VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1; + const VectorType &nDeformed = pMesh->vert[vi].cN(); + const VectorType referenceF1_deformed + = (referenceT1_deformed + - (node.initialNormal * (referenceT1_deformed * node.initialNormal))) + .Normalize(); - const VectorType referenceT1_initial - = computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation, - pMesh->nodes[node.referenceElement->cV(1)].initialLocation); - // const VectorType &n_initial = node.initialNormal; - const VectorType referenceF1_initial = (referenceT1_initial - - (node.initialNormal - * (referenceT1_initial * node.initialNormal))) - .Normalize(); - Eigen::Quaternion q_f1_nInit; //nr is with respect to f1 - q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector(), - referenceF1_deformed.ToEigenVector()); + const VectorType referenceT1_initial + = computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation, + pMesh->nodes[node.referenceElement->cV(1)].initialLocation); + // const VectorType &n_initial = node.initialNormal; + const VectorType referenceF1_initial = (referenceT1_initial + - (node.initialNormal + * (referenceT1_initial * node.initialNormal))) + .Normalize(); + Eigen::Quaternion q_f1_nInit; //nr is with respect to f1 + q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector(), + referenceF1_deformed.ToEigenVector()); - Eigen::Quaternion q_f1_nDeformed; //nr is with respect to f1 - // const VectorType &n_initial = node.initialNormal; - const VectorType referenceF1_initial_def - = (referenceT1_initial - (nDeformed * (referenceT1_initial * nDeformed))).Normalize(); - const VectorType referenceF1_deformed_def - = (referenceT1_deformed - (nDeformed * (referenceT1_deformed * nDeformed))).Normalize(); - q_f1_nDeformed.setFromTwoVectors(referenceF1_initial_def.ToEigenVector(), - referenceF1_deformed_def.ToEigenVector()); - const auto debug_qf1_nDef = (q_f1_nDeformed * q_nr_nDeformed) * nDeformed_eigen; - const auto debug_qf1_nInit = (q_f1_nInit * q_nr_nInit) * nInitial_eigen; - const auto debug_deformedNormal_f1Init = (q_normal * (q_f1_nInit * q_nr_nInit)) - * nInitial_eigen; - const auto debug_deformedNormal_f1Def = ((q_nr_nDeformed * q_f1_nDeformed) * q_normal) - * nInitial_eigen; - // Eigen::Quaternion q_t1; - // q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), - // referenceT1_deformed.ToEigenVector()); + Eigen::Quaternion q_f1_nDeformed; //nr is with respect to f1 + // const VectorType &n_initial = node.initialNormal; + const VectorType referenceF1_initial_def + = (referenceT1_initial - (nDeformed * (referenceT1_initial * nDeformed))).Normalize(); + const VectorType referenceF1_deformed_def = (referenceT1_deformed + - (nDeformed + * (referenceT1_deformed * nDeformed))) + .Normalize(); + q_f1_nDeformed + .setFromTwoVectors(referenceF1_initial_def.ToEigenVector(), + referenceF1_deformed_def.ToEigenVector()); + const auto debug_qf1_nDef = (q_f1_nDeformed * q_nr_nDeformed) * nDeformed_eigen; + const auto debug_qf1_nInit = (q_f1_nInit * q_nr_nInit) * nInitial_eigen; + const auto debug_deformedNormal_f1Init = (q_normal * (q_f1_nInit * q_nr_nInit)) + * nInitial_eigen; + const auto debug_deformedNormal_f1Def = ((q_nr_nDeformed * q_f1_nDeformed) * q_normal) + * nInitial_eigen; + // Eigen::Quaternion q_t1; + // q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), + // referenceT1_deformed.ToEigenVector()); - results.debug_q_f1[vi] = q_f1_nInit; - results.debug_q_normal[vi] = q_normal; - results.debug_q_nr[vi] = q_nr_nInit; - results.rotationalDisplacementQuaternion[vi] - = (q_normal * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; - //Update the displacement vector to contain the euler angles - const Eigen::Vector3d eulerAngles - = results.rotationalDisplacementQuaternion[vi].toRotationMatrix().eulerAngles(0, 1, 2); - results.displacements[vi][3] = eulerAngles[0]; - results.displacements[vi][4] = eulerAngles[1]; - results.displacements[vi][5] = eulerAngles[2]; + results.debug_q_f1[vi] = q_f1_nInit; + results.debug_q_normal[vi] = q_normal; + results.debug_q_nr[vi] = q_nr_nInit; + results.rotationalDisplacementQuaternion[vi] + = (q_normal + * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; + //Update the displacement vector to contain the euler angles + const Eigen::Vector3d eulerAngles = results.rotationalDisplacementQuaternion[vi] + .toRotationMatrix() + .eulerAngles(0, 1, 2); + results.displacements[vi][3] = eulerAngles[0]; + results.displacements[vi][4] = eulerAngles[1]; + results.displacements[vi][5] = eulerAngles[2]; - // Eigen::Quaterniond q_test = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX()) - // * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY()) - // * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ()); + // Eigen::Quaterniond q_test = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX()) + // * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY()) + // * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ()); - // const double dot_test = results.rotationalDisplacementQuaternion[vi].dot(q_test); - // assert(dot_test > 1 - 1e5); + // const double dot_test = results.rotationalDisplacementQuaternion[vi].dot(q_test); + // assert(dot_test > 1 - 1e5); - int i = 0; - i++; - } - } + int i = 0; + i++; + } + } - if (mSettings.shouldCreatePlots && results.converged) { - SimulationResultsReporter reporter; - reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); - } + if (!mSettings.isDebugMode && mSettings.shouldCreatePlots) { + SimulationResultsReporter reporter; + reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); + } return results; } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 392fd95..4c62eb3 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -31,8 +31,7 @@ public: bool beVerbose{false}; bool shouldCreatePlots{false}; int drawingStep{1}; - double totalTranslationalKineticEnergyThreshold{1e-10}; - double totalExternalForcesNormPercentageTermination{1e-3}; + double totalTranslationalKineticEnergyThreshold{1e-8}; double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; double residualForcesMovingAverageNormThreshold{1e-8}; double Dtini{0.1}; diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index 08b92d1..40878ce 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -254,7 +254,8 @@ struct Colors { std::string label{"EmptyLabel"}; double time{-1}; - int numberOfSimulationCrashes{0}; + bool wasSuccessful{true}; + // int numberOfSimulationCrashes{0}; Settings settings; struct ObjectiveValues { @@ -273,6 +274,7 @@ struct Colors PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern vcg::Triangle3 baseTriangle; + std::string notes; struct JsonKeys { @@ -294,15 +296,19 @@ struct Colors //Save optimal X nlohmann::json json_optimizationResults; json_optimizationResults[JsonKeys::Label] = label; - std::string jsonValue_optimizationVariables; - for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { - jsonValue_optimizationVariables.append(optimalXNameValuePair.first + ","); - } - jsonValue_optimizationVariables.pop_back(); // for deleting the last comma - json_optimizationResults[JsonKeys::optimizationVariables] = jsonValue_optimizationVariables; + if (wasSuccessful) { + std::string jsonValue_optimizationVariables; + for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { + jsonValue_optimizationVariables.append(optimalXNameValuePair.first + ","); + } + jsonValue_optimizationVariables.pop_back(); // for deleting the last comma + json_optimizationResults[JsonKeys::optimizationVariables] + = jsonValue_optimizationVariables; - for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { - json_optimizationResults[optimalXNameValuePair.first] = optimalXNameValuePair.second; + for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { + json_optimizationResults[optimalXNameValuePair.first] = optimalXNameValuePair + .second; + } } //base triangle json_optimizationResults[JsonKeys::baseTriangle] @@ -650,7 +656,8 @@ struct Colors os << nameValuePair.first; } os << "Time(s)"; - os << "#Crashes"; + // os << "#Crashes"; + // os << "notes"; } void writeResultsTo(const Settings &settings_optimization, csvFile &os) const @@ -669,11 +676,12 @@ struct Colors } os << time; - if (numberOfSimulationCrashes == 0) { - os << "-"; - } else { - os << numberOfSimulationCrashes; - } + // if (numberOfSimulationCrashes == 0) { + // os << "-"; + // } else { + // os << numberOfSimulationCrashes; + // } + // os << notes; } }; diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index ce317b9..371caec 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -63,7 +63,7 @@ struct SimulationResultsReporter { const std::string &graphSuffix = std::string()) { const auto simulationResultPath = std::filesystem::path(reportFolderPath).append(history.label); - std::filesystem::create_directory(simulationResultPath); + std::filesystem::create_directories(simulationResultPath); createPlots(history, simulationResultPath, graphSuffix); } diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index e8c4f99..dba4356 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -177,7 +177,12 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced auto perEdgeResultPath = std::filesystem::path(resultsPath) .append(std::to_string(numberOfEdges)); if (std::filesystem::exists(perEdgeResultPath)) { - continue; + if (debugIsOn) { + std::filesystem::remove_all(perEdgeResultPath); + + } else { + continue; + } } std::filesystem::create_directory(perEdgeResultPath); computeValidPatterns(numberOfNodesPerSlot, @@ -517,8 +522,7 @@ void TopologyEnumerator::computeValidPatterns( // Check dangling edges with vcg method // const bool vcg_tiledPatternHasDangling = // tiledPatternGeometry.hasUntiledDanglingEdges(); - if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents && - !hasArticulationPoints*/) { + if (tiledPatternHasDanglingEdges && !hasFloatingComponents /*&& !hasArticulationPoints*/) { statistics.numberOfPatternsWithADanglingEdgeOrNode++; if (debugIsOn) { if (savePlyFiles) { @@ -539,8 +543,7 @@ void TopologyEnumerator::computeValidPatterns( } } - if (hasFloatingComponents /*&& !hasArticulationPoints && - !tiledPatternHasDanglingEdges*/) { + if (hasFloatingComponents && !hasArticulationPoints && !tiledPatternHasDanglingEdges) { statistics.numberOfPatternsWithMoreThanASingleCC++; if (debugIsOn) { if (savePlyFiles) { diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp index 2f630b7..c1c384d 100755 --- a/vcgtrimesh.cpp +++ b/vcgtrimesh.cpp @@ -15,12 +15,12 @@ bool VCGTriMesh::load(const std::string &filename) { mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper::IO_FACEINDEX; - if (nanoply::NanoPlyWrapper::LoadModel( - std::filesystem::absolute(filename).string().c_str(), *this, mask) + // if (nanoply::NanoPlyWrapper::LoadModel( + // std::filesystem::absolute(filename).string().c_str(), *this, mask) + // != 0) { + if (vcg::tri::io::Importer::Open(*this, + std::filesystem::absolute(filename).string().c_str()) != 0) { - // if (vcg::tri::io::Importer::Open(*this, - // std::filesystem::absolute(filename).string().c_str()) - // != 0) { std::cout << "Could not load tri mesh" << std::endl; return false; }