From f8a98deb17af4f695fbe4cf92365a521146f83b2 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Mon, 15 Nov 2021 11:08:39 +0200 Subject: [PATCH] Refactoring --- CMakeLists.txt | 35 ++- csvfile.hpp | 3 +- drmsimulationmodel.cpp | 435 +++++++++++++++--------------- drmsimulationmodel.hpp | 12 +- edgemesh.cpp | 79 ++++-- edgemesh.hpp | 20 +- hexagonremesher.hpp | 128 +++++++++ linearsimulationmodel.hpp | 76 ++---- polymesh.hpp | 51 +++- reducedmodeloptimizer_structs.hpp | 254 +++++++++++------ simulation_structs.hpp | 156 ++++++++--- simulationmesh.cpp | 157 ++++++----- simulationmesh.hpp | 16 +- simulationmodel.hpp | 13 + topologyenumerator.cpp | 77 ++++-- topologyenumerator.hpp | 26 +- trianglepatterngeometry.cpp | 248 +++++++++-------- trianglepatterngeometry.hpp | 7 +- utilities.hpp | 14 + vcgtrimesh.cpp | 37 ++- vcgtrimesh.hpp | 1 + 21 files changed, 1191 insertions(+), 654 deletions(-) create mode 100644 hexagonremesher.hpp create mode 100644 simulationmodel.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 3bc7eec..596193f 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,32 @@ if(NOT EXISTS ${EXTERNAL_DEPS_DIR}) endif() ##Create directory for the external libraries file(MAKE_DIRECTORY ${EXTERNAL_DEPS_DIR}) +#message(${POLYSCOPE_ALREADY_COMPILED}) + +if(${USE_POLYSCOPE}) +download_project(PROJ POLYSCOPE + GIT_REPOSITORY https://github.com/nmwsharp/polyscope.git + GIT_TAG master + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +if(NOT EXISTS ${POLYSCOPE_BINARY_DIR}) +add_subdirectory(${POLYSCOPE_SOURCE_DIR} ${POLYSCOPE_BINARY_DIR}) +add_compile_definitions(POLYSCOPE_DEFINED) +endif() +endif() + +#dlib +set(DLIB_BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/dlib_bin) +file(MAKE_DIRECTORY ${DLIB_BIN_DIR}) +download_project(PROJ DLIB + GIT_REPOSITORY https://github.com/davisking/dlib.git + GIT_TAG master + BINARY_DIR ${DLIB_BIN_DIR} + PREFIX ${EXTERNAL_DEPS_DIR} + ${UPDATE_DISCONNECTED_IF_AVAILABLE} +) +add_subdirectory(${DLIB_SOURCE_DIR} ${DLIB_BINARY_DIR}) ##vcglib devel branch download_project(PROJ vcglib_devel @@ -55,7 +81,6 @@ option(TBB_BUILD_TESTS "Build TBB tests and enable testing infrastruct add_subdirectory(${TBB_SOURCE_DIR} ${TBB_BINARY_DIR}) link_directories(${TBB_BINARY_DIR}) - ###Eigen 3 NOTE: Eigen is required on the system the code is ran find_package(Eigen3 3.3 REQUIRED) @@ -75,9 +100,11 @@ target_include_directories(${PROJECT_NAME} PUBLIC ${MySourcesFiles} ) -if(${USE_POLYSCOPE}) - target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot polyscope glad ThreedBeamFEA tbb pthread) +if(${MYSOURCES_STATIC_LINK}) + target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot dlib::dlib ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) +elseif(${USE_POLYSCOPE}) + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread polyscope) else() - target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen matplot ThreedBeamFEA ${TBB_BINARY_DIR}/libtbb_static.a pthread) + target_link_libraries(${PROJECT_NAME} Eigen3::Eigen matplot dlib::dlib glad ThreedBeamFEA tbb pthread) endif() target_link_directories(MySources PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boost_graph/libs) diff --git a/csvfile.hpp b/csvfile.hpp index 7c9dffd..43d3c02 100755 --- a/csvfile.hpp +++ b/csvfile.hpp @@ -36,8 +36,7 @@ public: std::ofstream outfile(filename); outfile.close(); } - overwrite ? fs_.open(filename, std::ios::trunc) - : fs_.open(filename, std::ios::app); + overwrite ? fs_.open(filename, std::ios::trunc) : fs_.open(filename, std::ios::app); } } diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 03af7fe..cad27ec 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -218,22 +218,58 @@ void DRMSimulationModel::runUnitTests() // polyscope::show(); } -void DRMSimulationModel::reset() +void DRMSimulationModel::reset(const std::shared_ptr &pJob, const Settings &settings) { + mSettings = settings; mCurrentSimulationStep = 0; history.clear(); + history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); constrainedVertices.clear(); pMesh.reset(); + pMesh = std::make_unique(*pJob->pMesh); + vcg::tri::UpdateBounding::Box(*pMesh); plotYValues.clear(); plotHandle.reset(); checkedForMaximumMoment = false; externalMomentsNorm = 0; - Dt = mSettings.Dtini; + Dt = settings.Dtini; numOfDampings = 0; shouldTemporarilyDampForces = false; externalLoadStep = 1; isVertexConstrained.clear(); minTotalResidualForcesNorm = std::numeric_limits::max(); + + constrainedVertices = pJob->constrainedVertices; + if (!pJob->nodalForcedDisplacements.empty()) { + for (std::pair viDisplPair : pJob->nodalForcedDisplacements) { + const VertexIndex vi = viDisplPair.first; + constrainedVertices[vi].insert({0, 1, 2}); + } + } + computeRigidSupports(); + isVertexConstrained.resize(pMesh->VN(), false); + for (auto fixedVertex : constrainedVertices) { + isVertexConstrained[fixedVertex.first] = true; + } + +#ifdef POLYSCOPE_DEFINED + if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { + PolyscopeInterface::init(); + polyscope::registerCurveNetwork(meshPolyscopeLabel, + pMesh->getEigenVertices(), + pMesh->getEigenEdges()); + polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, + pMesh->getEigenVertices(), + pMesh->getEigenEdges()) + ->setRadius(0.002); + // registerWorldAxes(); + } +#endif + + if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { + shouldApplyInitialDistortion = true; + } + updateElementalFrames(); } VectorType DRMSimulationModel::computeDisplacementDifferenceDerivative( @@ -755,18 +791,18 @@ double DRMSimulationModel::computeTotalInternalPotentialEnergy() const double axialPE = pow(e_k, 2) * element.material.youngsModulus * element.A / (2 * element.initialLength); const double theta1Diff = theta1_jplus1 - theta1_j; - const double torsionalPE = element.G * element.J * pow(theta1Diff, 2) + const double torsionalPE = element.G * element.inertia.J * pow(theta1Diff, 2) / (2 * element.initialLength); const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) + 4 * theta2_j * theta2_jplus1 + 4 * pow(theta2_jplus1, 2); const double firstBendingPE = firstBendingPE_inBracketsTerm * element.material.youngsModulus - * element.I2 / (2 * element.initialLength); + * element.inertia.I2 / (2 * element.initialLength); const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) + 4 * theta3_j * theta3_jplus1 + 4 * pow(theta3_jplus1, 2); const double secondBendingPE = secondBendingPE_inBracketsTerm - * element.material.youngsModulus * element.I3 + * element.material.youngsModulus * element.inertia.I3 / (2 * element.initialLength); totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE; @@ -1220,11 +1256,12 @@ void DRMSimulationModel::updateNodalMasses() const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const long double SkRotational_I2 = elem.material.youngsModulus * elem.I2 + const long double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 / lengthToThe3; - const long double SkRotational_I3 = elem.material.youngsModulus * elem.I3 + const long double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 / lengthToThe3; - const long double SkRotational_J = elem.material.youngsModulus * elem.J / lengthToThe3; + const long double SkRotational_J = elem.material.youngsModulus * elem.inertia.J + / lengthToThe3; rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_J += SkRotational_J; @@ -1586,11 +1623,11 @@ void DRMSimulationModel::updatePositionsOnTheFly( const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const double SkRotational_I2 = elem.material.youngsModulus * elem.I2 + const double SkRotational_I2 = elem.material.youngsModulus * elem.inertia.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_I3 = elem.material.youngsModulus * elem.I3 + const double SkRotational_I3 = elem.material.youngsModulus * elem.inertia.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_J = elem.material.youngsModulus * elem.J + const double SkRotational_J = elem.material.youngsModulus * elem.inertia.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1686,9 +1723,129 @@ SimulationResults DRMSimulationModel::computeResults(const std::shared_ptrcurrentTotalKineticEnergy)) { + if (mSettings.beVerbose) { + std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; + } + results.converged = false; + } + // mesh.printVertexCoordinates(mesh.VN() / 2); +#ifdef POLYSCOPE_DEFINED + if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) { + draw(); + } +#endif + if (!std::isnan(pMesh->currentTotalKineticEnergy)) { + results.debug_drmDisplacements = results.displacements; + results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); + + results.rotationalDisplacementQuaternion.resize(pMesh->VN()); + results.debug_q_f1.resize(pMesh->VN()); + results.debug_q_normal.resize(pMesh->VN()); + results.debug_q_nr.resize(pMesh->VN()); + for (int vi = 0; vi < pMesh->VN(); vi++) { + const Node &node = pMesh->nodes[vi]; + const Eigen::Vector3d nInitial_eigen = node.initialNormal + .ToEigenVector(); + const Eigen::Vector3d nDeformed_eigen + = pMesh->vert[vi].cN().ToEigenVector(); + + Eigen::Quaternion q_normal; + q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); + Eigen::Quaternion q_nr_nDeformed; + q_nr_nDeformed = Eigen::AngleAxis(pMesh->nodes[vi].nR, nDeformed_eigen); + Eigen::Quaternion q_nr_nInit; + q_nr_nInit = Eigen::AngleAxis(pMesh->nodes[vi].nR, nInitial_eigen); + const auto w = q_nr_nDeformed.w(); + const auto theta = 2 * acos(q_nr_nDeformed.w()); + const auto nr = pMesh->nodes[vi].nR; + Eigen::Vector3d deformedNormal_debug(q_nr_nDeformed.x() * sin(theta / 2), + q_nr_nDeformed.y() * sin(theta / 2), + q_nr_nDeformed.z() * sin(theta / 2)); + deformedNormal_debug.normalize(); + // const double nr = nr_0To2pi <= M_PI ? nr_0To2pi : (nr_0To2pi - 2 * M_PI); + const double nr_debug = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? theta : -theta; + assert(pMesh->nodes[vi].nR - nr_debug < 1e-6); + VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1; + const VectorType &nDeformed = pMesh->vert[vi].cN(); + const VectorType referenceF1_deformed + = (referenceT1_deformed + - (node.initialNormal * (referenceT1_deformed * node.initialNormal))) + .Normalize(); + + const VectorType referenceT1_initial + = computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation, + pMesh->nodes[node.referenceElement->cV(1)].initialLocation); + // const VectorType &n_initial = node.initialNormal; + const VectorType referenceF1_initial = (referenceT1_initial + - (node.initialNormal + * (referenceT1_initial * node.initialNormal))) + .Normalize(); + Eigen::Quaternion q_f1_nInit; //nr is with respect to f1 + q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector(), + referenceF1_deformed.ToEigenVector()); + + Eigen::Quaternion q_f1_nDeformed; //nr is with respect to f1 + // const VectorType &n_initial = node.initialNormal; + const VectorType referenceF1_initial_def + = (referenceT1_initial - (nDeformed * (referenceT1_initial * nDeformed))).Normalize(); + const VectorType referenceF1_deformed_def = (referenceT1_deformed + - (nDeformed + * (referenceT1_deformed * nDeformed))) + .Normalize(); + q_f1_nDeformed + .setFromTwoVectors(referenceF1_initial_def.ToEigenVector(), + referenceF1_deformed_def.ToEigenVector()); + const auto debug_qf1_nDef = (q_f1_nDeformed * q_nr_nDeformed) * nDeformed_eigen; + const auto debug_qf1_nInit = (q_f1_nInit * q_nr_nInit) * nInitial_eigen; + const auto debug_deformedNormal_f1Init = (q_normal * (q_f1_nInit * q_nr_nInit)) + * nInitial_eigen; + const auto debug_deformedNormal_f1Def = ((q_nr_nDeformed * q_f1_nDeformed) * q_normal) + * nInitial_eigen; + // Eigen::Quaternion q_t1; + // q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), + // referenceT1_deformed.ToEigenVector()); + + results.debug_q_f1[vi] = q_f1_nInit; + results.debug_q_normal[vi] = q_normal; + results.debug_q_nr[vi] = q_nr_nInit; + results.rotationalDisplacementQuaternion[vi] + = (q_normal + * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; + //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]; + + // Eigen::Quaterniond q_test = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX()) + // * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY()) + // * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ()); + + // const double dot_test = results.rotationalDisplacementQuaternion[vi].dot(q_test); + // assert(dot_test > 1 - 1e5); + + int i = 0; + i++; + } + } + return results; } @@ -1844,9 +2001,9 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) for (const EdgeType &e : pMesh->edge) { const size_t ei = pMesh->getIndex(e); A[ei] = pMesh->elements[e].A; - J[ei] = pMesh->elements[e].J; - I2[ei] = pMesh->elements[e].I2; - I3[ei] = pMesh->elements[e].I3; + J[ei] = pMesh->elements[e].inertia.J; + I2[ei] = pMesh->elements[e].inertia.I2; + I3[ei] = pMesh->elements[e].inertia.I3; } polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); @@ -2105,93 +2262,9 @@ void DRMSimulationModel::applyKineticDamping(const std::shared_ptr &pJob, - const Settings &settings, - const SimulationResults &solutionGuess) +SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob) { - assert(pJob->pMesh.operator bool()); auto beginTime = std::chrono::high_resolution_clock::now(); - mSettings = settings; - reset(); - - history.label = pJob->pMesh->getLabel() + "_" + pJob->getLabel(); - - // if (!pJob->nodalExternalForces.empty()) { - // double externalForcesNorm = 0; - // for (const auto &externalForce : pJob->nodalExternalForces) { - // externalForcesNorm += externalForce.second.norm(); - // } - // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; - // } - - constrainedVertices = pJob->constrainedVertices; - if (!pJob->nodalForcedDisplacements.empty()) { - for (std::pair viDisplPair : pJob->nodalForcedDisplacements) { - const VertexIndex vi = viDisplPair.first; - constrainedVertices[vi].insert({0, 1, 2}); - } - } - // if (!pJob->nodalForcedNormals.empty()) { - // for (std::pair viNormalPair : - // pJob->nodalForcedDisplacements) { - // assert(viNormalPair3second[2] >= 0); - // } - // } - - pMesh.reset(); - pMesh = std::make_unique(*pJob->pMesh); - if (mSettings.beVerbose) { - std::cout << "Executing simulation for mesh with:" << pMesh->VN() << " nodes and " - << pMesh->EN() << " elements." << std::endl; - } - vcg::tri::UpdateBounding::Box(*pMesh); - computeRigidSupports(); - isVertexConstrained.resize(pMesh->VN(), false); - for (auto fixedVertex : constrainedVertices) { - isVertexConstrained[fixedVertex.first] = true; - } - -#ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw || mSettings.debugModeStep.has_value()) { - PolyscopeInterface::init(); - polyscope::registerCurveNetwork(meshPolyscopeLabel, - pMesh->getEigenVertices(), - pMesh->getEigenEdges()); - polyscope::registerCurveNetwork("Initial_" + meshPolyscopeLabel, - pMesh->getEigenVertices(), - pMesh->getEigenEdges()) - ->setRadius(0.002); - // registerWorldAxes(); - } -#endif - - if (!pJob->nodalForcedDisplacements.empty() && pJob->nodalExternalForces.empty()) { - shouldApplyInitialDistortion = true; - } - updateElementalFrames(); - if (!solutionGuess.displacements.empty()) { - assert(!mSettings.linearGuessForceScaleFactor.has_value()); - applySolutionGuess(solutionGuess, pJob); - shouldTemporarilyDampForces = true; - } - - if (mSettings.linearGuessForceScaleFactor.has_value()) { - const double forceScaleFactor = 1; - for (auto &externalForce : pJob->nodalExternalForces) { - externalForce.second = externalForce.second * forceScaleFactor; - } - LinearSimulationModel linearSimulationModel; - SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel - .executeSimulation(pJob); - // simulationResults_fullPatternLinearModel.save(fullResultsFolderPath); - for (auto &externalForce : pJob->nodalExternalForces) { - externalForce.second = externalForce.second / forceScaleFactor; - } - - applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); - shouldTemporarilyDampForces = true; - } - updateNodalMasses(); // std::unordered_map nodalExternalForces = pJob->nodalExternalForces; // double totalExternalForcesNorm = 0; @@ -2206,7 +2279,7 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces, constrainedVertices); if (!pJob->nodalExternalForces.empty()) { mSettings.totalResidualForcesNormThreshold - = settings.totalExternalForcesNormPercentageTermination + = mSettings.totalExternalForcesNormPercentageTermination * pMesh->totalExternalForcesNorm; } else { mSettings.totalResidualForcesNormThreshold = 1e-3; @@ -2220,13 +2293,18 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN() << " nodes and " + << pMesh->EN() << " elements." << std::endl; + } + const size_t movingAverageSampleSize = 200; std::queue residualForcesMovingAverageHistorySample; std::vector percentageOfConvergence; - // double residualForcesMovingAverageDerivativeMax = 0; - while (!settings.maxDRMIterations.has_value() - || mCurrentSimulationStep < settings.maxDRMIterations.value()) { + while (!mSettings.maxDRMIterations.has_value() + || mCurrentSimulationStep < mSettings.maxDRMIterations.value()) { if ((mSettings.debugModeStep.has_value() && mCurrentSimulationStep == 50000)) { // std::filesystem::create_directory("./PatternOptimizationNonConv"); // pJob->save("./PatternOptimizationNonConv"); @@ -2532,123 +2610,6 @@ mesh->currentTotalPotentialEnergykN*/ results.executionTime = std::chrono::duration_cast(endTime - beginTime) .count(); - if (settings.maxDRMIterations.has_value() && mCurrentSimulationStep == settings.maxDRMIterations - && mCurrentSimulationStep != 0) { - if (mSettings.beVerbose) { - std::cout << "Did not reach equilibrium before reaching the maximum number " - "of DRM steps (" - << settings.maxDRMIterations.value() << "). Breaking simulation" << std::endl; - } - results.converged = false; - } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { - if (mSettings.beVerbose) { - std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; - } - results.converged = false; - } - // mesh.printVertexCoordinates(mesh.VN() / 2); -#ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw && !mSettings.debugModeStep.has_value()) { - draw(); - } -#endif - if (!std::isnan(pMesh->currentTotalKineticEnergy)) { - results.debug_drmDisplacements = results.displacements; - results.internalPotentialEnergy = computeTotalInternalPotentialEnergy(); - - results.rotationalDisplacementQuaternion.resize(pMesh->VN()); - results.debug_q_f1.resize(pMesh->VN()); - results.debug_q_normal.resize(pMesh->VN()); - results.debug_q_nr.resize(pMesh->VN()); - for (int vi = 0; vi < pMesh->VN(); vi++) { - const Node &node = pMesh->nodes[vi]; - const Eigen::Vector3d nInitial_eigen = node.initialNormal - .ToEigenVector(); - const Eigen::Vector3d nDeformed_eigen - = pMesh->vert[vi].cN().ToEigenVector(); - - Eigen::Quaternion q_normal; - q_normal.setFromTwoVectors(nInitial_eigen, nDeformed_eigen); - Eigen::Quaternion q_nr_nDeformed; - q_nr_nDeformed = Eigen::AngleAxis(pMesh->nodes[vi].nR, nDeformed_eigen); - Eigen::Quaternion q_nr_nInit; - q_nr_nInit = Eigen::AngleAxis(pMesh->nodes[vi].nR, nInitial_eigen); - const auto w = q_nr_nDeformed.w(); - const auto theta = 2 * acos(q_nr_nDeformed.w()); - const auto nr = pMesh->nodes[vi].nR; - Eigen::Vector3d deformedNormal_debug(q_nr_nDeformed.x() * sin(theta / 2), - q_nr_nDeformed.y() * sin(theta / 2), - q_nr_nDeformed.z() * sin(theta / 2)); - deformedNormal_debug.normalize(); - // const double nr = nr_0To2pi <= M_PI ? nr_0To2pi : (nr_0To2pi - 2 * M_PI); - const double nr_debug = deformedNormal_debug.dot(nDeformed_eigen) > 0 ? theta : -theta; - assert(pMesh->nodes[vi].nR - nr_debug < 1e-6); - VectorType referenceT1_deformed = pMesh->elements[node.referenceElement].frame.t1; - const VectorType &nDeformed = pMesh->vert[vi].cN(); - const VectorType referenceF1_deformed - = (referenceT1_deformed - - (node.initialNormal * (referenceT1_deformed * node.initialNormal))) - .Normalize(); - - const VectorType referenceT1_initial - = computeT1Vector(pMesh->nodes[node.referenceElement->cV(0)].initialLocation, - pMesh->nodes[node.referenceElement->cV(1)].initialLocation); - // const VectorType &n_initial = node.initialNormal; - const VectorType referenceF1_initial = (referenceT1_initial - - (node.initialNormal - * (referenceT1_initial * node.initialNormal))) - .Normalize(); - Eigen::Quaternion q_f1_nInit; //nr is with respect to f1 - q_f1_nInit.setFromTwoVectors(referenceF1_initial.ToEigenVector(), - referenceF1_deformed.ToEigenVector()); - - Eigen::Quaternion q_f1_nDeformed; //nr is with respect to f1 - // const VectorType &n_initial = node.initialNormal; - const VectorType referenceF1_initial_def - = (referenceT1_initial - (nDeformed * (referenceT1_initial * nDeformed))).Normalize(); - const VectorType referenceF1_deformed_def = (referenceT1_deformed - - (nDeformed - * (referenceT1_deformed * nDeformed))) - .Normalize(); - q_f1_nDeformed - .setFromTwoVectors(referenceF1_initial_def.ToEigenVector(), - referenceF1_deformed_def.ToEigenVector()); - const auto debug_qf1_nDef = (q_f1_nDeformed * q_nr_nDeformed) * nDeformed_eigen; - const auto debug_qf1_nInit = (q_f1_nInit * q_nr_nInit) * nInitial_eigen; - const auto debug_deformedNormal_f1Init = (q_normal * (q_f1_nInit * q_nr_nInit)) - * nInitial_eigen; - const auto debug_deformedNormal_f1Def = ((q_nr_nDeformed * q_f1_nDeformed) * q_normal) - * nInitial_eigen; - // Eigen::Quaternion q_t1; - // q_t1.setFromTwoVectors(referenceT1_initial.ToEigenVector(), - // referenceT1_deformed.ToEigenVector()); - - results.debug_q_f1[vi] = q_f1_nInit; - results.debug_q_normal[vi] = q_normal; - results.debug_q_nr[vi] = q_nr_nInit; - results.rotationalDisplacementQuaternion[vi] - = (q_normal - * (q_f1_nInit * q_nr_nInit)); //q_f1_nDeformed * q_nr_nDeformed * q_normal; - //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]; - - // Eigen::Quaterniond q_test = Eigen::AngleAxisd(eulerAngles[0], Eigen::Vector3d::UnitX()) - // * Eigen::AngleAxisd(eulerAngles[1], Eigen::Vector3d::UnitY()) - // * Eigen::AngleAxisd(eulerAngles[2], Eigen::Vector3d::UnitZ()); - - // const double dot_test = results.rotationalDisplacementQuaternion[vi].dot(q_test); - // assert(dot_test > 1 - 1e5); - - int i = 0; - i++; - } - } - if (!mSettings.debugModeStep.has_value() && mSettings.shouldCreatePlots) { SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); @@ -2660,7 +2621,39 @@ mesh->currentTotalPotentialEnergykN*/ polyscope::removeCurveNetwork("Initial_" + meshPolyscopeLabel); } #endif + return results; +} +SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob, + const Settings &settings, + const SimulationResults &solutionGuess) +{ + reset(pJob, settings); + assert(pJob->pMesh != nullptr); + if (!solutionGuess.displacements.empty()) { + assert(!mSettings.linearGuessForceScaleFactor.has_value()); + applySolutionGuess(solutionGuess, pJob); + shouldTemporarilyDampForces = true; + } + + if (mSettings.linearGuessForceScaleFactor.has_value()) { + const double forceScaleFactor = mSettings.linearGuessForceScaleFactor.value(); + for (auto &externalForce : pJob->nodalExternalForces) { + externalForce.second = externalForce.second * forceScaleFactor; + } + LinearSimulationModel linearSimulationModel; + SimulationResults simulationResults_fullPatternLinearModel = linearSimulationModel + .executeSimulation(pJob); + // simulationResults_fullPatternLinearModel.save(fullResultsFolderPath); + for (auto &externalForce : pJob->nodalExternalForces) { + externalForce.second = externalForce.second / forceScaleFactor; + } + + applySolutionGuess(simulationResults_fullPatternLinearModel, pJob); + shouldTemporarilyDampForces = true; + } + + SimulationResults results = executeSimulation(pJob); //calling the virtual function return results; } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 31d31b0..e83c647 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -1,9 +1,9 @@ #ifndef BEAMFORMFINDER_HPP #define BEAMFORMFINDER_HPP -#include "simulationmesh.hpp" #include "matplot/matplot.h" -#include "simulation_structs.hpp" +#include "simulationmesh.hpp" +#include "simulationmodel.hpp" #include #include #include @@ -20,7 +20,7 @@ struct DifferentiateWithRespectTo { const DoFType &dofi; }; -class DRMSimulationModel +class DRMSimulationModel : public SimulationModel { public: struct Settings @@ -95,7 +95,7 @@ private: // std::unordered_map theta3Derivatives; bool shouldApplyInitialDistortion = false; - void reset(); + void reset(const std::shared_ptr &pJob, const Settings &settings); void updateNodalInternalForces( const std::unordered_map> &fixedVertices); void updateNodalExternalForces( @@ -239,10 +239,12 @@ private: void applyKineticDamping(const std::shared_ptr &pJob); + virtual SimulationResults executeSimulation(const std::shared_ptr &pJob); + public: DRMSimulationModel(); SimulationResults executeSimulation(const std::shared_ptr &pJob, - const Settings &settings = Settings(), + const Settings &settings, const SimulationResults &solutionGuess = SimulationResults()); static void runUnitTests(); diff --git a/edgemesh.cpp b/edgemesh.cpp index 87b87d7..b698300 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -1,9 +1,21 @@ #include "edgemesh.hpp" #include "vcg/simplex/face/topology.h" -#include +//#include +#include +#include Eigen::MatrixX2i VCGEdgeMesh::getEigenEdges() const { return eigenEdges; } +std::vector VCGEdgeMesh::getEdges() +{ + getEdges(eigenEdges); + std::vector edges(eigenEdges.rows()); + for (int ei = 0; ei < eigenEdges.rows(); ei++) { + edges[ei] = vcg::Point2i(eigenEdges(ei, 0), eigenEdges(ei, 1)); + } + return edges; +} + Eigen::MatrixX3d VCGEdgeMesh::getEigenVertices() { getVertices(eigenVertices); return eigenVertices; @@ -23,11 +35,12 @@ bool VCGEdgeMesh::save(const string &plyFilename) } 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; + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_EDGEINDEX; + mask |= vcg::tri::io::Mask::IOM_VERTNORMAL; mask |= vcg::tri::io::Mask::IOM_VERTCOLOR; - if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), *this, mask, false) != 0) { + // if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), *this, mask, false) != 0) { + if (vcg::tri::io::Exporter::Save(*this, filename.c_str(), mask) != 0) { return false; } return true; @@ -192,9 +205,13 @@ bool VCGEdgeMesh::load(const string &plyFilename) { } assert(std::filesystem::exists(usedPath)); this->Clear(); - if (!loadUsingNanoply(usedPath)) { - std::cerr << "Error: Unable to open " + usedPath << std::endl; - return false; + // if (!loadUsingNanoply(usedPath)) { + // std::cerr << "Error: Unable to open " + usedPath << std::endl; + // return false; + // } + if (!loadUsingDefaultLoader(usedPath)) { + std::cerr << "Error: Unable to open " + usedPath << std::endl; + return false; } getEdges(eigenEdges); @@ -212,21 +229,39 @@ bool VCGEdgeMesh::load(const string &plyFilename) { return true; } -bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) { +//bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) { - this->Clear(); - // assert(plyFileHasAllRequiredFields(plyFilename)); - // Load the ply file - unsigned int mask = 0; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; - mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; - mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; - if (nanoply::NanoPlyWrapper::LoadModel(plyFilename.c_str(), - *this, mask) != 0) { - return false; - } - return true; +// this->Clear(); +// // assert(plyFileHasAllRequiredFields(plyFilename)); +// // Load the ply file +// unsigned int mask = 0; +// mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; +// mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; +// mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; +// mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; +// if (nanoply::NanoPlyWrapper::LoadModel(plyFilename.c_str(), +// *this, mask) != 0) { +// return false; +// } +// return true; +//} + +bool VCGEdgeMesh::loadUsingDefaultLoader(const std::string &plyFilename) +{ + this->Clear(); + // assert(plyFileHasAllRequiredFields(plyFilename)); + // Load the ply file + int mask = 0; + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_VERTNORMAL; + mask |= vcg::tri::io::Mask::IOM_VERTCOLOR; + mask |= vcg::tri::io::Mask::IOM_EDGEINDEX; + // if (nanoply::NanoPlyWrapper::LoadModel(plyFilename.c_str(), + // *this, mask) != 0) { + if (vcg::tri::io::Importer::Open(*this, plyFilename.c_str(), mask) != 0) { + return false; + } + return true; } // bool VCGEdgeMesh::plyFileHasAllRequiredFields(const std::string diff --git a/edgemesh.hpp b/edgemesh.hpp index 65ef25f..1379318 100755 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -68,7 +68,7 @@ public: 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; @@ -78,13 +78,15 @@ public: void createSpiral(const float °reesOfArm, const size_t &numberOfSamples); Eigen::MatrixX2i getEigenEdges() const; + std::vector getEdges(); + 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 double &desiredRadius = 0.002, + const double &desiredRadius = 0.001, const bool &shouldEnable = true); void unregister() const; void drawInitialFrames(polyscope::CurveNetwork *polyscopeHandle_initialMesh) const; @@ -93,10 +95,24 @@ public: vcg::tri::Allocator::PointerUpdater &pu_vertices, vcg::tri::Allocator::PointerUpdater &pu_edges); + void centerMesh() + { + CoordType centerOfMass(0, 0, 0); + + for (const auto &v : vert) { + centerOfMass = centerOfMass + v.cP(); + } + centerOfMass = centerOfMass / VN(); + for (auto &v : vert) { + v.P() = v.cP() - centerOfMass; + } + } + private: void GeneratedRegularSquaredPattern(const double angleDeg, std::vector> &pattern, const size_t &desiredNumberOfSamples); + bool loadUsingDefaultLoader(const std::string &plyFilename); }; using VectorType = VCGEdgeMesh::CoordType; diff --git a/hexagonremesher.hpp b/hexagonremesher.hpp new file mode 100644 index 0000000..6d76967 --- /dev/null +++ b/hexagonremesher.hpp @@ -0,0 +1,128 @@ +#ifndef HEXAGONREMESHER_HPP +#define HEXAGONREMESHER_HPP +#include "polymesh.hpp" +#include "vcgtrimesh.hpp" +#include +#include +#include +#include +#include +#include + +namespace PolygonalRemeshing { +//std::shared_ptr +std::shared_ptr remeshWithPolygons(VCGTriMesh &surface) +{ + // const std::string tileIntoFilePath( + // "/home/iason/Coding/build/PatternTillingReducedModel/RelWithDebInfo/plane.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/hexagon.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/instantMeshes_plane_34.ply"); + // const std::string tileIntoFilePath("/home/iason/Coding/build/PatternTillingReducedModel/" + // "RelWithDebInfo/instantMeshes_plane.ply"); + // assert(std::filesystem::exists(tileIntoFilePath)); + // assert(tileInto_polyMesh.load(tileIntoFilePath)); + + ///Load triangle mesh + // VCGTriMesh tileInto_triMesh; + // assert(std::filesystem::exists(tileIntoFilePath)); + // tileInto_triMesh.load(tileIntoFilePath); + + ///Sample the triangle mesh + // std::vector pointVec; + // double radius = 1; + // int sampleNum = 50; + // // // vcg::tri::PoissonSampling(tileInto_triMesh, pointVec, sampleNum, radius); + // // // vcg::tri::TrivialSampler ps; + // // // vcg::tri::SurfaceSampling::VertexUniform(tileInto_triMesh, ps, sampleNum); + // // // pointVec = ps.SampleVec(); + // vcg::tri::MontecarloSampling(tileInto_triMesh, pointVec, sampleNum); + + // vcg::tri::VoronoiProcessingParameter vpp; + // vcg::tri::VoronoiProcessing::PreprocessForVoronoi(tileInto_triMesh, radius, vpp); + // std::vector seedVec; + // vcg::tri::VoronoiProcessing::SeedToVertexConversion(tileInto_triMesh, + // pointVec, + // seedVec); + + // vcg::tri::EuclideanDistance df; + // vpp.geodesicRelaxFlag = false; + // const int iterNum = 100; + // int actualIter = vcg::tri::VoronoiProcessing::VoronoiRelaxing(tileInto_triMesh, + // seedVec, + // iterNum, + // df, + // vpp); + // std::cout << "Voronoi relaxation iterations performed:" << actualIter << std::endl; + + // bool saveWasSuccessful + // = 0 + // == vcg::tri::io::ExporterOBJ::Save(tileInto_triMesh, + // "./tileInto_triMesh.obj", + // vcg::tri::io::Mask::IOM_ALL); + + // VCGTriMesh delaunay; + // vcg::tri::VoronoiProcessing::ConvertDelaunayTriangulationToMesh(tileInto_triMesh, + // delaunay, + // seedVec); + // delaunay.setLabel("Delaunay"); + // delaunay.save(delaunay.getLabel() + ".ply"); + // delaunay.registerForDrawing(); + + // ////Load surface + VCGPolyMesh tileInto_polyMesh, dual; + // vcg::tri::PolygonSupport::ImportFromTriMesh(tileInto_polyMesh, + // delaunay); + vcg::tri::PolygonSupport::ImportFromTriMesh(tileInto_polyMesh, surface); + vcg::tri::DualMeshing::MakeDual(tileInto_polyMesh, dual, false); + + // for(size_t i=0; i< Edges.size(); ++i) + // { + // std::vector fpVec; + // std::vector eiVec; + // face::EFStarFF(Edges[i].f,Edges[i].z,fpVec,eiVec); + // for(size_t j=0;jFEp(eiVec[j])=&(m.edge[i]); + // } + dual.setLabel(surface.getLabel() + "_polyDual"); + + // bool saveWasSuccessful + // = 0 + // == vcg::tri::io::ExporterOBJ::Save(dual, + // (dual.getLabel() + ".obj").c_str(), + // vcg::tri::io::Mask::IOM_BITPOLYGONAL); + // assert(saveWasSuccessful); + vcg::tri::UpdateNormal::PerFaceNormalized(dual); + + // vcg::PolygonalAlgorithm::SmoothReprojectPCA(dual, 1000, false, 0.5, 0, false, false); + // vcg::tri::io::ExporterOBJ::Save(dual, + // "./dual_optimized.obj", + // vcg::tri::io::Mask::IOM_BITPOLYGONAL); + + // if (vcg::tri::VoronoiProcessing::CheckVoronoiTopology(tileInto_triMesh, seedVec)) { + // VCGTriMesh voronoiMesh, voronoiPoly; + // vcg::tri::VoronoiProcessing::ConvertVoronoiDiagramToMesh(tileInto_triMesh, + // voronoiMesh, + // voronoiPoly, + // seedVec, + // vpp); + // vcg::tri::Allocator::CompactEveryVector(voronoiMesh); + // voronoiMesh.setLabel("Voro"); + // voronoiMesh.save(voronoiMesh.getLabel() + ".ply"); + // // voronoiMesh.registerForDrawing(); + // vcg::tri::Allocator::CompactEveryVector(voronoiPoly); + // voronoiPoly.setLabel("Poly"); + // voronoiPoly.save(voronoiPoly.getLabel() + ".ply"); + // // voronoiPoly.registerForDrawing(); + // // polyscope::show(); + // } + // dual.registerForDrawing(); + // polyscope::show(); + // assert(vcg::tri::IsValidPointer(dual, dual.face[0].cFEp(0))); + return std::make_shared(dual); +} + +} // namespace PolygonalRemeshing + +#endif // HEXAGONREMESHER_HPP diff --git a/linearsimulationmodel.hpp b/linearsimulationmodel.hpp index 1ca6621..d2c544f 100644 --- a/linearsimulationmodel.hpp +++ b/linearsimulationmodel.hpp @@ -1,55 +1,13 @@ #ifndef LINEARSIMULATIONMODEL_HPP #define LINEARSIMULATIONMODEL_HPP -//#include "beam.hpp" -#include "simulation_structs.hpp" +#include "simulationmodel.hpp" #include "threed_beam_fea.h" #include #include -// struct BeamSimulationProperties { -// float crossArea; -// float I2; -// float I3; -// float polarInertia; -// float G; -// // Properties used by fea -// float EA; -// float EIz; -// float EIy; -// float GJ; - -// BeamSimulationProperties(const BeamDimensions &dimensions, -// const BeamMaterial &material); -//}; - -// struct NodalForce { -// int index; -// int dof; -// double magnitude; -//}; - -// struct SimulationJob { -// Eigen::MatrixX3d nodes; -// Eigen::MatrixX2i elements; -// Eigen::MatrixX3d elementalNormals; -// Eigen::VectorXi fixedNodes; -// std::vector nodalForces; -// std::vector beamDimensions; -// std::vector beamMaterial; -//}; - -// struct SimulationResults { -// std::vector edgeForces; ///< Force values per force -// component -// ///< #force components x #edges -// Eigen::MatrixXd -// nodalDisplacements; ///< The displacement of each node #nodes x 3 -// SimulationResults(const fea::Summary &feaSummary); -// SimulationResults() {} -//}; - -class LinearSimulationModel { +class LinearSimulationModel : public SimulationModel +{ public: LinearSimulationModel(){ @@ -67,8 +25,11 @@ public: (evn0[2] + evn1[2]) / 2}; const Element &element = job->pMesh->elements[edgeIndex]; const double E = element.material.youngsModulus; - fea::Props feaProperties(E * element.A, E * element.I3, E * element.I2, - element.G * element.J, nAverageVector); + fea::Props feaProperties(E * element.A, + E * element.inertia.I3, + E * element.inertia.I2, + element.G * element.inertia.J, + nAverageVector); const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0)); const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1)); elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties); @@ -91,8 +52,8 @@ public: static std::vector getFeaFixedNodes(const std::shared_ptr &job) { std::vector boundaryConditions; - boundaryConditions.reserve(job->constrainedVertices.size() * 6 - + job->nodalForcedDisplacements.size() * 3); + // boundaryConditions.reserve(job->constrainedVertices.size() * 6 + // + job->nodalForcedDisplacements.size() * 3); for (auto fixedVertex : job->constrainedVertices) { const int vertexIndex = fixedVertex.first; for (int dofIndex : fixedVertex.second) { @@ -157,6 +118,7 @@ public: // results.rotationalDisplacementQuaternion[vi] = q_nz * q_ny * q_nx; // results.rotationalDisplacementQuaternion[vi] = q_nz * q_nx * q_ny; } + results.setLabelPrefix("Linear"); // // Convert forces // // Convert to vector of eigen matrices of the form force component-> per @@ -212,19 +174,19 @@ public: const double elementPotentialEnergy_bendingY_v0 = std::pow(elementForces[4], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.I2); + * element.inertia.I2); const double elementPotentialEnergy_bendingZ_v0 = std::pow(elementForces[5], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.I3); + * element.inertia.I3); const double elementPotentialEnergy_bendingY_v1 = std::pow(elementForces[10], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.I2); + * element.inertia.I2); const double elementPotentialEnergy_bendingZ_v1 = std::pow(elementForces[11], 2) * element.initialLength / (4 * element.material.youngsModulus - * element.I3); + * element.inertia.I3); const double elementPotentialEnergy_bending = elementPotentialEnergy_bendingY_v0 + elementPotentialEnergy_bendingZ_v0 + elementPotentialEnergy_bendingY_v1 @@ -232,10 +194,10 @@ public: //Torsion const double elementPotentialEnergy_torsion_v0 = std::pow(elementForces[3], 2) * element.initialLength - / (4 * element.J * element.G); + / (4 * element.inertia.J * element.G); const double elementPotentialEnergy_torsion_v1 = std::pow(elementForces[9], 2) * element.initialLength - / (4 * element.J * element.G); + / (4 * element.inertia.J * element.G); const double elementPotentialEnergy_torsion = elementPotentialEnergy_torsion_v0 + elementPotentialEnergy_torsion_v1; const double elementInternalPotentialEnergy = elementPotentialEnergy_axial @@ -244,7 +206,7 @@ public: + elementPotentialEnergy_torsion; results.internalPotentialEnergy += elementInternalPotentialEnergy; } - results.job = simulationJob; + results.pJob = simulationJob; return results; } @@ -285,7 +247,7 @@ public: fea::solve(job, fixedVertices, nodalForces, ties, equations, opts); SimulationResults results = getResults(feaResults, simulationJob); - results.job = simulationJob; + results.pJob = simulationJob; return results; } // SimulationResults getResults() const; diff --git a/polymesh.hpp b/polymesh.hpp index 8642dd7..aa50d47 100755 --- a/polymesh.hpp +++ b/polymesh.hpp @@ -3,7 +3,8 @@ #include "mesh.hpp" #include "vcg/complex/complex.h" #include -#include +#include +#include #ifdef POLYSCOPE_DEFINED #include @@ -54,7 +55,8 @@ class PFace : public vcg::Face + , + vcg::face::Color4b> {}; class VCGPolyMesh @@ -110,6 +112,49 @@ public: return faces; } + + bool saveOBJ(const std::filesystem::path &objFilePath = std::filesystem::path()) + { + const std::string extension = ".obj"; + std::filesystem::path filePath = objFilePath; + if (filePath.empty()) { + filePath = std::filesystem::current_path().append(getLabel() + extension).string(); + } else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) { + filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string(); + } + assert(std::filesystem::path(filePath).extension().string() == extension); + unsigned int mask = 0; + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_VERTNORMAL; + mask |= vcg::tri::io::Mask::IOM_FACEINDEX; + mask |= vcg::tri::io::Mask::IOM_FACECOLOR; + if (vcg::tri::io::ExporterOBJ::Save(*this, filePath.c_str(), mask) != 0) { + return false; + } + return true; + } + + bool savePLY(const std::filesystem::path &objFilePath = std::filesystem::path()) + { + const std::string extension = ".ply"; + std::filesystem::path filePath = objFilePath; + if (filePath.empty()) { + filePath = std::filesystem::current_path().append(getLabel() + extension).string(); + } else if (std::filesystem::is_directory(std::filesystem::path(objFilePath))) { + filePath = std::filesystem::path(objFilePath).append(getLabel() + extension).string(); + } + assert(std::filesystem::path(filePath).extension().string() == extension); + unsigned int mask = 0; + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_VERTNORMAL; + mask |= vcg::tri::io::Mask::IOM_FACEINDEX; + mask |= vcg::tri::io::Mask::IOM_FACECOLOR; + if (vcg::tri::io::ExporterPLY::Save(*this, filePath.c_str(), mask, false) != 0) { + return false; + } + return true; + } + #ifdef POLYSCOPE_DEFINED polyscope::SurfaceMesh *registerForDrawing( const std::optional> &desiredColor = std::nullopt, @@ -131,7 +176,7 @@ public: const glm::vec3 desiredColor_glm(desiredColor.value()[0], desiredColor.value()[1], desiredColor.value()[2]); - polyscopeHandle_mesh->setSurfaceColor(desiredColor_glm); + polyscopeHandle_mesh->setSurfaceColor(glm::normalize(desiredColor_glm)); } return polyscopeHandle_mesh; diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index 6d42852..84dc216 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -25,19 +25,23 @@ inline static std::vector baseSimulationScenarioNames{"Axial", struct Colors { - inline static std::array fullInitial{0.416666, 0.109804, 0.890196}; - inline static std::array fullDeformed{0.583333, 0.890196, 0.109804}; + inline static std::array fullInitial{0.094, 0.094, 0.094}; + // inline static std::array fullDeformed{0.583333, 0.890196, 0.109804}; + inline static std::array fullDeformed{0.094, 0.094, 0.094}; inline static std::array reducedInitial{0.890196, 0.109804, 0.193138}; - inline static std::array reducedDeformed{0.109804, 0.890196, 0.806863}; + inline static std::array reducedDeformed{0.262, 0.627, 0.910}; }; - struct xRange{ - std::string label; +struct xRange +{ + std::string label; double min; double max; - std::string toString() const { - return label + "=[" + std::to_string(min) + "," + std::to_string(max) + - "]"; + + inline bool operator<(const xRange &other) { return label < other.label; } + std::string toString() const + { + return label + "=[" + std::to_string(min) + "," + std::to_string(max) + "]"; } void fromString(const std::string &s) @@ -56,7 +60,8 @@ struct Colors { return label == xrange.label && min == xrange.min && max == xrange.max; } - }; +}; + struct Settings { enum NormalizationStrategy { NonNormalized, Epsilon }; @@ -65,7 +70,9 @@ struct Colors int numberOfFunctionCalls{100000}; double solverAccuracy{1e-2}; NormalizationStrategy normalizationStrategy{Epsilon}; - double normalizationParameter{0.003}; + double translationNormalizationParameter{0.003}; + double rotationNormalizationParameter{3}; + bool splitGeometryMaterialOptimization{false}; struct ObjectiveWeights { double translational{1}; @@ -168,6 +175,7 @@ struct Colors os << "Normalization strategy"; os << "Trans weight"; os << "Rot weight"; + os << "Splitted geo from mat opt"; // os << std::endl; } @@ -182,9 +190,10 @@ struct Colors os << numberOfFunctionCalls; os << solverAccuracy; os << normalizationStrategyStrings[normalizationStrategy] + "_" - + std::to_string(normalizationParameter); + + std::to_string(translationNormalizationParameter); os << objectiveWeights.translational; os << objectiveWeights.rotational; + os << (splitGeometryMaterialOptimization == true ? "true" : "false"); } }; @@ -193,8 +202,9 @@ struct Colors return settings1.numberOfFunctionCalls == settings2.numberOfFunctionCalls && settings1.xRanges == settings2.xRanges && settings1.solverAccuracy == settings2.solverAccuracy - && settings1.objectiveWeights.translational - == settings2.objectiveWeights.translational; + && settings1.objectiveWeights.translational == settings2.objectiveWeights.translational + && settings1.translationNormalizationParameter + == settings2.translationNormalizationParameter; } inline void updateMeshWithOptimalVariables(const std::vector &x, SimulationMesh &m) @@ -212,9 +222,9 @@ struct Colors Element &e = m.elements[ei]; e.setDimensions(CrossSectionType(beamWidth, beamHeight)); e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); - e.J = J; - e.I2 = I2; - e.I3 = I3; + e.inertia.J = J; + e.inertia.I2 = I2; + e.inertia.I3 = I3; } CoordType center_barycentric(1, 0, 0); @@ -257,6 +267,16 @@ struct Colors bool wasSuccessful{true}; // int numberOfSimulationCrashes{0}; Settings settings; + + std::vector> optimalXNameValuePairs; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + + PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern + vcg::Triangle3 baseTriangle; + std::string notes; + + //Data gathered for csv exporting struct ObjectiveValues { double totalRaw; @@ -268,15 +288,18 @@ struct Colors std::vector perSimulationScenario_total; } objectiveValue; std::vector perScenario_fullPatternPotentialEnergy; - // std::vector> optimalXNameValuePairs; - std::vector> optimalXNameValuePairs; - std::vector> fullPatternSimulationJobs; - std::vector> reducedPatternSimulationJobs; + std::vector objectiveValueHistory; + std::vector objectiveValueHistory_iteration; - PatternGeometry baseTriangleFullPattern; //non-fanned,non-tiled full pattern - vcg::Triangle3 baseTriangle; - std::string notes; + struct CSVExportingSettings + { + bool exportPE{false}; + bool exportIterationOfMinima{true}; + bool exportRawObjectiveValue{false}; + CSVExportingSettings() {} + }; + const CSVExportingSettings exportSettings; struct JsonKeys { inline static std::string filename{"OptimizationResults.json"}; @@ -290,10 +313,14 @@ struct Colors void save(const string &saveToPath, const bool shouldExportDebugFiles = false) { - std::filesystem::create_directories(saveToPath); //clear directory - for (const auto &entry : std::filesystem::directory_iterator(saveToPath)) - std::filesystem::remove_all(entry.path()); + if (std::filesystem::exists(saveToPath)) { + for (const auto &entry : std::filesystem::directory_iterator(saveToPath)) { + std::error_code ec; + std::filesystem::remove_all(entry.path(), ec); + } + } + std::filesystem::create_directories(saveToPath); //Save optimal X nlohmann::json json_optimizationResults; @@ -373,6 +400,12 @@ struct Colors } } + //save minima info + std::filesystem::path csvFilepathMinimaInfo = std::filesystem::path(saveToPath) + .append("minimaInfo.csv"); + csvFile csv_minimaInfo(csvFilepathMinimaInfo, false); + writeMinimaInfoTo(csv_minimaInfo); + settings.save(saveToPath); #ifdef POLYSCOPE_DEFINED @@ -460,7 +493,7 @@ struct Colors if (filepath.extension() == ".json") { SimulationJob job; job.load(filepath.string()); - applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh); + // applyOptimizationResults_innerHexagon(*this, baseTriangle, *job.pMesh); applyOptimizationResults_elements(*this, job.pMesh); reducedPatternSimulationJobs.push_back( std::make_shared(job)); @@ -511,6 +544,15 @@ struct Colors vcg::math::ToRad(hexAngle)) * reducedPattern.vert[0].cP(); } + // for (int rotationCounter = 0; + // rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { + // pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + // = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, + // vcg::math::ToRad(60.0 * rotationCounter)) + // * baseTriangleMovableVertexPosition; + // } + // reducedPattern.registerForDrawing(); + // polyscope::show(); // CoordType p0 = reducedPattern.vert[0].P(); // CoordType p1 = reducedPattern.vert[1].P(); // int i = 0; @@ -528,36 +570,53 @@ struct Colors reducedPattern_optimizationResults.optimalXNameValuePairs.end()); const std::string ALabel = "A"; - assert(optimalXVariables.contains(ALabel)); - const double A = optimalXVariables.at(ALabel); - const double beamWidth = std::sqrt(A); - const double beamHeight = beamWidth; - CrossSectionType elementDimensions(beamWidth, beamHeight); + if (optimalXVariables.contains(ALabel)) { + const double A = optimalXVariables.at(ALabel); + const double beamWidth = std::sqrt(A); + const double beamHeight = beamWidth; + CrossSectionType elementDimensions(beamWidth, beamHeight); + for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.setDimensions(elementDimensions); + } + } const double poissonsRatio = 0.3; const std::string ymLabel = "E"; - assert(optimalXVariables.contains(ymLabel)); - const double E = optimalXVariables.at(ymLabel); - const ElementMaterial elementMaterial(poissonsRatio, E); + if (optimalXVariables.contains(ymLabel)) { + const double E = optimalXVariables.at(ymLabel); + const ElementMaterial elementMaterial(poissonsRatio, E); + for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.setMaterial(elementMaterial); + } + } const std::string JLabel = "J"; - assert(optimalXVariables.contains(JLabel)); - const double J = optimalXVariables.at(JLabel); + if (optimalXVariables.contains(JLabel)) { + const double J = optimalXVariables.at(JLabel); + for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.inertia.J = J; + } + } const std::string I2Label = "I2"; - assert(optimalXVariables.contains(I2Label)); - const double I2 = optimalXVariables.at(I2Label); + if (optimalXVariables.contains(I2Label)) { + const double I2 = optimalXVariables.at(I2Label); + for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.inertia.I2 = I2; + } + } const std::string I3Label = "I3"; - assert(optimalXVariables.contains(I3Label)); - const double I3 = optimalXVariables.at(I3Label); - for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { - Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; - e.setDimensions(elementDimensions); - e.setMaterial(elementMaterial); - e.J = J; - e.I2 = I2; - e.I3 = I3; + if (optimalXVariables.contains(I3Label)) { + const double I3 = optimalXVariables.at(I3Label); + for (int ei = 0; ei < pTiledReducedPattern_simulationMesh->EN(); ei++) { + Element &e = pTiledReducedPattern_simulationMesh->elements[ei]; + e.inertia.I3 = I3; + } } pTiledReducedPattern_simulationMesh->reset(); } @@ -579,9 +638,10 @@ struct Colors = fullPatternSimulationJobs[simulationJobIndex]; pFullPatternSimulationJob->registerForDrawing( fullPatternSimulationJobs[0]->pMesh->getLabel()); - SimulationResults fullModelResults = drmSimulator.executeSimulation( - pFullPatternSimulationJob); - fullModelResults.registerForDrawing(Colors::fullDeformed, true, true); + DRMSimulationModel::Settings drmSettings; + SimulationResults fullModelResults + = drmSimulator.executeSimulation(pFullPatternSimulationJob, drmSettings); + fullModelResults.registerForDrawing(Colors::fullDeformed, true); // SimulationResults fullModelLinearResults = // linearSimulator.executeSimulation(pFullPatternSimulationJob); // fullModelLinearResults.setLabelPrefix("linear"); @@ -589,15 +649,20 @@ struct Colors // Drawing of reduced pattern results const std::shared_ptr &pReducedPatternSimulationJob = reducedPatternSimulationJobs[simulationJobIndex]; - // SimulationResults reducedModelResults = drmSimulator.executeSimulation( - // pReducedPatternSimulationJob); + // pReducedPatternSimulationJob->pMesh->registerForDrawing(); + // polyscope::show(); + // SimulationResults reducedModelResults = drmSimulator.executeSimulation( + // pReducedPatternSimulationJob); // reducedModelResults.registerForDrawing(Colors::reducedDeformed, false); - SimulationResults reducedModelLinearResults = linearSimulator.executeSimulation( + // SimulationResults reducedModelResults + // = drmSimulator.executeSimulation(pReducedPatternSimulationJob, + // DRMSimulationModel::Settings()); + SimulationResults reducedModelResults = linearSimulator.executeSimulation( pReducedPatternSimulationJob); - reducedModelLinearResults.setLabelPrefix("linear"); - reducedModelLinearResults.registerForDrawing(Colors::reducedDeformed, false); + reducedModelResults.setLabelPrefix("linear"); + reducedModelResults.registerForDrawing(Colors::reducedDeformed, true); polyscope::options::programName = fullPatternSimulationJobs[0]->pMesh->getLabel(); - polyscope::view::resetCameraToHomeView(); + // polyscope::view::resetCameraToHomeView(); polyscope::show(); // Save a screensh const std::string screenshotFilename @@ -606,18 +671,16 @@ struct Colors + fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" + pFullPatternSimulationJob->getLabel(); polyscope::screenshot(screenshotFilename, false); + pFullPatternSimulationJob->unregister(fullPatternSimulationJobs[0]->pMesh->getLabel()); fullModelResults.unregister(); // reducedModelResults.unregister(); - reducedModelLinearResults.unregister(); + reducedModelResults.unregister(); // fullModelLinearResults.unregister(); - // double error = computeError( - // reducedModelResults.displacements,fullModelResults.displacements, - // ); - // std::cout << "Error of simulation scenario " - // << - // simula simulationScenarioStrings[simulationScenarioIndex] - // << " is " - // << error << std::endl; + // double error = computeError(reducedModelLinearResults.displacements, + // fullModelResults.displacements); + // std::cout << "Error of simulation scenario " + // << simulationScenarioStrings[simulationScenarioIndex] << " is " << error + // << std::endl; } } #endif // POLYSCOPE_DEFINED @@ -638,8 +701,9 @@ struct Colors // Drawing of full pattern results const std::shared_ptr &pFullPatternSimulationJob = fullPatternSimulationJobs[simulationJobIndex]; - SimulationResults fullModelResults = simulator_drm.executeSimulation( - pFullPatternSimulationJob); + DRMSimulationModel::Settings drmSettings; + SimulationResults fullModelResults + = simulator_drm.executeSimulation(pFullPatternSimulationJob, drmSettings); fullModelResults.saveDeformedModel(); // Drawing of reduced pattern results @@ -653,8 +717,13 @@ struct Colors void writeHeaderTo(csvFile &os) { - os << "Total raw Obj value"; + if (exportSettings.exportRawObjectiveValue) { + os << "Total raw Obj value"; + } os << "Total Obj value"; + if (exportSettings.exportIterationOfMinima) { + os << "Iteration of minima"; + } for (int simulationScenarioIndex = 0; simulationScenarioIndex < fullPatternSimulationJobs.size(); simulationScenarioIndex++) { @@ -664,12 +733,15 @@ struct Colors os << "Obj value Rot " + simulationScenarioName; os << "Obj value Total " + simulationScenarioName; } - for (int simulationScenarioIndex = 0; - simulationScenarioIndex < fullPatternSimulationJobs.size(); - simulationScenarioIndex++) { - const std::string simulationScenarioName - = fullPatternSimulationJobs[simulationScenarioIndex]->getLabel(); - os << "PE " + simulationScenarioName; + + if (exportSettings.exportPE) { + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < fullPatternSimulationJobs.size(); + simulationScenarioIndex++) { + const std::string simulationScenarioName + = fullPatternSimulationJobs[simulationScenarioIndex]->getLabel(); + os << "PE " + simulationScenarioName; + } } for (const auto &nameValuePair : optimalXNameValuePairs) { os << nameValuePair.first; @@ -679,10 +751,15 @@ struct Colors // os << "notes"; } - void writeResultsTo(const Settings &settings_optimization, csvFile &os) const + void writeResultsTo(csvFile &os) const { - os << objectiveValue.totalRaw; + if (exportSettings.exportRawObjectiveValue) { + os << objectiveValue.totalRaw; + } os << objectiveValue.total; + if (exportSettings.exportIterationOfMinima) { + os << objectiveValueHistory_iteration.back(); + } for (int simulationScenarioIndex = 0; simulationScenarioIndex < fullPatternSimulationJobs.size(); simulationScenarioIndex++) { @@ -690,10 +767,12 @@ struct Colors os << objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex]; os << objectiveValue.perSimulationScenario_total[simulationScenarioIndex]; } - for (int simulationScenarioIndex = 0; - simulationScenarioIndex < fullPatternSimulationJobs.size(); - simulationScenarioIndex++) { - os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex]; + if (exportSettings.exportPE) { + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < fullPatternSimulationJobs.size(); + simulationScenarioIndex++) { + os << perScenario_fullPatternPotentialEnergy[simulationScenarioIndex]; + } } for (const auto &optimalXNameValuePair : optimalXNameValuePairs) { os << optimalXNameValuePair.second; @@ -707,7 +786,20 @@ struct Colors // } // os << notes; } + + void writeMinimaInfoTo(csvFile &outputCsv) + { + outputCsv << "Iteration"; + outputCsv << "Objective value"; + for (int objectiveValueIndex = 0; objectiveValueIndex < objectiveValueHistory.size(); + objectiveValueIndex++) { + outputCsv.endrow(); + outputCsv << objectiveValueHistory_iteration[objectiveValueIndex]; + outputCsv << objectiveValueHistory[objectiveValueIndex]; + } + } }; + enum SimulationModelTag { DRM, Linear }; } // namespace ReducedPatternOptimization #endif // REDUCEDMODELOPTIMIZER_STRUCTS_HPP diff --git a/simulation_structs.hpp b/simulation_structs.hpp index 344779d..835d8f7 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -23,6 +23,13 @@ void read_binary(const std::string &filename, Matrix &matrix) { in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar)); in.close(); } +const static IOFormat CSVFormat(StreamPrecision, DontAlignCols, ", ", "\n"); +template +void writeToCSV(const std::string &filename, Matrix &matrix) +{ + ofstream file(filename.c_str()); + file << matrix.format(CSVFormat); +} } // namespace Eigen #include "simulationmesh.hpp" @@ -195,7 +202,6 @@ public: bool load(const std::string &jsonFilename, const bool &shouldLoadMesh = true) { - std::cout << jsonFilename << std::endl; label = "empty_job"; constrainedVertices.clear(); nodalExternalForces.clear(); @@ -406,11 +412,11 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m } }); - if (!nodeColors.empty()) { - polyscope::getCurveNetwork(meshLabel) - ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors) - ->setEnabled(shouldEnable); - } + // if (!nodeColors.empty()) { + // polyscope::getCurveNetwork(meshLabel) + // ->addNodeColorQuantity("Boundary conditions_" + label, nodeColors) + // ->setEnabled(shouldEnable); + // } // per node external forces std::vector> externalForces(pMesh->VN()); @@ -421,9 +427,15 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m } if (!externalForces.empty()) { - polyscope::getCurveNetwork(meshLabel) - ->addNodeVectorQuantity("External force_" + label, externalForces) - ->setEnabled(shouldEnable); + polyscope::CurveNetworkNodeVectorQuantity *externalForcesVectors + = polyscope::getCurveNetwork(meshLabel)->addNodeVectorQuantity("External force_" + + label, + externalForces); + + const std::array color_loads{1.0, 0, 0}; + externalForcesVectors->setVectorColor( + glm::vec3(color_loads[0], color_loads[1], color_loads[2])); + externalForcesVectors->setEnabled(shouldEnable); } // per node external moments bool hasExternalMoments = false; @@ -454,6 +466,19 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m if (!constrainedVertices.empty()) { polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label); } + + // per node external moments + bool hasExternalMoments = false; + for (const auto &forcePair : nodalExternalForces) { + const Vector6d &load = forcePair.second; + if (load.getRotation().norm() != 0) { + hasExternalMoments = true; + break; + } + } + if (hasExternalMoments) { + polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment_" + label); + } } #endif // POLYSCOPE_DEFINED }; @@ -463,7 +488,7 @@ struct SimulationResults * vector contains the same info using euler angles */ bool converged{false}; - std::shared_ptr job; + std::shared_ptr pJob; SimulationHistory history; std::vector debug_drmDisplacements; std::vector> debug_q_f1; //per vertex @@ -475,7 +500,7 @@ struct SimulationResults double executionTime{0}; std::string labelPrefix{"deformed"}; inline static char deliminator{' '}; - SimulationResults() { job = std::make_shared(); } + SimulationResults() { pJob = std::make_shared(); } std::vector getTranslationalDisplacements() const { @@ -491,13 +516,13 @@ struct SimulationResults void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; } std::string getLabel() const { - return labelPrefix + deliminator + job->pMesh->getLabel() + deliminator + job->getLabel(); + return labelPrefix + deliminator + pJob->pMesh->getLabel() + deliminator + pJob->getLabel(); } bool saveDeformedModel(const std::string &outputFolder = std::string()) { VCGEdgeMesh m; - vcg::tri::Append::MeshCopy(m, *job->pMesh); + vcg::tri::Append::MeshCopy(m, *pJob->pMesh); for (int vi = 0; vi < m.VN(); vi++) { m.vert[vi].P() = m.vert[vi].P() @@ -513,10 +538,11 @@ struct SimulationResults const std::filesystem::path outputFolderPath = outputFolder.empty() ? std::filesystem::current_path() : std::filesystem::path(outputFolder); + std::cout << "Saving results to:" << outputFolderPath << std::endl; std::filesystem::path simulationJobOutputFolderPath = std::filesystem::path(outputFolderPath).append("SimulationJob"); std::filesystem::create_directories(simulationJobOutputFolderPath); - job->save(simulationJobOutputFolderPath.string()); + pJob->save(simulationJobOutputFolderPath.string()); const std::string filename(getLabel() + "_displacements.eigenBin"); Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); @@ -543,7 +569,7 @@ struct SimulationResults void load(const std::filesystem::path &loadFromPath, const std::filesystem::path &loadJobFrom) { //load job - job->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string()); + pJob->load(std::filesystem::path(loadJobFrom).append("SimulationJob.json").string()); //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") { @@ -554,8 +580,8 @@ struct SimulationResults } } - rotationalDisplacementQuaternion.resize(job->pMesh->VN()); - for (int vi = 0; vi < job->pMesh->VN(); vi++) { + rotationalDisplacementQuaternion.resize(pJob->pMesh->VN()); + for (int vi = 0; vi < pJob->pMesh->VN(); vi++) { rotationalDisplacementQuaternion[vi] = Eigen::AngleAxisd(displacements[vi][3], Eigen::Vector3d::UnitX()) * Eigen::AngleAxisd(displacements[vi][4], @@ -573,15 +599,18 @@ struct SimulationResults std::cerr << "Nothing to remove." << std::endl; return; } - job->unregister(getLabel()); + pJob->unregister(getLabel()); polyscope::removeCurveNetwork(getLabel()); } - void registerForDrawing(const std::optional> &desiredColor = std::nullopt, - const bool &shouldEnable = true, - const bool &shouldShowFrames = false) const + polyscope::CurveNetwork *registerForDrawing( + const std::optional> &desiredColor = std::nullopt, + const bool &shouldEnable = true, + const double &desiredRadius = 0.001, + // const double &desiredRadius = 0.0001, + const bool &shouldShowFrames = false) const { PolyscopeInterface::init(); - const std::shared_ptr &mesh = job->pMesh; + const std::shared_ptr &mesh = pJob->pMesh; polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh; if (!polyscope::hasStructure(getLabel())) { polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(), @@ -589,7 +618,7 @@ struct SimulationResults mesh->getEigenEdges()); } polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable); - polyscopeHandle_deformedEdmeMesh->setRadius(0.002, true); + polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true); if (desiredColor.has_value()) { const glm::vec3 desiredColor_glm(desiredColor.value()[0], desiredColor.value()[1], @@ -600,42 +629,94 @@ struct SimulationResults Eigen::MatrixX3d framesX(mesh->VN(), 3); Eigen::MatrixX3d framesY(mesh->VN(), 3); Eigen::MatrixX3d framesZ(mesh->VN(), 3); + + Eigen::MatrixX3d framesX_initial(mesh->VN(), 3); + Eigen::MatrixX3d framesY_initial(mesh->VN(), 3); + Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3); + + // std::unordered_set interfaceNodes{1, 3, 5, 7, 9, 11}; + std::unordered_set interfaceNodes{3, 7, 11, 15, 19, 23}; for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { const Vector6d &nodalDisplacement = displacements[vi]; nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0], nodalDisplacement[1], nodalDisplacement[2]); - auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1); - auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0); - auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0); - framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], deformedFrameX[1], deformedFrameX[2]); - framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], deformedFrameY[1], deformedFrameY[2]); - framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], deformedNormal[1], deformedNormal[2]); + if (interfaceNodes.contains(vi)) { + auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1); + auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0); + auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0); + framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], + deformedFrameX[1], + deformedFrameX[2]); + framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], + deformedFrameY[1], + deformedFrameY[2]); + framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], + deformedNormal[1], + deformedNormal[2]); + framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0); + framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0); + framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1); + } else { + framesX.row(vi) = Eigen::Vector3d(0, 0, 0); + framesY.row(vi) = Eigen::Vector3d(0, 0, 0); + framesZ.row(vi) = Eigen::Vector3d(0, 0, 0); + framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + } } polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices() + nodalDisplacements); + const double frameRadius_default = 0.035; + const double frameLength_default = 0.035; + const bool shouldEnable_default = true; //if(showFramesOn.contains(vi)){ auto polyscopeHandle_frameX = polyscopeHandle_deformedEdmeMesh ->addNodeVectorQuantity("FrameX", framesX); - polyscopeHandle_frameX->setVectorLengthScale(0.01); - polyscopeHandle_frameX->setVectorRadius(0.01); + polyscopeHandle_frameX->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameX->setVectorRadius(frameRadius_default); 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->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameY->setVectorRadius(frameRadius_default); 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->setVectorLengthScale(frameLength_default); + polyscopeHandle_frameZ->setVectorRadius(frameRadius_default); polyscopeHandle_frameZ->setVectorColor( /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); - //} - job->registerForDrawing(getLabel()); + + // auto polyscopeHandle_initialMesh = polyscope::getCurveNetwork(mesh->getLabel()); + // if (!polyscopeHandle_initialMesh) { + // polyscopeHandle_initialMesh = mesh->registerForDrawing(); + // } + + // auto polyscopeHandle_frameX_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameX", framesX_initial); + // polyscopeHandle_frameX_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameX_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameX_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(1, 0, 0)); + // auto polyscopeHandle_frameY_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameY", framesY_initial); + // polyscopeHandle_frameY_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameY_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameY_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 1, 0)); + // auto polyscopeHandle_frameZ_initial = polyscopeHandle_initialMesh + // ->addNodeVectorQuantity("FrameZ", framesZ_initial); + // polyscopeHandle_frameZ_initial->setVectorLengthScale(frameLength_default); + // polyscopeHandle_frameZ_initial->setVectorRadius(frameRadius_default); + // polyscopeHandle_frameZ_initial->setVectorColor( + // /*polyscope::getNextUniqueColor()*/ glm::vec3(0, 0, 1)); + // //} + pJob->registerForDrawing(getLabel()); static bool wasExecuted = false; if (!wasExecuted) { std::function callback = [&]() { @@ -651,6 +732,7 @@ struct SimulationResults PolyscopeInterface::addUserCallback(callback); wasExecuted = true; } + return polyscopeHandle_deformedEdmeMesh; } #endif }; diff --git a/simulationmesh.cpp b/simulationmesh.cpp index 1bd4076..65b0f2e 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -1,5 +1,5 @@ #include "simulationmesh.hpp" -#include +//#include SimulationMesh::SimulationMesh() { elements = vcg::tri::Allocator::GetPerEdgeAttribute( @@ -292,23 +292,23 @@ bool SimulationMesh::load(const string &plyFilename) { this->Clear(); // assert(plyFileHasAllRequiredFields(plyFilename)); // Load the ply file - VCGEdgeMesh::PerEdgeAttributeHandle handleBeamDimensions = - vcg::tri::Allocator::AddPerEdgeAttribute< - CrossSectionType>(*this, plyPropertyBeamDimensionsID); - VCGEdgeMesh::PerEdgeAttributeHandle handleBeamMaterial = - vcg::tri::Allocator::AddPerEdgeAttribute( - *this, plyPropertyBeamMaterialID); - nanoply::NanoPlyWrapper::CustomAttributeDescriptor - customAttrib; - customAttrib.GetMeshAttrib(plyFilename); - customAttrib.AddEdgeAttribDescriptor( - plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); - /*FIXME: Since I allow CrossSectionType to take two types I should export the - * type as well such that that when loaded the correct type of cross section - * is used. - */ - customAttrib.AddEdgeAttribDescriptor( - plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); + // VCGEdgeMesh::PerEdgeAttributeHandle handleBeamDimensions = + // vcg::tri::Allocator::AddPerEdgeAttribute< + // CrossSectionType>(*this, plyPropertyBeamDimensionsID); + // VCGEdgeMesh::PerEdgeAttributeHandle handleBeamMaterial = + // vcg::tri::Allocator::AddPerEdgeAttribute( + // *this, plyPropertyBeamMaterialID); + // nanoply::NanoPlyWrapper::CustomAttributeDescriptor + // customAttrib; + // customAttrib.GetMeshAttrib(plyFilename); + // customAttrib.AddEdgeAttribDescriptor( + // plyPropertyBeamDimensionsID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); + // /*FIXME: Since I allow CrossSectionType to take two types I should export the + // * type as well such that that when loaded the correct type of cross section + // * is used. + // */ + // customAttrib.AddEdgeAttribDescriptor( + // plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); // VCGEdgeMesh::PerEdgeAttributeHandle> // handleBeamProperties = @@ -322,21 +322,20 @@ bool SimulationMesh::load(const string &plyFilename) { // customAttrib.AddEdgeAttribDescriptor( // plyPropertyBeamRigidityConstantsID, nanoply::NNP_LIST_INT8_FLOAT32, // nullptr); - unsigned int mask = 0; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; - mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; - mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; - mask |= nanoply::NanoPlyWrapper::IO_MESHATTRIB; - if (nanoply::NanoPlyWrapper::LoadModel( - plyFilename.c_str(), *this, mask, customAttrib) != 0) { - return false; + // unsigned int mask = 0; + // mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; + // mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + // mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; + // mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; + // mask |= nanoply::NanoPlyWrapper::IO_MESHATTRIB; + if (!load(plyFilename)) { + return false; } - elements = vcg::tri::Allocator::GetPerEdgeAttribute( - *this, std::string("Elements")); - nodes = vcg::tri::Allocator::GetPerVertexAttribute( - *this, std::string("Nodes")); + // elements = vcg::tri::Allocator::GetPerEdgeAttribute( + // *this, std::string("Elements")); + // nodes = vcg::tri::Allocator::GetPerVertexAttribute( + // *this, std::string("Nodes")); vcg::tri::UpdateTopology::VertexEdge(*this); initializeNodes(); initializeElements(); @@ -349,11 +348,11 @@ bool SimulationMesh::load(const string &plyFilename) { // elements[ei].updateRigidity(); // } // } - for (size_t ei = 0; ei < EN(); ei++) { - elements[ei].setDimensions(handleBeamDimensions[ei]); - elements[ei].setMaterial(handleBeamMaterial[ei]); - elements[ei].updateRigidity(); - } + // for (size_t ei = 0; ei < EN(); ei++) { + // elements[ei].setDimensions(handleBeamDimensions[ei]); + // elements[ei].setMaterial(handleBeamMaterial[ei]); + // elements[ei].updateRigidity(); + // } bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; if (normalsAreAbsent) { @@ -375,28 +374,31 @@ bool SimulationMesh::save(const std::string &plyFilename) if (filename.empty()) { filename = std::filesystem::current_path().append(getLabel() + ".ply").string(); } - nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; - customAttrib.GetMeshAttrib(filename); + // nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; + // customAttrib.GetMeshAttrib(filename); - std::vector dimensions = getBeamDimensions(); - customAttrib.AddEdgeAttribDescriptor(plyPropertyBeamDimensionsID, - nanoply::NNP_LIST_INT8_FLOAT64, - dimensions.data()); - std::vector material = getBeamMaterial(); - customAttrib.AddEdgeAttribDescriptor(plyPropertyBeamMaterialID, - nanoply::NNP_LIST_INT8_FLOAT64, - material.data()); - unsigned int mask = 0; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; - mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; - mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; - mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; - if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), - *this, - mask, - customAttrib, - false) - != 1) { + // std::vector dimensions = getBeamDimensions(); + // customAttrib.AddEdgeAttribDescriptor(plyPropertyBeamDimensionsID, + // nanoply::NNP_LIST_INT8_FLOAT64, + // dimensions.data()); + // std::vector material = getBeamMaterial(); + // customAttrib.AddEdgeAttribDescriptor(plyPropertyBeamMaterialID, + // nanoply::NNP_LIST_INT8_FLOAT64, + // material.data()); + // unsigned int mask = 0; + // mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; + // mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; + // mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; + // mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + // if (nanoply::NanoPlyWrapper::SaveModel(filename.c_str(), + // *this, + // mask, + // customAttrib, + // false) + // != 1) { + // return false; + // } + if (!save(plyFilename)) { return false; } return true; @@ -449,31 +451,34 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, const VectorType cp = vector0 ^ vector1; if (cp.dot(normalVector) < 0) { - return -angle; + return -angle; } return angle; } void Element::computeMaterialProperties(const ElementMaterial &material) { - this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); + G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); +} + +void Element::computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A) +{ + A = dimensions.b * dimensions.h; } void Element::computeDimensionsProperties( const RectangularBeamDimensions &dimensions) { assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); - A = (dimensions.b * dimensions.h); - I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; - I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; - J = I2 + I3; + computeCrossSectionArea(dimensions, A); + computeMomentsOfInertia(dimensions, inertia); } void Element::computeDimensionsProperties( const CylindricalBeamDimensions &dimensions) { assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); - I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64; - I3 = I2; - J = I2 + I3; + inertia.I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64; + inertia.I3 = inertia.I2; + inertia.J = inertia.I2 + inertia.I3; } void Element::setDimensions(const CrossSectionType &dimensions) { @@ -488,9 +493,23 @@ void Element::setMaterial(const ElementMaterial &material) { updateRigidity(); } +double Element::getMass(const double &materialDensity) +{ + const double beamVolume = A * length; + return beamVolume * materialDensity; +} + +void Element::computeMomentsOfInertia(const RectangularBeamDimensions &dimensions, + Element::MomentsOfInertia &inertia) +{ + inertia.I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; + inertia.I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; + inertia.J = inertia.I2 + inertia.I3; +} + void Element::updateRigidity() { rigidity.axial = material.youngsModulus * A / initialLength; - rigidity.torsional = G * J / initialLength; - rigidity.firstBending = 2 * material.youngsModulus * I2 / initialLength; - rigidity.secondBending = 2 * material.youngsModulus * I3 / initialLength; + rigidity.torsional = G * inertia.J / initialLength; + rigidity.firstBending = 2 * material.youngsModulus * inertia.I2 / initialLength; + rigidity.secondBending = 2 * material.youngsModulus * inertia.I3 / initialLength; } diff --git a/simulationmesh.hpp b/simulationmesh.hpp index a7b85c9..35cc499 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -64,15 +64,23 @@ struct Element { ElementMaterial material; double G{0}; double A{0}; // cross sectional area - double I2{0}; // second moment of inertia - double I3{0}; // third moment of inertia - double J{0}; // torsional constant (polar moment of inertia) void computeMaterialProperties(const ElementMaterial &material); void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); void setDimensions(const CrossSectionType &dimensions); void setMaterial(const ElementMaterial &material); + double getMass(const double &matrialDensity); + + struct MomentsOfInertia + { + double I2{0}; // second moment of inertia + double I3{0}; // third moment of inertia + double J{0}; // torsional constant (polar moment of inertia) + } inertia; + + static void computeMomentsOfInertia(const RectangularBeamDimensions &dimensions, + MomentsOfInertia &inertia); struct LocalFrame { VectorType t1; @@ -127,6 +135,8 @@ struct Element { }; RotationalDisplacements rotationalDisplacements_j; RotationalDisplacements rotationalDisplacements_jplus1; + + static void computeCrossSectionArea(const RectangularBeamDimensions &dimensions, double &A); }; struct Node { diff --git a/simulationmodel.hpp b/simulationmodel.hpp new file mode 100644 index 0000000..6842260 --- /dev/null +++ b/simulationmodel.hpp @@ -0,0 +1,13 @@ +#ifndef SIMULATIONMODEL_HPP +#define SIMULATIONMODEL_HPP + +#include "simulation_structs.hpp" + +class SimulationModel +{ +public: + virtual SimulationResults executeSimulation(const std::shared_ptr &simulationJob) + = 0; +}; + +#endif // SIMULATIONMODEL_HPP diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index 4bad741..261f4e0 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -152,6 +152,9 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced } } + const std::unordered_set interfaceNodes = patternGeometryAllEdges.getInterfaceNodes( + numberOfNodesPerSlot); + // assert(validEdges.size() == allPossibleEdges.size() - // coincideEdges.size() - // duplicateEdges.size()); @@ -191,7 +194,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced perEdgeResultPath, patternGeometryAllEdges.getVertices(), intersectingEdges, - validEdges); + validEdges, + interfaceNodes); statistics.print(setupString, perEdgeResultPath); statistics.reset(); } @@ -211,7 +215,8 @@ void TopologyEnumerator::computeValidPatterns(const std::vector &reduced perEdgeResultPath, patternGeometryAllEdges.getVertices(), intersectingEdges, - validEdges); + validEdges, + interfaceNodes); statistics.print(setupString, perEdgeResultPath); } } @@ -395,7 +400,8 @@ void TopologyEnumerator::computeValidPatterns( const std::filesystem::path &resultsPath, const std::vector &allVertices, const std::unordered_map> &intersectingEdges, - const std::vector &validEdges) + const std::vector &validEdges, + const std::unordered_set &interfaceNodes) { assert(numberOfNodesPerSlot.size() == 7); // Iterate over all patterns which have numberOfDesiredEdges edges from @@ -435,7 +441,7 @@ void TopologyEnumerator::computeValidPatterns( // std::string previousPatternBinaryRepresentation(validEdges.size(),'0'); size_t patternIndex = 0; bool validPatternsExist = false; - const bool exportTilledPattern = debugIsOn; + const bool exportTilledPattern = true; const bool saveCompressedFormat = false; do { patternIndex++; @@ -459,6 +465,23 @@ void TopologyEnumerator::computeValidPatterns( patternGeometry.add(allVertices, patternEdges); patternGeometry.setLabel(patternName); + // Check if pattern contains intersecting edges + const bool isInterfaceConnected = patternGeometry.isInterfaceConnected(interfaceNodes); + // Export the tiled ply file if it contains intersecting edges + if (!isInterfaceConnected) { + // create the tiled geometry of the pattern + statistics.numberOfPatternViolatingInterfaceEnforcement++; + if (debugIsOn) { + if (savePlyFiles) { + exportPattern(std::filesystem::path(resultsPath).append("InterfaceEnforcement"), + patternGeometry, + exportTilledPattern); + } + } else { + continue; // should be uncommented in order to improve performance + } + } + // Check if pattern contains intersecting edges const bool patternContainsIntersectingEdges = patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, intersectingEdges); @@ -477,6 +500,7 @@ void TopologyEnumerator::computeValidPatterns( } } + // const bool shouldBreak = numberOfDesiredEdges == 4 && patternIndex == 53; const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 15); if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { @@ -512,7 +536,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges( numberOfNodesPerSlot); // marks the nodes with valence>=1 // Create the tiled geometry of the pattern - const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenTiled(); + const bool hasFloatingComponents = !patternGeometry.isFullyConnectedWhenFanned(); PatternGeometry fanPatternGeometry = PatternGeometry::createFan(patternGeometry); const int interfaceNodeVi = 3; @@ -559,10 +583,19 @@ void TopologyEnumerator::computeValidPatterns( std::vector articulationPoints; boost::articulation_points(fanPatternGraph, std::back_inserter(articulationPoints)); const bool hasArticulationPoints = !articulationPoints.empty(); - // if (hasArticulationPoints /*&& !patternContainsIntersectingEdges - // && !tiledPatternHasDanglingEdges && !hasFloatingComponents - // && !tiledPatternHasNodeWithValenceGreaterThanDesired - // && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) { + // if (!hasArticulationPoints && tiledPatternHasDanglingEdges) { + // PatternGeometry tilledPatternGeometry = PatternGeometry::createTile(patternGeometry); + // tilledPatternGeometry.updateEigenEdgeAndVertices(); + // tilledPatternGeometry.registerForDrawing(); + // // ->addNodeColorQuantity("de_noAp_tilled", fanVertexColors) + // // ->setEnabled(true); + // polyscope::show(); + // tilledPatternGeometry.unregister(); + // } + // if (hasArticulationPoints && !tiledPatternHasDanglingEdges/*&& !patternContainsIntersectingEdges + // && !hasFloatingComponents + // && !tiledPatternHasNodeWithValenceGreaterThanDesired + // && !tiledPatternHasEdgesWithAngleSmallerThanThreshold*/) { // for (PatternGeometry::VertexType &v : patternGeometry.vert) { // v.C() = vcg::Color4b::Yellow; // } @@ -582,12 +615,12 @@ void TopologyEnumerator::computeValidPatterns( // const auto vi = tilledPatternGeometry.getIndex(v); // fanVertexColors[vi] = vColor; // } - // // tilledPatternGeometry.updateEigenEdgeAndVertices(); - // // tilledPatternGeometry.registerForDrawing() - // // ->addNodeColorQuantity("ap_tilled", fanVertexColors) - // // ->setEnabled(true); - // // polyscope::show(); - // // tilledPatternGeometry.unregister(); + // tilledPatternGeometry.updateEigenEdgeAndVertices(); + // tilledPatternGeometry.registerForDrawing() + // ->addNodeColorQuantity("ap_tilled", fanVertexColors) + // ->setEnabled(true); + // polyscope::show(); + // tilledPatternGeometry.unregister(); // } // duplicated here // Check dangling edges with vcg method @@ -687,10 +720,20 @@ void TopologyEnumerator::computeValidPatterns( } const bool isValidPattern = !patternContainsIntersectingEdges - && !tiledPatternHasDanglingEdges && !hasFloatingComponents - && !hasArticulationPoints + && isInterfaceConnected + /*&& !tiledPatternHasDanglingEdges*/ + && !hasFloatingComponents && !hasArticulationPoints && !tiledPatternHasNodeWithValenceGreaterThanDesired && !tiledPatternHasEdgesWithAngleSmallerThanThreshold; + // if (!hasArticulationPoints && !patternContainsIntersectingEdges + // && !tiledPatternHasDanglingEdges && !hasFloatingComponents + // && !tiledPatternHasNodeWithValenceGreaterThanDesired + // && tiledPatternHasEdgesWithAngleSmallerThanThreshold && numberOfDesiredEdges > 4) { + // std::cout << "Pattern found:" << patternName << std::endl; + // patternGeometry.registerForDrawing(); + // polyscope::show(); + // patternGeometry.unregister(); + // } if (isValidPattern) { // if(patternName=='2055'){ // PatternGeometry tiledPatternGeometry = PatternGeometry::createTile( diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp index 0de1444..9a2157e 100755 --- a/topologyenumerator.hpp +++ b/topologyenumerator.hpp @@ -64,25 +64,27 @@ class TopologyEnumerator { size_t numberOfPatternsWithArticulationPoints{0}; size_t numberOfPatternsViolatingAngleThreshold{0}; size_t numberOfPatternsViolatingValenceThreshold{0}; + size_t numberOfPatternViolatingInterfaceEnforcement{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["numPossibleEdges"] = numberOfPossibleEdges; + // json["numCoincideEdges"] = numberOfCoincideEdges; + // json["numDuplicateEdges"] = numberOfDuplicateEdges; + // json["numValidEdges"] = numberOfValidEdges; + // json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs; json["numPatterns"] = numberOfPatterns; // json["numIntersectingEdgesOverAllPatterns"] = // numberOfIntersectingEdgesOverAllPatterns; json["numPatternsWithIntersectingEdges"] = numberOfPatternsWithIntersectingEdges; + json["numViolatingInterfaceEnforcement"] = numberOfPatternViolatingInterfaceEnforcement; json["numPatternsWithNotASingleCC"] = numberOfPatternsWithMoreThanASingleCC; - json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; json["numPatternsWithArticulationPoints"] = numberOfPatternsWithArticulationPoints; - json["numValidPatterns"] = numberOfValidPatterns; + json["numPatternsWithDangling"] = numberOfPatternsWithADanglingEdgeOrNode; json["violatingAngle"] = numberOfPatternsViolatingAngleThreshold; json["violatingValence"] = numberOfPatternsViolatingValenceThreshold; + json["numValidPatterns"] = numberOfValidPatterns; return json; } @@ -107,9 +109,10 @@ class TopologyEnumerator { // std::endl; std::cout << numberOfPatternsWithIntersectingEdges << " patterns found with intersecting edges" << std::endl; + std::cout << numberOfPatternViolatingInterfaceEnforcement + << " patterns found which violate interface enforcement" << std::endl; std::cout << numberOfPatternsWithMoreThanASingleCC - << " patterns found with more than one connected components" - << std::endl; + << " patterns found with more than one connected components" << std::endl; std::cout << numberOfPatternsWithADanglingEdgeOrNode << " patterns found with a dangling node or edge" << std::endl; std::cout << numberOfPatternsWithArticulationPoints @@ -141,6 +144,8 @@ class TopologyEnumerator { numberOfIntersectingEdgePairs = 0; numberOfPatterns = 0; numberOfPatternsWithIntersectingEdges = 0; + numberOfPatternViolatingInterfaceEnforcement = 0; + ; numberOfPatternsWithMoreThanASingleCC = 0; numberOfPatternsWithADanglingEdgeOrNode = 0; numberOfPatternsWithArticulationPoints = 0; @@ -195,7 +200,8 @@ private: const std::filesystem::path &resultsPath, const std::vector &vertices, const std::unordered_map> &intersectingEdges, - const std::vector &validEdges); + const std::vector &validEdges, + const std::unordered_set &interfaceNodes); void exportPattern(const std::filesystem::path &saveToPath, PatternGeometry &patternGeometry, const bool saveTilledPattern) const; diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index b8f4b22..74d67be 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -178,10 +178,11 @@ void PatternGeometry::add(const std::vector &edges) { } void PatternGeometry::add(const std::vector &vertices, - const std::vector &edges) { - add(vertices); - add(edges); - updateEigenEdgeAndVertices(); + const std::vector &edges) +{ + add(vertices); + add(edges); + updateEigenEdgeAndVertices(); } void PatternGeometry::add(const std::vector &numberOfNodesPerSlot, @@ -375,9 +376,9 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector &nu // 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)); + R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(-360.0 / fanSize)); } else { - R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(60.0)); + R = vcg::RotationMatrix(vcg::Point3d(0, 0, 1), vcg::math::ToRad(360.0 / fanSize)); } // const CoordType &correspondingVertexPosition_rotated = vert[vi].cP(); std::transform( @@ -452,6 +453,7 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector &nu } return false; } + bool PatternGeometry::hasValenceGreaterThan(const std::vector &numberOfNodesPerSlot, const size_t &valenceThreshold) { @@ -561,111 +563,111 @@ void PatternGeometry::constructCorrespondingNodeMap(const std::vector &n } } -bool PatternGeometry::isFullyConnectedWhenTiled() { - assert(vn != 0 /* && !correspondingNode.empty()*/); - // TrianglePatternGeometry copyOfPattern(*this); +bool PatternGeometry::isFullyConnectedWhenFanned() +{ + assert(vn != 0 /* && !correspondingNode.empty()*/); + // TrianglePatternGeometry copyOfPattern(*this); - // If bottom interface nodes have a valence of zero there definetely more than - // a single cc - bool bottomInterfaceIsConnected = false; - for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { - 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; - } + // // If bottom interface nodes have a valence of zero there definetely more than + // // a single cc + // bool bottomInterfaceIsConnected = false; + // for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + // 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; + // } - PatternGeometry fanedPattern = createFan(*this); - vcg::tri::Clean::MergeCloseVertex(fanedPattern, 0.000000005); - vcg::tri::Allocator::CompactEveryVector(fanedPattern); - vcg::tri::UpdateTopology::VertexEdge(fanedPattern); - vcg::tri::UpdateTopology::EdgeEdge(fanedPattern); - std::vector> eCC; - vcg::tri::Clean::edgeMeshConnectedComponents(fanedPattern, - eCC); + PatternGeometry fanedPattern = createFan(*this); + vcg::tri::Clean::MergeCloseVertex(fanedPattern, 0.000000005); + vcg::tri::Allocator::CompactEveryVector(fanedPattern); + vcg::tri::UpdateTopology::VertexEdge(fanedPattern); + vcg::tri::UpdateTopology::EdgeEdge(fanedPattern); + std::vector> eCC; + vcg::tri::Clean::edgeMeshConnectedComponents(fanedPattern, eCC); - const bool sideInterfaceIsConnected = 1 == eCC.size(); - if (!sideInterfaceIsConnected) { - return false; - } + const bool sideInterfaceIsConnected = 1 == eCC.size(); + if (!sideInterfaceIsConnected) { + return false; + } - // for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { - // if (nodeTiledValence[nodeIndex] == 0) - // continue; - // const size_t slotIndex = nodeSlot[nodeIndex]; - // // connect nodes belonging to first or third slot(nodes on the first - // // triangle edge) - // // if (nodeTiledValence[nodeIndex] == 0 && slotIndex == 3) { - // // continue; - // // } - // if (slotIndex == 1 || slotIndex == 3) { - // assert(correspondingNode.count(nodeIndex) != 0); - // const size_t correspondingNodeIndex = - // correspondingNode[nodeIndex]; - // std::vector starEdges; - // vcg::edge::VEStarVE(&vert[nodeIndex], starEdges); - // bool containsEdge = false; - // for (TrianglePatternGeometry::EdgeType *e : starEdges) { - // assert(vcg::tri::Index(*this, e->V(0)) == nodeIndex || - // vcg::tri::Index(*this, e->V(1)) == nodeIndex); - // if (vcg::tri::Index(*this, e->V(0)) == nodeIndex) { - // if (vcg::tri::Index(*this, e->V(1)) == correspondingNodeIndex) { - // containsEdge = true; - // break; - // } - // } else if (vcg::tri::Index(*this, e->V(1)) == nodeIndex) { - // if (vcg::tri::Index(*this, e->V(0)) == correspondingNodeIndex) { - // containsEdge = true; - // break; - // } - // } - // } - // if (!containsEdge) { - // vcg::tri::Allocator::AddEdge( - // copyOfPattern, nodeIndex, correspondingNodeIndex); - // } - // } else if (slotIndex == 2 || slotIndex == 5) { - // assert(correspondingNode.count(nodeIndex) != 0); - // } else { - // assert(correspondingNode.count(nodeIndex) == 0); - // } - // } + // for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + // if (nodeTiledValence[nodeIndex] == 0) + // continue; + // const size_t slotIndex = nodeSlot[nodeIndex]; + // // connect nodes belonging to first or third slot(nodes on the first + // // triangle edge) + // // if (nodeTiledValence[nodeIndex] == 0 && slotIndex == 3) { + // // continue; + // // } + // if (slotIndex == 1 || slotIndex == 3) { + // assert(correspondingNode.count(nodeIndex) != 0); + // const size_t correspondingNodeIndex = + // correspondingNode[nodeIndex]; + // std::vector starEdges; + // vcg::edge::VEStarVE(&vert[nodeIndex], starEdges); + // bool containsEdge = false; + // for (TrianglePatternGeometry::EdgeType *e : starEdges) { + // assert(vcg::tri::Index(*this, e->V(0)) == nodeIndex || + // vcg::tri::Index(*this, e->V(1)) == nodeIndex); + // if (vcg::tri::Index(*this, e->V(0)) == nodeIndex) { + // if (vcg::tri::Index(*this, e->V(1)) == correspondingNodeIndex) { + // containsEdge = true; + // break; + // } + // } else if (vcg::tri::Index(*this, e->V(1)) == nodeIndex) { + // if (vcg::tri::Index(*this, e->V(0)) == correspondingNodeIndex) { + // containsEdge = true; + // break; + // } + // } + // } + // if (!containsEdge) { + // vcg::tri::Allocator::AddEdge( + // copyOfPattern, nodeIndex, correspondingNodeIndex); + // } + // } else if (slotIndex == 2 || slotIndex == 5) { + // assert(correspondingNode.count(nodeIndex) != 0); + // } else { + // assert(correspondingNode.count(nodeIndex) == 0); + // } + // } - // std::vector> eCC; - // vcg::tri::Clean::edgeMeshConnectedComponents( - // copyOfPattern, eCC); - // size_t numberOfCC_edgeBased = eCC.size(); - // size_t numberOfCC_vertexBased = numberOfCC_edgeBased; - // if (numberOfCC_edgeBased == 1) { - // vcg::tri::UpdateTopology::VertexEdge( - // copyOfPattern); - // vcg::tri::UpdateTopology::EdgeEdge(copyOfPattern); - // vcg::tri::UpdateFlags::VertexSetV(copyOfPattern); - // vcg::tri::UpdateFlags::VertexClear(copyOfPattern); - // for (size_t ei = 0; ei < copyOfPattern.EN(); ei++) { - // copyOfPattern.edge[ei].V(0)->SetV(); - // copyOfPattern.edge[ei].V(1)->SetV(); - // } + // std::vector> eCC; + // vcg::tri::Clean::edgeMeshConnectedComponents( + // copyOfPattern, eCC); + // size_t numberOfCC_edgeBased = eCC.size(); + // size_t numberOfCC_vertexBased = numberOfCC_edgeBased; + // if (numberOfCC_edgeBased == 1) { + // vcg::tri::UpdateTopology::VertexEdge( + // copyOfPattern); + // vcg::tri::UpdateTopology::EdgeEdge(copyOfPattern); + // vcg::tri::UpdateFlags::VertexSetV(copyOfPattern); + // vcg::tri::UpdateFlags::VertexClear(copyOfPattern); + // for (size_t ei = 0; ei < copyOfPattern.EN(); ei++) { + // copyOfPattern.edge[ei].V(0)->SetV(); + // copyOfPattern.edge[ei].V(1)->SetV(); + // } - // for (size_t vi = 0; vi < copyOfPattern.VN(); vi++) { - // if (!copyOfPattern.vert[vi].IsV()) { - // numberOfCC_vertexBased++; - // } - // } - // return numberOfCC_vertexBased; - // } + // for (size_t vi = 0; vi < copyOfPattern.VN(); vi++) { + // if (!copyOfPattern.vert[vi].IsV()) { + // numberOfCC_vertexBased++; + // } + // } + // return numberOfCC_vertexBased; + // } - // return numberOfCC_edgeBased == 1; // TODO: not good - return true; + // return numberOfCC_edgeBased == 1; // TODO: not good + return true; } bool PatternGeometry::hasIntersectingEdges( @@ -1062,7 +1064,7 @@ void PatternGeometry::createFan(const std::vector &connectToNeighborsVi, co for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) { vcg::Matrix44d R; auto rotationAxis = vcg::Point3d(0, 0, 1); - R.SetRotateDeg(360 / fanSize, rotationAxis); + R.SetRotateDeg(360.0 / fanSize, rotationAxis); vcg::tri::UpdatePosition::Matrix(rotatedPattern, R); vcg::tri::Append::Mesh(*this, rotatedPattern); //Add edges for connecting the desired vertices @@ -1097,3 +1099,39 @@ double PatternGeometry::getBaseTriangleHeight() const { return baseTriangleHeight; } + +bool PatternGeometry::isInterfaceConnected(const std::unordered_set &interfaceNodes) +{ + std::unordered_set interfaceNodesConnected; + + for (int ei = 0; ei < EN(); ei++) { + const int evi0 = getIndex(edge[ei].cV(0)); + const int evi1 = getIndex(edge[ei].cV(1)); + if (interfaceNodes.contains(evi0) && !interfaceNodes.contains(evi1)) { + interfaceNodesConnected.insert(evi0); + } else if (!interfaceNodes.contains(evi0) && interfaceNodes.contains(evi1)) { + interfaceNodesConnected.insert(evi1); + } + } + + if (interfaceNodesConnected.size() != interfaceNodes.size()) { + return false; + } + return true; +} + +std::unordered_set PatternGeometry::getInterfaceNodes( + const std::vector &numberOfNodesPerSlot) +{ + if (nodeToSlotMap.empty()) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlotMap); + } + std::unordered_set interfaceNodes; + for (std::pair viSlotPair : nodeToSlotMap) { + if (viSlotPair.second == 4) { + interfaceNodes.insert(viSlotPair.first); + } + } + + return interfaceNodes; +} diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 8364830..ccce648 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -49,8 +49,7 @@ private: 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, - const std::vector &edges); + void add(const std::vector &vertices, const std::vector &edges); void add(const std::vector &numberOfNodesPerSlot, const std::vector &edges); static std::vector @@ -71,7 +70,7 @@ private: assert(n >= m); return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1)); } - bool isFullyConnectedWhenTiled(); + bool isFullyConnectedWhenFanned(); bool hasIntersectingEdges( const std::string &patternBinaryRepresentation, @@ -123,6 +122,8 @@ private: const std::vector &perSurfaceFacePatternIndices, std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexTiledPatternEdgeIndex); + std::unordered_set getInterfaceNodes(const std::vector &numberOfNodesPerSlot); + bool isInterfaceConnected(const std::unordered_set &interfaceNodes); }; #endif // FLATPATTERNGEOMETRY_HPP diff --git a/utilities.hpp b/utilities.hpp index aac3f0a..f388c3f 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -257,6 +257,7 @@ inline void init() polyscope::view::upDir = polyscope::view::UpDir::ZUp; polyscope::state::userCallback = &mainCallback; + // polyscope::options::autocenterStructures = true; } using PolyscopeLabel = std::string; inline std::pair getSelection() @@ -346,4 +347,17 @@ inline size_t computeHashOrdered(const std::vector &v) return std::hash{}(elementsString); } +//inline std::filesystem::path getFilepathWithExtension(const std::filesystem::path &folderPath, +// const std::string &extension) +//{ +// for (const std::filesystem::directory_entry &dirEntry : +// std::filesystem::directory_iterator(folderPath)) { +// if (dirEntry.is_regular_file() && std::filesystem::path(dirEntry).extension() == extension) { +// return std::filesystem::path(dirEntry); +// } +// } + +// return ""; +//} + #endif // UTILITIES_H diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp index c1c384d..486c619 100755 --- a/vcgtrimesh.cpp +++ b/vcgtrimesh.cpp @@ -3,17 +3,18 @@ #include "wrap/io_trimesh/import_obj.h" #include "wrap/io_trimesh/import_off.h" #include +#include #include -#include +//#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; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; - mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; - mask |= nanoply::NanoPlyWrapper::IO_FACEINDEX; + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_VERTNORMAL; + mask |= vcg::tri::io::Mask::IOM_VERTCOLOR; + mask |= vcg::tri::io::Mask::IOM_EDGEINDEX; + mask |= vcg::tri::io::Mask::IOM_FACEINDEX; // if (nanoply::NanoPlyWrapper::LoadModel( // std::filesystem::absolute(filename).string().c_str(), *this, mask) @@ -79,12 +80,11 @@ bool VCGTriMesh::save(const std::string plyFilename) assert(std::filesystem::path(plyFilename).extension() == ".ply"); // Load the ply file unsigned int mask = 0; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; - mask |= nanoply::NanoPlyWrapper::IO_VERTCOLOR; - mask |= nanoply::NanoPlyWrapper::IO_FACEINDEX; - mask |= nanoply::NanoPlyWrapper::IO_FACENORMAL; - if (nanoply::NanoPlyWrapper::SaveModel(plyFilename.c_str(), *this, mask, false) - != 0) { + mask |= vcg::tri::io::Mask::IOM_VERTCOORD; + mask |= vcg::tri::io::Mask::IOM_VERTCOLOR; + mask |= vcg::tri::io::Mask::IOM_FACEINDEX; + mask |= vcg::tri::io::Mask::IOM_FACENORMAL; + if (vcg::tri::io::Exporter::Save(*this, plyFilename.c_str(), mask) != 0) { return false; } return true; @@ -117,7 +117,7 @@ polyscope::SurfaceMesh *VCGTriMesh::registerForDrawing( const glm::vec3 desiredColor_glm(desiredColor.value()[0], desiredColor.value()[1], desiredColor.value()[2]); - polyscopeHandle_mesh->setSurfaceColor(desiredColor_glm); + polyscopeHandle_mesh->setSurfaceColor(glm::normalize(desiredColor_glm)); } return polyscopeHandle_mesh; } @@ -146,3 +146,14 @@ VCGTriMesh::VCGTriMesh(const std::string &filename) vcg::tri::UpdateTopology::VertexFace(*this); vcg::tri::UpdateNormal::PerVertexNormalized(*this); } + +void VCGTriMesh::moveToCenter() +{ + CoordType centerOfMass(0, 0, 0); + + for (int vi = 0; vi < VN(); vi++) { + centerOfMass += vert[vi].cP(); + } + centerOfMass /= VN(); + vcg::tri::UpdatePosition::Translate(*this, -centerOfMass); +} diff --git a/vcgtrimesh.hpp b/vcgtrimesh.hpp index 23ac9f1..6f79d02 100755 --- a/vcgtrimesh.hpp +++ b/vcgtrimesh.hpp @@ -54,6 +54,7 @@ public: #endif Eigen::MatrixX2i getEdges() const; void updateEigenEdgeAndVertices(); + void moveToCenter(); }; #endif // VCGTRIMESH_HPP