commit 9668d030eab6d4e28d06f557dd8b800c423c635c Author: Iason Date: Fri Nov 27 12:47:21 2020 +0200 Initial commit 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