manolas
/
MySources
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merged master

This commit is contained in:
iasonmanolas 2021-07-21 17:59:46 +03:00
parent f20a7fbe55 846289448b
commit af4449a9dd
2 changed files with 32 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
drmsimulationmodel.cpp
View File

@ -222,6 +222,7 @@ void DRMSimulationModel::reset()
mCurrentSimulationStep = 0; mCurrentSimulationStep = 0;
history.clear(); history.clear();
constrainedVertices.clear(); constrainedVertices.clear();
rigidSupports.clear();
pMesh.reset(); pMesh.reset();
plotYValues.clear(); plotYValues.clear();
plotHandle.reset(); plotHandle.reset();
@ -692,8 +693,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
derivativeF1 = derivativeT1_k - (n * (derivativeT1_k * n)); derivativeF1 = derivativeT1_k - (n * (derivativeT1_k * n));
const double f1Norm = f1.Norm(); const double f1Norm = f1.Norm();
const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1) const double derivativeF1Norm = f1 * derivativeF1 / f1Norm;
/ (f1Norm * f1Norm * f1Norm) const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm)
+ derivativeF1 / f1Norm; + derivativeF1 / f1Norm;
derivativeF2 = derivativeF1Normalized * cosRotationAngle derivativeF2 = derivativeF1Normalized * cosRotationAngle
@ -713,8 +714,8 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
const VectorType &derivativeOfNormal = node.derivativeOfNormal[dofi]; const VectorType &derivativeOfNormal = node.derivativeOfNormal[dofi];
derivativeF1 = -(n * (t1_k * derivativeOfNormal) + derivativeOfNormal * (t1_k * n)); derivativeF1 = -(n * (t1_k * derivativeOfNormal) + derivativeOfNormal * (t1_k * n));
const double f1Norm = f1.Norm(); const double f1Norm = f1.Norm();
const VectorType derivativeF1Normalized = -f1 * (f1 * derivativeF1) const double derivativeF1Norm = f1 * derivativeF1 / f1Norm;
/ (f1Norm * f1Norm * f1Norm) const VectorType derivativeF1Normalized = -f1 * derivativeF1Norm / (f1Norm * f1Norm)
+ derivativeF1 / f1Norm; + derivativeF1 / f1Norm;
derivativeF2 = derivativeF1Normalized * cosRotationAngle + derivativeF2 = derivativeF1Normalized * cosRotationAngle +
@ -805,9 +806,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
pMesh->EN(), std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()})); pMesh->EN(), std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
// omp_lock_t writelock; // omp_lock_t writelock;
// omp_init_lock(&writelock); // omp_init_lock(&writelock);
#ifdef ENABLE_OPENMP //#ifdef ENABLE_OPENMP
#pragma omp parallel for schedule(static) num_threads(5) //#pragma omp parallel for schedule(static) num_threads(8)
#endif //#endif
for (int ei = 0; ei < pMesh->EN(); ei++) { for (int ei = 0; ei < pMesh->EN(); ei++) {
const EdgeType &e = pMesh->edge[ei]; const EdgeType &e = pMesh->edge[ei];
const SimulationMesh::VertexType &ev_j = *e.cV(0); const SimulationMesh::VertexType &ev_j = *e.cV(0);
@ -1007,7 +1008,6 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
sumOfResidualForces = sumOfResidualForces + nodeResidualForce; sumOfResidualForces = sumOfResidualForces + nodeResidualForce;
const double residualForceNorm = nodeResidualForce.norm(); const double residualForceNorm = nodeResidualForce.norm();
// const double residualForceNorm = nodeResidualForce.getTranslation().norm(); // const double residualForceNorm = nodeResidualForce.getTranslation().norm();
const bool shouldTrigger = std::isnan(residualForceNorm);
totalResidualForcesNorm += residualForceNorm; totalResidualForcesNorm += residualForceNorm;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (std::isnan(force.residual.norm())) { if (std::isnan(force.residual.norm())) {
@ -1192,10 +1192,10 @@ void DRMSimulationModel::updateNodalMasses()
if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) { if (shouldTemporarilyDampForces && mCurrentSimulationStep < untilStep) {
gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep); gamma *= 1e6 * (1 - static_cast<double>(mCurrentSimulationStep) / untilStep);
} }
if (mCurrentSimulationStep == static_cast<size_t>(1.2 * untilStep) // if (mCurrentSimulationStep == static_cast<size_t>(1.4 * untilStep)
&& shouldTemporarilyDampForces) { // && shouldTemporarilyDampForces) {
Dt = mSettings.Dtini; // Dt = mSettings.Dtini;
} // }
for (VertexType &v : pMesh->vert) { for (VertexType &v : pMesh->vert) {
const size_t vi = pMesh->getIndex(v); const size_t vi = pMesh->getIndex(v);
double translationalSumSk = 0; double translationalSumSk = 0;
@ -1777,7 +1777,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0]; accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0];
} }
meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false); 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("Internal force", internalForces)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false); meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true); meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true);
@ -2089,7 +2089,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh); vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
computeRigidSupports(); computeRigidSupports();
isVertexConstrained.resize(pMesh->VN(), false); isVertexConstrained.resize(pMesh->VN(), false);
for (auto fixedVertex : pJob->constrainedVertices) { for (auto fixedVertex : constrainedVertices) {
isVertexConstrained[fixedVertex.first] = true; isVertexConstrained[fixedVertex.first] = true;
} }
@ -2285,19 +2285,19 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv", // "Residual Forces mov aver deriv",
// movingAveragesDerivatives_norm); // movingAveragesDerivatives_norm);
// SimulationResultsReporter::createPlot("Number of Steps", SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver", "Residual Forces mov aver",
// history.residualForcesMovingAverage); history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Residual Forces", // "Log of Residual Forces",
// history.logResidualForces, // history.logResidualForces,
// {}, // {},
// history.redMarks); // history.redMarks);
SimulationResultsReporter::createPlot("Number of Steps", // SimulationResultsReporter::createPlot("Number of Steps",
"Log of Kinetic energy", // "Log of Kinetic energy",
history.kineticEnergy, // history.kineticEnergy,
{}, // {},
history.redMarks); // history.redMarks);
// SimulationResultsReporter reporter; // SimulationResultsReporter reporter;
// reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel()); // reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel());
// SimulationResultsReporter::createPlot("Number of Iterations", // SimulationResultsReporter::createPlot("Number of Iterations",
@ -2349,21 +2349,14 @@ currentSimulationStep > maxDRMIterations*/
= mSettings.useAverage = mSettings.useAverage
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm && (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold; < mSettings.averageResidualForcesCriterionThreshold;
if (fullfillsResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion) {
if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
}
break;
}
// Residual forces norm convergence
if ((mSettings.useKineticDamping)
&& (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
// && mCurrentSimulationStep > movingAverageSampleSize // Residual forces norm convergence
&& (pJob->nodalForcedDisplacements.empty() if (((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)) || fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion)
// && mCurrentSimulationStep > movingAverageSampleSize
&& (pJob->nodalForcedDisplacements.empty()
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps))
/* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 /* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1
&& mCurrentSimulationStep > 1)*/ && mCurrentSimulationStep > 1)*/
/*|| /*||
@ -2388,7 +2381,7 @@ mesh->currentTotalPotentialEnergykN*/
// || fullfillsMovingAverageDerivativeNormTerminationCriterion // || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsAverageResidualForcesNormTerminationCriterion || fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion; || fullfillsResidualForcesNormTerminationCriterion;
if (shouldTerminate) { if (shouldTerminate && mCurrentSimulationStep > 100) {
if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) { if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
std::cout << "Simulation converged." << std::endl; std::cout << "Simulation converged." << std::endl;
printCurrentState(); printCurrentState();

4
drmsimulationmodel.hpp
View File

@ -161,9 +161,9 @@ private:
void updateElementalFrames(); 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;
static double computeDerivativeOfNorm(const VectorType &x, static double computeDerivativeOfNorm(const VectorType &x,
const VectorType &derivativeOfX); const VectorType &derivativeOfX);