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

Refactoring

This commit is contained in:
iasonmanolas 2022-05-06 16:16:43 +03:00
parent 436ece0d88
commit e9707e2cfb
2 changed files with 417 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

690
drmsimulationmodel.cpp
View File

@ -263,10 +263,10 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
PolyscopeInterface::init(); PolyscopeInterface::init();
polyscope::registerCurveNetwork(meshPolyscopeLabel, polyscope::registerCurveNetwork(meshPolyscopeLabel + "_" + pMesh->getLabel(),
pMesh->getEigenVertices(), pMesh->getEigenVertices(),
pMesh->getEigenEdges()); pMesh->getEigenEdges());
polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel(),
pMesh->getEigenVertices(), pMesh->getEigenVertices(),
pMesh->getEigenEdges()) pMesh->getEigenEdges())
->setRadius(0.002); ->setRadius(0.002);
@ -274,9 +274,9 @@ void DRMSimulationModel::reset(const std::shared_ptr<SimulationJob> &pJob, const
} }
#endif #endif
if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { // if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true; // shouldApplyInitialDistortion = true;
} // }
updateElementalFrames(); updateElementalFrames();
} }
@ -584,8 +584,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v)
const double &nR = node.nR; // TODO: This makes the function non-const. const double &nR = node.nR; // TODO: This makes the function non-const.
// Should be refactored in the future // Should be refactored in the future
double theta3; // const bool shouldBreak = mCurrentSimulationStep == 12970;
const bool shouldBreak = mCurrentSimulationStep == 12970;
if (&e == node.referenceElement) { if (&e == node.referenceElement) {
// Use nR as theta3 only for the first star edge // Use nR as theta3 only for the first star edge
return nR; return nR;
@ -631,7 +630,7 @@ double DRMSimulationModel::computeTheta3(const EdgeType &e, const VertexType &v)
element.cosRotationAngle_jplus1 = cosRotationAngle; element.cosRotationAngle_jplus1 = cosRotationAngle;
element.sinRotationAngle_jplus1 = sinRotationAngle; element.sinRotationAngle_jplus1 = sinRotationAngle;
} }
theta3 = f3.dot(n); const double theta3 = f3.dot(n);
return theta3; return theta3;
} }
@ -845,9 +844,18 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices) const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices)
{ {
// std::vector<Vector6d> internalForcesParallel(mesh->vert.size()); // std::vector<Vector6d> internalForcesParallel(mesh->vert.size());
std::for_each(pMesh->nodes._handle->data.begin(),
pMesh->nodes._handle->data.end(),
[](Node &node) {
node.force.internalAxial = 0;
node.force.internalTorsion = 0;
node.force.internalFirstBending = 0;
node.force.internalSecondBending = 0;
});
std::vector<std::vector<std::pair<int, Vector6d>>> internalForcesContributionsFromEachEdge( std::vector<std::vector<std::pair<int, Vector6d>>> internalForcesContributionsFromEachEdge(
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
@ -909,6 +917,7 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double e_k = element.length - element.initialLength; const double e_k = element.length - element.initialLength;
const double e_kDeriv = computeDerivativeElementLength(e, dui); const double e_kDeriv = computeDerivativeElementLength(e, dui);
const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial; const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial;
pMesh->nodes[edgeNode.vi].force.internalAxial[dofi] += axialForce_dofi;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (std::isnan(axialForce_dofi)) { if (std::isnan(axialForce_dofi)) {
@ -923,9 +932,11 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv;
const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff
* element.rigidity.torsional; * element.rigidity.torsional;
pMesh->nodes[edgeNode.vi].force.internalTorsion[dofi] += torsionalForce_dofi;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (std::isnan(torsionalForce_dofi)) { if (std::isnan(torsionalForce_dofi)) {
std::cerr << "nan in torsional" << evi << std::endl; std::cerr << "nan in torsional" << evi << std::endl;
std::terminate();
} }
#endif #endif
@ -954,6 +965,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
+ firstBendingForce_inBracketsTerm_3; + firstBendingForce_inBracketsTerm_3;
const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm
* element.rigidity.firstBending; * element.rigidity.firstBending;
pMesh->nodes[edgeNode.vi].force.internalFirstBending[dofi]
+= firstBendingForce_dofi;
// Second bending force computation // Second bending force computation
////theta2_j derivative ////theta2_j derivative
@ -982,6 +995,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
+ secondBendingForce_inBracketsTerm_3; + secondBendingForce_inBracketsTerm_3;
double secondBendingForce_dofi = secondBendingForce_inBracketsTerm double secondBendingForce_dofi = secondBendingForce_inBracketsTerm
* element.rigidity.secondBending; * element.rigidity.secondBending;
pMesh->nodes[edgeNode.vi].force.internalSecondBending[dofi]
+= secondBendingForce_dofi;
internalForcesContributionFromThisEdge[evi].second[dofi] internalForcesContributionFromThisEdge[evi].second[dofi]
= axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi
+ torsionalForce_dofi; + torsionalForce_dofi;
@ -1021,6 +1036,8 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
* element.rigidity.secondBending; * element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi + 2].second[dofi] internalForcesContributionFromThisEdge[evi + 2].second[dofi]
= secondBendingForce_dofi; = secondBendingForce_dofi;
pMesh->nodes[refElemOtherVertexNode.vi].force.internalSecondBending[dofi]
+= secondBendingForce_dofi;
} }
} }
} }
@ -1028,18 +1045,14 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge; internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge;
}); });
//#pragma omp parallel for schedule(static) num_threads(8) std::for_each(pMesh->nodes._handle->data.begin(),
pMesh->nodes._handle->data.end(),
[](Node &node) {
Node::Forces &force = node.force;
force.residual = force.external;
force.internal = 0;
});
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
Node::Forces &force = pMesh->nodes[vi].force;
// Vector6d ext = force.external;
// if (mCurrentSimulationStep <= 100) {
// ext[3] = 0;
// ext[4] = 0;
// }
force.residual = force.external;
force.internal = 0;
}
for (size_t ei = 0; ei < pMesh->EN(); ei++) { for (size_t ei = 0; ei < pMesh->EN(); ei++) {
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
std::pair<int, Vector6d> internalForcePair std::pair<int, Vector6d> internalForcePair
@ -1054,10 +1067,23 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (std::isnan(internalForcePair.second.norm())) { if (std::isnan(internalForcePair.second.norm())) {
std::cerr << "nan on edge" << ei << std::endl; std::cerr << "nan on edge" << ei << std::endl;
std::terminate();
} }
#endif #endif
} }
} }
#ifdef POLYSCOPE_DEFINED
std::for_each(pMesh->nodes._handle->data.begin(),
pMesh->nodes._handle->data.end(),
[](Node &node) {
Node::Forces &force = node.force;
const Vector6d debug_residual = force.external - force.internalAxial
- force.internalTorsion
- force.internalFirstBending
- force.internalSecondBending;
assert((force.residual - debug_residual).norm() < 1e-5);
});
#endif
double totalResidualForcesNorm = 0; double totalResidualForcesNorm = 0;
Vector6d sumOfResidualForces(0); Vector6d sumOfResidualForces(0);
for (size_t vi = 0; vi < pMesh->VN(); vi++) { for (size_t vi = 0; vi < pMesh->VN(); vi++) {
@ -1070,7 +1096,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly(
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (std::isnan(force.residual.norm())) { if (std::isnan(force.residual.norm())) {
std::cout << "residual " << vi << ":" << force.residual.toString() << std::endl; std::cout << "residual " << vi << ":" << force.residual.toString() << std::endl;
std::terminate();
} }
// std::cout << "residual " << vi << ":" << force.residual.norm() << std::endl;
#endif #endif
} }
pMesh->totalResidualForcesNorm = totalResidualForcesNorm; pMesh->totalResidualForcesNorm = totalResidualForcesNorm;
@ -1448,6 +1476,7 @@ void DRMSimulationModel::updateT2Derivatives()
{ {
for (const EdgeType &e : pMesh->edge) { for (const EdgeType &e : pMesh->edge) {
for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
const auto ei = pMesh->getIndex(e);
DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi};
pMesh->elements[e].derivativeT2_j[dofi] = computeDerivativeT2(e, dui_v0); pMesh->elements[e].derivativeT2_j[dofi] = computeDerivativeT2(e, dui_v0);
DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi};
@ -1830,8 +1859,11 @@ void DRMSimulationModel::applyForcedNormals(
void DRMSimulationModel::updateKineticEnergy() void DRMSimulationModel::updateKineticEnergy()
{ {
pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy; pMesh->previousTotalKineticEnergy = pMesh->currentTotalKineticEnergy;
pMesh->previousTranslationalKineticEnergy = pMesh->currentTotalTranslationalKineticEnergy;
pMesh->previousTotalRotationalKineticEnergy = pMesh->currentTotalRotationalKineticEnergy;
pMesh->currentTotalKineticEnergy = 0; pMesh->currentTotalKineticEnergy = 0;
pMesh->currentTotalTranslationalKineticEnergy = 0; pMesh->currentTotalTranslationalKineticEnergy = 0;
pMesh->currentTotalRotationalKineticEnergy = 0;
for (const VertexType &v : pMesh->vert) { for (const VertexType &v : pMesh->vert) {
Node &node = pMesh->nodes[v]; Node &node = pMesh->nodes[v];
node.kineticEnergy = 0; node.kineticEnergy = 0;
@ -1853,6 +1885,7 @@ void DRMSimulationModel::updateKineticEnergy()
pMesh->currentTotalKineticEnergy += node.kineticEnergy; pMesh->currentTotalKineticEnergy += node.kineticEnergy;
pMesh->currentTotalTranslationalKineticEnergy += nodeTranslationalKineticEnergy; pMesh->currentTotalTranslationalKineticEnergy += nodeTranslationalKineticEnergy;
pMesh->currentTotalRotationalKineticEnergy += nodeRotationalKineticEnergy;
} }
// assert(mesh->currentTotalKineticEnergy < 100000000000000); // assert(mesh->currentTotalKineticEnergy < 100000000000000);
} }
@ -2009,7 +2042,7 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
// mesh.printVertexCoordinates(mesh.VN() / 2); // mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) { if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) {
draw(); draw(pJob);
} }
#endif #endif
if (!std::isnan(pMesh->currentTotalKineticEnergy)) { if (!std::isnan(pMesh->currentTotalKineticEnergy)) {
@ -2018,46 +2051,46 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptr<Simul
#ifdef COMPUTE_INTERNAL_FORCES #ifdef COMPUTE_INTERNAL_FORCES
results.perVertexInternalForces = computeInternalForces(pJob->constrainedVertices); results.perVertexInternalForces = computeInternalForces(pJob->constrainedVertices);
std::vector<double> perVertexInternalForces_axial = [=]() { // std::vector<double> perVertexInternalForces_axial = [=]() {
std::vector<double> axialInternalForces; // std::vector<double> axialInternalForces;
axialInternalForces.reserve(results.perVertexInternalForces.size()); // axialInternalForces.reserve(results.perVertexInternalForces.size());
for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) { // for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) {
axialInternalForces.push_back(internalForces[0].norm()); // axialInternalForces.push_back(internalForces[0].norm());
} // }
return axialInternalForces; // return axialInternalForces;
}(); // }();
const auto pCurvNet = pMesh->registerForDrawing(); // const auto pCurvNet = pMesh->registerForDrawing();
pCurvNet->addNodeScalarQuantity("axial forces", perVertexInternalForces_axial); // pCurvNet->addNodeScalarQuantity("axial forces", perVertexInternalForces_axial);
std::vector<double> perVertexInternalForces_torsion = [=]() { // std::vector<double> perVertexInternalForces_torsion = [=]() {
std::vector<double> axialInternalForces; // std::vector<double> axialInternalForces;
axialInternalForces.reserve(results.perVertexInternalForces.size()); // axialInternalForces.reserve(results.perVertexInternalForces.size());
for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) { // for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) {
axialInternalForces.push_back(internalForces[1].norm()); // axialInternalForces.push_back(internalForces[1].norm());
} // }
return axialInternalForces; // return axialInternalForces;
}(); // }();
pCurvNet->addNodeScalarQuantity("torsional forces", perVertexInternalForces_torsion); // pCurvNet->addNodeScalarQuantity("torsional forces", perVertexInternalForces_torsion);
std::vector<double> perVertexInternalForces_firstBending = [=]() { // std::vector<double> perVertexInternalForces_firstBending = [=]() {
std::vector<double> axialInternalForces; // std::vector<double> axialInternalForces;
axialInternalForces.reserve(results.perVertexInternalForces.size()); // axialInternalForces.reserve(results.perVertexInternalForces.size());
for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) { // for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) {
axialInternalForces.push_back(internalForces[2].norm()); // axialInternalForces.push_back(internalForces[2].norm());
} // }
return axialInternalForces; // return axialInternalForces;
}(); // }();
pCurvNet->addNodeScalarQuantity("first bending forces", // pCurvNet->addNodeScalarQuantity("first bending forces",
perVertexInternalForces_firstBending); // perVertexInternalForces_firstBending);
std::vector<double> perVertexInternalForces_secondBending = [=]() { // std::vector<double> perVertexInternalForces_secondBending = [=]() {
std::vector<double> axialInternalForces; // std::vector<double> axialInternalForces;
axialInternalForces.reserve(results.perVertexInternalForces.size()); // axialInternalForces.reserve(results.perVertexInternalForces.size());
for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) { // for (const std::array<Vector6d, 4> &internalForces : results.perVertexInternalForces) {
axialInternalForces.push_back(internalForces[3].norm()); // axialInternalForces.push_back(internalForces[3].norm());
} // }
return axialInternalForces; // return axialInternalForces;
}(); // }();
pCurvNet->addNodeScalarQuantity("second bending forces", // pCurvNet->addNodeScalarQuantity("second bending forces",
perVertexInternalForces_secondBending); // perVertexInternalForces_secondBending);
polyscope::show(); // polyscope::show();
#endif #endif
results.rotationalDisplacementQuaternion.resize(pMesh->VN()); results.rotationalDisplacementQuaternion.resize(pMesh->VN());
results.debug_q_f1.resize(pMesh->VN()); results.debug_q_f1.resize(pMesh->VN());
@ -2191,11 +2224,14 @@ void DRMSimulationModel::printDebugInfo() const
} }
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
void DRMSimulationModel::draw(const std::string &screenshotsFolder) void DRMSimulationModel::draw(const std::shared_ptr<SimulationJob> &pJob,
const std::string &screenshotsFolder)
{ {
// update positions // update positions
// polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false); // polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false);
polyscope::CurveNetwork *meshPolyscopeHandle = polyscope::getCurveNetwork(meshPolyscopeLabel); auto deformedMeshPolyscopeLabel = meshPolyscopeLabel + "_" + pMesh->getLabel();
polyscope::CurveNetwork *meshPolyscopeHandle = polyscope::getCurveNetwork(
deformedMeshPolyscopeLabel);
meshPolyscopeHandle->updateNodePositions(pMesh->getEigenVertices()); meshPolyscopeHandle->updateNodePositions(pMesh->getEigenVertices());
// Vertex quantities // Vertex quantities
@ -2259,11 +2295,12 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false); meshPolyscopeHandle->addNodeScalarQuantity("Kinetic Energy", kineticEnergies)->setEnabled(false);
meshPolyscopeHandle->addNodeVectorQuantity("Node normals", nodeNormals)->setEnabled(true); 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 Moments", externalMoments)->setEnabled(true);
meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm) meshPolyscopeHandle->addNodeScalarQuantity("Internal force norm", internalForcesNorm)
->setEnabled(true); ->setEnabled(true);
meshPolyscopeHandle->setRadius(0.002); meshPolyscopeHandle->setRadius(0.002);
// meshPolyscopeHandle->addNodeVectorQuantity("External force", externalForces)->setEnabled(true);
// meshPolyscopeHandle->addNodeVectorQuantity("External Moments", externalMoments)->setEnabled(true);
pJob->registerForDrawing(deformedMeshPolyscopeLabel, true);
// polyscope::getCurveNetwork(meshPolyscopeLabel) // polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) // ->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
// ->setEnabled(false); // ->setEnabled(false);
@ -2284,7 +2321,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
// ->addNodeVectorQuantity("Second bending force-Rotation", // ->addNodeVectorQuantity("Second bending force-Rotation",
// internalSecondBendingRotationForces) // internalSecondBendingRotationForces)
// ->setEnabled(false); // ->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel) polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)
->addNodeScalarQuantity("nR", nRs) ->addNodeScalarQuantity("nR", nRs)
->setEnabled(false); ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel) // polyscope::getCurveNetwork(meshPolyscopeLabel)
@ -2293,10 +2330,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
// polyscope::getCurveNetwork(meshPolyscopeLabel) // polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeScalarQuantity("theta2", theta2) // ->addNodeScalarQuantity("theta2", theta2)
// ->setEnabled(false); // ->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel) polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)
->addNodeVectorQuantity("Residual force", residualForces) ->addNodeVectorQuantity("Residual force", residualForces)
->setEnabled(false); ->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel) polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)
->addNodeScalarQuantity("Residual force norm", residualForcesNorm) ->addNodeScalarQuantity("Residual force norm", residualForcesNorm)
->setEnabled(false); ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel) // polyscope::getCurveNetwork(meshPolyscopeLabel)
@ -2315,11 +2352,10 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
I3[ei] = pMesh->elements[e].dimensions.inertia.I3; I3[ei] = pMesh->elements[e].dimensions.inertia.I3;
} }
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("A", A);
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J); polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("J", J);
polyscope::getCurveNetwork(meshPolyscopeLabel) polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("I2", I2);
->addEdgeScalarQuantity("I2", I2); polyscope::getCurveNetwork(deformedMeshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3);
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("I3", I3);
// Specify the callback // Specify the callback
static bool calledOnce = false; static bool calledOnce = false;
@ -2574,7 +2610,19 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr<SimulationJob
SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob) SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<SimulationJob> &pJob)
{ {
auto beginTime = std::chrono::high_resolution_clock::now(); auto beginTime = std::chrono::high_resolution_clock::now();
updateNodalMasses(); updateNodalMasses();
// constexpr bool useDRM = true;
//#ifdef USE_ENSMALLEN
// if (!useDRM) {
// setJob(pJob);
// // ens::L_BFGS optimizer(20);
// ens::SA optimizer;
// arma::mat x(pJob->pMesh->VN() * NumDoF, 1);
// optimizer.Optimize(*this, x);
// // getD
// } else {
//#endif
// std::unordered_map<VertexIndex, Vector6d> nodalExternalForces = pJob->nodalExternalForces; // std::unordered_map<VertexIndex, Vector6d> nodalExternalForces = pJob->nodalExternalForces;
// double totalExternalForcesNorm = 0; // double totalExternalForcesNorm = 0;
// Vector6d sumOfExternalForces(0); // Vector6d sumOfExternalForces(0);
@ -2594,197 +2642,224 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
mSettings.totalResidualForcesNormThreshold = 1e-3; mSettings.totalResidualForcesNormThreshold = 1e-3;
std::cout << "No forces setted default residual forces norm threshold" << std::endl; std::cout << "No forces setted default residual forces norm threshold" << std::endl;
} }
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" std::cout << "totalResidualForcesNormThreshold:"
<< mSettings.totalResidualForcesNormThreshold << std::endl; << mSettings.totalResidualForcesNormThreshold << std::endl;
if (mSettings.averageResidualForcesCriterionThreshold.has_value()) { if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
std::cout << "average/extNorm threshold:" std::cout << "average/extNorm threshold:"
<< *mSettings.averageResidualForcesCriterionThreshold << std::endl; << *mSettings.averageResidualForcesCriterionThreshold << std::endl;
}
} }
}
if (mSettings.beVerbose) { if (mSettings.beVerbose) {
std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and "
<< pMesh->EN() << " elements." << std::endl; << pMesh->EN() << " elements." << std::endl;
} }
const size_t movingAverageSampleSize = 200; const size_t movingAverageSampleSize = 200;
std::queue<double> residualForcesMovingAverageHistorySample; std::queue<double> residualForcesMovingAverageHistorySample;
std::vector<double> percentageOfConvergence; std::vector<double> percentageOfConvergence;
// double residualForcesMovingAverageDerivativeMax = 0; // double residualForcesMovingAverageDerivativeMax = 0;
while (!mSettings.maxDRMIterations.has_value() while (!mSettings.maxDRMIterations.has_value()
|| mCurrentSimulationStep < mSettings.maxDRMIterations.value()) { || mCurrentSimulationStep < mSettings.maxDRMIterations.value()) {
if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { // if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) {
// std::filesystem::create_directory("./PatternOptimizationNonConv"); // std::filesystem::create_directory("./PatternOptimizationNonConv");
// pJob->save("./PatternOptimizationNonConv"); // pJob->save("./PatternOptimizationNonConv");
// Dt = mSettings.Dtini; // Dt = mSettings.Dtini;
} // }
// if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) { // if (mCurrentSimulationStep == 500 && shouldTemporarilyDampForces) {
// Dt = mSettings.Dtini; // Dt = mSettings.Dtini;
// } // }
// while (true) { // while (true) {
updateNormalDerivatives(); updateNormalDerivatives();
updateT1Derivatives(); updateT1Derivatives();
updateRDerivatives(); updateRDerivatives();
updateT2Derivatives(); updateT2Derivatives();
updateT3Derivatives(); updateT3Derivatives();
const bool shouldBreak = mCurrentSimulationStep == 12970; const bool shouldBreak = mCurrentSimulationStep == 3935;
updateResidualForcesOnTheFly(constrainedVertices); if (shouldBreak) {
int i = 0;
// TODO: write parallel function for updating positions i++;
// TODO: Make a single function out of updateResidualForcesOnTheFly
// updatePositionsOnTheFly
// updatePositionsOnTheFly(constrainedVertices);
updateNodalMasses();
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements);
}
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
mCurrentSimulationStep++;
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
std::cout << pMesh->currentTotalKineticEnergy << std::endl;
if (mSettings.beVerbose) {
std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep
<< ". Exiting.." << std::endl;
} }
break; updateResidualForcesOnTheFly(constrainedVertices);
}
if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) { // TODO: write parallel function for updating positions
minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm; // TODO: Make a single function out of updateResidualForcesOnTheFly
static int lastSavedStep = mCurrentSimulationStep; // updatePositionsOnTheFly
if (mSettings.saveIntermediateBestStates.has_value() // updatePositionsOnTheFly(constrainedVertices);
&& mSettings.saveIntermediateBestStates.value() updateNodalMasses();
&& mCurrentSimulationStep - lastSavedStep > 20000) { updateNodalAccelerations();
lastSavedStep = mCurrentSimulationStep; updateNodalVelocities();
computeResults(pJob).save("./IntermediateResults_" + pJob->getLabel()); updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
pJob->nodalForcedDisplacements);
} }
} // if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
//normalized sum of displacements // }
// double sumOfDisplacementsNorm = 0; updateElementalLengths();
// for (size_t vi = 0; vi < pMesh->VN(); vi++) { mCurrentSimulationStep++;
// sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); if (std::isnan(pMesh->currentTotalKineticEnergy)) {
// } std::cout << pMesh->currentTotalKineticEnergy << std::endl;
// sumOfDisplacementsNorm /= pMesh->bbox.Diag(); if (mSettings.beVerbose) {
// pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; std::cerr << "Infinite kinetic energy at step " << mCurrentSimulationStep
<< ". Exiting.." << std::endl;
//moving average }
// // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm break;
// // + 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->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;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
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<std::chrono::minutes>(t2 - beginTime).count();
// std::cout << "Execution time(min):" << executionTime_min << std::endl;
if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
std::cout << "Best percentage of target (average):"
<< 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;
} }
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// history.residualForcesMovingAverage);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
// draw();
// SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps",
// "Time/(#nodes*#iterations)",
// timePerNodePerIterationHistory);
}
if ((mSettings.shouldCreatePlots || mSettings.debugModeStep.has_value()) if (minTotalResidualForcesNorm > pMesh->totalResidualForcesNorm) {
&& mCurrentSimulationStep != 0) { minTotalResidualForcesNorm = pMesh->totalResidualForcesNorm;
history.stepPulse(*pMesh); static int lastSavedStep = mCurrentSimulationStep;
percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold if (mSettings.saveIntermediateBestStates.has_value()
/ pMesh->totalResidualForcesNorm); && mSettings.saveIntermediateBestStates.value()
} && mCurrentSimulationStep - lastSavedStep > 20000) {
lastSavedStep = mCurrentSimulationStep;
computeResults(pJob).save("./IntermediateResults_" + pJob->getLabel());
}
}
if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value() // pMesh->previousTotalPotentialEnergykN =
&& mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) { // pMesh->currentTotalPotentialEnergykN;
// SimulationResultsReporter::createPlot("Number of Steps", // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// "Residual Forces mov aver deriv", // timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
// movingAveragesDerivatives_norm); if (mSettings.beVerbose && mSettings.debugModeStep.has_value()
// SimulationResultsReporter::createPlot("Number of Steps", && mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) {
// "Residual Forces mov aver", printCurrentState();
// history.residualForcesMovingAverage,
// {}, // auto t2 = std::chrono::high_resolution_clock::now();
// history.redMarks); // const double executionTime_min
SimulationResultsReporter::createPlot("Number of Steps", // = std::chrono::duration_cast<std::chrono::minutes>(t2 - beginTime).count();
"Log of Residual Forces", // std::cout << "Execution time(min):" << executionTime_min << std::endl;
history.logResidualForces, if (mSettings.averageResidualForcesCriterionThreshold.has_value()) {
{}, std::cout << "Best percentage of target (average):"
history.redMarks); << 100 * mSettings.averageResidualForcesCriterionThreshold.value()
// SimulationResultsReporter::createPlot("Number of Steps", * pMesh->totalExternalForcesNorm * pMesh->VN()
// "Log of Kinetic energy", / minTotalResidualForcesNorm
// history.kineticEnergy, << "%" << std::endl;
// {}, } else {
// history.redMarks); std::cout << "Best percentage of target:"
// SimulationResultsReporter reporter; << 100 * mSettings.totalExternalForcesNormPercentageTermination
// reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel()); * pMesh->totalExternalForcesNorm / minTotalResidualForcesNorm
// SimulationResultsReporter::createPlot("Number of Iterations", << "%" << std::endl;
// "Sum of normalized displacement norms", }
// history.sumOfNormalizedDisplacementNorms /*, // SimulationResultsReporter::createPlot("Number of Steps",
// std::filesystem::path("./") // "Residual Forces mov aver",
// .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") // history.residualForcesMovingAverage);
// .string()*/ // SimulationResultsReporter::createPlot("Number of Steps",
// , // "Residual Forces mov aver deriv",
// {}, // movingAveragesDerivatives);
// history.redMarks); // draw();
// SimulationResultsReporter::createPlot("Number of Steps", // SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps",
// "Percentage of convergence", // "Time/(#nodes*#iterations)",
// percentageOfConvergence, // timePerNodePerIterationHistory);
// {}, }
// history.redMarks);
} if ((mSettings.shouldCreatePlots || mSettings.debugModeStep.has_value())
&& mCurrentSimulationStep != 0) {
//normalized sum of displacements
double sumOfDisplacementsNorm = 0;
for (size_t vi = 0; vi < pMesh->VN(); vi++) {
sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm();
}
pMesh->perVertexAverageNormalizedDisplacementNorm = sumOfDisplacementsNorm
/ (pMesh->bbox.Diag()
* pMesh->VN());
//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->residualForcesMovingAverageDerivativeNorm
= std::abs((residualForcesMovingAverageHistorySample.back()
- residualForcesMovingAverageHistorySample.front()))
/ (movingAverageSampleSize);
// residualForcesMovingAverageDerivativeMax
// = std::max(pMesh->residualForcesMovingAverageDerivativeNorm,
// residualForcesMovingAverageDerivativeMax);
// pMesh->residualForcesMovingAverageDerivativeNorm
// /= residualForcesMovingAverageDerivativeMax;
// std::cout << "Normalized derivative:"
// << residualForcesMovingAverageDerivativeNorm
// << std::endl;
}
history.stepPulse(*pMesh);
percentageOfConvergence.push_back(100 * mSettings.totalResidualForcesNormThreshold
/ pMesh->totalResidualForcesNorm);
}
if (mSettings.shouldCreatePlots && mSettings.debugModeStep.has_value()
&& mCurrentSimulationStep % mSettings.debugModeStep.value() == 0) {
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// history
// .residualForcesMovingAverageDerivativesLog);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// history.residualForcesMovingAverage,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Sum of axial forces norms",
// history.logOfSumOfAxialForcesNorm,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Sum of torsion forces norms",
// history.logOfSumOfTorsionForcesNorm,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Sum of 1stBending forces norms",
// history.logOfSumOfFirstBendingForcesNorm,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Sum of 2ndBending forces norms",
// history.logOfSumOfSecondBendingForcesNorm,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Residual Forces",
// history.logResidualForces,
// {},
// history.redMarks);
SimulationResultsReporter::createPlot("Number of Steps",
"Log of Kinetic energy",
history.kineticEnergy,
{},
history.redMarks);
// SimulationResultsReporter reporter;
// reporter.reportHistory(history, "IntermediateResults", pJob->pMesh->getLabel());
// SimulationResultsReporter::createPlot(
// "Number of Iterations",
// "Per vertex average normalized displacement norms",
// history.perVertexAverageNormalizedDisplacementNorm /*,
// std::filesystem::path("./")
// .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
// .string()*/
// ,
// {},
// history.redMarks);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Percentage of convergence",
// percentageOfConvergence,
// {},
// history.redMarks);
}
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && mSettings.debugModeStep.has_value() if (mSettings.shouldDraw && mSettings.debugModeStep.has_value()
@ -2796,7 +2871,7 @@ currentSimulationStep > maxDRMIterations*/
// .append("Debugging_files") // .append("Debugging_files")
// .append("Screenshots") // .append("Screenshots")
// .string(); // .string();
draw(); draw(pJob);
// yValues = history.kineticEnergy; // yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues); // plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy"; // label = "Log of Kinetic energy";
@ -2837,17 +2912,14 @@ currentSimulationStep > maxDRMIterations*/
} }
break; break;
} }
// Residual forces norm convergence
if ((pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy const bool isKineticEnergyPeak = (mSettings.useTotalRotationalKineticEnergyForKineticDamping
// && mCurrentSimulationStep > movingAverageSampleSize && pMesh->previousTotalRotationalKineticEnergy
) > pMesh->currentTotalRotationalKineticEnergy)
/* || (pMesh->totalResidualForcesNorm / mSettings.totalResidualForcesNormThreshold <= 1 || (mSettings.useTranslationalKineticEnergyForKineticDamping
&& mCurrentSimulationStep > 1)*/ && pMesh->previousTranslationalKineticEnergy
/*|| > pMesh->currentTotalTranslationalKineticEnergy);
mesh->previousTotalPotentialEnergykN > if (isKineticEnergyPeak) {
mesh->currentTotalPotentialEnergykN*/
/*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) {
// if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) {
const bool fullfillsKineticEnergyTerminationCriterion const bool fullfillsKineticEnergyTerminationCriterion
= mSettings.totalTranslationalKineticEnergyThreshold.has_value() = mSettings.totalTranslationalKineticEnergyThreshold.has_value()
&& pMesh->currentTotalTranslationalKineticEnergy && pMesh->currentTotalTranslationalKineticEnergy
@ -2899,7 +2971,7 @@ mesh->currentTotalPotentialEnergykN*/
// } // }
} }
if (mSettings.useKineticDamping) { if (mSettings.useTranslationalKineticEnergyForKineticDamping) {
applyKineticDamping(pJob); applyKineticDamping(pJob);
Dt *= mSettings.xi; Dt *= mSettings.xi;
if (mSettings.shouldCreatePlots) { if (mSettings.shouldCreatePlots) {
@ -2912,22 +2984,24 @@ mesh->currentTotalPotentialEnergykN*/
// std::cout << "Number of dampings:" << numOfDampings << endl; // std::cout << "Number of dampings:" << numOfDampings << endl;
// draw(); // draw();
} }
} }
auto endTime = std::chrono::high_resolution_clock::now(); //#ifdef USE_ENSMALLEN
// }
//#endif
auto endTime = std::chrono::high_resolution_clock::now();
SimulationResults results = computeResults(pJob);
results.executionTime
= std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime).count();
SimulationResults results = computeResults(pJob); if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) {
results.executionTime = std::chrono::duration_cast<std::chrono::seconds>(endTime - beginTime) SimulationResultsReporter reporter;
.count(); reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}
if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) {
SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) {
polyscope::removeCurveNetwork(meshPolyscopeLabel); polyscope::removeCurveNetwork(meshPolyscopeLabel + "_" + pMesh->getLabel());
polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel + "_" + pMesh->getLabel());
} }
#endif #endif
return results; return results;
@ -2936,6 +3010,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
const Settings &settings, const Settings &settings,
const SimulationResults &solutionGuess) const SimulationResults &solutionGuess)
{ {
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
reset(pJob, settings); reset(pJob, settings);
assert(pJob->pMesh != nullptr); assert(pJob->pMesh != nullptr);
@ -2953,16 +3028,26 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
LinearSimulationModel linearSimulationModel; LinearSimulationModel linearSimulationModel;
SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel
.executeSimulation(pJob); .executeSimulation(pJob);
// simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); if (simulationResults_fullPatternLinearModel.converged) {
shouldTemporarilyDampForces = true; // simulationResults_fullPatternLinearModel.save(fullResultsFolderPath);
for (auto &externalForce : pJob->nodalExternalForces) {
externalForce.second = externalForce.second / forceScaleFactor;
}
applySolutionGuess(simulationResults_fullPatternLinearModel, pJob);
shouldTemporarilyDampForces = true;
}
} }
SimulationResults results = executeSimulation(pJob); //calling the virtual function SimulationResults results = executeSimulation(pJob); //calling the virtual function
std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
if (settings.beVerbose) {
std::cout << "Simulation converged in "
<< std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count()
/ (60 * 1000.0)
<< " min" << std::endl;
}
return results; return results;
} }
@ -2983,7 +3068,8 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// // } // // }
//} //}
//SimulationMesh *DRMSimulationModel::getDeformedMesh(const arma::mat &x) //SimulationMesh *DRMSimulationModel::getDeformedMesh(const arma::mat &x,
// const std::shared_ptr<SimulationJob> &pJob)
//{ //{
// for (int vi = 0; vi < pMesh->VN(); vi++) { // for (int vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi]; // Node &node = pMesh->nodes[vi];
@ -3000,6 +3086,32 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// return pMesh.get(); // return pMesh.get();
//} //}
//double DRMSimulationModel::Evaluate(const arma::mat &x)
//{
// for (int vi = 0; vi < pMesh->VN(); vi++) {
// Node &node = pMesh->nodes[vi];
// for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) {
// node.displacements[dofi] = x(vi * DoF::NumDoF + dofi);
// }
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(pJob->nodalForcedDisplacements);
// }
// updateElementalLengths();
// updateNormalDerivatives();
// updateT1Derivatives();
// updateRDerivatives();
// updateT2Derivatives();
// updateT3Derivatives();
// updateResidualForcesOnTheFly(constrainedVertices);
// const double PE = computeTotalPotentialEnergy();
// // const double PE = computeTotalInternalPotentialEnergy();
// // std::cout << "PE:" << PE << std::endl;
// return pMesh->totalResidualForcesNorm;
//}
//double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g) //double DRMSimulationModel::EvaluateWithGradient(const arma::mat &x, arma::mat &g)
//{ //{
// for (int vi = 0; vi < pMesh->VN(); vi++) { // for (int vi = 0; vi < pMesh->VN(); vi++) {
@ -3028,7 +3140,13 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
// } // }
// mCurrentSimulationStep++; // mCurrentSimulationStep++;
// const double PE = computeTotalPotentialEnergy(); // const double PE = computeTotalPotentialEnergy();
// // const double PE = computeTotalInternalPotentialEnergy();
// std::cout << "PE:" << PE << std::endl; // std::cout << "PE:" << PE << std::endl;
// double sumOfDisplacements = 0;
// for (int xi = 0; xi < pJob->pMesh->VN() * NumDoF; xi++) {
// sumOfDisplacements += x(xi) * x(xi);
// }
// std::cout << "Sum of disp:" << std::sqrt(sumOfDisplacements) << std::endl;
// return PE; // return PE;
//} //}
@ -3043,7 +3161,8 @@ void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePat
json[GET_VARIABLE_NAME(Dtini)] = Dtini; json[GET_VARIABLE_NAME(Dtini)] = Dtini;
json[GET_VARIABLE_NAME(xi)] = xi; json[GET_VARIABLE_NAME(xi)] = xi;
json[GET_VARIABLE_NAME(gamma)] = gamma; json[GET_VARIABLE_NAME(gamma)] = gamma;
json[GET_VARIABLE_NAME(useKineticDamping)] = totalResidualForcesNormThreshold; json[GET_VARIABLE_NAME(useTranslationalKineticEnergyForKineticDamping)]
= useTranslationalKineticEnergyForKineticDamping;
if (maxDRMIterations.has_value()) { if (maxDRMIterations.has_value()) {
json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value(); json[GET_VARIABLE_NAME(jsonLabels.maxDRMIterations)] = maxDRMIterations.value();
} }
@ -3064,7 +3183,8 @@ void DRMSimulationModel::Settings::save(const std::filesystem::path &jsonFilePat
if (viscousDampingFactor.has_value()) { if (viscousDampingFactor.has_value()) {
json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value(); json[GET_VARIABLE_NAME(viscousDampingFactor)] = viscousDampingFactor.value();
} }
json[GET_VARIABLE_NAME(useKineticDamping)] = useKineticDamping; json[GET_VARIABLE_NAME(useTranslationalKineticEnergyForKineticDamping)]
= useTranslationalKineticEnergyForKineticDamping;
// if (intermediateResultsSaveStep.has_value()) { // if (intermediateResultsSaveStep.has_value()) {
// json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value(); // json[GET_VARIABLE_NAME(intermediateResultsSaveStep)] = intermediateResultsSaveStep.value();
// } // }

22
drmsimulationmodel.hpp
View File

@ -16,7 +16,7 @@
#include <ensmallen.hpp> #include <ensmallen.hpp>
#endif #endif
struct SimulationJob; class SimulationJob;
enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF }; enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF };
using DoFType = int; using DoFType = int;
@ -33,6 +33,8 @@ class DRMSimulationModel : public SimulationModel
public: public:
struct Settings struct Settings
{ {
bool useTranslationalKineticEnergyForKineticDamping{true};
bool useTotalRotationalKineticEnergyForKineticDamping{false};
bool shouldDraw{false}; bool shouldDraw{false};
bool beVerbose{false}; bool beVerbose{false};
bool shouldCreatePlots{false}; bool shouldCreatePlots{false};
@ -45,7 +47,7 @@ public:
int gradualForcedDisplacementSteps{50}; int gradualForcedDisplacementSteps{50};
// int desiredGradualExternalLoadsSteps{1}; // int desiredGradualExternalLoadsSteps{1};
double gamma{0.8}; double gamma{0.8};
double totalResidualForcesNormThreshold{1e-8}; double totalResidualForcesNormThreshold{1e-20};
double totalExternalForcesNormPercentageTermination{1e-5}; double totalExternalForcesNormPercentageTermination{1e-5};
std::optional<int> maxDRMIterations; std::optional<int> maxDRMIterations;
std::optional<int> debugModeStep; std::optional<int> debugModeStep;
@ -55,7 +57,6 @@ public:
// std::optional<int> intermediateResultsSaveStep; // std::optional<int> intermediateResultsSaveStep;
std::optional<bool> saveIntermediateBestStates; std::optional<bool> saveIntermediateBestStates;
std::optional<double> viscousDampingFactor; std::optional<double> viscousDampingFactor;
bool useKineticDamping{true};
Settings() {} Settings() {}
void save(const std::filesystem::path &jsonFilePath) const; void save(const std::filesystem::path &jsonFilePath) const;
bool load(const std::filesystem::path &filePath); bool load(const std::filesystem::path &filePath);
@ -153,7 +154,7 @@ private:
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices); const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> &fixedVertices);
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
void draw(const std::string &screenshotsFolder = {}); void draw(const std::shared_ptr<SimulationJob> &pJob, const std::string &screenshotsFolder = {});
#endif #endif
void void
updateNodalInternalForce(Vector6d &nodalInternalForce, updateNodalInternalForce(Vector6d &nodalInternalForce,
@ -262,12 +263,13 @@ private:
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob, SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings, const Settings &settings,
const SimulationResults &solutionGuess = SimulationResults()); const SimulationResults &solutionGuess = SimulationResults());
#ifdef USE_ENSMALLEN //#ifdef USE_ENSMALLEN
double EvaluateWithGradient(const arma::mat &x, arma::mat &g); // std::shared_ptr<SimulationJob> pJob;
void setJob(const std::shared_ptr<SimulationJob> &pJob); // double EvaluateWithGradient(const arma::mat &x, arma::mat &g);
SimulationMesh *getDeformedMesh(const arma::mat &x); // void setJob(const std::shared_ptr<SimulationJob> &pJob);
// SimulationMesh *getDeformedMesh(const arma::mat &x, const std::shared_ptr<SimulationJob> &pJob);
#endif // double Evaluate(const arma::mat &x);
//#endif
static void runUnitTests(); static void runUnitTests();
}; };