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

Refactoring

This commit is contained in:
iasonmanolas 2021-05-24 14:43:32 +03:00
parent 6e682e1e71
commit 8793dc5408
12 changed files with 511 additions and 267 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
csvfile.hpp
View File

@ -1,10 +1,13 @@
#ifndef CSVFILE_HPP
#define CSVFILE_HPP
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
class csvFile;
@ -60,7 +63,46 @@ public:
csvFile &operator<<(const std::string &val) { return write(escape(val)); }
template <typename T> csvFile &operator<<(const T &val) { return write(val); }
template<typename T>
csvFile &operator<<(const T &val)
{
return write(val);
}
template<typename T>
static std::vector<std::vector<T>> parse(const std::filesystem::path &csvFilepath)
{
std::vector<std::vector<T>> resultCSV;
if (!std::filesystem::exists(csvFilepath)) {
std::cerr << "The file does not exist:" << csvFilepath.string() << std::endl;
return resultCSV;
}
std::ifstream inputfile(csvFilepath.string().c_str());
if (!inputfile.is_open()) {
std::cerr << "Can't open file:" << csvFilepath.string() << std::endl;
return resultCSV;
}
std::vector<T> row;
std::string line;
using Tokenizer = boost::tokenizer<boost::escaped_list_separator<char>>;
while (std::getline(inputfile, line)) {
Tokenizer tokenizer(line);
row.resize(std::distance(tokenizer.begin(), tokenizer.end()));
std::transform(tokenizer.begin(), tokenizer.end(), row.begin(), [](const std::string &el) {
return boost::lexical_cast<T>(el);
});
// std::cout << std::endl;
// row.assign(tokenizer.begin(), tokenizer.end());
// for (const auto &el : row) {
// std::cout << el << " ";
// }
// std::cout << std::endl;
resultCSV.push_back(row);
}
return resultCSV;
}
private:
template <typename T> csvFile &write(const T &val) {

486
drmsimulationmodel.cpp
View File

@ -22,7 +22,7 @@ void DRMSimulationModel::runUnitTests()
// mesh.createSpanGrid(spanGridSize);
beam.load(std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string());
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3};
fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2};
fixedVertices[beam.VN() - 1] = std::unordered_set<DoFType>{1, 2};
std::unordered_map<VertexIndex, Vector6d> nodalForces{
{beam.VN() / 2, Vector6d({0, 1000, 1000, 0, 0, 0})}};
@ -45,9 +45,10 @@ void DRMSimulationModel::runUnitTests()
settings.Dtini = 0.1;
settings.xi = 0.9969;
settings.maxDRMIterations = 0;
totalResidualForcesNormThreshold = 1;
// settings.shouldDraw = true;
settings.beVerbose = true;
formFinder.totalResidualForcesNormThreshold = 1;
settings.shouldDraw = false;
settings.beVerbose = false;
// settings.shouldCreatePlots = true;
SimulationResults simpleBeam_simulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>(beamSimulationJob), settings);
simpleBeam_simulationResults.save();
@ -58,10 +59,19 @@ void DRMSimulationModel::runUnitTests()
assert(std::filesystem::exists(std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename)));
Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements;
Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename, simpleBeam_groundOfTruthDisplacements);
if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements)) {
std::cerr << "Error:Beam simulation produces wrong results." << std::endl;
double error;
if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements, error)) {
std::cerr << "Error:Beam simulation produces wrong results. Error:" << error << std::endl;
// return;
}
#ifdef POLYSCOPE_DEFINED
beam.registerForDrawing();
simpleBeam_simulationResults.registerForDrawing();
polyscope::show();
beam.unregister();
simpleBeam_simulationResults.unregister();
#endif
// Second example of the paper
VCGEdgeMesh shortSpanGrid;
@ -109,13 +119,21 @@ void DRMSimulationModel::runUnitTests()
// shortSpanGridshellSimulationResults.registerForDrawing(
// shortSpanGridshellSimulationJob);
// polyscope::show();
assert(
shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements));
if (!shortSpanGridshellSimulationResults.isEqual(
shortSpanGridshell_groundOfTruthDisplacements)) {
std::cerr << "Error:Short span simulation produces wrong results." << std::endl;
return;
assert(shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements,
error));
if (!shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements,
error)) {
std::cerr << "Error:Short span simulation produces wrong results. Error:" << error
<< std::endl;
// return;
}
#ifdef POLYSCOPE_DEFINED
shortSpanGrid.registerForDrawing();
shortSpanGridshellSimulationResults.registerForDrawing();
polyscope::show();
shortSpanGrid.unregister();
shortSpanGridshellSimulationResults.unregister();
#endif
// Third example of the paper
VCGEdgeMesh longSpanGrid;
@ -163,6 +181,7 @@ void DRMSimulationModel::runUnitTests()
<< std::endl;
}
longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
SimulationResults longSpanGridshell_simulationResults
= formFinder.executeSimulation(std::make_shared<SimulationJob>(
longSpanGridshell_simulationJob),
@ -177,12 +196,21 @@ void DRMSimulationModel::runUnitTests()
Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements;
Eigen::read_binary(longSpan_groundOfTruthBinaryFilename,
longSpanGridshell_groundOfTruthDisplacements);
assert(
longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements));
if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements)) {
std::cerr << "Error:Long span simulation produces wrong results." << std::endl;
return;
assert(longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements,
error));
if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements,
error)) {
std::cerr << "Error:Long span simulation produces wrong results. Error:" << error
<< std::endl;
// return;
}
#ifdef POLYSCOPE_DEFINED
longSpanGrid.registerForDrawing();
longSpanGridshell_simulationResults.registerForDrawing();
polyscope::show();
longSpanGrid.unregister();
longSpanGridshell_simulationResults.unregister();
#endif
std::cout << "Form finder unit tests succesufully passed." << std::endl;
@ -1130,14 +1158,11 @@ void DRMSimulationModel::updateNodalMasses()
const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3);
const long double SkRotational_I2
= elem.material.youngsModulus * elem.I2
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const long double SkRotational_I3
= elem.material.youngsModulus * elem.I3
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const long double SkRotational_J = elem.material.youngsModulus * elem.J
/ lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const long double SkRotational_I2 = elem.material.youngsModulus * elem.I2
/ lengthToThe3;
const long double SkRotational_I3 = elem.material.youngsModulus * elem.I3
/ lengthToThe3;
const long double SkRotational_J = elem.material.youngsModulus * elem.J / lengthToThe3;
rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3;
rotationalSumSk_J += SkRotational_J;
@ -1336,9 +1361,8 @@ void DRMSimulationModel::updateElementalFrames()
}
void DRMSimulationModel::applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements)
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements)
{
const int gradualDisplacementSteps = 500;
for (const std::pair<VertexIndex, Eigen::Vector3d> vertexIndexDisplacementPair :
nodalForcedDisplacements) {
const VertexIndex vi = vertexIndexDisplacementPair.first;
@ -1347,9 +1371,9 @@ void DRMSimulationModel::applyForcedDisplacements(
VectorType displacementVector(vertexDisplacement(0),
vertexDisplacement(1),
vertexDisplacement(2));
if (mCurrentSimulationStep < gradualDisplacementSteps) {
if (mCurrentSimulationStep < mSettings.gradualForcedDisplacementSteps) {
displacementVector *= mCurrentSimulationStep
/ static_cast<double>(gradualDisplacementSteps);
/ static_cast<double>(mSettings.gradualForcedDisplacementSteps);
}
pMesh->vert[vi].P() = node.initialLocation + displacementVector;
node.displacements = Vector6d({displacementVector[0],
@ -1553,6 +1577,13 @@ void DRMSimulationModel::printCurrentState() const
std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl;
std::cout << "Residual forces norm: " << pMesh->totalResidualForcesNorm << 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>(
std::chrono::steady_clock::now() - begin)
.count()
* 1e-6 / (mCurrentSimulationStep * pMesh->VN());
std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl;
std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << std::endl;
}
void DRMSimulationModel::printDebugInfo() const
@ -1694,22 +1725,26 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder)
->addEdgeScalarQuantity("I3", I3);
// Specify the callback
PolyscopeInterface::addUserCallback([&]() {
// Since options::openImGuiWindowForUserCallback == true by default,
// we can immediately start using ImGui commands to build a UI
static bool calledOnce = false;
if (!calledOnce) {
PolyscopeInterface::addUserCallback([&]() {
// Since options::openImGuiWindowForUserCallback == true by default,
// we can immediately start using ImGui commands to build a UI
ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide,
// instead of full width. Must have
// matching PopItemWidth() below.
ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide,
// instead of full width. Must have
// matching PopItemWidth() below.
ImGui::InputInt("Simulation drawing step",
&mSettings.drawingStep); // set a int variable
ImGui::Checkbox("Enable drawing",
&mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
ImGui::InputInt("Simulation drawing step",
&mSettings.drawingStep); // set a int variable
ImGui::Checkbox("Enable drawing",
&mSettings.shouldDraw); // set a int variable
ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep);
ImGui::PopItemWidth();
});
ImGui::PopItemWidth();
});
calledOnce = true;
}
if (!screenshotsFolder.empty()) {
static bool firstDraw = true;
@ -1742,16 +1777,27 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
for (const std::pair<VertexIndex, Vector6d> &externalForce : pJob->nodalExternalForces) {
totalExternalForcesNorm += externalForce.second.norm();
}
totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm;
// if (!pJob->nodalExternalForces.empty()) {
// double externalForcesNorm = 0;
// for (const auto &externalForce : pJob->nodalExternalForces) {
// externalForcesNorm += externalForce.second.norm();
// }
// mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2;
// }
if (!pJob->nodalExternalForces.empty()) {
totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination
* totalExternalForcesNorm;
} else {
totalResidualForcesNormThreshold = 1;
std::cout << "No forces setted default residual forces norm threshold" << std::endl;
}
if (mSettings.beVerbose) {
std::cout << "totalResidualForcesNormThreshold:" << totalResidualForcesNormThreshold
<< std::endl;
}
// if (!pJob->nodalExternalForces.empty()) {
// double externalForcesNorm = 0;
// for (const auto &externalForce : pJob->nodalExternalForces) {
// externalForcesNorm += externalForce.second.norm();
// }
// mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2;
// }
constrainedVertices = pJob->constrainedVertices;
if (!pJob->nodalForcedDisplacements.empty()) {
@ -1773,18 +1819,20 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and "
<< pMesh->EN() << " elements." << std::endl;
}
vcg::tri::UpdateBounding<SimulationMesh>::Box(*pMesh);
computeRigidSupports();
for (auto fixedVertex : pJob->constrainedVertices) {
assert(fixedVertex.first < pMesh->VN());
}
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw ) {
PolyscopeInterface::init();
polyscope::registerCurveNetwork(
meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges());
// registerWorldAxes();
}
if (mSettings.shouldDraw || mSettings.isDebugMode) {
PolyscopeInterface::init();
polyscope::registerCurveNetwork(meshPolyscopeLabel,
pMesh->getEigenVertices(),
pMesh->getEigenEdges());
// registerWorldAxes();
}
#endif
updateElementalFrames();
@ -1793,130 +1841,203 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr<Si
shouldApplyInitialDistortion = true;
}
updateNodalExternalForces(pJob->nodalExternalForces, constrainedVertices);
while (settings.maxDRMIterations == 0 ||
mCurrentSimulationStep < settings.maxDRMIterations) {
// while (true) {
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
updateResidualForcesOnTheFly(constrainedVertices);
const size_t movingAverageSampleSize = 50;
std::queue<double> residualForcesMovingAverageHistorySample;
// TODO: write parallel function for updating positions
// TODO: Make a single function out of updateResidualForcesOnTheFly
// updatePositionsOnTheFly
// updatePositionsOnTheFly(constrainedVertices);
updateNodalMasses();
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
double residualForcesMovingAverageDerivativeMax = 0;
while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) {
// while (true) {
updateNormalDerivatives();
updateT1Derivatives();
updateRDerivatives();
updateT2Derivatives();
updateT3Derivatives();
updateResidualForcesOnTheFly(constrainedVertices);
pJob->nodalForcedDisplacements);
}
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// TODO: write parallel function for updating positions
// TODO: Make a single function out of updateResidualForcesOnTheFly
// updatePositionsOnTheFly
// updatePositionsOnTheFly(constrainedVertices);
updateNodalMasses();
updateNodalAccelerations();
updateNodalVelocities();
updateNodalDisplacements();
applyDisplacements(constrainedVertices);
if (!pJob->nodalForcedDisplacements.empty()) {
applyForcedDisplacements(
// if (mCurrentSimulationStep != 0) {
// history.stepPulse(*pMesh);
// }
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
pJob->nodalForcedDisplacements);
}
// if (!pJob->nodalForcedNormals.empty()) {
// applyForcedNormals(pJob->nodalForcedNormals);
// }
updateElementalLengths();
mCurrentSimulationStep++;
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;
// 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;
}
// pMesh->previousTotalPotentialEnergykN =
// pMesh->currentTotalPotentialEnergykN;
// pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000;
// timePerNodePerIterationHistor.push_back(timePerNodePerIteration);
if (mSettings.beVerbose && mSettings.isDebugMode
&& mCurrentSimulationStep % mSettings.debugModeStep == 0) {
printCurrentState();
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver",
// movingAverages);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Residual Forces mov aver deriv",
// movingAveragesDerivatives);
// draw();
std::cout << "Mov aver deriv:" << pMesh->residualForcesMovingAverageDerivativeNorm
<< std::endl;
// SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps",
// "Time/(#nodes*#iterations)",
// timePerNodePerIterationHistory);
}
if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) {
history.stepPulse(*pMesh);
}
if (std::isnan(pMesh->currentTotalKineticEnergy)) {
if (mSettings.beVerbose) {
std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl;
}
break;
}
break;
}
#ifdef POLYSCOPE_DEFINED
if (mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) /* &&
if ((mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0)
|| (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0)) /* &&
currentSimulationStep > maxDRMIterations*/
{
// std::string saveTo = std::filesystem::current_path()
// .append("Debugging_files")
// .append("Screenshots")
// .string();
draw();
// yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy";
// plotHandle->legend_string(label);
{
// std::string saveTo = std::filesystem::current_path()
// .append("Debugging_files")
// .append("Screenshots")
// .string();
draw();
// yValues = history.kineticEnergy;
// plotHandle = matplot::scatter(xPlot, yValues);
// label = "Log of Kinetic energy";
// plotHandle->legend_string(label);
// shouldUseKineticEnergyThreshold = true;
}
// shouldUseKineticEnergyThreshold = true;
}
#endif
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;
// Kinetic energy convergence
if ((mSettings.shouldUseTranslationalKineticEnergyThreshold || mCurrentSimulationStep > 500000)
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold) {
if (mSettings.beVerbose) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm"
<< std::endl;
}
#ifdef POLYSCOPE_DEFINED
std::cout << "Warning: The kinetic energy of the system was "
" used as a convergence criterion"
<< std::endl;
#endif
break;
}
// Residual forces norm convergence
if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
/*||
if (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) {
// 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",
// "Residual Forces",
// history.residualForces);
// SimulationResultsReporter::createPlot("Number of Steps",
// "Log of Kinetic energy",
// history.kineticEnergy);
}
// t = t + Dt;
// std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy
// << std::endl;
// std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm
// << std::endl;
// Residual forces norm convergence
if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy
&& mCurrentSimulationStep > movingAverageSampleSize
&& (pJob->nodalForcedDisplacements.empty()
|| mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps)
/*||
mesh->previousTotalPotentialEnergykN >
mesh->currentTotalPotentialEnergykN*/
/*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) {
if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) {
if (mSettings.beVerbose) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm"
<< std::endl;
}
break;
// }
}
// printCurrentState();
// for (VertexType &v : mesh->vert) {
// Node &node = mesh->nodes[v];
// node.displacements = node.displacements - node.velocity * Dt;
// }
// applyDisplacements(constrainedVertices);
// if (!pJob->nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(
/*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) {
// if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) {
const bool fullfillsKineticEnergyTerminationCriterion
= mSettings.shouldUseTranslationalKineticEnergyThreshold
&& pMesh->currentTotalTranslationalKineticEnergy
< mSettings.totalTranslationalKineticEnergyThreshold
&& mCurrentSimulationStep > 20 && numOfDampings > 2;
const bool fullfillsResidualForcesNormTerminationCriterion
= pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold;
const bool fullfillsMovingAverageNormTerminationCriterion
= pMesh->residualForcesMovingAverage
< mSettings.residualForcesMovingAverageNormThreshold;
/*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/
const bool shouldTerminate
= fullfillsKineticEnergyTerminationCriterion
// || fullfillsMovingAverageNormTerminationCriterion
|| fullfillsResidualForcesNormTerminationCriterion;
if (shouldTerminate) {
if (mSettings.beVerbose) {
std::cout << "Simulation converged." << std::endl;
printCurrentState();
if (fullfillsResidualForcesNormTerminationCriterion) {
std::cout << "Converged using residual forces norm threshold criterion"
<< std::endl;
} else if (fullfillsKineticEnergyTerminationCriterion) {
std::cout << "Warning: The kinetic energy of the system was "
" used as a convergence criterion"
<< std::endl;
} else if (fullfillsMovingAverageNormTerminationCriterion) {
std::cout << "Converged using residual forces moving average derivative norm "
"threshold criterion"
<< std::endl;
}
}
break;
// }
}
// printCurrentState();
// 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();
Dt = Dt * mSettings.xi;
++numOfDampings;
}
if (mSettings.debugMessages) {
printDebugInfo();
}
// pJob->nodalForcedDisplacements);
// }
// updateElementalLengths();
resetVelocities();
Dt = Dt * mSettings.xi;
++numOfDampings;
// std::cout << "Number of dampings:" << numOfDampings << endl;
}
}
SimulationResults results = computeResults(pJob);
@ -1929,20 +2050,21 @@ mesh->currentTotalPotentialEnergykN*/
<< std::endl;
results.converged = false;
} else if (std::isnan(pMesh->currentTotalKineticEnergy)) {
results.converged = false;
std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl;
results.converged = false;
} else if (mSettings.beVerbose) {
auto t2 = std::chrono::high_resolution_clock::now();
double simulationDuration =
std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
simulationDuration /= 1000;
std::cout << "Simulation converged after " << simulationDuration << "s"
<< std::endl;
std::cout << "Time/(node*iteration) "
<< simulationDuration /
(static_cast<double>(mCurrentSimulationStep) * pMesh->VN())
<< "s" << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
auto t2 = std::chrono::high_resolution_clock::now();
double simulationDuration = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1)
.count();
simulationDuration /= 1000;
std::cout << "Simulation converged after " << simulationDuration << "s" << std::endl;
std::cout << "Time/(node*iteration) "
<< simulationDuration / (static_cast<double>(mCurrentSimulationStep) * pMesh->VN())
<< "s" << std::endl;
std::cout << "Number of dampings:" << numOfDampings << endl;
std::cout << "Number of steps:" << mCurrentSimulationStep << endl;
results.converged = true;
}
// mesh.printVertexCoordinates(mesh.VN() / 2);
#ifdef POLYSCOPE_DEFINED
@ -1994,7 +2116,7 @@ mesh->currentTotalPotentialEnergykN*/
}
}
if (mSettings.shouldCreatePlots) {
if (mSettings.shouldCreatePlots && results.converged) {
SimulationResultsReporter reporter;
reporter.reportResults({results}, "Results", pJob->pMesh->getLabel());
}

43
drmsimulationmodel.hpp
View File

@ -23,28 +23,33 @@ struct DifferentiateWithRespectTo {
class DRMSimulationModel
{
public:
struct Settings {
bool debugMessages{false};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
int drawingStep{1};
double totalTranslationalKineticEnergyThreshold{1e-10};
double totalExternalForcesNormPercentageTermination{1e-3};
double Dtini{0.1};
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
Settings() {}
};
struct Settings
{
bool isDebugMode{false};
int debugModeStep{10000};
bool shouldDraw{false};
bool beVerbose{false};
bool shouldCreatePlots{false};
int drawingStep{1};
double totalTranslationalKineticEnergyThreshold{1e-10};
double totalExternalForcesNormPercentageTermination{1e-3};
double residualForcesMovingAverageDerivativeNormThreshold{1e-8};
double residualForcesMovingAverageNormThreshold{1e-8};
double Dtini{0.1};
double xi{0.9969};
int maxDRMIterations{0};
bool shouldUseTranslationalKineticEnergyThreshold{false};
int gradualForcedDisplacementSteps{100};
Settings() {}
};
double totalResidualForcesNormThreshold{1};
private:
private:
Settings mSettings;
double Dt{mSettings.Dtini};
bool checkedForMaximumMoment{false};
double externalMomentsNorm{0};
size_t mCurrentSimulationStep{0};
double totalResidualForcesNormThreshold{1e-3};
matplot::line_handle plotHandle;
std::vector<double> plotYValues;
size_t numOfDampings{0};
@ -80,8 +85,7 @@ private:
void updateNodalDisplacements();
void applyForcedDisplacements(
const std::unordered_map<VertexIndex, Eigen::Vector3d>
nodalForcedDisplacements);
const std::unordered_map<VertexIndex, Eigen::Vector3d> &nodalForcedDisplacements);
Vector6d computeElementTorsionalForce(
const Element &element, const Vector6d &displacementDifference,
@ -208,8 +212,7 @@ private:
executeSimulation(const std::shared_ptr<SimulationJob> &pJob,
const Settings &settings = Settings());
void runUnitTests();
void setTotalResidualForcesNormThreshold(double value);
static void runUnitTests();
};
template <typename PointType> PointType Cross(PointType p1, PointType p2) {

3
edgemesh.cpp
View File

@ -200,7 +200,6 @@ bool VCGEdgeMesh::load(const string &plyFilename) {
getEdges(eigenEdges);
getVertices(eigenVertices);
vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
label = std::filesystem::path(plyFilename).stem().string();
const bool printDebugInfo = false;
if (printDebugInfo) {
@ -208,7 +207,7 @@ bool VCGEdgeMesh::load(const string &plyFilename) {
std::cout << "Mesh has " << EN() << " edges." << std::endl;
}
label=std::filesystem::path(plyFilename).stem();
label = std::filesystem::path(plyFilename).stem().string();
return true;
}

8
reducedmodeloptimizer_structs.hpp
View File

@ -126,12 +126,14 @@ struct Colors
}
xRange x;
x.fromString(json.at(jsonXRangeKey));
xRanges.push_back(x);
xRangeIndex++;
}
numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls);
solverAccuracy = json.at(JsonKeys::SolverAccuracy);
objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight);
objectiveWeights.rotational = 2 - objectiveWeights.translational;
return true;
}
@ -197,6 +199,7 @@ struct Colors
inline void updateMeshWithOptimalVariables(const std::vector<double> &x, SimulationMesh &m)
{
assert(CrossSectionType::name == RectangularBeamDimensions::name);
const double E = x[0];
const double A = x[1];
const double beamWidth = std::sqrt(A);
@ -207,7 +210,7 @@ struct Colors
const double I3 = x[4];
for (EdgeIndex ei = 0; ei < m.EN(); ei++) {
Element &e = m.elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
e.setDimensions(CrossSectionType(beamWidth, beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio, E));
e.J = J;
e.I2 = I2;
@ -494,6 +497,7 @@ struct Colors
const ReducedPatternOptimization::Results &reducedPattern_optimizationResults,
const shared_ptr<SimulationMesh> &pTiledReducedPattern_simulationMesh)
{
assert(CrossSectionType::name == RectangularBeamDimensions::name);
// Set reduced parameters fron the optimization results
std::unordered_map<std::string, double>
optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(),
@ -504,7 +508,7 @@ struct Colors
const double A = optimalXVariables.at(ALabel);
const double beamWidth = std::sqrt(A);
const double beamHeight = beamWidth;
RectangularBeamDimensions elementDimensions(beamWidth, beamHeight);
CrossSectionType elementDimensions(beamWidth, beamHeight);
const double poissonsRatio = 0.3;
const std::string ymLabel = "E";

42
simulation_structs.hpp
View File

@ -9,11 +9,14 @@ struct SimulationHistory {
size_t numberOfSteps{0};
std::string label;
std::vector<double> residualForces;
std::vector<double> residualForcesLog;
std::vector<double> kineticEnergy;
std::vector<double> potentialEnergies;
std::vector<size_t> redMarks;
std::vector<double> greenMarks;
std::vector<double> residualForcesMovingAverage;
std::vector<double> sumOfNormalizedDisplacementNorms;
// std::vector<double> residualForcesMovingAverageDerivativesLog;
void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); }
@ -22,13 +25,21 @@ struct SimulationHistory {
void stepPulse(const SimulationMesh &mesh) {
kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy));
// potentialEnergy.push_back(mesh.totalPotentialEnergykN);
residualForces.push_back(mesh.totalResidualForcesNorm);
residualForcesLog.push_back(std::log(mesh.totalResidualForcesNorm));
residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage);
sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms);
// residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
}
void clear() {
residualForces.clear();
kineticEnergy.clear();
potentialEnergies.clear();
void clear()
{
residualForcesLog.clear();
kineticEnergy.clear();
potentialEnergies.clear();
residualForcesMovingAverage.clear();
sumOfNormalizedDisplacementNorms.clear();
// residualForcesMovingAverageDerivativesLog.clear();
}
};
@ -399,6 +410,9 @@ public:
}
void unregister(const std::string &meshLabel)
{
if (polyscope::getCurveNetwork(meshLabel) == nullptr) {
return;
}
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
}
@ -478,23 +492,29 @@ struct SimulationResults
void save(const std::string &outputFolder = std::string())
{
const std::filesystem::path outputPath(outputFolder);
job->save(std::filesystem::path(outputPath).append("SimulationJob").string());
const std::filesystem::path outputFolderPath = outputFolder.empty()
? std::filesystem::current_path()
: std::filesystem::path(outputFolder);
std::filesystem::path simulationJobOutputFolderPath
= std::filesystem::path(outputFolderPath).append("SimulationJob");
std::filesystem::create_directories(simulationJobOutputFolderPath);
job->save(simulationJobOutputFolderPath);
const std::string filename(getLabel() + "_displacements.eigenBin");
Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements);
Eigen::write_binary(std::filesystem::path(outputPath).append(filename).string(), m);
Eigen::write_binary(std::filesystem::path(outputFolderPath).append(filename).string(), m);
saveDeformedModel(outputFolder);
saveDeformedModel(outputFolderPath);
}
// The comparison of the results happens comparing the 6-dof nodal
// displacements
bool isEqual(const Eigen::MatrixXd &nodalDisplacements)
bool isEqual(const Eigen::MatrixXd &nodalDisplacements, double &error)
{
assert(nodalDisplacements.cols() == 6);
Eigen::MatrixXd eigenDisplacements = Utilities::toEigenMatrix(this->displacements);
const double errorNorm = (eigenDisplacements - nodalDisplacements).norm();
error = errorNorm;
return errorNorm < 1e-10;
// return eigenDisplacements.isApprox(nodalDisplacements);
}

