diff --git a/CMakeLists.txt.user b/CMakeLists.txt.user new file mode 100644 index 0000000..d9cb88e --- /dev/null +++ b/CMakeLists.txt.user @@ -0,0 +1,355 @@ + + + + + + EnvironmentId + {4ab9621f-3cb4-4d5c-a5a5-8f7ea53a6f0f} + + + ProjectExplorer.Project.ActiveTarget + 0 + + + ProjectExplorer.Project.EditorSettings + + false + false + true + + Cpp + + CppGlobal + + + + QmlJS + + QmlJSGlobal + + + 2 + UTF-8 + false + 4 + false + 80 + true + true + 1 + true + false + 0 + true + true + 0 + 8 + true + 1 + true + true + true + *.md, *.MD, Makefile + false + true + + + + ProjectExplorer.Project.PluginSettings + + + true + true + true + true + true + + + 0 + true + + true + true + Builtin.DefaultTidyAndClazy + 6 + + + + true + + + + + ProjectExplorer.Project.Target.0 + + Desktop + Clang + Clang + {08d5b07f-8010-48d8-a202-51411d8039c8} + 0 + 0 + 0 + + -GUnix Makefiles +-DCMAKE_BUILD_TYPE:String=Debug +-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable} +-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX} +-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C} +-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx} + /home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-Debug + + + + all + + true + Build + CMakeProjectManager.MakeStep + + 1 + Build + Build + ProjectExplorer.BuildSteps.Build + + + + + clean + + true + Build + CMakeProjectManager.MakeStep + + 1 + Clean + Clean + ProjectExplorer.BuildSteps.Clean + + 2 + false + + + Debug + CMakeProjectManager.CMakeBuildConfiguration + + + -GUnix Makefiles +-DCMAKE_BUILD_TYPE:String=Release +-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable} +-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX} +-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C} +-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx} + /home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-Release + + + + all + + true + CMakeProjectManager.MakeStep + + 1 + Build + Build + ProjectExplorer.BuildSteps.Build + + + + + clean + + true + CMakeProjectManager.MakeStep + + 1 + Clean + Clean + ProjectExplorer.BuildSteps.Clean + + 2 + false + + + Release + CMakeProjectManager.CMakeBuildConfiguration + + + -GUnix Makefiles +-DCMAKE_BUILD_TYPE:String=RelWithDebInfo +-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable} +-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX} +-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C} +-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx} + /home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-RelWithDebInfo + + + + all + + true + CMakeProjectManager.MakeStep + + 1 + Build + Build + ProjectExplorer.BuildSteps.Build + + + + + clean + + true + CMakeProjectManager.MakeStep + + 1 + Clean + Clean + ProjectExplorer.BuildSteps.Clean + + 2 + false + + + Release with Debug Information + CMakeProjectManager.CMakeBuildConfiguration + + + -GUnix Makefiles +-DCMAKE_BUILD_TYPE:String=MinSizeRel +-DQT_QMAKE_EXECUTABLE:STRING=%{Qt:qmakeExecutable} +-DCMAKE_PREFIX_PATH:STRING=%{Qt:QT_INSTALL_PREFIX} +-DCMAKE_C_COMPILER:STRING=%{Compiler:Executable:C} +-DCMAKE_CXX_COMPILER:STRING=%{Compiler:Executable:Cxx} + /home/iason/Coding/Libraries/MySources_viscousDampingBranch/../build-MySources_viscousDampingBranch-Clang-MinSizeRel + + + + all + + true + CMakeProjectManager.MakeStep + + 1 + Build + Build + ProjectExplorer.BuildSteps.Build + + + + + clean + + true + CMakeProjectManager.MakeStep + + 1 + Clean + Clean + ProjectExplorer.BuildSteps.Clean + + 2 + false + + + Minimum Size Release + CMakeProjectManager.CMakeBuildConfiguration + + 4 + + + 0 + Deploy + Deploy + ProjectExplorer.BuildSteps.Deploy + + 1 + + false + ProjectExplorer.DefaultDeployConfiguration + + 1 + + dwarf + + cpu-cycles + + + 250 + + -e + cpu-cycles + --call-graph + dwarf,4096 + -F + 250 + + -F + true + 4096 + false + false + 1000 + + true + + false + false + false + false + true + 0.01 + 10 + true + kcachegrind + 1 + 25 + + 1 + true + false + true + valgrind + + 0 + 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14 + + + 2 + + ProjectExplorer.CustomExecutableRunConfiguration + + false + true + false + true + + 1 + + + + ProjectExplorer.Project.TargetCount + 1 + + + ProjectExplorer.Project.Updater.FileVersion + 22 + + + Version + 22 + + diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 4dd80cb..b2d6645 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); @@ -1215,6 +1228,14 @@ void DRMSimulationModel::updateNodalMasses() * rotationalSumSk_I3; pMesh->nodes[v].mass.rotationalJ = gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J; + //fill 6d mass vector + 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); @@ -1234,17 +1255,7 @@ void DRMSimulationModel::updateNodalAccelerations() for (VertexType &v : pMesh->vert) { Node &node = pMesh->nodes[v]; const VertexIndex vi = pMesh->getIndex(v); - for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { - if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.translational; - } else if (dofi == DoF::Nx) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalJ; - } else if (dofi == DoF::Ny) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI3; - } else if (dofi == DoF::Nr) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2; - } - } + node.acceleration = node.force.residual / node.mass_6d; #ifdef POLYSCOPE_DEFINED if (std::isnan(node.acceleration.norm())) { std::cout << "acceleration " << vi << ":" << node.acceleration.toString() << std::endl; @@ -1271,7 +1282,16 @@ void DRMSimulationModel::updateNodalVelocities() std::cout << "Velocity " << vi << ":" << node.velocity.toString() << std::endl; } #endif - node.velocity = node.velocity + node.acceleration * Dt; + if (mSettings.useViscousDamping) { + const Vector6d massOverDt = node.mass_6d / Dt; + const Vector6d denominator = massOverDt + Vector6d(viscuousDampingConstant / 2); + const Vector6d firstTermNominator = massOverDt - Vector6d(viscuousDampingConstant / 2); + const Vector6d firstTerm = node.velocity * firstTermNominator / denominator; + const Vector6d secondTerm = node.force.residual / denominator; + node.velocity = firstTerm + secondTerm; + } else { + node.velocity = node.velocity + node.acceleration * Dt; + } } updateKineticEnergy(); } @@ -1328,7 +1348,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 +1529,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 +1636,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 +1670,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 +1959,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 +2069,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 +2082,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 +2145,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 +2233,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 +2322,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 +2349,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 +2369,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 +2405,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 +2436,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..e2afd98 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -45,8 +45,9 @@ public: double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; bool useAverage{false}; - double averageResidualForcesCriterionThreshold{1e-3}; + double averageResidualForcesCriterionThreshold{1e-5}; Settings() {} + bool useViscousDamping{false}; }; private: @@ -61,6 +62,10 @@ private: std::vector plotYValues; size_t numOfDampings{0}; int externalLoadStep{1}; + const double viscuousDampingConstant{0.01}; + std::vector isVertexConstrained; + std::vector isRigidSupport; + double minTotalResidualForcesNorm{std::numeric_limits::max()}; const std::string meshPolyscopeLabel{"Simulation mesh"}; std::unique_ptr pMesh; @@ -69,7 +74,6 @@ private: // Eigen::Tensor theta3Derivatives; // std::unordered_map theta3Derivatives; bool shouldApplyInitialDistortion = false; - std::unordered_set rigidSupports; void reset(); void updateNodalInternalForces( @@ -156,10 +160,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..234a535 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -148,6 +148,7 @@ struct Node { }; Mass mass; + Vector6d mass_6d; VertexIndex vi; CoordType initialLocation; CoordType initialNormal; @@ -157,10 +158,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; diff --git a/utilities.hpp b/utilities.hpp index 863768a..1d2dd6e 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -57,6 +57,15 @@ struct Vector6d : public std::array { return result; } + Vector6d operator/(const Vector6d &v) const + { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) / v[i]; + } + return result; + } + Vector6d operator+(const Vector6d &v) const { Vector6d result; for (size_t i = 0; i < 6; i++) {