diff --git a/CMakeLists.txt b/CMakeLists.txt index c24ff9a..e43756a 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,6 +43,17 @@ download_project(PROJ threed-beam-fea ) add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR}) +##TBB +download_project(PROJ TBB + GIT_REPOSITORY https://github.com/wjakob/tbb.git + GIT_TAG master + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) +link_directories(${TBB_BINARY_DIR}) +message(${TBB_BINARY_DIR}) + ###Eigen 3 NOTE: Eigen is required on the system the code is ran find_package(Eigen3 3.3 REQUIRED) @@ -53,6 +64,7 @@ endif(MSVC) #link_directories(${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs) file(GLOB MySourcesFiles ${CMAKE_CURRENT_LIST_DIR}/*.hpp ${CMAKE_CURRENT_LIST_DIR}/*.cpp) +#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ltbb") add_library(${PROJECT_NAME} ${MySourcesFiles} ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp) @@ -64,9 +76,8 @@ target_include_directories(${PROJECT_NAME} ) if(${USE_POLYSCOPE}) - find_package(OpenMP REQUIRED) - target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA OpenMP::OpenMP_CXX) + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) else() - target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA) + target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) endif() target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs) diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index e33b7b4..9649614 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -222,6 +222,7 @@ void DRMSimulationModel::reset() mCurrentSimulationStep = 0; history.clear(); constrainedVertices.clear(); + isRigidSupport.clear(); pMesh.reset(); plotYValues.clear(); plotHandle.reset(); @@ -808,134 +809,111 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( //#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); - const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1); - const Element &element = pMesh->elements[e]; - const Element::LocalFrame &ef = element.frame; - const VectorType t3CrossN_j = Cross(ef.t3, ev_j.cN()); - const VectorType t3CrossN_jplus1 = Cross(ef.t3, ev_jplus1.cN()); - const double theta1_j = ef.t1.dot(t3CrossN_j); - const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1); - const double theta2_j = ef.t2.dot(t3CrossN_j); - const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1); - const bool shouldBreak = mCurrentSimulationStep == 12970; - const double theta3_j = computeTheta3(e, ev_j); - const double theta3_jplus1 = computeTheta3(e, ev_jplus1); - // element.rotationalDisplacements_j.theta1 = theta1_j; - // element.rotationalDisplacements_jplus1.theta1 = theta1_jplus1; - // element.rotationalDisplacements_j.theta2 = theta2_j; - // element.rotationalDisplacements_jplus1.theta2 = theta2_jplus1; - // element.rotationalDisplacements_j.theta3 = theta3_j; - // element.rotationalDisplacements_jplus1.theta3 = theta3_jplus1; - std::vector> internalForcesContributionFromThisEdge(4, - {-1, - Vector6d()}); - for (VertexIndex evi = 0; evi < 2; evi++) { - const SimulationMesh::VertexType &ev = *e.cV(evi); - const Node &edgeNode = pMesh->nodes[ev]; - internalForcesContributionFromThisEdge[evi].first = edgeNode.vi; + std::for_each( + std::execution::par_unseq, + pMesh->edge.begin(), + pMesh->edge.end(), + [&](const EdgeType &e) + // for (int ei = 0; ei < pMesh->EN(); ei++) + { + const int ei = pMesh->getIndex(e); + // const EdgeType &e = pMesh->edge[ei]; + const SimulationMesh::VertexType &ev_j = *e.cV(0); + const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1); + const Element &element = pMesh->elements[e]; + const Element::LocalFrame &ef = element.frame; + const VectorType t3CrossN_j = Cross(ef.t3, ev_j.cN()); + const VectorType t3CrossN_jplus1 = Cross(ef.t3, ev_jplus1.cN()); + const double theta1_j = ef.t1.dot(t3CrossN_j); + const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1); + const double theta2_j = ef.t2.dot(t3CrossN_j); + const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1); + const bool shouldBreak = mCurrentSimulationStep == 12970; + const double theta3_j = computeTheta3(e, ev_j); + const double theta3_jplus1 = computeTheta3(e, ev_jplus1); + // element.rotationalDisplacements_j.theta1 = theta1_j; + // element.rotationalDisplacements_jplus1.theta1 = theta1_jplus1; + // element.rotationalDisplacements_j.theta2 = theta2_j; + // element.rotationalDisplacements_jplus1.theta2 = theta2_jplus1; + // element.rotationalDisplacements_j.theta3 = theta3_j; + // element.rotationalDisplacements_jplus1.theta3 = theta3_jplus1; + std::vector> + internalForcesContributionFromThisEdge(4, {-1, Vector6d()}); + for (VertexIndex evi = 0; evi < 2; evi++) { + const SimulationMesh::VertexType &ev = *e.cV(evi); + const Node &edgeNode = pMesh->nodes[ev]; + internalForcesContributionFromThisEdge[evi].first = edgeNode.vi; - const VertexPointer &rev_j = edgeNode.referenceElement->cV(0); - const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1); - // refElemOtherVertex can be precomputed - const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j; - const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex]; - if (edgeNode.referenceElement != &e) { - internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi; - } - const size_t vi = edgeNode.vi; - // #pragma omp parallel for schedule(static) num_threads(6) - for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { - const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi] - && fixedVertices.at(edgeNode.vi).contains(dofi); - if (!isDofConstrainedFor_ev) { - DifferentiateWithRespectTo dui{ev, dofi}; - // Axial force computation - const double e_k = element.length - element.initialLength; - const double e_kDeriv = computeDerivativeElementLength(e, dui); - const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial; - -#ifdef POLYSCOPE_DEFINED - if (std::isnan(axialForce_dofi)) { - std::cerr << "nan in axial" << evi << std::endl; - } -#endif - - // Torsional force computation - const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi); - const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi); - const double theta1Diff = theta1_jplus1 - theta1_j; - const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; - const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff - * element.rigidity.torsional; -#ifdef POLYSCOPE_DEFINED - if (std::isnan(torsionalForce_dofi)) { - std::cerr << "nan in torsional" << evi << std::endl; - } -#endif - - // First bending force computation - ////theta2_j derivative - const double theta2_j_deriv = computeDerivativeTheta2(e, 0, evi, dofi); - ////theta2_jplus1 derivative - const double theta2_jplus1_deriv = computeDerivativeTheta2(e, 1, evi, dofi); - ////1st in bracket term - const double firstBendingForce_inBracketsTerm_0 = theta2_j_deriv * 2 * theta2_j; - ////2nd in bracket term - const double firstBendingForce_inBracketsTerm_1 = theta2_jplus1_deriv - * theta2_j; - ////3rd in bracket term - const double firstBendingForce_inBracketsTerm_2 = theta2_j_deriv - * theta2_jplus1; - ////4th in bracket term - const double firstBendingForce_inBracketsTerm_3 = 2 * theta2_jplus1_deriv - * theta2_jplus1; - // 3rd term computation - const double firstBendingForce_inBracketsTerm - = firstBendingForce_inBracketsTerm_0 + firstBendingForce_inBracketsTerm_1 - + firstBendingForce_inBracketsTerm_2 + firstBendingForce_inBracketsTerm_3; - const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm - * element.rigidity.firstBending; - - // Second bending force computation - ////theta2_j derivative - const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui); - ////theta2_jplus1 derivative - const double theta3_jplus1_deriv = computeDerivativeTheta3(e, ev_jplus1, dui); - ////1st in bracket term - const double secondBendingForce_inBracketsTerm_0 = theta3_j_deriv * 2 - * theta3_j; - ////2nd in bracket term - const double secondBendingForce_inBracketsTerm_1 = theta3_jplus1_deriv - * theta3_j; - ////3rd in bracket term - const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv - * theta3_jplus1; - ////4th in bracket term - const double secondBendingForce_inBracketsTerm_3 = 2 * theta3_jplus1_deriv - * theta3_jplus1; - // 3rd term computation - const double secondBendingForce_inBracketsTerm - = secondBendingForce_inBracketsTerm_0 + secondBendingForce_inBracketsTerm_1 - + secondBendingForce_inBracketsTerm_2 - + secondBendingForce_inBracketsTerm_3; - double secondBendingForce_dofi = secondBendingForce_inBracketsTerm - * element.rigidity.secondBending; - internalForcesContributionFromThisEdge[evi].second[dofi] - = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi - + torsionalForce_dofi; - } + const VertexPointer &rev_j = edgeNode.referenceElement->cV(0); + const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1); + // refElemOtherVertex can be precomputed + const VertexPointer &refElemOtherVertex = rev_j == &ev ? rev_jplus1 : rev_j; + const Node &refElemOtherVertexNode = pMesh->nodes[refElemOtherVertex]; if (edgeNode.referenceElement != &e) { - const bool isDofConstrainedFor_refElemOtherVertex - = isVertexConstrained[refElemOtherVertexNode.vi] - && fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi); - if (!isDofConstrainedFor_refElemOtherVertex) { - DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi}; - ////theta3_j derivative + internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi; + } + const size_t vi = edgeNode.vi; + // #pragma omp parallel for schedule(static) num_threads(6) + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + const bool isDofConstrainedFor_ev = isVertexConstrained[edgeNode.vi] + && fixedVertices.at(edgeNode.vi) + .contains(dofi); + if (!isDofConstrainedFor_ev) { + DifferentiateWithRespectTo dui{ev, dofi}; + // Axial force computation + const double e_k = element.length - element.initialLength; + const double e_kDeriv = computeDerivativeElementLength(e, dui); + const double axialForce_dofi = e_kDeriv * e_k * element.rigidity.axial; + +#ifdef POLYSCOPE_DEFINED + if (std::isnan(axialForce_dofi)) { + std::cerr << "nan in axial" << evi << std::endl; + } +#endif + + // Torsional force computation + const double theta1_j_deriv = computeDerivativeTheta1(e, 0, evi, dofi); + const double theta1_jplus1_deriv = computeDerivativeTheta1(e, 1, evi, dofi); + const double theta1Diff = theta1_jplus1 - theta1_j; + const double theta1DiffDerivative = theta1_jplus1_deriv - theta1_j_deriv; + const double torsionalForce_dofi = theta1DiffDerivative * theta1Diff + * element.rigidity.torsional; +#ifdef POLYSCOPE_DEFINED + if (std::isnan(torsionalForce_dofi)) { + std::cerr << "nan in torsional" << evi << std::endl; + } +#endif + + // First bending force computation + ////theta2_j derivative + const double theta2_j_deriv = computeDerivativeTheta2(e, 0, evi, dofi); + ////theta2_jplus1 derivative + const double theta2_jplus1_deriv = computeDerivativeTheta2(e, 1, evi, dofi); + ////1st in bracket term + const double firstBendingForce_inBracketsTerm_0 = theta2_j_deriv * 2 + * theta2_j; + ////2nd in bracket term + const double firstBendingForce_inBracketsTerm_1 = theta2_jplus1_deriv + * theta2_j; + ////3rd in bracket term + const double firstBendingForce_inBracketsTerm_2 = theta2_j_deriv + * theta2_jplus1; + ////4th in bracket term + const double firstBendingForce_inBracketsTerm_3 = 2 * theta2_jplus1_deriv + * theta2_jplus1; + // 3rd term computation + const double firstBendingForce_inBracketsTerm + = firstBendingForce_inBracketsTerm_0 + + firstBendingForce_inBracketsTerm_1 + + firstBendingForce_inBracketsTerm_2 + + firstBendingForce_inBracketsTerm_3; + const double firstBendingForce_dofi = firstBendingForce_inBracketsTerm + * element.rigidity.firstBending; + + // Second bending force computation + ////theta2_j derivative const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui); - ////theta3_jplus1 derivative + ////theta2_jplus1 derivative const double theta3_jplus1_deriv = computeDerivativeTheta3(e, ev_jplus1, dui); @@ -949,25 +927,61 @@ void DRMSimulationModel::updateResidualForcesOnTheFly( const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv * theta3_jplus1; ////4th in bracket term - const double secondBendingForce_inBracketsTerm_3 = theta3_jplus1_deriv * 2 + const double secondBendingForce_inBracketsTerm_3 = 2 * theta3_jplus1_deriv * theta3_jplus1; - - // 4th term computation + // 3rd term computation const double secondBendingForce_inBracketsTerm = secondBendingForce_inBracketsTerm_0 + secondBendingForce_inBracketsTerm_1 + secondBendingForce_inBracketsTerm_2 + secondBendingForce_inBracketsTerm_3; - const double secondBendingForce_dofi = secondBendingForce_inBracketsTerm - * element.rigidity.secondBending; - internalForcesContributionFromThisEdge[evi + 2].second[dofi] - = secondBendingForce_dofi; + double secondBendingForce_dofi = secondBendingForce_inBracketsTerm + * element.rigidity.secondBending; + internalForcesContributionFromThisEdge[evi].second[dofi] + = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi + + torsionalForce_dofi; + } + if (edgeNode.referenceElement != &e) { + const bool isDofConstrainedFor_refElemOtherVertex + = isVertexConstrained[refElemOtherVertexNode.vi] + && fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi); + if (!isDofConstrainedFor_refElemOtherVertex) { + DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi}; + ////theta3_j derivative + const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui); + ////theta3_jplus1 derivative + const double theta3_jplus1_deriv = computeDerivativeTheta3(e, + ev_jplus1, + dui); + ////1st in bracket term + const double secondBendingForce_inBracketsTerm_0 = theta3_j_deriv * 2 + * theta3_j; + ////2nd in bracket term + const double secondBendingForce_inBracketsTerm_1 = theta3_jplus1_deriv + * theta3_j; + ////3rd in bracket term + const double secondBendingForce_inBracketsTerm_2 = theta3_j_deriv + * theta3_jplus1; + ////4th in bracket term + const double secondBendingForce_inBracketsTerm_3 = theta3_jplus1_deriv + * 2 * theta3_jplus1; + + // 4th term computation + const double secondBendingForce_inBracketsTerm + = secondBendingForce_inBracketsTerm_0 + + secondBendingForce_inBracketsTerm_1 + + secondBendingForce_inBracketsTerm_2 + + secondBendingForce_inBracketsTerm_3; + const double secondBendingForce_dofi = secondBendingForce_inBracketsTerm + * element.rigidity.secondBending; + internalForcesContributionFromThisEdge[evi + 2].second[dofi] + = secondBendingForce_dofi; + } } } } - } - internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge; - } + internalForcesContributionsFromEachEdge[ei] = internalForcesContributionFromThisEdge; + }); //#pragma omp parallel for schedule(static) num_threads(8) @@ -1247,17 +1261,7 @@ void DRMSimulationModel::updateNodalAccelerations() for (VertexType &v : pMesh->vert) { Node &node = pMesh->nodes[v]; const VertexIndex vi = pMesh->getIndex(v); - for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { - if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.translational; - } else if (dofi == DoF::Nx) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalJ; - } else if (dofi == DoF::Ny) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI3; - } else if (dofi == DoF::Nr) { - node.acceleration[dofi] = node.force.residual[dofi] / node.mass.rotationalI2; - } - } + node.acceleration = node.force.residual / node.mass_6d; #ifdef POLYSCOPE_DEFINED if (std::isnan(node.acceleration.norm())) { std::cout << "acceleration " << vi << ":" << node.acceleration.toString() << std::endl; @@ -2578,3 +2582,19 @@ mesh->currentTotalPotentialEnergykN*/ return results; } + +bool DRMSimulationModel::Settings::save(const filesystem::__cxx11::path &folderPath) const +{ + bool returnValue = true; + std::filesystem::create_directories(folderPath); + nlohmann::json json; + json[jsonLabels.meshFilename] + = std::filesystem::relative(std::filesystem::path(meshFilename), + std::filesystem::path(jsonFilename).parent_path()) + .string(); + const std::string jsonFilename = "drmSettings.json"; + std::ofstream jsonFile(std::filesystem::path(folderPath).append(jsonFilename)); + jsonFile << json; +} + +bool DRMSimulationModel::Settings::load(const filesystem::__cxx11::path &filePath) const {} diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index f415b80..bb68a07 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -25,28 +25,32 @@ class DRMSimulationModel public: struct Settings { - bool isDebugMode{false}; - int debugModeStep{100000}; + // bool isDebugMode{false}; + std::optional debugModeStep{0}; bool shouldDraw{false}; bool beVerbose{false}; bool shouldCreatePlots{false}; - int drawingStep{1}; - double totalTranslationalKineticEnergyThreshold{1e-8}; - double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; - double residualForcesMovingAverageNormThreshold{1e-8}; + // int drawingStep{0}; + // double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; + // double residualForcesMovingAverageNormThreshold{1e-8}; double Dtini{0.1}; double xi{0.9969}; int maxDRMIterations{0}; - bool shouldUseTranslationalKineticEnergyThreshold{false}; - int gradualForcedDisplacementSteps{50}; - int desiredGradualExternalLoadsSteps{1}; + std::optional shouldUseTranslationalKineticEnergyThreshold; + // int gradualForcedDisplacementSteps{50}; + // int desiredGradualExternalLoadsSteps{1}; double gamma{0.8}; std::optional displacementCap; double totalResidualForcesNormThreshold{1e-3}; double totalExternalForcesNormPercentageTermination{1e-3}; - bool useAverage{false}; - double averageResidualForcesCriterionThreshold{1e-5}; + std::optional averageResidualForcesCriterionThreshold{1e-5}; Settings() {} + bool save(const std::filesystem::path &folderPath) const; + bool load(const std::filesystem::path &filePath) const; + bool useViscousDamping{false}; + bool useKineticDamping{true}; + struct JsonLabels + {}; }; private: @@ -158,9 +162,9 @@ private: void updateElementalFrames(); - VectorType computeDerivativeOfR(const EdgeType &e, const DifferentiateWithRespectTo &dui) const; + VectorType computeDerivativeOfR(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; - // bool isRigidSupport(const VertexType &v) const; static double computeDerivativeOfNorm(const VectorType &x, const VectorType &derivativeOfX); diff --git a/simulationmesh.hpp b/simulationmesh.hpp index cecadf8..a7b85c9 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -148,6 +148,8 @@ struct Node { }; Mass mass; + Vector6d mass_6d; + Vector6d damping_6d; VertexIndex vi; CoordType initialLocation; CoordType initialNormal; diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index dba4356..d62a085 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -191,7 +191,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced patternGeometryAllEdges.getVertices(), intersectingEdges, validEdges); - // statistics.print(setupString, perEdgeResultPath); + statistics.print(setupString, perEdgeResultPath); + statistics.reset(); } } else { std::cout << "Computing " + setupString << " with " << numberOfDesiredEdges << " edges." @@ -372,6 +373,20 @@ std::vector TopologyEnumerator::getValidEdges( return validEdges; } +void TopologyEnumerator::exportPattern(const std::filesystem::path &saveToPath, + PatternGeometry &patternGeometry, + const bool saveTilledPattern) const +{ + const std::string patternName = patternGeometry.getLabel(); + std::filesystem::create_directory(saveToPath); + patternGeometry.save(std::filesystem::path(saveToPath).append(patternName).string() + ".ply"); + if (saveTilledPattern) { + PatternGeometry tiledPatternGeometry = PatternGeometry::createTile(patternGeometry); + tiledPatternGeometry.save( + std::filesystem::path(saveToPath).append(patternName + "_tiled").string() + ".ply"); + } +} + void TopologyEnumerator::computeValidPatterns( const std::vector &numberOfNodesPerSlot, const size_t &numberOfDesiredEdges, @@ -418,6 +433,8 @@ void TopologyEnumerator::computeValidPatterns( // std::string previousPatternBinaryRepresentation(validEdges.size(),'0'); size_t patternIndex = 0; bool validPatternsExist = false; + const bool exportTilledPattern = false; + const bool saveCompressedFormat = false; do { patternIndex++; const std::string patternName = std::to_string(patternIndex); @@ -438,6 +455,7 @@ void TopologyEnumerator::computeValidPatterns( PatternGeometry patternGeometry; patternGeometry.add(allVertices, patternEdges); + patternGeometry.setLabel(patternName); // Check if pattern contains intersecting edges const bool patternContainsIntersectingEdges @@ -448,21 +466,13 @@ void TopologyEnumerator::computeValidPatterns( 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"); + exportPattern(std::filesystem::path(resultsPath).append("Intersecting"), + patternGeometry, + exportTilledPattern); } + } else { + continue; // should be uncommented in order to improve performance } - continue; // should be uncommented in order to improve performance } const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold @@ -470,18 +480,10 @@ void TopologyEnumerator::computeValidPatterns( 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"); + exportPattern(std::filesystem::path(resultsPath) + .append("ExceedingAngleThreshold"), + patternGeometry, + exportTilledPattern); } } else { continue; @@ -493,18 +495,9 @@ void TopologyEnumerator::computeValidPatterns( 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"); + auto highValencePath = std::filesystem::path(resultsPath) + .append("HighValencePatterns"); + exportPattern(highValencePath, patternGeometry, exportTilledPattern); } } else { continue; @@ -522,28 +515,19 @@ void TopologyEnumerator::computeValidPatterns( // Check dangling edges with vcg method // const bool vcg_tiledPatternHasDangling = // tiledPatternGeometry.hasUntiledDanglingEdges(); - if (tiledPatternHasDanglingEdges && !hasFloatingComponents /*&& !hasArticulationPoints*/) { + 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"); + exportPattern(danglingEdgesPath, patternGeometry, exportTilledPattern); } } else { continue; } } - if (hasFloatingComponents && !hasArticulationPoints && !tiledPatternHasDanglingEdges) { + if (hasFloatingComponents /*&& !hasArticulationPoints && !tiledPatternHasDanglingEdges */) { statistics.numberOfPatternsWithMoreThanASingleCC++; if (debugIsOn) { if (savePlyFiles) { @@ -599,34 +583,25 @@ void TopologyEnumerator::computeValidPatterns( assert(colorsRegistered == eCC[0].first); - tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath) - .append(patternName + "_tiled") - .string() - + ".ply"); + if (exportTilledPattern) { + tiledPatternGeometry.save(std::filesystem::path(moreThanOneCCPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } } } else { continue; } } - if (hasArticulationPoints && !hasFloatingComponents && !tiledPatternHasDanglingEdges) { + if (hasArticulationPoints /*&& !hasFloatingComponents && !tiledPatternHasDanglingEdges */) { 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; + exportPattern(articulationPointsPath, patternGeometry, exportTilledPattern); } } else { continue; @@ -649,39 +624,27 @@ void TopologyEnumerator::computeValidPatterns( 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"); + exportPattern(validPatternsPath, patternGeometry, exportTilledPattern); } - 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); + if (saveCompressedFormat) { + 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()) { + if (!patternSet.patterns.empty() && saveCompressedFormat) { PatternIO::save(compressedPatternsFilePath, patternSet); } diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp index a11855e..7871ce6 100755 --- a/topologyenumerator.hpp +++ b/topologyenumerator.hpp @@ -1,12 +1,12 @@ #ifndef TOPOLOGYENUMERATOR_HPP #define TOPOLOGYENUMERATOR_HPP +#include "nlohmann/json.hpp" #include "patternIO.hpp" #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" #include #include #include -//#include #include #include @@ -63,28 +63,26 @@ class TopologyEnumerator { size_t numberOfPatternsWithADanglingEdgeOrNode{0}; size_t numberOfPatternsWithArticulationPoints{0}; size_t numberOfValidPatterns{0}; - // nlohmann::json convertToJson() const { - // nlohmann::json json; - // json["numPossibleEdges"] = numberOfPossibleEdges; - // json["numCoincideEdges"] = numberOfCoincideEdges; - // json["numDuplicateEdges"] = numberOfDuplicateEdges; - // json["numValidEdges"] = numberOfValidEdges; - // json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs; - // json["numPatterns"] = numberOfPatterns; - // // json["numIntersectingEdgesOverAllPatterns"] = - // // numberOfIntersectingEdgesOverAllPatterns; - // json["numPatternsWithIntersectingEdges"] = - // numberOfPatternsWithIntersectingEdges; - // json["numPatternsWithNotASingleCC"] = - // numberOfPatternsWithMoreThanASingleCC; - // json["numPatternsWithDangling"] = - // numberOfPatternsWithADanglingEdgeOrNode; - // json["numPatternsWithArticulationPoints"] = - // numberOfPatternsWithArticulationPoints; - // json["numValidPatterns"] = numberOfValidPatterns; + nlohmann::json convertToJson() const + { + nlohmann::json json; + json["numPossibleEdges"] = numberOfPossibleEdges; + json["numCoincideEdges"] = numberOfCoincideEdges; + json["numDuplicateEdges"] = numberOfDuplicateEdges; + json["numValidEdges"] = numberOfValidEdges; + json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs; + json["numPatterns"] = numberOfPatterns; + // json["numIntersectingEdgesOverAllPatterns"] = + // numberOfIntersectingEdgesOverAllPatterns; + json["numPatternsWithIntersectingEdges"] = numberOfPatternsWithIntersectingEdges; + json["numPatternsWithNotASingleCC"] = numberOfPatternsWithMoreThanASingleCC; + json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; + json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints; + json["numValidPatterns"] = numberOfValidPatterns; + + return json; + } - // return json; - // } void print(const std::string &setupString, const std::filesystem::path &directoryPath) const { std::cout << "The setup " << setupString << std::endl; @@ -112,20 +110,37 @@ class TopologyEnumerator { << " patterns found with a dangling node or edge" << std::endl; std::cout << numberOfPatternsWithArticulationPoints << " patterns found with an articulation point" << std::endl; - std::cout << numberOfValidPatterns << " valid patterns were found" - << std::endl; - // if (!directoryPath.empty()) { - // auto json = convertToJson(); + std::cout << numberOfValidPatterns << " valid patterns were found" << std::endl; + if (!directoryPath.empty()) { + auto json = convertToJson(); - // std::ofstream file; - // file.open(std::filesystem::path(directoryPath) - // .append("statistics.csv") - // .string()); - // file << "setup," << setupString << "\n"; - // for (const auto &el : json.items()) { - // file << el.key() << "," << el.value() << "\n"; - // } - // } + std::ofstream file; + file.open(std::filesystem::path(directoryPath).append("statistics.csv").string()); + file << "setup," << setupString << "\n"; + for (const auto &el : json.items()) { + file << el.key() << ","; + } + file << "\n"; + for (const auto &el : json.items()) { + file << el.value() << ","; + } + file << "\n"; + } + } + + void reset() + { + numberOfPossibleEdges = 0; + numberOfCoincideEdges = 0; + numberOfDuplicateEdges = 0; + numberOfValidEdges = 0; + numberOfIntersectingEdgePairs = 0; + numberOfPatterns = 0; + numberOfPatternsWithIntersectingEdges = 0; + numberOfPatternsWithMoreThanASingleCC = 0; + numberOfPatternsWithADanglingEdgeOrNode = 0; + numberOfPatternsWithArticulationPoints = 0; + numberOfValidPatterns = 0; } }; @@ -173,8 +188,10 @@ private: const size_t &numberOfDesiredEdges, const std::filesystem::path &resultsPath, const std::vector &vertices, - const std::unordered_map> - &intersectingEdges, + const std::unordered_map> &intersectingEdges, const std::vector &validEdges); + void exportPattern(const std::filesystem::path &saveToPath, + PatternGeometry &patternGeometry, + const bool saveTilledPattern) const; }; #endif // TOPOLOGYENUMERATOR_HPP diff --git a/utilities.hpp b/utilities.hpp index 863768a..1d2dd6e 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -57,6 +57,15 @@ struct Vector6d : public std::array { return result; } + Vector6d operator/(const Vector6d &v) const + { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) / v[i]; + } + return result; + } + Vector6d operator+(const Vector6d &v) const { Vector6d result; for (size_t i = 0; i < 6; i++) {