105
simulationhistoryplotter.hpp
View File

@ -38,11 +38,11 @@ struct SimulationResultsReporter {
file << "\n";
}
if (!history.residualForces.empty()) {
if (!history.residualForcesLog.empty()) {
file << "Residual forces"
<< "\n";
for (size_t step = 0; step < numberOfSteps; step++) {
file << history.residualForces[step] << "\n";
file << history.residualForcesLog[step] << "\n";
}
file << "\n";
}
@ -92,27 +92,31 @@ struct SimulationResultsReporter {
}
}
static void createPlot(const std::string &xLabel, const std::string &yLabel,
static void createPlot(const std::string &xLabel,
const std::string &yLabel,
const std::vector<double> &YvaluesToPlot,
const std::string &saveTo = {}) {
auto x =
matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo);
// matplot::figure(true);
// matplot::hold(matplot::on);
const std::string &saveTo = {})
{
if (YvaluesToPlot.size() < 2) {
return;
}
auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size());
createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo);
// matplot::figure(true);
// matplot::hold(matplot::on);
// ->marker_indices(history.redMarks)
// ->marker_indices(truncatedRedMarks)
// .marker_color("r")
// ->marker_size(1)
;
// auto greenMarksY = matplot::transform(
// history.greenMarks, [&](auto x) { return history.kineticEnergy[x];
// });
// matplot::scatter(history.greenMarks, greenMarksY)
// ->color("green")
// .marker_size(10);
// matplot::hold(matplot::off);
// ->marker_indices(history.redMarks)
// ->marker_indices(truncatedRedMarks)
// .marker_color("r")
// ->marker_size(1)
;
// auto greenMarksY = matplot::transform(
// history.greenMarks, [&](auto x) { return history.kineticEnergy[x];
// });
// matplot::scatter(history.greenMarks, greenMarksY)
// ->color("green")
// .marker_size(10);
// matplot::hold(matplot::off);
}
void createPlots(const SimulationHistory &history,
@ -124,27 +128,54 @@ struct SimulationResultsReporter {
std::filesystem::create_directory(graphsFolder.string());
if (!history.kineticEnergy.empty()) {
createPlot("Number of Iterations", "Log of Kinetic Energy",
history.kineticEnergy,
std::filesystem::path(graphsFolder)
.append("KineticEnergy" + graphSuffix + ".png")
.string());
createPlot("Number of Iterations",
"Log of Kinetic Energy log",
history.kineticEnergy,
std::filesystem::path(graphsFolder)
.append("KineticEnergyLog_" + graphSuffix + ".png")
.string());
}
if (!history.residualForces.empty()) {
createPlot("Number of Iterations", "Residual Forces norm",
history.residualForces,
std::filesystem::path(graphsFolder)
.append("ResidualForces" + graphSuffix + ".png")
.string());
if (!history.residualForcesLog.empty()) {
createPlot("Number of Iterations",
"Residual Forces norm log",
history.residualForcesLog,
std::filesystem::path(graphsFolder)
.append("ResidualForcesLog_" + graphSuffix + ".png")
.string());
}
if (!history.potentialEnergies.empty()) {
createPlot("Number of Iterations", "Potential energy",
history.potentialEnergies,
std::filesystem::path(graphsFolder)
.append("PotentialEnergy" + graphSuffix + ".png")
.string());
createPlot("Number of Iterations",
"Potential energy",
history.potentialEnergies,
std::filesystem::path(graphsFolder)
.append("PotentialEnergy_" + graphSuffix + ".png")
.string());
}
if (!history.residualForcesMovingAverage.empty()) {
createPlot("Number of Iterations",
"Residual forces moving average",
history.residualForcesMovingAverage,
std::filesystem::path(graphsFolder)
.append("ResidualForcesMovingAverage_" + graphSuffix + ".png")
.string());
}
// if (!history.residualForcesMovingAverageDerivativesLog.empty()) {
// createPlot("Number of Iterations",
// "Residual forces moving average derivative log",
// history.residualForcesMovingAverageDerivativesLog,
// std::filesystem::path(graphsFolder)
// .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png")
// .string());
// }
if (!history.sumOfNormalizedDisplacementNorms.empty()) {
createPlot("Number of Iterations",
"Sum of normalized displacement norms",
history.sumOfNormalizedDisplacementNorms,
std::filesystem::path(graphsFolder)
.append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png")
.string());
}
}
};

5
simulationmesh.hpp
View File

@ -9,7 +9,7 @@
struct Element;
struct Node;
using CrossSectionType = RectangularBeamDimensions;
// using CrossSectionType = CylindricalBeamDimensions;
//using CrossSectionType = CylindricalBeamDimensions;
class SimulationMesh : public VCGEdgeMesh {
private:
@ -44,6 +44,9 @@ public:
double totalResidualForcesNorm{0};
double currentTotalPotentialEnergykN{0};
double previousTotalPotentialEnergykN{0};
double residualForcesMovingAverageDerivativeNorm{0};
double residualForcesMovingAverage{0};
double sumOfNormalizedDisplacementNorms{0};
bool save(const std::string &plyFilename = std::string());
void setBeamCrossSection(const CrossSectionType &beamDimensions);
void setBeamMaterial(const double &pr, const double &ym);

2
topologyenumerator.cpp
View File

@ -464,7 +464,7 @@ void TopologyEnumerator::computeValidPatterns(
const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold
= patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 5);
if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) {
if (debugIsOn || savePlyFiles) {
if (debugIsOn /*|| savePlyFiles*/) {
if (savePlyFiles) {
auto danglingEdgesPath = std::filesystem::path(resultsPath)
.append("ExceedingAngleThreshold");

4
trianglepatterngeometry.cpp
View File

@ -435,7 +435,7 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector<size_t> &nu
thetaAnglesOfIncidentVectors[thetaAngleIndex]
- thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]);
if (absAngleDifference < angleThreshold_rad
&& absAngleDifference > vcg::math::ToRad(0.01)) {
/*&& absAngleDifference > vcg::math::ToRad(0.01)*/) {
// std::cout << "Found angDiff:" << absAngleDifference << std::endl;
// vert[vi].C() = vcg::Color4b::Magenta;
// hasAngleSmallerThan = true;
@ -817,7 +817,7 @@ std::shared_ptr<PatternGeometry> PatternGeometry::tilePattern(
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices,
std::vector<int> &tileIntoEdgesToTiledVi,
std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexToTiledPatternEdgeIndex)
{
perPatternIndexToTiledPatternEdgeIndex.resize(patterns.size());

2
trianglepatterngeometry.hpp
View File

@ -121,7 +121,7 @@ private:
const std::vector<int> &connectToNeighborsVi,
const VCGPolyMesh &tileInto,
const std::vector<size_t> &perSurfaceFacePatternIndices,
std::vector<int> &tileIntoEdgesToTiledVi,
std::vector<size_t> &tileIntoEdgesToTiledVi,
std::vector<std::vector<size_t>> &perPatternIndexTiledPatternEdgeIndex);
};

34
utilities.hpp
View File

@ -130,13 +130,13 @@ inline void parseIntegers(const std::string &str, std::vector<size_t> &result) {
}
inline Eigen::MatrixXd toEigenMatrix(const std::vector<Vector6d> &v) {
Eigen::MatrixXd m(v.size(), 6);
Eigen::MatrixXd m(v.size(), 6);
for (size_t vi = 0; vi < v.size(); vi++) {
const Vector6d &vec = v[vi];
for (size_t i = 0; i < 6; i++) {
m(vi, i) = vec[i];
}
for (size_t vi = 0; vi < v.size(); vi++) {
const Vector6d &vec = v[vi];
for (size_t i = 0; i < 6; i++) {
m(vi, i) = vec[i];
}
}
return m;
@ -220,7 +220,7 @@ inline void deinitPolyscope()
return;
}
polyscope::shutdown();
polyscope::render::engine->shutdownImGui();
}
inline void init()
@ -302,4 +302,24 @@ std::string to_string_with_precision(const T a_value, const int n = 2)
return out.str();
}
template<typename T>
size_t computeHashUnordered(const std::vector<T> &v)
{
size_t hash = 0;
for (const auto &el : v) {
hash += std::hash<T>{}(el);
}
return hash;
}
inline size_t computeHashOrdered(const std::vector<size_t> &v)
{
std::string elementsString;
for (const auto &el : v) {
elementsString += std::to_string(el);
}
return std::hash<std::string>{}(elementsString);
}
#endif // UTILITIES_H