From 55def2cfd505aadeda00bbe2c7a481be277dcc8b Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Fri, 30 Apr 2021 13:13:58 +0300 Subject: [PATCH] Refactoring --- drmsimulationmodel.cpp | 140 ++-- drmsimulationmodel.hpp | 10 +- edgemesh.cpp | 37 +- edgemesh.hpp | 99 ++- linearsimulationmodel.hpp | 1 + mesh.hpp | 8 +- polymesh.hpp | 22 +- reducedmodeloptimizer_structs.hpp | 549 +++++++++------ simulation_structs.hpp | 677 +++++++++++------- simulationmesh.cpp | 8 +- simulationmesh.hpp | 2 +- topologyenumerator.cpp | 1090 +++++++++++++++-------------- trianglepatterngeometry.cpp | 612 ++++++++++++---- trianglepatterngeometry.hpp | 282 ++------ utilities.hpp | 142 +++- vcgtrimesh.cpp | 5 +- vcgtrimesh.hpp | 4 +- 17 files changed, 2156 insertions(+), 1532 deletions(-) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 83f49e0..33b40c5 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -4,7 +4,10 @@ #include #include #include + +#ifdef ENABLE_OPENMP #include +#endif void DRMSimulationModel::runUnitTests() { @@ -42,7 +45,7 @@ void DRMSimulationModel::runUnitTests() settings.Dtini = 0.1; settings.xi = 0.9969; settings.maxDRMIterations = 0; - settings.totalResidualForcesNormThreshold = 1; + totalResidualForcesNormThreshold = 1; // settings.shouldDraw = true; settings.beVerbose = true; SimulationResults simpleBeam_simulationResults @@ -320,6 +323,7 @@ VectorType DRMSimulationModel::computeDerivativeOfR(const EdgeType &e, const VectorType derivativeOfSumOfNormals = derivativeOfNormal_j + derivativeOfNormal_jplus1; const VectorType sumOfNormals = normal_j + normal_jplus1; const double normOfSumOfNormals = sumOfNormals.Norm(); + assert(normOfSumOfNormals != 0); const double derivativeOfNormOfSumOfNormals = computeDerivativeOfNorm(sumOfNormals, derivativeOfSumOfNormals); @@ -761,7 +765,9 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( pMesh->EN(), std::vector>(4, {-1, Vector6d()})); // omp_lock_t writelock; // omp_init_lock(&writelock); - //#pragma omp parallel for schedule(static) num_threads(8) +#ifdef ENABLE_OPENMP +#pragma omp parallel for schedule(static) num_threads(5) +#endif for (int ei = 0; ei < pMesh->EN(); ei++) { const EdgeType &e = pMesh->edge[ei]; const SimulationMesh::VertexType &ev_j = *e.cV(0); @@ -1272,11 +1278,11 @@ void DRMSimulationModel::updateNodeNormal( const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0; if (!checkedForMaximumMoment && viHasMoments) { mSettings.shouldUseTranslationalKineticEnergyThreshold = true; - if (mSettings.beVerbose) { - std::cout << "Maximum moment reached.The Kinetic energy of the system will " - "be used as a convergence criterion" - << std::endl; - } +#ifdef POLYSCOPE_DEFINED + std::cout << "Maximum moment reached.The Kinetic energy of the system will " + "be used as a convergence criterion" + << std::endl; +#endif checkedForMaximumMoment = true; } @@ -1534,8 +1540,12 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptrnodes[vi].displacements; } history.numberOfSteps = mCurrentSimulationStep; - - return SimulationResults{true, pJob, history, displacements}; + SimulationResults results; + results.converged = true; + results.job = pJob; + results.history = history; + results.displacements = displacements; + return results; } void DRMSimulationModel::printCurrentState() const @@ -1684,7 +1694,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) ->addEdgeScalarQuantity("I3", I3); // Specify the callback - polyscope::state::userCallback = [&]() { + PolyscopeInterface::addUserCallback([&]() { // Since options::openImGuiWindowForUserCallback == true by default, // we can immediately start using ImGui commands to build a UI @@ -1699,7 +1709,7 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); ImGui::PopItemWidth(); - }; + }); if (!screenshotsFolder.empty()) { static bool firstDraw = true; @@ -1728,6 +1738,12 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &externalForce : pJob->nodalExternalForces) { + totalExternalForcesNorm += externalForce.second.norm(); + } + totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination + * totalExternalForcesNorm; // if (!pJob->nodalExternalForces.empty()) { // double externalForcesNorm = 0; @@ -1764,7 +1780,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrgetEigenVertices(), pMesh->getEigenEdges()); // registerWorldAxes(); @@ -1809,9 +1825,9 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrcurrentTotalPotentialEnergykN; // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; - if (mCurrentSimulationStep != 0) { - history.stepPulse(*pMesh); - } + // if (mCurrentSimulationStep != 0) { + // history.stepPulse(*pMesh); + // } if (std::isnan(pMesh->currentTotalKineticEnergy)) { if (mSettings.beVerbose) { @@ -1820,6 +1836,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr maxDRMIterations*/ { @@ -1835,6 +1852,7 @@ currentSimulationStep > maxDRMIterations*/ // shouldUseKineticEnergyThreshold = true; } +#endif if (mSettings.shouldCreatePlots && mCurrentSimulationStep % 10 == 0) { printCurrentState(); SimulationResultsReporter::createPlot( @@ -1847,20 +1865,21 @@ currentSimulationStep > maxDRMIterations*/ // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm // << std::endl; // Kinetic energy convergence - if ((mSettings.shouldUseTranslationalKineticEnergyThreshold || - mCurrentSimulationStep > 100000) && - pMesh->currentTotalTranslationalKineticEnergy < - mSettings.totalTranslationalKineticEnergyThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN - << " kNm" << std::endl; + 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; - } - break; +#endif + break; } // Residual forces norm convergence if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy @@ -1868,32 +1887,31 @@ currentSimulationStep > maxDRMIterations*/ mesh->previousTotalPotentialEnergykN > mesh->currentTotalPotentialEnergykN*/ /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { - if (pMesh->totalResidualForcesNorm < - mSettings.totalResidualForcesNormThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" - << pMesh->currentTotalPotentialEnergykN << " kNm" - << std::endl; + if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + if (mSettings.beVerbose) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" + << std::endl; + } + break; + // } } - break; - // } - } - // for (VertexType &v : mesh->vert) { - // Node &node = mesh->nodes[v]; - // node.displacements = node.displacements - node.velocity * Dt; - // } - // applyDisplacements(constrainedVertices); - // if (!pJob->nodalForcedDisplacements.empty()) { - // applyForcedDisplacements( + // 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; + // pJob->nodalForcedDisplacements); + // } + // updateElementalLengths(); + resetVelocities(); + Dt = Dt * mSettings.xi; + ++numOfDampings; } if (mSettings.debugMessages) { @@ -1927,16 +1945,20 @@ mesh->currentTotalPotentialEnergykN*/ std::cout << "Number of dampings:" << numOfDampings << endl; } // mesh.printVertexCoordinates(mesh.VN() / 2); +#ifdef POLYSCOPE_DEFINED if (settings.shouldDraw) { draw(); } +#endif if (results.converged) { + results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); + results.rotationalDisplacementQuaternion.resize(pMesh->VN()); for (int vi = 0; vi < pMesh->VN(); vi++) { const Node &node = pMesh->nodes[vi]; const Eigen::Vector3d initialNormal = node.initialNormal.ToEigenVector(); - const Eigen::Vector3d deformedNormal= - pMesh->vert[vi].cN().ToEigenVector(); + const Eigen::Vector3d deformedNormal + = pMesh->vert[vi].cN().ToEigenVector(); Eigen::Quaternion q; q.setFromTwoVectors(initialNormal, deformedNormal); @@ -1963,15 +1985,19 @@ mesh->currentTotalPotentialEnergykN*/ referenceT1_deformed.ToEigenVector()); results.rotationalDisplacementQuaternion[vi] = q * q_f1 * q_nr; + //Update the displacement vector to contain the euler angles + const Eigen::Vector3d eulerAngles + = results.rotationalDisplacementQuaternion[vi].toRotationMatrix().eulerAngles(0, 1, 2); + results.displacements[vi][3] = eulerAngles[0]; + results.displacements[vi][4] = eulerAngles[1]; + results.displacements[vi][5] = eulerAngles[2]; } + } - results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); - } - - if (mSettings.shouldCreatePlots) { + if (mSettings.shouldCreatePlots) { SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); - } + } return results; } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 9e4fc50..928ef08 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -29,8 +29,8 @@ public: bool beVerbose{false}; bool shouldCreatePlots{false}; int drawingStep{1}; - double totalTranslationalKineticEnergyThreshold{1e-7}; - double totalResidualForcesNormThreshold{1e-3}; + double totalTranslationalKineticEnergyThreshold{1e-10}; + double totalExternalForcesNormPercentageTermination{1e-3}; double Dtini{0.1}; double xi{0.9969}; int maxDRMIterations{0}; @@ -44,6 +44,7 @@ private: bool checkedForMaximumMoment{false}; double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; + double totalResidualForcesNormThreshold{1e-3}; matplot::line_handle plotHandle; std::vector plotYValues; size_t numOfDampings{0}; @@ -109,8 +110,9 @@ private: const std::unordered_map> &fixedVertices); +#ifdef POLYSCOPE_DEFINED void draw(const string &screenshotsFolder= {}); - +#endif void updateNodalInternalForce(Vector6d &nodalInternalForce, const Vector6d &elementInternalForce, @@ -206,7 +208,7 @@ private: executeSimulation(const std::shared_ptr &pJob, const Settings &settings = Settings()); - static void runUnitTests(); + void runUnitTests(); void setTotalResidualForcesNormThreshold(double value); }; diff --git a/edgemesh.cpp b/edgemesh.cpp index 436f49a..d1f4f2d 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -15,13 +15,19 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const { bool VCGEdgeMesh::save(const string &plyFilename) { - assert(std::filesystem::path(plyFilename).extension() == ".ply"); + std::string filename = plyFilename; + if (filename.empty()) { + filename = std::filesystem::current_path().append(getLabel() + ".ply"); + } else if (std::filesystem::is_directory(std::filesystem::path(plyFilename))) { + filename = std::filesystem::path(plyFilename).append(getLabel() + ".ply"); + } + assert(std::filesystem::path(filename).extension().string() == ".ply"); unsigned int mask = 0; mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - if (nanoply::NanoPlyWrapper::SaveModel(plyFilename.c_str(), *this, mask, false) - != 0) { + mask |= vcg::tri::io::Mask::IOM_VERTCOLOR; + if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), *this, mask, false) != 0) { return false; } return true; @@ -295,9 +301,20 @@ bool VCGEdgeMesh::copy(VCGEdgeMesh &mesh) { return true; } -void VCGEdgeMesh::removeDuplicateVertices() { +void VCGEdgeMesh::removeDuplicateVertices() +{ vcg::tri::Clean::MergeCloseVertex(*this, 0.000061524); - vcg::tri::Allocator::CompactEveryVector(*this); + vcg::tri::Allocator::CompactVertexVector(*this); + vcg::tri::UpdateTopology::VertexEdge(*this); +} + +void VCGEdgeMesh::removeDuplicateVertices( + vcg::tri::Allocator::PointerUpdater &pu_vertices, + vcg::tri::Allocator::PointerUpdater &pu_edges) +{ + vcg::tri::Clean::MergeCloseVertex(*this, 0.000061524); + vcg::tri::Allocator::CompactVertexVector(*this, pu_vertices); + vcg::tri::Allocator::CompactEdgeVector(*this, pu_edges); vcg::tri::UpdateTopology::VertexEdge(*this); } @@ -334,10 +351,6 @@ void VCGEdgeMesh::getVertices(Eigen::MatrixX3d &vertices) { } } -std::string VCGEdgeMesh::getLabel() const { return label; } - -void VCGEdgeMesh::setLabel(const std::string &value) { label = value; } - void VCGEdgeMesh::getEdges(Eigen::MatrixX2i &edges) { edges = Eigen::MatrixX2i(); edges.resize(EN(), 2); @@ -365,10 +378,10 @@ void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const { polyscope::CurveNetwork *VCGEdgeMesh::registerForDrawing( const std::optional> &desiredColor, const bool &shouldEnable) { - initPolyscope(); + PolyscopeInterface::init(); const double drawingRadius = 0.002; - polyscope::CurveNetwork *polyscopeHandle_edgeMesh( - polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges())); + polyscope::CurveNetwork *polyscopeHandle_edgeMesh + = polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges()); polyscopeHandle_edgeMesh->setEnabled(shouldEnable); polyscopeHandle_edgeMesh->setRadius(drawingRadius, true); diff --git a/edgemesh.hpp b/edgemesh.hpp index 46d3a21..c373add 100755 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -25,73 +25,72 @@ class VCGEdgeMeshEdgeType vcg::edge::BitFlags, vcg::edge::EEAdj, vcg::edge::VEAdj> {}; -class VCGEdgeMesh : public vcg::tri::TriMesh, - std::vector>, - Mesh { +class VCGEdgeMesh + : public vcg::tri::TriMesh, std::vector>, + public Mesh +{ +protected: + Eigen::MatrixX2i eigenEdges; + Eigen::MatrixX3d eigenVertices; + Eigen::MatrixX3d eigenEdgeNormals; - protected: - Eigen::MatrixX2i eigenEdges; - Eigen::MatrixX3d eigenVertices; - Eigen::MatrixX3d eigenEdgeNormals; - std::string label{"No_name"}; + void getEdges(Eigen::MatrixX2i &edges); + void getVertices(Eigen::MatrixX3d &vertices); - void getEdges(Eigen::MatrixX2i &edges); - void getVertices(Eigen::MatrixX3d &vertices); - - public: - VCGEdgeMesh(); - template - size_t getIndex(const MeshElement &meshElement) { - return vcg::tri::Index(*this, meshElement); - } - void updateEigenEdgeAndVertices(); - /* +public: + VCGEdgeMesh(); + template + size_t getIndex(const MeshElement &meshElement) + { + return vcg::tri::Index(*this, meshElement); + } + void updateEigenEdgeAndVertices(); + /* * The copy function shold be a virtual function of the base interface Mesh class. * https://stackoverflow.com/questions/2354210/can-a-class-member-function-template-be-virtual * use type erasure (?) * */ - bool copy(VCGEdgeMesh &mesh); + bool copy(VCGEdgeMesh &mesh); -void removeDuplicateVertices(); -void deleteDanglingVertices(); -void deleteDanglingVertices(vcg::tri::Allocator::PointerUpdater &pu); + void removeDuplicateVertices(); + void deleteDanglingVertices(); + void deleteDanglingVertices(vcg::tri::Allocator::PointerUpdater &pu); + void getEdges(Eigen::MatrixX3d &edgeStartingPoints, Eigen::MatrixX3d &edgeEndingPoints) const; - void getEdges(Eigen::MatrixX3d &edgeStartingPoints, - Eigen::MatrixX3d &edgeEndingPoints) const; + Eigen::MatrixX3d getNormals() const; - Eigen::MatrixX3d getNormals() const; + bool plyFileHasAllRequiredFields(const std::string &plyFilename); - bool plyFileHasAllRequiredFields(const std::string &plyFilename); + bool loadUsingNanoply(const std::string &plyFilename); - bool loadUsingNanoply(const std::string &plyFilename); + bool load(const std::string &plyFilename) override; + bool save(const std::string &plyFilename = std::string()) override; - bool load(const std::string &plyFilename) override; - bool save(const std::string &plyFilename) override; + bool createSpanGrid(const size_t squareGridDimension); + bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight); + void createSpiral(const float °reesOfArm, const size_t &numberOfSamples); - bool createSpanGrid(const size_t squareGridDimension); - bool createSpanGrid(const size_t desiredWidth, const size_t desiredHeight); - void createSpiral(const float °reesOfArm, const size_t &numberOfSamples); - - Eigen::MatrixX2i getEigenEdges() const; - Eigen::MatrixX3d getEigenVertices(); - Eigen::MatrixX3d getEigenEdgeNormals() const; - void printVertexCoordinates(const size_t &vi) const; + Eigen::MatrixX2i getEigenEdges() const; + Eigen::MatrixX3d getEigenVertices(); + Eigen::MatrixX3d getEigenEdgeNormals() const; + void printVertexCoordinates(const size_t &vi) const; #ifdef POLYSCOPE_DEFINED - polyscope::CurveNetwork *registerForDrawing( - const std::optional> &desiredColor = std::nullopt, - const bool &shouldEnable = true); - void unregister() const; - void drawInitialFrames(polyscope::CurveNetwork *polyscopeHandle_initialMesh) const; + polyscope::CurveNetwork *registerForDrawing( + const std::optional> &desiredColor = std::nullopt, + const bool &shouldEnable = true); + void unregister() const; + void drawInitialFrames(polyscope::CurveNetwork *polyscopeHandle_initialMesh) const; #endif - std::string getLabel() const; - void setLabel(const std::string &value); + void removeDuplicateVertices( + vcg::tri::Allocator::PointerUpdater &pu_vertices, + vcg::tri::Allocator::PointerUpdater &pu_edges); - private: - void GeneratedRegularSquaredPattern( - const double angleDeg, std::vector> &pattern, - const size_t &desiredNumberOfSamples); - }; +private: + void GeneratedRegularSquaredPattern(const double angleDeg, + std::vector> &pattern, + const size_t &desiredNumberOfSamples); +}; using VectorType = VCGEdgeMesh::CoordType; using CoordType = VCGEdgeMesh::CoordType; diff --git a/linearsimulationmodel.hpp b/linearsimulationmodel.hpp index 26a5c94..5086bfb 100644 --- a/linearsimulationmodel.hpp +++ b/linearsimulationmodel.hpp @@ -244,6 +244,7 @@ public: + elementPotentialEnergy_torsion; results.internalPotentialEnergy += elementInternalPotentialEnergy; } + results.job = simulationJob; return results; } diff --git a/mesh.hpp b/mesh.hpp index f2d6b7d..c08ab86 100755 --- a/mesh.hpp +++ b/mesh.hpp @@ -7,7 +7,7 @@ class Mesh { protected: - std::string label; + std::string label{"empty_label"}; public: virtual ~Mesh() = default; @@ -16,10 +16,16 @@ public: std::string getLabel() const; void setLabel(const std::string &newLabel); + void prependToLabel(const std::string &text); }; inline std::string Mesh::getLabel() const { return label; } inline void Mesh::setLabel(const std::string &newLabel) { label = newLabel; } +inline void Mesh::prependToLabel(const std::string &text) +{ + label = text + label; +} + #endif // MESH_HPP diff --git a/polymesh.hpp b/polymesh.hpp index 843b2cc..931fca6 100755 --- a/polymesh.hpp +++ b/polymesh.hpp @@ -3,9 +3,11 @@ #include "mesh.hpp" #include "vcg/complex/complex.h" #include -#include #include +#ifdef POLYSCOPE_DEFINED +#include +#endif class PFace; class PVertex; class PEdge; @@ -29,10 +31,7 @@ class PEdge : public vcg::Edge + vcg::edge::EVAdj> {}; class PFace : public vcg::Face polyscopeHandle_mesh( polyscope::registerSurfaceMesh(label, vertices, faces)); @@ -191,6 +192,9 @@ public: // } vcg::tri::UpdateTopology::VertexEdge(*this); vcg::tri::UpdateTopology::EdgeEdge(*this); + vcg::tri::UpdateTopology::FaceFace(*this); + vcg::tri::UpdateTopology::VertexFace(*this); + vcg::tri::UpdateTopology::AllocateEdge(*this); // vcg::tri::UpdateTopology::VertexFace(*this); return true; @@ -205,4 +209,6 @@ public: } }; +using ConstVCGPolyMesh = VCGPolyMesh; + #endif // POLYMESH_HPP diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index add823c..fa881f5 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -39,65 +39,161 @@ struct Colors return label + "=[" + std::to_string(min) + "," + std::to_string(max) + "]"; } + + void fromString(const std::string &s) + { + const std::size_t equalPos = s.find("="); + label = s.substr(0, equalPos); + const std::size_t commaPos = s.find(","); + const size_t minBeginPos = equalPos + 2; + min = std::stod(s.substr(minBeginPos, commaPos - minBeginPos)); + const size_t maxBeginPos = commaPos + 1; + const std::size_t closingBracketPos = s.find("]"); + max = std::stod(s.substr(maxBeginPos, closingBracketPos - maxBeginPos)); + } + + bool operator==(const xRange &xrange) const + { + return label == xrange.label && min == xrange.min && max == xrange.max; + } + }; + struct Settings + { + enum NormalizationStrategy { NonNormalized, Epsilon }; + std::vector xRanges; + inline static vector normalizationStrategyStrings{"NonNormalized", "Epsilon"}; + int numberOfFunctionCalls{100000}; + double solverAccuracy{1e-2}; + NormalizationStrategy normalizationStrategy{Epsilon}; + double normalizationParameter{0.003}; + struct ObjectiveWeights + { + double translational{1}; + double rotational{1}; + } objectiveWeights; + + struct JsonKeys + { + inline static std::string filename{"OptimizationSettings.json"}; + inline static std::string OptimizationVariables{"OptimizationVariables"}; + inline static std::string NumberOfFunctionCalls{"NumberOfFunctionCalls"}; + inline static std::string SolverAccuracy{"SolverAccuracy"}; + inline static std::string TranslationalObjectiveWeight{"TransObjWeight"}; + }; + + void save(const std::filesystem::path &saveToPath) + { + assert(std::filesystem::is_directory(saveToPath)); + + nlohmann::json json; + //write x ranges + for (size_t xRangeIndex = 0; xRangeIndex < xRanges.size(); xRangeIndex++) { + const auto &xRange = xRanges[xRangeIndex]; + json[JsonKeys::OptimizationVariables + "_" + std::to_string(xRangeIndex)] + = xRange.toString(); + } + + json[JsonKeys::NumberOfFunctionCalls] = numberOfFunctionCalls; + json[JsonKeys::SolverAccuracy] = solverAccuracy; + json[JsonKeys::TranslationalObjectiveWeight] = objectiveWeights.translational; + + std::filesystem::path jsonFilePath( + std::filesystem::path(saveToPath).append(JsonKeys::filename)); + std::ofstream jsonFile(jsonFilePath.string()); + jsonFile << json; + } + + bool load(const std::filesystem::path &loadFromPath) + { + assert(std::filesystem::is_directory(loadFromPath)); + //Load optimal X + nlohmann::json json; + std::filesystem::path jsonFilepath( + std::filesystem::path(loadFromPath).append(JsonKeys::filename)); + if (!std::filesystem::exists(jsonFilepath)) { + return false; + } + std::ifstream ifs(jsonFilepath); + ifs >> json; + + //read x ranges + size_t xRangeIndex = 0; + while (true) { + const std::string jsonXRangeKey = JsonKeys::OptimizationVariables + "_" + + std::to_string(xRangeIndex); + if (!json.contains(jsonXRangeKey)) { + break; + } + xRange x; + x.fromString(json.at(jsonXRangeKey)); + xRangeIndex++; + } + + numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); + solverAccuracy = json.at(JsonKeys::SolverAccuracy); + objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); + return true; + } + + std::string toString() const + { + std::string settingsString; + if (!xRanges.empty()) { + std::string xRangesString; + for (const xRange &range : xRanges) { + xRangesString += range.toString() + " "; + } + settingsString += xRangesString; + } + settingsString += "FuncCalls=" + std::to_string(numberOfFunctionCalls) + + " Accuracy=" + std::to_string(solverAccuracy) + + " TransWeight=" + std::to_string(objectiveWeights.translational) + + " RotWeight=" + std::to_string(objectiveWeights.rotational); + + return settingsString; + } + + void writeHeaderTo(csvFile &os) const + { + if (!xRanges.empty()) { + for (const xRange &range : xRanges) { + os << range.label + " max"; + os << range.label + " min"; + } + } + os << "Function Calls"; + os << "Solution Accuracy"; + os << "Normalization strategy"; + os << "Trans weight"; + os << "Rot weight"; + // os << std::endl; + } + + void writeSettingsTo(csvFile &os) const + { + if (!xRanges.empty()) { + for (const xRange &range : xRanges) { + os << range.max; + os << range.min; + } + } + os << numberOfFunctionCalls; + os << solverAccuracy; + os << normalizationStrategyStrings[normalizationStrategy] + "_" + + std::to_string(normalizationParameter); + os << objectiveWeights.translational; + os << objectiveWeights.rotational; + } }; - struct Settings { - enum NormalizationStrategy { - NonNormalized, - Epsilon, - MaxDisplacement, - EqualDisplacements - }; - inline static vector normalizationStrategyStrings{ - "NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"}; - std::vector xRanges; - int numberOfFunctionCalls{100}; - double solutionAccuracy{1e-2}; - NormalizationStrategy normalizationStrategy{Epsilon}; - double normalizationParameter{0.003}; - std::string toString() const { - std::string settingsString; - if (!xRanges.empty()) { - std::string xRangesString; - for (const xRange &range : xRanges) { - xRangesString += range.toString() + " "; - } - settingsString += xRangesString; - } - settingsString += - "FuncCalls=" + std::to_string(numberOfFunctionCalls) + - " Accuracy=" + std::to_string(solutionAccuracy) + - " Norm=" + normalizationStrategyStrings[normalizationStrategy]; - - return settingsString; - } - - void writeHeaderTo(csvFile &os) const { - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - os << range.label + " max"; - os << range.label + " min"; - } - } - os << "Function Calls"; - os << "Solution Accuracy"; - os << "Normalization strategy"; - // os << std::endl; - } - - void writeSettingsTo(csvFile &os) const { - if (!xRanges.empty()) { - for (const xRange &range : xRanges) { - os << range.max; - os << range.min; - } - } - os << numberOfFunctionCalls; - os << solutionAccuracy; - os << normalizationStrategyStrings[normalizationStrategy] + "_" + - std::to_string(normalizationParameter); - } - }; + inline bool operator==(const Settings &settings1, const Settings &settings2) noexcept + { + return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls + && settings1.xRanges == settings2.xRanges + && settings1.solverAccuracy == settings2.solverAccuracy + && settings1.objectiveWeights.translational + == settings2.objectiveWeights.translational; + } inline void updateMeshWithOptimalVariables(const std::vector &x, SimulationMesh &m) { @@ -153,9 +249,10 @@ struct Colors struct Results { + std::string label{"EmptyLabel"}; double time{-1}; int numberOfSimulationCrashes{0}; - + Settings settings; struct ObjectiveValues { double totalRaw; @@ -166,28 +263,34 @@ struct Colors std::vector perSimulationScenario_rotational; std::vector perSimulationScenario_total; } objectiveValue; - // std::vector> optimalXNameValuePairs; std::vector> optimalXNameValuePairs; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; + + PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern vcg::Triangle3 baseTriangle; + struct JsonKeys { inline static std::string filename{"OptimizationResults.json"}; inline static std::string optimizationVariables{"OptimizationVariables"}; inline static std::string baseTriangle{"BaseTriangle"}; + inline static std::string Label{"Label"}; + inline static std::string FullPatternLabel{"FullPatternLabel"}; + inline static std::string Settings{"OptimizationSettings"}; }; - void save(const string &saveToPath) const + void save(const string &saveToPath) { - assert(std::filesystem::is_directory(saveToPath)); + std::filesystem::create_directories(saveToPath); //clear directory for (const auto &entry : std::filesystem::directory_iterator(saveToPath)) std::filesystem::remove_all(entry.path()); //Save optimal X nlohmann::json json_optimizationResults; + json_optimizationResults[JsonKeys::Label] = label; std::string jsonValue_optimizationVariables; for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { jsonValue_optimizationVariables.append(optimalXNameValuePair.first + ","); @@ -209,6 +312,8 @@ struct Colors baseTriangle.cP2(0)[0], baseTriangle.cP2(0)[1], baseTriangle.cP2(0)[2]}; + baseTriangleFullPattern.save(std::filesystem::path(saveToPath)); + json_optimizationResults[JsonKeys::FullPatternLabel] = baseTriangleFullPattern.getLabel(); ////Save to json file std::filesystem::path jsonFilePath( std::filesystem::path(saveToPath).append(JsonKeys::filename)); @@ -219,15 +324,16 @@ struct Colors //Save jobs and meshes for each simulation scenario //Save the reduced and full patterns - + const std::filesystem::path simulationJobsPath( + std::filesystem::path(saveToPath).append("SimulationJobs")); const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationJobs; simulationScenarioIndex++) { const std::shared_ptr &pFullPatternSimulationJob = fullPatternSimulationJobs[simulationScenarioIndex]; std::filesystem::path simulationJobFolderPath( - std::filesystem::path(saveToPath).append(pFullPatternSimulationJob->label)); - std::filesystem::create_directory(simulationJobFolderPath); + std::filesystem::path(simulationJobsPath).append(pFullPatternSimulationJob->label)); + std::filesystem::create_directories(simulationJobFolderPath); const auto fullPatternDirectoryPath = std::filesystem::path(simulationJobFolderPath) .append("Full"); std::filesystem::create_directory(fullPatternDirectoryPath); @@ -242,7 +348,102 @@ struct Colors pReducedPatternSimulationJob->save(reducedPatternDirectoryPath.string()); } + settings.save(saveToPath); + +#ifdef POLYSCOPE_DEFINED std::cout << "Saved optimization results to:" << saveToPath << std::endl; +#endif + } + + bool load(const string &loadFromPath) + { + assert(std::filesystem::is_directory(loadFromPath)); + std::filesystem::path jsonFilepath( + std::filesystem::path(loadFromPath).append(JsonKeys::filename)); + if (!std::filesystem::exists(jsonFilepath)) { + return false; + } + //Load optimal X + nlohmann::json json_optimizationResults; + std::ifstream ifs(std::filesystem::path(loadFromPath).append(JsonKeys::filename)); + ifs >> json_optimizationResults; + + label = json_optimizationResults.at(JsonKeys::Label); + std::string optimizationVariablesString = *json_optimizationResults.find( + JsonKeys::optimizationVariables); + std::string optimizationVariablesDelimeter = ","; + size_t pos = 0; + std::vector optimizationVariablesNames; + while ((pos = optimizationVariablesString.find(optimizationVariablesDelimeter)) + != std::string::npos) { + const std::string variableName = optimizationVariablesString.substr(0, pos); + optimizationVariablesNames.push_back(variableName); + optimizationVariablesString.erase(0, pos + optimizationVariablesDelimeter.length()); + } + optimizationVariablesNames.push_back(optimizationVariablesString); //add last variable name + + optimalXNameValuePairs.resize(optimizationVariablesNames.size()); + for (int xVariable_index = 0; xVariable_index < optimizationVariablesNames.size(); + xVariable_index++) { + const std::string xVariable_name = optimizationVariablesNames[xVariable_index]; + const double xVariable_value = *json_optimizationResults.find(xVariable_name); + optimalXNameValuePairs[xVariable_index] = std::make_pair(xVariable_name, + xVariable_value); + } + + const std::string fullPatternLabel = json_optimizationResults.at( + JsonKeys::FullPatternLabel); + baseTriangleFullPattern.load( + std::filesystem::path(loadFromPath).append(fullPatternLabel + ".ply")); + + std::vector baseTriangleVertices = json_optimizationResults.at( + JsonKeys::baseTriangle); + baseTriangle.P0(0) = CoordType(baseTriangleVertices[0], + baseTriangleVertices[1], + baseTriangleVertices[2]); + baseTriangle.P1(0) = CoordType(baseTriangleVertices[3], + baseTriangleVertices[4], + baseTriangleVertices[5]); + baseTriangle.P2(0) = CoordType(baseTriangleVertices[6], + baseTriangleVertices[7], + baseTriangleVertices[8]); + + const std::filesystem::path simulationJobsFolderPath( + std::filesystem::path(loadFromPath).append("SimulationJobs")); + for (const auto &directoryEntry : + filesystem::directory_iterator(simulationJobsFolderPath)) { + const auto simulationScenarioPath = directoryEntry.path(); + if (!std::filesystem::is_directory(simulationScenarioPath)) { + continue; + } + // Load full pattern files + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath).append("Full"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + + job.load(filepath.string()); + fullPatternSimulationJobs.push_back(std::make_shared(job)); + } + } + + // Load reduced pattern files and apply optimal parameters + for (const auto &fileEntry : filesystem::directory_iterator( + std::filesystem::path(simulationScenarioPath).append("Reduced"))) { + const auto filepath = fileEntry.path(); + if (filepath.extension() == ".json") { + SimulationJob job; + job.load(filepath.string()); + applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh); + applyOptimizationResults_elements(*this, job.pMesh); + reducedPatternSimulationJobs.push_back(std::make_shared(job)); + } + } + } + settings.load(loadFromPath); + + return true; } template @@ -255,19 +456,38 @@ struct Colors optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(), reducedPattern_optimizationResults.optimalXNameValuePairs.end()); assert(optimalXVariables.contains("HexSize") && optimalXVariables.contains("HexAngle")); + applyOptimizationResults_innerHexagon(optimalXVariables.at("HexSize"), + optimalXVariables.at("HexAngle"), + patternBaseTriangle, + reducedPattern); + } + template + static void applyOptimizationResults_innerHexagon( + const double &hexSize, + const double &hexAngle, + const vcg::Triangle3 &patternBaseTriangle, + MeshType &reducedPattern) + { //Set optimal geometrical params of the reduced pattern - CoordType center_barycentric(1, 0, 0); - CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); - CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) - * optimalXVariables.at("HexSize")); - reducedPattern.vert[1].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0] + // CoordType center_barycentric(1, 0, 0); + // CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); + // CoordType movableVertex_barycentric( + // (center_barycentric * (1 - hexSize) + interfaceEdgeMiddle_barycentric)); + CoordType movableVertex_barycentric(1 - hexSize, hexSize / 2, hexSize / 2); + + reducedPattern.vert[0].P() = patternBaseTriangle.cP(0) * movableVertex_barycentric[0] + patternBaseTriangle.cP(1) * movableVertex_barycentric[1] + patternBaseTriangle.cP(2) * movableVertex_barycentric[2]; - reducedPattern.vert[1].P() = vcg::RotationMatrix(CoordType{0, 0, 1}, - vcg::math::ToRad( - optimalXVariables.at("HexAngle"))) - * reducedPattern.vert[1].cP(); + if (hexAngle != 0) { + reducedPattern.vert[0].P() = vcg::RotationMatrix(CoordType{0, 0, 1}, + vcg::math::ToRad(hexAngle)) + * reducedPattern.vert[0].cP(); + } + // CoordType p0 = reducedPattern.vert[0].P(); + // CoordType p1 = reducedPattern.vert[1].P(); + // int i = 0; + // i++; } static void applyOptimizationResults_elements( @@ -314,139 +534,62 @@ struct Colors pTiledReducedPattern_simulationMesh->reset(); } - void load(const string &loadFromPath) - { - assert(std::filesystem::is_directory(loadFromPath)); - //Load optimal X - nlohmann::json json_optimizationResults; - std::ifstream ifs(std::filesystem::path(loadFromPath).append(JsonKeys::filename)); - ifs >> json_optimizationResults; - - std::string optimizationVariablesString = *json_optimizationResults.find( - JsonKeys::optimizationVariables); - std::string optimizationVariablesDelimeter = ","; - size_t pos = 0; - std::vector optimizationVariablesNames; - while ((pos = optimizationVariablesString.find(optimizationVariablesDelimeter)) - != std::string::npos) { - const std::string variableName = optimizationVariablesString.substr(0, pos); - optimizationVariablesNames.push_back(variableName); - optimizationVariablesString.erase(0, pos + optimizationVariablesDelimeter.length()); - } - optimizationVariablesNames.push_back(optimizationVariablesString); //add last variable name - - optimalXNameValuePairs.resize(optimizationVariablesNames.size()); - for (int xVariable_index = 0; xVariable_index < optimizationVariablesNames.size(); - xVariable_index++) { - const std::string xVariable_name = optimizationVariablesNames[xVariable_index]; - const double xVariable_value = *json_optimizationResults.find(xVariable_name); - optimalXNameValuePairs[xVariable_index] = std::make_pair(xVariable_name, - xVariable_value); - } - - std::vector baseTriangleVertices = json_optimizationResults.at( - JsonKeys::baseTriangle); - baseTriangle.P0(0) = CoordType(baseTriangleVertices[0], - baseTriangleVertices[1], - baseTriangleVertices[2]); - baseTriangle.P1(0) = CoordType(baseTriangleVertices[3], - baseTriangleVertices[4], - baseTriangleVertices[5]); - baseTriangle.P2(0) = CoordType(baseTriangleVertices[6], - baseTriangleVertices[7], - baseTriangleVertices[8]); - - for (const auto &directoryEntry : filesystem::directory_iterator(loadFromPath)) { - const auto simulationScenarioPath = directoryEntry.path(); - if (!std::filesystem::is_directory(simulationScenarioPath)) { - continue; - } - // Load reduced pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath).append("Full"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - - job.load(filepath.string()); - fullPatternSimulationJobs.push_back(std::make_shared(job)); - } - } - - // Load full pattern files - for (const auto &fileEntry : filesystem::directory_iterator( - std::filesystem::path(simulationScenarioPath).append("Reduced"))) { - const auto filepath = fileEntry.path(); - if (filepath.extension() == ".json") { - SimulationJob job; - job.load(filepath.string()); - applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh); - applyOptimizationResults_elements(*this, job.pMesh); - reducedPatternSimulationJobs.push_back(std::make_shared(job)); - } - } - } - } #if POLYSCOPE_DEFINED void draw() const { - initPolyscope(); - DRMSimulationModel drmSimulator; - LinearSimulationModel linearSimulator; - assert(fullPatternSimulationJobs.size() == - reducedPatternSimulationJobs.size()); - fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - Colors::fullInitial); - reducedPatternSimulationJobs[0]->pMesh->registerForDrawing( - Colors::reducedInitial, false); + PolyscopeInterface::init(); + DRMSimulationModel drmSimulator; + LinearSimulationModel linearSimulator; + assert(fullPatternSimulationJobs.size() == reducedPatternSimulationJobs.size()); + fullPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::fullInitial); + reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial, false); - const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); - for (int simulationJobIndex = 0; - simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) { - // Drawing of full pattern results - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationJobIndex]; - pFullPatternSimulationJob->registerForDrawing( - fullPatternSimulationJobs[0]->pMesh->getLabel()); - SimulationResults fullModelResults = drmSimulator.executeSimulation( - pFullPatternSimulationJob); - fullModelResults.registerForDrawing(Colors::fullDeformed, true, true); - // SimulationResults fullModelLinearResults = - // linearSimulator.executeSimulation(pFullPatternSimulationJob); - // fullModelLinearResults.setLabelPrefix("linear"); - // fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false); - // Drawing of reduced pattern results - const std::shared_ptr &pReducedPatternSimulationJob - = reducedPatternSimulationJobs[simulationJobIndex]; - // SimulationResults reducedModelResults = drmSimulator.executeSimulation( - // pReducedPatternSimulationJob); - // reducedModelResults.registerForDrawing(Colors::reducedDeformed, false); - SimulationResults reducedModelLinearResults = - linearSimulator.executeSimulation(pReducedPatternSimulationJob); - reducedModelLinearResults.setLabelPrefix("linear"); - reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, true, true); - polyscope::options::programName = - fullPatternSimulationJobs[0]->pMesh->getLabel(); - polyscope::view::resetCameraToHomeView(); - polyscope::show(); - // Save a screensh - const std::string screenshotFilename = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/Results/Images/" + - fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" + - pFullPatternSimulationJob->getLabel(); - polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - // reducedModelResults.unregister(); - reducedModelLinearResults.unregister(); - // fullModelLinearResults.unregister(); - // double error = computeError( - // reducedModelResults.displacements,fullModelResults.displacements, - // ); - // std::cout << "Error of simulation scenario " - // << - // simula simulationScenarioStrings[simulationScenarioIndex] - // << " is " - // << error << std::endl; + const int numberOfSimulationJobs = fullPatternSimulationJobs.size(); + for (int simulationJobIndex = 0; simulationJobIndex < numberOfSimulationJobs; + simulationJobIndex++) { + // Drawing of full pattern results + const std::shared_ptr &pFullPatternSimulationJob + = fullPatternSimulationJobs[simulationJobIndex]; + pFullPatternSimulationJob->registerForDrawing( + fullPatternSimulationJobs[0]->pMesh->getLabel()); + SimulationResults fullModelResults = drmSimulator.executeSimulation( + pFullPatternSimulationJob); + fullModelResults.registerForDrawing(Colors::fullDeformed, true, true); + // SimulationResults fullModelLinearResults = + // linearSimulator.executeSimulation(pFullPatternSimulationJob); + // fullModelLinearResults.setLabelPrefix("linear"); + // fullModelLinearResults.registerForDrawing(Colors::fullDeformed,false); + // Drawing of reduced pattern results + const std::shared_ptr &pReducedPatternSimulationJob + = reducedPatternSimulationJobs[simulationJobIndex]; + // SimulationResults reducedModelResults = drmSimulator.executeSimulation( + // pReducedPatternSimulationJob); + // reducedModelResults.registerForDrawing(Colors::reducedDeformed, false); + SimulationResults reducedModelLinearResults = linearSimulator.executeSimulation( + pReducedPatternSimulationJob); + reducedModelLinearResults.setLabelPrefix("linear"); + reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, true, true); + polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel(); + polyscope::view::resetCameraToHomeView(); + polyscope::show(); + // Save a screensh + const std::string screenshotFilename + = "/home/iason/Coding/Projects/Approximating shapes with flat " + "patterns/RodModelOptimizationForPatterns/Results/Images/" + + fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" + + pFullPatternSimulationJob->getLabel(); + polyscope::screenshot(screenshotFilename, false); + fullModelResults.unregister(); + // reducedModelResults.unregister(); + reducedModelLinearResults.unregister(); + // fullModelLinearResults.unregister(); + // double error = computeError( + // reducedModelResults.displacements,fullModelResults.displacements, + // ); + // std::cout << "Error of simulation scenario " + // << + // simula simulationScenarioStrings[simulationScenarioIndex] + // << " is " + // << error << std::endl; } } #endif // POLYSCOPE_DEFINED diff --git a/simulation_structs.hpp b/simulation_structs.hpp index 159dae5..fd75566 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -1,4 +1,4 @@ -#ifndef SIMULATIONSTRUCTS_HPP +#ifndef SIMULATIONSTRUCTS_HPP #define SIMULATIONSTRUCTS_HPP #include "simulationmesh.hpp" @@ -57,6 +57,7 @@ class SimulationJob { // json labels struct JSONLabels { inline static std::string meshFilename{"mesh filename"}; + inline static std::string forcedDisplacements{"forced displacements"}; inline static std::string constrainedVertices{"fixed vertices"}; inline static std::string nodalForces{"forces"}; inline static std::string label{"label"}; @@ -81,6 +82,33 @@ public: SimulationJob() {} SimulationJob(const std::string &jsonFilename) { load(jsonFilename); } + void remap(const std::unordered_map &thisToOtherViMap, + SimulationJob &simulationJobMapped) const + { + assert(simulationJobMapped.pMesh->VN() != 0); + std::unordered_map> reducedModelFixedVertices; + for (auto fullModelFixedVertex : this->constrainedVertices) { + reducedModelFixedVertices[thisToOtherViMap.at(fullModelFixedVertex.first)] + = fullModelFixedVertex.second; + } + + std::unordered_map reducedModelNodalForces; + for (auto fullModelNodalForce : this->nodalExternalForces) { + reducedModelNodalForces[thisToOtherViMap.at(fullModelNodalForce.first)] + = fullModelNodalForce.second; + } + + std::unordered_map reducedNodalForcedDisplacements; + for (auto fullForcedDisplacement : this->nodalForcedDisplacements) { + reducedNodalForcedDisplacements[thisToOtherViMap.at(fullForcedDisplacement.first)] + = fullForcedDisplacement.second; + } + + simulationJobMapped.constrainedVertices = reducedModelFixedVertices; + simulationJobMapped.nodalExternalForces = reducedModelNodalForces; + simulationJobMapped.label = this->getLabel(); + simulationJobMapped.nodalForcedDisplacements = reducedNodalForcedDisplacements; + } SimulationJob getCopy() const { SimulationJob jobCopy; jobCopy.pMesh = std::make_shared(); @@ -109,198 +137,270 @@ public: return json.dump(); } - bool load(const std::string &jsonFilename) { - const bool beVerbose = false; - 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; - } - - if (beVerbose) { - 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->load(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>>(); - if (beVerbose) { - std::cout << "Loaded constrained vertices. Number of constrained " - "vertices found:" - << constrainedVertices.size() << std::endl; + void clear() + { + label = "empty_job"; + constrainedVertices.clear(); + nodalExternalForces.clear(); + nodalForcedDisplacements.clear(); + if (pMesh.use_count() == 1) { + std::cout << "Job mesh is deleted" << 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); - } - if (beVerbose) { - 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; + pMesh.reset(); } - 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->save(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; + bool load(const std::string &jsonFilename, const bool &shouldLoadMesh = true) + { + label = "empty_job"; + constrainedVertices.clear(); + nodalExternalForces.clear(); + nodalForcedDisplacements.clear(); + const bool beVerbose = false; + 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; } - 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; + 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; + } - return returnValue; + if (beVerbose) { + std::cout << "Loading json file:" << jsonFilename << std::endl; + } + nlohmann::json json; + std::ifstream ifs(jsonFilename); + ifs >> json; + + if (shouldLoadMesh) { + pMesh.reset(); + 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->load(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.meshLabel)) { + pMesh->setLabel(json[jsonLabels.meshLabel]); + } + } + + if (json.contains(jsonLabels.constrainedVertices)) { + constrainedVertices = + // auto conV = + json[jsonLabels.constrainedVertices] + .get>>(); + if (beVerbose) { + 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); + } + if (beVerbose) { + std::cout << "Loaded forces. Number of forces found:" << nodalExternalForces.size() + << std::endl; + } + } + + if (json.contains(jsonLabels.forcedDisplacements)) { + // auto conV = + std::unordered_map> forcedDisp + = json[jsonLabels.forcedDisplacements] + .get>>(); + + for (const auto &fd : forcedDisp) { + nodalForcedDisplacements[fd.first] = Eigen::Vector3d(fd.second[0], + fd.second[1], + fd.second[2]); + } + if (beVerbose) { + std::cout << "Loaded forced displacements. Number of forced displaced" + "vertices found:" + << nodalForcedDisplacements.size() << std::endl; + } + } + + if (json.contains(jsonLabels.label)) { + label = json[jsonLabels.label]; + } + + 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") + .append("SimulationJob.json") + .string()); + + const std::string meshFilename = std::filesystem::absolute( + std::filesystem::canonical( + std::filesystem::path(pathFolderDirectory))) + .append(pMesh->getLabel() + ".ply") + .string(); + returnValue = pMesh->save(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 (!nodalForcedDisplacements.empty()) { + std::unordered_map> forcedDisp; + for (const auto &fd : nodalForcedDisplacements) { + forcedDisp[fd.first] = {fd.second[0], fd.second[1], fd.second[2]}; + } + + json[jsonLabels.forcedDisplacements] = forcedDisp; + } + + 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()) { - assert(false); - std::cerr << "Expects a mesh label on which to draw the simulation job." - << std::endl; - return; - } - - auto structs = polyscope::state::structures; - if (!polyscope::hasCurveNetwork(meshLabel)) { - assert(false); - std::cerr << "Expects mesh already being registered to draw the " - "simulation job. No struct named " + - meshLabel - << std::endl; - return; - } - std::vector> nodeColors(pMesh->VN()); - for (auto fixedVertex : constrainedVertices) { - nodeColors[fixedVertex.first] = {0, 0, 1}; - } - if (!nodalForcedDisplacements.empty()) { - for (std::pair viDisplPair : - nodalForcedDisplacements) { - const VertexIndex vi = viDisplPair.first; - nodeColors[vi][0] += 1; - nodeColors[vi][0] /= 2; - nodeColors[vi][1] += 0; - nodeColors[vi][1] /= 2; - nodeColors[vi][2] += 0; - nodeColors[vi][2] /= 2; + void registerForDrawing(const std::string &meshLabel, const bool &shouldEnable = false) const + { + PolyscopeInterface::init(); + if (meshLabel.empty()) { + assert(false); + std::cerr << "Expects a mesh label on which to draw the simulation job." << std::endl; + return; } - } - std::for_each(nodeColors.begin(), nodeColors.end(), - [](std::array &color) { - const double norm = - sqrt(std::pow(color[0], 2) + std::pow(color[1], 2) + - std::pow(color[2], 2)); - if (norm > std::pow(10, -7)) { - color[0] /= norm; - color[1] /= norm; - color[2] /= norm; - } - }); - if (!nodeColors.empty()) { - polyscope::getCurveNetwork(meshLabel) - ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors) - ->setEnabled(false); - } + if (!polyscope::hasCurveNetwork(meshLabel)) { + assert(false); + std::cerr << "Expects mesh already being registered to draw the " + "simulation job. No struct named " + + meshLabel + << std::endl; + return; + } + std::vector> nodeColors(pMesh->VN()); + for (auto fixedVertex : constrainedVertices) { + const bool hasRotationalDoFConstrained = fixedVertex.second.contains(3) + || fixedVertex.second.contains(4) + || fixedVertex.second.contains(5); + const bool hasTranslationalDoFConstrained = fixedVertex.second.contains(0) + || fixedVertex.second.contains(1) + || fixedVertex.second.contains(2); + if (hasTranslationalDoFConstrained && !hasRotationalDoFConstrained) { + nodeColors[fixedVertex.first] = {0, 0, 1}; + } else if (!hasTranslationalDoFConstrained && hasRotationalDoFConstrained) { + nodeColors[fixedVertex.first] = {0, 1, 0}; + } else { + nodeColors[fixedVertex.first] = {0, 1, 1}; + } + } + if (!nodalForcedDisplacements.empty()) { + for (std::pair viDisplPair : nodalForcedDisplacements) { + const VertexIndex vi = viDisplPair.first; + nodeColors[vi][0] += 1; + nodeColors[vi][0] /= 2; + nodeColors[vi][1] += 0; + nodeColors[vi][1] /= 2; + nodeColors[vi][2] += 0; + nodeColors[vi][2] /= 2; + } + } + std::for_each(nodeColors.begin(), nodeColors.end(), [](std::array &color) { + const double norm = sqrt(std::pow(color[0], 2) + std::pow(color[1], 2) + + std::pow(color[2], 2)); + if (norm > std::pow(10, -7)) { + color[0] /= norm; + color[1] /= norm; + color[2] /= norm; + } + }); - // per node external forces - std::vector> externalForces(pMesh->VN()); - for (const auto &forcePair : nodalExternalForces) { - auto index = forcePair.first; - auto force = forcePair.second; - externalForces[index] = {force[0], force[1], force[2]}; - } + if (!nodeColors.empty()) { + polyscope::getCurveNetwork(meshLabel) + ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors) + ->setEnabled(shouldEnable); + } - if (!externalForces.empty()) { - polyscope::getCurveNetwork(meshLabel) - ->addNodeVectorQuantity("External force_" + label, externalForces) - ->setEnabled(false); - } + // per node external forces + std::vector> externalForces(pMesh->VN()); + for (const auto &forcePair : nodalExternalForces) { + auto index = forcePair.first; + auto force = forcePair.second; + externalForces[index] = {force[0], force[1], force[2]}; + } + + if (!externalForces.empty()) { + polyscope::getCurveNetwork(meshLabel) + ->addNodeVectorQuantity("External force_" + label, externalForces) + ->setEnabled(shouldEnable); + } + // per node external moments + std::vector> externalMoments(pMesh->VN()); + for (const auto &forcePair : nodalExternalForces) { + auto index = forcePair.first; + auto load = forcePair.second; + externalMoments[index] = {load[3], load[4], load[5]}; + } + + if (!externalMoments.empty()) { + polyscope::getCurveNetwork(meshLabel) + ->addNodeVectorQuantity("External moment_" + label, externalMoments) + ->setEnabled(shouldEnable); + } + } + void unregister(const std::string &meshLabel) + { + polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label); + polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label); } #endif // POLYSCOPE_DEFINED }; @@ -329,64 +429,100 @@ void read_binary(const std::string &filename, Matrix &matrix) { } } // namespace Eigen -struct SimulationResults { - bool converged{true}; - std::shared_ptr job; - SimulationHistory history; - std::vector displacements; - std::vector> rotationalDisplacementQuaternion; //per vertex - double internalPotentialEnergy{0}; - double executionTime{0}; - std::string labelPrefix{"deformed"}; - inline static char deliminator{' '}; +struct SimulationResults +{ + /*TODO: remove rotationalDisplacementQuaternion since the last three components of the displacments + * vector contains the same info using euler angles +*/ + bool converged{true}; + std::shared_ptr job; + SimulationHistory history; + std::vector displacements; + std::vector> rotationalDisplacementQuaternion; //per vertex + double internalPotentialEnergy{0}; + double executionTime{0}; + std::string labelPrefix{"deformed"}; + inline static char deliminator{' '}; + SimulationResults() { job = std::make_shared(); } - 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]); - }); + 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; - } - std::string getLabel() const { - return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + - job->getLabel(); - } - - void saveDeformedModel() { - VCGEdgeMesh m; - vcg::tri::Append::MeshCopy(m, *job->pMesh); - - for (int vi = 0; vi < m.VN(); vi++) { - m.vert[vi].P() = - m.vert[vi].P() + CoordType(displacements[vi][0], displacements[vi][1], - displacements[vi][2]); + return translationalDisplacements; } - m.save(getLabel() + ".ply"); - } - void save() { - const std::string filename(getLabel() + "_displacements.eigenBin"); + void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; } + std::string getLabel() const + { + return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel(); + } - Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); - Eigen::write_binary(filename, m); - } + void saveDeformedModel(const std::string &outputFolder = std::string()) + { + VCGEdgeMesh m; + vcg::tri::Append::MeshCopy(m, *job->pMesh); - // The comparison of the results happens comparing the 6-dof nodal - // displacements - bool isEqual(const Eigen::MatrixXd &nodalDisplacements) { - assert(nodalDisplacements.cols() == 6); - Eigen::MatrixXd eigenDisplacements = - Utilities::toEigenMatrix(this->displacements); - const double errorNorm = (eigenDisplacements - nodalDisplacements).norm(); - return errorNorm < 1e-10; - // return eigenDisplacements.isApprox(nodalDisplacements); - } + for (int vi = 0; vi < m.VN(); vi++) { + m.vert[vi].P() = m.vert[vi].P() + + CoordType(displacements[vi][0], + displacements[vi][1], + displacements[vi][2]); + } + m.save(std::filesystem::path(outputFolder).append(getLabel() + ".ply").string()); + } + + 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::string filename(getLabel() + "_displacements.eigenBin"); + + Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); + Eigen::write_binary(std::filesystem::path(outputPath).append(filename).string(), m); + + saveDeformedModel(outputFolder); + } + + // The comparison of the results happens comparing the 6-dof nodal + // displacements + bool isEqual(const Eigen::MatrixXd &nodalDisplacements) + { + assert(nodalDisplacements.cols() == 6); + Eigen::MatrixXd eigenDisplacements = Utilities::toEigenMatrix(this->displacements); + const double errorNorm = (eigenDisplacements - nodalDisplacements).norm(); + return errorNorm < 1e-10; + // return eigenDisplacements.isApprox(nodalDisplacements); + } + + void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom) + { + //load job + job->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json")); + //Use the first .eigenBin file for loading the displacements + for (auto const &entry : std::filesystem::recursive_directory_iterator(loadFromPath)) { + if (filesystem::is_regular_file(entry) && entry.path().extension() == ".eigenBin") { + Eigen::MatrixXd displacements_eigen; + Eigen::read_binary(entry.path().string(), displacements_eigen); + displacements = Utilities::fromEigenMatrix(displacements_eigen); + break; + } + } + + rotationalDisplacementQuaternion.resize(job->pMesh->VN()); + for (int vi = 0; vi < job->pMesh->VN(); vi++) { + rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3], + Eigen::Vector3d::UnitX()) + * Eigen::AngleAxisd(displacements[vi][4], + Eigen::Vector3d::UnitY()) + * Eigen::AngleAxisd(displacements[vi][5], + Eigen::Vector3d::UnitZ()); + } + } #ifdef POLYSCOPE_DEFINED void unregister() const { @@ -396,29 +532,21 @@ struct SimulationResults { std::cerr << "Nothing to remove." << std::endl; return; } + job->unregister(getLabel()); polyscope::removeCurveNetwork(getLabel()); } void registerForDrawing(const std::optional> &desiredColor = std::nullopt, const bool &shouldEnable = true, - const bool &showFrames = false) const + const bool &shouldShowFrames = false) const { - polyscope::options::groundPlaneEnabled = false; - polyscope::view::upDir = polyscope::view::UpDir::ZUp; - if (!polyscope::state::initialized) { - polyscope::init(); - } /* else { - polyscope::removeAllStructures(); - }*/ - // const std::string undeformedMeshName = "Undeformed_" + label; - // polyscope::registerCurveNetwork( - // undeformedMeshName, mesh->getEigenVertices(), - // mesh->getEigenEdges()); - + PolyscopeInterface::init(); const std::shared_ptr &mesh = job->pMesh; - auto polyscopeHandle_deformedEdmeMesh - = polyscope::registerCurveNetwork(getLabel(), - mesh->getEigenVertices(), - mesh->getEigenEdges()); + polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh; + if (!polyscope::hasStructure(getLabel())) { + polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(), + mesh->getEigenVertices(), + mesh->getEigenEdges()); + } polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable); polyscopeHandle_deformedEdmeMesh->setRadius(0.002, true); if (desiredColor.has_value()) { @@ -446,23 +574,42 @@ struct SimulationResults { polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices() + nodalDisplacements); - polyscopeHandle_deformedEdmeMesh->addNodeVectorQuantity("FrameX", framesX) - ->setVectorLengthScale(0.1) - ->setVectorRadius(0.01) - ->setVectorColor(/*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)) - ->setEnabled(showFrames); - polyscopeHandle_deformedEdmeMesh->addNodeVectorQuantity("FrameY", framesY) - ->setVectorLengthScale(0.1) - ->setVectorRadius(0.01) - ->setVectorColor(/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)) - ->setEnabled(showFrames); - polyscopeHandle_deformedEdmeMesh->addNodeVectorQuantity("FrameZ", framesZ) - ->setVectorLengthScale(0.1) - ->setVectorRadius(0.01) - ->setVectorColor(/*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)) - ->setEnabled(showFrames); - + //if(showFramesOn.contains(vi)){ + auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameX", framesX); + polyscopeHandle_frameX->setVectorLengthScale(0.01); + polyscopeHandle_frameX->setVectorRadius(0.01); + polyscopeHandle_frameX->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); + auto polyscopeHandle_frameY = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameY", framesY); + polyscopeHandle_frameY->setVectorLengthScale(0.01); + polyscopeHandle_frameY->setVectorRadius(0.01); + polyscopeHandle_frameY->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); + auto polyscopeHandle_frameZ = polyscopeHandle_deformedEdmeMesh + ->addNodeVectorQuantity("FrameZ", framesZ); + polyscopeHandle_frameZ->setVectorLengthScale(0.01); + polyscopeHandle_frameZ->setVectorRadius(0.01); + polyscopeHandle_frameZ->setVectorColor( + /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); + //} job->registerForDrawing(getLabel()); + static bool wasExecuted = false; + if (!wasExecuted) { + std::function callback = [&]() { + static bool showFrames = shouldShowFrames; + + if (ImGui::Checkbox("Show Frames", &showFrames) && showFrames) { + polyscopeHandle_frameX->setEnabled(showFrames); + polyscopeHandle_frameY->setEnabled(showFrames); + polyscopeHandle_frameZ->setEnabled(showFrames); + } + }; + + PolyscopeInterface::addUserCallback(callback); + wasExecuted = true; + } } #endif }; diff --git a/simulationmesh.cpp b/simulationmesh.cpp index ec3ad3e..4892939 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -343,8 +343,12 @@ bool SimulationMesh::load(const string &plyFilename) { bool SimulationMesh::save(const std::string &plyFilename) { + std::string filename = plyFilename; + if (filename.empty()) { + filename = std::filesystem::current_path().append(getLabel() + ".ply"); + } nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; - customAttrib.GetMeshAttrib(plyFilename); + customAttrib.GetMeshAttrib(filename); std::vector dimensions = getBeamDimensions(); customAttrib.AddEdgeAttribDescriptor(plyPropertyBeamDimensionsID, @@ -359,7 +363,7 @@ bool SimulationMesh::save(const std::string &plyFilename) mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - if (nanoply::NanoPlyWrapper::SaveModel(plyFilename.c_str(), + if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), *this, mask, customAttrib, diff --git a/simulationmesh.hpp b/simulationmesh.hpp index e815e51..95e02b6 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -44,7 +44,7 @@ public: double totalResidualForcesNorm{0}; double currentTotalPotentialEnergykN{0}; double previousTotalPotentialEnergykN{0}; - bool save(const std::string &plyFilename); + bool save(const std::string &plyFilename = std::string()); void setBeamCrossSection(const CrossSectionType &beamDimensions); void setBeamMaterial(const double &pr, const double &ym); void reset(); diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index a8efcec..42f242f 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -6,8 +6,7 @@ #include const bool debugIsOn{false}; -const bool exportArticulationPointsPatterns{false}; -const bool savePlyFiles{false}; +const bool savePlyFiles{true}; // size_t binomialCoefficient(size_t n, size_t m) { // assert(n > m); @@ -34,351 +33,338 @@ const bool savePlyFiles{false}; // .string()); //} -size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const { - if (ni1 <= ni0) { - std::swap(ni0, ni1); - } - assert(ni1 > ni0); - const size_t &n = numberOfNodes; - return (n * (n - 1) / 2) - (n - ni0) * ((n - ni0) - 1) / 2 + ni1 - ni0 - 1; +size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const +{ + if (ni1 <= ni0) { + std::swap(ni0, ni1); + } + assert(ni1 > ni0); + const size_t &n = numberOfNodes; + return (n * (n - 1) / 2) - (n - ni0) * ((n - ni0) - 1) / 2 + ni1 - ni0 - 1; } TopologyEnumerator::TopologyEnumerator() {} -void TopologyEnumerator::computeValidPatterns(const std::vector &reducedNumberOfNodesPerSlot,const std::string& desiredResultsPath, const int &numberOfDesiredEdges) { - assert(reducedNumberOfNodesPerSlot.size() == 5); - assert(reducedNumberOfNodesPerSlot[0] == 0 || - reducedNumberOfNodesPerSlot[0] == 1); - assert(reducedNumberOfNodesPerSlot[1] == 0 || - reducedNumberOfNodesPerSlot[1] == 1); - std::vector numberOfNodesPerSlot{ - reducedNumberOfNodesPerSlot[0], reducedNumberOfNodesPerSlot[1], - reducedNumberOfNodesPerSlot[1], reducedNumberOfNodesPerSlot[2], - reducedNumberOfNodesPerSlot[3], reducedNumberOfNodesPerSlot[2], - reducedNumberOfNodesPerSlot[4]}; - // Generate an edge mesh wih all possible edges - numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.end(), 0); - const size_t numberOfAllPossibleEdges = - numberOfNodes * (numberOfNodes - 1) / 2; +void TopologyEnumerator::computeValidPatterns(const std::vector &reducedNumberOfNodesPerSlot, + const std::string &desiredResultsPath, + const int &numberOfDesiredEdges) +{ + assert(reducedNumberOfNodesPerSlot.size() == 5); + assert(reducedNumberOfNodesPerSlot[0] == 0 || reducedNumberOfNodesPerSlot[0] == 1); + assert(reducedNumberOfNodesPerSlot[1] == 0 || reducedNumberOfNodesPerSlot[1] == 1); + std::vector numberOfNodesPerSlot{reducedNumberOfNodesPerSlot[0], + reducedNumberOfNodesPerSlot[1], + reducedNumberOfNodesPerSlot[1], + reducedNumberOfNodesPerSlot[2], + reducedNumberOfNodesPerSlot[3], + reducedNumberOfNodesPerSlot[2], + reducedNumberOfNodesPerSlot[4]}; + // Generate an edge mesh wih all possible edges + numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(), numberOfNodesPerSlot.end(), 0); + const size_t numberOfAllPossibleEdges = numberOfNodes * (numberOfNodes - 1) / 2; - std::vector allPossibleEdges(numberOfAllPossibleEdges); - const int &n = numberOfNodes; - for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges; - edgeIndex++) { - const int ni0 = - n - 2 - - std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); - const int ni1 = - edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2; - allPossibleEdges[edgeIndex] = vcg::Point2i(ni0, ni1); - } - PatternGeometry patternGeometryAllEdges; - patternGeometryAllEdges.add(numberOfNodesPerSlot, allPossibleEdges); - // Create Results path - auto resultPath=std::filesystem::path(desiredResultsPath); - assert(std::filesystem::exists(resultPath)); + std::vector allPossibleEdges(numberOfAllPossibleEdges); + const int &n = numberOfNodes; + for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges; edgeIndex++) { + const int ni0 = n - 2 + - std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); + const int ni1 = edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2; + allPossibleEdges[edgeIndex] = vcg::Point2i(ni0, ni1); + } + PatternGeometry patternGeometryAllEdges; + patternGeometryAllEdges.add(numberOfNodesPerSlot, allPossibleEdges); + // Create Results path + auto resultPath = std::filesystem::path(desiredResultsPath); + assert(std::filesystem::exists(resultPath)); - auto allResultsPath = resultPath.append("Results"); - std::filesystem::create_directory(allResultsPath); - std::string setupString; - // for (size_t numberOfNodes : reducedNumberOfNodesPerSlot) { - for (size_t numberOfNodesPerSlotIndex = 0; - numberOfNodesPerSlotIndex < reducedNumberOfNodesPerSlot.size(); - numberOfNodesPerSlotIndex++) { - std::string elemID; - if (numberOfNodesPerSlotIndex == 0 || numberOfNodesPerSlotIndex == 1) { - elemID = "v"; - } else if (numberOfNodesPerSlotIndex == 2 || - numberOfNodesPerSlotIndex == 3) { - elemID = "e"; + auto allResultsPath = resultPath.append("Results"); + std::filesystem::create_directory(allResultsPath); + std::string setupString; + // for (size_t numberOfNodes : reducedNumberOfNodesPerSlot) { + for (size_t numberOfNodesPerSlotIndex = 0; + numberOfNodesPerSlotIndex < reducedNumberOfNodesPerSlot.size(); + numberOfNodesPerSlotIndex++) { + std::string elemID; + if (numberOfNodesPerSlotIndex == 0 || numberOfNodesPerSlotIndex == 1) { + elemID = "v"; + } else if (numberOfNodesPerSlotIndex == 2 || numberOfNodesPerSlotIndex == 3) { + elemID = "e"; + } else { + elemID = "c"; + } + setupString += std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex]) + + elemID + "_"; + } + setupString += std::to_string(PatternGeometry().getFanSize()) + "fan"; + if (debugIsOn) { + setupString += "_debug"; + } + auto resultsPath = std::filesystem::path(allResultsPath).append(setupString); + // std::filesystem::remove_all(resultsPath); // delete previous results + std::filesystem::create_directory(resultsPath); + + if (debugIsOn) { + patternGeometryAllEdges.save( + std::filesystem::path(resultsPath).append("allPossibleEdges.ply").string()); + } + // statistics.numberOfPossibleEdges = numberOfAllPossibleEdges; + + std::vector validEdges = getValidEdges(numberOfNodesPerSlot, + resultsPath, + patternGeometryAllEdges, + allPossibleEdges); + PatternGeometry patternAllValidEdges; + patternAllValidEdges.add(patternGeometryAllEdges.getVertices(), validEdges); + if (debugIsOn) { + // Export all valid edges in a ply + patternAllValidEdges.save( + std::filesystem::path(resultsPath).append("allValidEdges.ply").string()); + } + // statistics.numberOfValidEdges = validEdges.size(); + + // Find pairs of intersecting edges + std::unordered_map> intersectingEdges + = patternAllValidEdges.getIntersectingEdges(statistics.numberOfIntersectingEdgePairs); + if (debugIsOn) { + auto intersectingEdgesPath = std::filesystem::path(resultsPath) + .append("All_intersecting_edge_pairs"); + std::filesystem::create_directory(intersectingEdgesPath); + // Export intersecting pairs in ply files + for (auto mapIt = intersectingEdges.begin(); mapIt != intersectingEdges.end(); mapIt++) { + for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end(); setIt++) { + PatternGeometry intersectingEdgePair; + const size_t ei0 = mapIt->first; + const size_t ei1 = *setIt; + vcg::tri::Allocator::AddEdge( + intersectingEdgePair, + patternGeometryAllEdges.getVertices()[validEdges[ei0][0]], + patternGeometryAllEdges.getVertices()[validEdges[ei0][1]]); + vcg::tri::Allocator::AddEdge( + intersectingEdgePair, + patternGeometryAllEdges.getVertices()[validEdges[ei1][0]], + patternGeometryAllEdges.getVertices()[validEdges[ei1][1]]); + intersectingEdgePair.save(std::filesystem::path(intersectingEdgesPath) + .append(std::to_string(mapIt->first) + "_" + + std::to_string(*setIt) + ".ply") + .string()); + } + } + } + + // assert(validEdges.size() == allPossibleEdges.size() - + // coincideEdges.size() - + // duplicateEdges.size()); + + // PatternSet patternSet; + // const std::vector nodes = + // patternGeometryAllEdges.getVertices(); const size_t numberOfNodes = + // nodes.size(); patternSet.nodes.resize(numberOfNodes); for (size_t + // nodeIndex = 0; nodeIndex < numberOfNodes; nodeIndex++) { + // patternSet.nodes[nodeIndex] = + // vcg::Point3d(nodes[nodeIndex][0], nodes[nodeIndex][1],0); + // } + // if (std::filesystem::exists(std::filesystem::path(resultsPath) + // .append("patterns.patt") + // .string())) { + // std::filesystem::remove( + // std::filesystem::path(resultsPath).append("patterns.patt")); + // } + if (numberOfDesiredEdges == -1) { + for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); numberOfEdges++) { + std::cout << "Computing " + setupString << " with " << numberOfEdges << " edges." + << std::endl; + + auto perEdgeResultPath = std::filesystem::path(resultsPath) + .append(std::to_string(numberOfEdges)); + if (std::filesystem::exists(perEdgeResultPath)) { + continue; + } + std::filesystem::create_directory(perEdgeResultPath); + computeValidPatterns(numberOfNodesPerSlot, + numberOfEdges, + perEdgeResultPath, + patternGeometryAllEdges.getVertices(), + intersectingEdges, + validEdges); + // statistics.print(setupString, perEdgeResultPath); + } } else { - elemID = "c"; + std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges << " edges." + << std::endl; + + auto perEdgeResultPath = std::filesystem::path(resultsPath) + .append(std::to_string(numberOfDesiredEdges)); + if (std::filesystem::exists(perEdgeResultPath)) { + // return; + std::filesystem::remove_all(perEdgeResultPath); + } + std::filesystem::create_directory(perEdgeResultPath); + computeValidPatterns(numberOfNodesPerSlot, + numberOfDesiredEdges, + perEdgeResultPath, + patternGeometryAllEdges.getVertices(), + intersectingEdges, + validEdges); } - setupString += - std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex]) + - elemID + "_"; - } - setupString += std::to_string(PatternGeometry().getFanSize()) + "fan"; - if (debugIsOn) { - setupString += "_debug"; - } - auto resultsPath = std::filesystem::path(allResultsPath).append(setupString); - // std::filesystem::remove_all(resultsPath); // delete previous results - std::filesystem::create_directory(resultsPath); - - if (debugIsOn) { - patternGeometryAllEdges.save(std::filesystem::path(resultsPath) - .append("allPossibleEdges.ply") - .string()); - } - // statistics.numberOfPossibleEdges = numberOfAllPossibleEdges; - - std::vector validEdges = - getValidEdges(numberOfNodesPerSlot, resultsPath, patternGeometryAllEdges, - allPossibleEdges); - PatternGeometry patternAllValidEdges; - patternAllValidEdges.add(patternGeometryAllEdges.getVertices(), validEdges); - if (debugIsOn) { - // Export all valid edges in a ply - patternAllValidEdges.save(std::filesystem::path(resultsPath) - .append("allValidEdges.ply") - .string()); - } - // statistics.numberOfValidEdges = validEdges.size(); - - // Find pairs of intersecting edges - std::unordered_map> intersectingEdges = - patternAllValidEdges.getIntersectingEdges( - statistics.numberOfIntersectingEdgePairs); - if (debugIsOn) { - auto intersectingEdgesPath = std::filesystem::path(resultsPath) - .append("All_intersecting_edge_pairs"); - std::filesystem::create_directory(intersectingEdgesPath); - // Export intersecting pairs in ply files - for (auto mapIt = intersectingEdges.begin(); - mapIt != intersectingEdges.end(); mapIt++) { - for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end(); - setIt++) { - PatternGeometry intersectingEdgePair; - const size_t ei0 = mapIt->first; - const size_t ei1 = *setIt; - vcg::tri::Allocator::AddEdge( - intersectingEdgePair, - patternGeometryAllEdges.getVertices()[validEdges[ei0][0]], - patternGeometryAllEdges.getVertices()[validEdges[ei0][1]]); - vcg::tri::Allocator::AddEdge( - intersectingEdgePair, - patternGeometryAllEdges.getVertices()[validEdges[ei1][0]], - patternGeometryAllEdges.getVertices()[validEdges[ei1][1]]); - intersectingEdgePair.save( - std::filesystem::path(intersectingEdgesPath) - .append(std::to_string(mapIt->first) + "_" + - std::to_string(*setIt) + ".ply") - .string()); - } - } - } - - // assert(validEdges.size() == allPossibleEdges.size() - - // coincideEdges.size() - - // duplicateEdges.size()); - - // PatternSet patternSet; - // const std::vector nodes = - // patternGeometryAllEdges.getVertices(); const size_t numberOfNodes = - // nodes.size(); patternSet.nodes.resize(numberOfNodes); for (size_t - // nodeIndex = 0; nodeIndex < numberOfNodes; nodeIndex++) { - // patternSet.nodes[nodeIndex] = - // vcg::Point3d(nodes[nodeIndex][0], nodes[nodeIndex][1],0); - // } - // if (std::filesystem::exists(std::filesystem::path(resultsPath) - // .append("patterns.patt") - // .string())) { - // std::filesystem::remove( - // std::filesystem::path(resultsPath).append("patterns.patt")); - // } - if(numberOfDesiredEdges==-1){ - for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); - numberOfEdges++) { - std::cout << "Computing " + setupString << " with " << numberOfEdges - << " edges." << std::endl; - - auto perEdgeResultPath = std::filesystem::path(resultsPath) - .append(std::to_string(numberOfEdges)); - if (std::filesystem::exists(perEdgeResultPath)) { - continue; - } - std::filesystem::create_directory(perEdgeResultPath); - computeValidPatterns(numberOfNodesPerSlot, numberOfEdges, perEdgeResultPath, - patternGeometryAllEdges.getVertices(), - intersectingEdges, validEdges); - // statistics.print(setupString, perEdgeResultPath); - } - } - else{ - std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges - << " edges." << std::endl; - - auto perEdgeResultPath = std::filesystem::path(resultsPath) - .append(std::to_string(numberOfDesiredEdges)); - if (std::filesystem::exists(perEdgeResultPath)) { - return; - } - std::filesystem::create_directory(perEdgeResultPath); - computeValidPatterns(numberOfNodesPerSlot,numberOfDesiredEdges, perEdgeResultPath, - patternGeometryAllEdges.getVertices(), - intersectingEdges, validEdges); - - } } -void TopologyEnumerator::computeEdgeNodes( - const std::vector &numberOfNodesPerSlot, - std::vector &nodesEdge0, std::vector &nodesEdge1, - std::vector &nodesEdge2) { - // Create vectors holding the node indices of each pattern node of each - // triangle edge - size_t nodeIndex = 0; - if (numberOfNodesPerSlot[0] != 0) { - nodesEdge0.push_back(nodeIndex++); - } - if (numberOfNodesPerSlot[1] != 0) - nodesEdge1.push_back(nodeIndex++); - if (numberOfNodesPerSlot[2] != 0) - nodesEdge2.push_back(nodeIndex++); - - if (numberOfNodesPerSlot[3] != 0) { - for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3]; - edgeNodeIndex++) { - nodesEdge0.push_back(nodeIndex++); +void TopologyEnumerator::computeEdgeNodes(const std::vector &numberOfNodesPerSlot, + std::vector &nodesEdge0, + std::vector &nodesEdge1, + std::vector &nodesEdge2) +{ + // Create vectors holding the node indices of each pattern node of each + // triangle edge + size_t nodeIndex = 0; + if (numberOfNodesPerSlot[0] != 0) { + nodesEdge0.push_back(nodeIndex++); } - } - if (numberOfNodesPerSlot[4] != 0) { - for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4]; - edgeNodeIndex++) { - nodesEdge1.push_back(nodeIndex++); - } - } + if (numberOfNodesPerSlot[1] != 0) + nodesEdge1.push_back(nodeIndex++); + if (numberOfNodesPerSlot[2] != 0) + nodesEdge2.push_back(nodeIndex++); - if (numberOfNodesPerSlot[5] != 0) { - for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5]; - edgeNodeIndex++) { - nodesEdge2.push_back(nodeIndex++); + if (numberOfNodesPerSlot[3] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3]; edgeNodeIndex++) { + nodesEdge0.push_back(nodeIndex++); + } + } + if (numberOfNodesPerSlot[4] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4]; edgeNodeIndex++) { + nodesEdge1.push_back(nodeIndex++); + } } - } - if (numberOfNodesPerSlot[1] != 0) { - assert(numberOfNodesPerSlot[2]); - nodesEdge0.push_back(1); - nodesEdge1.push_back(2); - } - if (numberOfNodesPerSlot[0] != 0) { - nodesEdge2.push_back(0); - } + if (numberOfNodesPerSlot[5] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5]; edgeNodeIndex++) { + nodesEdge2.push_back(nodeIndex++); + } + } + if (numberOfNodesPerSlot[1] != 0) { + assert(numberOfNodesPerSlot[2]); + nodesEdge0.push_back(1); + nodesEdge1.push_back(2); + } + + if (numberOfNodesPerSlot[0] != 0) { + nodesEdge2.push_back(0); + } } std::unordered_set TopologyEnumerator::computeCoincideEdges( - const std::vector &numberOfNodesPerSlot) { - /* + const std::vector &numberOfNodesPerSlot) +{ + /* * A coincide edge is defined as an edge connection between two nodes that lay * on a triangle edge and which have another node in between * */ - std::vector nodesEdge0; // left edge - std::vector nodesEdge1; // bottom edge - std::vector nodesEdge2; // right edge - computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); + std::vector nodesEdge0; // left edge + std::vector nodesEdge1; // bottom edge + std::vector nodesEdge2; // right edge + computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); - std::vector coincideEdges0 = getCoincideEdges(nodesEdge0); - std::vector coincideEdges1 = getCoincideEdges(nodesEdge1); - std::vector coincideEdges2 = getCoincideEdges(nodesEdge2); - std::unordered_set coincideEdges{coincideEdges0.begin(), - coincideEdges0.end()}; - std::copy(coincideEdges1.begin(), coincideEdges1.end(), - std::inserter(coincideEdges, coincideEdges.end())); - std::copy(coincideEdges2.begin(), coincideEdges2.end(), - std::inserter(coincideEdges, coincideEdges.end())); + std::vector coincideEdges0 = getCoincideEdges(nodesEdge0); + std::vector coincideEdges1 = getCoincideEdges(nodesEdge1); + std::vector coincideEdges2 = getCoincideEdges(nodesEdge2); + std::unordered_set coincideEdges{coincideEdges0.begin(), coincideEdges0.end()}; + std::copy(coincideEdges1.begin(), + coincideEdges1.end(), + std::inserter(coincideEdges, coincideEdges.end())); + std::copy(coincideEdges2.begin(), + coincideEdges2.end(), + std::inserter(coincideEdges, coincideEdges.end())); - if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[1]) { - coincideEdges.insert(getEdgeIndex(0, 2)); - } + if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[1]) { + coincideEdges.insert(getEdgeIndex(0, 2)); + } - if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[2]) { - assert(numberOfNodesPerSlot[1]); - coincideEdges.insert(getEdgeIndex(0, 3)); - } + if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[2]) { + assert(numberOfNodesPerSlot[1]); + coincideEdges.insert(getEdgeIndex(0, 3)); + } - return coincideEdges; + return coincideEdges; } std::unordered_set TopologyEnumerator::computeDuplicateEdges( - const std::vector &numberOfNodesPerSlot) { - /* + const std::vector &numberOfNodesPerSlot) +{ + /* * A duplicate edges are all edges the "right" edge since due to rotational * symmetry "left" edge=="right" edge * */ - std::unordered_set duplicateEdges; - std::vector nodesEdge0; // left edge - std::vector nodesEdge1; // bottom edge - std::vector nodesEdge2; // right edge - computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); - if (numberOfNodesPerSlot[5]) { - for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1; - edge2NodeIndex++) { - const size_t nodeIndex = nodesEdge2[edge2NodeIndex]; - const size_t nextNodeIndex = nodesEdge2[edge2NodeIndex + 1]; - duplicateEdges.insert(getEdgeIndex(nodeIndex, nextNodeIndex)); + std::unordered_set duplicateEdges; + std::vector nodesEdge0; // left edge + std::vector nodesEdge1; // bottom edge + std::vector nodesEdge2; // right edge + computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); + if (numberOfNodesPerSlot[5]) { + for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1; edge2NodeIndex++) { + const size_t nodeIndex = nodesEdge2[edge2NodeIndex]; + const size_t nextNodeIndex = nodesEdge2[edge2NodeIndex + 1]; + duplicateEdges.insert(getEdgeIndex(nodeIndex, nextNodeIndex)); + } } - } - return duplicateEdges; + return duplicateEdges; } std::vector TopologyEnumerator::getValidEdges( const std::vector &numberOfNodesPerSlot, const std::filesystem::path &resultsPath, const PatternGeometry &patternGeometryAllEdges, - const std::vector &allPossibleEdges) { - - std::unordered_set coincideEdges = - computeCoincideEdges(numberOfNodesPerSlot); - // Export each coincide edge into a ply file - if (!coincideEdges.empty() && debugIsOn) { - auto coincideEdgesPath = - std::filesystem::path(resultsPath).append("Coincide_edges"); - std::filesystem::create_directories(coincideEdgesPath); - for (auto coincideEdgeIndex : coincideEdges) { - PatternGeometry::EdgeType e = - patternGeometryAllEdges.edge[coincideEdgeIndex]; - PatternGeometry singleEdgeMesh; - vcg::Point3d p0 = e.cP(0); - vcg::Point3d p1 = e.cP(1); - std::vector edgeVertices; - edgeVertices.push_back(p0); - edgeVertices.push_back(p1); - singleEdgeMesh.add(edgeVertices); - singleEdgeMesh.add(std::vector{vcg::Point2i{0, 1}}); - singleEdgeMesh.save(std::filesystem::path(coincideEdgesPath) - .append(std::to_string(coincideEdgeIndex)) - .string() + - ".ply"); + const std::vector &allPossibleEdges) +{ + std::unordered_set coincideEdges = computeCoincideEdges(numberOfNodesPerSlot); + // Export each coincide edge into a ply file + if (!coincideEdges.empty() && debugIsOn) { + auto coincideEdgesPath = std::filesystem::path(resultsPath).append("Coincide_edges"); + std::filesystem::create_directories(coincideEdgesPath); + for (auto coincideEdgeIndex : coincideEdges) { + PatternGeometry::EdgeType e = patternGeometryAllEdges.edge[coincideEdgeIndex]; + PatternGeometry singleEdgeMesh; + vcg::Point3d p0 = e.cP(0); + vcg::Point3d p1 = e.cP(1); + std::vector edgeVertices; + edgeVertices.push_back(p0); + edgeVertices.push_back(p1); + singleEdgeMesh.add(edgeVertices); + singleEdgeMesh.add(std::vector{vcg::Point2i{0, 1}}); + singleEdgeMesh.save(std::filesystem::path(coincideEdgesPath) + .append(std::to_string(coincideEdgeIndex)) + .string() + + ".ply"); + } } - } - statistics.numberOfCoincideEdges = coincideEdges.size(); + statistics.numberOfCoincideEdges = coincideEdges.size(); - // Compute duplicate edges - std::unordered_set duplicateEdges = - computeDuplicateEdges(numberOfNodesPerSlot); - if (!duplicateEdges.empty() && debugIsOn) { - // Export duplicate edges in a single ply file - auto duplicateEdgesPath = - std::filesystem::path(resultsPath).append("duplicate"); - std::filesystem::create_directory(duplicateEdgesPath); - PatternGeometry patternDuplicateEdges; - for (auto duplicateEdgeIndex : duplicateEdges) { - PatternGeometry::EdgeType e = - patternGeometryAllEdges.edge[duplicateEdgeIndex]; - vcg::Point3d p0 = e.cP(0); - vcg::Point3d p1 = e.cP(1); - vcg::tri::Allocator::AddEdge(patternDuplicateEdges, p0, - p1); + // Compute duplicate edges + std::unordered_set duplicateEdges = computeDuplicateEdges(numberOfNodesPerSlot); + if (!duplicateEdges.empty() && debugIsOn) { + // Export duplicate edges in a single ply file + auto duplicateEdgesPath = std::filesystem::path(resultsPath).append("duplicate"); + std::filesystem::create_directory(duplicateEdgesPath); + PatternGeometry patternDuplicateEdges; + for (auto duplicateEdgeIndex : duplicateEdges) { + PatternGeometry::EdgeType e = patternGeometryAllEdges.edge[duplicateEdgeIndex]; + vcg::Point3d p0 = e.cP(0); + vcg::Point3d p1 = e.cP(1); + vcg::tri::Allocator::AddEdge(patternDuplicateEdges, p0, p1); + } + patternDuplicateEdges.save( + std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply").string()); } - patternDuplicateEdges.save(std::filesystem::path(duplicateEdgesPath) - .append("duplicateEdges.ply") - .string()); - } - statistics.numberOfDuplicateEdges = duplicateEdges.size(); + statistics.numberOfDuplicateEdges = duplicateEdges.size(); - // Create the set of all possible edges without coincide and duplicate edges - std::vector validEdges; - for (size_t edgeIndex = 0; edgeIndex < allPossibleEdges.size(); edgeIndex++) { - if (coincideEdges.count(edgeIndex) == 0 && - duplicateEdges.count(edgeIndex) == 0) { - validEdges.push_back(allPossibleEdges[edgeIndex]); + // Create the set of all possible edges without coincide and duplicate edges + std::vector validEdges; + for (size_t edgeIndex = 0; edgeIndex < allPossibleEdges.size(); edgeIndex++) { + if (coincideEdges.count(edgeIndex) == 0 && duplicateEdges.count(edgeIndex) == 0) { + validEdges.push_back(allPossibleEdges[edgeIndex]); + } } - } - return validEdges; + return validEdges; } void TopologyEnumerator::computeValidPatterns( @@ -386,252 +372,302 @@ void TopologyEnumerator::computeValidPatterns( const size_t &numberOfDesiredEdges, const std::filesystem::path &resultsPath, const std::vector &allVertices, - const std::unordered_map> - &intersectingEdges, - const std::vector &validEdges) { - assert(numberOfNodesPerSlot.size() == 7); - // Iterate over all patterns which have numberOfDesiredEdges edges from - // from the validEdges Identify patterns that contain dangling edges - const bool enoughValidEdgesExist = validEdges.size() >= numberOfDesiredEdges; - if (!enoughValidEdgesExist) { - std::filesystem::remove_all(resultsPath); // delete previous results folder - return; - } - assert(enoughValidEdgesExist); + const std::unordered_map> &intersectingEdges, + const std::vector &validEdges) +{ + assert(numberOfNodesPerSlot.size() == 7); + // Iterate over all patterns which have numberOfDesiredEdges edges from + // from the validEdges Identify patterns that contain dangling edges + const bool enoughValidEdgesExist = validEdges.size() >= numberOfDesiredEdges; + if (!enoughValidEdgesExist) { + std::filesystem::remove_all(resultsPath); // delete previous results folder + return; + } + assert(enoughValidEdgesExist); - // Create pattern result paths - auto validPatternsPath = std::filesystem::path(resultsPath).append("Valid"); - std::filesystem::create_directory(validPatternsPath); - // std::ofstream validPatternsFileStream; - // validPatternsFileStream.open( - // validPatternsPath.append("patterns.patt").string()); - const std::string validPatternsFilePath = - validPatternsPath.append("patterns.patt").string(); - PatternIO::PatternSet patternSet; - patternSet.nodes = allVertices; - const int patternSetBufferSize = 10000; + // Create pattern result paths + const auto validPatternsPath = std::filesystem::path(resultsPath).append("Valid"); + const bool validPathCreatedSuccesfully = std::filesystem::create_directories(validPatternsPath); + assert(validPathCreatedSuccesfully && std::filesystem::exists(validPatternsPath)); + // std::ofstream validPatternsFileStream; + // validPatternsFileStream.open( + // validPatternsPath.append("patterns.patt").string()); + const std::string compressedPatternsFilePath + = std::filesystem::path(validPatternsPath).append("patterns.patt").string(); + PatternIO::PatternSet patternSet; + patternSet.nodes = allVertices; + const int patternSetBufferSize = 10000; - const size_t numberOfPatterns = PatternGeometry::binomialCoefficient( - validEdges.size(), numberOfDesiredEdges); - statistics.numberOfPatterns = numberOfPatterns; + const size_t numberOfPatterns = PatternGeometry::binomialCoefficient(validEdges.size(), + numberOfDesiredEdges); + statistics.numberOfPatterns = numberOfPatterns; - // Initialize pattern binary representation - std::string patternBinaryRepresentation; - patternBinaryRepresentation = std::string(numberOfDesiredEdges, '1'); - patternBinaryRepresentation += - std::string(validEdges.size() - numberOfDesiredEdges, '0'); - std::sort(patternBinaryRepresentation.begin(), - patternBinaryRepresentation.end()); - /*TODO: Performance could be improved by changing the patternGeometry with + // Initialize pattern binary representation + std::string patternBinaryRepresentation; + patternBinaryRepresentation = std::string(numberOfDesiredEdges, '1'); + patternBinaryRepresentation += std::string(validEdges.size() - numberOfDesiredEdges, '0'); + std::sort(patternBinaryRepresentation.begin(), patternBinaryRepresentation.end()); + /*TODO: Performance could be improved by changing the patternGeometry with * respect to the previous one. Maybe I could xor the binaryRepresentation * to the previous one.*/ - // std::string previousPatternBinaryRepresentation(validEdges.size(),'0'); - size_t patternIndex = 0; - do { - patternIndex++; - const std::string patternName = std::to_string(patternIndex); - // std::cout << "Pattern name:" + patternBinaryRepresentation << - // std::endl; isValidPattern(patternBinaryRepresentation, validEdges, - // numberOfDesiredEdges); - // Create the geometry of the pattern - // Compute the pattern edges from the binary representation - std::vector patternEdges(numberOfDesiredEdges); - size_t patternEdgeIndex = 0; - for (size_t validEdgeIndex = 0; - validEdgeIndex < patternBinaryRepresentation.size(); - validEdgeIndex++) { - if (patternBinaryRepresentation[validEdgeIndex] == '1') { - assert(patternEdgeIndex < numberOfDesiredEdges); - patternEdges[patternEdgeIndex++] = validEdges[validEdgeIndex]; - } - } - - PatternGeometry patternGeometry; - patternGeometry.add(allVertices, patternEdges); - - // Check if pattern contains intersecting edges - const bool patternContainsIntersectingEdges = - patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, - intersectingEdges); - // Export the tiled ply file if it contains intersecting edges - if (patternContainsIntersectingEdges) { - // create the tiled geometry of the pattern - statistics.numberOfPatternsWithIntersectingEdges++; - if (debugIsOn) { - if (savePlyFiles) { - PatternGeometry tiledPatternGeometry = - PatternGeometry::createTile(patternGeometry); - auto intersectingPatternsPath = - std::filesystem::path(resultsPath).append("Intersecting"); - std::filesystem::create_directory(intersectingPatternsPath); - patternGeometry.save( - std::filesystem::path(intersectingPatternsPath) - .append(patternName) - .string() + - ".ply"); - tiledPatternGeometry.save( - std::filesystem::path(intersectingPatternsPath) - .append(patternName + "_tiled") - .string() + - ".ply"); + // std::string previousPatternBinaryRepresentation(validEdges.size(),'0'); + size_t patternIndex = 0; + bool validPatternsExist = false; + do { + patternIndex++; + const std::string patternName = std::to_string(patternIndex); + // std::cout << "Pattern name:" + patternBinaryRepresentation << + // std::endl; isValidPattern(patternBinaryRepresentation, validEdges, + // numberOfDesiredEdges); + // Create the geometry of the pattern + // Compute the pattern edges from the binary representation + std::vector patternEdges(numberOfDesiredEdges); + size_t patternEdgeIndex = 0; + for (size_t validEdgeIndex = 0; validEdgeIndex < patternBinaryRepresentation.size(); + validEdgeIndex++) { + if (patternBinaryRepresentation[validEdgeIndex] == '1') { + assert(patternEdgeIndex < numberOfDesiredEdges); + patternEdges[patternEdgeIndex++] = validEdges[validEdgeIndex]; + } } - } else { - continue; // should be uncommented in order to improve performance - } - } - // Compute the tiled valence - const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges( - numberOfNodesPerSlot); // marks the nodes with valence>=1 - // Create the tiled geometry of the pattern - const bool hasFloatingComponents = - !patternGeometry.isFullyConnectedWhenTiled(); - FlatPatternTopology topology(numberOfNodesPerSlot, patternEdges); - const bool hasArticulationPoints = topology.containsArticulationPoints(); - // duplicated here - // Check dangling edges with vcg method - // const bool vcg_tiledPatternHasDangling = - // tiledPatternGeometry.hasUntiledDanglingEdges(); - if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents && + PatternGeometry patternGeometry; + patternGeometry.add(allVertices, patternEdges); + + // Check if pattern contains intersecting edges + const bool patternContainsIntersectingEdges + = patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, intersectingEdges); + // Export the tiled ply file if it contains intersecting edges + if (patternContainsIntersectingEdges) { + // create the tiled geometry of the pattern + statistics.numberOfPatternsWithIntersectingEdges++; + if (debugIsOn) { + if (savePlyFiles) { + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); + auto intersectingPatternsPath = std::filesystem::path(resultsPath) + .append("Intersecting"); + std::filesystem::create_directory(intersectingPatternsPath); + patternGeometry.save( + std::filesystem::path(intersectingPatternsPath).append(patternName).string() + + ".ply"); + tiledPatternGeometry.save(std::filesystem::path(intersectingPatternsPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + } else { + continue; // should be uncommented in order to improve performance + } + } + + const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold + = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 5); + if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { + if (debugIsOn || savePlyFiles) { + if (savePlyFiles) { + auto danglingEdgesPath = std::filesystem::path(resultsPath) + .append("ExceedingAngleThreshold"); + std::filesystem::create_directory(danglingEdgesPath); + patternGeometry.save( + std::filesystem::path(danglingEdgesPath).append(patternName).string() + + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + } else { + continue; + } + } + + const bool tiledPatternHasNodeWithValenceGreaterThanDesired + = patternGeometry.hasValenceGreaterThan(numberOfNodesPerSlot, 8); + if (tiledPatternHasNodeWithValenceGreaterThanDesired) { + if (debugIsOn) { + if (savePlyFiles) { + auto danglingEdgesPath = std::filesystem::path(resultsPath) + .append("HighValencePatterns"); + std::filesystem::create_directory(danglingEdgesPath); + patternGeometry.save( + std::filesystem::path(danglingEdgesPath).append(patternName).string() + + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + } else { + continue; + } + } + + // Compute the tiled valence + const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges( + numberOfNodesPerSlot); // marks the nodes with valence>=1 + // Create the tiled geometry of the pattern + const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenTiled(); + FlatPatternTopology topology(numberOfNodesPerSlot, patternEdges); + const bool hasArticulationPoints = topology.containsArticulationPoints(); + // duplicated here + // Check dangling edges with vcg method + // const bool vcg_tiledPatternHasDangling = + // tiledPatternGeometry.hasUntiledDanglingEdges(); + if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents && !hasArticulationPoints*/) { - statistics.numberOfPatternsWithADanglingEdgeOrNode++; - if (debugIsOn) { - if (savePlyFiles) { - auto danglingEdgesPath = - std::filesystem::path(resultsPath).append("Dangling"); - std::filesystem::create_directory(danglingEdgesPath); - patternGeometry.save(std::filesystem::path(danglingEdgesPath) - .append(patternName) - .string() + - ".ply"); - PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( - patternGeometry); // the marked nodes of hasDanglingEdges are - tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath) - .append(patternName + "_tiled") - .string() + - ".ply"); + statistics.numberOfPatternsWithADanglingEdgeOrNode++; + if (debugIsOn) { + if (savePlyFiles) { + auto danglingEdgesPath = std::filesystem::path(resultsPath).append("Dangling"); + std::filesystem::create_directory(danglingEdgesPath); + patternGeometry.save( + std::filesystem::path(danglingEdgesPath).append(patternName).string() + + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save(std::filesystem::path(danglingEdgesPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + } else { + continue; + } } - } else { - continue; - } - } - if (hasFloatingComponents /*&& !hasArticulationPoints && + if (hasFloatingComponents /*&& !hasArticulationPoints && !tiledPatternHasDanglingEdges*/) { - statistics.numberOfPatternsWithMoreThanASingleCC++; - if (debugIsOn) { - if (savePlyFiles) { - auto moreThanOneCCPath = - std::filesystem::path(resultsPath).append("MoreThanOneCC"); - std::filesystem::create_directory(moreThanOneCCPath); - patternGeometry.save(std::filesystem::path(moreThanOneCCPath) - .append(patternName) - .string() + - ".ply"); - PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( - patternGeometry); // the marked nodes of hasDanglingEdges are - tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath) - .append(patternName + "_tiled") - .string() + - ".ply"); + statistics.numberOfPatternsWithMoreThanASingleCC++; + if (debugIsOn) { + if (savePlyFiles) { + auto moreThanOneCCPath = std::filesystem::path(resultsPath) + .append("MoreThanOneCC"); + std::filesystem::create_directory(moreThanOneCCPath); + patternGeometry.save( + std::filesystem::path(moreThanOneCCPath).append(patternName).string() + + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + } else { + continue; + } } - } else { - continue; - } - } - if (hasArticulationPoints /*&& !hasFloatingComponents && + if (hasArticulationPoints /*&& !hasFloatingComponents && !tiledPatternHasDanglingEdges*/) { - statistics.numberOfPatternsWithArticulationPoints++; - if (exportArticulationPointsPatterns || debugIsOn) { - if (savePlyFiles) { - auto articulationPointsPath = - std::filesystem::path(resultsPath).append("ArticulationPoints"); - std::filesystem::create_directory(articulationPointsPath); - patternGeometry.save(std::filesystem::path(articulationPointsPath) - .append(patternName) - .string() + - ".ply"); - PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( - patternGeometry); // the marked nodes of hasDanglingEdges are - tiledPatternGeometry.save( - std::filesystem::path(articulationPointsPath) - .append(patternName + "_tiled") - .string() + - ".ply"); + statistics.numberOfPatternsWithArticulationPoints++; + if (debugIsOn) { + if (savePlyFiles) { + auto articulationPointsPath = std::filesystem::path(resultsPath) + .append("ArticulationPoints"); + std::filesystem::create_directory(articulationPointsPath); + patternGeometry.save( + std::filesystem::path(articulationPointsPath).append(patternName).string() + + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save(std::filesystem::path(articulationPointsPath) + .append(patternName + "_tiled") + .string() + + ".ply"); - // std::cout << "Pattern:" << patternName << std::endl; + // std::cout << "Pattern:" << patternName << std::endl; + } + } else { + continue; + } } - } else { - continue; - } + + const bool isValidPattern = !patternContainsIntersectingEdges + && !tiledPatternHasDanglingEdges && !hasFloatingComponents + && !hasArticulationPoints + && !tiledPatternHasNodeWithValenceGreaterThanDesired + && !tiledPatternHasEdgesWithAngleSmallerThanThreshold; + if (isValidPattern) { + // if(patternName=='2055'){ + // PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + // patternGeometry); // the marked nodes of hasDanglingEdges are + // tiledPatternGeometry.registerForDrawing(); + // polyscope::show(); + // tiledPatternGeometry.unregister(); + // } + statistics.numberOfValidPatterns++; + validPatternsExist = true; + if (savePlyFiles) { + // if (numberOfDesiredEdges == 4) { + // std::cout << "Saving:" + // << std::filesystem::path(validPatternsPath) + // .append(patternName) + // .string() + + // ".ply" + // << std::endl; + // } + patternGeometry.save( + std::filesystem::path(validPatternsPath).append(patternName).string() + ".ply"); + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + tiledPatternGeometry.save( + std::filesystem::path(validPatternsPath).append(patternName + "_tiled").string() + + ".ply"); + } + PatternIO::Pattern pattern; + pattern.edges = patternEdges; + pattern.name = patternIndex; + patternSet.patterns.emplace_back(pattern); + // Save valid patterns + // if (patternIndex% patternSetBufferSize == 0) { + if (statistics.numberOfValidPatterns % patternSetBufferSize == 0) { + PatternIO::save(compressedPatternsFilePath, patternSet); + patternSet.patterns.clear(); + patternSet.patterns.reserve(patternSetBufferSize); + } + } + + // assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges); + } while (std::next_permutation(patternBinaryRepresentation.begin(), + patternBinaryRepresentation.end())); + if (!patternSet.patterns.empty()) { + PatternIO::save(compressedPatternsFilePath, patternSet); } - const bool isValidPattern = - !patternContainsIntersectingEdges && !tiledPatternHasDanglingEdges && - !hasFloatingComponents && !hasArticulationPoints; - if (isValidPattern) { - statistics.numberOfValidPatterns++; - if (savePlyFiles) { - // if (numberOfDesiredEdges == 4) { - // std::cout << "Saving:" - // << std::filesystem::path(validPatternsPath) - // .append(patternName) - // .string() + - // ".ply" - // << std::endl; - // } - patternGeometry.save(std::filesystem::path(validPatternsPath) - .append(patternName) - .string() + - ".ply"); - PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( - patternGeometry); // the marked nodes of hasDanglingEdges are - tiledPatternGeometry.save(std::filesystem::path(validPatternsPath) - .append(patternName + "_tiled") - .string() + - ".ply"); - } - PatternIO::Pattern pattern; - pattern.edges = patternEdges; - pattern.name = patternIndex; - patternSet.patterns.emplace_back(pattern); - // Save valid patterns - if (patternIndex % patternSetBufferSize == 0) { - PatternIO::save(validPatternsFilePath, patternSet); - patternSet.patterns.clear(); - patternSet.patterns.reserve(patternSetBufferSize); - } + if (!validPatternsExist) { + std::filesystem::remove_all(validPatternsPath); + if (!debugIsOn) { + std::filesystem::remove_all(resultsPath); + } } - - // assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges); - } while (std::next_permutation(patternBinaryRepresentation.begin(), - patternBinaryRepresentation.end())); - if (!patternSet.patterns.empty()) { - PatternIO::save(validPatternsFilePath, patternSet); - } } - std::vector TopologyEnumerator::getCoincideEdges( - const std::vector &edgeNodeIndices) const { - std::vector coincideEdges; - if (edgeNodeIndices.size() < 3) - return coincideEdges; - for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2; - edgeNodeIndex++) { - const size_t &firstNodeIndex = edgeNodeIndices[edgeNodeIndex]; - for (size_t secondEdgeNodeIndex = edgeNodeIndex + 2; - secondEdgeNodeIndex < edgeNodeIndices.size(); secondEdgeNodeIndex++) { - const size_t &secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex]; - coincideEdges.push_back(getEdgeIndex(firstNodeIndex, secondNodeIndex)); + const std::vector &edgeNodeIndices) const +{ + std::vector coincideEdges; + if (edgeNodeIndices.size() < 3) + return coincideEdges; + for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2; edgeNodeIndex++) { + const size_t &firstNodeIndex = edgeNodeIndices[edgeNodeIndex]; + for (size_t secondEdgeNodeIndex = edgeNodeIndex + 2; + secondEdgeNodeIndex < edgeNodeIndices.size(); + secondEdgeNodeIndex++) { + const size_t &secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex]; + coincideEdges.push_back(getEdgeIndex(firstNodeIndex, secondNodeIndex)); + } } - } - return coincideEdges; + return coincideEdges; } -bool TopologyEnumerator::isValidPattern( - const std::string &patternBinaryRepresentation, - const std::vector &validEdges, - const size_t &numberOfDesiredEdges) const { - return true; +bool TopologyEnumerator::isValidPattern(const std::string &patternBinaryRepresentation, + const std::vector &validEdges, + const size_t &numberOfDesiredEdges) const +{ + return true; } diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index 564ca84..9eb7e61 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -7,6 +7,8 @@ #include #include #include +/* include the support for half edges */ +#include size_t PatternGeometry::computeTiledValence( const size_t &nodeIndex, @@ -14,8 +16,8 @@ size_t PatternGeometry::computeTiledValence( std::vector connectedEdges; vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges); const size_t nodeValence = connectedEdges.size(); - assert(nodeSlot.count(nodeIndex) != 0); - const size_t nodeSlotIndex = nodeSlot.at(nodeIndex); + assert(nodeToSlotMap.count(nodeIndex) != 0); + const size_t nodeSlotIndex = nodeToSlotMap.at(nodeIndex); if (nodeSlotIndex == 0) { return nodeValence * fanSize; } else if (nodeSlotIndex == 1 || nodeSlotIndex == 2) { @@ -136,23 +138,23 @@ PatternGeometry PatternGeometry::createFan(PatternGeometry &pattern) { PatternGeometry::PatternGeometry(PatternGeometry &other) { vcg::tri::Append::MeshCopy(*this, other); this->vertices = other.getVertices(); + baseTriangle = other.getBaseTriangle(); + baseTriangleHeight = computeBaseTriangleHeight(); vcg::tri::UpdateTopology::VertexEdge(*this); vcg::tri::UpdateTopology::EdgeEdge(*this); } -bool PatternGeometry::save(const std::string plyFilename) { +bool PatternGeometry::load(const string &plyFilename) +{ + if (!VCGEdgeMesh::load(plyFilename)) { + return false; + } + vcg::tri::UpdateTopology::VertexEdge(*this); + baseTriangleHeight = computeBaseTriangleHeight(); + baseTriangle = computeBaseTriangle(); - int returnValue = vcg::tri::io::ExporterPLY::Save( - *this, plyFilename.c_str(), - vcg::tri::io::Mask::IOM_EDGEINDEX | vcg::tri::io::Mask::IOM_VERTCOLOR, - false); - if (returnValue != 0) { - std::cerr << vcg::tri::io::ExporterPLY::ErrorMsg( - returnValue) - << std::endl; - } - - return static_cast(returnValue); + updateEigenEdgeAndVertices(); + return true; } void PatternGeometry::add(const std::vector &vertices) { @@ -296,17 +298,199 @@ void PatternGeometry::constructNodeToTiledValenceMap( } } +std::vector PatternGeometry::getEdgeVectorsWithVertexAsOrigin( + std::vector &edgePointers, const int &vi) +{ + std::vector incidentElementsVectors(edgePointers.size()); + for (int incidentElementsIndex = 0; incidentElementsIndex < edgePointers.size(); + incidentElementsIndex++) { + assert(vi == getIndex(edgePointers[incidentElementsIndex]->cV(0)) + || vi == getIndex(edgePointers[incidentElementsIndex]->cV(1))); + + incidentElementsVectors[incidentElementsIndex] + = vi == getIndex(edgePointers[incidentElementsIndex]->cV(0)) + ? edgePointers[incidentElementsIndex]->cP(1) + - edgePointers[incidentElementsIndex]->cP(0) + : edgePointers[incidentElementsIndex]->cP(0) + - edgePointers[incidentElementsIndex]->cP(1); + } + + return incidentElementsVectors; +} + +bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector &numberOfNodesPerSlot, + const double &angleThreshold_degrees) +{ + assert(numberOfNodesPerSlot.size() == 7); + //Initialize helping structs + if (nodeToSlotMap.empty()) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap); + } + if (correspondingNode.empty()) { + constructCorrespondingNodeMap(numberOfNodesPerSlot); + } + + // const double theta0 = vcg::math::ToDeg(vcg::Point2d(1, 0).Angle()); + // const double theta1 = vcg::math::ToDeg(vcg::Point2d(1, 1).Angle()); + // const double theta2 = vcg::math::ToDeg(vcg::Point2d(0, 1).Angle()); + // const double theta3 = vcg::math::ToDeg(vcg::Point2d(-1, 1).Angle()); + // const double theta4 = vcg::math::ToDeg(vcg::Point2d(-1, 0).Angle()); + // const double theta5 = vcg::math::ToDeg(vcg::Point2d(-1, -1).Angle() + 2 * M_PI); + // const double theta6 = vcg::math::ToDeg(vcg::Point2d(0, -1).Angle() + 2 * M_PI); + // const double theta7 = vcg::math::ToDeg(vcg::Point2d(1, -1).Angle() + 2 * M_PI); + + // std::set test{theta0, theta1, theta2, theta3, theta4, theta5, theta6, theta7}; + + bool hasAngleSmallerThan = false; + //find tiled incident edges for each node + using EdgeIndex = int; + using ThetaWithXAxis = double; + for (size_t vi = 0; vi < VN(); vi++) { + std::vector incidentElementPointers; + vcg::edge::VEStarVE(&vert[vi], incidentElementPointers); + const size_t numberOfIncidentEdges = incidentElementPointers.size(); + if (numberOfIncidentEdges == 0) { + continue; + } + std::vector incidentEdgeVectors + = getEdgeVectorsWithVertexAsOrigin(incidentElementPointers, vi); + + std::vector tiledIncidentVectors = incidentEdgeVectors; + const size_t slotIndex = nodeToSlotMap[vi]; + if (slotIndex == 2) { + //NOTE:I assume that base triangle slots 1,2(bottom triangle nodes) are not used + std::cerr + << "Slot of bottom base triangle nodes detected.This case is not handled.Exiting.." + << std::endl; + std::terminate(); + } else if (slotIndex == 3 || slotIndex == 5) { + //NOTE: I don't need to check both triangle edges + std::vector correspondingVertexIncidentElementPointers; + vcg::edge::VEStarVE(&vert[correspondingNode[vi]], + correspondingVertexIncidentElementPointers); + std::vector correspondingVertexIncidentVectors + = getEdgeVectorsWithVertexAsOrigin(correspondingVertexIncidentElementPointers, + correspondingNode[vi]); + // const CoordType &correspondingVertexPosition = vert[correspondingNode[vi]].cP(); + vcg::Matrix33d R; + if (slotIndex == 3) { + R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(-60.0)); + } else { + R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(60.0)); + } + // const CoordType &correspondingVertexPosition_rotated = vert[vi].cP(); + std::transform( + correspondingVertexIncidentVectors.begin(), + correspondingVertexIncidentVectors.end(), + correspondingVertexIncidentVectors.begin(), + [&](const VectorType &v) { + // CoordType rotatedTarget = v * R; + // v = rotatedTarget - correspondingVertexPosition_rotated; + return R * v; + }); + tiledIncidentVectors.insert(tiledIncidentVectors.end(), + correspondingVertexIncidentVectors.begin(), + correspondingVertexIncidentVectors.end()); + } else if (slotIndex == 4) { + std::vector reversedIncidentElementVectors(incidentEdgeVectors.size()); + std::transform(incidentEdgeVectors.begin(), + incidentEdgeVectors.end(), + reversedIncidentElementVectors.begin(), + [&](const VectorType &v) { return -v; }); + //NOTE: for checking the angles not all opposite incident element vectors are needed but rather the two "extreme" ones of slot 4 + //here I simply add them all + tiledIncidentVectors.insert(tiledIncidentVectors.end(), + reversedIncidentElementVectors.begin(), + reversedIncidentElementVectors.end()); + } + // std::vector> edgePoints; + // for (int tiledVectorIndex = 0; tiledVectorIndex < tiledIncidentVectors.size(); + // tiledVectorIndex++) { + // edgePoints.push_back( + // std::array{vert[vi].cP()[0], vert[vi].cP()[1], vert[vi].cP()[2]}); + // edgePoints.push_back( + // std::array{(vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[0], + // (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[1], + // (vert[vi].cP() + tiledIncidentVectors[tiledVectorIndex])[2]}); + // } + // polyscope::registerCurveNetworkLine("temp", edgePoints); + // polyscope::show(); + if (tiledIncidentVectors.size() == 1) { + continue; + } + + std::vector thetaAnglesOfIncidentVectors(tiledIncidentVectors.size()); + std::transform(tiledIncidentVectors.begin(), + tiledIncidentVectors.end(), + thetaAnglesOfIncidentVectors.begin(), + [](const VectorType &v) { return vcg::Point2d(v[0], v[1]).Angle(); }); + //sort them using theta angles + std::sort(thetaAnglesOfIncidentVectors.begin(), thetaAnglesOfIncidentVectors.end()); + + //find nodes that contain incident edges with relative angles less than the threshold + const double angleThreshold_rad = vcg::math::ToRad(angleThreshold_degrees); + for (int thetaAngleIndex = 1; thetaAngleIndex < thetaAnglesOfIncidentVectors.size(); + thetaAngleIndex++) { + const double absAngleDifference = std::abs( + thetaAnglesOfIncidentVectors[thetaAngleIndex] + - thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]); + if (absAngleDifference < angleThreshold_rad + && absAngleDifference > vcg::math::ToRad(0.01)) { + // std::cout << "Found angDiff:" << absAngleDifference << std::endl; + // vert[vi].C() = vcg::Color4b::Magenta; + // hasAngleSmallerThan = true; + return true; + } + } + const double firstLastPairAngleDiff = std::abs( + thetaAnglesOfIncidentVectors[0] + - thetaAnglesOfIncidentVectors[thetaAnglesOfIncidentVectors.size() - 1]); + if (firstLastPairAngleDiff < angleThreshold_rad && firstLastPairAngleDiff > 0.01) { + return true; + } + } + return false; +} +bool PatternGeometry::hasValenceGreaterThan(const std::vector &numberOfNodesPerSlot, + const size_t &valenceThreshold) +{ + //Initialize helping structs + if (nodeToSlotMap.empty()) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap); + } + if (correspondingNode.empty()) { + constructCorrespondingNodeMap(numberOfNodesPerSlot); + } + if (nodeTiledValence.empty()) { + constructNodeToTiledValenceMap(numberOfNodesPerSlot); + } + + //Check tiled valence of the pattern's nodes + bool hasNodeWithValenceGreaterThanDesired = false; + for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + const size_t tiledValence = nodeTiledValence[nodeIndex]; + if (tiledValence > valenceThreshold) { + vert[nodeIndex].C() = vcg::Color4b::LightRed; + hasNodeWithValenceGreaterThanDesired = true; + } + } + return hasNodeWithValenceGreaterThanDesired; +} + bool PatternGeometry::hasDanglingEdges( const std::vector &numberOfNodesPerSlot) { - if (nodeSlot.empty()) { - FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeSlot); - } - if (correspondingNode.empty()) { - constructCorresponginNodeMap(numberOfNodesPerSlot); - } + //Initialize helping structs + if (nodeToSlotMap.empty()) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap); + } + if (correspondingNode.empty()) { + constructCorrespondingNodeMap(numberOfNodesPerSlot); + } if (nodeTiledValence.empty()) { constructNodeToTiledValenceMap(numberOfNodesPerSlot); } + + //Check tiled valence of the pattern's nodes bool hasDanglingEdges = false; for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { const size_t tiledValence = nodeTiledValence[nodeIndex]; @@ -343,36 +527,37 @@ bool PatternGeometry::hasUntiledDanglingEdges() { // TODO: The function expects that the numberOfNodesPerSlot follows a specific // format and that the vertex container was populated in a particular order. -void PatternGeometry::constructCorresponginNodeMap( - const std::vector &numberOfNodesPerSlot) { - assert(vn != 0 && !nodeSlot.empty() && correspondingNode.empty() && - numberOfNodesPerSlot.size() == 7); +void PatternGeometry::constructCorrespondingNodeMap(const std::vector &numberOfNodesPerSlot) +{ + assert(vn != 0 && !nodeToSlotMap.empty() && correspondingNode.empty() + && numberOfNodesPerSlot.size() == 7); - for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { - const size_t slotIndex = nodeSlot[nodeIndex]; - if (slotIndex == 1) { - correspondingNode[nodeIndex] = nodeIndex + 1; - } else if (slotIndex == 2) { - correspondingNode[nodeIndex] = nodeIndex - 1; - } else if (slotIndex == 3) { - const size_t numberOfNodesBefore = - nodeIndex - std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.begin() + 3, 0); - correspondingNode[nodeIndex] = - std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.begin() + 6, 0) - - 1 - numberOfNodesBefore; - } else if (slotIndex == 5) { - const size_t numberOfNodesAfter = - std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.begin() + 6, 0) - - 1 - nodeIndex; - correspondingNode[nodeIndex] = - numberOfNodesAfter + std::accumulate(numberOfNodesPerSlot.begin(), - numberOfNodesPerSlot.begin() + 3, - 0); + for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + const size_t slotIndex = nodeToSlotMap[nodeIndex]; + if (slotIndex == 1) { + correspondingNode[nodeIndex] = nodeIndex + 1; + } else if (slotIndex == 2) { + correspondingNode[nodeIndex] = nodeIndex - 1; + } else if (slotIndex == 3) { + const size_t numberOfNodesBefore = nodeIndex + - std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, + 0); + correspondingNode[nodeIndex] = std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, + 0) + - 1 - numberOfNodesBefore; + } else if (slotIndex == 5) { + const size_t numberOfNodesAfter = std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, + 0) + - 1 - nodeIndex; + correspondingNode[nodeIndex] = numberOfNodesAfter + + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, + 0); + } } - } } bool PatternGeometry::isFullyConnectedWhenTiled() { @@ -383,16 +568,16 @@ bool PatternGeometry::isFullyConnectedWhenTiled() { // a single cc bool bottomInterfaceIsConnected = false; for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { - if (nodeSlot[nodeIndex] == 1 || nodeSlot[nodeIndex] == 4 || - nodeSlot[nodeIndex] == 2) { - std::vector connectedEdges; - vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges); - const size_t nodeValence = connectedEdges.size(); - if (nodeValence != 0) { - bottomInterfaceIsConnected = true; - break; + if (nodeToSlotMap[nodeIndex] == 1 || nodeToSlotMap[nodeIndex] == 4 + || nodeToSlotMap[nodeIndex] == 2) { + std::vector connectedEdges; + vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges); + const size_t nodeValence = connectedEdges.size(); + if (nodeValence != 0) { + bottomInterfaceIsConnected = true; + break; + } } - } } if (!bottomInterfaceIsConnected) { return false; @@ -580,6 +765,7 @@ PatternGeometry::PatternGeometry(const std::string &filename, vcg::tri::UpdateTopology::VertexEdge(*this); baseTriangleHeight = computeBaseTriangleHeight(); + baseTriangle = computeBaseTriangle(); updateEigenEdgeAndVertices(); } @@ -589,18 +775,25 @@ double PatternGeometry::computeBaseTriangleHeight() const return vcg::Distance(vert[0].cP(), vert[interfaceNodeIndex].cP()); } -void PatternGeometry::addNormals() { - bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; - if (normalsAreAbsent) { - for (auto &v : vert) { - v.N() = CoordType(0, 0, 1); - } - } -} - vcg::Triangle3 PatternGeometry::getBaseTriangle() const { - double bottomEdgeHalfSize = computeBaseTriangleHeight() / std::tan(M_PI / 3); + return baseTriangle; +} + +void PatternGeometry::addNormals() +{ + bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; + if (normalsAreAbsent) { + for (auto &v : vert) { + v.N() = CoordType(0, 0, 1); + } + } +} + +vcg::Triangle3 PatternGeometry::computeBaseTriangle() const +{ + const double baseTriangleHeight = computeBaseTriangleHeight(); + double bottomEdgeHalfSize = baseTriangleHeight / std::tan(M_PI / 3); CoordType patternCoord0 = vert[0].cP(); CoordType patternBottomRight = vert[interfaceNodeIndex].cP() + CoordType(bottomEdgeHalfSize, 0, 0); @@ -615,69 +808,245 @@ PatternGeometry::PatternGeometry( add(numberOfNodesPerSlot, edges); addNormals(); baseTriangleHeight = computeBaseTriangleHeight(); + baseTriangle = computeBaseTriangle(); updateEigenEdgeAndVertices(); } -bool PatternGeometry::createHoneycombAtom() { - VCGEdgeMesh honeycombQuarter; - const VCGEdgeMesh::CoordType n(0, 0, 1); - const double H = 0.2; - const double height = 1.5 * H; - const double width = 0.2; - const double theta = 70; - const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2; - vcg::tri::Allocator::AddVertex( - honeycombQuarter, VCGEdgeMesh::CoordType(0, height / 2, 0), n); - vcg::tri::Allocator::AddVertex( - honeycombQuarter, VCGEdgeMesh::CoordType(0, H / 2 - dy, 0), n); - vcg::tri::Allocator::AddVertex( - honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, H / 2, 0), n); - vcg::tri::Allocator::AddVertex( - honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, 0, 0), n); - vcg::tri::Allocator::AddEdge(honeycombQuarter, 0, 1); - vcg::tri::Allocator::AddEdge(honeycombQuarter, 1, 2); - vcg::tri::Allocator::AddEdge(honeycombQuarter, 2, 3); +std::shared_ptr PatternGeometry::tilePattern( + std::vector &patterns, + const std::vector &connectToNeighborsVi, + const VCGPolyMesh &tileInto, + const std::vector &perSurfaceFacePatternIndices, + std::vector &tileIntoEdgesToTiledVi, + std::vector> &perPatternIndexToTiledPatternEdgeIndex) +{ + perPatternIndexToTiledPatternEdgeIndex.resize(patterns.size()); + std::shared_ptr pTiledPattern(new PatternGeometry); + std::vector> tileIntoEdgeToInterfaceVi( + tileInto.EN()); //ith element contains all the interface nodes that lay on the ith edge of tileInto + //Compute the barycentric coords of the verts in the base triangle pattern + std::vector> barycentricCoordinates(patterns.size()); + for (size_t patternIndex = 0; patternIndex < patterns.size(); patternIndex++) { + const PatternGeometry &pattern = patterns[patternIndex]; + const vcg::Triangle3 &baseTriangle = pattern.getBaseTriangle(); + barycentricCoordinates[patternIndex].resize(pattern.VN()); + for (int vi = 0; vi < pattern.VN(); vi++) { + CoordType barycentricCoords_vi; + vcg::InterpolationParameters, double>(baseTriangle, + pattern.vert[vi].cP(), + barycentricCoords_vi); + barycentricCoordinates[patternIndex][vi] = barycentricCoords_vi; + } + } + VCGTriMesh tileIntoEdgeMesh; - VCGEdgeMesh honeycombAtom; - // Top right - vcg::tri::Append::MeshCopy(honeycombAtom, - honeycombQuarter); - // Bottom right - vcg::Matrix44d rotM; - rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); - vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); - vcg::tri::Append::Mesh(honeycombAtom, - honeycombQuarter); - // Bottom left - rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0)); - vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); - vcg::tri::Append::Mesh(honeycombAtom, - honeycombQuarter); - // Top left - rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); - vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); - vcg::tri::Append::Mesh(honeycombAtom, - honeycombQuarter); + assert(vcg::tri::HasFEAdjacency(tileInto)); + assert(vcg::tri::HasFVAdjacency(tileInto)); + for (const VCGPolyMesh::FaceType &f : tileInto.face) { + CoordType centerOfFace(0, 0, 0); + for (size_t vi = 0; vi < f.VN(); vi++) { + centerOfFace = centerOfFace + f.cP(vi); + } + centerOfFace /= f.VN(); + vcg::tri::Allocator::AddVertex(tileIntoEdgeMesh, + centerOfFace, + vcg::Color4b::Yellow); - for (VertexType &v : honeycombAtom.vert) { - v.P()[2] = 0; - } + std::vector firstInFanConnectToNeighbor_vi(connectToNeighborsVi); + for (int &vi : firstInFanConnectToNeighbor_vi) { + vi += pTiledPattern->VN(); + } + const size_t facePatternIndex = perSurfaceFacePatternIndices[tileInto.getIndex(f)]; + ConstPatternGeometry &pattern = patterns[facePatternIndex]; - return true; + for (size_t vi = 0; vi < f.VN(); vi++) { + auto ep = f.FEp(vi); + assert(vcg::tri::IsValidPointer(tileInto, ep)); + std::vector meshTrianglePoints{centerOfFace, + f.cP(vi), + vi + 1 == f.VN() ? f.cP(0) : f.cP(vi + 1)}; + // std::vector meshTrianglePoints{centerOfFace, ep->cP(0), ep->cP(1)}; + auto fit = vcg::tri::Allocator::AddFace(tileIntoEdgeMesh, + meshTrianglePoints[0], + meshTrianglePoints[1], + meshTrianglePoints[2]); + // CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) + // ^ (meshTrianglePoints[2] - meshTrianglePoints[0])) + // .Normalize(); + //TODO: in planar surfaces I should compute the face of triangle + CoordType faceNormal = f.cN(); + + fit->N() = faceNormal; + PatternGeometry transformedPattern; + transformedPattern.copy(pattern); + //Transform the base triangle nodes to the mesh triangle using barycentric coords + for (int vi = 0; vi < transformedPattern.VN(); vi++) { + transformedPattern.vert[vi].P() = CoordType( + meshTrianglePoints[0] * barycentricCoordinates[facePatternIndex][vi][0] + + meshTrianglePoints[1] * barycentricCoordinates[facePatternIndex][vi][1] + + meshTrianglePoints[2] * barycentricCoordinates[facePatternIndex][vi][2]); + } + + for (VertexType &v : transformedPattern.vert) { + v.N() = faceNormal; + } + + vcg::tri::Append::Remap remap; + vcg::tri::Append::Mesh(*pTiledPattern, + transformedPattern, + remap); + for (size_t ei = 0; ei < pattern.EN(); ei++) { + perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back(remap.edge[ei]); + } + const size_t ei = tileInto.getIndex(ep); + tileIntoEdgeToInterfaceVi[ei].push_back(remap.vert[pattern.interfaceNodeIndex]); + //Add edges for connecting the desired vertices + if (!connectToNeighborsVi.empty()) { + if (vi + 1 == f.VN()) { + for (int connectToNeighborIndex = 0; + connectToNeighborIndex < connectToNeighborsVi.size(); + connectToNeighborIndex++) { + auto eIt = vcg::tri::Allocator::AddEdge( + *pTiledPattern, + firstInFanConnectToNeighbor_vi[connectToNeighborIndex], + pTiledPattern->VN() - pattern.VN() + + connectToNeighborsVi[connectToNeighborIndex]); + perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back( + pTiledPattern->getIndex(*eIt)); + } + } + if (vi != 0) { + for (int connectToNeighborIndex = 0; + connectToNeighborIndex < connectToNeighborsVi.size(); + connectToNeighborIndex++) { + auto eIt = vcg::tri::Allocator::AddEdge( + *pTiledPattern, + pTiledPattern->VN() - 2 * pattern.VN() + + connectToNeighborsVi[connectToNeighborIndex], + pTiledPattern->VN() - pattern.VN() + + connectToNeighborsVi[connectToNeighborIndex]); + perPatternIndexToTiledPatternEdgeIndex[facePatternIndex].push_back( + pTiledPattern->getIndex(*eIt)); + } + } + } + // tiledPattern.updateEigenEdgeAndVertices(); + // tiledPattern.registerForDrawing(); + + // polyscope::show(); + } + } + vcg::tri::Allocator::PointerUpdater pu_vertices; + vcg::tri::Allocator::PointerUpdater pu_edges; + pTiledPattern->removeDuplicateVertices(pu_vertices, pu_edges); + //Update perPattern + // for (std::vector &tiledPatternEdges : perPatternIndexToTiledPatternEdgeIndex) { + // for (size_t &ei : tiledPatternEdges) { + // ei = pu_edges.remap[ei]; + // } + + // auto end = tiledPatternEdges.end(); + // for (auto it = tiledPatternEdges.begin(); it != end; ++it) { + // end = std::remove(it + 1, end, *it); + // } + + // tiledPatternEdges.erase(end, tiledPatternEdges.end()); + // } + const size_t sumOfEdgeIndices = std::accumulate(perPatternIndexToTiledPatternEdgeIndex.begin(), + perPatternIndexToTiledPatternEdgeIndex.end(), + 0, + [](const size_t &sum, + const std::vector &v) { + return sum + v.size(); + }); + + assert(pTiledPattern->EN() == sumOfEdgeIndices); + + tileIntoEdgesToTiledVi.clear(); + tileIntoEdgesToTiledVi.resize(tileInto.EN()); + for (int ei = 0; ei < tileInto.EN(); ei++) { + const std::vector interfaceVis = tileIntoEdgeToInterfaceVi[ei]; + assert(interfaceVis.size() == 1 || interfaceVis.size() == 2); + assert( + interfaceVis.size() == 1 + || (interfaceVis.size() == 2 + && (pu_vertices.remap[interfaceVis[0]] == std::numeric_limits::max() + || pu_vertices.remap[interfaceVis[1]] == std::numeric_limits::max()))); + tileIntoEdgesToTiledVi[ei] = pu_vertices.remap[interfaceVis[0]]; + } + + pTiledPattern->deleteDanglingVertices(); + vcg::tri::Allocator::CompactEveryVector(*pTiledPattern); + pTiledPattern->updateEigenEdgeAndVertices(); + pTiledPattern->save(); + return pTiledPattern; } -void PatternGeometry::copy(PatternGeometry ©From) { - VCGEdgeMesh::copy(copyFrom); - baseTriangleHeight = copyFrom.getBaseTriangleHeight(); +bool PatternGeometry::createHoneycombAtom() +{ + VCGEdgeMesh honeycombQuarter; + const VCGEdgeMesh::CoordType n(0, 0, 1); + const double H = 0.2; + const double height = 1.5 * H; + const double width = 0.2; + const double theta = 70; + const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2; + vcg::tri::Allocator::AddVertex(honeycombQuarter, + VCGEdgeMesh::CoordType(0, height / 2, 0), + n); + vcg::tri::Allocator::AddVertex(honeycombQuarter, + VCGEdgeMesh::CoordType(0, H / 2 - dy, 0), + n); + vcg::tri::Allocator::AddVertex(honeycombQuarter, + VCGEdgeMesh::CoordType(width / 2, H / 2, 0), + n); + vcg::tri::Allocator::AddVertex(honeycombQuarter, + VCGEdgeMesh::CoordType(width / 2, 0, 0), + n); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 0, 1); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 1, 2); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 2, 3); + + VCGEdgeMesh honeycombAtom; + // Top right + vcg::tri::Append::MeshCopy(honeycombAtom, honeycombQuarter); + // Bottom right + vcg::Matrix44d rotM; + rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, honeycombQuarter); + // Bottom left + rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, honeycombQuarter); + // Top left + rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, honeycombQuarter); + + for (VertexType &v : honeycombAtom.vert) { + v.P()[2] = 0; + } + + return true; } -void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize,const int& interfaceNodeIndex) { - this->baseTriangleHeight = desiredBaseTriangleCentralEdgeSize; - const double baseTriangleCentralEdgeSize = getBaseTriangleHeight(); - const double scaleRatio = - desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize; +void PatternGeometry::copy(PatternGeometry ©From) +{ + VCGEdgeMesh::copy(copyFrom); + baseTriangleHeight = copyFrom.getBaseTriangleHeight(); + baseTriangle = copyFrom.getBaseTriangle(); + interfaceNodeIndex = copyFrom.interfaceNodeIndex; +} - vcg::tri::UpdatePosition::Scale(*this, scaleRatio); +void PatternGeometry::scale(const double &desiredBaseTriangleCentralEdgeSize, + const int &interfaceNodeIndex) +{ + const double baseTriangleCentralEdgeSize = computeBaseTriangleHeight(); + const double scaleRatio = desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize; + vcg::tri::UpdatePosition::Scale(*this, scaleRatio); + baseTriangle = computeBaseTriangle(); + baseTriangleHeight = computeBaseTriangleHeight(); } void PatternGeometry::createFan(const std::vector &connectToNeighborsVi, const size_t &fanSize) @@ -718,6 +1087,7 @@ void PatternGeometry::createFan(const std::vector &connectToNeighborsVi, co } } -double PatternGeometry::getBaseTriangleHeight() const { - return baseTriangleHeight; +double PatternGeometry::getBaseTriangleHeight() const +{ + return baseTriangleHeight; } diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 2902044..4b19f44 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -7,7 +7,11 @@ #include "vcgtrimesh.hpp" #include "polymesh.hpp" -class PatternGeometry : public VCGEdgeMesh { +class PatternGeometry; +using ConstPatternGeometry = PatternGeometry; + +class PatternGeometry : public VCGEdgeMesh +{ private: size_t computeTiledValence(const size_t &nodeIndex, @@ -19,28 +23,30 @@ private: double baseTriangleHeight; double computeBaseTriangleHeight() const; - const int interfaceNodeIndex{3}; - const size_t fanSize{6}; + inline static size_t fanSize{6}; std::vector vertices; const double triangleEdgeSize{1}; // radius edge - std::unordered_map nodeSlot; + std::unordered_map nodeToSlotMap; std::unordered_map correspondingNode; std::unordered_map nodeTiledValence; + vcg::Triangle3 baseTriangle; void - constructCorresponginNodeMap(const std::vector &numberOfNodesPerSlot); + constructCorrespondingNodeMap(const std::vector &numberOfNodesPerSlot); - void constructNodeToTiledValenceMap( - const std::vector &numberOfNodesPerSlot); + void constructNodeToTiledValenceMap(const std::vector &numberOfNodesPerSlot); -public: + std::vector getEdgeVectorsWithVertexAsOrigin(std::vector &edgePointers, + const int &vi); + + public: PatternGeometry(); /*The following function should be a copy constructor with * a const argument but this is impossible due to the * vcglib interface. * */ PatternGeometry(PatternGeometry &other); - bool save(const std::string plyFilename); + bool load(const std::string &plyFilename) override; void add(const std::vector &vertices); void add(const std::vector &edges); void add(const std::vector &vertices, @@ -51,6 +57,8 @@ public: constructVertexVector(const std::vector &numberOfNodesPerSlot, const size_t &fanSize, const double &triangleEdgeSize); bool hasDanglingEdges(const std::vector &numberOfNodesPerSlot); + bool hasValenceGreaterThan(const std::vector &numberOfNodesPerSlot, + const size_t &valenceThreshold); std::vector getVertices() const; static PatternGeometry createFan(PatternGeometry &pattern); static PatternGeometry createTile(PatternGeometry &pattern); @@ -81,7 +89,6 @@ public: bool createHoneycombAtom(); void copy(PatternGeometry ©From); - void tilePattern(PatternGeometry &pattern, VCGTriMesh &tileInto); void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto, const int &interfaceNodeIndex, @@ -90,243 +97,32 @@ public: void scale(const double &desiredBaseTriangleCentralEdgeSize, const int &interfaceNodeIndex); double getBaseTriangleHeight() const; - vcg::Triangle3 getBaseTriangle() const; + vcg::Triangle3 computeBaseTriangle() const; PatternGeometry(const std::vector &vertices, const std::vector &edges); - - // static std::shared_ptr tilePattern(PatternGeometry &pattern, - // VCGPolyMesh &tileInto, - // const int &interfaceNodeIndex, - // const bool &shouldDeleteDanglingEdges) - // { - // PatternGeometry tiledPattern; - // double bottomEdgeHalfSize = vcg::Distance(pattern.vert[0].cP(), - // pattern.vert[interfaceNodeIndex].cP()) - // / std::tan(M_PI / 3); - // CoordType patternCoord0 = pattern.vert[0].cP(); - // CoordType patternBottomRight = pattern.vert[interfaceNodeIndex].cP() - // + CoordType(bottomEdgeHalfSize, 0, 0); - // CoordType patternBottomLeft = pattern.vert[interfaceNodeIndex].cP() - // - CoordType(bottomEdgeHalfSize, 0, 0); - // std::vector patternTrianglePoints{patternCoord0, - // patternBottomRight, - // patternBottomLeft}; - // CoordType pointOnPattern = patternCoord0 + (patternBottomLeft - patternCoord0) - // ^ (patternBottomRight - patternCoord0); - // // patternTrianglePoints.push_back(pointOnPattern); - - // VCGTriMesh tileIntoEdgeMesh; - - // for (VCGPolyMesh::FaceType &f : tileInto.face) { - // size_t numberOfNodes = 0; - // CoordType centerOfFace(0, 0, 0); - // for (size_t vi = 0; vi < f.VN(); vi++) { - // numberOfNodes++; - // centerOfFace = centerOfFace + f.cP(vi); - // } - // centerOfFace /= f.VN(); - // vcg::tri::Allocator::AddVertex(tileIntoEdgeMesh, - // centerOfFace, - // vcg::Color4b::Yellow); - // for (size_t vi = 0; vi < f.VN(); vi++) { - // // Compute transformation matrix M - // // vcg::Matrix44d M; - // std::vector meshTrianglePoints{centerOfFace, - // f.cP(vi), - // vi + 1 == f.VN() ? f.cP(0) - // : f.cP(vi + 1)}; - // CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) - // ^ (meshTrianglePoints[2] - meshTrianglePoints[0])) - // .Normalize(); - // auto fit = vcg::tri::Allocator::AddFace(tileIntoEdgeMesh, - // meshTrianglePoints[0], - // meshTrianglePoints[1], - // meshTrianglePoints[2]); - // fit->N() = faceNormal; - // CoordType pointOnTriMesh = meshTrianglePoints[0] - // + (meshTrianglePoints[1] - meshTrianglePoints[0]) - // ^ (meshTrianglePoints[2] - meshTrianglePoints[0]); - // vcg::Matrix44d M; - - // // meshTrianglePoints.push_back(pointOnTriMesh); - // // vcg::ComputeRigidMatchMatrix(meshTrianglePoints, patternTrianglePoints, M); - - // vcg::Matrix44d A_prime; - // A_prime[0][0] = meshTrianglePoints[0][0]; - // A_prime[1][0] = meshTrianglePoints[0][1]; - // A_prime[2][0] = meshTrianglePoints[0][2]; - // A_prime[3][0] = 1; - // A_prime[0][1] = meshTrianglePoints[1][0]; - // A_prime[1][1] = meshTrianglePoints[1][1]; - // A_prime[2][1] = meshTrianglePoints[1][2]; - // A_prime[3][1] = 1; - // A_prime[0][2] = meshTrianglePoints[2][0]; - // A_prime[1][2] = meshTrianglePoints[2][1]; - // A_prime[2][2] = meshTrianglePoints[2][2]; - // A_prime[3][2] = 1; - // A_prime[0][3] = pointOnTriMesh[0]; - // A_prime[1][3] = pointOnTriMesh[1]; - // A_prime[2][3] = pointOnTriMesh[2]; - // A_prime[3][3] = 1; - // vcg::Matrix44d A; - // A[0][0] = patternTrianglePoints[0][0]; - // A[1][0] = patternTrianglePoints[0][1]; - // A[2][0] = patternTrianglePoints[0][2]; - // A[3][0] = 1; - // A[0][1] = patternTrianglePoints[1][0]; - // A[1][1] = patternTrianglePoints[1][1]; - // A[2][1] = patternTrianglePoints[1][2]; - // A[3][1] = 1; - // A[0][2] = patternTrianglePoints[2][0]; - // A[1][2] = patternTrianglePoints[2][1]; - // A[2][2] = patternTrianglePoints[2][2]; - // A[3][2] = 1; - // A[0][3] = pointOnPattern[0]; - // A[1][3] = pointOnPattern[1]; - // A[2][3] = pointOnPattern[2]; - // A[3][3] = 1; - // M = A_prime * vcg::Inverse(A); - - // PatternGeometry transformedPattern; - // transformedPattern.copy(pattern); - // vcg::tri::UpdatePosition::Matrix(transformedPattern, M); - // for (VertexType &v : transformedPattern.vert) { - // v.N() = faceNormal; - // } - // // const double innerHexagonInitialRotationAngle = 30; - // // vcg::Matrix44d Rlocal; - // // Rlocal.SetRotateDeg( - // // (0 - // // // optimizationResults.optimalXNameValuePairs["HexAngle"] - // // - innerHexagonInitialRotationAngle), - // // VectorType(0, 0, 1)); - // // vcg::Matrix44d T; - // // T.SetTranslate(-transformedPattern.vert[0].cP()); - // // vcg::Matrix44d Trev; - // // Trev.SetTranslate(transformedPattern.vert[0].cP()); - // // auto R = T * Rlocal * Trev; - // // transformedPattern.vert[1].P() - // // = R * meshTrianglePoints[2] - // // * 0.5; //optimizationResults.optimalXNameValuePairs["HexSize"]; - // // transformedPattern.vert[5].P() - // // = R * meshTrianglePoints[1] - // // * 0.5; //optimizationResults.optimalXNameValuePairs["HexSize"]; - - // vcg::tri::Append::Mesh(tiledPattern, transformedPattern); - // tiledPattern.updateEigenEdgeAndVertices(); - // tiledPattern.registerForDrawing(); - - // polyscope::show(); - // } - // } - // tiledPattern.removeDuplicateVertices(); - // tiledPattern.deleteDanglingVertices(); - // vcg::tri::Allocator::CompactEveryVector(tiledPattern); - // tiledPattern.updateEigenEdgeAndVertices(); - // return std::make_shared(tiledPattern); - // } + // static std::shared_ptr tilePattern( + // std::vector &pattern, + // const std::vector &connectToNeighborsVi, + // const VCGPolyMesh &tileInto, + // std::vector &tileIntoEdgesToTiledVi, + // std::vector> &perPatternEdges); + void createFan(const std::vector &connectToNeighborsVi = std::vector(), + const size_t &fanSize = 6); + int interfaceNodeIndex{3}; //TODO: Fix this. This should be automatically computed + bool hasAngleSmallerThanThreshold(const std::vector &numberOfNodesPerSlot, + const double &angleThreshold_degrees); + vcg::Triangle3 getBaseTriangle() const; static std::shared_ptr tilePattern(PatternGeometry &pattern, - const vcg::Triangle3 &baseTriangle, const std::vector &connectToNeighborsVi, - const VCGPolyMesh &tileInto) - { - std::shared_ptr pTiledPattern(new PatternGeometry); - - //Compute the barycentric coords of the verts in the base triangle pattern - std::vector barycentricCoordinates(pattern.VN()); - for (int vi = 0; vi < pattern.VN(); vi++) { - CoordType barycentricCoords_vi; - vcg::InterpolationParameters, double>(baseTriangle, - pattern.vert[vi].cP(), - barycentricCoords_vi); - barycentricCoordinates[vi] = barycentricCoords_vi; - } - VCGTriMesh tileIntoEdgeMesh; - - for (const VCGPolyMesh::FaceType &f : tileInto.face) { - CoordType centerOfFace(0, 0, 0); - for (size_t vi = 0; vi < f.VN(); vi++) { - centerOfFace = centerOfFace + f.cP(vi); - } - centerOfFace /= f.VN(); - vcg::tri::Allocator::AddVertex(tileIntoEdgeMesh, - centerOfFace, - vcg::Color4b::Yellow); - - std::vector firstInFanConnectToNeighbor_vi(connectToNeighborsVi); - for (int &vi : firstInFanConnectToNeighbor_vi) { - vi += pTiledPattern->VN(); - } - - for (size_t vi = 0; vi < f.VN(); vi++) { - std::vector meshTrianglePoints{centerOfFace, - f.cP(vi), - vi + 1 == f.VN() ? f.cP(0) - : f.cP(vi + 1)}; - auto fit = vcg::tri::Allocator::AddFace(tileIntoEdgeMesh, - meshTrianglePoints[0], - meshTrianglePoints[1], - meshTrianglePoints[2]); - CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) - ^ (meshTrianglePoints[2] - meshTrianglePoints[0])) - .Normalize(); - fit->N() = faceNormal; - - PatternGeometry transformedPattern; - transformedPattern.copy(pattern); - //Transform the base triangle nodes to the mesh triangle using barycentric coords - for (int vi = 0; vi < transformedPattern.VN(); vi++) { - transformedPattern.vert[vi].P() = CoordType( - meshTrianglePoints[0] * barycentricCoordinates[vi][0] - + meshTrianglePoints[1] * barycentricCoordinates[vi][1] - + meshTrianglePoints[2] * barycentricCoordinates[vi][2]); - } - - for (VertexType &v : transformedPattern.vert) { - v.N() = faceNormal; - } - - vcg::tri::Append::Mesh(*pTiledPattern, - transformedPattern); - //Add edges for connecting the desired vertices - if (!connectToNeighborsVi.empty()) { - if (vi + 1 == f.VN()) { - for (int connectToNeighborIndex = 0; - connectToNeighborIndex < connectToNeighborsVi.size(); - connectToNeighborIndex++) { - vcg::tri::Allocator::AddEdge( - *pTiledPattern, - firstInFanConnectToNeighbor_vi[connectToNeighborIndex], - pTiledPattern->VN() - pattern.VN() - + connectToNeighborsVi[connectToNeighborIndex]); - } - } - if (vi != 0) { - for (int connectToNeighborIndex = 0; - connectToNeighborIndex < connectToNeighborsVi.size(); - connectToNeighborIndex++) { - vcg::tri::Allocator::AddEdge( - *pTiledPattern, - pTiledPattern->VN() - 2 * pattern.VN() - + connectToNeighborsVi[connectToNeighborIndex], - pTiledPattern->VN() - pattern.VN() - + connectToNeighborsVi[connectToNeighborIndex]); - } - } - } - // tiledPattern.updateEigenEdgeAndVertices(); - // tiledPattern.registerForDrawing(); - - // polyscope::show(); - } - } - pTiledPattern->removeDuplicateVertices(); - pTiledPattern->deleteDanglingVertices(); - vcg::tri::Allocator::CompactEveryVector(*pTiledPattern); - pTiledPattern->updateEigenEdgeAndVertices(); - return pTiledPattern; - } - void createFan(const std::vector &connectToNeighborsVi=std::vector(), const size_t &fanSize=6); + const VCGPolyMesh &tileInto, + std::vector &tileIntoEdgesToTiledVi); + static std::shared_ptr tilePattern( + std::vector &patterns, + const std::vector &connectToNeighborsVi, + const VCGPolyMesh &tileInto, + const std::vector &perSurfaceFacePatternIndices, + std::vector &tileIntoEdgesToTiledVi, + std::vector> &perPatternIndexTiledPatternEdgeIndex); }; #endif // FLATPATTERNGEOMETRY_HPP diff --git a/utilities.hpp b/utilities.hpp index 2424397..8a1be83 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -108,6 +108,11 @@ struct Vector6d : public std::array { return Eigen::Vector3d(this->operator[](0), this->operator[](1), this->operator[](2)); } + + Eigen::Vector3d getRotation() const + { + return Eigen::Vector3d(this->operator[](3), this->operator[](4), this->operator[](5)); + } }; namespace Utilities { @@ -137,6 +142,20 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { return m; } +inline std::vector fromEigenMatrix(const Eigen::MatrixXd &m) +{ + std::vector v(m.rows()); + + for (size_t vi = 0; vi < m.rows(); vi++) { + const Eigen::RowVectorXd &row = m.row(vi); + for (size_t i = 0; i < 6; i++) { + v[vi][i] = row(i); + } + } + + return v; +} + // std::string convertToLowercase(const std::string &s) { // std::string lowercase; // std::transform(s.begin(), s.end(), lowercase.begin(), @@ -167,51 +186,93 @@ inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { #ifdef POLYSCOPE_DEFINED #include "polyscope/curve_network.h" +#include "polyscope/pick.h" #include "polyscope/polyscope.h" +#include -inline void deinitPolyscope() { - if (!polyscope::state::initialized) { - return; - } +namespace PolyscopeInterface { +inline struct GlobalPolyscopeData +{ + std::vector> userCallbacks; +} globalPolyscopeData; - polyscope::shutdown(); +inline void mainCallback() +{ + ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, + // instead of full width. Must have + // matching PopItemWidth() below. + + for (std::function userCallback : globalPolyscopeData.userCallbacks) { + userCallback(); + } + + ImGui::PopItemWidth(); } -inline void initPolyscope() { - if (polyscope::state::initialized) { - return; - } - polyscope::init(); - polyscope::options::groundPlaneEnabled = false; - polyscope::view::upDir = polyscope::view::UpDir::ZUp; +inline void addUserCallback(const std::function &userCallback) +{ + globalPolyscopeData.userCallbacks.push_back(userCallback); } -inline void registerWorldAxes() { - initPolyscope(); +inline void deinitPolyscope() +{ + if (!polyscope::state::initialized) { + return; + } - Eigen::MatrixX3d axesPositions(4, 3); - axesPositions.row(0) = Eigen::Vector3d(0, 0, 0); - axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0); - axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0); - axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale); - - Eigen::MatrixX2i axesEdges(3, 2); - axesEdges.row(0) = Eigen::Vector2i(0, 1); - axesEdges.row(1) = Eigen::Vector2i(0, 2); - axesEdges.row(2) = Eigen::Vector2i(0, 3); - Eigen::MatrixX3d axesColors(3, 3); - axesColors.row(0) = Eigen::Vector3d(1, 0, 0); - axesColors.row(1) = Eigen::Vector3d(0, 1, 0); - axesColors.row(2) = Eigen::Vector3d(0, 0, 1); - - const std::string worldAxesName = "World Axes"; - polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges); - polyscope::getCurveNetwork(worldAxesName)->setRadius(0.0001, false); - const std::string worldAxesColorName = worldAxesName + " Color"; - polyscope::getCurveNetwork(worldAxesName) - ->addEdgeColorQuantity(worldAxesColorName, axesColors) - ->setEnabled(true); + polyscope::shutdown(); } + +inline void init() +{ + if (polyscope::state::initialized) { + return; + } + polyscope::init(); + polyscope::options::groundPlaneEnabled = false; + polyscope::view::upDir = polyscope::view::UpDir::ZUp; + + polyscope::state::userCallback = &mainCallback; +} +using PolyscopeLabel = std::string; +inline std::pair getSelection() +{ + std::pair selection = polyscope::pick::getSelection(); + if (selection.first == nullptr) { + return std::make_pair(std::string(), 0); + } + return std::make_pair(selection.first->name, selection.second); +} + +inline void registerWorldAxes() +{ + PolyscopeInterface::init(); + + Eigen::MatrixX3d axesPositions(4, 3); + axesPositions.row(0) = Eigen::Vector3d(0, 0, 0); + axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0); + axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0); + axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale); + + Eigen::MatrixX2i axesEdges(3, 2); + axesEdges.row(0) = Eigen::Vector2i(0, 1); + axesEdges.row(1) = Eigen::Vector2i(0, 2); + axesEdges.row(2) = Eigen::Vector2i(0, 3); + Eigen::MatrixX3d axesColors(3, 3); + axesColors.row(0) = Eigen::Vector3d(1, 0, 0); + axesColors.row(1) = Eigen::Vector3d(0, 1, 0); + axesColors.row(2) = Eigen::Vector3d(0, 0, 1); + + const std::string worldAxesName = "World Axes"; + polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges); + polyscope::getCurveNetwork(worldAxesName)->setRadius(0.0001, false); + const std::string worldAxesColorName = worldAxesName + " Color"; + polyscope::getCurveNetwork(worldAxesName) + ->addEdgeColorQuantity(worldAxesColorName, axesColors) + ->setEnabled(true); +} +} // namespace PolyscopeInterface + #endif // namespace ConfigurationFile { @@ -232,4 +293,13 @@ template std::string toString(const T &v) { std::to_string(v[2]) + ")"; } +template +std::string to_string_with_precision(const T a_value, const int n = 2) +{ + std::ostringstream out; + out.precision(n); + out << std::fixed << a_value; + return out.str(); +} + #endif // UTILITIES_H diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp index dc938cb..53de708 100755 --- a/vcgtrimesh.cpp +++ b/vcgtrimesh.cpp @@ -7,6 +7,7 @@ #include bool VCGTriMesh::load(const std::string &filename) { + assert(std::filesystem::exists(filename)); // unsigned int mask = 0; // mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; // mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; @@ -85,12 +86,13 @@ bool VCGTriMesh::save(const std::string plyFilename) return true; } +#ifdef POLYSCOPE_DEFINED polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing( const std::optional> &desiredColor, const bool &shouldEnable) const { auto vertices = getVertices(); auto faces = getFaces(); - initPolyscope(); + PolyscopeInterface::init(); if (faces.rows() == 0) { auto edges = getEdges(); @@ -114,6 +116,7 @@ polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing( } return polyscopeHandle_mesh; } +#endif VCGTriMesh::VCGTriMesh() {} diff --git a/vcgtrimesh.hpp b/vcgtrimesh.hpp index 582a074..02a180e 100755 --- a/vcgtrimesh.hpp +++ b/vcgtrimesh.hpp @@ -1,9 +1,11 @@ #ifndef VCGTRIMESH_HPP #define VCGTRIMESH_HPP #include "mesh.hpp" -#include #include +#ifdef POLYSCOPE_DEFINED +#include +#endif using VertexIndex = size_t; class VCGTriMeshVertex; class VCGTriMeshEdge;