From a1f155c0f73b060cf2872f9f412324e696c3f78e Mon Sep 17 00:00:00 2001 From: Iason Date: Thu, 17 Dec 2020 20:06:17 +0200 Subject: [PATCH] Added unit tests for form finder. Refactoring. --- beamformfinder.cpp | 726 +++++++++++++++++------------------- beamformfinder.hpp | 14 +- edgemesh.hpp | 2 +- flatpattern.cpp | 8 + trianglepatterngeometry.cpp | 12 +- 5 files changed, 364 insertions(+), 398 deletions(-) diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 853855e..1441476 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -8,6 +8,197 @@ #include #include +void FormFinder::runUnitTests() { + const std::filesystem::path groundOfTruthFolder{ + "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"}; + + FormFinder formFinder; + formFinder.setTotalResidualForcesNormThreshold(1); + + // First example of the paper + VCGEdgeMesh beam; + // const size_t spanGridSize = 11; + // mesh.createSpanGrid(spanGridSize); + beam.loadFromPly( + std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply")); + std::unordered_map> fixedVertices; + fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; + fixedVertices[beam.VN() - 1] = std::unordered_set{1, 2}; + std::unordered_map nodalForces{ + {2, Vector6d({0, 1000, 1000, 0, 0, 0})}}; + // Forced displacements + std::unordered_map nodalForcedDisplacements; + nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)}); + + SimulationJob beamSimulationJob{ + std::make_shared(beam), + // SimulationJob::constructFixedVerticesSpanGrid(spanGridSize, + // mesh.VN()), + fixedVertices, nodalForces, nodalForcedDisplacements}; + beamSimulationJob.mesh->setBeamMaterial(0.3, 200); + assert(CrossSectionType::name == CylindricalBeamDimensions::name); + + beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + SimulationResults simpleBeam_simulationResults = + formFinder.executeSimulation(beamSimulationJob, false, true); + simpleBeam_simulationResults.label = "SimpleBeam"; + simpleBeam_simulationResults.save(); + const std::string simpleBeamGroundOfTruthBinaryFilename = + std::filesystem::path(groundOfTruthFolder) + .append("SimpleBeam_displacements.eigenBin"); + assert(std::filesystem::exists( + std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename))); + Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements; + Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename, + simpleBeam_groundOfTruthDisplacements); + if (!simpleBeam_simulationResults.isEqual( + simpleBeam_groundOfTruthDisplacements)) { + std::cerr << "Error:Beam simulation produces wrong results." << std::endl; + return; + } + + // Second example of the paper + VCGEdgeMesh shortSpanGrid; + // const size_t spanGridSize = 11; + // mesh.createSpanGrid(spanGridSize); + shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) + .append("shortSpanGridshell.ply") + .string()); + + fixedVertices.clear(); + //// Corner nodes + fixedVertices[0] = std::unordered_set{2}; + fixedVertices[4] = std::unordered_set{2}; + fixedVertices[16] = std::unordered_set{2}; + fixedVertices[20] = std::unordered_set{2}; + //// Center node + fixedVertices[10] = std::unordered_set{0, 1, 3, 4, 5}; + + nodalForcedDisplacements.clear(); + nodalForcedDisplacements.insert({{0, Eigen::Vector3d(0.1, 0.1, 0)}, + {4, Eigen::Vector3d(-0.1, 0.1, 0)}, + {16, Eigen::Vector3d(0.1, -0.1, 0)}, + {20, Eigen::Vector3d(-0.1, -0.1, 0)}}); + + SimulationJob shortSpanGridshellSimulationJob{ + std::make_shared(shortSpanGrid), + fixedVertices, + {}, + nodalForcedDisplacements}; + shortSpanGridshellSimulationJob.mesh->setBeamMaterial(0.3, 200); + assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); + shortSpanGridshellSimulationJob.mesh->setBeamCrossSection( + CrossSectionType{0.03, 0.026}); + SimulationResults shortSpanGridshellSimulationResults = + formFinder.executeSimulation(shortSpanGridshellSimulationJob, false, + false); + shortSpanGridshellSimulationResults.label = "ShortSpanGridshell"; + shortSpanGridshellSimulationResults.save(); + + const std::string groundOfTruthBinaryFilename = + std::filesystem::path(groundOfTruthFolder) + .append("ShortSpanGridshell_displacements.eigenBin") + .string(); + assert(std::filesystem::exists( + std::filesystem::path(groundOfTruthBinaryFilename))); + Eigen::MatrixXd shortSpanGridshell_groundOfTruthDisplacements; + Eigen::read_binary(groundOfTruthBinaryFilename, + shortSpanGridshell_groundOfTruthDisplacements); + // shortSpanGridshellSimulationResults.registerForDrawing( + // shortSpanGridshellSimulationJob); + // polyscope::show(); + assert(shortSpanGridshellSimulationResults.isEqual( + shortSpanGridshell_groundOfTruthDisplacements)); + if (!shortSpanGridshellSimulationResults.isEqual( + shortSpanGridshell_groundOfTruthDisplacements)) { + std::cerr << "Error:Short span simulation produces wrong results." + << std::endl; + return; + } + + // Third example of the paper + VCGEdgeMesh longSpanGrid; + longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) + .append("longSpanGridshell.ply") + .string()); + const size_t spanGridSize = 11; + + fixedVertices.clear(); + for (size_t vi = 0; vi < spanGridSize - 1; vi++) { + fixedVertices[vi] = {0, 2}; + } + for (size_t vi = longSpanGrid.VN() - 1 - (spanGridSize - 2); + vi < longSpanGrid.VN(); vi++) { + fixedVertices[vi] = {0, 2}; + } + for (size_t vi = spanGridSize - 1; + vi < longSpanGrid.VN() - 1 - (spanGridSize - 2) - spanGridSize + 1; + vi += spanGridSize + 1) { + fixedVertices[vi] = {1, 2}; + fixedVertices[vi + spanGridSize] = {1, 2}; + } + + nodalForcedDisplacements.clear(); + const size_t horizontalOffset = std::floor((spanGridSize - 2) / 2); + nodalForcedDisplacements.insert( + {{horizontalOffset, Eigen::Vector3d(0, 0.3, 0)}, + {horizontalOffset + 1, Eigen::Vector3d(0, 0.3, 0)}, + {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset, + Eigen::Vector3d(0, -0.3, 0)}, + {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset + 1, + Eigen::Vector3d(0, -0.3, 0)}, + {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2, + Eigen::Vector3d(0.3, 0, 0)}, + {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2, + Eigen::Vector3d(0.3, 0, 0)}, + {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2 + spanGridSize, + Eigen::Vector3d(-0.3, 0, 0)}, + {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2 + + spanGridSize, + Eigen::Vector3d(-0.3, 0, 0)}}); + + SimulationJob longSpanGridshell_simulationJob{ + std::make_shared(longSpanGrid), + fixedVertices, + {}, + nodalForcedDisplacements}; + longSpanGridshell_simulationJob.mesh->setBeamMaterial(0.3, 200); + if (typeid(CrossSectionType) != typeid(CylindricalBeamDimensions)) { + std::cerr << "A cylindrical cross section is expected for running the " + "paper examples." + << std::endl; + } + longSpanGridshell_simulationJob.mesh->setBeamCrossSection( + CrossSectionType{0.03, 0.026}); + SimulationResults longSpanGridshell_simulationResults = + formFinder.executeSimulation(longSpanGridshell_simulationJob, false, + false); + longSpanGridshell_simulationResults.label = "LongSpanGridshell"; + longSpanGridshell_simulationResults.save(); + + const std::string longSpan_groundOfTruthBinaryFilename = + std::filesystem::path(groundOfTruthFolder) + .append("LongSpanGridshell_displacements.eigenBin") + .string(); + assert(std::filesystem::exists( + std::filesystem::path(longSpan_groundOfTruthBinaryFilename))); + Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements; + Eigen::read_binary(longSpan_groundOfTruthBinaryFilename, + longSpanGridshell_groundOfTruthDisplacements); + assert(longSpanGridshell_simulationResults.isEqual( + longSpanGridshell_groundOfTruthDisplacements)); + if (!longSpanGridshell_simulationResults.isEqual( + longSpanGridshell_groundOfTruthDisplacements)) { + std::cerr << "Error:Long span simulation produces wrong results." + << std::endl; + return; + } + + std::cout << "Form finder unit tests succesufully passed." << std::endl; + + // polyscope::show(); +} + void FormFinder::setTotalResidualForcesNormThreshold(double value) { totalResidualForcesNormThreshold = value; } @@ -21,6 +212,7 @@ void FormFinder::reset() { mesh.release(); plotYValues.clear(); plotHandle.reset(); + checkedForMaximumMoment = false; shouldUseKineticEnergyThreshold = false; externalMomentsNorm = 0; } @@ -791,6 +983,7 @@ void FormFinder::updateResidualForcesOnTheFly( const Vector6d &nodeResidualForce = force.residual; // sumOfResidualForces = sumOfResidualForces + nodeResidualForce; const double residualForceNorm = nodeResidualForce.norm(); + const bool shouldTrigger = std::isnan(residualForceNorm); totalResidualForcesNorm += residualForceNorm; } mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm; @@ -1094,16 +1287,15 @@ void FormFinder::updateNodeNormal( * algorithm if the kinetic energy of the system drops below the set * threshold. * */ - static bool wasExecuted{false}; const bool viHasMoments = node.force.external[3] != 0 || node.force.external[4] != 0; - if (!wasExecuted && viHasMoments) { + if (!checkedForMaximumMoment && viHasMoments) { shouldUseKineticEnergyThreshold = true; std::cout << "Maximum moment reached.The Kinetic energy of the system will " "be used as a convergence criterion" << std::endl; - wasExecuted = true; + checkedForMaximumMoment = true; } } else { @@ -1359,92 +1551,17 @@ FormFinder::computeResults(const SimulationMesh &initialMesh) { void FormFinder::draw(const std::string &screenshotsFolder = {}) { // update positions // polyscope::getCurveNetwork("Undeformed edge mesh")->setEnabled(false); - polyscope::getCurveNetwork("Deformed edge mesh") + polyscope::getCurveNetwork(meshPolyscopeLabel) ->updateNodePositions(mesh->getEigenVertices()); - // Per node kinetic energies + + // Vertex quantities 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> externalMoments(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]}; - externalMoments[mesh->getIndex(v)] = {nodeExternalForce[3], - nodeExternalForce[4], 0}; - } - 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") - ->addNodeVectorQuantity("External Moments", externalMoments) - ->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( @@ -1456,7 +1573,26 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { std::vector nRs(mesh->VN()); std::vector theta2(mesh->VN()); std::vector theta3(mesh->VN()); + std::vector> residualForces(mesh->VN()); + std::vector residualForcesNorm(mesh->VN()); + std::vector accelerationX(mesh->VN()); for (const VertexType &v : mesh->vert) { + kineticEnergies[mesh->getIndex(v)] = mesh->nodes[v].kineticEnergy; + const VectorType n = v.cN(); + nodeNormals[mesh->getIndex(v)] = {n[0], n[1], n[2]}; + // 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]}; + externalMoments[mesh->getIndex(v)] = {nodeExternalForce[3], + nodeExternalForce[4], 0}; + internalAxialForces[mesh->getIndex(v)] = {nodeforce[0], nodeforce[1], + nodeforce[2]}; const Node &node = mesh->nodes[v]; const Vector6d nodeForceFirst = node.force.internalFirstBending * (-1); internalFirstBendingTranslationForces[mesh->getIndex(v)] = { @@ -1470,116 +1606,89 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { 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); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("Kinetic Energy", kineticEnergies) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Node normals", nodeNormals) + ->setEnabled(true); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Internal force", internalForces) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("External force", externalForces) + ->setEnabled(true); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("External Moments", externalMoments) + ->setEnabled(true); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("Internal force norm", internalForcesNorm) + ->setEnabled(true); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Internal Axial force", internalAxialForces) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("First bending force-Translation", + internalFirstBendingTranslationForces) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("First bending force-Rotation", + internalFirstBendingRotationForces) + ->setEnabled(false); - // 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") + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Second bending force-Translation", + internalSecondBendingTranslationForces) ->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") + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Second bending force-Rotation", + internalSecondBendingRotationForces); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->getQuantity("Second bending force-Rotation") ->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") + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("nR", nRs) ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("theta3", theta3) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("theta2", theta2) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeVectorQuantity("Residual force", residualForces) + ->setEnabled(false); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addNodeScalarQuantity("Residual force norm", residualForcesNorm) + ->setEnabled(false); + polyscope::getCurveNetwork()->addNodeScalarQuantity("Node acceleration x", + accelerationX); + + // Edge quantities + std::vector A(mesh->EN()); + std::vector J(mesh->EN()); + std::vector I2(mesh->EN()); + std::vector I3(mesh->EN()); + for (const EdgeType &e : mesh->edge) { + const size_t ei = mesh->getIndex(e); + A[ei] = mesh->elements[e].properties.A; + J[ei] = mesh->elements[e].properties.J; + I2[ei] = mesh->elements[e].properties.I2; + I3[ei] = mesh->elements[e].properties.I3; + } + + polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); + polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("J", J); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addEdgeScalarQuantity("I2", I2); + polyscope::getCurveNetwork(meshPolyscopeLabel) + ->addEdgeScalarQuantity("I3", I3); // Specify the callback polyscope::state::userCallback = [&]() { @@ -1618,15 +1727,23 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { } } -SimulationResults -FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, - const bool &shouldDraw, - const size_t &maxDRMIterations) { +SimulationResults FormFinder::executeSimulation( + const SimulationJob &job, const bool &beVerbose, const bool &shouldDraw, + const bool &shouldCreatePlots, const size_t &maxDRMIterations) { assert(job.mesh.operator bool()); auto t1 = std::chrono::high_resolution_clock::now(); reset(); mShouldDraw = shouldDraw; + // if(!job.nodalExternalForces.empty()){ + // double externalForcesNorm=0; + // for(const auto& externalForce:job.nodalExternalForces) + // { + // externalForcesNorm+=externalForce.norm(); + // } + // setTotalResidualForcesNormThreshold(externalForcesNorm); + // } + constrainedVertices = job.fixedVertices; if (!job.nodalForcedDisplacements.empty()) { for (std::pair viDisplPair : @@ -1649,13 +1766,12 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, } computeRigidSupports(); - const std::string deformedMeshName = "Deformed edge mesh"; - if (mShouldDraw) { + if (mShouldDraw || true) { if (!polyscope::state::initialized) { initPolyscope(); } - polyscope::registerCurveNetwork(deformedMeshName, mesh->getEigenVertices(), - mesh->getEigenEdges()); + polyscope::registerCurveNetwork( + meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges()); // registerWorldAxes(); } for (auto fixedVertex : job.fixedVertices) { @@ -1670,7 +1786,7 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, } updateNodalExternalForces(job.nodalExternalForces, constrainedVertices); - while (mCurrentSimulationStep < maxDRMIterations) { + while (maxDRMIterations == 0 || mCurrentSimulationStep < maxDRMIterations) { // while (true) { updateNormalDerivatives(); updateT1Derivatives(); @@ -1701,24 +1817,40 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, if (mCurrentSimulationStep != 0) { history.stepPulse(*mesh); } - if (mShouldDraw && - mCurrentSimulationStep % mDrawingStep == 0 /* && + if (mCurrentSimulationStep > 100000) { + shouldUseKineticEnergyThreshold = true; + } + if (mCurrentSimulationStep>500000 ||(mShouldDraw && + mCurrentSimulationStep % mDrawingStep == 0) /* && currentSimulationStep > maxDRMIterations*/) { // std::string saveTo = std::filesystem::current_path() // .append("Debugging_files") // .append("Screenshots") // .string(); // draw(saveTo); + draw(); + static bool wasExecuted = false; + if (mCurrentSimulationStep > 500000) { + FormFinder debug; + debug.executeSimulation(job, true, true, true); + wasExecuted = true; + } std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm << std::endl; std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy << std::endl; - const auto yValues = history.residualForces; + auto yValues = history.residualForces; auto xPlot = matplot::linspace(0, yValues.size(), yValues.size()); plotHandle = matplot::scatter(xPlot, yValues); - const std::string label = "Residual forces"; + std::string label = "Residual forces"; plotHandle->legend_string(label); - draw(); + + // yValues = history.kineticEnergy; + // plotHandle = matplot::scatter(xPlot, yValues); + // label = "Log of Kinetic energy"; + // plotHandle->legend_string(label); + + // shouldUseKineticEnergyThreshold = true; } // t = t + Dt; mCurrentSimulationStep++; @@ -1728,8 +1860,7 @@ currentSimulationStep > maxDRMIterations*/) { // << std::endl; if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy) { if (/*mesh.totalPotentialEnergykN < 10 ||*/ - // mesh->totalResidualForcesNorm < - // totalResidualForcesNormThreshold || + mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold || (shouldUseKineticEnergyThreshold && mesh->currentTotalKineticEnergy < totalKineticEnergyThreshold)) { if (beVerbose) { @@ -1744,8 +1875,7 @@ currentSimulationStep > maxDRMIterations*/) { << " kNm" << std::endl; } if (shouldUseKineticEnergyThreshold) { - std::cout << "Warning: Since maximum applied external moment reached " - "the kinetic energy of the system was " + std::cout << "Warning: The kinetic energy of the system was " " used as a convergence criterion" << std::endl; } @@ -1775,197 +1905,15 @@ currentSimulationStep > maxDRMIterations*/) { // mesh.printVertexCoordinates(mesh.VN() / 2); if (mShouldDraw) { draw(); - // if (!polyscope::hasCurveNetwork(deformedMeshName)) { - polyscope::removeCurveNetwork(deformedMeshName); - // } + } + if (polyscope::hasCurveNetwork(meshPolyscopeLabel)) { + polyscope::removeCurveNetwork(meshPolyscopeLabel); } SimulationResults results = computeResults(*job.mesh); + results.label = job.mesh->getLabel(); + if (shouldCreatePlots) { + SimulationResultsReporter reporter; + reporter.reportResults({results}, "Results", job.mesh->getLabel()); + } return results; } - -void FormFinder::runUnitTests() { - const std::filesystem::path groundOfTruthFolder{ - "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"}; - - FormFinder formFinder; - formFinder.setTotalResidualForcesNormThreshold(1); - - // First example of the paper - VCGEdgeMesh beam; - // const size_t spanGridSize = 11; - // mesh.createSpanGrid(spanGridSize); - beam.loadFromPly( - std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply")); - std::unordered_map> fixedVertices; - fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; - fixedVertices[beam.VN() - 1] = std::unordered_set{1, 2}; - std::unordered_map nodalForces{ - {2, Vector6d({0, 1000, 1000, 0, 0, 0})}}; - // Forced displacements - std::unordered_map nodalForcedDisplacements; - nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)}); - - SimulationJob beamSimulationJob{ - std::make_shared(beam), - // SimulationJob::constructFixedVerticesSpanGrid(spanGridSize, - // mesh.VN()), - fixedVertices, nodalForces, nodalForcedDisplacements}; - beamSimulationJob.mesh->setBeamMaterial(0.3, 200); - assert(CrossSectionType::name == CylindricalBeamDimensions::name); - - beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); - SimulationResults simpleBeam_simulationResults = - formFinder.executeSimulation(beamSimulationJob, false, true); - simpleBeam_simulationResults.label = "SimpleBeam"; - simpleBeam_simulationResults.save(); - const std::string simpleBeamGroundOfTruthBinaryFilename = - std::filesystem::path(groundOfTruthFolder) - .append("SimpleBeam_displacements.eigenBin"); - assert(std::filesystem::exists( - std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename))); - Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements; - Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename, - simpleBeam_groundOfTruthDisplacements); - if (!simpleBeam_simulationResults.isEqual( - simpleBeam_groundOfTruthDisplacements)) { - std::cerr << "Error:Beam simulation produces wrong results." << std::endl; - return; - } - - // Second example of the paper - VCGEdgeMesh shortSpanGrid; - // const size_t spanGridSize = 11; - // mesh.createSpanGrid(spanGridSize); - shortSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) - .append("shortSpanGridshell.ply") - .string()); - - fixedVertices.clear(); - //// Corner nodes - fixedVertices[0] = std::unordered_set{2}; - fixedVertices[4] = std::unordered_set{2}; - fixedVertices[16] = std::unordered_set{2}; - fixedVertices[20] = std::unordered_set{2}; - //// Center node - fixedVertices[10] = std::unordered_set{0, 1, 3, 4, 5}; - - nodalForcedDisplacements.clear(); - nodalForcedDisplacements.insert({{0, Eigen::Vector3d(0.1, 0.1, 0)}, - {4, Eigen::Vector3d(-0.1, 0.1, 0)}, - {16, Eigen::Vector3d(0.1, -0.1, 0)}, - {20, Eigen::Vector3d(-0.1, -0.1, 0)}}); - - SimulationJob shortSpanGridshellSimulationJob{ - std::make_shared(shortSpanGrid), - fixedVertices, - {}, - nodalForcedDisplacements}; - shortSpanGridshellSimulationJob.mesh->setBeamMaterial(0.3, 200); - assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); - shortSpanGridshellSimulationJob.mesh->setBeamCrossSection( - CrossSectionType{0.03, 0.026}); - SimulationResults shortSpanGridshellSimulationResults = - formFinder.executeSimulation(shortSpanGridshellSimulationJob, false, - false); - shortSpanGridshellSimulationResults.label = "ShortSpanGridshell"; - shortSpanGridshellSimulationResults.save(); - - const std::string groundOfTruthBinaryFilename = - std::filesystem::path(groundOfTruthFolder) - .append("ShortSpanGridshell_displacements.eigenBin") - .string(); - assert(std::filesystem::exists( - std::filesystem::path(groundOfTruthBinaryFilename))); - Eigen::MatrixXd shortSpanGridshell_groundOfTruthDisplacements; - Eigen::read_binary(groundOfTruthBinaryFilename, - shortSpanGridshell_groundOfTruthDisplacements); - // shortSpanGridshellSimulationResults.registerForDrawing( - // shortSpanGridshellSimulationJob); - // polyscope::show(); - assert(shortSpanGridshellSimulationResults.isEqual( - shortSpanGridshell_groundOfTruthDisplacements)); - if (!shortSpanGridshellSimulationResults.isEqual( - shortSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Short span simulation produces wrong results." - << std::endl; - return; - } - - // Third example of the paper - VCGEdgeMesh longSpanGrid; - longSpanGrid.loadFromPly(std::filesystem::path(groundOfTruthFolder) - .append("longSpanGridshell.ply") - .string()); - const size_t spanGridSize = 11; - - fixedVertices.clear(); - for (size_t vi = 0; vi < spanGridSize - 1; vi++) { - fixedVertices[vi] = {0, 2}; - } - for (size_t vi = longSpanGrid.VN() - 1 - (spanGridSize - 2); - vi < longSpanGrid.VN(); vi++) { - fixedVertices[vi] = {0, 2}; - } - for (size_t vi = spanGridSize - 1; - vi < longSpanGrid.VN() - 1 - (spanGridSize - 2) - spanGridSize + 1; - vi += spanGridSize + 1) { - fixedVertices[vi] = {1, 2}; - fixedVertices[vi + spanGridSize] = {1, 2}; - } - - nodalForcedDisplacements.clear(); - const size_t horizontalOffset = std::floor((spanGridSize - 2) / 2); - nodalForcedDisplacements.insert( - {{horizontalOffset, Eigen::Vector3d(0, 0.3, 0)}, - {horizontalOffset + 1, Eigen::Vector3d(0, 0.3, 0)}, - {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset, - Eigen::Vector3d(0, -0.3, 0)}, - {spanGridSize * (spanGridSize + 1) - 2 + horizontalOffset + 1, - Eigen::Vector3d(0, -0.3, 0)}, - {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2, - Eigen::Vector3d(0.3, 0, 0)}, - {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2, - Eigen::Vector3d(0.3, 0, 0)}, - {std::floor(spanGridSize / 2) * (spanGridSize + 1) - 2 + spanGridSize, - Eigen::Vector3d(-0.3, 0, 0)}, - {(std::floor(spanGridSize / 2) + 1) * (spanGridSize + 1) - 2 + - spanGridSize, - Eigen::Vector3d(-0.3, 0, 0)}}); - - SimulationJob longSpanGridshell_simulationJob{ - std::make_shared(longSpanGrid), - fixedVertices, - {}, - nodalForcedDisplacements}; - longSpanGridshell_simulationJob.mesh->setBeamMaterial(0.3, 200); - assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); - longSpanGridshell_simulationJob.mesh->setBeamCrossSection( - CrossSectionType{0.03, 0.026}); - SimulationResults longSpanGridshell_simulationResults = - formFinder.executeSimulation(longSpanGridshell_simulationJob, false, - false); - longSpanGridshell_simulationResults.label = "LongSpanGridshell"; - longSpanGridshell_simulationResults.save(); - - const std::string longSpan_groundOfTruthBinaryFilename = - std::filesystem::path(groundOfTruthFolder) - .append("LongSpanGridshell_displacements.eigenBin") - .string(); - assert(std::filesystem::exists( - std::filesystem::path(longSpan_groundOfTruthBinaryFilename))); - Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements; - Eigen::read_binary(longSpan_groundOfTruthBinaryFilename, - longSpanGridshell_groundOfTruthDisplacements); - assert(longSpanGridshell_simulationResults.isEqual( - longSpanGridshell_groundOfTruthDisplacements)); - if (!longSpanGridshell_simulationResults.isEqual( - longSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Long span simulation produces wrong results." - << std::endl; - return; - } - - std::cout << "Form finder unit tests succesufully passed." << std::endl; - - // polyscope::show(); -} diff --git a/beamformfinder.hpp b/beamformfinder.hpp index e7d31c5..70397ff 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -27,16 +27,18 @@ private: const double Dtini{0.1}; double Dt{Dtini}; const double xi{0.9969}; - double totalResidualForcesNormThreshold{0.001}; - bool shouldUseKineticEnergyThreshold{true}; - const double totalKineticEnergyThreshold{1e-11}; + double totalResidualForcesNormThreshold{1e-5}; + bool shouldUseKineticEnergyThreshold{false}; + bool checkedForMaximumMoment{false}; + const double totalKineticEnergyThreshold{1e-9}; double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; - int mDrawingStep{5}; + int mDrawingStep{1}; bool mShouldDraw{false}; matplot::line_handle plotHandle; std::vector plotYValues; + const std::string meshPolyscopeLabel{"Simulation mesh"}; std::unique_ptr mesh; std::unordered_map> constrainedVertices; @@ -186,9 +188,9 @@ public: FormFinder(); SimulationResults executeSimulation( const SimulationJob &job, const bool &beVerbose = false, - const bool &shouldDraw = false, + const bool &shouldDraw = false, const bool &createPlots = false, const size_t &maxDRMIterations = FormFinder::maxDRMIterations); - inline static const size_t maxDRMIterations{500000}; + inline static const size_t maxDRMIterations{0}; static void runUnitTests(); void setTotalResidualForcesNormThreshold(double value); diff --git a/edgemesh.hpp b/edgemesh.hpp index 9f8ab96..18588f5 100644 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -21,7 +21,7 @@ struct VCGEdgeMeshUsedTypes class VCGEdgeMeshVertexType : public vcg::Vertex {}; + vcg::vertex::Color4b, vcg::vertex::VEAdj> {}; class VCGEdgeMeshEdgeType : public vcg::Edge FlatPatternGeometry::getVertices() const {} +std::vector FlatPatternGeometry::getVertices() const { + std::vector verts(VN()); + for (size_t vi = 0; vi < VN(); vi++) { + verts[vi] = vert[vi].cP(); + } + return verts; +} FlatPatternGeometry FlatPatternGeometry::createTile(FlatPatternGeometry &pattern) { @@ -235,7 +241,9 @@ std::vector FlatPatternGeometry::constructVertexVector( } } - if (numberOfNodesPerSlot[4] != 0) { + if (numberOfNodesPerSlot[4] == 1) { + vertices.push_back(vcg::Point3d(0, -triangleHeight, 0)); + } else 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];