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

Refactoring

This commit is contained in:
Iason 2021-01-12 14:40:25 +02:00
parent 2e787689fe
commit d933f92eeb
9 changed files with 273 additions and 212 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
beam.hpp
View File

@ -12,7 +12,7 @@ struct RectangularBeamDimensions {
: b(width), h(height) {
assert(width > 0 && height > 0);
}
RectangularBeamDimensions() : b(1), h(1) {}
RectangularBeamDimensions() : b(0.002), h(0.002) {}
};
struct CylindricalBeamDimensions {

377
beamformfinder.cpp
View File

@ -8,18 +8,19 @@
#include <execution>
#include <omp.h>
const bool debug = true;
void FormFinder::runUnitTests() {
const std::filesystem::path groundOfTruthFolder{
"/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"};
FormFinder formFinder;
formFinder.setTotalResidualForcesNormThreshold(1);
// First example of the paper
VCGEdgeMesh beam;
// const size_t spanGridSize = 11;
// mesh.createSpanGrid(spanGridSize);
beam.loadFromPly(
beam.loadPly(
std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply"));
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3};
@ -39,8 +40,13 @@ void FormFinder::runUnitTests() {
assert(CrossSectionType::name == CylindricalBeamDimensions::name);
beamSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
Settings settings;
settings.Dtini = 0.1;
settings.xi = 0.9969;
settings.maxDRMIterations = 0;
settings.totalResidualForcesNormThreshold = 1;
SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(
std::make_shared<SimulationJob>(beamSimulationJob), false, true);
std::make_shared<SimulationJob>(beamSimulationJob), settings);
simpleBeam_simulationResults.save();
const std::string simpleBeamGroundOfTruthBinaryFilename =
std::filesystem::path(groundOfTruthFolder)
@ -60,9 +66,9 @@ void FormFinder::runUnitTests() {
VCGEdgeMesh shortSpanGrid;
// const size_t spanGridSize = 11;
// mesh.createSpanGrid(spanGridSize);
shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
.append("shortSpanGridshell.ply")
.string());
shortSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder)
.append("shortSpanGridshell.ply")
.string());
fixedVertices.clear();
//// Corner nodes
@ -92,7 +98,7 @@ void FormFinder::runUnitTests() {
SimulationResults shortSpanGridshellSimulationResults =
formFinder.executeSimulation(
std::make_shared<SimulationJob>(shortSpanGridshellSimulationJob),
false, false);
settings);
shortSpanGridshellSimulationResults.save();
const std::string groundOfTruthBinaryFilename =
@ -118,9 +124,9 @@ void FormFinder::runUnitTests() {
// Third example of the paper
VCGEdgeMesh longSpanGrid;
longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder)
.append("longSpanGridshell.ply")
.string());
longSpanGrid.loadPly(std::filesystem::path(groundOfTruthFolder)
.append("longSpanGridshell.ply")
.string());
const size_t spanGridSize = 11;
fixedVertices.clear();
@ -174,7 +180,7 @@ void FormFinder::runUnitTests() {
SimulationResults longSpanGridshell_simulationResults =
formFinder.executeSimulation(
std::make_shared<SimulationJob>(longSpanGridshell_simulationJob),
false, false);
settings);
longSpanGridshell_simulationResults.save();
const std::string longSpan_groundOfTruthBinaryFilename =
@ -200,12 +206,7 @@ void FormFinder::runUnitTests() {
// polyscope::show();
}
void FormFinder::setTotalResidualForcesNormThreshold(double value) {
settings.totalResidualForcesNormThreshold = value;
}
void FormFinder::reset() {
settings.Dt = settings.Dtini;
mCurrentSimulationStep = 0;
history.clear();
constrainedVertices.clear();
@ -214,9 +215,10 @@ void FormFinder::reset() {
plotYValues.clear();
plotHandle.reset();
checkedForMaximumMoment = false;
shouldUseKineticEnergyThreshold = false;
mSettings.shouldUseTranslationalKineticEnergyThreshold = false;
externalMomentsNorm = 0;
settings.drawingStep = 1;
mSettings.drawingStep = 1;
Dt = mSettings.Dtini;
}
VectorType FormFinder::computeDisplacementDifferenceDerivative(
@ -841,7 +843,7 @@ void FormFinder::updateResidualForcesOnTheFly(
const double e_k = element.length - element.initialLength;
const double e_kDeriv = computeDerivativeElementLength(e, dui);
const double axialForce_dofi =
e_kDeriv * e_k * element.axialConstFactor;
e_kDeriv * e_k * element.rigidity.axial;
// Torsional force computation
const double theta1_j_deriv =
@ -852,7 +854,7 @@ void FormFinder::updateResidualForcesOnTheFly(
const double theta1DiffDerivative =
theta1_jplus1_deriv - theta1_j_deriv;
const double torsionalForce_dofi =
theta1DiffDerivative * theta1Diff * element.torsionConstFactor;
theta1DiffDerivative * theta1Diff * element.rigidity.torsional;
// First bending force computation
////theta2_j derivative
@ -880,8 +882,7 @@ void FormFinder::updateResidualForcesOnTheFly(
firstBendingForce_inBracketsTerm_2 +
firstBendingForce_inBracketsTerm_3;
const double firstBendingForce_dofi =
firstBendingForce_inBracketsTerm *
element.firstBendingConstFactor;
firstBendingForce_inBracketsTerm * element.rigidity.firstBending;
// Second bending force computation
////theta2_j derivative
@ -908,7 +909,7 @@ void FormFinder::updateResidualForcesOnTheFly(
secondBendingForce_inBracketsTerm_2 +
secondBendingForce_inBracketsTerm_3;
double secondBendingForce_dofi = secondBendingForce_inBracketsTerm *
element.secondBendingConstFactor;
element.rigidity.secondBending;
const bool shouldBreak =
mCurrentSimulationStep == 118 && edgeNode.vi == 1 && dofi == 3;
internalForcesContributionFromThisEdge[evi].second[dofi] =
@ -947,7 +948,7 @@ void FormFinder::updateResidualForcesOnTheFly(
secondBendingForce_inBracketsTerm_3;
const double secondBendingForce_dofi =
secondBendingForce_inBracketsTerm *
element.secondBendingConstFactor;
element.rigidity.secondBending;
internalForcesContributionFromThisEdge[evi + 2].second[dofi] =
secondBendingForce_dofi;
}
@ -1162,24 +1163,24 @@ void FormFinder::updateNodalMasses() {
assert(rotationalSumSk_J != 0);
}
mesh->nodes[v].translationalMass =
gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk;
gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk;
mesh->nodes[v].rotationalMass_I2 =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2;
mesh->nodes[v].rotationalMass_I3 =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3;
mesh->nodes[v].rotationalMass_J =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J;
assert(std::pow(settings.Dtini, 2.0) * translationalSumSk /
assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk /
mesh->nodes[v].translationalMass <
2);
assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I2 /
assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I2 /
mesh->nodes[v].rotationalMass_I2 <
2);
assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_I3 /
assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I3 /
mesh->nodes[v].rotationalMass_I3 <
2);
assert(std::pow(settings.Dtini, 2.0) * rotationalSumSk_J /
assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_J /
mesh->nodes[v].rotationalMass_J <
2);
}
@ -1211,7 +1212,7 @@ void FormFinder::updateNodalVelocities() {
for (VertexType &v : mesh->vert) {
const VertexIndex vi = mesh->getIndex(v);
Node &node = mesh->nodes[v];
node.velocity = node.velocity + node.acceleration * settings.Dt;
node.velocity = node.velocity + node.acceleration * Dt;
}
updateKineticEnergy();
}
@ -1219,14 +1220,17 @@ void FormFinder::updateNodalVelocities() {
void FormFinder::updateNodalDisplacements() {
for (VertexType &v : mesh->vert) {
Node &node = mesh->nodes[v];
node.displacements = node.displacements + node.velocity * settings.Dt;
if (settings.beVerbose) {
std::cout << "Node " << node.vi << ":" << endl;
std::cout << node.displacements[0] << " " << node.displacements[0] << " "
<< node.displacements[1] << " " << node.displacements[2] << " "
<< node.displacements[3] << " " << node.displacements[4] << " "
<< node.displacements[5] << std::endl;
}
node.displacements = node.displacements + node.velocity * Dt;
// if (mSettings.beVerbose) {
// std::cout << "Node " << node.vi << ":" << endl;
// std::cout << node.displacements[0] << " " << node.displacements[0]
// << " "
// << node.displacements[1] << " " << node.displacements[2]
// << " "
// << node.displacements[3] << " " << node.displacements[4]
// << " "
// << node.displacements[5] << std::endl;
// }
}
}
@ -1299,7 +1303,7 @@ void FormFinder::updateNodeNormal(
const bool viHasMoments =
node.force.external[3] != 0 || node.force.external[4] != 0;
if (!checkedForMaximumMoment && viHasMoments) {
shouldUseKineticEnergyThreshold = true;
mSettings.shouldUseTranslationalKineticEnergyThreshold = true;
std::cout
<< "Maximum moment reached.The Kinetic energy of the system will "
"be used as a convergence criterion"
@ -1340,7 +1344,7 @@ void FormFinder::applyDisplacements(
updateNodeNormal(v, fixedVertices);
}
updateElementalFrames();
if (settings.shouldDraw || true) {
if (mSettings.shouldDraw || true) {
mesh->updateEigenEdgeAndVertices();
}
}
@ -1392,6 +1396,7 @@ void FormFinder::applyForcedNormals(
void FormFinder::updateKineticEnergy() {
mesh->previousTotalKineticEnergy = mesh->currentTotalKineticEnergy;
mesh->currentTotalKineticEnergy = 0;
mesh->currentTotalTranslationalKineticEnergy = 0;
for (const VertexType &v : mesh->vert) {
Node &node = mesh->nodes[v];
node.kineticEnergy = 0;
@ -1399,19 +1404,21 @@ void FormFinder::updateKineticEnergy() {
const double translationalVelocityNorm = std::sqrt(
std::pow(node.velocity[0], 2) + std::pow(node.velocity[1], 2) +
std::pow(node.velocity[2], 2));
node.kineticEnergy +=
const double nodeTranslationalKineticEnergy =
0.5 * node.translationalMass * pow(translationalVelocityNorm, 2);
const double nodeRotationalKineticEnergy =
0.5 * (node.rotationalMass_J * pow(node.velocity[3], 2) +
+node.rotationalMass_I3 * pow(node.velocity[4], 2) +
+node.rotationalMass_I2 * pow(node.velocity[5], 2));
node.kineticEnergy +=
nodeTranslationalKineticEnergy /*+ nodeRotationalKineticEnergy*/;
assert(node.kineticEnergy < 10000000000000000000);
// const double rotationalVelocityNorm = std::sqrt(
// std::pow(node.velocity[3], 2) + std::pow(node.velocity[4], 2) +
// std::pow(node.velocity[5], 2));
// node.kineticEnergy +=
// 0.5 * node.rotationalMass_J * pow(node.velocity[3], 2) +
// 0.5 * node.rotationalMass_I3 * pow(node.velocity[4], 2) +
// 0.5 * node.rotationalMass_I2 * pow(node.velocity[5], 2);
mesh->currentTotalKineticEnergy += node.kineticEnergy;
mesh->currentTotalTranslationalKineticEnergy +=
nodeTranslationalKineticEnergy;
}
// assert(mesh->currentTotalKineticEnergy < 100000000000000);
}
@ -1460,24 +1467,24 @@ void FormFinder::updatePositionsOnTheFly(
// assert(rotationalSumSk_J != 0);
}
mesh->nodes[v].translationalMass =
gamma * pow(settings.Dtini, 2) * 2 * translationalSumSk;
gamma * pow(mSettings.Dtini, 2) * 2 * translationalSumSk;
mesh->nodes[v].rotationalMass_I2 =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I2;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I2;
mesh->nodes[v].rotationalMass_I3 =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_I3;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_I3;
mesh->nodes[v].rotationalMass_J =
gamma * pow(settings.Dtini, 2) * 8 * rotationalSumSk_J;
gamma * pow(mSettings.Dtini, 2) * 8 * rotationalSumSk_J;
// assert(std::pow(Dtini, 2.0) * translationalSumSk /
// assert(std::pow(mSettings.Dtini, 2.0) * translationalSumSk /
// mesh->nodes[v].translationalMass <
// 2);
// assert(std::pow(Dtini, 2.0) * rotationalSumSk_I2 /
// assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I2 /
// mesh->nodes[v].rotationalMass_I2 <
// 2);
// assert(std::pow(Dtini, 2.0) * rotationalSumSk_I3 /
// assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_I3 /
// mesh->nodes[v].rotationalMass_I3 <
// 2);
// assert(std::pow(Dtini, 2.0) * rotationalSumSk_J /
// assert(std::pow(mSettings.Dtini, 2.0) * rotationalSumSk_J /
// mesh->nodes[v].rotationalMass_J <
// 2);
}
@ -1503,13 +1510,13 @@ void FormFinder::updatePositionsOnTheFly(
for (VertexType &v : mesh->vert) {
Node &node = mesh->nodes[v];
node.velocity = node.velocity + node.acceleration * settings.Dt;
node.velocity = node.velocity + node.acceleration * Dt;
}
updateKineticEnergy();
for (VertexType &v : mesh->vert) {
Node &node = mesh->nodes[v];
node.displacements = node.displacements + node.velocity * settings.Dt;
node.displacements = node.displacements + node.velocity * Dt;
}
for (VertexType &v : mesh->vert) {
@ -1636,45 +1643,43 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("Internal force norm", internalForcesNorm)
->setEnabled(true);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("First bending force-Translation",
internalFirstBendingTranslationForces)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("First bending force-Rotation",
internalFirstBendingRotationForces)
->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Internal Axial force", internalAxialForces)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("First bending force-Translation",
// internalFirstBendingTranslationForces)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("First bending force-Rotation",
// internalFirstBendingRotationForces)
// ->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("Second bending force-Translation",
internalSecondBendingTranslationForces)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("Second bending force-Rotation",
internalSecondBendingRotationForces);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->getQuantity("Second bending force-Rotation")
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("nR", nRs)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("theta3", theta3)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("theta2", theta2)
->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Second bending force-Translation",
// internalSecondBendingTranslationForces)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeVectorQuantity("Second bending force-Rotation",
// internalSecondBendingRotationForces)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeScalarQuantity("nR", nRs)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeScalarQuantity("theta3", theta3)
// ->setEnabled(false);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeScalarQuantity("theta2", theta2)
// ->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeVectorQuantity("Residual force", residualForces)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("Residual force norm", residualForcesNorm)
->setEnabled(false);
polyscope::getCurveNetwork(meshPolyscopeLabel)
->addNodeScalarQuantity("Node acceleration x", accelerationX);
// polyscope::getCurveNetwork(meshPolyscopeLabel)
// ->addNodeScalarQuantity("Node acceleration x", accelerationX);
// Edge quantities
std::vector<double> A(mesh->EN());
@ -1706,9 +1711,9 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
// matching PopItemWidth() below.
ImGui::InputInt("Simulation drawing step",
&settings.drawingStep); // set a int variable
&mSettings.drawingStep); // set a int variable
ImGui::Checkbox("Enable drawing",
&settings.shouldDraw); // set a int variable
&mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
ImGui::PopItemWidth();
@ -1733,25 +1738,35 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
}
}
void FormFinder::printCurrentState() {
std::cout << "Simulation steps executed:" << mCurrentSimulationStep
<< std::endl;
std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm
<< std::endl;
std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
<< std::endl;
}
void FormFinder::printDebugInfo() const {
std::cout << mesh->elements[0].rigidity.toString() << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
}
SimulationResults
FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const bool &shouldDraw, const bool &beVerbose,
const bool &shouldCreatePlots,
const size_t &maxDRMIterations) {
const Settings &settings) {
assert(pJob->pMesh.operator bool());
auto t1 = std::chrono::high_resolution_clock::now();
reset();
settings.shouldDraw = shouldDraw;
settings.beVerbose = beVerbose;
mSettings = settings;
// if(!job.nodalExternalForces.empty()){
// double externalForcesNorm=0;
// for(const auto& externalForce:job.nodalExternalForces)
// {
// externalForcesNorm+=externalForce.norm();
// }
// setTotalResidualForcesNormThreshold(externalForcesNorm);
// }
if (!pJob->nodalExternalForces.empty()) {
double externalForcesNorm = 0;
for (const auto &externalForce : pJob->nodalExternalForces) {
externalForcesNorm += externalForce.second.norm();
}
mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-3;
}
constrainedVertices = pJob->constrainedVertices;
if (!pJob->nodalForcedDisplacements.empty()) {
@ -1764,24 +1779,24 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
// if (!pJob->nodalForcedNormals.empty()) {
// for (std::pair<VertexIndex, Eigen::Vector3d> viNormalPair :
// pJob->nodalForcedDisplacements) {
// assert(viNormalPair.second[2] >= 0);
// assert(viNormalPair3second[2] >= 0);
// }
// }
mesh = std::make_unique<SimulationMesh>(*pJob->pMesh);
if (beVerbose) {
if (mSettings.beVerbose) {
std::cout << "Executing simulation for mesh with:" << mesh->VN()
<< " nodes and " << mesh->EN() << " elements." << std::endl;
}
computeRigidSupports();
if (settings.shouldDraw || true) {
if (mSettings.shouldDraw || true) {
if (!polyscope::state::initialized) {
initPolyscope();
}
polyscope::registerCurveNetwork(
meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges());
// registerWorldAxes();
registerWorldAxes();
}
for (auto fixedVertex : pJob->constrainedVertices) {
assert(fixedVertex.first < mesh->VN());
@ -1821,35 +1836,38 @@ FormFinder::executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
mesh->previousTotalPotentialEnergykN = mesh->currentTotalPotentialEnergykN;
mesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
if (mCurrentSimulationStep != 0) {
history.stepPulse(*mesh);
}
if (mCurrentSimulationStep > 100000) {
shouldUseKineticEnergyThreshold = true;
}
if (mCurrentSimulationStep > 500000) {
if (std::isnan(mesh->currentTotalKineticEnergy)) {
std::time_t now = time(0);
std::tm *ltm = std::localtime(&now);
std::string dir =
"../ProblematicSimulationJobs/" + std::to_string(ltm->tm_hour) + ":" +
std::to_string(ltm->tm_min) + "_" + std::to_string(ltm->tm_mday) +
"_" + std::to_string(1 + ltm->tm_mon) + "_" +
std::to_string(1900 + ltm->tm_year);
std::filesystem::create_directories(dir);
pJob->save(dir);
std::string dir = "../ProblematicSimulationJobs/" +
std::to_string(ltm->tm_mday) + "_" +
std::to_string(1 + ltm->tm_mon) + "_" +
std::to_string(1900 + ltm->tm_year);
const std::string subDir = std::filesystem::path(dir)
.append(std::to_string(ltm->tm_hour) +
":" + std::to_string(ltm->tm_min))
.string();
std::filesystem::create_directories(subDir);
pJob->save(subDir);
std::cout << "Non terminating simulation found. Saved simulation job to:"
<< dir << std::endl;
std::cout << "Exiting.." << std::endl;
FormFinder debug;
debug.executeSimulation(pJob, true, true, true);
std::terminate();
// FormFinder debug;
// debug.executeSimulation(pJob, true, true, true);
// std::terminate();
break;
}
if (settings.shouldDraw &&
mCurrentSimulationStep % settings.drawingStep == 0) /* &&
if (mSettings.shouldDraw &&
mCurrentSimulationStep % mSettings.drawingStep == 0) /* &&
currentSimulationStep > maxDRMIterations*/
{
// std::string saveTo = std::filesystem::current_path()
@ -1866,63 +1884,75 @@ currentSimulationStep > maxDRMIterations*/
// shouldUseKineticEnergyThreshold = true;
}
// if (mCurrentSimulationStep % 100 == 0 && mCurrentSimulationStep >
// 100000) {
// std::cout << "Current simulation step:" << mCurrentSimulationStep
// << std::endl;
// std::cout << "Residual forces norm: " <<
// mesh->totalResidualForcesNorm
// << std::endl;
// std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
// << std::endl;
// // auto yValues = history.residualForces;
// // auto xPlot = matplot::linspace(0, yValues.size(),
// yValues.size());
// // plotHandle = matplot::scatter(xPlot, yValues);
// // std::string label = "Residual forces";
// // plotHandle->legend_string(label);
// }
if (mSettings.shouldCreatePlots && mCurrentSimulationStep % 10 == 0) {
printCurrentState();
SimulationResultsReporter::createPlot(
"Number of Steps", "Log of Kinetic energy", history.kineticEnergy);
}
// t = t + Dt;
mCurrentSimulationStep++;
// std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
// << std::endl;
// std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm
// << std::endl;
if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) {
if (/*mesh.totalPotentialEnergykN < 10 ||*/
mesh->totalResidualForcesNorm <
settings.totalResidualForcesNormThreshold ||
(shouldUseKineticEnergyThreshold &&
mesh->currentTotalKineticEnergy <
settings.totalKineticEnergyThreshold)) {
if (beVerbose) {
std::cout << "Convergence after " << mCurrentSimulationStep
<< " steps" << std::endl;
std::cout << "Residual force of magnitude:"
<< mesh->previousTotalResidualForcesNorm << std::endl;
std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy
<< std::endl;
mesh->totalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
std::cout << "Total potential:" << mesh->totalPotentialEnergykN
// Kinetic energy convergence
if ((mSettings.shouldUseTranslationalKineticEnergyThreshold ||
mCurrentSimulationStep > 100000) &&
mesh->currentTotalTranslationalKineticEnergy <
mSettings.totalTranslationalKineticEnergyThreshold) {
if (mSettings.beVerbose) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
std::cout << "Total potential:" << mesh->currentTotalPotentialEnergykN
<< " kNm" << std::endl;
}
std::cout << "Warning: The kinetic energy of the system was "
" used as a convergence criterion"
<< std::endl;
break;
}
// Residual forces norm convergence
if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy
/*||
mesh->previousTotalPotentialEnergykN >
mesh->currentTotalPotentialEnergykN*/
/*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) {
if (mesh->totalResidualForcesNorm <
mSettings.totalResidualForcesNormThreshold) {
if (mSettings.beVerbose) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
std::cout << "Total potential:" << mesh->currentTotalPotentialEnergykN
<< " kNm" << std::endl;
}
if (shouldUseKineticEnergyThreshold) {
std::cout << "Warning: The kinetic energy of the system was "
" used as a convergence criterion"
<< std::endl;
}
break;
// }
}
// for (VertexType &v : mesh->vert) {
// Node &node = mesh->nodes[v];
// node.displacements = node.displacements - node.velocity * Dt;
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(
// pJob->nodalForcedDisplacements);
// }
// updateElementalLengths();
resetVelocities();
settings.Dt = settings.Dt * settings.xi;
Dt = Dt * mSettings.xi;
++numOfDampings;
}
if (debug) {
printDebugInfo();
}
}
if (mCurrentSimulationStep == maxDRMIterations) {
std::cout << "Did not reach equilibrium before reaching the maximum number "
"of DRM steps ("
<< maxDRMIterations << "). Breaking simulation" << std::endl;
} else if (beVerbose) {
} else if (mSettings.beVerbose) {
auto t2 = std::chrono::high_resolution_clock::now();
double simulationDuration =
std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
@ -1933,18 +1963,19 @@ currentSimulationStep > maxDRMIterations*/
<< simulationDuration /
(static_cast<double>(mCurrentSimulationStep) * mesh->VN())
<< "s" << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
}
// mesh.printVertexCoordinates(mesh.VN() / 2);
if (settings.shouldDraw) {
SimulationResults results = computeResults(pJob);
if (mSettings.shouldCreatePlots) {
SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}
if (mSettings.shouldDraw) {
draw();
}
if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) {
polyscope::removeCurveNetwork(meshPolyscopeLabel);
}
SimulationResults results = computeResults(pJob);
if (shouldCreatePlots) {
SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}
return results;
}

33
beamformfinder.hpp
View File

@ -21,27 +21,32 @@ struct DifferentiateWithRespectTo {
const VertexType &v;
const DoFType &dofi;
};
class FormFinder {
public:
struct Settings {
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
int drawingStep{1};
const double totalKineticEnergyThreshold{1e-9};
double totalTranslationalKineticEnergyThreshold{1e-7};
double totalResidualForcesNormThreshold{1e-3};
const double Dtini{0.000001};
double Dt{Dtini};
const double xi{0.9969};
bool beVerbose;
double Dtini{0.1};
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
Settings() {}
};
Settings settings;
private:
bool shouldUseKineticEnergyThreshold{false};
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
const std::string meshPolyscopeLabel{"Simulation mesh"};
std::unique_ptr<SimulationMesh> mesh;
@ -189,13 +194,15 @@ private:
void applyForcedNormals(
const std::unordered_map<VertexIndex, VectorType> nodalForcedRotations);
void printCurrentState();
void printDebugInfo() const;
public:
FormFinder();
SimulationResults executeSimulation(
const std::shared_ptr<SimulationJob> &pJob,
const bool &shouldDraw = false, const bool &beVerbose = false,
const bool &createPlots = false,
const size_t &maxDRMIterations = FormFinder::maxDRMIterations);
SimulationResults
executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings());
inline static const size_t maxDRMIterations{0};
static void runUnitTests();

3
edgemesh.cpp
View File

@ -174,7 +174,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth,
return true;
}
bool VCGEdgeMesh::loadFromPly(const std::string plyFilename) {
bool VCGEdgeMesh::loadPly(const std::string plyFilename) {
std::string usedPath = plyFilename;
if (std::filesystem::path(plyFilename).is_relative()) {
@ -352,5 +352,4 @@ void VCGEdgeMesh::registerForDrawing() const {
const double drawingRadius = 0.0007;
polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges())
->setRadius(drawingRadius, false);
std::cout << "Registered:" << label << std::endl;
}

2
edgemesh.hpp
View File

@ -66,7 +66,7 @@ public:
bool loadUsingNanoply(const std::string &plyFilename);
bool loadFromPly(const std::string plyFilename);
virtual bool loadPly(const std::string plyFilename);
bool savePly(const std::string plyFilename);

20
elementalmesh.cpp
View File

@ -141,7 +141,7 @@ void SimulationMesh::reset() {
const vcg::Segment3<double> s(p0, p1);
element.initialLength = s.Length();
element.length = element.initialLength;
element.updateConstFactors();
element.updateRigidity();
}
for (const VertexType &v : vert) {
@ -194,7 +194,7 @@ void SimulationMesh::initializeElements() {
element.initialLength = s.Length();
element.length = element.initialLength;
// Initialize const factors
element.updateConstFactors();
element.updateRigidity();
element.derivativeT1.resize(
2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
element.derivativeT2.resize(
@ -233,14 +233,14 @@ void SimulationMesh::setBeamCrossSection(
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties.dimensions = beamDimensions;
elements[ei].properties.computeDimensionsProperties(beamDimensions);
elements[ei].updateConstFactors();
elements[ei].updateRigidity();
}
}
void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties.setMaterial(ElementMaterial{pr, ym});
elements[ei].updateConstFactors();
elements[ei].updateRigidity();
}
}
@ -316,7 +316,7 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
if (!handleBeamProperties._handle->data.empty()) {
for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties = Element::Properties(handleBeamProperties[ei]);
elements[ei].updateConstFactors();
elements[ei].updateRigidity();
}
}
@ -469,9 +469,9 @@ std::array<double, 6> Element::Properties::toArray() const {
return std::array<double, 6>({E, G, A, I2, I3, J});
}
void Element::updateConstFactors() {
axialConstFactor = properties.E * properties.A / initialLength;
torsionConstFactor = properties.G * properties.J / initialLength;
firstBendingConstFactor = 2 * properties.E * properties.I2 / initialLength;
secondBendingConstFactor = 2 * properties.E * properties.I3 / initialLength;
void Element::updateRigidity() {
rigidity.axial = properties.E * properties.A / initialLength;
rigidity.torsional = properties.G * properties.J / initialLength;
rigidity.firstBending = 2 * properties.E * properties.I2 / initialLength;
rigidity.secondBending = 2 * properties.E * properties.I3 / initialLength;
}

30
elementalmesh.hpp
View File

@ -33,15 +33,17 @@ public:
std::vector<VCGEdgeMesh::EdgePointer>
getIncidentElements(const VCGEdgeMesh::VertexType &v);
bool loadPly(const string &plyFilename);
virtual bool loadPly(const string &plyFilename);
std::vector<CrossSectionType> getBeamDimensions();
std::vector<ElementMaterial> getBeamMaterial();
double previousTotalKineticEnergy{0};
double previousTotalResidualForcesNorm{0};
double currentTotalKineticEnergy{0};
double currentTotalTranslationalKineticEnergy{0};
double totalResidualForcesNorm{0};
double totalPotentialEnergykN{0};
double currentTotalPotentialEnergykN{0};
double previousTotalPotentialEnergykN{0};
bool savePly(const std::string &plyFilename);
void setBeamCrossSection(const CrossSectionType &beamDimensions);
void setBeamMaterial(const double &pr, const double &ym);
@ -50,6 +52,7 @@ public:
};
struct Element {
struct Properties {
CrossSectionType dimensions;
ElementMaterial material;
@ -81,17 +84,28 @@ struct Element {
VectorType t3;
};
void updateConstFactors();
EdgeIndex ei;
double length{0};
Properties properties;
double initialLength;
LocalFrame frame;
double axialConstFactor;
double torsionConstFactor;
double firstBendingConstFactor;
double secondBendingConstFactor;
struct Rigidity {
double axial;
double torsional;
double firstBending;
double secondBending;
std::string toString() const {
return std::string("Rigidity:") + std::string("\nAxial=") +
std::to_string(axial) + std::string("\nTorsional=") +
std::to_string(torsional) + std::string("\nFirstBending=") +
std::to_string(firstBending) + std::string("\nSecondBending=") +
std::to_string(secondBending);
}
};
Rigidity rigidity;
void updateRigidity();
VectorType f1_j;
VectorType f1_jplus1;
VectorType f2_j;

12
flatpattern.cpp
View File

@ -5,8 +5,16 @@
FlatPattern::FlatPattern() {}
FlatPattern::FlatPattern(const std::string &filename, bool addNormalsIfAbsent) {
assert(std::filesystem::exists(std::filesystem::path(filename)));
loadFromPly(filename);
if (!std::filesystem::exists(std::filesystem::path(filename))) {
assert(false);
std::cerr << "No flat pattern with name " << filename << std::endl;
return;
}
if (!loadPly(filename)) {
assert(false);
std::cerr << "File could not be loaded " << filename << std::endl;
return;
}
if (addNormalsIfAbsent) {
bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001;
if (normalsAreAbsent) {

6
simulationhistoryplotter.hpp
View File

@ -85,15 +85,17 @@ struct SimulationResultsReporter {
const std::string &saveTo = {}) {
matplot::xlabel(xLabel);
matplot::ylabel(yLabel);
matplot::figure(true);
// matplot::figure(true);
// matplot::hold(matplot::on);
matplot::grid(matplot::on);
auto x = matplot::linspace(0, YvaluesToPlot.size(), YvaluesToPlot.size());
matplot::scatter(x, YvaluesToPlot)
// ->marker_indices(history.redMarks)
// ->marker_indices(truncatedRedMarks)
// .marker_color("r")
->marker_size(1);
// ->marker_size(1)
;
// auto greenMarksY = matplot::transform(
// history.greenMarks, [&](auto x) { return history.kineticEnergy[x];
// });