From 5e1fae8b245f839c6067f98caa14f920978ff095 Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 1 Mar 2021 15:44:35 +0200 Subject: [PATCH] Refactoring --- beamformfinder.cpp | 233 ++------------------------- beamformfinder.hpp | 4 +- csvfile.hpp | 50 +++--- edgemesh.cpp | 41 +++-- edgemesh.hpp | 8 +- elementalmesh.cpp | 32 +--- optimizationstructs.hpp | 5 + simulationresult.hpp | 338 +++++++++++++++++++++++----------------- utilities.hpp | 15 +- 9 files changed, 296 insertions(+), 430 deletions(-) create mode 100644 optimizationstructs.hpp diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 87e2a2a..018d6f7 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -1,7 +1,4 @@ #include "beamformfinder.hpp" -#include "polyscope/curve_network.h" -#include "polyscope/histogram.h" -#include "polyscope/polyscope.h" #include "simulationhistoryplotter.hpp" #include #include @@ -18,8 +15,9 @@ void FormFinder::runUnitTests() { VCGEdgeMesh beam; // const size_t spanGridSize = 11; // mesh.createSpanGrid(spanGridSize); - beam.loadPly( - std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string()); + beam.loadPly(std::filesystem::path(groundOfTruthFolder) + .append("simpleBeam.ply") + .string()); std::unordered_map> fixedVertices; fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; fixedVertices[beam.VN() - 1] = std::unordered_set{1, 2}; @@ -50,7 +48,8 @@ void FormFinder::runUnitTests() { simpleBeam_simulationResults.save(); const std::string simpleBeamGroundOfTruthBinaryFilename = std::filesystem::path(groundOfTruthFolder) - .append("SimpleBeam_displacements.eigenBin").string(); + .append("SimpleBeam_displacements.eigenBin") + .string(); assert(std::filesystem::exists( std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename))); Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements; @@ -219,7 +218,7 @@ void FormFinder::reset() { externalMomentsNorm = 0; mSettings.drawingStep = 1; Dt = mSettings.Dtini; - numOfDampings=0; + numOfDampings = 0; } VectorType FormFinder::computeDisplacementDifferenceDerivative( @@ -794,7 +793,7 @@ void FormFinder::updateResidualForcesOnTheFly( std::vector>(4, {-1, Vector6d()})); // omp_lock_t writelock; // omp_init_lock(&writelock); -//#pragma omp parallel for //schedule(static) num_threads(8) + //#pragma omp parallel for //schedule(static) num_threads(8) for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; const SimulationMesh::VertexType &ev_j = *e.cV(0); @@ -1568,183 +1567,6 @@ FormFinder::computeResults(const std::shared_ptr &pJob) { return SimulationResults{true, pJob, history, displacements}; } -void FormFinder::draw(const std::string &screenshotsFolder = {}) { - // update positions - // polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->updateNodePositions(pMesh->getEigenVertices()); - - // Vertex quantities - std::vector kineticEnergies(pMesh->VN()); - std::vector> nodeNormals(pMesh->VN()); - std::vector> internalForces(pMesh->VN()); - std::vector> externalForces(pMesh->VN()); - std::vector> externalMoments(pMesh->VN()); - std::vector internalForcesNorm(pMesh->VN()); - std::vector> internalAxialForces(pMesh->VN()); - std::vector> internalFirstBendingTranslationForces( - pMesh->VN()); - std::vector> internalFirstBendingRotationForces( - pMesh->VN()); - std::vector> internalSecondBendingTranslationForces( - pMesh->VN()); - std::vector> internalSecondBendingRotationForces( - pMesh->VN()); - std::vector nRs(pMesh->VN()); - std::vector theta2(pMesh->VN()); - std::vector theta3(pMesh->VN()); - std::vector> residualForces(pMesh->VN()); - std::vector residualForcesNorm(pMesh->VN()); - std::vector accelerationX(pMesh->VN()); - for (const VertexType &v : pMesh->vert) { - kineticEnergies[pMesh->getIndex(v)] = pMesh->nodes[v].kineticEnergy; - const VectorType n = v.cN(); - nodeNormals[pMesh->getIndex(v)] = {n[0], n[1], n[2]}; - // per node internal forces - const Vector6d nodeforce = pMesh->nodes[v].force.internal * (-1); - internalForces[pMesh->getIndex(v)] = {nodeforce[0], nodeforce[1], - nodeforce[2]}; - internalForcesNorm[pMesh->getIndex(v)] = nodeforce.norm(); - // External force - const Vector6d nodeExternalForce = pMesh->nodes[v].force.external; - externalForces[pMesh->getIndex(v)] = { - nodeExternalForce[0], nodeExternalForce[1], nodeExternalForce[2]}; - externalMoments[pMesh->getIndex(v)] = {nodeExternalForce[3], - nodeExternalForce[4], 0}; - internalAxialForces[pMesh->getIndex(v)] = {nodeforce[0], nodeforce[1], - nodeforce[2]}; - const Node &node = pMesh->nodes[v]; - const Vector6d nodeForceFirst = node.force.internalFirstBending * (-1); - internalFirstBendingTranslationForces[pMesh->getIndex(v)] = { - nodeForceFirst[0], nodeForceFirst[1], nodeForceFirst[2]}; - internalFirstBendingRotationForces[pMesh->getIndex(v)] = { - nodeForceFirst[3], nodeForceFirst[4], 0}; - - const Vector6d nodeForceSecond = node.force.internalSecondBending * (-1); - internalSecondBendingTranslationForces[pMesh->getIndex(v)] = { - nodeForceSecond[0], nodeForceSecond[1], nodeForceSecond[2]}; - internalSecondBendingRotationForces[pMesh->getIndex(v)] = { - nodeForceSecond[3], nodeForceSecond[4], 0}; - nRs[pMesh->getIndex(v)] = node.nR; - const Vector6d nodeResidualForce = pMesh->nodes[v].force.residual; - residualForces[pMesh->getIndex(v)] = { - nodeResidualForce[0], nodeResidualForce[1], nodeResidualForce[2]}; - residualForcesNorm[pMesh->getIndex(v)] = nodeResidualForce.norm(); - accelerationX[pMesh->getIndex(v)] = pMesh->nodes[v].acceleration[0]; - } - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeScalarQuantity("Kinetic Energy", kineticEnergies) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("Node normals", nodeNormals) - ->setEnabled(true); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("Internal force", internalForces) - ->setEnabled(false); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("External force", externalForces) - ->setEnabled(true); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addNodeVectorQuantity("External Moments", externalMoments) - ->setEnabled(true); - 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("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); - - // Edge quantities - std::vector A(pMesh->EN()); - std::vector J(pMesh->EN()); - std::vector I2(pMesh->EN()); - std::vector I3(pMesh->EN()); - for (const EdgeType &e : pMesh->edge) { - const size_t ei = pMesh->getIndex(e); - A[ei] = pMesh->elements[e].A; - J[ei] = pMesh->elements[e].J; - I2[ei] = pMesh->elements[e].I2; - I3[ei] = pMesh->elements[e].I3; - } - - polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); - polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addEdgeScalarQuantity("I2", I2); - polyscope::getCurveNetwork(meshPolyscopeLabel) - ->addEdgeScalarQuantity("I3", I3); - - // Specify the callback - polyscope::state::userCallback = [&]() { - // 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::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(); - }; - - if (!screenshotsFolder.empty()) { - static bool firstDraw = true; - if (firstDraw) { - for (const auto &entry : - std::filesystem::directory_iterator(screenshotsFolder)) - std::filesystem::remove_all(entry.path()); - // polyscope::view::processZoom(5); - firstDraw = false; - } - polyscope::screenshot( - std::filesystem::path(screenshotsFolder) - .append(std::to_string(mCurrentSimulationStep) + ".png") - .string(), - false); - } else { - polyscope::show(); - } -} - void FormFinder::printCurrentState() const { std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl; @@ -1793,18 +1615,11 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, pMesh.reset(); pMesh = std::make_unique(*pJob->pMesh); - if (mSettings.beVerbose ) { + if (mSettings.beVerbose) { std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " << pMesh->EN() << " elements." << std::endl; } computeRigidSupports(); - - if (mSettings.shouldDraw ) { - initPolyscope(); - polyscope::registerCurveNetwork( - meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges()); - // registerWorldAxes(); - } for (auto fixedVertex : pJob->constrainedVertices) { assert(fixedVertex.first < pMesh->VN()); } @@ -1853,25 +1668,9 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, } if (std::isnan(pMesh->currentTotalKineticEnergy)) { - std::time_t now = time(0); - std::tm *ltm = std::localtime(&now); - 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 << "Infinite kinetic energy detected.Exiting.." << std::endl; - FormFinder debug; - FormFinder::Settings settings; - settings.shouldDraw = true; - settings.beVerbose = true; - debug.executeSimulation(pJob, settings); - std::terminate(); + if (mSettings.beVerbose) { + std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; + } break; } @@ -1884,8 +1683,6 @@ currentSimulationStep > maxDRMIterations*/ // .append("Screenshots") // .string(); // draw(saveTo); - draw(); - // yValues = history.kineticEnergy; // plotHandle = matplot::scatter(xPlot, yValues); // label = "Log of Kinetic energy"; @@ -1971,7 +1768,7 @@ mesh->currentTotalPotentialEnergykN*/ } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { results.converged = false; - } else if (mSettings.beVerbose ) { + } else if (mSettings.beVerbose) { auto t2 = std::chrono::high_resolution_clock::now(); double simulationDuration = std::chrono::duration_cast(t2 - t1).count(); @@ -1990,12 +1787,6 @@ mesh->currentTotalPotentialEnergykN*/ SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); } - if (mSettings.shouldDraw) { - draw(); - } - } - if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) { - polyscope::removeCurveNetwork(meshPolyscopeLabel); } return results; } diff --git a/beamformfinder.hpp b/beamformfinder.hpp index f1b0b7c..79797fc 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -3,8 +3,8 @@ #include "elementalmesh.hpp" #include "matplot/matplot.h" -#include "polyscope/curve_network.h" -#include "polyscope/polyscope.h" +//#include "polyscope/curve_network.h" +//#include "polyscope/polyscope.h" #include "simulationresult.hpp" #include #include diff --git a/csvfile.hpp b/csvfile.hpp index 9bd0511..303829f 100644 --- a/csvfile.hpp +++ b/csvfile.hpp @@ -5,12 +5,12 @@ #include #include -class csvfile; +class csvFile; -inline static csvfile &endrow(csvfile &file); -inline static csvfile &flush(csvfile &file); +inline static csvFile &endrow(csvFile &file); +inline static csvFile &flush(csvFile &file); -class csvfile { +class csvFile { std::ofstream fs_; bool is_first_; const std::string separator_; @@ -18,24 +18,30 @@ class csvfile { const std::string special_chars_; public: - csvfile(const std::string filename, const bool &overwrite, - const std::string separator = ";") + csvFile(const std::string filename, const bool &overwrite, + const std::string separator = ",") : fs_(), is_first_(true), separator_(separator), escape_seq_("\""), special_chars_("\"") { fs_.exceptions(std::ios::failbit | std::ios::badbit); - if (!std::filesystem::exists(std::filesystem::path("../OptimizationResults") - .append("statistics.csv") - )) { + if (filename.empty()) { + fs_.copyfmt(std::cout); + fs_.clear(std::cout.rdstate()); + fs_.basic_ios::rdbuf(std::cout.rdbuf()); + } else { + if (!std::filesystem::exists( + std::filesystem::path("../OptimizationResults") + .append("statistics.csv"))) { std::ofstream outfile(std::filesystem::path("../OptimizationResults") - .append("statistics.csv") - .string()); + .append("statistics.csv") + .string()); outfile.close(); -} - overwrite ? fs_.open(filename, std::ios::trunc) - : fs_.open(filename, std::ios::app); + } + overwrite ? fs_.open(filename, std::ios::trunc) + : fs_.open(filename, std::ios::app); + } } - ~csvfile() { + ~csvFile() { flush(); fs_.close(); } @@ -47,16 +53,16 @@ public: is_first_ = true; } - csvfile &operator<<(csvfile &(*val)(csvfile &)) { return val(*this); } + csvFile &operator<<(csvFile &(*val)(csvFile &)) { return val(*this); } - csvfile &operator<<(const char *val) { return write(escape(val)); } + csvFile &operator<<(const char *val) { return write(escape(val)); } - csvfile &operator<<(const std::string &val) { return write(escape(val)); } + csvFile &operator<<(const std::string &val) { return write(escape(val)); } - template csvfile &operator<<(const T &val) { return write(val); } + template csvFile &operator<<(const T &val) { return write(val); } private: - template csvfile &write(const T &val) { + template csvFile &write(const T &val) { if (!is_first_) { fs_ << separator_; } else { @@ -80,12 +86,12 @@ private: } }; -inline static csvfile &endrow(csvfile &file) { +inline static csvFile &endrow(csvFile &file) { file.endrow(); return file; } -inline static csvfile &flush(csvfile &file) { +inline static csvFile &flush(csvFile &file) { file.flush(); return file; } diff --git a/edgemesh.cpp b/edgemesh.cpp index a32257a..976c2b3 100644 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -183,18 +183,21 @@ bool VCGEdgeMesh::loadPly(const std::string plyFilename) { assert(std::filesystem::exists(usedPath)); this->Clear(); const bool useDefaultImporter = false; - if (!loadUsingNanoply(usedPath)) { - std::cerr << "Error: Unable to open " + usedPath << std::endl; - return false; - } + if (!loadUsingNanoply(usedPath)) { + std::cerr << "Error: Unable to open " + usedPath << std::endl; + return false; + } getEdges(eigenEdges); getVertices(eigenVertices); vcg::tri::UpdateTopology::VertexEdge(*this); label = std::filesystem::path(plyFilename).stem().string(); - std::cout << plyFilename << " was loaded successfuly." << std::endl; - std::cout << "Mesh has " << EN() << " edges." << std::endl; + const bool printDebugInfo = false; + if (printDebugInfo) { + std::cout << plyFilename << " was loaded successfuly." << std::endl; + std::cout << "Mesh has " << EN() << " edges." << std::endl; + } return true; } @@ -337,9 +340,27 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const { << " " << vert[vi].cP()[2] << std::endl; } -void VCGEdgeMesh::registerForDrawing() const { +#ifdef POLYSCOPE_DEFINED +void VCGEdgeMesh::registerForDrawing( + const std::optional &desiredColor) const { initPolyscope(); - const double drawingRadius = 0.0007; - polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges()) - ->setRadius(drawingRadius, false); + const double drawingRadius = 0.002; + auto polyscopeHandle_edgeMesh = + polyscope::registerCurveNetwork(label, getEigenVertices(), + getEigenEdges()) + ->setRadius(drawingRadius, true); + if (desiredColor.has_value()) { + polyscopeHandle_edgeMesh->setColor(desiredColor.value()); + } } + +void VCGEdgeMesh::unregister() const { + if (!polyscope::hasCurveNetwork(label)) { + std::cerr << "No curve network registered with a name: " << getLabel() + << std::endl; + std::cerr << "Nothing to remove." << std::endl; + return; + } + polyscope::removeCurveNetwork(label); +} +#endif diff --git a/edgemesh.hpp b/edgemesh.hpp index e420a4c..7995507 100644 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -1,6 +1,7 @@ #ifndef EDGEMESH_HPP #define EDGEMESH_HPP #include "beam.hpp" +#include "externvariables.hpp" #include "polymesh.hpp" #include "utilities.hpp" #include "vcgtrimesh.hpp" @@ -77,8 +78,11 @@ public: Eigen::MatrixX3d getEigenVertices() const; Eigen::MatrixX3d getEigenEdgeNormals() const; void printVertexCoordinates(const size_t &vi) const; - void registerForDrawing() const; - +#ifdef POLYSCOPE_DEFINED + void registerForDrawing( + const std::optional &desiredColor = std::nullopt) const; + void unregister() const; +#endif std::string getLabel() const; void setLabel(const std::string &value); diff --git a/elementalmesh.cpp b/elementalmesh.cpp index 1647c19..0cc48b2 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -38,12 +38,13 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) { eigenVertices = mesh.getEigenVertices(); } -SimulationMesh::~SimulationMesh() -{ - vcg::tri::Allocator::DeletePerEdgeAttribute(*this, elements); - vcg::tri::Allocator::DeletePerVertexAttribute(*this,nodes); -} - +SimulationMesh::~SimulationMesh() { + vcg::tri::Allocator::DeletePerEdgeAttribute(*this, + elements); + vcg::tri::Allocator::DeletePerVertexAttribute(*this, + nodes); +} + SimulationMesh::SimulationMesh(FlatPattern &pattern) { vcg::tri::MeshAssert::VertexNormalNormalized(pattern); @@ -305,6 +306,7 @@ bool SimulationMesh::loadPly(const string &plyFilename) { mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; + mask |= nanoply::NanoPlyWrapper::IO_MESHATTRIB; if (nanoply::NanoPlyWrapper::LoadModel( plyFilename.c_str(), *this, mask, customAttrib) != 0) { return false; @@ -358,19 +360,6 @@ bool SimulationMesh::savePly(const std::string &plyFilename) { customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, material.data()); - // std::vector> beamProperties(EN()); - // for (size_t ei = 0; ei < EN(); ei++) { - // auto props = elements[ei].toArray(); - // for (auto p : props) { - // std::cout << p << " "; - // } - // std::cout << std::endl; - // beamProperties[ei] = props; - // } - // customAttrib.AddEdgeAttribDescriptor, double, 6>( - // plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, - // beamProperties.data()); - // Load the ply file unsigned int mask = 0; mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; @@ -386,11 +375,6 @@ bool SimulationMesh::savePly(const std::string &plyFilename) { SimulationMesh::EdgePointer SimulationMesh::getReferenceElement(const VCGEdgeMesh::VertexType &v) { const VertexIndex vi = getIndex(v); - // return nodes[v].incidentElements[0]; - // if (vi == 0 || vi == 1) { - // return nodes[0].incidentElements[0]; - // } - return nodes[v].incidentElements[0]; } diff --git a/optimizationstructs.hpp b/optimizationstructs.hpp new file mode 100644 index 0000000..80077d0 --- /dev/null +++ b/optimizationstructs.hpp @@ -0,0 +1,5 @@ +#ifndef OPTIMIZATIONSTRUCTS_HPP +#define OPTIMIZATIONSTRUCTS_HPP +#include "beamformfinder.hpp" + +#endif // OPTIMIZATIONSTRUCTS_HPP diff --git a/simulationresult.hpp b/simulationresult.hpp index 1e60842..f1a7533 100644 --- a/simulationresult.hpp +++ b/simulationresult.hpp @@ -60,8 +60,8 @@ class SimulationJob { inline static std::string nodalForces{"forces"}; inline static std::string label{"label"}; inline static std::string meshLabel{ - "label"}; // TODO: should be in the savePly function of the simulation - // mesh class + "meshLabel"}; // TODO: should be in the savePly function of the + // simulation mesh class } jsonLabels; public: @@ -79,24 +79,155 @@ public: nodalExternalForces(ef), nodalForcedDisplacements(fd) {} SimulationJob() {} SimulationJob(const std::string &jsonFilename) { load(jsonFilename); } + + SimulationJob getCopy() const { + SimulationJob jobCopy; + jobCopy.pMesh = std::make_shared(); + jobCopy.pMesh->copy(*pMesh); + jobCopy.label = label; + jobCopy.constrainedVertices = constrainedVertices; + jobCopy.nodalExternalForces = nodalExternalForces; + jobCopy.nodalForcedDisplacements = nodalForcedDisplacements; + + return jobCopy; + } std::string getLabel() const { return label; } std::string toString() const { nlohmann::json json; if (!constrainedVertices.empty()) { - json[jsonLabel_constrainedVertices] = constrainedVertices; + json[jsonLabels.constrainedVertices] = constrainedVertices; } if (!nodalExternalForces.empty()) { std::unordered_map> arrForces; for (const auto &f : nodalExternalForces) { arrForces[f.first] = f.second; } - json[jsonLabel_nodalForces] = arrForces; + json[jsonLabels.nodalForces] = arrForces; } return json.dump(); } + bool load(const std::string &jsonFilename) { + if (std::filesystem::path(jsonFilename).extension() != ".json") { + std::cerr << "A json file is expected as input. The given file has the " + "following extension:" + << std::filesystem::path(jsonFilename).extension() << std::endl; + assert(false); + return false; + } + if (!std::filesystem::exists(std::filesystem::path(jsonFilename))) { + std::cerr << "The json file does not exist. Json file provided:" + << jsonFilename << std::endl; + assert(false); + return false; + } + + std::cout << "Loading json file:" << jsonFilename << std::endl; + nlohmann::json json; + std::ifstream ifs(jsonFilename); + ifs >> json; + + pMesh = std::make_shared(); + if (json.contains(jsonLabels.meshFilename)) { + const std::string relativeFilepath = json[jsonLabels.meshFilename]; + const auto meshFilepath = + std::filesystem::path( + std::filesystem::path(jsonFilename).parent_path()) + .append(relativeFilepath); + pMesh->loadPly(meshFilepath.string()); + pMesh->setLabel( + json[jsonLabels.meshLabel]); // FIXME: This should be exported using + // nanoply but nanoply might not be able + // to write a string(??) + } + + if (json.contains(jsonLabels.constrainedVertices)) { + constrainedVertices = + // auto conV = + json[jsonLabels.constrainedVertices] + .get>>(); + std::cout << "Loaded constrained vertices. Number of constrained " + "vertices found:" + << constrainedVertices.size() << std::endl; + } + + if (json.contains(jsonLabels.nodalForces)) { + auto f = std::unordered_map>( + json[jsonLabels.nodalForces]); + for (const auto &forces : f) { + nodalExternalForces[forces.first] = Vector6d(forces.second); + } + std::cout << "Loaded forces. Number of forces found:" + << nodalExternalForces.size() << std::endl; + } + + if (json.contains(jsonLabels.label)) { + label = json[jsonLabels.label]; + } + + if (json.contains(jsonLabels.meshLabel)) { + pMesh->setLabel(json[jsonLabels.meshLabel]); + } + + return true; + } + + bool save(const std::string &folderDirectory) const { + const std::filesystem::path pathFolderDirectory(folderDirectory); + if (!std::filesystem::is_directory(pathFolderDirectory)) { + std::cerr << "A folder directory is expected for saving the simulation " + "job. Exiting.." + << std::endl; + return false; + } + + bool returnValue = true; + std::string jsonFilename( + std::filesystem::path(pathFolderDirectory) + .append(label + "_" + pMesh->getLabel() + "_simulationJob.json") + .string()); + + const std::string meshFilename = + std::filesystem::absolute( + std::filesystem::canonical( + std::filesystem::path(pathFolderDirectory))) + .append(pMesh->getLabel() + ".ply") + .string(); + returnValue = pMesh->savePly(meshFilename); + nlohmann::json json; + json[jsonLabels.meshFilename] = std::filesystem::relative( + std::filesystem::path(meshFilename), + std::filesystem::path( + std::filesystem::path(jsonFilename).parent_path())); + json[jsonLabels.meshLabel] = + pMesh->getLabel(); // FIXME: This should be exported using nanoply but + // nanoply might not be able to write a string(??) + if (!constrainedVertices.empty()) { + json[jsonLabels.constrainedVertices] = constrainedVertices; + } + if (!nodalExternalForces.empty()) { + std::unordered_map> arrForces; + for (const auto &f : nodalExternalForces) { + arrForces[f.first] = f.second; + } + json[jsonLabels.nodalForces] = arrForces; + } + + if (!label.empty()) { + json[jsonLabels.label] = label; + } + if (!pMesh->getLabel().empty()) { + json[jsonLabels.meshLabel] = pMesh->getLabel(); + } + std::ofstream jsonFile(jsonFilename); + jsonFile << json; + // std::cout << "Saved simulation job as:" << jsonFilename << std::endl; + + return returnValue; + } +#ifdef POLYSCOPE_DEFINED void registerForDrawing(const std::string &meshLabel) const { initPolyscope(); if (meshLabel.empty()) { @@ -163,108 +294,7 @@ public: ->setEnabled(false); } } - - bool load(const std::string &jsonFilename) { - if (std::filesystem::path(jsonFilename).extension() != ".json") { - std::cerr << "A json file is expected as input. The given file has the " - "following extension:" - << std::filesystem::path(jsonFilename).extension() << std::endl; - assert(false); - return false; - } - - if (!std::filesystem::exists(std::filesystem::path(jsonFilename))) { - std::cerr << "The json file does not exist. Json file provided:" - << jsonFilename << std::endl; - assert(false); - return false; - } - - std::cout << "Loading json file:" << jsonFilename << std::endl; - nlohmann::json json; - std::ifstream ifs(jsonFilename); - json << ifs; - - pMesh = std::make_shared(); - if (json.contains(jsonLabels.meshFilename)) { - pMesh->loadPly(json[jsonLabels.meshFilename]); - } - - if (json.contains(jsonLabels.constrainedVertices)) { - constrainedVertices = - // auto conV = - json[jsonLabels.constrainedVertices].get>>(); - std::cout << "Loaded constrained vertices. Number of constrained " - "vertices found:" - << constrainedVertices.size() << std::endl; - } - - if (json.contains(jsonLabels.nodalForces)) { - auto f = std::unordered_map>( - json[jsonLabels.nodalForces]); - for (const auto &forces : f) { - nodalExternalForces[forces.first] = Vector6d(forces.second); - } - std::cout << "Loaded forces. Number of forces found:" - << nodalExternalForces.size() << std::endl; - } - - if (json.contains(jsonLabels.label)) { - label = json[jsonLabels.label]; - } - - if (json.contains(jsonLabels.meshLabel)) { - pMesh->setLabel(json[jsonLabels.meshLabel]); - } - - return true; - } - - bool save(const std::string &folderDirectory) const { - const std::filesystem::path pathFolderDirectory(folderDirectory); - if (!std::filesystem::is_directory(pathFolderDirectory)) { - std::cerr << "A folder directory is expected for saving the simulation " - "job. Exiting.." - << std::endl; - return false; - } - - bool returnValue = true; - const std::string meshFilename = - std::filesystem::absolute( - std::filesystem::canonical( - std::filesystem::path(pathFolderDirectory))) - .append(pMesh->getLabel() + ".ply").string(); - returnValue = pMesh->savePly(meshFilename); - nlohmann::json json; - json[jsonLabels.meshFilename] = meshFilename; - if (!constrainedVertices.empty()) { - json[jsonLabels.constrainedVertices] = constrainedVertices; - } - if (!nodalExternalForces.empty()) { - std::unordered_map> arrForces; - for (const auto &f : nodalExternalForces) { - arrForces[f.first] = f.second; - } - json[jsonLabels.nodalForces] = arrForces; - } - - if (!label.empty()) { - json[jsonLabels.label] = label; - } - if (!pMesh->getLabel().empty()) { - json[jsonLabels.meshLabel] = pMesh->getLabel(); - } - - std::string jsonFilename( - std::filesystem::path(pathFolderDirectory) - .append(label + "_" + pMesh->getLabel() + "_simulationJob.json").string()); - std::ofstream jsonFile(jsonFilename); - jsonFile << json; - std::cout << "Saved simulation job as:" << jsonFilename << std::endl; - - return returnValue; - } +#endif // POLYSCOPE_DEFINED }; namespace Eigen { @@ -300,6 +330,16 @@ struct SimulationResults { std::string labelPrefix{"deformed"}; inline static char deliminator{' '}; + std::vector getTranslationalDisplacements() const { + std::vector translationalDisplacements(displacements.size()); + std::transform(displacements.begin(), displacements.end(), + translationalDisplacements.begin(), [&](const Vector6d &d) { + return VectorType(d[0], d[1], d[2]); + }); + + return translationalDisplacements; + } + void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; } @@ -307,45 +347,6 @@ struct SimulationResults { return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel(); } - void unregister() const { - if (!polyscope::hasCurveNetwork(getLabel())) { - std::cerr << "No curve network registered with a name: " << getLabel() - << std::endl; - std::cerr << "Nothing to remove." << std::endl; - return; - } - polyscope::removeCurveNetwork(getLabel()); - } - void registerForDrawing() const { - polyscope::options::groundPlaneEnabled = false; - polyscope::view::upDir = polyscope::view::UpDir::ZUp; - const std::string branchName = "Branch:Polyscope"; - polyscope::options::programName = branchName; - if (!polyscope::state::initialized) { - polyscope::init(); - } /* else { - polyscope::removeAllStructures(); - }*/ - // const std::string undeformedMeshName = "Undeformed_" + label; - // polyscope::registerCurveNetwork( - // undeformedMeshName, mesh->getEigenVertices(), - // mesh->getEigenEdges()); - - const std::shared_ptr &mesh = job->pMesh; - polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(), - mesh->getEigenEdges()) - ->setRadius(0.0007, false); - Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3); - for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { - const Vector6d &nodalDisplacement = displacements[vi]; - nodalDisplacements.row(vi) = Eigen::Vector3d( - nodalDisplacement[0], nodalDisplacement[1], nodalDisplacement[2]); - } - polyscope::getCurveNetwork(getLabel()) - ->updateNodePositions(mesh->getEigenVertices() + nodalDisplacements); - - job->registerForDrawing(getLabel()); - } void saveDeformedModel() { VCGEdgeMesh m; @@ -376,6 +377,53 @@ struct SimulationResults { return errorNorm < 1e-10; // return eigenDisplacements.isApprox(nodalDisplacements); } + +#ifdef POLYSCOPE_DEFINED + void unregister() const { + if (!polyscope::hasCurveNetwork(getLabel())) { + std::cerr << "No curve network registered with a name: " << getLabel() + << std::endl; + std::cerr << "Nothing to remove." << std::endl; + return; + } + polyscope::removeCurveNetwork(getLabel()); + } + void registerForDrawing( + const std::optional &desiredColor = std::nullopt) const { + polyscope::options::groundPlaneEnabled = false; + polyscope::view::upDir = polyscope::view::UpDir::ZUp; + const std::string branchName = "Branch:Polyscope"; + polyscope::options::programName = branchName; + if (!polyscope::state::initialized) { + polyscope::init(); + } /* else { + polyscope::removeAllStructures(); + }*/ + // const std::string undeformedMeshName = "Undeformed_" + label; + // polyscope::registerCurveNetwork( + // undeformedMeshName, mesh->getEigenVertices(), + // mesh->getEigenEdges()); + + const std::shared_ptr &mesh = job->pMesh; + auto polyscopeHandle_deformedEdmeMesh = + polyscope::registerCurveNetwork(getLabel(), mesh->getEigenVertices(), + mesh->getEigenEdges()) + ->setRadius(0.002, true); + if (desiredColor.has_value()) { + polyscopeHandle_deformedEdmeMesh->setColor(desiredColor.value()); + } + Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3); + for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { + const Vector6d &nodalDisplacement = displacements[vi]; + nodalDisplacements.row(vi) = Eigen::Vector3d( + nodalDisplacement[0], nodalDisplacement[1], nodalDisplacement[2]); + } + polyscopeHandle_deformedEdmeMesh->updateNodePositions( + mesh->getEigenVertices() + nodalDisplacements); + + job->registerForDrawing(getLabel()); + } +#endif }; #endif // SIMULATIONHISTORY_HPP diff --git a/utilities.hpp b/utilities.hpp index 7918a0d..b1e5088 100644 --- a/utilities.hpp +++ b/utilities.hpp @@ -98,6 +98,11 @@ struct Vector6d : public std::array { return true; }); } + + Eigen::Vector3d getTranslation() const { + return Eigen::Vector3d(this->operator[](0), this->operator[](1), + this->operator[](2)); + } }; namespace Utilities { @@ -155,10 +160,7 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { } // namespace Utilities -// namespace ConfigurationFile { - -//} -//} // namespace ConfigurationFile +#ifdef POLYSCOPE_DEFINED #include "polyscope/curve_network.h" #include "polyscope/polyscope.h" @@ -204,7 +206,12 @@ inline void registerWorldAxes() { ->addEdgeColorQuantity(worldAxesColorName, axesColors) ->setEnabled(true); } +#endif +// namespace ConfigurationFile { + +//} +//} // namespace ConfigurationFile template void constructInverseMap(const T1 &map, T2 &oppositeMap) { assert(!map.empty());