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

Refactoring

This commit is contained in:
iasonmanolas 2021-07-12 23:53:04 +03:00
parent 2509b0a795
commit 5ebf354dcf
4 changed files with 156 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

236
drmsimulationmodel.cpp
View File

@ -640,6 +640,9 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
? node.derivativeOfNormal[dui.dofi]
: VectorType(0, 0, 0);
derivativeTheta3_dofi = derivativeG1 * n + g1 * derivativeNormal;
if (std::isnan(derivativeTheta3_dofi)) {
std::cerr << "nan derivative theta3 rigid" << std::endl;
}
return derivativeTheta3_dofi;
}
EdgeType &refElem = *node.referenceElement;
@ -718,6 +721,9 @@ double DRMSimulationModel::computeDerivativeTheta3(const EdgeType &e,
derivativeTheta3_dofi = derivativeF3 * n + f3 * derivativeOfNormal;
}
}
if (std::isnan(derivativeTheta3_dofi)) {
std::cerr << "nan derivative theta3" << std::endl;
}
return derivativeTheta3_dofi;
}
@ -794,9 +800,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
pMesh->EN(), std::vector<std::pair<int, Vector6d>>(4, {-1, Vector6d()}));
// omp_lock_t writelock;
// omp_init_lock(&writelock);
#ifdef ENABLE_OPENMP
#pragma omp parallel for schedule(static) num_threads(6)
#endif
//#ifdef ENABLE_OPENMP
//#pragma omp parallel for schedule(static) num_threads(8)
//#endif
for (int ei = 0; ei < pMesh->EN(); ei++) {
const EdgeType &e = pMesh->edge[ei];
const SimulationMesh::VertexType &ev_j = *e.cV(0);
@ -846,6 +852,12 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double e_kDeriv = computeDerivativeElementLength(e, dui);
const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(axialForce_dofi)) {
std::cerr << "nan in axial" << evi << std::endl;
}
#endif
// Torsional force computation
const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi);
const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi);
@ -853,6 +865,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv;
const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff
* element.rigidity.torsional;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(torsionalForce_dofi)) {
std::cerr << "nan in torsional" << evi << std::endl;
}
#endif
// First bending force computation
////theta2_j derivative
@ -970,6 +987,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
Node::Forces &force = pMesh->nodes[vi].force;
force.internal = force.internal + internalForcePair.second;
force.residual = force.residual + (internalForcePair.second * -1);
#ifdef POLYSCOPE_DEFINED
if (std::isnan(internalForcePair.second.norm())) {
std::cerr << "nan on edge" << ei << std::endl;
}
#endif
}
}
double totalResidualForcesNorm = 0;
@ -982,6 +1004,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
// const double residualForceNorm = nodeResidualForce.getTranslation().norm();
const bool shouldTrigger = std::isnan(residualForceNorm);
totalResidualForcesNorm += residualForceNorm;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(force.residual.norm())) {
std::cout << "residual " << vi << ":" << force.residual.toString() << std::endl;
}
#endif
}
pMesh->previousTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
@ -1218,6 +1245,11 @@ void DRMSimulationModel::updateNodalAccelerations()
node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2;
}
}
#ifdef POLYSCOPE_DEFINED
if (std::isnan(node.acceleration.norm())) {
std::cout << "acceleration " << vi << ":" << node.acceleration.toString() << std::endl;
}
#endif
// if (vi == 10) {
// std::cout << "Acceleration:" << node.acceleration[0] << " " << node.acceleration[1]
// << " " << node.acceleration[2] << std::endl;
@ -1234,10 +1266,12 @@ void DRMSimulationModel::updateNodalVelocities()
for (VertexType &v : pMesh->vert) {
const VertexIndex vi = pMesh->getIndex(v);
Node &node = pMesh->nodes[v];
node.velocity = node.velocity + node.acceleration * Dt;
#ifdef POLYSCOPE_DEFINED
if (std::isnan(node.velocity.norm())) {
std::cout << "Velocity:" << node.velocity.toString() << std::endl;
std::cout << "Velocity " << vi << ":" << node.velocity.toString() << std::endl;
}
#endif
node.velocity = node.velocity + node.acceleration * Dt;
}
updateKineticEnergy();
}
@ -1609,6 +1643,11 @@ void DRMSimulationModel::printCurrentState() const
{
std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl;
std::cout << "Residual forces norm: " << pMesh->totalResidualForcesNorm << std::endl;
std::cout << "Average Residual forces norm/extForcesNorm: "
<< pMesh->totalResidualForcesNorm / pMesh->VN() / pMesh->totalExternalForcesNorm
<< std::endl;
std::cout << "Moving average residual forces:" << pMesh->residualForcesMovingAverage
<< 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::microseconds>(
@ -1758,8 +1797,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addEdgeScalarQuantity("I2", I2);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addEdgeScalarQuantity("I3", I3);
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3);
// Specify the callback
static bool calledOnce = false;
@ -1798,9 +1836,8 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
.append(std::to_string(mCurrentSimulationStep) + ".png")
.string(),
false);
} else {
polyscope::show();
}
polyscope::show();
}
#endif
@ -1873,6 +1910,9 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue
Eigen::Quaternion<double> q_normal;
q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen);
Eigen::Quaternion<double> q_nr = q_f1.inverse() * q_normal.inverse() * q;
q_nr.w() = q_nr.w() >= 1 ? 1 : q_nr.w();
q_nr.w() = q_nr.w() <= -1 ? -1 : q_nr.w();
const double nr_0To2pi_pos = 2 * std::acos(q_nr.w());
// const double nr_0To2pi_groundOfTruth = 2 * std::acos(solutionGuess.debug_q_nr[vi].w());
const double nr_0To2pi = nr_0To2pi_pos <= M_PI ? nr_0To2pi_pos : (nr_0To2pi_pos - 2 * M_PI);
@ -1880,6 +1920,16 @@ void DRMSimulationModel::applySolutionGuess(const SimulationResults &solutionGue
q_nr.y() * sin(nr_0To2pi_pos / 2),
q_nr.z() * sin(nr_0To2pi_pos / 2));
/*deformedNormal_debug =*/deformedNormal_debug.normalize();
#ifdef POLYSCOPE_DEFINED
if (std::isnan(deformedNormal_debug.norm())) {
std::cerr << "nr_0To2pi_pos:" << nr_0To2pi_pos << std::endl;
std::cerr << "q_nrx:" << q_nr.x() << std::endl;
std::cerr << "q_nry:" << q_nr.y() << std::endl;
std::cerr << "q_nrz:" << q_nr.z() << std::endl;
std::cerr << "q_nrw:" << q_nr.w() << std::endl;
std::cerr << "nan deformedNormal in guess" << std::endl;
}
#endif
const double nr = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? nr_0To2pi : -nr_0To2pi;
if (!pJob->constrainedVertices.contains(vi)
|| !pJob->constrainedVertices.at(vi).contains(5)) {
@ -2006,6 +2056,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
polyscope::registerCurveNetwork(meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
// registerWorldAxes();
}
#endif
@ -2039,15 +2092,11 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
totalExternalForcesNorm += nodalForce.second.norm();
// sumOfExternalForces = sumOfExternalForces + nodalForce.second;
}
pMesh->totalExternalForcesNorm = totalExternalForcesNorm;
updateNodalExternalForces(nodalExternalForces, constrainedVertices);
double averageExternalForcesNorm =
// sumOfExternalForces.norm() / pMesh->VN();
totalExternalForcesNorm / pMesh->VN();
if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold
= std::min(settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm,
1e-3);
= settings.totalExternalForcesNormPercentageTermination * totalExternalForcesNorm;
} else {
mSettings.totalResidualForcesNormThreshold = 1e-3;
std::cout << "No forces setted default residual forces norm threshold" << std::endl;
@ -2055,9 +2104,13 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl;
if (mSettings.useAverage) {
std::cout << "average/extNorm threshold:"
<< mSettings.averageResidualForcesCriterionThreshold << std::endl;
}
}
// const size_t movingAverageSampleSize = 200;
// std::queue<double> residualForcesMovingAverageHistorySample;
const size_t movingAverageSampleSize = 200;
std::queue<double> residualForcesMovingAverageHistorySample;
std::vector<double> percentageOfConvergence;
// double residualForcesMovingAverageDerivativeMax = 0;
@ -2107,37 +2160,40 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
break;
}
//normalized sum of displacements
// double sumOfDisplacementsNorm = 0;
// for (size_t vi = 0; vi < pMesh->VN(); vi++) {
// sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
// }
// sumOfDisplacementsNorm /= pMesh->bbox.Diag();
// pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm;
//moving average
// // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm
// // + mCurrentSimulationStep
// // * pMesh->residualForcesMovingAverage)
// // / (1 + mCurrentSimulationStep);
// pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
// / movingAverageSampleSize;
// residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
// if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
// const double firstElementValue = residualForcesMovingAverageHistorySample.front();
// pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
// residualForcesMovingAverageHistorySample.pop();
pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm
/ movingAverageSampleSize;
residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage);
if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) {
const double firstElementValue = residualForcesMovingAverageHistorySample.front();
pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize;
residualForcesMovingAverageHistorySample.pop();
// pMesh->residualForcesMovingAverageDerivativeNorm
// = std::abs((residualForcesMovingAverageHistorySample.back()
// - residualForcesMovingAverageHistorySample.front()))
// / (movingAverageSampleSize);
// residualForcesMovingAverageDerivativeMax
// = std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
// residualForcesMovingAverageDerivativeMax);
// pMesh->residualForcesMovingAverageDerivativeNorm
// /= residualForcesMovingAverageDerivativeMax;
// // std::cout << "Normalized derivative:"
// // << residualForcesMovingAverageDerivativeNorm
// // << std::endl;
// }
// pMesh->residualForcesMovingAverageDerivativeNorm
// = std::abs((residualForcesMovingAverageHistorySample.back()
// - residualForcesMovingAverageHistorySample.front()))
// / (movingAverageSampleSize);
// residualForcesMovingAverageDerivativeMax
// = std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
// residualForcesMovingAverageDerivativeMax);
// pMesh->residualForcesMovingAverageDerivativeNorm
// /= residualForcesMovingAverageDerivativeMax;
// // std::cout << "Normalized derivative:"
// // << residualForcesMovingAverageDerivativeNorm
// // << std::endl;
}
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
@ -2151,19 +2207,20 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// const double executionTime_min
// = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count();
// std::cout << "Execution time(min):" << executionTime_min << std::endl;
if (mSettings.useAverageResidualForcesNorm) {
const double percOfConv = averageExternalForcesNorm
/ pMesh->averageResidualForcesNorm;
std::cout << "Percentage of target:" << percOfConv << "%" << std::endl;
} else {
std::cout << "Percentage of target:"
<< 100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm
if (mSettings.useAverage) {
std::cout << "Percentage of target (average):"
<< 100 * mSettings.averageResidualForcesCriterionThreshold
* totalExternalForcesNorm
/ (pMesh->totalResidualForcesNorm / pMesh->VN())
<< "%" << std::endl;
}
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// movingAverages);
std::cout << "Percentage of target:"
<< 100 * mSettings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm / pMesh->totalResidualForcesNorm
<< "%" << std::endl;
SimulationResultsReporter::createPlot("Number of Steps",
"Residual Forces mov aver",
history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
@ -2184,14 +2241,14 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// 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);
"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,
@ -2202,9 +2259,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// SimulationResultsReporter::createPlot("Number of Iterations",
// "Sum of normalized displacement norms",
// history.sumOfNormalizedDisplacementNorms /*,
// std::filesystem::path("./")
// .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
// .string()*/
// std::filesystem::path("./")
// .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
// .string()*/
// ,
// {},
// history.redMarks);
@ -2257,12 +2314,13 @@ mesh->currentTotalPotentialEnergykN*/
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 0;
const bool fullfillsAverageResidualForcesNormCriterion
= mSettings.useAverageResidualForcesNorm
&& pMesh->averageResidualForcesNorm / averageExternalForcesNorm < 1e-2;
const bool fullfillsResidualForcesNormTerminationCriterion
= mSettings.useResidualForcesNorm
&& pMesh->totalResidualForcesNorm < mSettings.totalResidualForcesNormThreshold;
= pMesh->totalResidualForcesNorm / totalExternalForcesNorm
< mSettings.totalExternalForcesNormPercentageTermination;
const bool fullfillsAverageResidualForcesNormTerminationCriterion
= mSettings.useAverage
&& (pMesh->totalResidualForcesNorm / pMesh->VN()) / totalExternalForcesNorm
< mSettings.averageResidualForcesCriterionThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold;
@ -2273,10 +2331,10 @@ mesh->currentTotalPotentialEnergykN*/
= fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion
// || fullfillsMovingAverageDerivativeNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion
|| fullfillsAverageResidualForcesNormCriterion;
|| fullfillsAverageResidualForcesNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion;
if (shouldTerminate) {
if (mSettings.beVerbose && !mSettings.isDebugMode) {
if (mSettings.beVerbose /*&& !mSettings.isDebugMode*/) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
if (fullfillsResidualForcesNormTerminationCriterion) {
@ -2308,7 +2366,6 @@ mesh->currentTotalPotentialEnergykN*/
}
updateNodalExternalForces(nodalExternalForces, constrainedVertices);
averageExternalForcesNorm = totalExternalForcesNorm / pMesh->VN();
if (!nodalExternalForces.empty()) {
mSettings.totalResidualForcesNormThreshold = 1e-2 * totalExternalForcesNorm;
}
@ -2320,25 +2377,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 * 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;
@ -2491,5 +2548,12 @@ mesh->currentTotalPotentialEnergykN*/
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}
return results;
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.isDebugMode) {
polyscope::removeCurveNetwork(meshPolyscopeLabel);
polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel);
}
#endif
return results;
}

