From 9668d030eab6d4e28d06f557dd8b800c423c635c Mon Sep 17 00:00:00 2001 From: Iason Date: Fri, 27 Nov 2020 12:47:21 +0200 Subject: [PATCH] Initial commit --- beam.hpp | 57 ++ beamformfinder.cpp | 1796 ++++++++++++++++++++++++++++++++++ beamformfinder.hpp | 224 +++++ edgemesh.cpp | 382 ++++++++ edgemesh.hpp | 105 ++ elementalmesh.cpp | 209 ++++ elementalmesh.hpp | 164 ++++ flatpattern.cpp | 397 ++++++++ flatpattern.hpp | 33 + patternIO.cpp | 116 +++ patternIO.hpp | 40 + polymesh.hpp | 64 ++ simulationhistoryplotter.hpp | 145 +++ simulationresult.hpp | 199 ++++ topologyenumerator.cpp | 609 ++++++++++++ topologyenumerator.hpp | 180 ++++ trianglepatterngeometry.cpp | 560 +++++++++++ trianglepatterngeometry.hpp | 68 ++ trianglepattterntopology.cpp | 298 ++++++ trianglepattterntopology.hpp | 49 + utilities.hpp | 272 +++++ vcgtrimesh.cpp | 78 ++ vcgtrimesh.hpp | 40 + 23 files changed, 6085 insertions(+) create mode 100644 beam.hpp create mode 100644 beamformfinder.cpp create mode 100644 beamformfinder.hpp create mode 100644 edgemesh.cpp create mode 100644 edgemesh.hpp create mode 100644 elementalmesh.cpp create mode 100644 elementalmesh.hpp create mode 100644 flatpattern.cpp create mode 100644 flatpattern.hpp create mode 100644 patternIO.cpp create mode 100644 patternIO.hpp create mode 100644 polymesh.hpp create mode 100644 simulationhistoryplotter.hpp create mode 100644 simulationresult.hpp create mode 100644 topologyenumerator.cpp create mode 100644 topologyenumerator.hpp create mode 100644 trianglepatterngeometry.cpp create mode 100644 trianglepatterngeometry.hpp create mode 100644 trianglepattterntopology.cpp create mode 100644 trianglepattterntopology.hpp create mode 100644 utilities.hpp create mode 100644 vcgtrimesh.cpp create mode 100644 vcgtrimesh.hpp diff --git a/beam.hpp b/beam.hpp new file mode 100644 index 0000000..3657070 --- /dev/null +++ b/beam.hpp @@ -0,0 +1,57 @@ +#ifndef BEAM_HPP +#define BEAM_HPP +#include +#include +#include +#include + +struct RectangularBeamDimensions { + float b; + float h; + RectangularBeamDimensions(const float &width, const float &height) + : b(width), h(height) { + assert(width > 0 && height > 0); + } + RectangularBeamDimensions() : b(1), h(1) {} +}; + +struct CylindricalElementDimensions { + float od; // Cylinder outside diameter + float + id; // Cylinder inside diameter + // https://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html + CylindricalElementDimensions(const float &outsideDiameter, + const float &insideDiameter) + : od(outsideDiameter), id(insideDiameter) { + assert(outsideDiameter > 0 && insideDiameter > 0 && + outsideDiameter > insideDiameter); + } + CylindricalElementDimensions() : od(0.03), id(0.026) {} +}; + +struct ElementMaterial { + float poissonsRatio; + float youngsModulusGPascal; + ElementMaterial(const float &poissonsRatio, const float &youngsModulusGPascal) + : poissonsRatio(poissonsRatio), + youngsModulusGPascal(youngsModulusGPascal) { + assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); + } + ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(7.5) {} +}; + +struct BeamProperties { + float crossArea; + float I2; + float I3; + float polarInertia; + BeamProperties(const RectangularBeamDimensions &dimensions, + const ElementMaterial &material) { + crossArea = (dimensions.b * dimensions.h); + I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; + I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; + polarInertia = (I2 + I3); + } +}; + +#endif // BEAM_HPP diff --git a/beamformfinder.cpp b/beamformfinder.cpp new file mode 100644 index 0000000..157837a --- /dev/null +++ b/beamformfinder.cpp @@ -0,0 +1,1796 @@ +#include "beamformfinder.hpp" +#include "polyscope/curve_network.h" +#include "polyscope/histogram.h" +#include "polyscope/polyscope.h" +#include "simulationhistoryplotter.hpp" +#include +#include +#include +#include + +void FormFinder::reset() { + Dt = Dtini; + t = 0; + currentSimulationStep = 0; + history.clear(); + // theta3Derivatives = + // Eigen::Tensor(DoF::NumDoFType, mesh->VN(), mesh->EN(), + // mesh->VN()); + // theta3Derivatives.setZero(); +} + +VectorType FormFinder::computeDisplacementDifferenceDerivative( + const EdgeType &e, const DifferentiateWithRespectTo &dui) const { + VectorType displacementDiffDeriv(0, 0, 0); + const DoFType &dofi = dui.dofi; + const bool differentiateWithRespectToANonEdgeNode = + e.cV(0) != &dui.v && e.cV(1) != &dui.v; + if (differentiateWithRespectToANonEdgeNode || dofi > 2) { + return displacementDiffDeriv; + } + + if (e.cV(0) == &dui.v) { + displacementDiffDeriv[dofi] = -1; + } else if (e.cV(1) == &dui.v) { + displacementDiffDeriv[dofi] = 1; + } + + return displacementDiffDeriv; +} + +VectorType FormFinder::computeDerivativeOfNormal( + const VertexType &v, const DifferentiateWithRespectTo &dui) const { + VectorType normalDerivative(0, 0, 0); + if (&dui.v != &v || + (dui.dofi == 0 || dui.dofi == 1 || dui.dofi == 2 || dui.dofi == 5)) { + return normalDerivative; + } + const VectorType &n = v.cN(); + const double &nx = n[0], ny = n[1]; + const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); + + if (dui.dofi == 3) { + if (nxnyMagnitude >= 1) { + const double normalDerivativeX = + 1 / sqrt(nxnyMagnitude) - + std::pow(nx, 2) / std::pow(nxnyMagnitude, 1.5); + const double normalDerivativeY = -nx * ny / std::pow(nxnyMagnitude, 1.5); + const double normalDerivativeZ = 0; + normalDerivative = + VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); + } else { + const double normalDerivativeX = 1; + const double normalDerivativeY = 0; + const double normalDerivativeZ = -nx / std::sqrt(1 - nxnyMagnitude); + normalDerivative = + VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); + } + } else if (dui.dofi == 4) { + if (nxnyMagnitude >= 1) { + const double normalDerivativeX = -nx * ny / std::pow(nxnyMagnitude, 1.5); + const double normalDerivativeY = + 1 / sqrt(nxnyMagnitude) - + std::pow(ny, 2) / std::pow(nxnyMagnitude, 1.5); + const double normalDerivativeZ = 0; + normalDerivative = + VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); + } else { + const double normalDerivativeX = 0; + const double normalDerivativeY = 1; + const double normalDerivativeZ = -ny / std::sqrt(1 - nxnyMagnitude); + normalDerivative = + VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); + } + } + + const bool normalDerivativeIsFinite = std::isfinite(normalDerivative[0]) && + std::isfinite(normalDerivative[1]) && + std::isfinite(normalDerivative[2]); + assert(normalDerivativeIsFinite); + + return normalDerivative; +} + +double FormFinder::computeDerivativeElementLength( + const EdgeType &e, const DifferentiateWithRespectTo &dui) const { + if (e.cV(0) != &dui.v && e.cV(1) != &dui.v) { + return 0; + } + + const VectorType &X_j = e.cP(0); + const VectorType &X_jplus1 = e.cP(1); + const VectorType positionVectorDiff = X_jplus1 - X_j; + const VectorType displacementDiffDeriv = + computeDisplacementDifferenceDerivative(e, dui); + const double edgeLength = mesh->elements[e].length; + const double L_kDeriv = + positionVectorDiff * displacementDiffDeriv / edgeLength; + return L_kDeriv; +} + +double FormFinder::computeDerivativeOfNorm(const VectorType &x, + const VectorType &derivativeOfX) { + return x.dot(derivativeOfX) / x.Norm(); +} + +VectorType FormFinder::computeDerivativeOfCrossProduct( + const VectorType &a, const VectorType &derivativeOfA, const VectorType &b, + const VectorType &derivativeOfB) { + const auto firstTerm = Cross(derivativeOfA, b); + const auto secondTerm = Cross(a, derivativeOfB); + return firstTerm + secondTerm; +} + +VectorType +FormFinder::computeDerivativeOfR(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const { + + const VertexType &v_j = *e.cV(0); + const VertexType &v_jplus1 = *e.cV(1); + const VectorType normal_j = v_j.cN(); + const VectorType normal_jplus1 = v_jplus1.cN(); + const VectorType derivativeOfNormal_j = + &v_j == &dui.v && dui.dofi > 2 + ? mesh->nodes[v_j].derivativeOfNormal[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeOfNormal_jplus1 = + &v_jplus1 == &dui.v && dui.dofi > 2 + ? mesh->nodes[v_jplus1].derivativeOfNormal[dui.dofi] + : VectorType(0, 0, 0); + + const VectorType derivativeOfSumOfNormals = + derivativeOfNormal_j + derivativeOfNormal_jplus1; + const VectorType sumOfNormals = normal_j + normal_jplus1; + const double normOfSumOfNormals = sumOfNormals.Norm(); + const double derivativeOfNormOfSumOfNormals = + computeDerivativeOfNorm(sumOfNormals, derivativeOfSumOfNormals); + + const VectorType derivativeOfR_firstTerm = -sumOfNormals * + derivativeOfNormOfSumOfNormals / + std::pow(normOfSumOfNormals, 2); + const VectorType derivativeOfR_secondTerm = + derivativeOfSumOfNormals / normOfSumOfNormals; + const VectorType derivativeOfR = + derivativeOfR_firstTerm + derivativeOfR_secondTerm; + + assert(std::isfinite(derivativeOfR[0]) && std::isfinite(derivativeOfR[1]) && + std::isfinite(derivativeOfR[2])); + + return derivativeOfR; +} + +VectorType +FormFinder::computeDerivativeT1(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const { + + const VectorType &X_j = e.cP(0); + const VectorType &X_jplus1 = e.cP(1); + const VectorType edgeVector = X_jplus1 - X_j; + const double L_kDerivative = computeDerivativeElementLength(e, dui); + const double edgeLength = mesh->elements[e].length; + const VectorType firstTerm = + -edgeVector * L_kDerivative / std::pow(edgeLength, 2); + const VectorType secondTerm = + computeDisplacementDifferenceDerivative(e, dui) / edgeLength; + const VectorType t1Derivative = firstTerm + secondTerm; + + return t1Derivative; +} + +VectorType +FormFinder::computeDerivativeT2(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const { + + const DoFType dofi = dui.dofi; + + const VertexType &v_j = *e.cV(0); + const VertexType &v_jplus1 = *e.cV(1); + + const VectorType r = (v_j.cN() + v_jplus1.cN()).Normalize(); + const VectorType derivativeR_j = + dofi > 2 && &v_j == &dui.v ? mesh->elements[e].derivativeR_j[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeR_jplus1 = + dofi > 2 && &v_jplus1 == &dui.v + ? mesh->elements[e].derivativeR_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeR = derivativeR_j + derivativeR_jplus1; + + const VectorType &t1 = mesh->elements[e].frame.t1; + const VectorType derivativeT1_j = + dofi < 3 && &v_j == &dui.v ? mesh->elements[e].derivativeT1_j[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_jplus1 = + dofi < 3 && &v_jplus1 == &dui.v + ? mesh->elements[e].derivativeT1_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1 = derivativeT1_j + derivativeT1_jplus1; + + const VectorType derivativeOfRCrossT1 = + computeDerivativeOfCrossProduct(r, derivativeR, t1, derivativeT1); + const VectorType rCrossT1 = Cross(r, t1); + const double normOfRCrossT1 = rCrossT1.Norm(); + const double derivativeNormRCrossT1 = + computeDerivativeOfNorm(rCrossT1, derivativeOfRCrossT1); + + const VectorType t2Deriv_firstTerm = + -(rCrossT1 * derivativeNormRCrossT1) / std::pow(normOfRCrossT1, 2); + const VectorType t2Deriv_secondTerm = derivativeOfRCrossT1 / normOfRCrossT1; + const VectorType t2Deriv = t2Deriv_firstTerm + t2Deriv_secondTerm; + + const double t2DerivNorm = t2Deriv.Norm(); + assert(std::isfinite(t2DerivNorm)); + return t2Deriv; +} + +VectorType +FormFinder::computeDerivativeT3(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const { + const Element &element = mesh->elements[e]; + const VectorType &t1 = element.frame.t1; + const VectorType &t2 = element.frame.t2; + const VectorType t1CrossT2 = Cross(t1, t2); + const VertexType &v_j = *e.cV(0); + const VertexType &v_jplus1 = *e.cV(1); + const VectorType derivativeT1_j = + dui.dofi < 3 && &v_j == &dui.v + ? mesh->elements[e].derivativeT1_j[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_jplus1 = + dui.dofi < 3 && &v_jplus1 == &dui.v + ? mesh->elements[e].derivativeT1_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1 = derivativeT1_j + derivativeT1_jplus1; + + // const VectorType derivativeOfT2 = computeDerivativeT2(e, dui); + // const VectorType derivativeT2_j = + // &v_j == &dui.v + // ? mesh->elements[e].derivativeT2_j[dui.dofi] + // : VectorType(0, 0, 0); + // const VectorType derivativeT2_jplus1 = + // &v_jplus1 == &dui.v + // ? mesh->elements[e].derivativeT2_jplus1[dui.dofi] + // : VectorType(0, 0, 0); + VectorType derivativeT2(0, 0, 0); + if (&v_j == &dui.v) { + derivativeT2 = mesh->elements[e].derivativeT2_j[dui.dofi]; + } else if (&v_jplus1 == &dui.v) { + derivativeT2 = mesh->elements[e].derivativeT2_jplus1[dui.dofi]; + } + + const VectorType derivativeT1CrossT2 = + computeDerivativeOfCrossProduct(t1, derivativeT1, t2, derivativeT2); + const double derivativeOfNormT1CrossT2 = + computeDerivativeOfNorm(t1CrossT2, derivativeT1CrossT2); + const double normT1CrossT2 = t1CrossT2.Norm(); + + const VectorType t3Deriv_firstTerm = + -(t1CrossT2 * derivativeOfNormT1CrossT2) / std::pow(normT1CrossT2, 2); + const VectorType t3Deriv_secondTerm = derivativeT1CrossT2 / normT1CrossT2; + const VectorType t3Deriv = t3Deriv_firstTerm + t3Deriv_secondTerm; + + assert(std::isfinite(t3Deriv[0]) && std::isfinite(t3Deriv[1]) && + std::isfinite(t3Deriv[2])); + return t3Deriv; +} + +double FormFinder::computeDerivativeTheta1(const EdgeType &e, + const VertexIndex &evi, + const VertexIndex &dwrt_evi, + const DoFType &dwrt_dofi) const { + const VertexType &v = *e.cV(evi); + const Element &element = mesh->elements[e]; + const VectorType derivativeT1 = element.derivativeT1[dwrt_evi][dwrt_dofi]; + const VectorType derivativeT3 = element.derivativeT3[dwrt_evi][dwrt_dofi]; + const VectorType nDerivative = + evi != dwrt_evi ? VectorType(0, 0, 0) + : mesh->nodes[v].derivativeOfNormal[dwrt_dofi]; + const VectorType n = v.cN(); + const VectorType &t1 = element.frame.t1; + const VectorType &t3 = element.frame.t3; + const double theta1Derivative = + derivativeT1 * Cross(t3, n) + + t1 * (Cross(derivativeT3, n) + Cross(t3, nDerivative)); + + return theta1Derivative; +} + +double FormFinder::computeDerivativeTheta2(const EdgeType &e, + const VertexIndex &evi, + const VertexIndex &dwrt_evi, + const DoFType &dwrt_dofi) const { + const VertexType &v = *e.cV(evi); + + const Element &element = mesh->elements[e]; + const VectorType derivativeT2 = element.derivativeT2[dwrt_evi][dwrt_dofi]; + const VectorType derivativeT3 = element.derivativeT3[dwrt_evi][dwrt_dofi]; + + const VectorType n = v.cN(); + const VectorType &t2 = element.frame.t2; + const VectorType &t3 = element.frame.t3; + const VectorType nDerivative = + dwrt_evi == evi ? mesh->nodes[v].derivativeOfNormal[dwrt_dofi] + : VectorType(0, 0, 0); + const double theta2Derivative = + derivativeT2 * Cross(t3, n) + + t2 * (Cross(derivativeT3, n) + Cross(t3, nDerivative)); + + return theta2Derivative; +} + +double FormFinder::computeTheta3(const EdgeType &e, const VertexType &v) { + const VertexIndex &vi = mesh->nodes[v].vi; + // assert(e.cV(0) == &v || e.cV(1) == &v); + + const Element &elem = mesh->elements[e]; + const EdgeIndex &ei = elem.ei; + const Element::LocalFrame &ef = elem.frame; + const VectorType &t1 = ef.t1; + const VectorType &n = v.cN(); + const Node &node = mesh->nodes[v]; + const double &nR = node.nR; // TODO: This makes the function non-const. + // Should be refactored in the future + + double theta3; + if (&e == node.referenceElement) { + // Use nR as theta3 only for the first star edge + return nR; + } + // assert(node.alphaAngles.contains(mesh->getIndex(e))); + double alphaAngle = node.alphaAngles.find(elem.ei)->second; + const EdgeType &refElem = *node.referenceElement; + const double rotationAngle = nR + alphaAngle; + + // const VectorType &t1_k = computeT1Vector(refElem); + const VectorType &t1_k = mesh->elements[refElem].frame.t1; + const VectorType f1 = (t1_k - (n * (t1_k * n))).Normalize(); + vcg::Matrix44 rotationMatrix; + rotationMatrix.SetRotateRad(rotationAngle, n); + const double cosRotationAngle = cos(rotationAngle); + const double sinRotationAngle = sin(rotationAngle); + const VectorType f2 = + (f1 * cosRotationAngle + Cross(n, f1) * sinRotationAngle).Normalize(); + const VectorType &t1Current = t1; + const VectorType f3 = Cross(t1Current, f2); + + Element &element = mesh->elements[e]; + // Save for computing theta3Derivative + if (&v == e.cV(0)) { + element.f1_j = f1; + element.f2_j = f2; + element.f3_j = f3; + element.cosRotationAngle_j = cosRotationAngle; + element.sinRotationAngle_j = sinRotationAngle; + } else { + element.f1_jplus1 = f1; + element.f2_jplus1 = f2; + element.f3_jplus1 = f3; + element.cosRotationAngle_jplus1 = cosRotationAngle; + element.sinRotationAngle_jplus1 = sinRotationAngle; + } + theta3 = f3.dot(n); + + return theta3; +} + +double FormFinder::computeDerivativeTheta3( + const EdgeType &e, const VertexType &v, + const DifferentiateWithRespectTo &dui) const { + const Node &node = mesh->nodes[v]; + const VertexIndex &vi = mesh->nodes[v].vi; + const bool isRigidSupport = rigidSupports.contains(vi); + if (&e == node.referenceElement && !isRigidSupport) { + if (dui.dofi == DoF::Nr && &dui.v == &v) { + return 1; + } else { + return 0; + } + } + // assert(e.cV(0) == &v || e.cV(1) == &v); + + const Element &element = mesh->elements[e]; + const Element::LocalFrame &ef = element.frame; + const VectorType &t1 = ef.t1; + const VectorType &n = v.cN(); + const DoFType &dofi = dui.dofi; + const VertexPointer &vp_j = e.cV(0); + const VertexPointer &vp_jplus1 = e.cV(1); + + double derivativeTheta3_dofi = 0; + if (isRigidSupport) { + const VectorType &t1Initial = + computeT1Vector(mesh->nodes[vp_j].initialLocation, + mesh->nodes[vp_jplus1].initialLocation); + VectorType g1 = Cross(t1, t1Initial); + + const VectorType derivativeInitialT1_dofi(0, 0, 0); + const VectorType derivativeT1_j = dui.dofi < 3 && vp_j == &dui.v + ? element.derivativeT1_j[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_jplus1 = + dui.dofi < 3 && vp_jplus1 == &dui.v + ? element.derivativeT1_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_dofi = derivativeT1_j + derivativeT1_jplus1; + // VectorType derivativeT1_dofi(0, 0, 0); + // if (dui.dofi < 3) { + // if (&v_j == &dui.v) { + // derivativeT1_dofi = mesh->elements[e].derivativeT1_j[dui.dofi]; + // } else if (&v_jplus1 == &dui.v) { + // derivativeT1_dofi = + // mesh->elements[e].derivativeT1_jplus1[dui.dofi]; + // } + // } + + const VectorType derivativeG1_firstTerm = + Cross(derivativeT1_dofi, t1Initial); + const VectorType derivativeG1_secondTerm = + Cross(t1, derivativeInitialT1_dofi); + const VectorType derivativeG1 = + derivativeG1_firstTerm + derivativeG1_secondTerm; + const VectorType derivativeNormal = &v == &dui.v && dui.dofi > 2 + ? node.derivativeOfNormal[dui.dofi] + : VectorType(0, 0, 0); + derivativeTheta3_dofi = derivativeG1 * n + g1 * derivativeNormal; + return derivativeTheta3_dofi; + } + EdgeType &refElem = *node.referenceElement; + const VectorType &t1_k = mesh->elements[refElem].frame.t1; + VectorType f1, f2, f3; + double cosRotationAngle, sinRotationAngle; + if (e.cV(0) == &v) { + f1 = element.f1_j; + cosRotationAngle = element.cosRotationAngle_j; + sinRotationAngle = element.sinRotationAngle_j; + f2 = element.f2_j; + f3 = element.f3_j; + } else { + f1 = element.f1_jplus1; + cosRotationAngle = element.cosRotationAngle_jplus1; + sinRotationAngle = element.sinRotationAngle_jplus1; + f2 = element.f2_jplus1; + f3 = element.f3_jplus1; + } + const VectorType &t1_kplus1 = t1; + const VectorType f1Normalized = f1 / f1.Norm(); + + VectorType derivativeF1(0, 0, 0); + VectorType derivativeF2(0, 0, 0); + VectorType derivativeF3(0, 0, 0); + if (dui.dofi < 3) { + + const VectorType derivativeT1_kplus1_j = + vp_j == &dui.v ? element.derivativeT1_j[dui.dofi] : VectorType(0, 0, 0); + const VectorType derivativeT1_kplus1_jplus1 = + vp_jplus1 == &dui.v ? element.derivativeT1_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_kplus1 = + derivativeT1_kplus1_j + derivativeT1_kplus1_jplus1; + + const VectorType derivativeT1_k_j = + refElem.cV(0) == &dui.v + ? mesh->elements[refElem].derivativeT1_j[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_k_jplus1 = + refElem.cV(1) == &dui.v + ? mesh->elements[refElem].derivativeT1_jplus1[dui.dofi] + : VectorType(0, 0, 0); + const VectorType derivativeT1_k = derivativeT1_k_j + derivativeT1_k_jplus1; + + derivativeF1 = derivativeT1_k - (n * (derivativeT1_k * n)); + + const double f1Norm = f1.Norm(); + const VectorType derivativeF1Normalized = + -f1 * (f1 * derivativeF1) / (f1Norm * f1Norm * f1Norm) + + derivativeF1 / f1Norm; + + derivativeF2 = derivativeF1Normalized * cosRotationAngle + + Cross(n, derivativeF1Normalized) * sinRotationAngle; + const VectorType derivativeF3_firstTerm = Cross(derivativeT1_kplus1, f2); + const VectorType derivativeF3_secondTerm = Cross(t1_kplus1, derivativeF2); + derivativeF3 = derivativeF3_firstTerm + derivativeF3_secondTerm; + derivativeTheta3_dofi = derivativeF3 * n; + + } else if (dui.dofi == DoF::Nr && &dui.v == &v) { + derivativeF2 = f1Normalized * (-sinRotationAngle) + + Cross(n, f1Normalized) * cosRotationAngle; + derivativeF3 = Cross(t1_kplus1, derivativeF2); + derivativeTheta3_dofi = derivativeF3 * n; + } else { // 2vert) { + const Node &node = mesh->nodes[v]; + if (!node.force.hasExternalForce()) { + continue; + } + const auto nodeDisplacement = v.cP() - node.initialLocation; + const SimulationMesh::CoordType externalForce( + node.force.external[0], node.force.external[1], node.force.external[2]); + totalExternalPotentialEnergy += externalForce.dot(nodeDisplacement); + } + + double totalInternalPotentialEnergy = 0; + for (const SimulationMesh::EdgeType &e : mesh->edge) { + + const Element &element = mesh->elements[e]; + const EdgeIndex ei = mesh->getIndex(e); + const double e_k = element.length - element.initialLength; + const double axialPE = pow(e_k, 2) * element.properties.E * + element.properties.A / (2 * element.initialLength); + const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 - + element.rotationalDisplacements_j.theta1; + const double torsionalPE = element.properties.G * element.properties.J * + pow(theta1Diff, 2) / (2 * element.initialLength); + const double &theta2_j = element.rotationalDisplacements_j.theta2; + const double &theta2_jplus1 = element.rotationalDisplacements_jplus1.theta2; + 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.properties.E * element.properties.I2 / + (2 * element.initialLength); + const double &theta3_j = element.rotationalDisplacements_j.theta3; + const double &theta3_jplus1 = element.rotationalDisplacements_jplus1.theta3; + const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) + + 4 * theta3_j * theta3_jplus1 + + 4 * pow(theta3_jplus1, 2); + const double secondBendingPE = + secondBendingPE_inBracketsTerm * 2 * element.properties.E * + element.properties.I3 / element.initialLength; + + totalInternalPotentialEnergy += + axialPE + torsionalPE + firstBendingPE + secondBendingPE; + } + return totalInternalPotentialEnergy - totalExternalPotentialEnergy; +} + +void FormFinder::updateResidualForcesOnTheFly( + const std::unordered_map> + &fixedVertices) { + + // std::vector internalForcesParallel(mesh->vert.size()); + + std::vector>> + internalForcesContributionsFromEachEdge( + mesh->EN(), + std::vector>(4, {-1, Vector6d()})); + // omp_lock_t writelock; + // omp_init_lock(&writelock); +#pragma omp parallel for schedule(static) num_threads(8) + for (size_t ei = 0; ei < mesh->EN(); ei++) { + const EdgeType &e = mesh->edge[ei]; + const SimulationMesh::VertexType &ev_j = *e.cV(0); + const SimulationMesh::VertexType &ev_jplus1 = *e.cV(1); + const Element &element = mesh->elements[e]; + const Element::LocalFrame &ef = element.frame; + const VectorType t3CrossN_j = Cross(ef.t3, ev_j.cN()); + const VectorType t3CrossN_jplus1 = Cross(ef.t3, ev_jplus1.cN()); + const double theta1_j = ef.t1.dot(t3CrossN_j); + const double theta1_jplus1 = ef.t1.dot(t3CrossN_jplus1); + const double theta2_j = ef.t2.dot(t3CrossN_j); + const double theta2_jplus1 = ef.t2.dot(t3CrossN_jplus1); + const double theta3_j = computeTheta3(e, ev_j); + const double theta3_jplus1 = computeTheta3(e, ev_jplus1); + // element.rotationalDisplacements_j.theta1 = theta1_j; + // element.rotationalDisplacements_jplus1.theta1 = theta1_jplus1; + // element.rotationalDisplacements_j.theta2 = theta2_j; + // element.rotationalDisplacements_jplus1.theta2 = theta2_jplus1; + // element.rotationalDisplacements_j.theta3 = theta3_j; + // element.rotationalDisplacements_jplus1.theta3 = theta3_jplus1; + std::vector> + internalForcesContributionFromThisEdge(4, {-1, Vector6d()}); + for (VertexIndex evi = 0; evi < 2; evi++) { + const SimulationMesh::VertexType &ev = *e.cV(evi); + const Node &edgeNode = mesh->nodes[ev]; + internalForcesContributionFromThisEdge[evi].first = edgeNode.vi; + + const VertexPointer &rev_j = edgeNode.referenceElement->cV(0); + const VertexPointer &rev_jplus1 = edgeNode.referenceElement->cV(1); + // refElemOtherVertex can be precomputed + const VertexPointer &refElemOtherVertex = + rev_j == &ev ? rev_jplus1 : rev_j; + const Node &refElemOtherVertexNode = mesh->nodes[refElemOtherVertex]; + if (edgeNode.referenceElement != &e) { + internalForcesContributionFromThisEdge[evi + 2].first = + refElemOtherVertexNode.vi; + } + // #pragma omp parallel for schedule(static) num_threads(6) + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + const bool isDofConstrainedFor_ev = + fixedVertices.contains(edgeNode.vi) && + fixedVertices.at(edgeNode.vi).contains(dofi); + if (!isDofConstrainedFor_ev) { + DifferentiateWithRespectTo dui{ev, dofi}; + // Axial force computation + const double e_k = element.length - element.initialLength; + const double e_kDeriv = computeDerivativeElementLength(e, dui); + const double axialForce_dofi = + e_kDeriv * e_k * element.axialConstFactor; + + // Torsional force computation + const double theta1_j_deriv = + computeDerivativeTheta1(e, 0, evi, dofi); + const double theta1_jplus1_deriv = + computeDerivativeTheta1(e, 1, evi, dofi); + const double theta1Diff = theta1_jplus1 - theta1_j; + const double theta1DiffDerivative = + theta1_jplus1_deriv - theta1_j_deriv; + const double torsionalForce_dofi = + theta1DiffDerivative * theta1Diff * element.torsionConstFactor; + + // First bending force computation + ////theta2_j derivative + const double theta2_j_deriv = + computeDerivativeTheta2(e, 0, evi, dofi); + ////theta2_jplus1 derivative + const double theta2_jplus1_deriv = + computeDerivativeTheta2(e, 1, evi, dofi); + ////1st in bracket term + const double firstBendingForce_inBracketsTerm_0 = + theta2_j_deriv * 2 * theta2_j; + ////2nd in bracket term + const double firstBendingForce_inBracketsTerm_1 = + theta2_jplus1_deriv * theta2_j; + ////3rd in bracket term + const double firstBendingForce_inBracketsTerm_2 = + theta2_j_deriv * theta2_jplus1; + ////4th in bracket term + const double firstBendingForce_inBracketsTerm_3 = + 2 * theta2_jplus1_deriv * theta2_jplus1; + // 3rd term computation + const double firstBendingForce_inBracketsTerm = + firstBendingForce_inBracketsTerm_0 + + firstBendingForce_inBracketsTerm_1 + + firstBendingForce_inBracketsTerm_2 + + firstBendingForce_inBracketsTerm_3; + const double firstBendingForce_dofi = + firstBendingForce_inBracketsTerm * + element.firstBendingConstFactor; + + // Second bending force computation + ////theta2_j derivative + const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui); + ////theta2_jplus1 derivative + const double theta3_jplus1_deriv = + computeDerivativeTheta3(e, ev_jplus1, dui); + ////1st in bracket term + const double secondBendingForce_inBracketsTerm_0 = + theta3_j_deriv * 2 * theta3_j; + ////2nd in bracket term + const double secondBendingForce_inBracketsTerm_1 = + theta3_jplus1_deriv * theta3_j; + ////3rd in bracket term + const double secondBendingForce_inBracketsTerm_2 = + theta3_j_deriv * theta3_jplus1; + ////4th in bracket term + const double secondBendingForce_inBracketsTerm_3 = + 2 * theta3_jplus1_deriv * theta3_jplus1; + // 3rd term computation + const double secondBendingForce_inBracketsTerm = + secondBendingForce_inBracketsTerm_0 + + secondBendingForce_inBracketsTerm_1 + + secondBendingForce_inBracketsTerm_2 + + secondBendingForce_inBracketsTerm_3; + double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * + element.secondBendingConstFactor; + internalForcesContributionFromThisEdge[evi].second[dofi] = + axialForce_dofi + firstBendingForce_dofi + + secondBendingForce_dofi + torsionalForce_dofi; + } + if (edgeNode.referenceElement != &e) { + const bool isDofConstrainedFor_refElemOtherVertex = + fixedVertices.contains(refElemOtherVertexNode.vi) && + fixedVertices.at(refElemOtherVertexNode.vi).contains(dofi); + if (!isDofConstrainedFor_refElemOtherVertex) { + DifferentiateWithRespectTo dui{*refElemOtherVertex, dofi}; + ////theta3_j derivative + const double theta3_j_deriv = computeDerivativeTheta3(e, ev_j, dui); + ////theta3_jplus1 derivative + const double theta3_jplus1_deriv = + computeDerivativeTheta3(e, ev_jplus1, dui); + ////1st in bracket term + const double secondBendingForce_inBracketsTerm_0 = + theta3_j_deriv * 2 * theta3_j; + ////2nd in bracket term + const double secondBendingForce_inBracketsTerm_1 = + theta3_jplus1_deriv * theta3_j; + ////3rd in bracket term + const double secondBendingForce_inBracketsTerm_2 = + theta3_j_deriv * theta3_jplus1; + ////4th in bracket term + const double secondBendingForce_inBracketsTerm_3 = + theta3_jplus1_deriv * 2 * theta3_jplus1; + + // 4th term computation + const double secondBendingForce_inBracketsTerm = + secondBendingForce_inBracketsTerm_0 + + secondBendingForce_inBracketsTerm_1 + + secondBendingForce_inBracketsTerm_2 + + secondBendingForce_inBracketsTerm_3; + const double secondBendingForce_dofi = + secondBendingForce_inBracketsTerm * + element.secondBendingConstFactor; + internalForcesContributionFromThisEdge[evi + 2].second[dofi] = + secondBendingForce_dofi; + } + } + } + } + internalForcesContributionsFromEachEdge[ei] = + internalForcesContributionFromThisEdge; + } + + //#pragma omp parallel for schedule(static) num_threads(8) + + for (size_t vi = 0; vi < mesh->VN(); vi++) { + Node::Forces &force = mesh->nodes[vi].force; + force.residual = force.external; + force.internal = 0; + } + double totalResidualForcesNorm = 0; + for (size_t ei = 0; ei < mesh->EN(); ei++) { + for (int i = 0; i < 4; i++) { + std::pair internalForcePair = + internalForcesContributionsFromEachEdge[ei][i]; + int vi = internalForcePair.first; + if (i > 1 && vi == -1) { + continue; + } + Node::Forces &force = mesh->nodes[vi].force; + force.internal = force.internal + internalForcePair.second; + force.residual = force.residual + (internalForcePair.second * -1); + } + } + for (size_t vi = 0; vi < mesh->VN(); vi++) { + const double residualForceNorm = (mesh->nodes[vi].force.residual).norm(); + totalResidualForcesNorm += residualForceNorm; + } + mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm; + mesh->totalResidualForcesNorm = totalResidualForcesNorm; +} + +void FormFinder::updateNodalExternalForces( + const std::unordered_map &nodalForces, + const std::unordered_map> + &fixedVertices) { + + for (const std::pair &nodalForce : nodalForces) { + const VertexIndex nodeIndex = nodalForce.first; + const bool isNodeConstrained = fixedVertices.contains(nodeIndex); + Node &node = mesh->nodes[nodeIndex]; + Vector6d nodalExternalForce(0); + for (DoFType dofi = 0; dofi < 6; dofi++) { + const bool isDofConstrained = + isNodeConstrained && fixedVertices.at(nodeIndex).contains(dofi); + if (isDofConstrained) { + continue; + } + nodalExternalForce[dofi] = nodalForce.second[dofi]; + } + node.force.external = nodalExternalForce; + } +} + +void FormFinder::updateResidualForces() { + + mesh->totalResidualForcesNorm = 0; + for (const VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.force.residual = node.force.external - node.force.internal; + const double residualForceNorm = (node.force.residual).norm(); + mesh->totalResidualForcesNorm += residualForceNorm; + } +} + +void FormFinder::computeRigidSupports() { + + for (const VertexType &v : mesh->vert) { + const VertexIndex vi = mesh->nodes[v].vi; + const bool isVertexConstrained = constrainedVertices.contains(vi); + if (isVertexConstrained) { + auto constrainedDoFType = constrainedVertices.at(vi); + const bool hasAllDoFTypeConstrained = + constrainedDoFType.contains(DoF::Ux) && + constrainedDoFType.contains(DoF::Uy) && + constrainedDoFType.contains(DoF::Uz) && + constrainedDoFType.contains(DoF::Nx) && + constrainedDoFType.contains(DoF::Ny) && + constrainedDoFType.contains(DoF::Nr); + if (hasAllDoFTypeConstrained) { + rigidSupports.insert(vi); + } + } + } +} + +void FormFinder::updateNormalDerivatives() { + for (const VertexType &v : mesh->vert) { + for (DoFType dofi = DoF::Nx; dofi < DoF::NumDoF; dofi++) { + DifferentiateWithRespectTo dui{v, dofi}; + mesh->nodes[v].derivativeOfNormal[dofi] = + computeDerivativeOfNormal(v, dui); + } + } +} + +void FormFinder::updateT1Derivatives() { + for (const EdgeType &e : mesh->edge) { + for (DoFType dofi = DoF::Ux; dofi < DoF::Nx; dofi++) { + DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; + Element &element = mesh->elements[e]; + element.derivativeT1_j[dofi] = computeDerivativeT1(e, dui_v0); + DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; + element.derivativeT1_jplus1[dofi] = computeDerivativeT1(e, dui_v1); + + element.derivativeT1[0][dofi] = element.derivativeT1_j[dofi]; + element.derivativeT1[1][dofi] = element.derivativeT1_jplus1[dofi]; + } + } +} + +void FormFinder::updateT2Derivatives() { + for (const EdgeType &e : mesh->edge) { + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; + mesh->elements[e].derivativeT2_j[dofi] = computeDerivativeT2(e, dui_v0); + DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; + mesh->elements[e].derivativeT2_jplus1[dofi] = + computeDerivativeT2(e, dui_v1); + + Element &element = mesh->elements[e]; + element.derivativeT2[0][dofi] = element.derivativeT2_j[dofi]; + element.derivativeT2[1][dofi] = element.derivativeT2_jplus1[dofi]; + } + } +} + +void FormFinder::updateT3Derivatives() { + for (const EdgeType &e : mesh->edge) { + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; + Element &element = mesh->elements[e]; + element.derivativeT3_j[dofi] = computeDerivativeT3(e, dui_v0); + DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; + element.derivativeT3_jplus1[dofi] = computeDerivativeT3(e, dui_v1); + + element.derivativeT3[0][dofi] = element.derivativeT3_j[dofi]; + element.derivativeT3[1][dofi] = element.derivativeT3_jplus1[dofi]; + } + } +} + +void FormFinder::updateRDerivatives() { + for (const EdgeType &e : mesh->edge) { + for (DoFType dofi = DoF::Nx; dofi < DoF::NumDoF; dofi++) { + DifferentiateWithRespectTo dui_v0{*e.cV(0), dofi}; + mesh->elements[e].derivativeR_j[dofi] = computeDerivativeOfR(e, dui_v0); + DifferentiateWithRespectTo dui_v1{*e.cV(1), dofi}; + mesh->elements[e].derivativeR_jplus1[dofi] = + computeDerivativeOfR(e, dui_v1); + } + } +} + +void FormFinder::updateElementalLengths() { mesh->updateElementalLengths(); } + +FormFinder::FormFinder() {} + +void FormFinder::updateNodalMasses() { + const double gamma = 0.8; + for (VertexType &v : mesh->vert) { + double translationalSumSk = 0; + double rotationalSumSk_I2 = 0; + double rotationalSumSk_I3 = 0; + double rotationalSumSk_J = 0; + for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { + const size_t ei = mesh->getIndex(ep); + const Element &elem = mesh->elements[ep]; + const double SkTranslational = + elem.properties.E * elem.properties.A / elem.length; + translationalSumSk += SkTranslational; + const double lengthToThe3 = std::pow(elem.length, 3); + const double SkRotational_I2 = + elem.properties.E * elem.properties.I2 / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const double SkRotational_I3 = + elem.properties.E * elem.properties.I3 / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const double SkRotational_J = + elem.properties.E * elem.properties.J / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + rotationalSumSk_I2 += SkRotational_I2; + rotationalSumSk_I3 += SkRotational_I3; + rotationalSumSk_J += SkRotational_J; + assert(rotationalSumSk_I2 != 0); + assert(rotationalSumSk_I3 != 0); + assert(rotationalSumSk_J != 0); + } + mesh->nodes[v].translationalMass = + gamma * pow(Dtini, 2) * 2 * translationalSumSk; + mesh->nodes[v].rotationalMass_I2 = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I2; + mesh->nodes[v].rotationalMass_I3 = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I3; + mesh->nodes[v].rotationalMass_J = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_J; + + assert(std::pow(Dtini, 2.0) * translationalSumSk / + mesh->nodes[v].translationalMass < + 2); + assert(std::pow(Dtini, 2.0) * rotationalSumSk_I2 / + mesh->nodes[v].rotationalMass_I2 < + 2); + assert(std::pow(Dtini, 2.0) * rotationalSumSk_I3 / + mesh->nodes[v].rotationalMass_I3 < + 2); + assert(std::pow(Dtini, 2.0) * rotationalSumSk_J / + mesh->nodes[v].rotationalMass_J < + 2); + } +} + +void FormFinder::updateNodalAccelerations() { + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.translationalMass; + } else if (dofi == DoF::Nx) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_J; + } else if (dofi == DoF::Ny) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_I2; + } else if (dofi == DoF::Nr) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_I3; + } + } + } +} + +void FormFinder::updateNodalVelocities() { + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.velocity = node.velocity + node.acceleration * Dt; + } + updateKineticEnergy(); +} + +void FormFinder::updateNodalDisplacements() { + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.displacements = node.displacements + node.velocity * Dt; + } +} + +void FormFinder::updateNodePosition( + VertexType &v, + const std::unordered_map> + &fixedVertices) { + Node &node = mesh->nodes[v]; + CoordType previousLocation = v.cP(); + const VertexIndex &vi = mesh->nodes[v].vi; + + VectorType displacementVector(0, 0, 0); + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(0)) { + displacementVector += VectorType(node.displacements[0], 0, 0); + } + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(1)) { + displacementVector += VectorType(0, node.displacements[1], 0); + } + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(2)) { + displacementVector += VectorType(0, 0, node.displacements[2]); + } + node.previousLocation = previousLocation; + v.P() = node.initialLocation + displacementVector; + if (shouldApplyInitialDistortion && currentSimulationStep < 40) { + VectorType desiredInitialDisplacement(0, 0, 0.1); + v.P() = v.P() + desiredInitialDisplacement; + } +} + +void FormFinder::applyDisplacements( + const std::unordered_map> + &fixedVertices) { + for (VertexType &v : mesh->vert) { + updateNodePosition(v, fixedVertices); + Node &node = mesh->nodes[v]; + const VertexIndex &vi = node.vi; + VectorType normalDisplacementVector(0, 0, 0); + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(3)) { + normalDisplacementVector += VectorType(node.displacements[3], 0, 0); + } + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(4)) { + normalDisplacementVector += VectorType(0, node.displacements[4], 0); + } + v.N() = node.initialNormal + normalDisplacementVector; + const double &nx = v.N()[0], ny = v.N()[1]; + const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); + if (nxnyMagnitude > 1) { + v.N() = VectorType(nx / std::sqrt(nxnyMagnitude), + ny / std::sqrt(nxnyMagnitude), 0); + } else { + const double nzSquared = 1.0 - nxnyMagnitude; + const double nz = std::sqrt(nzSquared); + VectorType newNormal(nx, ny, nz); + v.N() = newNormal; + } + if (!rigidSupports.contains(vi)) { + node.nR = node.displacements[5]; + } else { + const EdgePointer &refElem = node.referenceElement; + const VectorType &refT1 = mesh->elements[refElem].frame.t1; + + const VectorType &t1Initial = + computeT1Vector(mesh->nodes[refElem->cV(0)].initialLocation, + mesh->nodes[refElem->cV(1)].initialLocation); + VectorType g1 = Cross(refT1, t1Initial); + node.nR = g1.dot(v.cN()); + } + } + updateElementalFrames(); +} + +void FormFinder::updateElementalFrames() { + for (const EdgeType &e : mesh->edge) { + const VectorType elementNormal = + (e.cV(0)->cN() + e.cV(1)->cN()).Normalize(); + mesh->elements[e].frame = + computeElementFrame(e.cP(0), e.cP(1), elementNormal); + } +} + +void FormFinder::applyForcedDisplacements( + const std::unordered_map + nodalForcedDisplacements) { + for (const std::pair + vertexIndexDisplacementPair : nodalForcedDisplacements) { + const VertexIndex vi = vertexIndexDisplacementPair.first; + const Eigen::Vector3d vertexDisplacement = + vertexIndexDisplacementPair.second; + Node &node = mesh->nodes[vi]; + VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1), + vertexDisplacement(2)); + mesh->vert[vi].P() = node.initialLocation + displacementVector; + } +} + +// NOTE: Is this correct? Should the kinetic energy be computed like that? +void FormFinder::updateKineticEnergy() { + mesh->previousTotalKineticEnergy = mesh->currentTotalKineticEnergy; + mesh->currentTotalKineticEnergy = 0; + for (const VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.kineticEnergy = 0; + + const double translationalVelocityNorm = std::sqrt( + std::pow(node.velocity[0], 2) + std::pow(node.velocity[1], 2) + + std::pow(node.velocity[2], 2)); + node.kineticEnergy += + 0.5 * node.translationalMass * pow(translationalVelocityNorm, 2); + + assert(node.kineticEnergy < 10000000000000000000); + // const double rotationalVelocityNorm = std::sqrt( + // std::pow(node.velocity[3], 2) + std::pow(node.velocity[4], 2) + + // std::pow(node.velocity[5], 2)); + node.kineticEnergy += + 0.5 * node.rotationalMass_J * pow(node.velocity[3], 2) + + 0.5 * node.rotationalMass_I3 * pow(node.velocity[4], 2) + + 0.5 * node.rotationalMass_I2 * pow(node.velocity[5], 2); + + mesh->currentTotalKineticEnergy += node.kineticEnergy; + } + // assert(mesh->currentTotalKineticEnergy < 100000000000000); +} + +void FormFinder::resetVelocities() { + for (const VertexType &v : mesh->vert) { + mesh->nodes[v].velocity = 0; + // mesh->nodes[v].force.residual * 0.5 * Dt / + // mesh->nodes[v].mass; // NOTE: Do I reset the previous + // velocity? + // reset + // current to 0 or this? + } + updateKineticEnergy(); +} + +void FormFinder::updatePositionsOnTheFly( + const std::unordered_map> + &fixedVertices) { + const double gamma = 0.8; + for (VertexType &v : mesh->vert) { + double translationalSumSk = 0; + double rotationalSumSk_I2 = 0; + double rotationalSumSk_I3 = 0; + double rotationalSumSk_J = 0; + for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { + const Element &elem = mesh->elements[ep]; + const double SkTranslational = + elem.properties.E * elem.properties.A / elem.length; + translationalSumSk += SkTranslational; + const double lengthToThe3 = std::pow(elem.length, 3); + const double SkRotational_I2 = + elem.properties.E * elem.properties.I2 / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const double SkRotational_I3 = + elem.properties.E * elem.properties.I3 / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const double SkRotational_J = + elem.properties.E * elem.properties.J / + lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + rotationalSumSk_I2 += SkRotational_I2; + rotationalSumSk_I3 += SkRotational_I3; + rotationalSumSk_J += SkRotational_J; + // assert(rotationalSumSk_I2 != 0); + // assert(rotationalSumSk_I3 != 0); + // assert(rotationalSumSk_J != 0); + } + mesh->nodes[v].translationalMass = + gamma * pow(Dtini, 2) * 2 * translationalSumSk; + mesh->nodes[v].rotationalMass_I2 = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I2; + mesh->nodes[v].rotationalMass_I3 = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_I3; + mesh->nodes[v].rotationalMass_J = + gamma * pow(Dtini, 2) * 8 * rotationalSumSk_J; + + // assert(std::pow(Dtini, 2.0) * translationalSumSk / + // mesh->nodes[v].translationalMass < + // 2); + // assert(std::pow(Dtini, 2.0) * rotationalSumSk_I2 / + // mesh->nodes[v].rotationalMass_I2 < + // 2); + // assert(std::pow(Dtini, 2.0) * rotationalSumSk_I3 / + // mesh->nodes[v].rotationalMass_I3 < + // 2); + // assert(std::pow(Dtini, 2.0) * rotationalSumSk_J / + // mesh->nodes[v].rotationalMass_J < + // 2); + } + + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { + if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.translationalMass; + } else if (dofi == DoF::Nx) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_J; + } else if (dofi == DoF::Ny) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_I3; + } else if (dofi == DoF::Nr) { + node.acceleration[dofi] = + node.force.residual[dofi] / node.rotationalMass_I2; + } + } + } + + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.velocity = node.velocity + node.acceleration * Dt; + } + updateKineticEnergy(); + + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.displacements = node.displacements + node.velocity * Dt; + } + + for (VertexType &v : mesh->vert) { + updateNodePosition(v, fixedVertices); + Node &node = mesh->nodes[v]; + const VertexIndex &vi = node.vi; + VectorType normalDisplacementVector(0, 0, 0); + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(3)) { + normalDisplacementVector += VectorType(node.displacements[3], 0, 0); + } + if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(4)) { + normalDisplacementVector += VectorType(0, node.displacements[4], 0); + } + v.N() = node.initialNormal + normalDisplacementVector; + const double &nx = v.N()[0], ny = v.N()[1]; + const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); + if (nxnyMagnitude > 1) { + v.N() = VectorType(nx / std::sqrt(nxnyMagnitude), + ny / std::sqrt(nxnyMagnitude), 0); + } else { + const double nzSquared = 1.0 - nxnyMagnitude; + const double nz = std::sqrt(nzSquared); + VectorType newNormal(nx, ny, nz); + v.N() = newNormal; + } + if (!rigidSupports.contains(vi)) { + node.nR = node.displacements[5]; + } else { + } + } + updateElementalFrames(); +} + +SimulationResults +FormFinder::computeResults(const SimulationMesh &initialMesh) { + const size_t numberOfVertices = initialMesh.VN(); + std::vector displacements(numberOfVertices); + + for (size_t vi = 0; vi < numberOfVertices; vi++) { + displacements[vi] = mesh->nodes[vi].displacements; + } + + history.numberOfSteps = currentSimulationStep; + + return SimulationResults{history, displacements}; +} + +void FormFinder::draw(const std::string &screenshotsFolder = {}) { + // update positions + // polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->updateNodePositions(mesh->getEigenVertices()); + // Per node kinetic energies + std::vector kineticEnergies(mesh->VN()); + for (const VertexType &v : mesh->vert) { + kineticEnergies[mesh->getIndex(v)] = mesh->nodes[v].kineticEnergy; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("Kinetic Energy", kineticEnergies); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Kinetic Energy") + ->setEnabled(false); + + // Per node normals + std::vector> nodeNormals(mesh->VN()); + for (const VertexType &v : mesh->vert) { + const VectorType n = v.cN(); + nodeNormals[mesh->getIndex(v)] = {n[0], n[1], n[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Node normals", nodeNormals) + ->setEnabled(true); + + // per node external forces + // std::vector> externalForces(mesh->VN()); + // for (const VertexType &v : mesh->vert) { + // const Vector6d nodeForce = + // mesh->nodes[v].force.external.value_or(Vector6d(0)); + // externalForces[mesh->getIndex(v)] = {nodeForce[0], nodeForce[1], + // nodeForce[2]}; + // } + // polyscope::getCurveNetwork("Deformed edge mesh") + // ->addNodeVectorQuantity("External force", externalForces); + // polyscope::getCurveNetwork("Deformed edge mesh") + // ->getQuantity("External force") + // ->setEnabled(true); + + std::vector> internalForces(mesh->VN()); + std::vector> externalForces(mesh->VN()); + std::vector internalForcesNorm(mesh->VN()); + for (const VertexType &v : mesh->vert) { + // per node internal forces + const Vector6d nodeforce = mesh->nodes[v].force.internal * (-1); + internalForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1], + nodeforce[2]}; + internalForcesNorm[mesh->getIndex(v)] = nodeforce.norm(); + // External force + const Vector6d nodeExternalForce = mesh->nodes[v].force.external; + externalForces[mesh->getIndex(v)] = { + nodeExternalForce[0], nodeExternalForce[1], nodeExternalForce[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Internal force", internalForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Internal force") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("External force", externalForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("External force") + ->setEnabled(true); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("Internal force norm", internalForcesNorm); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Internal force norm") + ->setEnabled(true); + // per node internal forces + std::vector> internalAxialForces(mesh->VN()); + for (const VertexType &v : mesh->vert) { + const Vector6d nodeforce = mesh->nodes[v].force.internalAxial * (-1); + internalAxialForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1], + nodeforce[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Internal Axial force", internalAxialForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Internal Axial force") + ->setEnabled(false); + // per node internal first bending force + std::vector> internalFirstBendingTranslationForces( + mesh->VN()); + std::vector> internalFirstBendingRotationForces( + mesh->VN()); + std::vector> internalSecondBendingTranslationForces( + mesh->VN()); + std::vector> internalSecondBendingRotationForces( + mesh->VN()); + std::vector nRs(mesh->VN()); + std::vector theta2(mesh->VN()); + std::vector theta3(mesh->VN()); + for (const VertexType &v : mesh->vert) { + const Node &node = mesh->nodes[v]; + const Vector6d nodeForceFirst = node.force.internalFirstBending * (-1); + internalFirstBendingTranslationForces[mesh->getIndex(v)] = { + nodeForceFirst[0], nodeForceFirst[1], nodeForceFirst[2]}; + internalFirstBendingRotationForces[mesh->getIndex(v)] = { + nodeForceFirst[3], nodeForceFirst[4], 0}; + + const Vector6d nodeForceSecond = node.force.internalSecondBending * (-1); + internalSecondBendingTranslationForces[mesh->getIndex(v)] = { + nodeForceSecond[0], nodeForceSecond[1], nodeForceSecond[2]}; + internalSecondBendingRotationForces[mesh->getIndex(v)] = { + nodeForceSecond[3], nodeForceSecond[4], 0}; + nRs[mesh->getIndex(v)] = node.nR; + const std::vector incidentEdges = node.incidentElements; + const EdgeIndex ei = mesh->getIndex(incidentEdges.back()); + // theta2[mesh->getIndex(v)] = + // node.rotationalDisplacements.at(ei).theta2; + // theta3[mesh->getIndex(v)] = + // node.rotationalDisplacements.at(ei).theta3; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("First bending force-Translation", + internalFirstBendingTranslationForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("First bending force-Translation") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("First bending force-Rotation", + internalFirstBendingRotationForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("First bending force-Rotation") + ->setEnabled(false); + + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Second bending force-Translation", + internalSecondBendingTranslationForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Second bending force-Translation") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Second bending force-Rotation", + internalSecondBendingRotationForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Second bending force-Rotation") + ->setEnabled(false); + + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("nR", nRs); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("nR") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("theta3", theta3); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("theta3") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("theta2", theta2); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("theta2") + ->setEnabled(false); + + // per node residual forces + std::vector> residualForces(mesh->VN()); + std::vector residualForcesNorm(mesh->VN()); + for (const VertexType &v : mesh->vert) { + const Vector6d nodeResidualForce = mesh->nodes[v].force.residual; + residualForces[mesh->getIndex(v)] = { + nodeResidualForce[0], nodeResidualForce[1], nodeResidualForce[2]}; + residualForcesNorm[mesh->getIndex(v)] = nodeResidualForce.norm(); + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeVectorQuantity("Residual force", residualForces); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Residual force") + ->setEnabled(false); + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("Residual force norm", residualForcesNorm); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("Residual force norm") + ->setEnabled(false); + + // per node x acceleration + std::vector accelerationX(mesh->VN()); + for (const VertexType &v : mesh->vert) { + accelerationX[mesh->getIndex(v)] = mesh->nodes[v].acceleration[0]; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addNodeScalarQuantity("Node acceleration x", accelerationX); + + // per element t1 + std::vector> t1s(mesh->EN()); + for (const EdgeType &e : mesh->edge) { + const VectorType &t1 = mesh->elements[e].frame.t1; + t1s[mesh->getIndex(e)] = {t1[0], t1[1], t1[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addEdgeVectorQuantity("t1Check", t1s); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("t1Check") + ->setEnabled(false); + // per element t2 + std::vector> t2s(mesh->EN()); + for (const EdgeType &e : mesh->edge) { + const VectorType &t2 = mesh->elements[e].frame.t2; + t2s[mesh->getIndex(e)] = {t2[0], t2[1], t2[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addEdgeVectorQuantity("t2", t2s); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("t2") + ->setEnabled(false); + // per element t3 + std::vector> t3s(mesh->EN()); + for (const EdgeType &e : mesh->edge) { + const VectorType &t3 = mesh->elements[e].frame.t3; + t3s[mesh->getIndex(e)] = {t3[0], t3[1], t3[2]}; + } + polyscope::getCurveNetwork("Deformed edge mesh") + ->addEdgeVectorQuantity("t3", t3s); + polyscope::getCurveNetwork("Deformed edge mesh") + ->getQuantity("t3") + ->setEnabled(false); + + // Specify the callback + polyscope::state::userCallback = [&]() { + // Since options::openImGuiWindowForUserCallback == true by default, + // we can immediately start using ImGui commands to build a UI + + ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, + // instead of full width. Must have + // matching PopItemWidth() below. + + ImGui::InputInt("Simulation drawing step", + &mDrawingStep); // set a int variable + + ImGui::PopItemWidth(); + }; + + if (!screenshotsFolder.empty()) { + static bool firstDraw = true; + if (firstDraw) { + for (const auto &entry : + std::filesystem::directory_iterator(screenshotsFolder)) + std::filesystem::remove_all(entry.path()); + // polyscope::view::processZoom(5); + firstDraw = false; + } + polyscope::screenshot( + std::filesystem::path(screenshotsFolder) + .append(std::to_string(currentSimulationStep) + ".png") + .string(), + false); + } else { + polyscope::show(); + } +} + +SimulationResults +FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, + const bool &shouldDraw, + const size_t &maxDRMIterations) { + assert(job.mesh.operator bool()); + auto t1 = std::chrono::high_resolution_clock::now(); + + constrainedVertices = job.fixedVertices; + if (!job.nodalForcedDisplacements.empty()) { + for (std::pair viDisplPair : + job.nodalForcedDisplacements) { + const VertexIndex vi = viDisplPair.first; + constrainedVertices[vi].insert({0, 1, 2}); + } + } + + // VCGEdgeMesh *jobEdgeMesh = job.mesh.get(); + mesh = std::make_unique(*job.mesh); + // ElementalMesh *jobElementalMesh = job.mesh.get(); + // vcg::tri::Append::MeshCopy( + // *(this->mesh), *jobElementalMesh); + if (beVerbose) { + std::cout << "Executing simulation for mesh with:" << mesh->VN() + << " nodes and " << mesh->EN() << " elements." << std::endl; + } + reset(); + computeRigidSupports(); + if (shouldDraw) { + if (!polyscope::state::initialized) { + initPolyscope(); + } + const std::string deformedMeshName = "Deformed edge mesh"; + if (!polyscope::hasCurveNetwork(deformedMeshName)) { + polyscope::registerCurveNetwork( + deformedMeshName, mesh->getEigenVertices(), mesh->getEigenEdges()); + } + + registerWorldAxes(); + } + // const double beamRadius = mesh.getBeamDimensions()[0].od / 2; + // std::cout << "BeamRadius:" << beamRadius << std::endl; + // polyscope::getCurveNetwork("Deformed edge mesh") + // // + // // ->setRadius(beamRadius); + // ->setRadius(0.01); + + // const bool hasExternalLoads = + // !job.nodalExternalForces.empty() || + // !job.nodalForcedDisplacements.empty(); + // assert(hasExternalLoads); + for (auto fixedVertex : job.fixedVertices) { + assert(fixedVertex.first < mesh->VN()); + } + + updateElementalFrames(); + updateNodalMasses(); + if (job.nodalExternalForces.empty()) { + shouldApplyInitialDistortion = true; + } + // std::unordered_map temporaryForces{ + // // // {2, Eigen::Vector3d(0, 0, 1000)}, + // // // {12, Eigen::Vector3d(0, 0, 1000)}, + // // // {18, Eigen::Vector3d(0, 0, 1000)}, + // // // {8, Eigen::Vector3d(0, 0, 1000)}}; + // {10, Eigen::Vector3d(0, 0, 10000)}}; + // std::unordered_map temporaryForces; + // for (VertexIndex vi = 0; vi < mesh.VN(); vi++) { + // for (VertexType &v : mesh.vert) { + // const VertexIndex vi = mesh.getIndex(v); + // VectorType allowedDoFType(1, 1, 1); + // if (constrainedVertices.contains(vi)) { + // std::unordered_set constrainedDof = + // constrainedVertices.at(vi); + // if (constrainedDof.contains(0)) { + // allowedDoFType[0] = 0; + // } else if (constrainedDof.contains(1)) { + // allowedDoFType[1] = 0; + // } else if (constrainedDof.contains(2)) { + // allowedDoFType[2] = 0; + // } + // } + // VectorType desiredDisplacement = VectorType(0, 0, 0.2); + // VectorType displacement(allowedDoFType[0] * desiredDisplacement[0], + // allowedDoFType[1] * desiredDisplacement[1], + // allowedDoFType[2] * desiredDisplacement[2]); + // v.P() = v.P() + displacement; + // temporaryForces.insert({vi, Eigen::Vector3d(0, 0, 100000)}); + // } + // updateNodalExternalForces(temporaryForces, job.fixedVertices); + // appliedTemporaryForce = true; + // } else { + // std::unordered_map appliedLoad = + // job.nodalExternalForces; + // const size_t numberOfStepsForApplyingExternalLoads = 10; + // int externalLoadStep = 1; + // std::for_each(appliedLoad.begin(), appliedLoad.end(), [](auto &pair) { + // pair.second /= numberOfStepsForApplyingExternalLoads; + // }); + // updateNodalExternalForces(appliedLoad, constrainedVertices); + updateNodalExternalForces(job.nodalExternalForces, constrainedVertices); + // } + + while (currentSimulationStep < maxDRMIterations) { + // while (true) { + updateNormalDerivatives(); + updateT1Derivatives(); + updateRDerivatives(); + updateT2Derivatives(); + updateT3Derivatives(); + // updateRotationalDisplacements(); + // if (currentSimulationStep > lastPulseStepIndex) { + // const std::vector internalForces = + updateResidualForcesOnTheFly(constrainedVertices); + //#pragma omp parallel for schedule(static) num_threads(8) + // double totalResidualForcesNorm = 0; + // for (size_t vi = 0; vi < mesh.VN(); vi++) { + // Node::Forces &force = mesh.nodes[vi].force; + // // const Vector6d ¶llelForce = internalForcesParallel[vi]; + // // const Vector6d &serialForce = internalForces[vi]; + // // const double normOfDifference = (serialForce + + // (parallelForce + // * + // // -1)).norm(); assert(normOfDifference < 0.000001); + // force.residual = force.external - internalForces[vi]; + // const double residualForceNorm = (force.residual).norm(); + // totalResidualForcesNorm += residualForceNorm; + // // assert(residualForceNorm < std::pow(10, 20)); + // } + // mesh.totalResidualForcesNorm = totalResidualForcesNorm; + // } else { + // if (currentSimulationStep == lastPulseStepIndex && + // appliedTemporaryForce) { + // for (const VertexType &v : mesh.vert) { + // Node &node = mesh.nodes[v]; + // node.force.external = std::optional(); + // } + // } + // updateNodalInternalForces(job.fixedVertices); + // } + + // TODO: write parallel function for updating positions + // TODO: Make a single function out of updateResidualForcesOnTheFly + // updatePositionsOnTheFly + // updatePositionsOnTheFly(constrainedVertices); + updateNodalMasses(); + updateNodalAccelerations(); + updateNodalVelocities(); + updateNodalDisplacements(); + applyDisplacements(constrainedVertices); + if (!job.nodalForcedDisplacements.empty()) { + applyForcedDisplacements( + + job.nodalForcedDisplacements); + } + updateElementalLengths(); + + // if (!std::isfinite(mesh.currentTotalKineticEnergy)) { + // std::cerr << "Infinite kinetic energy. Breaking simulation.." + // << std::endl; + // break; + // } + // assert(std::isfinite(mesh.currentTotalKineticEnergy)); + // mesh.printVertexCoordinates(mesh.VN() / 2); + // draw(); + if (shouldDraw && + currentSimulationStep % mDrawingStep == 0 /* && +currentSimulationStep > maxDRMIterations*/) { + // std::string saveTo = std::filesystem::current_path() + // .append("Debugging_files") + // .append("Screenshots") + // .string(); + // draw(saveTo); + // SimulationResultsReporter::createPlot( + // "Number of iterations", "Log of Kinetic Energy", + // history.kineticEnergy, + // std::filesystem::path(saveTo).append( + // std::to_string(currentSimulationStep) + "_kin.png")); + + // draw(); + // auto t2 = std::chrono::high_resolution_clock::now(); + // auto timePerNodePerIteration = + // std::chrono::duration_cast(t2 - + // t1) + // .count() / + // (mesh.VN() * (currentSimulationStep + 1)); + // std::cout << "Time/(node*iteration) " + // << timePerNodePerIteration * std::pow(10, -6) << "s" + // << std::endl; + draw(); + } + if (currentSimulationStep != 0) { + history.stepPulse(*mesh); + } + // t = t + Dt; + currentSimulationStep++; + // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy + // << std::endl; + // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm + // << std::endl; + if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) { + if (/*mesh.totalPotentialEnergykN < 10 ||*/ + mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + if (beVerbose) { + std::cout << "Convergence after " << currentSimulationStep << " steps" + << std::endl; + std::cout << "Residual force of magnitude:" + << mesh->previousTotalResidualForcesNorm << std::endl; + std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy + << std::endl; + mesh->totalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + std::cout << "Total potential:" << mesh->totalPotentialEnergykN + << " kNm" << std::endl; + } + // if (externalLoadStep != + // numberOfStepsForApplyingExternalLoads) + // { + // std::for_each( + // appliedLoad.begin(), appliedLoad.end(), [&](auto + // &pair) + // { + // pair.second += + // job.nodalExternalForces.at(pair.first) + // / + // numberOfStepsForApplyingExternalLoads; + // }); + // updateNodalExternalForces(appliedLoad, + // constrainedVertices); + // } else { + break; + // } + } + // history.markRed(currentSimulationStep); + // std::cout << "Total potential:" << mesh.totalPotentialEnergykN + // << " kNm" + // << std::endl; + // reset displacements to previous position where the peak occured + for (VertexType &v : mesh->vert) { + Node &node = mesh->nodes[v]; + node.displacements = node.displacements - node.velocity * Dt; + } + applyDisplacements(constrainedVertices); + if (!job.nodalForcedDisplacements.empty()) { + applyForcedDisplacements( + + job.nodalForcedDisplacements); + } + updateElementalLengths(); + resetVelocities(); + Dt = Dt * xi; + // std::cout << "Residual forces norm:" << + // mesh.totalResidualForcesNorm + // << std::endl; + } + } + if (currentSimulationStep == maxDRMIterations) { + std::cout << "Did not reach equilibrium before reaching the maximum number " + "of DRM steps (" + << maxDRMIterations << "). Breaking simulation" << std::endl; + } else if (beVerbose) { + auto t2 = std::chrono::high_resolution_clock::now(); + double simulationDuration = + std::chrono::duration_cast(t2 - t1).count(); + simulationDuration /= 1000; + std::cout << "Simulation converged after " << simulationDuration << "s" + << std::endl; + std::cout << "Time/(node*iteration) " + << simulationDuration / + (static_cast(currentSimulationStep) * mesh->VN()) + << "s" << std::endl; + } + // mesh.printVertexCoordinates(mesh.VN() / 2); + if (shouldDraw) { + draw(); + } + // if (!polyscope::hasCurveNetwork(deformedMeshName)) { + // polyscope::removeCurveNetwork(deformedMeshName); + // } + // debugger.createPlots(); + SimulationResults results = computeResults(*job.mesh); + // Eigen::write_binary("optimizedResult.eigenBin", + // results.nodalDisplacements); + + // const std::string groundOfTruthBinaryFilename = + // "../../groundOfTruths/grid6_groundOfTruth.eigenBin"; + // // "../../groundOfTruths/longSpanGridshell_groundOfTruth.eigenBin"; + // assert(std::filesystem::exists( + // std::filesystem::path(groundOfTruthBinaryFilename))); + // Eigen::MatrixX3d groundOfTruthMatrix; + // Eigen::read_binary(groundOfTruthBinaryFilename, groundOfTruthMatrix); + // assert(results.nodalDisplacements.isApprox(groundOfTruthMatrix)); + return results; + // return history; +} diff --git a/beamformfinder.hpp b/beamformfinder.hpp new file mode 100644 index 0000000..bd268da --- /dev/null +++ b/beamformfinder.hpp @@ -0,0 +1,224 @@ +#ifndef BEAMFORMFINDER_HPP +#define BEAMFORMFINDER_HPP + +#include "elementalmesh.hpp" +#include "polyscope/curve_network.h" +#include "polyscope/polyscope.h" +#include "simulationresult.hpp" +#include +#include +#include + +struct SimulationJob; + +enum DoF { Ux = 0, Uy, Uz, Nx, Ny, Nr, NumDoF }; +using DoFType = int; +using EdgeType = VCGEdgeMesh::EdgeType; +using VertexType = VCGEdgeMesh::VertexType; + +struct DifferentiateWithRespectTo { + const VertexType &v; + const DoFType &dofi; +}; +class FormFinder { + +private: + const double Dtini{0.1}; + double Dt{Dtini}; + double t{0}; + const double xi{0.9969}; + const double totalResidualForcesNormThreshold{300}; + size_t currentSimulationStep{0}; + int mDrawingStep{1}; + std::unique_ptr mesh; + std::unordered_map> + constrainedVertices; + SimulationHistory history; + // Eigen::Tensor theta3Derivatives; + // std::unordered_map theta3Derivatives; + bool shouldApplyInitialDistortion = false; + std::unordered_set rigidSupports; + + void reset(); + void updateNodalInternalForces( + const std::unordered_map> + &fixedVertices); + void updateNodalExternalForces( + const std::unordered_map &nodalForces, + const std::unordered_map> + &fixedVertices); + void updateResidualForces(); + void updateRotationalDisplacements(); + void updateElementalLengths(); + + void updateNodalMasses(); + + void updateNodalAccelerations(); + + void updateNodalVelocities(); + + void updateNodalDisplacements(); + + void applyForcedDisplacements( + const std::unordered_map + nodalForcedDisplacements); + + ::Vector6d computeElementTorsionalForce( + const Element &element, const Vector6d &displacementDifference, + const std::unordered_set &constrainedDof); + + // BeamFormFinder::Vector6d computeElementFirstBendingForce( + // const Element &element, const Vector6d &displacementDifference, + // const std::unordered_set &constrainedDof); + + // BeamFormFinder::Vector6d computeElementSecondBendingForce( + // const Element &element, const Vector6d &displacementDifference, + // const std::unordered_set &constrainedDof); + + void updateKineticEnergy(); + + void resetVelocities(); + + SimulationResults computeResults(const SimulationMesh &initialiMesh); + + void updateNodePosition( + VertexType &v, + const std::unordered_map> + &fixedVertices); + + void applyDisplacements( + const std::unordered_map> + &fixedVertices); + + void draw(const string &screenshotsFolder); + + void + updateNodalInternalForce(Vector6d &nodalInternalForce, + const Vector6d &elementInternalForce, + const std::unordered_set &nodalFixedDof); + + ::Vector6d computeElementInternalForce( + const Element &elem, const Node &n0, const Node &n1, + const std::unordered_set &n0ConstrainedDof, + const std::unordered_set &n1ConstrainedDof); + + ::Vector6d computeElementAxialForce(const ::EdgeType &e) const; + VectorType computeDisplacementDifferenceDerivative( + const EdgeType &e, const DifferentiateWithRespectTo &dui) const; + double + computeDerivativeElementLength(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; + + VectorType computeDerivativeT1(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; + + VectorType + computeDerivativeOfNormal(const VertexType &v, + const DifferentiateWithRespectTo &dui) const; + + VectorType computeDerivativeT3(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; + + VectorType computeDerivativeT2(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; + + double computeDerivativeTheta2(const EdgeType &e, const VertexIndex &evi, + const VertexIndex &dwrt_evi, + const DoFType &dwrt_dofi) const; + + void updateElementalFrames(); + + VectorType computeDerivativeOfR(const EdgeType &e, + const DifferentiateWithRespectTo &dui) const; + + bool isRigidSupport(const VertexType &v) const; + + static double computeDerivativeOfNorm(const VectorType &x, + const VectorType &derivativeOfX); + static VectorType computeDerivativeOfCrossProduct( + const VectorType &a, const VectorType &derivativeOfA, const VectorType &b, + const VectorType &derivativeOfB); + + double computeTheta3(const EdgeType &e, const VertexType &v); + double computeDerivativeTheta3(const EdgeType &e, const VertexType &v, + const DifferentiateWithRespectTo &dui) const; + double computeTotalPotentialEnergy() const; + void computeRigidSupports(); + void updateNormalDerivatives(); + void updateT1Derivatives(); + void updateT2Derivatives(); + void updateT3Derivatives(); + void updateRDerivatives(); + + double computeDerivativeTheta1(const EdgeType &e, const VertexIndex &evi, + const VertexIndex &dwrt_evi, + const DoFType &dwrt_dofi) const; + + // void updatePositionsOnTheFly( + // const std::unordered_map> + // &fixedVertices); + + void updateResidualForcesOnTheFly( + const std::unordered_map> + &fixedVertices); + + void updatePositionsOnTheFly( + const std::unordered_map> + &fixedVertices); + +public: + FormFinder(); + SimulationResults executeSimulation( + const SimulationJob &job, const bool &beVerbose = false, + const bool &shouldDraw = false, + const size_t &maxDRMIterations = FormFinder::maxDRMIterations); + inline static const size_t maxDRMIterations{50000}; +}; + +template PointType Cross(PointType p1, PointType p2) { + return p1 ^ p2; +} + +inline size_t currentStep{0}; +inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei, + const DoFType &dofi) { + const size_t numberOfVertices = 10; + const VertexIndex middleNodeIndex = numberOfVertices / 2; + // return vi == middleNodeIndex && dofi == 1; + return dofi == 1 && ((vi == 1 && ei == 0) || (vi == 9 && ei == 9)); +} + +inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei) { + const size_t numberOfVertices = 10; + const VertexIndex middleNodeIndex = numberOfVertices / 2; + return (vi == middleNodeIndex); + // return (vi == 0 || vi == numberOfVertices - 1) && currentStep == 1; + return (vi == 1 && ei == 0) || (vi == 9 && ei == 9); +} + +namespace Eigen { +template +void write_binary(const std::string &filename, const Matrix &matrix) { + std::ofstream out(filename, + std::ios::out | std::ios::binary | std::ios::trunc); + typename Matrix::Index rows = matrix.rows(), cols = matrix.cols(); + out.write((char *)(&rows), sizeof(typename Matrix::Index)); + out.write((char *)(&cols), sizeof(typename Matrix::Index)); + out.write((char *)matrix.data(), + rows * cols * sizeof(typename Matrix::Scalar)); + out.close(); +} +template +void read_binary(const std::string &filename, Matrix &matrix) { + std::ifstream in(filename, std::ios::in | std::ios::binary); + typename Matrix::Index rows = 0, cols = 0; + in.read((char *)(&rows), sizeof(typename Matrix::Index)); + in.read((char *)(&cols), sizeof(typename Matrix::Index)); + matrix.resize(rows, cols); + in.read((char *)matrix.data(), rows * cols * sizeof(typename Matrix::Scalar)); + in.close(); +} +} // namespace Eigen + +#endif // BEAMFORMFINDER_HPP diff --git a/edgemesh.cpp b/edgemesh.cpp new file mode 100644 index 0000000..9c3d5b6 --- /dev/null +++ b/edgemesh.cpp @@ -0,0 +1,382 @@ +#include "edgemesh.hpp" +#include "vcg/simplex/face/topology.h" + +Eigen::MatrixX2i VCGEdgeMesh::getEigenEdges() const { return eigenEdges; } + +Eigen::MatrixX3d VCGEdgeMesh::getEigenVertices() const { + // getVertices(eigenVertices); + return eigenVertices; +} + +Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const { + return eigenEdgeNormals; +} + +std::vector +VCGEdgeMesh::getBeamDimensions() const { + return handleBeamDimensions._handle->data; +} + +std::vector VCGEdgeMesh::getBeamMaterial() const { + return handleBeamMaterial._handle->data; +} + +bool VCGEdgeMesh::savePly(const std::string plyFilename) { + nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; + customAttrib.GetMeshAttrib(plyFilename); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, + &handleBeamDimensions[0]); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, + &handleBeamMaterial[0]); + // Load the ply file + 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( + plyFilename.c_str(), *this, mask, customAttrib, false) != 0) { + return false; + } + return true; +} + +void VCGEdgeMesh::GeneratedRegularSquaredPattern( + const double angleDeg, std::vector> &pattern, + const size_t &desiredNumberOfSamples) { + static const size_t piSamples = 10; + + // generate a pattern in a 1x1 quad + const vcg::Point2d offset(0, 0); + + const size_t samplesNo = desiredNumberOfSamples; + // std::max(desiredNumberOfSamples, size_t(piSamples * (angleDeg / + // 180))); + const double angle = vcg::math::ToRad(angleDeg); + + pattern.clear(); + + // first arm + std::vector arm; + { + for (int k = 0; k <= samplesNo; k++) { + const double t = double(k) / samplesNo; + const double a = (1 - t) * angle; + // const double r = vcg::math::Sin(t*M_PI_2) /*(1-((1-t)*(1-t)))*/ * 0.5; + const double r = t * 0.5; // linear + + vcg::Point2d p(vcg::math::Cos(a), vcg::math::Sin(a)); + + arm.push_back((p * r)); + } + pattern.push_back(arm); + } + + // other arms + for (int i = 0; i < 3; i++) { + for (vcg::Point2d &p : arm) { + p = vcg::Point2d(-p.Y(), p.X()); + } + pattern.push_back(arm); + } + + assert(pattern.size() == 4); + + // offset all + for (auto &arm : pattern) { + for (vcg::Point2d &p : arm) { + p += offset; + } + } +} + +void VCGEdgeMesh::createSpiral(const float °reesOfArm, + const size_t &numberOfSamples) { + std::vector> spiralPoints; + GeneratedRegularSquaredPattern(degreesOfArm, spiralPoints, numberOfSamples); + + for (size_t armIndex = 0; armIndex < spiralPoints.size(); armIndex++) { + for (size_t pointIndex = 0; pointIndex < spiralPoints[armIndex].size() - 1; + pointIndex++) { + const vcg::Point2d p0 = spiralPoints[armIndex][pointIndex]; + const vcg::Point2d p1 = spiralPoints[armIndex][pointIndex + 1]; + CoordType n(0, 0, 1); + auto ei = vcg::tri::Allocator::AddEdge( + *this, VCGEdgeMesh::CoordType(p0.X(), p0.Y(), 0), + VCGEdgeMesh::CoordType(p1.X(), p1.Y(), 0)); + ei->cV(0)->N() = n; + ei->cV(1)->N() = n; + } + } + + // setDefaultAttributes(); +} + +bool VCGEdgeMesh::createGrid(const size_t squareGridDimension) { + return createGrid(squareGridDimension, squareGridDimension); +} + +bool VCGEdgeMesh::createGrid(const size_t desiredWidth, + const size_t desiredHeight) { + std::cout << "Grid of dimensions:" << desiredWidth << "," << desiredHeight + << std::endl; + const VCGEdgeMesh::CoordType n(0, 0, 1); + int x = 0; + int y = 0; + // for (size_t vi = 0; vi < numberOfVertices; vi++) { + while (y <= desiredHeight) { + // std::cout << x << " " << y << std::endl; + auto p = VCGEdgeMesh::CoordType(x, y, 0); + vcg::tri::Allocator::AddVertex(*this, p, n); + x++; + if (x > desiredWidth) { + x = 0; + y++; + } + } + + for (size_t vi = 0; vi < VN(); vi++) { + int x = vi % (desiredWidth + 1); + int y = vi / (desiredWidth + 1); + const bool isCornerNode = (y == 0 && x == 0) || + (y == 0 && x == desiredWidth) || + (y == desiredHeight && x == 0) || + (y == desiredHeight && x == desiredWidth); + if (isCornerNode) { + continue; + } + if (y == 0) { // row 0.Connect with node above + vcg::tri::Allocator::AddEdge(*this, vi, + vi + desiredWidth + 1); + continue; + } else if (x == 0) { // col 0.Connect with node to the right + vcg::tri::Allocator::AddEdge(*this, vi, vi + 1); + continue; + } else if (y == desiredHeight) { // row 0.Connect with node below + // vcg::tri::Allocator::AddEdge(*this, vi, + // vi - (desiredWidth + + // 1)); + continue; + } else if (x == desiredWidth) { // row 0.Connect with node to the left + // vcg::tri::Allocator::AddEdge(*this, vi, vi - 1); + continue; + } + + vcg::tri::Allocator::AddEdge(*this, vi, vi + desiredWidth + 1); + vcg::tri::Allocator::AddEdge(*this, vi, vi + 1); + // vcg::tri::Allocator::AddEdge(*this, vi, + // vi - (desiredWidth + 1)); + // vcg::tri::Allocator::AddEdge(*this, vi, vi - 1); + } + + vcg::tri::Allocator::DeleteVertex(*this, vert[0]); + vcg::tri::Allocator::DeleteVertex(*this, vert[desiredWidth]); + vcg::tri::Allocator::DeleteVertex( + *this, vert[desiredHeight * (desiredWidth + 1)]); + vcg::tri::Allocator::DeleteVertex( + *this, vert[(desiredHeight + 1) * (desiredWidth + 1) - 1]); + vcg::tri::Allocator::CompactVertexVector(*this); + getEdges(eigenEdges); + getVertices(eigenVertices); + // vcg::tri::Allocator::CompactEdgeVector(*this); + + // const size_t numberOfEdges = + // desiredHeight * (desiredWidth - 1) + desiredWidth * (desiredHeight - + // 1); + // handleBeamDimensions._handle->data.resize( + // numberOfEdges, CylindricalElementDimensions(0.03, 0.026)); + // handleBeamMaterial._handle->data.resize(numberOfEdges, + // ElementMaterial(0.3, 200)); + + return true; +} + +bool VCGEdgeMesh::loadFromPly(const std::string plyFilename) { + std::string usedPath = plyFilename; + if (std::filesystem::path(plyFilename).is_relative()) { + usedPath = std::filesystem::absolute(plyFilename).string(); + } + assert(std::filesystem::exists(usedPath)); + this->Clear(); + const bool useDefaultImporter = false; + if (useDefaultImporter) { + if (!loadUsingDefaultLoader(usedPath)) { + return false; + } + + eigenEdgeNormals.resize(EN(), 3); + for (int i = 0; i < EN(); i++) { + eigenEdgeNormals.row(i) = Eigen::Vector3d(0, 1, 0); + } + } else { + if (!loadUsingNanoply(usedPath)) { + std::cerr << "Error: Unable to open " + usedPath << std::endl; + return false; + } + } + getEdges(eigenEdges); + getVertices(eigenVertices); + vcg::tri::UpdateTopology::VertexEdge(*this); + + std::cout << plyFilename << " was loaded successfuly." << std::endl; + std::cout << "Mesh has " << EN() << " edges." << std::endl; + return true; +} + +bool VCGEdgeMesh::loadUsingNanoply(const std::string &plyFilename) { + + this->Clear(); + // assert(plyFileHasAllRequiredFields(plyFilename)); + nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; + customAttrib.GetMeshAttrib(plyFilename); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); + // Load the ply file + unsigned int mask = 0; + mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; + mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; + mask |= nanoply::NanoPlyWrapper::IO_EDGEATTRIB; + if (nanoply::NanoPlyWrapper::LoadModel( + plyFilename.c_str(), *this, mask, customAttrib) != 0) { + return false; + } + return true; +} + +bool VCGEdgeMesh::plyFileHasAllRequiredFields(const std::string &plyFilename) { + const nanoply::Info info(plyFilename); + const std::vector::const_iterator edgeElemIt = + std::find_if(info.elemVec.begin(), info.elemVec.end(), + [&](const nanoply::PlyElement &plyElem) { + return plyElem.plyElem == nanoply::NNP_EDGE_ELEM; + }); + if (edgeElemIt == info.elemVec.end()) { + std::cerr << "Ply file is missing edge elements." << std::endl; + return false; + } + + const std::vector &edgePropertyVector = + edgeElemIt->propVec; + return hasPlyEdgeProperty(plyFilename, edgePropertyVector, + plyPropertyBeamDimensionsID) && + hasPlyEdgeProperty(plyFilename, edgePropertyVector, + plyPropertyBeamMaterialID); +} + +bool VCGEdgeMesh::hasPlyEdgeProperty( + const std::string &plyFilename, + const std::vector &edgeProperties, + const std::string &plyEdgePropertyName) { + const bool hasEdgeProperty = hasProperty(edgeProperties, plyEdgePropertyName); + if (!hasEdgeProperty) { + std::cerr << "Ply file " + plyFilename + + " is missing the propertry:" + plyEdgePropertyName + << std::endl; + return false; + } + return true; +} + +bool VCGEdgeMesh::hasProperty(const std::vector &v, + const std::string &propertyName) { + return v.end() != std::find_if(v.begin(), v.end(), + [&](const nanoply::PlyProperty &plyProperty) { + return plyProperty.name == propertyName; + }); +} + +Eigen::MatrixX3d VCGEdgeMesh::getNormals() const { + Eigen::MatrixX3d vertexNormals; + vertexNormals.resize(VN(), 3); + + for (int vertexIndex = 0; vertexIndex < VN(); vertexIndex++) { + VCGEdgeMesh::CoordType vertexNormal = + vert[static_cast(vertexIndex)].cN(); + vertexNormals.row(vertexIndex) = + vertexNormal.ToEigenVector(); + } + + return vertexNormals; +} +void VCGEdgeMesh::getEdges(Eigen::MatrixX3d &edgeStartingPoints, + Eigen::MatrixX3d &edgeEndingPoints) const { + edgeStartingPoints.resize(EN(), 3); + edgeEndingPoints.resize(EN(), 3); + for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) { + const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex]; + edgeStartingPoints.row(edgeIndex) = + edge.cP(0).ToEigenVector(); + edgeEndingPoints.row(edgeIndex) = + edge.cP(1).ToEigenVector(); + } +} + +VCGEdgeMesh::VCGEdgeMesh() { + handleBeamDimensions = vcg::tri::Allocator::AddPerEdgeAttribute< + CylindricalElementDimensions>(*this, plyPropertyBeamDimensionsID); + handleBeamMaterial = + vcg::tri::Allocator::AddPerEdgeAttribute( + *this, plyPropertyBeamMaterialID); +} + +void VCGEdgeMesh::updateEigenEdgeAndVertices() { + getEdges(eigenEdges); + getVertices(eigenVertices); +} + +void VCGEdgeMesh::copy(VCGEdgeMesh &mesh) { + vcg::tri::Append::MeshCopy(*this, mesh); + label = mesh.getLabel(); + eigenEdges = mesh.getEigenEdges(); + if (eigenEdges.rows() == 0) { + getEdges(eigenEdges); + } + eigenVertices = mesh.getEigenVertices(); + if (eigenVertices.rows() == 0) { + getVertices(eigenVertices); + } + vcg::tri::UpdateTopology::VertexEdge(*this); +} + +void VCGEdgeMesh::getVertices(Eigen::MatrixX3d &vertices) { + vertices = Eigen::MatrixX3d(); + vertices.resize(VN(), 3); + for (int vi = 0; vi < VN(); vi++) { + if (vert[vi].IsD()) { + continue; + } + VCGEdgeMesh::CoordType vertexCoordinates = + vert[static_cast(vi)].cP(); + vertices.row(vi) = vertexCoordinates.ToEigenVector(); + } +} + +std::string VCGEdgeMesh::getLabel() const { return label; } + +void VCGEdgeMesh::setLabel(const std::string &value) { label = value; } + +void VCGEdgeMesh::getEdges(Eigen::MatrixX2i &edges) { + edges = Eigen::MatrixX2i(); + edges.resize(EN(), 2); + for (int edgeIndex = 0; edgeIndex < EN(); edgeIndex++) { + const VCGEdgeMesh::EdgeType &edge = this->edge[edgeIndex]; + const size_t nodeIndex0 = vcg::tri::Index(*this, edge.cV(0)); + const size_t nodeIndex1 = vcg::tri::Index(*this, edge.cV(1)); + edges.row(edgeIndex) = Eigen::Vector2i(nodeIndex0, nodeIndex1); + } +} + +void VCGEdgeMesh::printVertexCoordinates(const size_t &vi) const { + std::cout << "vi:" << vi << " " << vert[vi].cP()[0] << " " << vert[vi].cP()[1] + << " " << vert[vi].cP()[2] << std::endl; +} + +void VCGEdgeMesh::registerForDrawing() const { + initPolyscope(); + polyscope::registerCurveNetwork(label, getEigenVertices(), getEigenEdges()); +} diff --git a/edgemesh.hpp b/edgemesh.hpp new file mode 100644 index 0000000..af2fadd --- /dev/null +++ b/edgemesh.hpp @@ -0,0 +1,105 @@ +#ifndef EDGEMESH_HPP +#define EDGEMESH_HPP +#include "beam.hpp" +#include "polymesh.hpp" +#include "utilities.hpp" +#include "vcgtrimesh.hpp" +#include +#include +#include +#include + +using EdgeIndex = size_t; + +class VCGEdgeMeshEdgeType; +class VCGEdgeMeshVertexType; + +struct VCGEdgeMeshUsedTypes + : public vcg::UsedTypes::AsVertexType, + vcg::Use::AsEdgeType> {}; + +class VCGEdgeMeshVertexType + : public vcg::Vertex {}; +class VCGEdgeMeshEdgeType + : public vcg::Edge {}; + +class VCGEdgeMesh : public vcg::tri::TriMesh, + std::vector> { + const std::string plyPropertyBeamDimensionsID{"beam_dimensions"}; + const std::string plyPropertyBeamMaterialID{"beam_material"}; + VCGEdgeMesh::PerEdgeAttributeHandle + handleBeamDimensions; + VCGEdgeMesh::PerEdgeAttributeHandle handleBeamMaterial; + +protected: + Eigen::MatrixX2i eigenEdges; + Eigen::MatrixX3d eigenVertices; + Eigen::MatrixX3d eigenEdgeNormals; + std::string label{"No_name"}; + + void getEdges(Eigen::MatrixX2i &edges); + void getVertices(Eigen::MatrixX3d &vertices); + +public: + VCGEdgeMesh(); + template + size_t getIndex(const MeshElement &meshElement) { + return vcg::tri::Index(*this, meshElement); + } + void updateEigenEdgeAndVertices(); + void copy(VCGEdgeMesh &mesh); + + void getEdges(Eigen::MatrixX3d &edgeStartingPoints, + Eigen::MatrixX3d &edgeEndingPoints) const; + + Eigen::MatrixX3d getNormals() const; + + bool loadUsingDefaultLoader(const std::string &plyFilename); + bool hasProperty(const std::vector &v, + const std::string &propertyName); + + bool + hasPlyEdgeProperty(const std::string &plyFilename, + const std::vector &edgeProperties, + const std::string &plyEdgePropertyName); + + bool plyFileHasAllRequiredFields(const std::string &plyFilename); + + bool loadUsingNanoply(const std::string &plyFilename); + + bool loadFromPly(const std::string plyFilename); + + bool savePly(const std::string plyFilename); + + bool createGrid(const size_t squareGridDimension); + bool createGrid(const size_t desiredWidth, const size_t desiredHeight); + void createSpiral(const float °reesOfArm, const size_t &numberOfSamples); + + Eigen::MatrixX2i getEigenEdges() const; + Eigen::MatrixX3d getEigenVertices() const; + Eigen::MatrixX3d getEigenEdgeNormals() const; + std::vector getBeamDimensions() const; + std::vector getBeamMaterial() const; + void printVertexCoordinates(const size_t &vi) const; + void registerForDrawing() const; + + std::string getLabel() const; + void setLabel(const std::string &value); + +private: + void GeneratedRegularSquaredPattern( + const double angleDeg, std::vector> &pattern, + const size_t &desiredNumberOfSamples); +}; + +using VectorType = VCGEdgeMesh::CoordType; +using CoordType = VCGEdgeMesh::CoordType; +using VertexPointer = VCGEdgeMesh::VertexPointer; +using EdgePointer = VCGEdgeMesh::EdgePointer; +using ConstVCGEdgeMesh = VCGEdgeMesh; + +#endif // EDGEMESH_HPP diff --git a/elementalmesh.cpp b/elementalmesh.cpp new file mode 100644 index 0000000..cb91033 --- /dev/null +++ b/elementalmesh.cpp @@ -0,0 +1,209 @@ +#include "elementalmesh.hpp" + +SimulationMesh::SimulationMesh(FlatPattern &pattern) { + vcg::tri::MeshAssert::VertexNormalNormalized(pattern); + + vcg::tri::Append::MeshCopy(*this, pattern); + 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(); + label = pattern.getLabel(); + eigenEdges = pattern.getEigenEdges(); + eigenVertices = pattern.getEigenVertices(); +} + +SimulationMesh::SimulationMesh(SimulationMesh &elementalMesh) { + vcg::tri::Append::MeshCopy(*this, + elementalMesh); + elements = vcg::tri::Allocator::GetPerEdgeAttribute( + *this, std::string("Elements")); + nodes = vcg::tri::Allocator::GetPerVertexAttribute( + *this, std::string("Nodes")); + vcg::tri::UpdateTopology::VertexEdge(*this); + initializeNodes(); + + for (size_t ei = 0; ei < EN(); ei++) { + elements[ei] = elementalMesh.elements[ei]; + } + + label = label; + eigenEdges = elementalMesh.getEigenEdges(); + eigenVertices = elementalMesh.getEigenVertices(); +} + +void SimulationMesh::computeElementalProperties() { + const std::vector elementalDimensions = + getBeamDimensions(); + const std::vector elementalMaterials = getBeamMaterial(); + assert(EN() == elementalDimensions.size() && + elementalDimensions.size() == elementalMaterials.size()); + + for (const EdgeType &e : edge) { + const EdgeIndex ei = getIndex(e); + elements[e].properties = + Element::Properties{elementalDimensions[ei], elementalMaterials[ei]}; + } +} + +void SimulationMesh::initializeNodes() { + // set initial and previous locations, + for (const VertexType &v : vert) { + const VertexIndex vi = getIndex(v); + Node &node = nodes[v]; + node.vi = vi; + node.initialLocation = v.cP(); + node.previousLocation = v.cP(); + node.initialNormal = v.cN(); + node.derivativeOfNormal.resize(6, VectorType(0, 0, 0)); + + node.displacements[3] = + v.cN()[0]; // initialize nx diplacement with vertex normal x + // component + node.displacements[4] = + v.cN()[1]; // initialize ny(in the paper) diplacement with vertex + // normal + // y component. + + // Initialize incident elements + std::vector incidentElements; + vcg::edge::VEStarVE(&v, incidentElements); + assert( + vcg::tri::IsValidPointer(*this, incidentElements[0]) && + incidentElements.size() > 0); + nodes[v].incidentElements = std::move(incidentElements); + node.referenceElement = getReferenceElement(v); + // Initialze alpha angles + + const EdgeType &referenceElement = *getReferenceElement(v); + const VectorType t01 = + computeT1Vector(referenceElement.cP(0), referenceElement.cP(1)); + const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize(); + + for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) { + assert(referenceElement.cV(0) == ep->cV(0) || + referenceElement.cV(0) == ep->cV(1) || + referenceElement.cV(1) == ep->cV(0) || + referenceElement.cV(1) == ep->cV(1)); + const VectorType t1 = computeT1Vector(*ep); + const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize(); + const EdgeIndex ei = getIndex(ep); + const double alphaAngle = computeAngle(f01, f1, v.cN()); + node.alphaAngles[ei] = alphaAngle; + } + } +} + +void SimulationMesh::initializeElements() { + computeElementalProperties(); + for (const EdgeType &e : edge) { + Element &element = elements[e]; + element.ei = getIndex(e); + // Initialize lengths + const VCGEdgeMesh::CoordType p0 = e.cP(0); + const VCGEdgeMesh::CoordType p1 = e.cP(1); + const vcg::Segment3 s(p0, p1); + element.initialLength = s.Length(); + element.length = element.initialLength; + // Initialize const factors + element.axialConstFactor = + element.properties.E * element.properties.A / element.initialLength; + element.torsionConstFactor = + element.properties.G * element.properties.J / element.initialLength; + element.firstBendingConstFactor = 2 * element.properties.E * + element.properties.I2 / + element.initialLength; + element.secondBendingConstFactor = 2 * element.properties.E * + element.properties.I3 / + element.initialLength; + element.derivativeT1.resize( + 2, std::vector(6, VectorType(0, 0, 0))); + element.derivativeT2.resize( + 2, std::vector(6, VectorType(0, 0, 0))); + element.derivativeT3.resize( + 2, std::vector(6, VectorType(0, 0, 0))); + element.derivativeT1_jplus1.resize(6); + element.derivativeT1_j.resize(6); + element.derivativeT1_jplus1.resize(6); + element.derivativeT2_j.resize(6); + element.derivativeT2_jplus1.resize(6); + element.derivativeT3_j.resize(6); + element.derivativeT3_jplus1.resize(6); + element.derivativeR_j.resize(6); + element.derivativeR_jplus1.resize(6); + } +} + +void SimulationMesh::updateElementalLengths() { + for (const EdgeType &e : edge) { + const EdgeIndex ei = getIndex(e); + const VertexIndex vi0 = getIndex(e.cV(0)); + const VCGEdgeMesh::CoordType p0 = e.cP(0); + const VertexIndex vi1 = getIndex(e.cV(1)); + const VCGEdgeMesh::CoordType p1 = e.cP(1); + const vcg::Segment3 s(p0, p1); + const double elementLength = s.Length(); + elements[e].length = elementLength; + int i = 0; + i++; + } +} + +SimulationMesh::EdgePointer +SimulationMesh::getReferenceElement(const VCGEdgeMesh::VertexType &v) { + const VertexIndex vi = getIndex(v); + // return nodes[v].incidentElements[0]; + // if (vi == 0 || vi == 1) { + // return nodes[0].incidentElements[0]; + // } + + return nodes[v].incidentElements[0]; +} + +VectorType computeT1Vector(const SimulationMesh::EdgeType &e) { + return computeT1Vector(e.cP(0), e.cP(1)); +} + +VectorType computeT1Vector(const CoordType &p0, const CoordType &p1) { + const VectorType t1 = (p1 - p0).Normalize(); + return t1; +} + +Element::LocalFrame computeElementFrame(const CoordType &p0, + const CoordType &p1, + const VectorType &elementNormal) { + const VectorType t1 = computeT1Vector(p0, p1); + const VectorType t2 = (elementNormal ^ t1).Normalize(); + const VectorType t3 = (t1 ^ t2).Normalize(); + + return Element::LocalFrame{t1, t2, t3}; +} + +double computeAngle(const VectorType &vector0, const VectorType &vector1, + const VectorType &normalVector) { + double cosAngle = vector0.dot(vector1); + const double epsilon = std::pow(10, -6); + if (abs(cosAngle) > 1 && abs(cosAngle) < 1 + epsilon) { + if (cosAngle > 0) { + cosAngle = 1; + + } else { + cosAngle = -1; + } + } + assert(abs(cosAngle) <= 1); + const double angle = + acos(cosAngle); // NOTE: I compute the alpha angle not between + // two consecutive elements but rather between + // the first and the ith. Is this correct? + assert(!std::isnan(angle)); + + const VectorType cp = vector0 ^ vector1; + if (cp.dot(normalVector) < 0) { + return -angle; + } + return angle; +} diff --git a/elementalmesh.hpp b/elementalmesh.hpp new file mode 100644 index 0000000..682277b --- /dev/null +++ b/elementalmesh.hpp @@ -0,0 +1,164 @@ +#ifndef ELEMENTALMESH_HPP +#define ELEMENTALMESH_HPP + +#include "Eigen/Dense" +#include "edgemesh.hpp" +#include "flatpattern.hpp" + +struct Element; +struct Node; + +class SimulationMesh : public VCGEdgeMesh { +private: + void computeElementalProperties(); + void initializeNodes(); + void initializeElements(); + EdgePointer getReferenceElement(const VertexType &v); + +public: + PerEdgeAttributeHandle elements; + PerVertexAttributeHandle nodes; + SimulationMesh(FlatPattern &pattern); + SimulationMesh(ConstVCGEdgeMesh &edgeMesh); + SimulationMesh(SimulationMesh &elementalMesh); + void updateElementalLengths(); + + std::vector + getIncidentElements(const VCGEdgeMesh::VertexType &v); + + double previousTotalKineticEnergy{0}; + double previousTotalResidualForcesNorm{0}; + double currentTotalKineticEnergy{0}; + double totalResidualForcesNorm{0}; + double totalPotentialEnergykN{0}; +}; + +struct Element { + struct Properties { + double E{0}; // youngs modulus in pascal + double G{0}; // shear modulus + 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) { + E = material.youngsModulusGPascal * std::pow(10, 9); + G = E / (2 * (1 + material.poissonsRatio)); + } + void + computeDimensionsProperties(const RectangularBeamDimensions &dimensions) { + A = (dimensions.b * dimensions.h); + I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; + I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; + } + void computeDimensionsProperties( + const CylindricalElementDimensions &dimensions) { + 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; + } + Properties(const RectangularBeamDimensions &dimensions, + const ElementMaterial &material) { + computeDimensionsProperties(dimensions); + computeMaterialProperties(material); + J = I2 + I3; + } + Properties(const CylindricalElementDimensions &dimensions, + const ElementMaterial &material) { + computeDimensionsProperties(dimensions); + computeMaterialProperties(material); + J = I2 + I3; + } + Properties() {} + }; + + struct LocalFrame { + VectorType t1; + VectorType t2; + VectorType t3; + }; + + EdgeIndex ei; + double length{0}; + Properties properties; + double initialLength; + LocalFrame frame; + double axialConstFactor; + double torsionConstFactor; + double firstBendingConstFactor; + double secondBendingConstFactor; + VectorType f1_j; + VectorType f1_jplus1; + VectorType f2_j; + VectorType f2_jplus1; + VectorType f3_j; + VectorType f3_jplus1; + double cosRotationAngle_j; + double cosRotationAngle_jplus1; + double sinRotationAngle_j; + double sinRotationAngle_jplus1; + std::vector> derivativeT1; + std::vector> derivativeT2; + std::vector> derivativeT3; + std::vector derivativeT1_j; + std::vector derivativeT1_jplus1; + std::vector derivativeT2_j; + std::vector derivativeT2_jplus1; + std::vector derivativeT3_j; + std::vector derivativeT3_jplus1; + std::vector derivativeR_j; + std::vector derivativeR_jplus1; + struct RotationalDisplacements { + double theta1{0}, theta2{0}, theta3{0}; + }; + RotationalDisplacements rotationalDisplacements_j; + RotationalDisplacements rotationalDisplacements_jplus1; +}; + +struct Node { + struct Forces { + Vector6d external{0}; + Vector6d internal{0}; + Vector6d residual{0}; + Vector6d internalAxial{0}; + Vector6d internalFirstBending{0}; + Vector6d internalSecondBending{0}; + bool hasExternalForce() const { return external.isZero(); } + }; + + VertexIndex vi; + CoordType initialLocation; + CoordType previousLocation; + CoordType initialNormal; + double translationalMass; + double rotationalMass_I2; + double rotationalMass_I3; + double rotationalMass_J; + Vector6d acceleration{0}; + Forces force; + Vector6d velocity{0}; + double kineticEnergy{0}; + Vector6d displacements{0}; + double nR{0}; + std::unordered_map + alphaAngles; // contains the initial angles between the first star element + // incident to this node and the other elements of the star + // has size equal to the valence of the vertex + + std::vector incidentElements; + std::vector derivativeOfNormal; + SimulationMesh::EdgePointer referenceElement; +}; + +Element::LocalFrame computeElementFrame(const CoordType &p0, + const CoordType &p1, + const VectorType &elementNormal); +VectorType computeT1Vector(const ::SimulationMesh::EdgeType &e); + +VectorType computeT1Vector(const CoordType &p0, const CoordType &p1); +double computeAngle(const VectorType &vector0, const VectorType &vector1, + const VectorType &normalVector); + +#endif // ELEMENTALMESH_HPP diff --git a/flatpattern.cpp b/flatpattern.cpp new file mode 100644 index 0000000..381bbbf --- /dev/null +++ b/flatpattern.cpp @@ -0,0 +1,397 @@ +#include "flatpattern.hpp" +#include "trianglepatterngeometry.hpp" +#include + +FlatPattern::FlatPattern() {} + +FlatPattern::FlatPattern(const string &filename, bool addNormalsIfAbsent) { + assert(std::filesystem::exists(std::filesystem::path(filename))); + loadFromPly(filename); + if (addNormalsIfAbsent) { + bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; + if (normalsAreAbsent) { + for (auto &v : vert) { + v.N() = CoordType(0, 0, 1); + } + } + } + + vcg::tri::UpdateTopology::VertexEdge(*this); + + // scale(); +} + +FlatPattern::FlatPattern(const std::vector &numberOfNodesPerSlot, + const std::vector &edges) { + add(numberOfNodesPerSlot, edges); + // add normals + bool normalsAreAbsent = vert[0].cN().Norm() < 0.000001; + if (normalsAreAbsent) { + for (auto &v : vert) { + v.N() = CoordType(0, 0, 1); + } + } + updateEigenEdgeAndVertices(); +} + +bool FlatPattern::createHoneycombAtom() { + VCGEdgeMesh honeycombQuarter; + const VCGEdgeMesh::CoordType n(0, 0, 1); + const double H = 0.2; + const double height = 1.5 * H; + const double width = 0.2; + const double theta = 70; + const double dy = tan(vcg::math::ToRad(90 - theta)) * width / 2; + vcg::tri::Allocator::AddVertex( + honeycombQuarter, VCGEdgeMesh::CoordType(0, height / 2, 0), n); + vcg::tri::Allocator::AddVertex( + honeycombQuarter, VCGEdgeMesh::CoordType(0, H / 2 - dy, 0), n); + vcg::tri::Allocator::AddVertex( + honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, H / 2, 0), n); + vcg::tri::Allocator::AddVertex( + honeycombQuarter, VCGEdgeMesh::CoordType(width / 2, 0, 0), n); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 0, 1); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 1, 2); + vcg::tri::Allocator::AddEdge(honeycombQuarter, 2, 3); + + VCGEdgeMesh honeycombAtom; + // Top right + vcg::tri::Append::MeshCopy(honeycombAtom, + honeycombQuarter); + // Bottom right + vcg::Matrix44d rotM; + rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, + honeycombQuarter); + // Bottom left + rotM.SetRotateDeg(180, vcg::Point3d(0, 1, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, + honeycombQuarter); + // Top left + rotM.SetRotateDeg(180, vcg::Point3d(1, 0, 0)); + vcg::tri::UpdatePosition::Matrix(honeycombQuarter, rotM); + vcg::tri::Append::Mesh(honeycombAtom, + honeycombQuarter); + + for (VertexType &v : honeycombAtom.vert) { + v.P()[2] = 0; + } + + return true; +} + +void FlatPattern::deleteDanglingEdges() { + for (VertexType &v : vert) { + std::vector incidentElements; + vcg::edge::VEStarVE(&v, incidentElements); + if (incidentElements.size() == 1) { + vcg::tri::Allocator::DeleteEdge(*this, *incidentElements[0]); + } + if (incidentElements.size() == 1) { + vcg::tri::Allocator::DeleteVertex(*this, v); + } + } + + vcg::tri::Clean::RemoveDegenerateVertex(*this); + vcg::tri::Clean::RemoveDegenerateEdge(*this); + vcg::tri::Allocator::CompactEveryVector(*this); +} + +void FlatPattern::scale() { + const double baseTriangleCentralEdgeSize = + (vert[0].cP() - vert[3].cP()).Norm(); + const double scaleRatio = + desiredBaseTriangleCentralEdgeSize / baseTriangleCentralEdgeSize; + + vcg::tri::UpdatePosition::Scale(*this, scaleRatio); +} + +void FlatPattern::deleteDanglingVertices() { + vcg::tri::Allocator::PointerUpdater pu; + deleteDanglingVertices(pu); +} + +void FlatPattern::deleteDanglingVertices( + vcg::tri::Allocator::PointerUpdater &pu) { + for (VertexType &v : vert) { + std::vector incidentElements; + vcg::edge::VEStarVE(&v, incidentElements); + if (incidentElements.size() == 0) { + vcg::tri::Allocator::DeleteVertex(*this, v); + } + } + + vcg::tri::Allocator::CompactVertexVector(*this, pu); + vcg::tri::Allocator::CompactEdgeVector(*this); + + updateEigenEdgeAndVertices(); +} + +void FlatPattern::tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto, + const bool &shouldDeleteDanglingEdges) { + const size_t middleIndex = 3; + double xOffset = + vcg::Distance(pattern.vert[0].cP(), pattern.vert[middleIndex].cP()) / + std::tan(M_PI / 3); + CoordType patternCoord0 = pattern.vert[0].cP(); + CoordType patternBottomRight = + pattern.vert[middleIndex].cP() + CoordType(xOffset, 0, 0); + CoordType patternBottomLeft = + pattern.vert[middleIndex].cP() - CoordType(xOffset, 0, 0); + std::vector patternTrianglePoints{ + patternCoord0, patternBottomRight, patternBottomLeft}; + CoordType pointOnPattern = + patternCoord0 + (patternBottomLeft - patternCoord0) ^ + (patternBottomRight - patternCoord0); + + std::vector faceCenters(FN()); + VCGTriMesh tileIntoEdgeMesh; + + for (VCGPolyMesh::FaceType &f : tileInto.face) { + std::vector incidentVertices; + vcg::face::VFIterator vfi( + &f, 0); // initialize the iterator to the first face + // vcg::face::Pos p(vfi.F(), f.cV(0)); + // vcg::face::VVOrderedStarFF(p, incidentVertices); + size_t numberOfNodes = 0; + CoordType centerOfFace(0, 0, 0); + for (size_t vi = 0; vi < f.VN(); vi++) { + numberOfNodes++; + centerOfFace = centerOfFace + f.cP(vi); + } + centerOfFace /= f.VN(); + vcg::tri::Allocator::AddVertex(tileIntoEdgeMesh, centerOfFace, + vcg::Color4b::Yellow); + + // const size_t vi = vcg::tri::Index(tileInto, v); + // std::cout << "vertex " << vi << " has incident vertices:" << + // std::endl; + for (size_t vi = 0; vi < f.VN(); vi++) { + // size_t f = 0; + // std::cout << vcg::tri::Index(tileInto, + // / incidentVertices[f]) / << std::endl; + + // Compute transformation matrix M + // vcg::Matrix44d M; + std::vector meshTrianglePoints{ + centerOfFace, f.cP(vi), vi + 1 == f.VN() ? f.cP(0) : f.cP(vi + 1)}; + CoordType faceNormal = ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^ + (meshTrianglePoints[2] - meshTrianglePoints[0])) + .Normalize(); + auto fit = vcg::tri::Allocator::AddFace( + tileIntoEdgeMesh, meshTrianglePoints[0], meshTrianglePoints[1], + meshTrianglePoints[2]); + fit->N() = faceNormal; + CoordType pointOnTriMesh = + meshTrianglePoints[0] + + (meshTrianglePoints[1] - meshTrianglePoints[0]) ^ + (meshTrianglePoints[2] - meshTrianglePoints[0]); + vcg::Matrix44d M; + // vcg::ComputeRigidMatchMatrix(meshTrianglePoints, + // patternTrianglePoints, + // M); + vcg::Matrix44d A_prime; + A_prime[0][0] = meshTrianglePoints[0][0]; + A_prime[1][0] = meshTrianglePoints[0][1]; + A_prime[2][0] = meshTrianglePoints[0][2]; + A_prime[3][0] = 1; + A_prime[0][1] = meshTrianglePoints[1][0]; + A_prime[1][1] = meshTrianglePoints[1][1]; + A_prime[2][1] = meshTrianglePoints[1][2]; + A_prime[3][1] = 1; + A_prime[0][2] = meshTrianglePoints[2][0]; + A_prime[1][2] = meshTrianglePoints[2][1]; + A_prime[2][2] = meshTrianglePoints[2][2]; + A_prime[3][2] = 1; + A_prime[0][3] = pointOnTriMesh[0]; + A_prime[1][3] = pointOnTriMesh[1]; + A_prime[2][3] = pointOnTriMesh[2]; + A_prime[3][3] = 1; + vcg::Matrix44d A; + A[0][0] = patternTrianglePoints[0][0]; + A[1][0] = patternTrianglePoints[0][1]; + A[2][0] = patternTrianglePoints[0][2]; + A[3][0] = 1; + A[0][1] = patternTrianglePoints[1][0]; + A[1][1] = patternTrianglePoints[1][1]; + A[2][1] = patternTrianglePoints[1][2]; + A[3][1] = 1; + A[0][2] = patternTrianglePoints[2][0]; + A[1][2] = patternTrianglePoints[2][1]; + A[2][2] = patternTrianglePoints[2][2]; + A[3][2] = 1; + A[0][3] = pointOnPattern[0]; + A[1][3] = pointOnPattern[1]; + A[2][3] = pointOnPattern[2]; + A[3][3] = 1; + M = A_prime * vcg::Inverse(A); + + VCGEdgeMesh transformedPattern; + vcg::tri::Append::MeshCopy(transformedPattern, + pattern); + vcg::tri::UpdatePosition::Matrix(transformedPattern, M); + for (VertexType &v : transformedPattern.vert) { + v.N() = faceNormal; + } + + vcg::tri::Append::Mesh(*this, + transformedPattern); + } + } + + // vcg::tri::Clean::MergeCloseVertex(*this, 0.0000000001); + // vcg::tri::Clean::RemoveDegenerateVertex(*this); + // vcg::tri::Clean::RemoveDegenerateEdge(*this); + // vcg::tri::Allocator::CompactEveryVector(*this); + + vcg::tri::UpdateTopology::VertexEdge(*this); + + vcg::tri::Clean::MergeCloseVertex(*this, 0.0000000001); + deleteDanglingVertices(); + deleteDanglingEdges(); + vcg::tri::Clean::RemoveDegenerateVertex(*this); + vcg::tri::Clean::RemoveDegenerateEdge(*this); + vcg::tri::Allocator::CompactEveryVector(*this); + updateEigenEdgeAndVertices(); + savePly("tiledPattern.ply"); + + vcg::tri::Clean::MergeCloseVertex(tileIntoEdgeMesh, 0.0000000001); + vcg::tri::Clean::RemoveDegenerateVertex(tileIntoEdgeMesh); + vcg::tri::Clean::RemoveDegenerateEdge(tileIntoEdgeMesh); + tileIntoEdgeMesh.savePly("tileIntoTriMesh.ply"); +} + +void FlatPattern::createFan(const size_t &fanSize) { + FlatPattern rotatedPattern; + vcg::tri::Append::MeshCopy(rotatedPattern, *this); + for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) { + vcg::Matrix44d R; + auto rotationAxis = vcg::Point3d(0, 0, 1); + R.SetRotateDeg(360 / fanSize, rotationAxis); + vcg::tri::UpdatePosition::Matrix(rotatedPattern, R); + vcg::tri::Append::Mesh(*this, rotatedPattern); + } + + removeDuplicateVertices(); + updateEigenEdgeAndVertices(); +} + +void FlatPattern::removeDuplicateVertices() { + vcg::tri::Clean::MergeCloseVertex(*this, 0.0000000001); + vcg::tri::Clean::RemoveDegenerateVertex(*this); + vcg::tri::Clean::RemoveDegenerateEdge(*this); + vcg::tri::Allocator::CompactEveryVector(*this); + vcg::tri::UpdateTopology::VertexEdge(*this); +} + +void FlatPattern::tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto) { + + const size_t middleIndex = 3; + double xOffset = + vcg::Distance(pattern.vert[0].cP(), pattern.vert[middleIndex].cP()) / + std::tan(M_PI / 3); + CoordType patternCoord0 = pattern.vert[0].cP(); + CoordType patternBottomRight = + pattern.vert[middleIndex].cP() + CoordType(xOffset, 0, 0); + CoordType patternBottomLeft = + pattern.vert[middleIndex].cP() - CoordType(xOffset, 0, 0); + std::vector patternTrianglePoints{ + patternCoord0, patternBottomRight, patternBottomLeft}; + CoordType pointOnPattern = + patternCoord0 + (patternBottomLeft - patternCoord0) ^ + (patternBottomRight - patternCoord0); + + for (VCGTriMesh::VertexType &v : tileInto.vert) { + const auto centralNodeColor = vcg::Color4(64, 64, 64, 255); + const bool isCentralNode = v.cC() == centralNodeColor; + if (isCentralNode) { + std::vector incidentVertices; + vcg::face::VFIterator vfi( + &v); // initialize the iterator tohe first face + vcg::face::Pos p(vfi.F(), &v); + vcg::face::VVOrderedStarFF(p, incidentVertices); + + const size_t vi = vcg::tri::Index(tileInto, v); + std::cout << "vertex " << vi << " has incident vertices:" << std::endl; + for (size_t f = 0; f < incidentVertices.size(); f++) { + // size_t f = 0; + std::cout << vcg::tri::Index(tileInto, incidentVertices[f]) + << std::endl; + + // Compute transformation matrix M + // vcg::Matrix44d M; + std::vector meshTrianglePoints{ + v.cP(), incidentVertices[f]->cP(), + f + 1 == incidentVertices.size() ? incidentVertices[0]->cP() + : incidentVertices[f + 1]->cP()}; + CoordType faceNormal = + ((meshTrianglePoints[1] - meshTrianglePoints[0]) ^ + (meshTrianglePoints[2] - meshTrianglePoints[0])) + .Normalize(); + CoordType pointOnTriMesh = + meshTrianglePoints[0] + + (meshTrianglePoints[1] - meshTrianglePoints[0]) ^ + (meshTrianglePoints[2] - meshTrianglePoints[0]); + vcg::Matrix44d M; + // vcg::ComputeRigidMatchMatrix(meshTrianglePoints, + // patternTrianglePoints, + // M); + vcg::Matrix44d A_prime; + A_prime[0][0] = meshTrianglePoints[0][0]; + A_prime[1][0] = meshTrianglePoints[0][1]; + A_prime[2][0] = meshTrianglePoints[0][2]; + A_prime[3][0] = 1; + A_prime[0][1] = meshTrianglePoints[1][0]; + A_prime[1][1] = meshTrianglePoints[1][1]; + A_prime[2][1] = meshTrianglePoints[1][2]; + A_prime[3][1] = 1; + A_prime[0][2] = meshTrianglePoints[2][0]; + A_prime[1][2] = meshTrianglePoints[2][1]; + A_prime[2][2] = meshTrianglePoints[2][2]; + A_prime[3][2] = 1; + A_prime[0][3] = pointOnTriMesh[0]; + A_prime[1][3] = pointOnTriMesh[1]; + A_prime[2][3] = pointOnTriMesh[2]; + A_prime[3][3] = 1; + vcg::Matrix44d A; + A[0][0] = patternTrianglePoints[0][0]; + A[1][0] = patternTrianglePoints[0][1]; + A[2][0] = patternTrianglePoints[0][2]; + A[3][0] = 1; + A[0][1] = patternTrianglePoints[1][0]; + A[1][1] = patternTrianglePoints[1][1]; + A[2][1] = patternTrianglePoints[1][2]; + A[3][1] = 1; + A[0][2] = patternTrianglePoints[2][0]; + A[1][2] = patternTrianglePoints[2][1]; + A[2][2] = patternTrianglePoints[2][2]; + A[3][2] = 1; + A[0][3] = pointOnPattern[0]; + A[1][3] = pointOnPattern[1]; + A[2][3] = pointOnPattern[2]; + A[3][3] = 1; + M = A_prime * vcg::Inverse(A); + + VCGEdgeMesh transformedPattern; + vcg::tri::Append::MeshCopy(transformedPattern, + pattern); + vcg::tri::UpdatePosition::Matrix(transformedPattern, M); + for (VertexType &v : transformedPattern.vert) { + v.N() = faceNormal; + } + + vcg::tri::Append::Mesh(*this, + transformedPattern); + } + } + } + + vcg::tri::UpdateTopology::VertexEdge(*this); + deleteDanglingVertices(); + deleteDanglingEdges(); + vcg::tri::Allocator::CompactEveryVector(*this); + + updateEigenEdgeAndVertices(); +} diff --git a/flatpattern.hpp b/flatpattern.hpp new file mode 100644 index 0000000..47b42d5 --- /dev/null +++ b/flatpattern.hpp @@ -0,0 +1,33 @@ +#ifndef FLATPATTERN_HPP +#define FLATPATTERN_HPP + +#include "trianglepatterngeometry.hpp" + +class FlatPattern : public FlatPatternGeometry { +public: + FlatPattern(); + FlatPattern(const std::string &filename, bool addNormalsIfAbsent = true); + FlatPattern(const std::vector &numberOfNodesPerSlot, + const std::vector &edges); + + bool createHoneycombAtom(); + + void tilePattern(VCGEdgeMesh &pattern, VCGTriMesh &tileInto); + void tilePattern(VCGEdgeMesh &pattern, VCGPolyMesh &tileInto, + const bool &shouldDeleteDanglingEdges); + + void createFan(const size_t &fanSize = 6); + + void deleteDanglingVertices( + vcg::tri::Allocator::PointerUpdater &pu); + void deleteDanglingVertices(); + +private: + void deleteDanglingEdges(); + void removeDuplicateVertices(); + void scale(); + + const double desiredBaseTriangleCentralEdgeSize{0.25 / std::tan(M_PI / 6)}; +}; + +#endif // FLATPATTERN_HPP diff --git a/patternIO.cpp b/patternIO.cpp new file mode 100644 index 0000000..10bb10e --- /dev/null +++ b/patternIO.cpp @@ -0,0 +1,116 @@ +#include "patternIO.hpp" +#include +#include +#include +#include + +PatternIO::PatternIO() {} + +void PatternIO::save(const std::string &filepath, + const PatternSet &patternSet) { + std::ofstream fileStream; + if (std::filesystem::exists(filepath)) { + fileStream.open(filepath, std::ios_base::app); + } else { + fileStream.open(filepath); + fileStream << "#Nodes" << std::endl; + for (vcg::Point2d node : patternSet.nodes) { + fileStream << "n " << node.X() << " " << node.Y() << std::endl; + } + fileStream << "#Patterns" << std::endl; + } + + for (const Pattern &pattern : patternSet.patterns) { + fileStream << "p " << pattern.labels.size() << " " << pattern.edges.size() + << " "; + for (size_t labelIndex = 0; labelIndex < pattern.labels.size() - 1; + labelIndex++) { + fileStream << pattern.labels[labelIndex] << " "; + } + fileStream << pattern.labels[pattern.labels.size() - 1] << " "; + + for (const vcg::Point2i &edge : pattern.edges) { + fileStream << " " << edge[0] << " " << edge[1] << " "; + } + fileStream << std::endl; + } +} + +void PatternIO::save(const std::string &filepath, const Pattern &pattern) { + std::ofstream fileStream; + if (std::filesystem::exists(filepath)) { + fileStream.open(filepath, std::ios_base::app); + } else { + fileStream.open(filepath); + fileStream << "#Nodes" << std::endl; + fileStream << "#Patterns" << std::endl; + } + + fileStream << "p " << pattern.labels.size() << " " << pattern.edges.size() + << " "; + for (size_t labelIndex = 0; labelIndex < pattern.labels.size() - 1; + labelIndex++) { + fileStream << pattern.labels[labelIndex] << " "; + } + fileStream << pattern.labels[pattern.labels.size() - 1] << " "; + + for (const vcg::Point2i &edge : pattern.edges) { + fileStream << " " << edge[0] << " " << edge[1] << " "; + } + fileStream << std::endl; +} + +void PatternIO::load(const std::string &filepath, PatternSet &patternSet, + const std::vector &targetLabels) { + std::ifstream fileStream(filepath); + std::string line; + std::vector sortedTargetPatternLabels(targetLabels); + std::sort(sortedTargetPatternLabels.begin(), sortedTargetPatternLabels.end()); + while (std::getline(fileStream, line)) { + std::cout << line << std::endl; + std::istringstream iss(line); + char lineID; + iss >> lineID; + if (lineID == 'n') { + double x, y; + iss >> x >> y; + std::cout << x << " " << y << std::endl; + } else if (lineID == 'p') { + Pattern pattern; + size_t numberOfLabels, numberOfEdges; + iss >> numberOfLabels >> numberOfEdges; + std::cout << "NumberOfLabels:" << numberOfLabels << std::endl; + std::cout << "NumberOfEdges:" << numberOfEdges << std::endl; + std::vector patternLabels; + for (size_t labelIndex = 0; labelIndex < numberOfLabels; labelIndex++) { + size_t patternLabel; + iss >> patternLabel; + PatternLabel pl = static_cast(patternLabel); + std::cout << "Pattern label read:" << patternLabel << std::endl; + patternLabels.push_back(pl); + } + if (!targetLabels.empty()) { + std::sort(patternLabels.begin(), patternLabels.end()); + std::vector labelIntersection; + std::set_intersection(patternLabels.begin(), patternLabels.end(), + sortedTargetPatternLabels.begin(), + sortedTargetPatternLabels.end(), + labelIntersection.begin()); + if (!labelIntersection.empty()) { + pattern.labels = patternLabels; + } else { + continue; + } + } else { + pattern.labels = patternLabels; + } + for (size_t edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) { + size_t ni0, ni1; + iss >> ni0 >> ni1; + vcg::Point2i edge(ni0, ni1); + pattern.edges.push_back(edge); + std::cout << "Node indices read:" << ni0 << "," << ni1 << std::endl; + } + } + } +} diff --git a/patternIO.hpp b/patternIO.hpp new file mode 100644 index 0000000..daef17a --- /dev/null +++ b/patternIO.hpp @@ -0,0 +1,40 @@ +#ifndef PATTERNEXPORTER_HPP +#define PATTERNEXPORTER_HPP + +#include +#include +#include +#include + +enum PatternLabel { + Valid = 0, + MultipleCC, + DanglingEdge, + IntersectingEdges, + ArticulationPoints +}; + +struct Pattern { + std::vector edges; + std::vector labels; +}; + +/* + * A set of planar patterns using the same nodes + * */ +struct PatternSet { + std::vector nodes; + std::vector patterns; +}; + +class PatternIO { + +public: + PatternIO(); + static void save(const std::string &filepath, const Pattern &pattern); + static void save(const std::string &filepath, const PatternSet &patterns); + static void load(const std::string &filepath, PatternSet &patternSet, + const std::vector &targetLabels = {}); +}; + +#endif // PATTERNEXPORTER_HPP diff --git a/polymesh.hpp b/polymesh.hpp new file mode 100644 index 0000000..91256d8 --- /dev/null +++ b/polymesh.hpp @@ -0,0 +1,64 @@ +#ifndef POLYMESH_HPP +#define POLYMESH_HPP +#include "vcg/complex/complex.h" +#include +#include +//#include + +class PFace; +class PVertex; + +struct PUsedTypes : public vcg::UsedTypes::AsVertexType, + vcg::Use::AsFaceType> {}; + +class PVertex : public vcg::Vertex {}; + +class PFace + : public vcg::Face< + PUsedTypes, + vcg::face::PolyInfo // this is necessary if you use component in + // vcg/simplex/face/component_polygon.h + // It says "this class is a polygon and the memory for its components + // (e.g. pointer to its vertices will be allocated dynamically") + , + vcg::face::PFVAdj // Pointer to the vertices (just like FVAdj ) + , + vcg::face::PFVAdj, + vcg::face::PFFAdj // Pointer to edge-adjacent face (just like FFAdj ) + , + vcg::face::BitFlags // bit flags + , + vcg::face::Qualityf // quality + , + vcg::face::Normal3f // normal + > {}; + +class VCGPolyMesh + : public vcg::tri::TriMesh, // the vector of vertices + std::vector // the vector of faces + > { +public: + void loadFromPlyFile(const std::string &filename) { + vcg::tri::io::ImporterOBJ::Info info; + vcg::tri::io::ImporterOBJ::Open(*this, filename.c_str(), info); + // 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; + // if (nanoply::NanoPlyWrapper::LoadModel( + // std::filesystem::absolute(filename).c_str(), *this, mask) != + // 0) { + // std::cout << "Could not load tri mesh" << std::endl; + // } + vcg::tri::UpdateTopology::FaceFace(*this); + // vcg::tri::UpdateTopology::VertexFace(*this); + vcg::tri::UpdateNormal::PerVertexNormalized(*this); + } +}; + +#endif // POLYMESH_HPP diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp new file mode 100644 index 0000000..6d5ffc5 --- /dev/null +++ b/simulationhistoryplotter.hpp @@ -0,0 +1,145 @@ +#ifndef SIMULATIONHISTORYPLOTTER_HPP +#define SIMULATIONHISTORYPLOTTER_HPP + +#include "elementalmesh.hpp" +#include "simulationresult.hpp" +#include "utilities.hpp" +#include +#include + +struct SimulationResultsReporter { + using VertexType = VCGEdgeMesh::VertexType; + using CoordType = VCGEdgeMesh::CoordType; + using Vector6d = Vector6d; + + SimulationResultsReporter() {} + + void writeStatistics(const SimulationResults &results, + const std::string &reportFolderPath) { + + ofstream file; + file.open( + std::filesystem::path(reportFolderPath).append("results.txt").string()); + + const size_t numberOfSteps = results.history.numberOfSteps; + file << "Number of steps " << numberOfSteps << "\n"; + + // file << "Force threshold used " << 1000 << "\n"; + + // assert(numberOfSteps == results.history.potentialEnergy.size() && + // numberOfSteps == results.history.residualForces.size()); + // Write kinetic energies + const SimulationHistory &history = results.history; + if (!history.kineticEnergy.empty()) { + file << "Kinetic energies" + << "\n"; + for (size_t step = 0; step < numberOfSteps; step++) { + file << history.kineticEnergy[step] << "\n"; + } + file << "\n"; + } + + if (!history.residualForces.empty()) { + file << "Residual forces" + << "\n"; + for (size_t step = 0; step < numberOfSteps; step++) { + file << history.residualForces[step] << "\n"; + } + file << "\n"; + } + + if (!history.potentialEnergies.empty()) { + file << "Potential energies" + << "\n"; + for (size_t step = 0; step < numberOfSteps; step++) { + file << history.potentialEnergies[step] << "\n"; + } + file << "\n"; + } + file.close(); + } + + void reportResults( + const std::vector &results, + const std::string &reportFolderPath, + const std::string &graphSuffix = std::string(), + const SimulationResults &groundOfTruthResults = SimulationResults()) { + if (results.empty()) { + return; + } + + // std::filesystem::remove_all(debuggerFolder); + std::filesystem::create_directory(reportFolderPath); + for (const SimulationResults &simulationResult : results) { + const auto simulationResultPath = + std::filesystem::path(reportFolderPath) + .append(simulationResult.simulationLabel); + std::filesystem::create_directory(simulationResultPath.string()); + + createPlots(simulationResult.history, simulationResultPath.string(), + graphSuffix); + writeStatistics(simulationResult, simulationResultPath); + } + } + + static void createPlot(const std::string &xLabel, const std::string &yLabel, + const std::vector &YvaluesToPlot, + const std::string &saveTo = {}) { + matplot::xlabel(xLabel); + matplot::ylabel(yLabel); + matplot::figure(true); + // matplot::hold(matplot::on); + matplot::grid(matplot::on); + auto x = matplot::linspace(0, YvaluesToPlot.size(), YvaluesToPlot.size()); + matplot::scatter(x, YvaluesToPlot) + // ->marker_indices(history.redMarks) + // ->marker_indices(truncatedRedMarks) + // .marker_color("r") + ->marker_size(1); + // auto greenMarksY = matplot::transform( + // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; + // }); + // matplot::scatter(history.greenMarks, greenMarksY) + // ->color("green") + // .marker_size(10); + // matplot::hold(matplot::off); + if (!saveTo.empty()) { + matplot::save(saveTo); + } + } + + void createPlots(const SimulationHistory &history, + const std::string &reportFolderPath, + const std::string &graphSuffix) { + const auto graphsFolder = + std::filesystem::path(reportFolderPath).append("Graphs"); + std::filesystem::remove_all(graphsFolder); + std::filesystem::create_directory(graphsFolder.string()); + + if (!history.kineticEnergy.empty()) { + createPlot("Number of Iterations", "Log of Kinetic Energy", + history.kineticEnergy, + std::filesystem::path(graphsFolder) + .append("KineticEnergy" + graphSuffix + ".png") + .string()); + } + + if (!history.residualForces.empty()) { + createPlot("Number of Iterations", "Residual Forces norm", + history.residualForces, + std::filesystem::path(graphsFolder) + .append("ResidualForces" + graphSuffix + ".png") + .string()); + } + + if (!history.potentialEnergies.empty()) { + createPlot("Number of Iterations", "Potential energy", + history.potentialEnergies, + std::filesystem::path(graphsFolder) + .append("PotentialEnergy" + graphSuffix + ".png") + .string()); + } + } +}; + +#endif // SIMULATIONHISTORYPLOTTER_HPP diff --git a/simulationresult.hpp b/simulationresult.hpp new file mode 100644 index 0000000..33c990c --- /dev/null +++ b/simulationresult.hpp @@ -0,0 +1,199 @@ +#ifndef SIMULATIONHISTORY_HPP +#define SIMULATIONHISTORY_HPP + +#include "elementalmesh.hpp" + +struct SimulationHistory { + SimulationHistory() {} + + size_t numberOfSteps{0}; + std::string label; + std::vector residualForces; + std::vector kineticEnergy; + std::vector potentialEnergies; + std::vector redMarks; + std::vector greenMarks; + + void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); } + + void markGreen(const size_t &stepNumber) { greenMarks.push_back(stepNumber); } + + void stepPulse(const SimulationMesh &mesh) { + kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy)); + // potentialEnergy.push_back(mesh.totalPotentialEnergykN); + residualForces.push_back(mesh.totalResidualForcesNorm); + } + + void clear() { + residualForces.clear(); + kineticEnergy.clear(); + potentialEnergies.clear(); + } +}; + +struct SimulationJob { + std::shared_ptr mesh; + std::unordered_map> fixedVertices; + std::unordered_map nodalExternalForces; + std::unordered_map nodalForcedDisplacements; + + void registerForDrawing() const { + initPolyscope(); + const std::string polyscopeName = mesh->getLabel() + "_Simulation Job"; + polyscope::registerCurveNetwork(polyscopeName, mesh->getEigenVertices(), + mesh->getEigenEdges()); + // polyscope::registerPointCloud("Undeformed edge mesh PC", + // job.edgeMesh.getEigenVertices()); + + std::vector> nodeColors(mesh->VN()); + for (auto fixedVertex : fixedVertices) { + nodeColors[fixedVertex.first] = {0, 0, 1}; + } + if (!nodalForcedDisplacements.empty()) { + for (std::pair viDisplPair : + nodalForcedDisplacements) { + const VertexIndex vi = viDisplPair.first; + nodeColors[vi][0] += 1; + nodeColors[vi][0] /= 2; + nodeColors[vi][1] += 0; + nodeColors[vi][1] /= 2; + nodeColors[vi][2] += 0; + nodeColors[vi][2] /= 2; + } + } + std::for_each(nodeColors.begin(), nodeColors.end(), + [](std::array &color) { + const double norm = + sqrt(std::pow(color[0], 2) + std::pow(color[1], 2) + + std::pow(color[2], 2)); + if (norm > std::pow(10, -7)) { + color[0] /= norm; + color[1] /= norm; + color[2] /= norm; + } + }); + + if (!nodeColors.empty()) { + polyscope::getCurveNetwork(polyscopeName) + ->addNodeColorQuantity("Boundary conditions", nodeColors); + polyscope::getCurveNetwork(polyscopeName) + ->getQuantity("Boundary conditions") + ->setEnabled(true); + } + + // per node external forces + std::vector> externalForces(mesh->VN()); + for (const auto &forcePair : nodalExternalForces) { + auto index = forcePair.first; + auto force = forcePair.second; + externalForces[index] = {force[0], force[1], force[2]}; + } + + if (!externalForces.empty()) { + polyscope::getCurveNetwork(polyscopeName) + ->addNodeVectorQuantity("External force", externalForces); + polyscope::getCurveNetwork(polyscopeName) + ->getQuantity("External force") + ->setEnabled(true); + } + } +}; + +struct SimulationResults { + SimulationHistory history; + std::vector displacements; + double executionTime{0}; + std::string simulationLabel{"NoLabel"}; + + void registerForDrawing(const SimulationJob &job) const { + polyscope::options::groundPlaneEnabled = false; + polyscope::view::upDir = polyscope::view::UpDir::ZUp; + const std::string branchName = "Branch:Polyscope"; + polyscope::options::programName = branchName; + if (!polyscope::state::initialized) { + polyscope::init(); + } /* else { + polyscope::removeAllStructures(); + }*/ + const std::string undeformedMeshName = "Undeformed_" + simulationLabel; + polyscope::registerCurveNetwork(undeformedMeshName, + job.mesh->getEigenVertices(), + job.mesh->getEigenEdges()); + + const std::string deformedMeshName = "Deformed_" + simulationLabel; + polyscope::registerCurveNetwork(deformedMeshName, + job.mesh->getEigenVertices(), + job.mesh->getEigenEdges()); + Eigen::MatrixX3d nodalDisplacements(job.mesh->VN(), 3); + for (VertexIndex vi = 0; vi < job.mesh->VN(); vi++) { + const Vector6d &nodalDisplacement = displacements[vi]; + nodalDisplacements.row(vi) = Eigen::Vector3d( + nodalDisplacement[0], nodalDisplacement[1], nodalDisplacement[2]); + } + polyscope::getCurveNetwork(deformedMeshName) + ->updateNodePositions(job.mesh->getEigenVertices() + + nodalDisplacements); + + std::vector> nodeColors(job.mesh->VN()); + for (auto fixedVertex : job.fixedVertices) { + nodeColors[fixedVertex.first] = {0, 0, 1}; + } + if (!job.nodalForcedDisplacements.empty()) { + for (std::pair viDisplPair : + job.nodalForcedDisplacements) { + const VertexIndex vi = viDisplPair.first; + nodeColors[vi][0] += 1; + nodeColors[vi][0] /= 2; + nodeColors[vi][1] += 0; + nodeColors[vi][1] /= 2; + nodeColors[vi][2] += 0; + nodeColors[vi][2] /= 2; + } + } + std::for_each(nodeColors.begin(), nodeColors.end(), + [](std::array &color) { + const double norm = + sqrt(std::pow(color[0], 2) + std::pow(color[1], 2) + + std::pow(color[2], 2)); + if (norm > std::pow(10, -7)) { + color[0] /= norm; + color[1] /= norm; + color[2] /= norm; + } + }); + // per node external forces + std::vector> externalForces(job.mesh->VN()); + for (const auto &forcePair : job.nodalExternalForces) { + auto index = forcePair.first; + auto force = forcePair.second; + externalForces[index] = {force[0], force[1], force[2]}; + } + + polyscope::getCurveNetwork(undeformedMeshName) + ->addNodeColorQuantity("Boundary conditions", nodeColors); + polyscope::getCurveNetwork(undeformedMeshName) + ->getQuantity("Boundary conditions") + ->setEnabled(true); + + polyscope::getCurveNetwork(undeformedMeshName) + ->addNodeVectorQuantity("External force", externalForces); + polyscope::getCurveNetwork(undeformedMeshName) + ->getQuantity("External force") + ->setEnabled(true); + polyscope::getCurveNetwork(deformedMeshName) + ->addNodeColorQuantity("Boundary conditions", nodeColors); + polyscope::getCurveNetwork(deformedMeshName) + ->getQuantity("Boundary conditions") + ->setEnabled(false); + + polyscope::getCurveNetwork(deformedMeshName) + ->addNodeVectorQuantity("External force", externalForces); + polyscope::getCurveNetwork(deformedMeshName) + ->getQuantity("External force") + ->setEnabled(true); + + // polyscope::show(); + } +}; + +#endif // SIMULATIONHISTORY_HPP diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp new file mode 100644 index 0000000..04a85ff --- /dev/null +++ b/topologyenumerator.cpp @@ -0,0 +1,609 @@ +#include "topologyenumerator.hpp" +#include +#include +#include +#include +#include + +const bool debugIsOn{false}; +const bool exportArticulationPointsPatterns{false}; +const bool savePlyFiles{true}; + +// size_t binomialCoefficient(size_t n, size_t m) { +// assert(n > m); +// return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1)); +//} + +// void TopologyEnumerator::createLabelMesh( +// const std::vector vertices, +// const std::filesystem::path &savePath) const { +// const std::string allOnes(patternTopology.getNumberOfPossibleEdges(), '1'); +// const std::vector allEdges = +// TrianglePatternTopology::convertToEdges(allOnes, vertices.size()); +// TrianglePatternGeometry labelMesh; +// std::vector labelVertices(allEdges.size()); +// for (size_t edgeIndex = 0; edgeIndex < allEdges.size(); edgeIndex++) { +// const vcg::Point3d edgeMidpoint = +// (vertices[allEdges[edgeIndex][0]] + vertices[allEdges[edgeIndex][1]]) +// / 2; +// labelVertices[edgeIndex] = edgeMidpoint; +// } +// labelMesh.set(labelVertices); +// labelMesh.savePly(std::filesystem::path(savePath) +// .append(std::string("labelMesh.ply")) +// .string()); +//} + +size_t TopologyEnumerator::getEdgeIndex(size_t ni0, size_t ni1) const { + if (ni1 <= ni0) { + std::swap(ni0, ni1); + } + assert(ni1 > ni0); + const size_t &n = numberOfNodes; + return (n * (n - 1) / 2) - (n - ni0) * ((n - ni0) - 1) / 2 + ni1 - ni0 - 1; +} + +TopologyEnumerator::TopologyEnumerator() {} + +void TopologyEnumerator::computeValidPatterns( + const std::vector &reducedNumberOfNodesPerSlot) { + assert(reducedNumberOfNodesPerSlot.size() == 5); + assert(reducedNumberOfNodesPerSlot[0] == 0 || + reducedNumberOfNodesPerSlot[0] == 1); + assert(reducedNumberOfNodesPerSlot[1] == 0 || + reducedNumberOfNodesPerSlot[1] == 1); + std::vector numberOfNodesPerSlot{ + reducedNumberOfNodesPerSlot[0], reducedNumberOfNodesPerSlot[1], + reducedNumberOfNodesPerSlot[1], reducedNumberOfNodesPerSlot[2], + reducedNumberOfNodesPerSlot[3], reducedNumberOfNodesPerSlot[2], + reducedNumberOfNodesPerSlot[4]}; + // Generate an edge mesh wih all possible edges + numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.end(), 0); + const size_t numberOfAllPossibleEdges = + numberOfNodes * (numberOfNodes - 1) / 2; + + std::vector allPossibleEdges(numberOfAllPossibleEdges); + const int &n = numberOfNodes; + for (size_t edgeIndex = 0; edgeIndex < numberOfAllPossibleEdges; + edgeIndex++) { + const int ni0 = + n - 2 - + std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); + const int ni1 = + edgeIndex + ni0 + 1 - n * (n - 1) / 2 + (n - ni0) * ((n - ni0) - 1) / 2; + allPossibleEdges[edgeIndex] = vcg::Point2i(ni0, ni1); + } + FlatPatternGeometry patternGeometryAllEdges; + patternGeometryAllEdges.add(numberOfNodesPerSlot, allPossibleEdges); + // Create Results path + auto resultPath = + // std::filesystem::path("/home/iason/Documents/PhD/Research/Enumerating\\ + // " + // "2d\\ connections\\ of\\ nodes"); + std::filesystem::current_path() + .parent_path() + .parent_path() + .parent_path() + .parent_path(); + assert(std::filesystem::exists(resultPath)); + + auto allResultsPath = resultPath.append("Results"); + std::filesystem::create_directory(allResultsPath); + std::string setupString; + // for (size_t numberOfNodes : reducedNumberOfNodesPerSlot) { + for (size_t numberOfNodesPerSlotIndex = 0; + numberOfNodesPerSlotIndex < reducedNumberOfNodesPerSlot.size(); + numberOfNodesPerSlotIndex++) { + std::string elemID; + if (numberOfNodesPerSlotIndex == 0 || numberOfNodesPerSlotIndex == 1) { + elemID = "v"; + } else if (numberOfNodesPerSlotIndex == 2 || + numberOfNodesPerSlotIndex == 3) { + elemID = "e"; + } else { + elemID = "c"; + } + setupString += + std::to_string(reducedNumberOfNodesPerSlot[numberOfNodesPerSlotIndex]) + + elemID + "_"; + } + setupString += std::to_string(FlatPatternGeometry().getFanSize()) + "fan"; + if (debugIsOn) { + setupString += "_debug"; + } + auto resultsPath = std::filesystem::path(allResultsPath).append(setupString); + // std::filesystem::remove_all(resultsPath); // delete previous results + std::filesystem::create_directory(resultsPath); + if (debugIsOn) { + patternGeometryAllEdges.savePly(std::filesystem::path(resultsPath) + .append("allPossibleEdges.ply") + .string()); + } + // statistics.numberOfPossibleEdges = numberOfAllPossibleEdges; + + std::vector validEdges = + getValidEdges(numberOfNodesPerSlot, resultsPath, patternGeometryAllEdges, + allPossibleEdges); + FlatPatternGeometry patternAllValidEdges; + patternAllValidEdges.add(patternGeometryAllEdges.getVertices(), validEdges); + if (debugIsOn) { + // Export all valid edges in a ply + patternAllValidEdges.savePly( + std::filesystem::path(resultsPath).append("allValidEdges.ply")); + } + // statistics.numberOfValidEdges = validEdges.size(); + + // Find pairs of intersecting edges + std::unordered_map> intersectingEdges = + patternAllValidEdges.getIntersectingEdges( + statistics.numberOfIntersectingEdgePairs); + if (debugIsOn) { + auto intersectingEdgesPath = std::filesystem::path(resultsPath) + .append("All_intersecting_edge_pairs"); + std::filesystem::create_directory(intersectingEdgesPath); + // Export intersecting pairs in ply files + for (auto mapIt = intersectingEdges.begin(); + mapIt != intersectingEdges.end(); mapIt++) { + for (auto setIt = mapIt->second.begin(); setIt != mapIt->second.end(); + setIt++) { + FlatPatternGeometry intersectingEdgePair; + const size_t ei0 = mapIt->first; + const size_t ei1 = *setIt; + vcg::tri::Allocator::AddEdge( + intersectingEdgePair, + patternGeometryAllEdges.getVertices()[validEdges[ei0][0]], + patternGeometryAllEdges.getVertices()[validEdges[ei0][1]]); + vcg::tri::Allocator::AddEdge( + intersectingEdgePair, + patternGeometryAllEdges.getVertices()[validEdges[ei1][0]], + patternGeometryAllEdges.getVertices()[validEdges[ei1][1]]); + intersectingEdgePair.savePly( + std::filesystem::path(intersectingEdgesPath) + .append(std::to_string(mapIt->first) + "_" + + std::to_string(*setIt) + ".ply") + .string()); + } + } + } + + // assert(validEdges.size() == allPossibleEdges.size() - + // coincideEdges.size() - + // duplicateEdges.size()); + + PatternSet patternSet; + const std::vector nodes = patternGeometryAllEdges.getVertices(); + const size_t numberOfNodes = nodes.size(); + patternSet.nodes.resize(numberOfNodes); + for (size_t nodeIndex = 0; nodeIndex < numberOfNodes; nodeIndex++) { + patternSet.nodes[nodeIndex] = + vcg::Point2d(nodes[nodeIndex][0], nodes[nodeIndex][1]); + } + if (std::filesystem::exists(std::filesystem::path(resultsPath) + .append("patterns.patt") + .string())) { + std::filesystem::remove( + std::filesystem::path(resultsPath).append("patterns.patt")); + } + for (size_t numberOfEdges = 2; numberOfEdges < validEdges.size(); + numberOfEdges++) { + // for (size_t numberOfEdges = 1; numberOfEdges < 3; numberOfEdges++) { + std::cout << "Computing " + setupString << " with " << numberOfEdges + << " edges." << std::endl; + auto perEdgeResultPath = std::filesystem::path(resultsPath) + .append(std::to_string(numberOfEdges)); + // if (std::filesystem::exists(perEdgeResultPath)) { + // continue; + // } + std::filesystem::create_directory(perEdgeResultPath); + computeValidPatterns(numberOfNodesPerSlot, numberOfEdges, perEdgeResultPath, + patternGeometryAllEdges.getVertices(), + intersectingEdges, validEdges, patternSet); + // statistics.print(setupString, perEdgeResultPath); + PatternIO::save( + std::filesystem::path(resultsPath).append("patterns.patt").string(), + patternSet); + } +} + +void TopologyEnumerator::computeEdgeNodes( + const std::vector &numberOfNodesPerSlot, + std::vector &nodesEdge0, std::vector &nodesEdge1, + std::vector &nodesEdge2) { + // Create vectors holding the node indices of each pattern node of each + // triangle edge + size_t nodeIndex = 0; + if (numberOfNodesPerSlot[0] != 0) { + nodesEdge0.push_back(nodeIndex++); + } + if (numberOfNodesPerSlot[1] != 0) + nodesEdge1.push_back(nodeIndex++); + if (numberOfNodesPerSlot[2] != 0) + nodesEdge2.push_back(nodeIndex++); + + if (numberOfNodesPerSlot[3] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[3]; + edgeNodeIndex++) { + nodesEdge0.push_back(nodeIndex++); + } + } + if (numberOfNodesPerSlot[4] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[4]; + edgeNodeIndex++) { + nodesEdge1.push_back(nodeIndex++); + } + } + + if (numberOfNodesPerSlot[5] != 0) { + for (size_t edgeNodeIndex = 0; edgeNodeIndex < numberOfNodesPerSlot[5]; + edgeNodeIndex++) { + nodesEdge2.push_back(nodeIndex++); + } + } + if (numberOfNodesPerSlot[1] != 0) { + assert(numberOfNodesPerSlot[2]); + nodesEdge0.push_back(1); + nodesEdge1.push_back(2); + } + + if (numberOfNodesPerSlot[0] != 0) { + nodesEdge2.push_back(0); + } +} + +std::unordered_set TopologyEnumerator::computeCoincideEdges( + const std::vector &numberOfNodesPerSlot) { + /* + * A coincide edge is defined as an edge connection between two nodes that lay + * on a triangle edge and which have another node in between + * */ + std::vector nodesEdge0; // left edge + std::vector nodesEdge1; // bottom edge + std::vector nodesEdge2; // right edge + computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); + + std::vector coincideEdges0 = getCoincideEdges(nodesEdge0); + std::vector coincideEdges1 = getCoincideEdges(nodesEdge1); + std::vector coincideEdges2 = getCoincideEdges(nodesEdge2); + std::unordered_set coincideEdges{coincideEdges0.begin(), + coincideEdges0.end()}; + std::copy(coincideEdges1.begin(), coincideEdges1.end(), + std::inserter(coincideEdges, coincideEdges.end())); + std::copy(coincideEdges2.begin(), coincideEdges2.end(), + std::inserter(coincideEdges, coincideEdges.end())); + + if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[1]) { + coincideEdges.insert(getEdgeIndex(0, 2)); + } + + if (numberOfNodesPerSlot[0] && numberOfNodesPerSlot[2]) { + assert(numberOfNodesPerSlot[1]); + coincideEdges.insert(getEdgeIndex(0, 3)); + } + + return coincideEdges; +} + +std::unordered_set TopologyEnumerator::computeDuplicateEdges( + const std::vector &numberOfNodesPerSlot) { + /* + * A duplicate edges are all edges the "right" edge since due to rotational + * symmetry "left" edge=="right" edge + * */ + std::unordered_set duplicateEdges; + std::vector nodesEdge0; // left edge + std::vector nodesEdge1; // bottom edge + std::vector nodesEdge2; // right edge + computeEdgeNodes(numberOfNodesPerSlot, nodesEdge0, nodesEdge1, nodesEdge2); + if (numberOfNodesPerSlot[5]) { + for (size_t edge2NodeIndex = 0; edge2NodeIndex < nodesEdge2.size() - 1; + edge2NodeIndex++) { + const size_t nodeIndex = nodesEdge2[edge2NodeIndex]; + const size_t nextNodeIndex = nodesEdge2[edge2NodeIndex + 1]; + duplicateEdges.insert(getEdgeIndex(nodeIndex, nextNodeIndex)); + } + } + + return duplicateEdges; +} + +std::vector TopologyEnumerator::getValidEdges( + const std::vector &numberOfNodesPerSlot, + const std::filesystem::path &resultsPath, + const FlatPatternGeometry &patternGeometryAllEdges, + const std::vector &allPossibleEdges) { + + std::unordered_set coincideEdges = + computeCoincideEdges(numberOfNodesPerSlot); + // Export each coincide edge into a ply file + if (!coincideEdges.empty() && debugIsOn) { + auto coincideEdgesPath = + std::filesystem::path(resultsPath).append("Coincide_edges"); + std::filesystem::create_directories(coincideEdgesPath); + for (auto coincideEdgeIndex : coincideEdges) { + FlatPatternGeometry::EdgeType e = + patternGeometryAllEdges.edge[coincideEdgeIndex]; + FlatPatternGeometry singleEdgeMesh; + vcg::Point3d p0 = e.cP(0); + vcg::Point3d p1 = e.cP(1); + std::vector edgeVertices; + edgeVertices.push_back(p0); + edgeVertices.push_back(p1); + singleEdgeMesh.add(edgeVertices); + singleEdgeMesh.add(std::vector{vcg::Point2i{0, 1}}); + singleEdgeMesh.savePly(std::filesystem::path(coincideEdgesPath) + .append(std::to_string(coincideEdgeIndex)) + .string() + + ".ply"); + } + } + statistics.numberOfCoincideEdges = coincideEdges.size(); + + // Compute duplicate edges + std::unordered_set duplicateEdges = + computeDuplicateEdges(numberOfNodesPerSlot); + if (!duplicateEdges.empty() && debugIsOn) { + // Export duplicate edges in a single ply file + auto duplicateEdgesPath = + std::filesystem::path(resultsPath).append("duplicate"); + std::filesystem::create_directory(duplicateEdgesPath); + FlatPatternGeometry patternDuplicateEdges; + for (auto duplicateEdgeIndex : duplicateEdges) { + FlatPatternGeometry::EdgeType e = + patternGeometryAllEdges.edge[duplicateEdgeIndex]; + vcg::Point3d p0 = e.cP(0); + vcg::Point3d p1 = e.cP(1); + vcg::tri::Allocator::AddEdge( + patternDuplicateEdges, p0, p1); + } + patternDuplicateEdges.savePly( + std::filesystem::path(duplicateEdgesPath).append("duplicateEdges.ply")); + } + statistics.numberOfDuplicateEdges = duplicateEdges.size(); + + // Create the set of all possible edges without coincide and duplicate edges + std::vector validEdges; + for (size_t edgeIndex = 0; edgeIndex < allPossibleEdges.size(); edgeIndex++) { + if (coincideEdges.count(edgeIndex) == 0 && + duplicateEdges.count(edgeIndex) == 0) { + validEdges.push_back(allPossibleEdges[edgeIndex]); + } + } + + return validEdges; +} + +void TopologyEnumerator::computeValidPatterns( + const std::vector &numberOfNodesPerSlot, + const size_t &numberOfDesiredEdges, + const std::filesystem::path &resultsPath, + const std::vector &allVertices, + const std::unordered_map> + &intersectingEdges, + const std::vector &validEdges, PatternSet &patternsSet) { + assert(numberOfNodesPerSlot.size() == 7); + + // Iterate over all patterns which have numberOfDesiredEdges edges from + // from the validEdges Identify patterns that contain dangling edges + const bool enoughValidEdgesExist = validEdges.size() >= numberOfDesiredEdges; + if (!enoughValidEdgesExist) { + std::filesystem::remove_all(resultsPath); // delete previous results folder + return; + } + assert(enoughValidEdgesExist); + + // Create pattern result paths + auto validPatternsPath = std::filesystem::path(resultsPath).append("Valid"); + std::filesystem::create_directory(validPatternsPath); + + const size_t numberOfPatterns = FlatPatternGeometry::binomialCoefficient( + validEdges.size(), numberOfDesiredEdges); + statistics.numberOfPatterns = numberOfPatterns; + + // Initialize pattern binary representation + std::string patternBinaryRepresentation; + patternBinaryRepresentation = std::string(numberOfDesiredEdges, '1'); + patternBinaryRepresentation += + std::string(validEdges.size() - numberOfDesiredEdges, '0'); + std::sort(patternBinaryRepresentation.begin(), + patternBinaryRepresentation.end()); + size_t patternIndex = 0; + do { + patternIndex++; + const std::string patternName = std::to_string(patternIndex); + // std::cout << "Pattern name:" + patternBinaryRepresentation << + // std::endl; isValidPattern(patternBinaryRepresentation, validEdges, + // numberOfDesiredEdges); + // Create the geometry of the pattern + // Compute the pattern edges from the binary representation + std::vector patternEdges(numberOfDesiredEdges); + size_t patternEdgeIndex = 0; + for (size_t validEdgeIndex = 0; + validEdgeIndex < patternBinaryRepresentation.size(); + validEdgeIndex++) { + if (patternBinaryRepresentation[validEdgeIndex] == '1') { + assert(patternEdgeIndex < numberOfDesiredEdges); + patternEdges[patternEdgeIndex++] = validEdges[validEdgeIndex]; + } + } + Pattern pattern; + pattern.edges = patternEdges; + + FlatPatternGeometry patternGeometry; + patternGeometry.add(allVertices, patternEdges); + + // Check if pattern contains intersecting edges + const bool patternContainsIntersectingEdges = + patternGeometry.hasIntersectingEdges(patternBinaryRepresentation, + intersectingEdges); + // Export the tiled ply file if it contains intersecting edges + if (patternContainsIntersectingEdges) { + // create the tiled geometry of the pattern + statistics.numberOfPatternsWithIntersectingEdges++; + if (debugIsOn) { + if (savePlyFiles) { + FlatPatternGeometry tiledPatternGeometry = + FlatPatternGeometry::createTile(patternGeometry); + auto intersectingPatternsPath = + std::filesystem::path(resultsPath).append("Intersecting"); + std::filesystem::create_directory(intersectingPatternsPath); + patternGeometry.savePly( + std::filesystem::path(intersectingPatternsPath) + .append(patternName) + .string() + + ".ply"); + tiledPatternGeometry.savePly( + std::filesystem::path(intersectingPatternsPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + pattern.labels.push_back(PatternLabel::IntersectingEdges); + } else { + continue; // should be uncommented in order to improve performance + } + } + + // Compute the tiled valence + const bool tiledPatternHasDanglingEdges = patternGeometry.hasDanglingEdges( + numberOfNodesPerSlot); // marks the nodes with valence>=1 + // Create the tiled geometry of the pattern + const bool hasFloatingComponents = + !patternGeometry.isFullyConnectedWhenTiled(); + FlatPatternTopology topology(numberOfNodesPerSlot, patternEdges); + const bool hasArticulationPoints = topology.containsArticulationPoints(); + FlatPatternGeometry tiledPatternGeometry = + FlatPatternGeometry::createTile( + patternGeometry); // the marked nodes of hasDanglingEdges are + // duplicated here + // Check dangling edges with vcg method + // const bool vcg_tiledPatternHasDangling = + // tiledPatternGeometry.hasUntiledDanglingEdges(); + if (tiledPatternHasDanglingEdges /*&& !hasFloatingComponents && + !hasArticulationPoints*/) { + statistics.numberOfPatternsWithADanglingEdgeOrNode++; + if (debugIsOn) { + if (savePlyFiles) { + auto danglingEdgesPath = + std::filesystem::path(resultsPath).append("Dangling"); + std::filesystem::create_directory(danglingEdgesPath); + patternGeometry.savePly(std::filesystem::path(danglingEdgesPath) + .append(patternName) + .string() + + ".ply"); + tiledPatternGeometry.savePly(std::filesystem::path(danglingEdgesPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + pattern.labels.push_back(PatternLabel::DanglingEdge); + } else { + continue; + } + } + + if (hasFloatingComponents /*&& !hasArticulationPoints && + !tiledPatternHasDanglingEdges*/) { + statistics.numberOfPatternsWithMoreThanASingleCC++; + if (debugIsOn) { + if (savePlyFiles) { + auto moreThanOneCCPath = + std::filesystem::path(resultsPath).append("MoreThanOneCC"); + std::filesystem::create_directory(moreThanOneCCPath); + patternGeometry.savePly(std::filesystem::path(moreThanOneCCPath) + .append(patternName) + .string() + + ".ply"); + tiledPatternGeometry.savePly(std::filesystem::path(moreThanOneCCPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + pattern.labels.push_back(PatternLabel::MultipleCC); + } else { + continue; + } + } + + if (hasArticulationPoints /*&& !hasFloatingComponents && + !tiledPatternHasDanglingEdges*/) { + statistics.numberOfPatternsWithArticulationPoints++; + if (exportArticulationPointsPatterns || debugIsOn) { + if (savePlyFiles) { + auto articulationPointsPath = + std::filesystem::path(resultsPath).append("ArticulationPoints"); + std::filesystem::create_directory(articulationPointsPath); + patternGeometry.savePly(std::filesystem::path(articulationPointsPath) + .append(patternName) + .string() + + ".ply"); + tiledPatternGeometry.savePly( + std::filesystem::path(articulationPointsPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + + // std::cout << "Pattern:" << patternName << std::endl; + } + pattern.labels.push_back(PatternLabel::ArticulationPoints); + } else { + continue; + } + } + + const bool isValidPattern = + !patternContainsIntersectingEdges && !tiledPatternHasDanglingEdges && + !hasFloatingComponents && !hasArticulationPoints; + if (isValidPattern) { + statistics.numberOfValidPatterns++; + if (savePlyFiles) { + // if (numberOfDesiredEdges == 4) { + // std::cout << "Saving:" + // << std::filesystem::path(validPatternsPath) + // .append(patternName) + // .string() + + // ".ply" + // << std::endl; + // } + patternGeometry.savePly(std::filesystem::path(validPatternsPath) + .append(patternName) + .string() + + ".ply"); + tiledPatternGeometry.savePly(std::filesystem::path(validPatternsPath) + .append(patternName + "_tiled") + .string() + + ".ply"); + } + pattern.labels.push_back(PatternLabel::Valid); + } + + assert(!pattern.labels.empty()); + patternsSet.patterns.push_back(pattern); + // assert(vcg_tiledPatternHasDangling == tiledPatternHasDanglingEdges); + } while (std::next_permutation(patternBinaryRepresentation.begin(), + patternBinaryRepresentation.end())); +} + +std::vector TopologyEnumerator::getCoincideEdges( + const std::vector &edgeNodeIndices) const { + std::vector coincideEdges; + if (edgeNodeIndices.size() < 3) + return coincideEdges; + for (size_t edgeNodeIndex = 0; edgeNodeIndex < edgeNodeIndices.size() - 2; + edgeNodeIndex++) { + const size_t &firstNodeIndex = edgeNodeIndices[edgeNodeIndex]; + for (size_t secondEdgeNodeIndex = edgeNodeIndex + 2; + secondEdgeNodeIndex < edgeNodeIndices.size(); secondEdgeNodeIndex++) { + const size_t &secondNodeIndex = edgeNodeIndices[secondEdgeNodeIndex]; + coincideEdges.push_back(getEdgeIndex(firstNodeIndex, secondNodeIndex)); + } + } + return coincideEdges; +} + +bool TopologyEnumerator::isValidPattern( + const std::string &patternBinaryRepresentation, + const std::vector &validEdges, + const size_t &numberOfDesiredEdges) const { + return true; +} diff --git a/topologyenumerator.hpp b/topologyenumerator.hpp new file mode 100644 index 0000000..47d8cff --- /dev/null +++ b/topologyenumerator.hpp @@ -0,0 +1,180 @@ +#ifndef TOPOLOGYENUMERATOR_HPP +#define TOPOLOGYENUMERATOR_HPP +#include "patternIO.hpp" +#include "trianglepatterngeometry.hpp" +#include "trianglepattterntopology.hpp" +#include +#include +#include +//#include +#include +#include + +using GraphEdge = std::pair; +/* + * Represents a graph formed by slots on a triangle that can either be filled or + * not. The slots are three and are: 1) one slot per vertex. Each slot can hold + * a signle node. 2) one slot per edge. Each slot can hold many nodes. 3) one + * slot in the inside of the triangle. Each slot can hold multiple nodes. + * */ +class TopologyEnumerator { + /* + * Holds the in a CCW order the vertex and the edge slots and then the face + * slot. [0],[1],[2] can either be 0 or 1 [3],[4],[5] are integers in [0,n] + * [6] an integer [0,n] + * */ + + /* + reduced nodes per slot + 0 + /\ + / \ + 2 2 + / 4 \ + / \ + 1----3-----1 + + nodes per slot + 0 + /\ + / \ + 3 5 + / 6 \ + / \ + 1----4-----2 + + slot=0 -> vi=0 + slot=1 -> vi=1 + slot=2 -> vi=2 + ... + see TrianglePatternGeometry::constructVertexVector + + */ + struct TopologicalStatistics { + size_t numberOfPossibleEdges{0}; + size_t numberOfCoincideEdges{0}; + size_t numberOfDuplicateEdges{0}; + size_t numberOfValidEdges{0}; + size_t numberOfIntersectingEdgePairs{0}; + size_t numberOfPatterns{0}; + // size_t numberOfIntersectingEdgesOverAllPatterns{0}; + size_t numberOfPatternsWithIntersectingEdges{0}; + size_t numberOfPatternsWithMoreThanASingleCC{0}; + size_t numberOfPatternsWithADanglingEdgeOrNode{0}; + size_t numberOfPatternsWithArticulationPoints{0}; + size_t numberOfValidPatterns{0}; + // nlohmann::json convertToJson() const { + // nlohmann::json json; + // json["numPossibleEdges"] = numberOfPossibleEdges; + // json["numCoincideEdges"] = numberOfCoincideEdges; + // json["numDuplicateEdges"] = numberOfDuplicateEdges; + // json["numValidEdges"] = numberOfValidEdges; + // json["numIntersectingEdgePairs"] = numberOfIntersectingEdgePairs; + // json["numPatterns"] = numberOfPatterns; + // // json["numIntersectingEdgesOverAllPatterns"] = + // // numberOfIntersectingEdgesOverAllPatterns; + // json["numPatternsWithIntersectingEdges"] = + // numberOfPatternsWithIntersectingEdges; + // json["numPatternsWithNotASingleCC"] = + // numberOfPatternsWithMoreThanASingleCC; + // json["numPatternsWithDangling"] = + // numberOfPatternsWithADanglingEdgeOrNode; + // json["numPatternsWithArticulationPoints"] = + // numberOfPatternsWithArticulationPoints; + // json["numValidPatterns"] = numberOfValidPatterns; + + // return json; + // } + void print(const std::string &setupString, + const std::filesystem::path &directoryPath) const { + std::cout << "The setup " << setupString << std::endl; + std::cout << "Has following statistics:" << std::endl; + std::cout << numberOfPossibleEdges + << " possible edges with the given arrangement of nodes" + << std::endl; + std::cout << numberOfCoincideEdges << " coincide edges" << std::endl; + std::cout << numberOfDuplicateEdges << " duplicate edges" << std::endl; + std::cout << numberOfValidEdges << " valid edges" << std::endl; + std::cout << numberOfIntersectingEdgePairs << "intersecting edge pairs " + << std::endl; + std::cout << numberOfPatterns + << " patterns can be generated with the given setup" + << std::endl; + // std::cout << numberOfIntersectingEdgesOverAllPatterns + // << " intersecting edges found over all patterns" << + // std::endl; + std::cout << numberOfPatternsWithIntersectingEdges + << " patterns found with intersecting edges" << std::endl; + std::cout << numberOfPatternsWithMoreThanASingleCC + << " 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 + << " patterns found with an articulation point" << std::endl; + std::cout << numberOfValidPatterns << " valid patterns were found" + << std::endl; + // if (!directoryPath.empty()) { + // auto json = convertToJson(); + + // std::ofstream file; + // file.open(std::filesystem::path(directoryPath) + // .append("statistics.csv") + // .string()); + // file << "setup," << setupString << "\n"; + // for (const auto &el : json.items()) { + // file << el.key() << "," << el.value() << "\n"; + // } + // } + } + }; + + TopologicalStatistics statistics; + std::vector numberOfNodesPerSlot; + size_t numberOfNodes; + + size_t + computeCoincideEdges(const std::vector &numberOfNodesPerSlot) const; + size_t computeNumberOfPossibleEdges( + const std::vector &numberOfNodesPerSlot) const; + + void createLabelMesh(const std::vector vertices, + const std::filesystem::path &savePath) const; + size_t getEdgeIndex(size_t ni0, size_t ni1) const; + +public: + TopologyEnumerator(); + + void + computeValidPatterns(const std::vector &reducedNumberOfNodesPerSlot); + +private: + std::vector + getCoincideEdges(const std::vector &edgeNodeIndices) const; + bool isValidPattern(const std::string &patternIndex, + const std::vector &validEdges, + const size_t &numberOfDesiredEdges) const; + std::vector + getValidEdges(const std::vector &numberOfNodesPerSlot, + const std::filesystem::__cxx11::path &resultsPath, + const FlatPatternGeometry &patternGeometryAllEdges, + const std::vector &allPossibleEdges); + std::unordered_set computeDuplicateEdges(); + std::unordered_set + computeCoincideEdges(const std::vector &numberOfNodesPerSlot); + void computeEdgeNodes(const std::vector &numberOfNodesPerSlot, + std::vector &nodesEdge0, + std::vector &nodesEdge1, + std::vector &nodesEdge2); + std::unordered_set + computeDuplicateEdges(const std::vector &numberOfNodesPerSlot); + void computeValidPatterns( + const std::vector &numberOfNodesPerSlot, + const size_t &numberOfDesiredEdges, + const std::filesystem::path &resultsPath, + const std::vector &vertices, + const std::unordered_map> + &intersectingEdges, + const std::vector &validEdges, PatternSet &patternsSet); +}; +#endif // TOPOLOGYENUMERATOR_HPP diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp new file mode 100644 index 0000000..1c58b29 --- /dev/null +++ b/trianglepatterngeometry.cpp @@ -0,0 +1,560 @@ +#include "trianglepatterngeometry.hpp" +#include "trianglepattterntopology.hpp" +#include +#include +#include +#include +#include +#include +#include + +size_t FlatPatternGeometry::computeTiledValence( + const size_t &nodeIndex, + const std::vector &numberOfNodesPerSlot) const { + std::vector connectedEdges; + vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges); + const size_t nodeValence = connectedEdges.size(); + assert(nodeSlot.count(nodeIndex) != 0); + const size_t nodeSlotIndex = nodeSlot.at(nodeIndex); + if (nodeSlotIndex == 0) { + return nodeValence * fanSize; + } else if (nodeSlotIndex == 1 || nodeSlotIndex == 2) { + size_t correspondingNodeIndex; + if (nodeSlotIndex == 1) { + correspondingNodeIndex = nodeIndex + 1; + } else { + correspondingNodeIndex = nodeIndex - 1; + } + std::vector + connectedEdgesCorrespondingNode; + vcg::edge::VEStarVE(&vert[correspondingNodeIndex], + connectedEdgesCorrespondingNode); + size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size(); + return fanSize / 2 * nodeValence + fanSize / 2 * correspondingNodeValence; + } else if (nodeSlotIndex == 3 || nodeSlotIndex == 5) { + size_t correspondingNodeIndex; + size_t numberOfNodesBefore; + size_t numberOfNodesAfter; + if (nodeSlotIndex == 3) { + numberOfNodesBefore = + nodeIndex - std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, 0); + correspondingNodeIndex = + + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - numberOfNodesBefore; + } else { + numberOfNodesAfter = + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - nodeIndex; + correspondingNodeIndex = + numberOfNodesAfter + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, + 0); + } + assert(correspondingNodeIndex < vn); + std::vector + connectedEdgesCorrespondingNode; + vcg::edge::VEStarVE(&vert[correspondingNodeIndex], + connectedEdgesCorrespondingNode); + size_t correspondingNodeValence = connectedEdgesCorrespondingNode.size(); + return nodeValence + correspondingNodeValence; + } else if (nodeSlotIndex == 4) { + return 2 * nodeValence; + } else if (nodeSlotIndex == 6) { + return nodeValence; + } else { + std::cerr << "Error in slot index computation" << std::endl; + } + assert(false); + return 0; +} + +size_t FlatPatternGeometry::getFanSize() const { return fanSize; } + +double FlatPatternGeometry::getTriangleEdgeSize() const { + return triangleEdgeSize; +} + +FlatPatternGeometry::FlatPatternGeometry() {} + +std::vector FlatPatternGeometry::getVertices() const {} + +FlatPatternGeometry +FlatPatternGeometry::createTile(FlatPatternGeometry &pattern) { + + const size_t fanSize = FlatPatternGeometry().getFanSize(); + FlatPatternGeometry fan(createFan(pattern)); + FlatPatternGeometry tile(fan); + + if (fanSize % 2 == 1) { + vcg::Matrix44d R; + auto rotationAxis = vcg::Point3d(0, 0, 1); + R.SetRotateDeg(180, rotationAxis); + vcg::tri::UpdatePosition::Matrix(fan, R); + } + vcg::Matrix44d T; + const double centerAngle = 2 * M_PI / fanSize; + const double triangleHeight = std::sin((M_PI - centerAngle) / 2) * + FlatPatternGeometry().triangleEdgeSize; + T.SetTranslate(0, -2 * triangleHeight, 0); + vcg::tri::UpdatePosition::Matrix(fan, T); + + FlatPatternGeometry fanOfFan = createFan(fan); + vcg::tri::Append::Mesh(tile, + fanOfFan); + vcg::tri::Clean::MergeCloseVertex(tile, 0.0000005); + vcg::tri::Allocator::CompactEveryVector(tile); + vcg::tri::UpdateTopology::VertexEdge(tile); + vcg::tri::UpdateTopology::EdgeEdge(tile); + + for (size_t vi = 0; vi < pattern.vn; vi++) { + tile.vert[vi].C() = vcg::Color4b::Blue; + } + + return tile; +} + +FlatPatternGeometry +FlatPatternGeometry::createFan(FlatPatternGeometry &pattern) { + const size_t fanSize = FlatPatternGeometry().getFanSize(); + FlatPatternGeometry fan(pattern); + FlatPatternGeometry rotatedPattern(pattern); + for (int rotationCounter = 1; rotationCounter < fanSize; rotationCounter++) { + vcg::Matrix44d R; + auto rotationAxis = vcg::Point3d(0, 0, 1); + R.SetRotateDeg(360 / fanSize, rotationAxis); + vcg::tri::UpdatePosition::Matrix(rotatedPattern, R); + vcg::tri::Append::Mesh( + fan, rotatedPattern); + } + return fan; +} + +FlatPatternGeometry::FlatPatternGeometry(FlatPatternGeometry &other) { + vcg::tri::Append::MeshCopy(*this, + other); + this->vertices = other.getVertices(); + vcg::tri::UpdateTopology::VertexEdge(*this); + vcg::tri::UpdateTopology::EdgeEdge(*this); +} + +bool FlatPatternGeometry::savePly(const std::string plyFilename) { + + int returnValue = vcg::tri::io::ExporterPLY::Save( + *this, plyFilename.c_str(), + vcg::tri::io::Mask::IOM_EDGEINDEX | vcg::tri::io::Mask::IOM_VERTCOLOR, + false); + if (returnValue != 0) { + std::cerr << vcg::tri::io::ExporterPLY::ErrorMsg( + returnValue) + << std::endl; + } + + return static_cast(returnValue); +} + +void FlatPatternGeometry::add(const std::vector &vertices) { + this->vertices = vertices; + std::for_each(vertices.begin(), vertices.end(), [&](const vcg::Point3d &p) { + vcg::tri::Allocator::AddVertex(*this, p); + }); + vcg::tri::UpdateTopology::VertexEdge(*this); + vcg::tri::UpdateTopology::EdgeEdge(*this); +} + +void FlatPatternGeometry::add(const std::vector &edges) { + std::for_each(edges.begin(), edges.end(), [&](const vcg::Point2i &e) { + vcg::tri::Allocator::AddEdge(*this, e[0], e[1]); + }); + vcg::tri::UpdateTopology::VertexEdge(*this); + vcg::tri::UpdateTopology::EdgeEdge(*this); +} + +void FlatPatternGeometry::add(const std::vector &vertices, + const std::vector &edges) { + add(vertices); + add(edges); +} + +void FlatPatternGeometry::add(const std::vector &numberOfNodesPerSlot, + const std::vector &edges) { + assert(numberOfNodesPerSlot.size() == 7); + auto vertices = + constructVertexVector(numberOfNodesPerSlot, fanSize, triangleEdgeSize); + add(vertices, edges); + vcg::tri::UpdateTopology::VertexEdge(*this); + vcg::tri::UpdateTopology::EdgeEdge(*this); + updateEigenEdgeAndVertices(); +} + +std::vector FlatPatternGeometry::constructVertexVector( + const std::vector &numberOfNodesPerSlot, const size_t &fanSize, + const double &triangleEdgeSize) { + + std::vector vertices; + const double centerAngle = 2 * M_PI / fanSize; + const double triangleHeight = + std::sin((M_PI - centerAngle) / 2) * triangleEdgeSize; + const double baseEdgeSize = + 2 * triangleEdgeSize * std::cos((M_PI - centerAngle) / 2); + const vcg::Point3d triangleV0 = vcg::Point3d(0, 0, 0); + const vcg::Point3d triangleV1 = + vcg::Point3d(-baseEdgeSize / 2, -triangleHeight, 0); + const vcg::Point3d triangleV2 = + vcg::Point3d(baseEdgeSize / 2, -triangleHeight, 0); + + // Nodes + if (numberOfNodesPerSlot[0] == 1) { + // vertices[0] = triangleV0; + vertices.push_back(triangleV0); + } + if (numberOfNodesPerSlot[1] == 1) { + // vertices[1] = triangleV1; + vertices.push_back(triangleV1); + } + if (numberOfNodesPerSlot[2] == 1) { + // vertices[2] = triangleV2; + vertices.push_back(triangleV2); + } + + // Edges + if (numberOfNodesPerSlot[3] != 0) { + const double offset0 = 1.0 / (numberOfNodesPerSlot[3] + 1); + const vcg::Point3d edgeVector0(triangleV1 - triangleV0); + for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[3]; + vertexIndex++) { + + // vertices[std::accumulate(numberOfNodesPerSlot.begin(), + // numberOfNodesPerSlot.begin() + 2, 0) + + // vertexIndex] = + vertices.push_back(triangleV0 + + edgeVector0.operator*((vertexIndex + 1) * offset0)); + } + } + + if (numberOfNodesPerSlot[4] != 0) { + const double offset1 = 1.0 / (numberOfNodesPerSlot[4] + 1); + const vcg::Point3d edgeVector1(triangleV2 - triangleV1); + for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[4]; + vertexIndex++) { + + // vertices[std::accumulate(numberOfNodesPerSlot.begin(), + // numberOfNodesPerSlot.begin() + 3, 0) + + // vertexIndex] = + vertices.push_back(triangleV1 + + edgeVector1.operator*((vertexIndex + 1) * offset1)); + } + } + + if (numberOfNodesPerSlot[5] != 0) { + const double offset2 = 1.0 / (numberOfNodesPerSlot[5] + 1); + const vcg::Point3d edgeVector2(triangleV0 - triangleV2); + for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[5]; + vertexIndex++) { + // vertices[std::accumulate(numberOfNodesPerSlot.begin(), + // numberOfNodesPerSlot.begin() + 4, 0) + + // vertexIndex] = + vertices.push_back(triangleV2 + + edgeVector2.operator*((vertexIndex + 1) * offset2)); + } + } + + // Face + if (numberOfNodesPerSlot[6] != 0) { + const vcg::Point3d triangleCenter((triangleV0 + triangleV1 + triangleV2) / + 3); + const double radius = 0.1; + const double offsetRad = 2 * M_PI / numberOfNodesPerSlot[6]; + for (size_t vertexIndex = 0; vertexIndex < numberOfNodesPerSlot[6]; + vertexIndex++) { + const double pX = + triangleCenter[0] + radius * std::cos(vertexIndex * offsetRad); + const double pY = + triangleCenter[1] + radius * std::sin(vertexIndex * offsetRad); + /*vertices[std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 5, 0) + + vertexIndex] =*/ + vertices.push_back(vcg::Point3d(pX, pY, 0)); + } + } + return vertices; +} + +void FlatPatternGeometry::constructNodeToTiledValenceMap( + const std::vector &numberOfNodesPerSlot) { + for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + const size_t tiledValence = + computeTiledValence(nodeIndex, numberOfNodesPerSlot); + nodeTiledValence[nodeIndex] = tiledValence; + } +} + +bool FlatPatternGeometry::hasDanglingEdges( + const std::vector &numberOfNodesPerSlot) { + if (nodeSlot.empty()) { + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeSlot); + } + if (correspondingNode.empty()) { + constructCorresponginNodeMap(numberOfNodesPerSlot); + } + if (nodeTiledValence.empty()) { + constructNodeToTiledValenceMap(numberOfNodesPerSlot); + } + bool hasDanglingEdges = false; + for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + const size_t tiledValence = nodeTiledValence[nodeIndex]; + if (tiledValence == 1) { + vert[nodeIndex].C() = vcg::Color4b::Red; + hasDanglingEdges = true; + } + } + return hasDanglingEdges; +} + +bool FlatPatternGeometry::hasUntiledDanglingEdges() { + // vcg::tri::Clean::MergeCloseVertex(*this, + // 0.0000005); + // vcg::tri::Allocator::CompactEveryVector(*this); + // vcg::tri::UpdateTopology::VertexEdge(*this); + // vcg::tri::UpdateTopology::EdgeEdge(*this); + bool hasDanglingEdges = false; + for (size_t vi = 0; vi < vn; vi++) { + std::vector connectedEdges; + vcg::edge::VEStarVE(&vert[vi], connectedEdges); + const size_t nodeValence = connectedEdges.size(); + if (nodeValence == 1) { + if (vert[vi].C().operator==(vcg::Color4b(vcg::Color4b::Red))) { + + } else { + vert[vi].C() = vcg::Color4b::Blue; + } + hasDanglingEdges = true; + } + } + return hasDanglingEdges; +} + +// TODO: The function expects that the numberOfNodesPerSlot follows a specific +// format and that the vertex container was populated in a particular order. +void FlatPatternGeometry::constructCorresponginNodeMap( + const std::vector &numberOfNodesPerSlot) { + assert(vn != 0 && !nodeSlot.empty() && correspondingNode.empty() && + numberOfNodesPerSlot.size() == 7); + + for (size_t nodeIndex = 0; nodeIndex < vn; nodeIndex++) { + const size_t slotIndex = nodeSlot[nodeIndex]; + if (slotIndex == 1) { + correspondingNode[nodeIndex] = nodeIndex + 1; + } else if (slotIndex == 2) { + correspondingNode[nodeIndex] = nodeIndex - 1; + } else if (slotIndex == 3) { + const size_t numberOfNodesBefore = + nodeIndex - std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, 0); + correspondingNode[nodeIndex] = + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - numberOfNodesBefore; + } else if (slotIndex == 5) { + const size_t numberOfNodesAfter = + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - nodeIndex; + correspondingNode[nodeIndex] = + numberOfNodesAfter + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, + 0); + } + } +} + +bool FlatPatternGeometry::isFullyConnectedWhenTiled() { + 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 (nodeSlot[nodeIndex] == 1 || nodeSlot[nodeIndex] == 4 || + nodeSlot[nodeIndex] == 2) { + std::vector connectedEdges; + vcg::edge::VEStarVE(&vert[nodeIndex], connectedEdges); + const size_t nodeValence = connectedEdges.size(); + if (nodeValence != 0) { + bottomInterfaceIsConnected = true; + break; + } + } + } + if (!bottomInterfaceIsConnected) { + return false; + } + + FlatPatternGeometry 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; + } + + // 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(); + // } + + // 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; +} + +bool FlatPatternGeometry::hasIntersectingEdges( + const std::string &patternBinaryRepresentation, + const std::unordered_map> + &intersectingEdges) { + bool containsIntersectingEdges = false; + for (size_t validEdgeIndex = 0; + validEdgeIndex < patternBinaryRepresentation.size(); validEdgeIndex++) { + if (patternBinaryRepresentation[validEdgeIndex] == '1' && + intersectingEdges.count(validEdgeIndex) != 0) { + for (auto edgeIndexIt = + intersectingEdges.find(validEdgeIndex)->second.begin(); + edgeIndexIt != intersectingEdges.find(validEdgeIndex)->second.end(); + edgeIndexIt++) { + if (patternBinaryRepresentation[*edgeIndexIt] == '1') { + containsIntersectingEdges = true; + // statistics.numberOfIntersectingEdgesOverAllPatterns++; + break; // should be uncommented in order to improve + // performance + } + } + if (containsIntersectingEdges) { + break; // should be uncommented in order to improve performance + } + } + } + + return containsIntersectingEdges; +} + +std::unordered_map> +FlatPatternGeometry::getIntersectingEdges( + size_t &numberOfIntersectingEdgePairs) const { + std::unordered_map> intersectingEdges; + numberOfIntersectingEdgePairs = 0; + size_t numberOfEdgePairs; + if (en == 0 || en == 1) { + numberOfEdgePairs = 0; + } else { + numberOfEdgePairs = binomialCoefficient(en, 2); + } + + for (size_t edgePairIndex = 0; edgePairIndex < numberOfEdgePairs; + edgePairIndex++) { + size_t ei0 = + en - 2 - + floor(sqrt(-8 * edgePairIndex + 4 * en * (en - 1) - 7) / 2.0 - 0.5); + size_t ei1 = edgePairIndex + ei0 + 1 - en * (en - 1) / 2 + + (en - ei0) * ((en - ei0) - 1) / 2; + const vcg::Point2d p0(edge[ei0].cP(0)[0], edge[ei0].cP(0)[1]); + const float epsilon = 0.000005; + vcg::Box2d bb0(p0 - vcg::Point2d(epsilon, epsilon), + p0 + vcg::Point2d(epsilon, epsilon)); + const vcg::Point2d p1(edge[ei0].cP(1)[0], edge[ei0].cP(1)[1]); + vcg::Box2d bb1(p1 - vcg::Point2d(epsilon, epsilon), + p1 + vcg::Point2d(epsilon, epsilon)); + const vcg::Point2d p2(edge[ei1].cP(0)[0], edge[ei1].cP(0)[1]); + vcg::Box2d bb2(p2 - vcg::Point2d(epsilon, epsilon), + p2 + vcg::Point2d(epsilon, epsilon)); + if (bb2.Collide(bb1) || bb2.Collide(bb0)) + continue; + const vcg::Point2d p3(edge[ei1].cP(1)[0], edge[ei1].cP(1)[1]); + vcg::Box2d bb3(p3 - vcg::Point2d(epsilon, epsilon), + p3 + vcg::Point2d(epsilon, epsilon)); + if (bb3.Collide(bb1) || bb3.Collide(bb0)) + continue; + const vcg::Segment2d s0(p0, p1); + const vcg::Segment2d s1(p2, p3); + + vcg::Point2d intersectionPoint; + const bool edgesIntersect = + vcg::SegmentSegmentIntersection(s0, s1, intersectionPoint); + if (edgesIntersect) { + numberOfIntersectingEdgePairs++; + intersectingEdges[ei0].insert(ei1); + intersectingEdges[ei1].insert(ei0); + } + } + return intersectingEdges; +} diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp new file mode 100644 index 0000000..73da326 --- /dev/null +++ b/trianglepatterngeometry.hpp @@ -0,0 +1,68 @@ +#ifndef FLATPATTERNGEOMETRY_HPP +#define FLATPATTERNGEOMETRY_HPP +#include "edgemesh.hpp" +#include +#include +#include + +class FlatPatternGeometry : public VCGEdgeMesh { +private: + size_t + computeTiledValence(const size_t &nodeIndex, + const std::vector &numberOfNodesPerSlot) const; + size_t getNodeValence(const size_t &nodeIndex) const; + size_t getNodeSlot(const size_t &nodeIndex) const; + + const size_t fanSize{6}; + std::vector vertices; + const double triangleEdgeSize{1}; // radius edge + std::unordered_map nodeSlot; + std::unordered_map correspondingNode; + std::unordered_map nodeTiledValence; + + void + constructCorresponginNodeMap(const std::vector &numberOfNodesPerSlot); + + void constructNodeToTiledValenceMap( + const std::vector &numberOfNodesPerSlot); + +public: + FlatPatternGeometry(); + /*The following function should be a copy constructor with + * a const argument but this is impossible due to the + * vcglib interface. + * */ + FlatPatternGeometry(FlatPatternGeometry &other); + bool savePly(const std::string plyFilename); + 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 &numberOfNodesPerSlot, + const std::vector &edges); + static std::vector + constructVertexVector(const std::vector &numberOfNodesPerSlot, + const size_t &fanSize, const double &triangleEdgeSize); + bool hasDanglingEdges(const std::vector &numberOfNodesPerSlot); + std::vector getVertices() const; + static FlatPatternGeometry createFan(FlatPatternGeometry &pattern); + static FlatPatternGeometry createTile(FlatPatternGeometry &pattern); + double getTriangleEdgeSize() const; + bool hasUntiledDanglingEdges(); + std::unordered_map> + getIntersectingEdges(size_t &numberOfIntersectingEdgePairs) const; + + static size_t binomialCoefficient(size_t n, size_t m) { + assert(n >= m); + return tgamma(n + 1) / (tgamma(m + 1) * tgamma(n - m + 1)); + } + bool isFullyConnectedWhenTiled(); + + bool hasIntersectingEdges( + const std::string &patternBinaryRepresentation, + const std::unordered_map> + &intersectingEdges); + bool isPatternValid(); + size_t getFanSize() const; +}; +#endif // FLATPATTERNGEOMETRY_HPP diff --git a/trianglepattterntopology.cpp b/trianglepattterntopology.cpp new file mode 100644 index 0000000..9686396 --- /dev/null +++ b/trianglepattterntopology.cpp @@ -0,0 +1,298 @@ +#include "trianglepattterntopology.hpp" +#include +#include +#include + +FlatPatternTopology::FlatPatternTopology() {} + +FlatPatternTopology::FlatPatternTopology( + const std::vector &numberOfNodesPerSlot, + const std::vector &edges) + : numberOfNodesPerSlot(numberOfNodesPerSlot) { + pattern = BoostGraph(std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.end(), 0)); + for (const vcg::Point2i &e : edges) { + boost::add_edge(e[0], e[1], pattern); + } + + constructNodeToSlotMap(); + constructSlotToNodeMap(); + constructCorresponginNodeMap(); +} + +bool FlatPatternTopology::containsArticulationPoints() const { + assert(numberOfNodesPerSlot.size() == 7 && numberOfNodesPerSlot[4] < 2 && + !nodeToSlot.empty() && !correspondingNode.empty()); + BoostGraph copyOfPattern(pattern); + // std::cout << std::endl; + std::vector componentsBefore(boost::num_vertices(copyOfPattern)); + size_t num_components = + boost::connected_components(copyOfPattern, &componentsBefore[0]); + // std::cout << "Number of cc before:" << num_components << std::endl; + // printGraph(copyOfPattern); + + copyOfPattern = constructRotationallySymmetricPattern( + copyOfPattern, slotToNode, nodeToSlot, correspondingNode); + // // Remove edges connected to the bottom edge node + // assert(slotToNode.find(4) != slotToNode.end()); + // std::unordered_set bottomEdgeNodeSet = + // (slotToNode.find(4))->second; + // size_t bottomEdgeNodeIndex = *bottomEdgeNodeSet.begin(); + // boost::clear_vertex(bottomEdgeNodeIndex, copyOfPattern); + + std::vector componentsAfter(boost::num_vertices(copyOfPattern)); + num_components = + boost::connected_components(copyOfPattern, &componentsAfter[0]); + // std::cout << "Number of cc after:" << num_components << std::endl; + // printGraph(copyOfPattern); + + // Compute articulation points on the edited graph + std::vector articulationPoints; + boost::articulation_points(copyOfPattern, + std::back_inserter(articulationPoints)); + // std::cout << "Found " << articulationPoints.size() + // << " articulation points.\n"; + // size_t numberOfNonValidArticulationPoints = 0; + for (std::size_t i = 0; i < articulationPoints.size(); ++i) { + // std::cout << articulationPoints[i] << std::endl; + // if (boost::out_degree(articulationPoints[i], copyOfPattern) < 3) { + // numberOfNonValidArticulationPoints++; + // } + } + // if (numberOfNonValidArticulationPoints == articulationPoints.size()) { + // return false; + // } + return !articulationPoints.empty(); +} + +void FlatPatternTopology::constructNodeToSlotMap( + const std::vector &numberOfNodesPerSlot, + std::unordered_map &nodeToSlot) { + const size_t numberOfNodes = std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.end(), 0); + assert(numberOfNodes != 0 && nodeToSlot.empty() && + numberOfNodesPerSlot.size() == 7); + std::vector maxNodeIndexPerSlot(numberOfNodesPerSlot.size()); + for (size_t i = 0; i < maxNodeIndexPerSlot.size(); ++i) { + maxNodeIndexPerSlot[i] = std::accumulate( + numberOfNodesPerSlot.begin(), numberOfNodesPerSlot.begin() + i + 1, 0); + } + + for (size_t nodeIndex = 0; nodeIndex < numberOfNodes; nodeIndex++) { + const size_t slotIndex = + std::distance(maxNodeIndexPerSlot.begin(), + std::upper_bound(maxNodeIndexPerSlot.begin(), + maxNodeIndexPerSlot.end(), nodeIndex)); + nodeToSlot[nodeIndex] = slotIndex; + } +} + +void FlatPatternTopology::constructNodeToSlotMap() { + constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); +} + +void FlatPatternTopology::constructSlotToNodeMap() { + constructSlotToNodeMap(nodeToSlot, slotToNode); +} + +void FlatPatternTopology::constructSlotToNodeMap( + const std::unordered_map &nodeToSlot, + std::unordered_map> &slotToNode) { + assert(!nodeToSlot.empty()); + for (size_t nodeIndex = 0; nodeIndex < nodeToSlot.size(); nodeIndex++) { + slotToNode[nodeToSlot.at(nodeIndex)].insert(nodeIndex); + } +} + +// TODO: The function expects that the numberOfNodesPerSlot follows a +// specific format and that the vertex container was populated in a +// particular order. +void FlatPatternTopology::constructCorresponginNodeMap() { + assert(!nodeToSlot.empty() && correspondingNode.empty() && + numberOfNodesPerSlot.size() == 7); + + for (size_t nodeIndex = 0; nodeIndex < boost::num_vertices(pattern); + nodeIndex++) { + const size_t slotIndex = nodeToSlot[nodeIndex]; + if (slotIndex == 1) { + correspondingNode[nodeIndex] = nodeIndex + 1; + } else if (slotIndex == 2) { + correspondingNode[nodeIndex] = nodeIndex - 1; + } else if (slotIndex == 3) { + const size_t numberOfNodesBefore = + nodeIndex - std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, 0); + correspondingNode[nodeIndex] = + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - numberOfNodesBefore; + } else if (slotIndex == 5) { + const size_t numberOfNodesAfter = + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 6, 0) - + 1 - nodeIndex; + correspondingNode[nodeIndex] = + numberOfNodesAfter + std::accumulate(numberOfNodesPerSlot.begin(), + numberOfNodesPerSlot.begin() + 3, + 0); + } + } +} + +/* + * In this function I create an "extended" pattern of the one in the base + * triangle by: + * 1)copying the edges from left edge to the right edge and vice versa + * 2)I label all nodes that lay on the boarder of the base triangle as + * "external" and add edges connecting them to each other (this is wrong in the + * case in which all "external" nodes are connected via a single node!)) + * */ + +BoostGraph FlatPatternTopology::constructRotationallySymmetricPattern( + const BoostGraph &pattern, + const std::unordered_map> &slotToNodes, + const std::unordered_map &nodeToSlot, + const std::unordered_map &correspondingNode) { + BoostGraph rotationallySymmetricPattern(pattern); + + boost::graph_traits::out_edge_iterator ei, ei_end; + // Copy edges that lay on the left edge to the right edge + const auto slot3NodesPairIt = slotToNodes.find(3); + if (slot3NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot3NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); + ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot == 3 || vtSlot == 1 || vtSlot == 0) { + // Connect the corresponding nodes on the opposite edge + auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); + assert(correspondingNodeIndexIt != correspondingNode.end()); + auto correspondingVtIt = correspondingNode.find(vt); + assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); + const size_t &correspondingNodeIndex = + correspondingNodeIndexIt->second; + size_t correspondingVt = 0; + if (correspondingVtIt != correspondingNode.end()) { + correspondingVt = correspondingVtIt->second; + } + boost::add_edge(correspondingNodeIndex, correspondingVt, + rotationallySymmetricPattern); + } + } + } + } + + // Copy edges that lay on the right edge to the left edge + const auto slot5NodesPairIt = slotToNodes.find(5); + if (slot5NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot5NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); + ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot == 5 || vtSlot == 2 || vtSlot == 0) { + // Connect the corresponding nodes on the opposite edge + auto correspondingNodeIndexIt = correspondingNode.find(nodeIndex); + assert(correspondingNodeIndexIt != correspondingNode.end()); + auto correspondingVtIt = correspondingNode.find(vt); + assert(correspondingVtIt != correspondingNode.end() || vtSlot == 0); + const size_t &correspondingNodeIndex = + correspondingNodeIndexIt->second; + size_t correspondingVt = 0; + if (correspondingVtIt != correspondingNode.end()) { + correspondingVt = correspondingVtIt->second; + } + boost::add_edge(correspondingNodeIndex, correspondingVt, + rotationallySymmetricPattern); + } + } + } + } + + // NOTE: The problem with that approach is that I connect !all! "external" + // nodes with each other, which might not be entirely true. If the number of + // cc of the tilled configuration is not 1 this might not label patterns as + // having an articulation node although they might have one. Create set of + // nodes connected with "external" edges + std::unordered_set externallyConnectedNodes; + // Mark the star node as external + const auto slot0NodesPairIt = slotToNodes.find(0); + if (slot0NodesPairIt != slotToNodes.end()) { + externallyConnectedNodes.insert(slot0NodesPairIt->second.begin(), + slot0NodesPairIt->second.end()); + } + // Mark all bottom nodes as external since they are allways connected to the + // south-neighbouring pattern + const auto slot4NodesPairIt = slotToNodes.find(4); + if (slot4NodesPairIt != slotToNodes.end()) { + externallyConnectedNodes.insert(slot4NodesPairIt->second.begin(), + slot4NodesPairIt->second.end()); + } + // Add all slot3 nodes that have a connection to the "inside" + if (slot3NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot3NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); + ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot != 3) { + auto correspondingNodePairIt = correspondingNode.find(nodeIndex); + assert(correspondingNodePairIt != correspondingNode.end()); + externallyConnectedNodes.insert(correspondingNodePairIt->second); + } + } + } + } + // Add all slot5 nodes that have a connection to the "inside" + if (slot5NodesPairIt != slotToNodes.end()) { + for (const size_t &nodeIndex : slot5NodesPairIt->second) { + for (boost::tie(ei, ei_end) = boost::out_edges(nodeIndex, pattern); + ei != ei_end; ++ei) { + auto vt = boost::target(*ei, pattern); + const auto vtNodeSlotPairIt = nodeToSlot.find(vt); + assert(vtNodeSlotPairIt != nodeToSlot.end()); + const size_t vtSlot = vtNodeSlotPairIt->second; + if (vtSlot != 5) { + auto correspondingNodePairIt = correspondingNode.find(nodeIndex); + assert(correspondingNodePairIt != correspondingNode.end()); + externallyConnectedNodes.insert(correspondingNodePairIt->second); + } + } + } + } + + // connecting all is wrong. Maybe I should check whether the external nodes + // are connected via a single node? If so this node is an articulation point. + // I could test this by checking if it filters out only the falsely labeled + // pattern 2367 + const size_t &n = externallyConnectedNodes.size(); + const size_t numberOfExternalEdges = n * (n - 1) / 2; + // Connect all external nodes with each other + for (size_t edgeIndex = 0; edgeIndex < numberOfExternalEdges; edgeIndex++) { + const int sei0 = + n - 2 - + std::floor(std::sqrt(-8 * edgeIndex + 4 * n * (n - 1) - 7) / 2.0 - 0.5); + const int sei1 = edgeIndex + sei0 + 1 - n * (n - 1) / 2 + + (n - sei0) * ((n - sei0) - 1) / 2; + const size_t ni0 = *std::next(externallyConnectedNodes.begin(), sei0); + const size_t ni1 = *std::next(externallyConnectedNodes.begin(), sei1); + boost::add_edge(ni0, ni1, rotationallySymmetricPattern); + } + + return rotationallySymmetricPattern; +} + +void FlatPatternTopology::printGraph(const BoostGraph &g) const { + boost::graph_traits::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(g); ei != ei_end; ++ei) + std::cout << (char)(boost::source(*ei, g) + 'A') << " -- " + << (char)(boost::target(*ei, g) + 'A'); + std::cout << std::endl; +} diff --git a/trianglepattterntopology.hpp b/trianglepattterntopology.hpp new file mode 100644 index 0000000..81844ea --- /dev/null +++ b/trianglepattterntopology.hpp @@ -0,0 +1,49 @@ +#ifndef FLATPATTTERNTOPOLOGY_HPP +#define FLATPATTTERNTOPOLOGY_HPP +#include +#include +#include +#include +#include + +using BoostGraph = + boost::adjacency_list; +using vertex_t = boost::graph_traits::vertex_descriptor; + +class FlatPatternTopology { + +public: + bool containsArticulationPoints() const; + FlatPatternTopology(const std::vector &numberOfNodesPerSlot, + const std::vector &edges); + static void + constructNodeToSlotMap(const std::vector &numberOfNodesPerSlot, + std::unordered_map &nodeToSlot); + static void constructSlotToNodeMap( + const std::unordered_map &nodeToSlot, + std::unordered_map> &slotToNode); + + FlatPatternTopology(); + +private: + BoostGraph pattern; + std::vector numberOfNodesPerSlot; + std::unordered_map nodeToSlot; + std::unordered_map> slotToNode; + std::unordered_map correspondingNode; + void constructCorresponginNodeMap(); + /* + * Creates a pattern which is a copy of the input pattern but with added edges + * that result + * */ + void printGraph(const BoostGraph &g) const; + static BoostGraph constructRotationallySymmetricPattern( + const BoostGraph &pattern, + const std::unordered_map> &slotToNodes, + const std::unordered_map &nodeToSlot, + const std::unordered_map &correspondingNode); + void constructNodeToSlotMap(); + void constructSlotToNodeMap(); +}; + +#endif // FLATPATTTERNTOPOLOGY_HPP diff --git a/utilities.hpp b/utilities.hpp new file mode 100644 index 0000000..73d5107 --- /dev/null +++ b/utilities.hpp @@ -0,0 +1,272 @@ +#ifndef UTILITIES_H +#define UTILITIES_H + +#include +#include +#include +#include + +namespace Utilities { +inline void parseIntegers(const std::string &str, std::vector &result) { + typedef std::regex_iterator re_iterator; + typedef re_iterator::value_type re_iterated; + + std::regex re("(\\d+)"); + + re_iterator rit(str.begin(), str.end(), re); + re_iterator rend; + + std::transform(rit, rend, std::back_inserter(result), + [](const re_iterated &it) { return std::stoi(it[1]); }); +} + +// std::string convertToLowercase(const std::string &s) { +// std::string lowercase; +// std::transform(s.begin(), s.end(), lowercase.begin(), +// [](unsigned char c) { return std::tolower(c); }); +// return lowercase; +//} +// bool hasExtension(const std::string &filename, const std::string &extension) +// { +// const std::filesystem::path path(filename); +// if (!path.has_extension()) { +// std::cerr << "Error: No file extension found in " << filename << +// std::endl; return false; +// } + +// const std::string detectedExtension = path.extension().string(); + +// if (convertToLowercase(detectedExtension) != convertToLowercase(extension)) +// { +// std::cerr << "Error: detected extension is " + detectedExtension + +// " and not " + extension +// << std::endl; +// return false; +// } +// return true; +//} + +} // namespace Utilities + +struct NodalForce { + size_t index; + size_t dof; + double magnitude; +}; +struct Vector6d : public std::array { + Vector6d() { + for (size_t i = 0; i < 6; i++) { + this->operator[](i) = 0; + } + } + + Vector6d(const double &d) { + for (size_t i = 0; i < 6; i++) { + this->operator[](i) = d; + } + } + + Vector6d(const std::initializer_list &initList) { + std::copy(initList.begin(), initList.end(), this->begin()); + } + + Vector6d operator*(const double &d) const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) * d; + } + return result; + } + + Vector6d operator*(const Vector6d &v) const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) * v[i]; + } + return result; + } + + Vector6d operator/(const double &d) const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) / d; + } + return result; + } + + Vector6d operator+(const Vector6d &v) const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) + v[i]; + } + return result; + } + + Vector6d operator-(const Vector6d &v) const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + result[i] = this->operator[](i) - v[i]; + } + return result; + } + + Vector6d inverted() const { + Vector6d result; + for (size_t i = 0; i < 6; i++) { + assert(this->operator[](i) != 0); + result[i] = 1 / this->operator[](i); + } + return result; + } + bool isZero() const { + for (size_t i = 0; i < 6; i++) { + if (this->operator[](i) != 0) + return false; + } + return true; + } + + double squaredNorm() const { + double squaredNorm = 0; + std::for_each(begin(), end(), + [&](const double &v) { squaredNorm += pow(v, 2); }); + return squaredNorm; + } + + double norm() const { return sqrt(squaredNorm()); } + + bool isFinite() const { + return std::any_of(begin(), end(), [](const double &v) { + if (!std::isfinite(v)) { + return false; + } + return true; + }); + } +}; + +// namespace ConfigurationFile { + +// inline void getPlyFilename(const std::string jsonFilepath, +// std::string &plyFilename) { +// std::ifstream inFile(jsonFilepath); +// std::string jsonContents((std::istreambuf_iterator(inFile)), +// std::istreambuf_iterator()); +// nlohmann::json jsonFile(nlohmann::json::parse(jsonContents)); + +// if (jsonFile.contains("plyFilename")) { +// plyFilename = jsonFile["plyFilename"]; +// } +//} + +// inline void +// getNodalForces(const std::string jsonFilepath, +// std::unordered_map &nodalForces) { +// std::ifstream inFile(jsonFilepath); +// std::string jsonContents((std::istreambuf_iterator(inFile)), +// std::istreambuf_iterator()); +// nlohmann::json jsonFile(nlohmann::json::parse(jsonContents)); + +// nodalForces.clear(); +// if (jsonFile.contains("forces")) { +// std::vector> forces = jsonFile["forces"]; +// for (size_t forceIndex = 0; forceIndex < forces.size(); forceIndex++) { +// const BeamFormFinder::NodalForce nf{ +// static_cast(forces[forceIndex][0]), +// static_cast(forces[forceIndex][1]), forces[forceIndex][2]}; +// const size_t vertexIndex = forces[forceIndex][0]; +// const Eigen::Vector3d forceVector( +// forces[forceIndex][1], forces[forceIndex][2], +// forces[forceIndex][3]); +// assert(forceIndex >= 0 && forceVector.norm() >= 0); +// nodalForces[vertexIndex] = forceVector; +// } +// } +//} + +// inline void getFixedVertices(const std::string jsonFilepath, +// std::vector &fixedVertices) { +// std::ifstream inFile(jsonFilepath); +// std::string jsonContents((std::istreambuf_iterator(inFile)), +// std::istreambuf_iterator()); +// nlohmann::json jsonFile(nlohmann::json::parse(jsonContents)); + +// fixedVertices.clear(); +// if (jsonFile.contains("fixedVertices")) { +// fixedVertices = +// static_cast>(jsonFile["fixedVertices"]); +// } +//} + +// struct SimulationScenario { +// std::string edgeMeshFilename; +// std::vector fixedVertices; +// std::unordered_map nodalForces; +//}; + +// void to_json(nlohmann::json &json, const SimulationScenario &scenario) { +// json["plyFilename"] = scenario.edgeMeshFilename; +// if (!scenario.fixedVertices.empty()) { +// json["fixedVertices"] = scenario.fixedVertices; +// } +// if (!scenario.nodalForces.empty()) { +// std::vector> forces; +// std::transform(scenario.nodalForces.begin(), scenario.nodalForces.end(), +// std::back_inserter(forces), +// [](const std::pair &f) { +// return std::tuple{ +// f.first, f.second[0], f.second[1], f.second[2]}; +// }); +// json["forces"] = forces; +// } +//} +//} // namespace ConfigurationFile +#include "polyscope/curve_network.h" +#include "polyscope/polyscope.h" + +inline void initPolyscope() { + if (polyscope::state::initialized) { + return; + } + polyscope::init(); + polyscope::options::groundPlaneEnabled = false; + polyscope::view::upDir = polyscope::view::UpDir::ZUp; +} + +inline void registerWorldAxes() { + if (!polyscope::state::initialized) { + initPolyscope(); + } + Eigen::MatrixX3d axesPositions(4, 3); + axesPositions.row(0) = Eigen::Vector3d(0, 0, 0); + axesPositions.row(1) = Eigen::Vector3d(1, 0, 0); + axesPositions.row(2) = Eigen::Vector3d(0, 1, 0); + axesPositions.row(3) = Eigen::Vector3d(0, 0, 1); + Eigen::MatrixX2i axesEdges(3, 2); + axesEdges.row(0) = Eigen::Vector2i(0, 1); + axesEdges.row(1) = Eigen::Vector2i(0, 2); + axesEdges.row(2) = Eigen::Vector2i(0, 3); + Eigen::MatrixX3d axesColors(3, 3); + axesColors.row(0) = Eigen::Vector3d(1, 0, 0); + axesColors.row(1) = Eigen::Vector3d(0, 1, 0); + axesColors.row(2) = Eigen::Vector3d(0, 0, 1); + + const std::string worldAxesName = "World Axes"; + polyscope::registerCurveNetwork(worldAxesName, axesPositions, axesEdges); + polyscope::getCurveNetwork(worldAxesName)->setRadius(0.001); + const std::string worldAxesColorName = worldAxesName + " Color"; + polyscope::getCurveNetwork(worldAxesName) + ->addEdgeColorQuantity(worldAxesColorName, axesColors) + ->setEnabled(true); +} + +template +void constructInverseMap(const T1 &map, T2 &oppositeMap) { + assert(!map.empty()); + oppositeMap.clear(); + for (const auto &mapIt : map) { + oppositeMap[mapIt.second] = mapIt.first; + } +} + +#endif // UTILITIES_H diff --git a/vcgtrimesh.cpp b/vcgtrimesh.cpp new file mode 100644 index 0000000..e68b312 --- /dev/null +++ b/vcgtrimesh.cpp @@ -0,0 +1,78 @@ +#include "vcgtrimesh.hpp" +#include "wrap/io_trimesh/import_obj.h" +#include "wrap/io_trimesh/import_off.h" +#include + +void VCGTriMesh::loadFromPlyFile(const std::string &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; + if (nanoply::NanoPlyWrapper::LoadModel( + std::filesystem::absolute(filename).c_str(), *this, mask) != 0) { + std::cout << "Could not load tri mesh" << std::endl; + } + vcg::tri::UpdateTopology::FaceFace(*this); + vcg::tri::UpdateTopology::VertexFace(*this); + vcg::tri::UpdateNormal::PerVertexNormalized(*this); +} + +Eigen::MatrixX3d VCGTriMesh::getVertices() const { + Eigen::MatrixX3d vertices(VN(), 3); + for (size_t vi = 0; vi < VN(); vi++) { + VCGTriMesh::CoordType vertexCoordinates = vert[vi].cP(); + vertices.row(vi) = vertexCoordinates.ToEigenVector(); + } + return vertices; +} + +Eigen::MatrixX3i VCGTriMesh::getFaces() const { + Eigen::MatrixX3i faces(FN(), 3); + for (int fi = 0; fi < FN(); fi++) { + const VCGTriMesh::FaceType &face = this->face[fi]; + const size_t v0 = vcg::tri::Index(*this, face.cV(0)); + const size_t v1 = vcg::tri::Index(*this, face.cV(1)); + const size_t v2 = vcg::tri::Index(*this, face.cV(2)); + faces.row(fi) = Eigen::Vector3i(v0, v1, v2); + } + return faces; +} + +bool VCGTriMesh::savePly(const std::string plyFilename) { + // 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) { + return false; + } + return true; +} + +VCGTriMesh::VCGTriMesh() {} + +VCGTriMesh::VCGTriMesh(const std::string &filename) { + const std::string extension = std::filesystem::path(filename).extension(); + if (extension == ".ply") { + loadFromPlyFile(filename); + } else if (extension == ".obj") { + vcg::tri::io::ImporterOBJ::Info info; + vcg::tri::io::ImporterOBJ::Open(*this, filename.c_str(), info); + } else if (extension == ".off") { + vcg::tri::io::ImporterOFF::Open(*this, filename.c_str()); + + } else { + std::cerr << "Uknown file extension " << extension << ". Could not open " + << filename << std::endl; + assert(false); + } + vcg::tri::UpdateTopology::AllocateEdge(*this); + vcg::tri::UpdateTopology::FaceFace(*this); + vcg::tri::UpdateTopology::VertexFace(*this); + vcg::tri::UpdateNormal::PerVertexNormalized(*this); +} diff --git a/vcgtrimesh.hpp b/vcgtrimesh.hpp new file mode 100644 index 0000000..245031c --- /dev/null +++ b/vcgtrimesh.hpp @@ -0,0 +1,40 @@ +#ifndef VCGTRIMESH_HPP +#define VCGTRIMESH_HPP +#include +#include + +using VertexIndex = size_t; +class VCGTriMeshVertex; +class VCGTriMeshEdge; +class VCGTriMeshFace; +struct VCGTriMeshUsedTypes + : public vcg::UsedTypes::AsVertexType, + vcg::Use::AsEdgeType, + vcg::Use::AsFaceType> {}; +class VCGTriMeshVertex + : public vcg::Vertex {}; +class VCGTriMeshFace + : public vcg::Face {}; +class VCGTriMeshEdge + : public vcg::Edge {}; + +class VCGTriMesh : public vcg::tri::TriMesh, + std::vector, + std::vector> { +public: + VCGTriMesh(); + VCGTriMesh(const std::string &filename); + void loadFromPlyFile(const std::string &filename); + Eigen::MatrixX3d getVertices() const; + Eigen::MatrixX3i getFaces() const; + bool savePly(const std::string plyFilename); + template size_t getIndex(const MeshElement &element) { + return vcg::tri::Index(*this, element); + } +}; + +#endif // VCGTRIMESH_HPP