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

refactoring. Added functionality to export internal forces

This commit is contained in:
iasonmanolas 2022-05-06 16:36:47 +03:00
parent f531b16b19
commit 7f11ea8a47
1 changed files with 193 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

233
simulation_structs.hpp
View File

@ -1,6 +1,12 @@
#ifndef SIMULATIONSTRUCTS_HPP
#define SIMULATIONSTRUCTS_HPP
#include "csvfile.hpp"
#include "nlohmann/json.hpp"
#include "simulationmesh.hpp"
#include "utilities.hpp"
#include <fstream>
#include <string>
#include <vector>
namespace Eigen {
template <class Matrix>
@ -33,10 +39,6 @@ void read_binary(const std::string &filename, Matrix &matrix) {
//}
} // namespace Eigen
#include "simulationmesh.hpp"
#include "nlohmann/json.hpp"
#include <string>
#include <vector>
struct SimulationHistory {
SimulationHistory() {}
@ -49,8 +51,13 @@ struct SimulationHistory {
std::vector<size_t> redMarks;
std::vector<double> greenMarks;
std::vector<double> residualForcesMovingAverage;
std::vector<double> sumOfNormalizedDisplacementNorms;
// std::vector<double> residualForcesMovingAverageDerivativesLog;
std::vector<double> perVertexAverageNormalizedDisplacementNorm;
std::vector<double> residualForcesMovingAverageDerivativesLog;
//internal forces
std::vector<double> logOfSumOfAxialForcesNorm;
std::vector<double> logOfSumOfTorsionForcesNorm;
std::vector<double> logOfSumOfFirstBendingForcesNorm;
std::vector<double> logOfSumOfSecondBendingForcesNorm;
void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); }
@ -62,9 +69,47 @@ struct SimulationHistory {
// potentialEnergy.push_back(mesh.totalPotentialEnergykN);
logResidualForces.push_back(std::log10(mesh.totalResidualForcesNorm));
residualForcesMovingAverage.push_back(std::log10(mesh.residualForcesMovingAverage));
sumOfNormalizedDisplacementNorms.push_back(std::log10(mesh.sumOfNormalizedDisplacementNorms));
// residualForcesMovingAverageDerivativesLog.push_back(
// std::log(mesh.residualForcesMovingAverageDerivativeNorm));
perVertexAverageNormalizedDisplacementNorm.push_back(
mesh.perVertexAverageNormalizedDisplacementNorm);
residualForcesMovingAverageDerivativesLog.push_back(
std::log(mesh.residualForcesMovingAverageDerivativeNorm));
//Internal forces
const double axialSumOfNorms = std::accumulate(mesh.nodes._handle->data.begin(),
mesh.nodes._handle->data.end(),
0.0,
[](const double &sum, const Node &node) {
return sum
+ node.force.internalAxial.norm();
});
logOfSumOfAxialForcesNorm.push_back(std::log(axialSumOfNorms));
const double torsionSumOfNorms
= std::accumulate(mesh.nodes._handle->data.begin(),
mesh.nodes._handle->data.end(),
0.0,
[](const double &sum, const Node &node) {
return sum + node.force.internalTorsion.norm();
});
logOfSumOfTorsionForcesNorm.push_back(std::log(torsionSumOfNorms));
const double firstBendingSumOfNorms
= std::accumulate(mesh.nodes._handle->data.begin(),
mesh.nodes._handle->data.end(),
0.0,
[](const double &sum, const Node &node) {
return sum + node.force.internalFirstBending.norm();
});
logOfSumOfFirstBendingForcesNorm.push_back(std::log(firstBendingSumOfNorms));
const double secondBendingSumOfNorms
= std::accumulate(mesh.nodes._handle->data.begin(),
mesh.nodes._handle->data.end(),
0.0,
[](const double &sum, const Node &node) {
return sum + node.force.internalSecondBending.norm();
});
logOfSumOfSecondBendingForcesNorm.push_back(std::log(secondBendingSumOfNorms));
}
void clear()
@ -73,7 +118,7 @@ struct SimulationHistory {
kineticEnergy.clear();
potentialEnergies.clear();
residualForcesMovingAverage.clear();
sumOfNormalizedDisplacementNorms.clear();
perVertexAverageNormalizedDisplacementNorm.clear();
// residualForcesMovingAverageDerivativesLog.clear();
}
};
@ -190,6 +235,10 @@ public:
return json.dump();
}
bool operator==(const SimulationJob &otherSimulationJob)
{
return this->toString() == otherSimulationJob.toString();
}
void clear()
{
@ -416,11 +465,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
});
// if (!nodeColors.empty()) {
// polyscope::getCurveNetwork(meshLabel)
// ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors)
// ->setEnabled(shouldEnable);
// }
if (!nodeColors.empty()) {
polyscope::getCurveNetwork(meshLabel)
->addNodeColorQuantity("Boundary conditions", nodeColors)
->setEnabled(shouldEnable);
}
// per node external forces
std::vector<std::array<double, 3>> externalForces(pMesh->VN());
@ -431,10 +480,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
if (!externalForces.empty()) {
const std::string polyscopeLabel_externalForces = "External force";
polyscope::getCurveNetwork(meshLabel)->removeQuantity(polyscopeLabel_externalForces);
polyscope::CurveNetworkNodeVectorQuantity *externalForcesVectors
= polyscope::getCurveNetwork(meshLabel)->addNodeVectorQuantity("External force_"
+ label,
externalForces);
= polyscope::getCurveNetwork(meshLabel)
->addNodeVectorQuantity(polyscopeLabel_externalForces, externalForces);
const std::array<double, 3> color_loads{1.0, 0, 0};
externalForcesVectors->setVectorColor(
@ -455,7 +505,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
if (hasExternalMoments) {
polyscope::getCurveNetwork(meshLabel)
->addNodeVectorQuantity("External moment_" + label, externalMoments)
->addNodeVectorQuantity("External moment", externalMoments)
->setEnabled(shouldEnable);
}
}
@ -465,10 +515,10 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
return;
}
if (!nodalExternalForces.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label);
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force");
}
if (!constrainedVertices.empty()) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label);
polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions");
}
// per node external moments
@ -481,7 +531,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m
}
}
if (hasExternalMoments) {
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment_" + label);
polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment");
}
}
#endif // POLYSCOPE_DEFINED
@ -503,6 +553,7 @@ struct SimulationResults
std::vector<Eigen::Quaternion<double>> rotationalDisplacementQuaternion; //per vertex
double internalPotentialEnergy{0};
double executionTime{0};
std::vector<std::array<Vector6d, 4>> perVertexInternalForces; //axial,torsion,bending1,bending2
std::string labelPrefix{"deformed"};
inline static char deliminator{' '};
SimulationResults() { pJob = std::make_shared<SimulationJob>(); }
@ -538,6 +589,53 @@ struct SimulationResults
return m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string());
}
void saveInternalForces(const std::filesystem::path &outputDirPath)
{
std::cout << "out to:" << outputDirPath << std::endl;
const std::filesystem::path internalForcesDirPath = std::filesystem::path(outputDirPath);
std::filesystem::create_directories(internalForcesDirPath);
csvFile csv_axial6d(std::filesystem::path(internalForcesDirPath)
.append("forces_axial_6d.csv"),
true);
csvFile csv_axialMagn(std::filesystem::path(internalForcesDirPath)
.append("forces_axial_magn.csv"),
true);
csvFile csv_torsion6d(std::filesystem::path(internalForcesDirPath)
.append("forces_torsion_6d.csv"),
true);
csvFile csv_torsionMagn(std::filesystem::path(internalForcesDirPath)
.append("forces_torsion_magn.csv"),
true);
csvFile csv_firstBending6d(std::filesystem::path(internalForcesDirPath)
.append("forces_firstBending_6d.csv"),
true);
csvFile csv_firstBendingMagn(std::filesystem::path(internalForcesDirPath)
.append("forces_firstBending_magn.csv"),
true);
csvFile csv_secondBending6d(std::filesystem::path(internalForcesDirPath)
.append("forces_secondBending_6d.csv"),
true);
csvFile csv_secondBendingMagn(std::filesystem::path(internalForcesDirPath)
.append("forces_secondBending_magn.csv"),
true);
for (const std::array<Vector6d, 4> &internalForce : perVertexInternalForces) {
for (int dofi = 0; dofi < 6; dofi++) {
csv_axial6d << internalForce[0][dofi];
csv_torsion6d << internalForce[1][dofi];
csv_firstBending6d << internalForce[2][dofi];
csv_secondBending6d << internalForce[3][dofi];
}
csv_axial6d << endrow;
csv_torsion6d << endrow;
csv_firstBending6d << endrow;
csv_secondBending6d << endrow;
csv_axialMagn << internalForce[0].norm() << endrow;
csv_torsionMagn << internalForce[1].norm() << endrow;
csv_firstBendingMagn << internalForce[2].norm() << endrow;
csv_secondBendingMagn << internalForce[3].norm() << endrow;
}
}
void save(const std::string &outputFolder = std::string())
{
const std::filesystem::path outputFolderPath = outputFolder.empty()
@ -559,7 +657,27 @@ struct SimulationResults
nlohmann::json json;
json[GET_VARIABLE_NAME(internalPotentialEnergy)] = internalPotentialEnergy;
//Write internal forces
if (!perVertexInternalForces.empty()) {
std::vector<Vector6d> internalForces_axial(perVertexInternalForces.size());
std::vector<Vector6d> internalForces_torsion(perVertexInternalForces.size());
std::vector<Vector6d> internalForces_firstBending(perVertexInternalForces.size());
std::vector<Vector6d> internalForces_secondBending(perVertexInternalForces.size());
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
internalForces_axial[vi] = perVertexInternalForces[vi][0];
internalForces_torsion[vi] = perVertexInternalForces[vi][1];
internalForces_firstBending[vi] = perVertexInternalForces[vi][2];
internalForces_secondBending[vi] = perVertexInternalForces[vi][3];
}
json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial"]
= internalForces_axial;
json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_torsion"]
= internalForces_torsion;
json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_firstBending"]
= internalForces_firstBending;
json[std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_secondBending"]
= internalForces_secondBending;
}
std::filesystem::path jsonFilePath(
std::filesystem::path(resultsFolderPath).append(defaultJsonFilename));
std::ofstream jsonFile(jsonFilePath.string());
@ -649,7 +767,7 @@ struct SimulationResults
const glm::vec3 desiredColor_glm(desiredColor.value()[0],
desiredColor.value()[1],
desiredColor.value()[2]);
polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm);
polyscopeHandle_deformedEdmeMesh->setColor(/*glm::normalize(*/ desiredColor_glm /*)*/);
}
Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3);
Eigen::MatrixX3d framesX(mesh->VN(), 3);
@ -723,10 +841,10 @@ struct SimulationResults
polyscopeHandle_frameZ->setVectorColor(
/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel());
if (!polyscopeHandle_initialMesh) {
polyscopeHandle_initialMesh = mesh->registerForDrawing();
}
// if (!polyscope::hasCurveNetwork(mesh->getLabel())) {
// const std::array<double, 3> initialColor({0, 0, 0});
// /*auto polyscopeHandle_initialMesh =*/mesh->registerForDrawing(initialColor);
// }
// auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh
// ->addNodeVectorQuantity("FrameX", framesX_initial);
@ -747,7 +865,7 @@ struct SimulationResults
// polyscopeHandle_frameZ_initial->setVectorColor(
// /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1));
// //}
pJob->registerForDrawing(getLabel());
pJob->registerForDrawing(getLabel(), false);
// static bool wasExecuted =false;
// if (!wasExecuted) {
// std::function<void()> callback = [&]() {
@ -794,17 +912,52 @@ private:
const std::filesystem::path jsonFilepath = std::filesystem::path(loadFromPath)
.append(defaultJsonFilename);
if (!std::filesystem::exists(jsonFilepath)) {
std::cerr << "Simulation results could not be loaded because filepath does "
"not exist:"
<< jsonFilepath << std::endl;
return false;
}
std::ifstream ifs(jsonFilepath);
nlohmann::json json;
ifs >> json;
if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
if (std::filesystem::exists(jsonFilepath)) {
std::ifstream ifs(jsonFilepath);
nlohmann::json json;
ifs >> json;
// if (json.contains(GET_VARIABLE_NAME(internalPotentialEnergy))) {
// internalPotentialEnergy = json.at(GET_VARIABLE_NAME(internalPotentialEnergy));
// }
if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial")) {
perVertexInternalForces.resize(pJob->pMesh->VN());
std::vector<Vector6d> perVertexInternalForces_axial
= static_cast<std::vector<Vector6d>>(json.at(
std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_axial"));
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
}
}
if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+ "_torsion")) {
perVertexInternalForces.resize(pJob->pMesh->VN());
std::vector<Vector6d> perVertexInternalForces_axial
= static_cast<std::vector<Vector6d>>(json.at(
std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_torsion"));
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
}
}
if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+ "_firstBending")) {
perVertexInternalForces.resize(pJob->pMesh->VN());
std::vector<Vector6d> perVertexInternalForces_axial
= static_cast<std::vector<Vector6d>>(json.at(
std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_firstBending"));
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
}
}
if (json.contains(std::string(GET_VARIABLE_NAME(perVertexInternalForces))
+ "_secondBending")) {
perVertexInternalForces.resize(pJob->pMesh->VN());
std::vector<Vector6d> perVertexInternalForces_axial
= static_cast<std::vector<Vector6d>>(json.at(
std::string(GET_VARIABLE_NAME(perVertexInternalForces)) + "_secondBending"));
for (int vi = 0; vi < pJob->pMesh->VN(); vi++) {
perVertexInternalForces[vi][0] = perVertexInternalForces_axial[vi];
}
}
}
return true;
}