diff --git a/CMakeLists.txt b/CMakeLists.txt index e43756a..a807c52 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,6 +50,7 @@ download_project(PROJ TBB PREFIX ${EXTERNAL_DEPS_DIR} ${UPDATE_DISCONNECTED_IF_AVAILABLE} ) +option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastructure" OFF) add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) link_directories(${TBB_BINARY_DIR}) message(${TBB_BINARY_DIR}) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 9649614..0de885c 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -1,4 +1,5 @@ #include "drmsimulationmodel.hpp" +#include "linearsimulationmodel.hpp" #include "simulationhistoryplotter.hpp" #include #include @@ -810,7 +811,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( //#pragma omp parallel for schedule(static) num_threads(5) //#endif std::for_each( +#ifdef ENABLE_PARALLEL_DRM std::execution::par_unseq, +#endif pMesh->edge.begin(), pMesh->edge.end(), [&](const EdgeType &e) @@ -1049,6 +1052,7 @@ void DRMSimulationModel::updateNodalExternalForces( const std::unordered_map> &fixedVertices) { externalMomentsNorm = 0; + double totalExternalForcesNorm = 0; for (const std::pair &nodalForce : nodalForces) { const VertexIndex nodeIndex = nodalForce.first; const bool isNodeConstrained = fixedVertices.contains(nodeIndex); @@ -1086,7 +1090,10 @@ void DRMSimulationModel::updateNodalExternalForces( } node.force.external = nodalExternalForce; + totalExternalForcesNorm += node.force.external.norm(); } + + pMesh->totalExternalForcesNorm = totalExternalForcesNorm; } void DRMSimulationModel::updateResidualForces() @@ -1242,6 +1249,13 @@ void DRMSimulationModel::updateNodalMasses() * rotationalSumSk_I3; pMesh->nodes[v].mass.rotationalJ = gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J; + pMesh->nodes[v].mass_6d[DoF::Ux] = pMesh->nodes[v].mass.translational; + pMesh->nodes[v].mass_6d[DoF::Uy] = pMesh->nodes[v].mass.translational; + pMesh->nodes[v].mass_6d[DoF::Uz] = pMesh->nodes[v].mass.translational; + pMesh->nodes[v].mass_6d[DoF::Nx] = pMesh->nodes[v].mass.rotationalJ; + pMesh->nodes[v].mass_6d[DoF::Ny] = pMesh->nodes[v].mass.rotationalI3; + pMesh->nodes[v].mass_6d[DoF::Nr] = pMesh->nodes[v].mass.rotationalI2; + assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk / pMesh->nodes[v].mass.translational < 2); @@ -1295,17 +1309,17 @@ void DRMSimulationModel::updateNodalVelocities() void DRMSimulationModel::updateNodalDisplacements() { - const bool shouldCapDisplacements = mSettings.displacementCap.has_value(); + // 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()); - std::cout << "Displacement capped" << std::endl; - } + // if (shouldCapDisplacements + // && stepDisplacement.getTranslation().norm() > mSettings.displacementCap) { + // stepDisplacement = stepDisplacement + // * (*mSettings.displacementCap + // / stepDisplacement.getTranslation().norm()); + // std::cout << "Displacement capped" << std::endl; + // } node.displacements = node.displacements + stepDisplacement; // if (mSettings.isDebugMode && mSettings.beVerbose && pMesh->getIndex(v) == 40) { // std::cout << "Node " << node.vi << ":" << endl; @@ -1397,7 +1411,7 @@ void DRMSimulationModel::updateNodeNormal( * */ const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0; if (!checkedForMaximumMoment && viHasMoments) { - mSettings.shouldUseTranslationalKineticEnergyThreshold = true; + mSettings.totalTranslationalKineticEnergyThreshold = 1e-8; #ifdef POLYSCOPE_DEFINED std::cout << "Maximum moment reached.The Kinetic energy of the system will " "be used as a convergence criterion" @@ -1830,9 +1844,14 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) // instead of full width. Must have // matching PopItemWidth() below. - ImGui::InputInt("Simulation debug step", - &mSettings.debugModeStep); // set a int variable - mSettings.drawingStep = mSettings.debugModeStep; + static int debugModeStep = 0; + if (ImGui::InputInt("Simulation debug step", + &debugModeStep)) { // set a int variable + if (debugModeStep != 0) { + *mSettings.debugModeStep = debugModeStep; + } + } + // mSettings.drawingStep = mSettings.debugModeStep; ImGui::Checkbox("Enable drawing", &mSettings.shouldDraw); // set a int variable ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); @@ -2024,6 +2043,42 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue // //} // polyscope::show(); } +void DRMSimulationModel::applyKineticDamping(const std::shared_ptr &pJob) +{ + // 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(); + + // 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); + // } + resetVelocities(); + ++numOfDampings; +} SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob, const Settings &settings, @@ -2072,7 +2127,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrgetEigenVertices(), @@ -2090,35 +2145,44 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces) { + externalForce.second = externalForce.second * forceScaleFactor; + } + LinearSimulationModel linearSimulationModel; + SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel + .executeSimulation(pJob); + // simulationResults_fullPatternLinearModel.save(fullResultsFolderPath); + for (auto &externalForce : pJob->nodalExternalForces) { + externalForce.second = externalForce.second / forceScaleFactor; + } + + applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); + shouldTemporarilyDampForces = true; } updateNodalMasses(); - std::unordered_map nodalExternalForces = pJob->nodalExternalForces; - double totalExternalForcesNorm = 0; + // std::unordered_map nodalExternalForces = pJob->nodalExternalForces; + // double totalExternalForcesNorm = 0; // Vector6d sumOfExternalForces(0); - for (auto &nodalForce : nodalExternalForces) { - const double percentageOfExternalLoads = double(externalLoadStep) - / mSettings.desiredGradualExternalLoadsSteps; - nodalForce.second = nodalForce.second * percentageOfExternalLoads; - totalExternalForcesNorm += nodalForce.second.norm(); - // sumOfExternalForces = sumOfExternalForces + nodalForce.second; - } - pMesh->totalExternalForcesNorm = totalExternalForcesNorm; - updateNodalExternalForces(nodalExternalForces, constrainedVertices); - if (!nodalExternalForces.empty()) { + // for (auto &nodalForce : nodalExternalForces) { + // const double percentageOfExternalLoads = double(externalLoadStep) + // / mSettings.desiredGradualExternalLoadsSteps; + // nodalForce.second = nodalForce.second * percentageOfExternalLoads; + // totalExternalForcesNorm += nodalForce.second.norm(); + // // sumOfExternalForces = sumOfExternalForces + nodalForce.second; + // } + updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices); + if (!pJob->nodalExternalForces.empty()) { mSettings.totalResidualForcesNormThreshold - = settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm; + = settings.totalExternalForcesNormPercentageTermination + * pMesh->totalExternalForcesNorm; } else { mSettings.totalResidualForcesNormThreshold = 1e-3; std::cout << "No forces setted default residual forces norm threshold" << std::endl; @@ -2126,9 +2190,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr percentageOfConvergence; // double residualForcesMovingAverageDerivativeMax = 0; - while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { - if ((mSettings.isDebugMode && mCurrentSimulationStep == 50000)) { + while (!settings.maxDRMIterations.has_value() + || mCurrentSimulationStep < settings.maxDRMIterations.value()) { + if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { // std::filesystem::create_directory("./PatternOptimizationNonConv"); // pJob->save("./PatternOptimizationNonConv"); // Dt = mSettings.Dtini; @@ -2221,25 +2286,26 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrcurrentTotalPotentialEnergykN; // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; // timePerNodePerIterationHistor.push_back(timePerNodePerIteration); - if (mSettings.beVerbose && mSettings.isDebugMode - && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + 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.useAverage) { + if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { std::cout << "Best percentage of target (average):" - << 100 * mSettings.averageResidualForcesCriterionThreshold - * totalExternalForcesNorm - / (minTotalResidualForcesNorm / pMesh->VN()) + << 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; } - std::cout << "Best percentage of target:" - << 100 * mSettings.totalExternalForcesNormPercentageTermination - * totalExternalForcesNorm / minTotalResidualForcesNorm - << "%" << std::endl; // SimulationResultsReporter::createPlot("Number of Steps", // "Residual Forces mov aver", // history.residualForcesMovingAverage); @@ -2252,14 +2318,15 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrtotalResidualForcesNorm); } - if (mSettings.shouldCreatePlots && mSettings.isDebugMode - && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() + && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { // SimulationResultsReporter::createPlot("Number of Steps", // "Residual Forces mov aver deriv", // movingAveragesDerivatives_norm); @@ -2295,8 +2362,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr maxDRMIterations*/ { @@ -2320,20 +2387,35 @@ currentSimulationStep > maxDRMIterations*/ // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // << std::endl; const bool fullfillsResidualForcesNormTerminationCriterion - = !mSettings.useAverage - && pMesh->totalResidualForcesNorm / totalExternalForcesNorm + = !mSettings.averageResidualForcesCriterionThreshold.has_value() + && pMesh->totalResidualForcesNorm / pMesh->totalExternalForcesNorm < mSettings.totalExternalForcesNormPercentageTermination; const bool fullfillsAverageResidualForcesNormTerminationCriterion - = mSettings.useAverage - && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm - < mSettings.averageResidualForcesCriterionThreshold; + = mSettings.averageResidualForcesCriterionThreshold.has_value() + && (pMesh->totalResidualForcesNorm / pMesh->VN()) / pMesh->totalExternalForcesNorm + < mSettings.averageResidualForcesCriterionThreshold.value() + && numOfDampings > 2; + if ((fullfillsAverageResidualForcesNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion) + && (pJob->nodalForcedDisplacements.empty() + || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) { + 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 (fullfillsAverageResidualForcesNormTerminationCriterion) { + std::cout << "Converged using average residual forces norm threshold criterion" + << std::endl; + } + } + break; + } // Residual forces norm convergence - if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - || fullfillsAverageResidualForcesNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion) + if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy // && mCurrentSimulationStep > movingAverageSampleSize - && (pJob->nodalForcedDisplacements.empty() - || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) + ) /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 && mCurrentSimulationStep > 1)*/ /*|| @@ -2342,103 +2424,61 @@ mesh->currentTotalPotentialEnergykN*/ /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { const bool fullfillsKineticEnergyTerminationCriterion - = mSettings.shouldUseTranslationalKineticEnergyThreshold + = mSettings.totalTranslationalKineticEnergyThreshold.has_value() && pMesh->currentTotalTranslationalKineticEnergy - < mSettings.totalTranslationalKineticEnergyThreshold + < mSettings.totalTranslationalKineticEnergyThreshold.value() && mCurrentSimulationStep > 20 && numOfDampings > 0; - const bool fullfillsMovingAverageNormTerminationCriterion - = pMesh->residualForcesMovingAverage - < mSettings.residualForcesMovingAverageNormThreshold; + // const bool fullfillsMovingAverageNormTerminationCriterion + // = pMesh->residualForcesMovingAverage + // < mSettings.residualForcesMovingAverageNormThreshold; /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ - const bool fullfillsMovingAverageDerivativeNormTerminationCriterion - = pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4; - const bool shouldTerminate - = fullfillsKineticEnergyTerminationCriterion - // || fullfillsMovingAverageNormTerminationCriterion - // || fullfillsMovingAverageDerivativeNormTerminationCriterion - || fullfillsAverageResidualForcesNormTerminationCriterion - || fullfillsResidualForcesNormTerminationCriterion; - if (shouldTerminate && mCurrentSimulationStep > 100) { + // const bool fullfillsMovingAverageDerivativeNormTerminationCriterion + // = pMesh->residualForcesMovingAverageDerivativeNorm < 1e-4; + const bool shouldTerminate = fullfillsKineticEnergyTerminationCriterion + // || fullfillsMovingAverageNormTerminationCriterion + // || fullfillsMovingAverageDerivativeNormTerminationCriterion + ; + 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) { + if (fullfillsKineticEnergyTerminationCriterion) { std::cout << "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; } } - if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) { - break; - } else { - externalLoadStep++; - std::unordered_map nodalExternalForces - = pJob->nodalExternalForces; - double totalExternalForcesNorm = 0; - for (auto &nodalForce : nodalExternalForces) { - const double percentageOfExternalLoads - = double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps; - nodalForce.second = nodalForce.second * percentageOfExternalLoads; - totalExternalForcesNorm += nodalForce.second.norm(); - } - updateNodalExternalForces(nodalExternalForces, constrainedVertices); + // if (mSettings.desiredGradualExternalLoadsSteps == externalLoadStep) { + break; + // } else { + // externalLoadStep++; + // std::unordered_map nodalExternalForces + // = pJob->nodalExternalForces; + // double totalExternalForcesNorm = 0; + // for (auto &nodalForce : nodalExternalForces) { + // const double percentageOfExternalLoads + // = double(externalLoadStep) / mSettings.desiredGradualExternalLoadsSteps; + // nodalForce.second = nodalForce.second * percentageOfExternalLoads; + // totalExternalForcesNorm += nodalForce.second.norm(); + // } + // updateNodalExternalForces(nodalExternalForces, constrainedVertices); - if (!nodalExternalForces.empty()) { - mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm; - } - if (mSettings.beVerbose) { - std::cout << "totalResidualForcesNormThreshold:" - << mSettings.totalResidualForcesNormThreshold << std::endl; - } - } + // if (!nodalExternalForces.empty()) { + // mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm; + // } + // if (mSettings.beVerbose) { + // std::cout << "totalResidualForcesNormThreshold:" + // << mSettings.totalResidualForcesNormThreshold << std::endl; + // } + // } // } } - 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(); - - // 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); - // } - resetVelocities(); + applyKineticDamping(pJob); Dt *= mSettings.xi; // if (mSettings.isDebugMode) { // std::cout << Dt << std::endl; // } - ++numOfDampings; if (mSettings.shouldCreatePlots) { history.redMarks.push_back(mCurrentSimulationStep); } @@ -2452,11 +2492,12 @@ mesh->currentTotalPotentialEnergykN*/ results.executionTime = std::chrono::duration_cast(endTime - beginTime) .count(); - if (mCurrentSimulationStep == settings.maxDRMIterations && mCurrentSimulationStep != 0) { + if (settings.maxDRMIterations.has_value() && 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; + << settings.maxDRMIterations.value() << "). Breaking simulation" << std::endl; } results.converged = false; } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { @@ -2467,7 +2508,7 @@ mesh->currentTotalPotentialEnergykN*/ } // mesh.printVertexCoordinates(mesh.VN() / 2); #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw && !mSettings.isDebugMode) { + if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) { draw(); } #endif @@ -2568,13 +2609,13 @@ mesh->currentTotalPotentialEnergykN*/ } } - if (!mSettings.isDebugMode && mSettings.shouldCreatePlots) { + if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) { SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); } #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw || mSettings.isDebugMode) { + if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { polyscope::removeCurveNetwork(meshPolyscopeLabel); polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); } @@ -2583,18 +2624,103 @@ mesh->currentTotalPotentialEnergykN*/ return results; } -bool DRMSimulationModel::Settings::save(const filesystem::__cxx11::path &folderPath) const +void DRMSimulationModel::Settings::save(const std::filesystem::path &folderPath) const { - bool returnValue = true; std::filesystem::create_directories(folderPath); nlohmann::json json; - json[jsonLabels.meshFilename] - = std::filesystem::relative(std::filesystem::path(meshFilename), - std::filesystem::path(jsonFilename).parent_path()) - .string(); + + json[jsonLabels.gamma] = gamma; + json[jsonLabels.shouldDraw] = shouldDraw ? "true" : "false"; + json[jsonLabels.shouldCreatePlots] = shouldCreatePlots ? "true" : "false"; + json[jsonLabels.beVerbose] = beVerbose ? "true" : "false"; + json[jsonLabels.Dtini] = Dtini; + json[jsonLabels.xi] = xi; + if (maxDRMIterations.has_value()) { + json[jsonLabels.maxDRMIterations] = maxDRMIterations.value(); + } + if (debugModeStep.has_value()) { + json[jsonLabels.debugModeStep] = debugModeStep.value(); + } + if (totalTranslationalKineticEnergyThreshold.has_value()) { + json[jsonLabels.totalTranslationalKineticEnergyThreshold] + = totalTranslationalKineticEnergyThreshold.value(); + } + + if (averageResidualForcesCriterionThreshold.has_value()) { + json[jsonLabels.averageResidualForcesCriterionThreshold] + = averageResidualForcesCriterionThreshold.value(); + } + + if (averageResidualForcesCriterionThreshold.has_value()) { + json[jsonLabels.averageResidualForcesCriterionThreshold] + = averageResidualForcesCriterionThreshold.value(); + } + if (linearGuessForceScaleFactor.has_value()) { + json[jsonLabels.linearGuessForceScaleFactor] = linearGuessForceScaleFactor.value(); + } + const std::string jsonFilename = "drmSettings.json"; std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename)); jsonFile << json; } -bool DRMSimulationModel::Settings::load(const filesystem::__cxx11::path &filePath) const {} +bool DRMSimulationModel::Settings::load(const std::filesystem::path &jsonFilePath) +{ + if (!std::filesystem::exists(std::filesystem::path(jsonFilePath))) { + std::cerr << "The json file does not exist. Json file provided:" << jsonFilePath.string() + << std::endl; + assert(false); + return false; + } + + if (std::filesystem::path(jsonFilePath).extension() != ".json") { + std::cerr << "A json file is expected as input. The given file has the " + "following extension:" + << std::filesystem::path(jsonFilePath).extension() << std::endl; + assert(false); + return false; + } + + nlohmann::json json; + std::ifstream ifs(jsonFilePath.string()); + ifs >> json; + + if (json.contains(jsonLabels.shouldDraw)) { + shouldDraw = json.at(jsonLabels.shouldDraw) == "true" ? true : false; + } + if (json.contains(jsonLabels.beVerbose)) { + beVerbose = json.at(jsonLabels.beVerbose) == "true" ? true : false; + } + if (json.contains(jsonLabels.shouldCreatePlots)) { + shouldCreatePlots = json.at(jsonLabels.shouldCreatePlots) == "true" ? true : false; + } + if (json.contains(jsonLabels.Dtini)) { + Dtini = json.at(jsonLabels.Dtini); + } + if (json.contains(jsonLabels.xi)) { + xi = json.at(jsonLabels.xi); + } + if (json.contains(jsonLabels.maxDRMIterations)) { + maxDRMIterations = json.at(jsonLabels.maxDRMIterations); + } + if (json.contains(jsonLabels.debugModeStep)) { + debugModeStep = json.at(jsonLabels.debugModeStep); + } + if (json.contains(jsonLabels.totalTranslationalKineticEnergyThreshold)) { + totalTranslationalKineticEnergyThreshold = json.at( + jsonLabels.totalTranslationalKineticEnergyThreshold); + } + if (json.contains(jsonLabels.gamma)) { + gamma = json.at(jsonLabels.gamma); + } + if (json.contains(jsonLabels.averageResidualForcesCriterionThreshold)) { + averageResidualForcesCriterionThreshold = json.at( + jsonLabels.averageResidualForcesCriterionThreshold); + } + + if (json.contains(jsonLabels.linearGuessForceScaleFactor)) { + linearGuessForceScaleFactor = json.at(jsonLabels.linearGuessForceScaleFactor); + } + + return true; +} diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index bb68a07..ccba95a 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -25,32 +25,40 @@ class DRMSimulationModel public: struct Settings { - // bool isDebugMode{false}; - std::optional debugModeStep{0}; bool shouldDraw{false}; bool beVerbose{false}; bool shouldCreatePlots{false}; - // int drawingStep{0}; - // double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; - // double residualForcesMovingAverageNormThreshold{1e-8}; double Dtini{0.1}; double xi{0.9969}; - int maxDRMIterations{0}; - std::optional shouldUseTranslationalKineticEnergyThreshold; - // int gradualForcedDisplacementSteps{50}; - // int desiredGradualExternalLoadsSteps{1}; + int gradualForcedDisplacementSteps{50}; double gamma{0.8}; - std::optional displacementCap; double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; - std::optional averageResidualForcesCriterionThreshold{1e-5}; + std::optional maxDRMIterations; + std::optional debugModeStep; + std::optional totalTranslationalKineticEnergyThreshold; + std::optional averageResidualForcesCriterionThreshold; + std::optional linearGuessForceScaleFactor; Settings() {} - bool save(const std::filesystem::path &folderPath) const; - bool load(const std::filesystem::path &filePath) const; - bool useViscousDamping{false}; - bool useKineticDamping{true}; + void save(const std::filesystem::path &folderPath) const; + bool load(const std::filesystem::path &filePath); struct JsonLabels - {}; + { + const std::string shouldDraw{"shouldDraw"}; + const std::string beVerbose{"beVerbose"}; + const std::string shouldCreatePlots{"shouldCreatePlots"}; + const std::string Dtini{"DtIni"}; + const std::string xi{"xi"}; + const std::string maxDRMIterations{"maxIterations"}; + const std::string debugModeStep{"debugModeStep"}; + const std::string totalTranslationalKineticEnergyThreshold{ + "totalTranslationaKineticEnergyThreshold"}; + const std::string gamma{"gamma"}; + const std::string averageResidualForcesCriterionThreshold{ + "averageResidualForcesThreshold"}; + const std::string linearGuessForceScaleFactor{"linearGuessForceScaleFactor"}; + }; + static JsonLabels jsonLabels; }; private: @@ -219,6 +227,8 @@ private: void updateNodeNr(VertexType &v); + void applyKineticDamping(const std::shared_ptr &pJob); + public: DRMSimulationModel(); SimulationResults executeSimulation(const std::shared_ptr &pJob, @@ -228,6 +238,8 @@ private: static void runUnitTests(); }; +inline DRMSimulationModel::Settings::JsonLabels DRMSimulationModel::Settings::jsonLabels; + template PointType Cross(PointType p1, PointType p2) { return p1 ^ p2; } diff --git a/polymesh.hpp b/polymesh.hpp index ea11c57..3e7a4ad 100755 --- a/polymesh.hpp +++ b/polymesh.hpp @@ -192,10 +192,10 @@ public: // getVertices(); // } vcg::tri::UpdateTopology::VertexEdge(*this); - vcg::tri::UpdateTopology::EdgeEdge(*this); vcg::tri::UpdateTopology::FaceFace(*this); vcg::tri::UpdateTopology::VertexFace(*this); vcg::tri::UpdateTopology::AllocateEdge(*this); + vcg::tri::UpdateTopology::EdgeEdge(*this); // vcg::tri::UpdateTopology::VertexFace(*this); return true; diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index d62a085..6d02176 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -170,7 +170,7 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced // std::filesystem::path(resultsPath).append("patterns.patt")); // } if (numberOfDesiredEdges == -1) { - for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); numberOfEdges++) { + for (size_t numberOfEdges = 2; numberOfEdges <= validEdges.size(); numberOfEdges++) { std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges." << std::endl; @@ -478,6 +478,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15); if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { + statistics.numberOfPatternsViolatingAngleThreshold++; if (debugIsOn /*|| savePlyFiles*/) { if (savePlyFiles) { exportPattern(std::filesystem::path(resultsPath) @@ -493,6 +494,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasNodeWithValenceGreaterThanDesired = patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 6); if (tiledPatternHasNodeWithValenceGreaterThanDesired) { + statistics.numberOfPatternsViolatingValenceThreshold++; if (debugIsOn) { if (savePlyFiles) { auto highValencePath = std::filesystem::path(resultsPath) @@ -648,6 +650,9 @@ void TopologyEnumerator::computeValidPatterns( PatternIO::save(compressedPatternsFilePath, patternSet); } + std::cout << "patternIndex:" << patternIndex << std::endl; + std::cout << "numberOfPatterns:" << numberOfPatterns << std::endl; + if (!validPatternsExist) { std::filesystem::remove_all(validPatternsPath); if (!debugIsOn) { @@ -655,6 +660,7 @@ void TopologyEnumerator::computeValidPatterns( } } } + std::vector TopologyEnumerator::getCoincideEdges( const std::vector &edgeNodeIndices) const { diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp index 7871ce6..0de1444 100755 --- a/topologyenumerator.hpp +++ b/topologyenumerator.hpp @@ -62,6 +62,8 @@ class TopologyEnumerator { size_t numberOfPatternsWithMoreThanASingleCC{0}; size_t numberOfPatternsWithADanglingEdgeOrNode{0}; size_t numberOfPatternsWithArticulationPoints{0}; + size_t numberOfPatternsViolatingAngleThreshold{0}; + size_t numberOfPatternsViolatingValenceThreshold{0}; size_t numberOfValidPatterns{0}; nlohmann::json convertToJson() const { @@ -79,6 +81,8 @@ class TopologyEnumerator { json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints; json["numValidPatterns"] = numberOfValidPatterns; + json["violatingAngle"] = numberOfPatternsViolatingAngleThreshold; + json["violatingValence"] = numberOfPatternsViolatingValenceThreshold; return json; } @@ -141,6 +145,8 @@ class TopologyEnumerator { numberOfPatternsWithADanglingEdgeOrNode = 0; numberOfPatternsWithArticulationPoints = 0; numberOfValidPatterns = 0; + numberOfPatternsViolatingAngleThreshold = 0; + numberOfPatternsViolatingValenceThreshold = 0; } };