6
drmsimulationmodel.hpp
View File

@ -38,14 +38,14 @@ public:
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{100};
int gradualForcedDisplacementSteps{50};
int desiredGradualExternalLoadsSteps{1};
double gamma{0.8};
std::optional<double> displacementCap;
double totalResidualForcesNormThreshold{1e-3};
double totalExternalForcesNormPercentageTermination{1e-3};
bool useResidualForcesNorm{true};
bool useAverageResidualForcesNorm{false};
bool useAverage{false};
double averageResidualForcesCriterionThreshold{1e-3};
Settings() {}
};

4
simulation_structs.hpp
View File

@ -53,8 +53,8 @@ struct SimulationHistory {
kineticEnergy.push_back(std::log10(mesh.currentTotalKineticEnergy));
// potentialEnergy.push_back(mesh.totalPotentialEnergykN);
logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage);
sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms);
residualForcesMovingAverage.push_back(std::log10(mesh.residualForcesMovingAverage));
sumOfNormalizedDisplacementNorms.push_back(std::log10(mesh.sumOfNormalizedDisplacementNorms));
// residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
}

1
simulationmesh.hpp
View File

@ -43,6 +43,7 @@ public:
double currentTotalKineticEnergy{0};
double currentTotalTranslationalKineticEnergy{0};
double totalResidualForcesNorm{0};
double totalExternalForcesNorm{0};
double averageResidualForcesNorm{0};
double currentTotalPotentialEnergykN{0};
double previousTotalPotentialEnergykN{0};