diff --git a/simulation_structs.hpp b/simulation_structs.hpp index f1e267e..aa6cbfb 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -1,6 +1,12 @@ #ifndef SIMULATIONSTRUCTS_HPP #define SIMULATIONSTRUCTS_HPP +#include "csvfile.hpp" +#include "nlohmann/json.hpp" +#include "simulationmesh.hpp" +#include "utilities.hpp" #include +#include +#include namespace Eigen { template @@ -33,10 +39,6 @@ void read_binary(const std::string &filename, Matrix &matrix) { //} } // namespace Eigen -#include "simulationmesh.hpp" -#include "nlohmann/json.hpp" -#include -#include struct SimulationHistory { SimulationHistory() {} @@ -49,8 +51,13 @@ struct SimulationHistory { std::vector redMarks; std::vector greenMarks; std::vector residualForcesMovingAverage; - std::vector sumOfNormalizedDisplacementNorms; - // std::vector residualForcesMovingAverageDerivativesLog; + std::vector perVertexAverageNormalizedDisplacementNorm; + std::vector residualForcesMovingAverageDerivativesLog; + //internal forces + std::vector logOfSumOfAxialForcesNorm; + std::vector logOfSumOfTorsionForcesNorm; + std::vector logOfSumOfFirstBendingForcesNorm; + std::vector 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> 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 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> rotationalDisplacementQuaternion; //per vertex double internalPotentialEnergy{0}; double executionTime{0}; + std::vector> perVertexInternalForces; //axial,torsion,bending1,bending2 std::string labelPrefix{"deformed"}; inline static char deliminator{' '}; SimulationResults() { pJob = std::make_shared(); } @@ -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 &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 internalForces_axial(perVertexInternalForces.size()); + std::vector internalForces_torsion(perVertexInternalForces.size()); + std::vector internalForces_firstBending(perVertexInternalForces.size()); + std::vector 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 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 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 perVertexInternalForces_axial + = static_cast>(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 perVertexInternalForces_axial + = static_cast>(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 perVertexInternalForces_axial + = static_cast>(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 perVertexInternalForces_axial + = static_cast>(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; }