From db500efd056cdf665f838872f0aa8ffeeed33c17 Mon Sep 17 00:00:00 2001 From: Iason Date: Mon, 14 Dec 2020 18:02:54 +0200 Subject: [PATCH] Added tests for the form finder for the 3 cases of the paper as a static function of the form finder. --- beam.hpp | 2 + beamformfinder.cpp | 262 ++++++++++-- beamformfinder.hpp | 30 +- edgemesh.cpp | 14 +- elementalmesh.hpp | 4 +- .../LongSpanGridshell_displacements.eigenBin | Bin 0 -> 6736 bytes .../ShortSpanGridshell_displacements.eigenBin | Bin 0 -> 1024 bytes .../SimpleBeam_displacements.eigenBin | Bin 0 -> 256 bytes .../longSpanGridshell.ply | 374 ++++++++++++++++++ .../shortSpanGridshell.ply | 69 ++++ formFinder_unitTestFiles/simpleBeam.ply | 27 ++ simulationhistoryplotter.hpp | 2 +- simulationresult.hpp | 68 ++-- utilities.hpp | 102 ++--- 14 files changed, 830 insertions(+), 124 deletions(-) create mode 100644 formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin create mode 100644 formFinder_unitTestFiles/ShortSpanGridshell_displacements.eigenBin create mode 100644 formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin create mode 100644 formFinder_unitTestFiles/longSpanGridshell.ply create mode 100644 formFinder_unitTestFiles/shortSpanGridshell.ply create mode 100644 formFinder_unitTestFiles/simpleBeam.ply diff --git a/beam.hpp b/beam.hpp index 4851276..e965c9c 100644 --- a/beam.hpp +++ b/beam.hpp @@ -6,6 +6,7 @@ #include struct RectangularBeamDimensions { + inline static std::string name{"Rectangular"}; float b; float h; RectangularBeamDimensions(const float &width, const float &height) @@ -16,6 +17,7 @@ struct RectangularBeamDimensions { }; struct CylindricalBeamDimensions { + inline static std::string name{"Cylindrical"}; float od; // Cylinder outside diameter float id; // Cylinder inside diameter diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 1d39e20..853855e 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -8,6 +8,10 @@ #include #include +void FormFinder::setTotalResidualForcesNormThreshold(double value) { + totalResidualForcesNormThreshold = value; +} + void FormFinder::reset() { Dt = Dtini; mCurrentSimulationStep = 0; @@ -17,6 +21,8 @@ void FormFinder::reset() { mesh.release(); plotYValues.clear(); plotHandle.reset(); + shouldUseKineticEnergyThreshold = false; + externalMomentsNorm = 0; } VectorType FormFinder::computeDisplacementDifferenceDerivative( @@ -779,13 +785,17 @@ void FormFinder::updateResidualForcesOnTheFly( force.residual = force.residual + (internalForcePair.second * -1); } } + Vector6d sumOfResidualForces(0); for (size_t vi = 0; vi < mesh->VN(); vi++) { - const Vector6d &nodeResidualForce = mesh->nodes[vi].force.residual; + Node::Forces &force = mesh->nodes[vi].force; + const Vector6d &nodeResidualForce = force.residual; + // sumOfResidualForces = sumOfResidualForces + nodeResidualForce; const double residualForceNorm = nodeResidualForce.norm(); totalResidualForcesNorm += residualForceNorm; } mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm; mesh->totalResidualForcesNorm = totalResidualForcesNorm; + // mesh->totalResidualForcesNorm = sumOfResidualForces.norm(); // plotYValues.push_back(totalResidualForcesNorm); // auto xPlot = matplot::linspace(0, plotYValues.size(), plotYValues.size()); @@ -797,6 +807,7 @@ void FormFinder::updateNodalExternalForces( const std::unordered_map> &fixedVertices) { + externalMomentsNorm = 0; for (const std::pair &nodalForce : nodalForces) { const VertexIndex nodeIndex = nodalForce.first; const bool isNodeConstrained = fixedVertices.contains(nodeIndex); @@ -810,6 +821,12 @@ void FormFinder::updateNodalExternalForces( } nodalExternalForce[dofi] = nodalForce.second[dofi]; } + externalMomentsNorm += std::sqrt(pow(nodalExternalForce[3], 2) + + pow(nodalExternalForce[4], 2)); + // CoordType v = (mesh->vert[nodeIndex].cP() - + // mesh->vert[0].cP()).Normalize(); const double forceMagnitude = 0.1; + // nodalExternalForce[3] = v[0] * forceMagnitude; + // nodalExternalForce[4] = v[1] * forceMagnitude; node.force.external = nodalExternalForce; } } @@ -1032,7 +1049,7 @@ void FormFinder::updateNodePosition( node.previousLocation = previousLocation; v.P() = node.initialLocation + displacementVector; if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) { - VectorType desiredInitialDisplacement(0, 0, 0.001); + VectorType desiredInitialDisplacement(0, 0, 0.01); v.P() = v.P() + desiredInitialDisplacement; } } @@ -1061,11 +1078,34 @@ void FormFinder::updateNodeNormal( const double &nx = v.N()[0], ny = v.N()[1], nz = v.N()[2]; const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); VectorType n = v.N(); - const bool shouldBreak = mCurrentSimulationStep == 118 && vi == 1; + const bool shouldBreak = mCurrentSimulationStep == 118 && vi == 3; if (nxnyMagnitude > 1) { VectorType newNormal(nx / std::sqrt(nxnyMagnitude), ny / std::sqrt(nxnyMagnitude), 0); v.N() = newNormal; + + /*If an external moment caused the normal to lay on the xy plane this means + * that in order to disable its effect a greater internal force is needed + * than what is possible (the constraint on the z of the normals imposes a + * constraint on the maximum internal force). Because of that the + * totalResidualForcesNorm can't drop below the magnitude of external moment + * applied on vertex vi. In order to allow termination of the simulation + * when the described phenomenon happens we allow the termination of the + * 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) { + shouldUseKineticEnergyThreshold = true; + std::cout + << "Maximum moment reached.The Kinetic energy of the system will " + "be used as a convergence criterion" + << std::endl; + wasExecuted = true; + } + } else { const double nzSquared = 1.0 - nxnyMagnitude; const double nz = std::sqrt(nzSquared); @@ -1358,6 +1398,7 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { 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 @@ -1369,6 +1410,8 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { 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); @@ -1380,6 +1423,10 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { 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") @@ -1651,6 +1698,9 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, } updateElementalLengths(); + if (mCurrentSimulationStep != 0) { + history.stepPulse(*mesh); + } if (mShouldDraw && mCurrentSimulationStep % mDrawingStep == 0 /* && currentSimulationStep > maxDRMIterations*/) { @@ -1661,11 +1711,15 @@ currentSimulationStep > maxDRMIterations*/) { // draw(saveTo); std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm << std::endl; + std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy + << std::endl; + const auto yValues = history.residualForces; + auto xPlot = matplot::linspace(0, yValues.size(), yValues.size()); + plotHandle = matplot::scatter(xPlot, yValues); + const std::string label = "Residual forces"; + plotHandle->legend_string(label); draw(); } - if (mCurrentSimulationStep != 0) { - history.stepPulse(*mesh); - } // t = t + Dt; mCurrentSimulationStep++; // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy @@ -1674,7 +1728,10 @@ 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) { std::cout << "Convergence after " << mCurrentSimulationStep << " steps" << std::endl; @@ -1686,6 +1743,12 @@ currentSimulationStep > maxDRMIterations*/) { std::cout << "Total potential:" << mesh->totalPotentialEnergykN << " kNm" << std::endl; } + if (shouldUseKineticEnergyThreshold) { + std::cout << "Warning: Since maximum applied external moment reached " + "the kinetic energy of the system was " + " used as a convergence criterion" + << std::endl; + } break; // } } @@ -1721,37 +1784,188 @@ currentSimulationStep > maxDRMIterations*/) { } void FormFinder::runUnitTests() { + const std::filesystem::path groundOfTruthFolder{ + "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"}; + FormFinder formFinder; - VCGEdgeMesh mesh; + formFinder.setTotalResidualForcesNormThreshold(1); + + // First example of the paper + VCGEdgeMesh beam; // const size_t spanGridSize = 11; // mesh.createSpanGrid(spanGridSize); - mesh.loadFromPly("/home/iason/Models/simple_beam_paper_example.ply"); + beam.loadFromPly( + std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply")); std::unordered_map> fixedVertices; fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; - fixedVertices[mesh.VN() - 1] = std::unordered_set{1, 2}; + 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({mesh.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)}); + nodalForcedDisplacements.insert({beam.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)}); SimulationJob beamSimulationJob{ - std::make_shared(mesh), + std::make_shared(beam), // SimulationJob::constructFixedVerticesSpanGrid(spanGridSize, // mesh.VN()), fixedVertices, nodalForces, nodalForcedDisplacements}; beamSimulationJob.mesh->setBeamMaterial(0.3, 200); - assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); - beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); - SimulationResults beamSimulationResults = - formFinder.executeSimulation(beamSimulationJob, true, false); + assert(CrossSectionType::name == CylindricalBeamDimensions::name); - const bool testSuccessful = - beamSimulationResults.displacements[2][0] - (-0.09556) < 1e-5 && - beamSimulationResults.displacements[2][1] - 0.40666 < 1e-5 && - beamSimulationResults.displacements[2][2] - 0.40666 < 1e-5; - assert(testSuccessful); - beamSimulationResults.simulationLabel = "Beam"; - beamSimulationResults.registerForDrawing(beamSimulationJob); - polyscope::show(); + 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 2e6c08f..e7d31c5 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -27,7 +27,10 @@ private: const double Dtini{0.1}; double Dt{Dtini}; const double xi{0.9969}; - const double totalResidualForcesNormThreshold{0.01}; + double totalResidualForcesNormThreshold{0.001}; + bool shouldUseKineticEnergyThreshold{true}; + const double totalKineticEnergyThreshold{1e-11}; + double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; int mDrawingStep{5}; bool mShouldDraw{false}; @@ -188,6 +191,7 @@ public: inline static const size_t maxDRMIterations{500000}; static void runUnitTests(); + void setTotalResidualForcesNormThreshold(double value); }; template PointType Cross(PointType p1, PointType p2) { @@ -211,28 +215,4 @@ inline bool TriggerBreakpoint(const VertexIndex &vi, const EdgeIndex &ei) { 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 index 4ca26fd..fa1fa44 100644 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -12,7 +12,17 @@ Eigen::MatrixX3d VCGEdgeMesh::getEigenEdgeNormals() const { return eigenEdgeNormals; } -bool VCGEdgeMesh::savePly(const std::string plyFilename) {} +bool VCGEdgeMesh::savePly(const std::string plyFilename) { + unsigned int mask = 0; + mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; + mask |= nanoply::NanoPlyWrapper::IO_EDGEINDEX; + mask |= nanoply::NanoPlyWrapper::IO_VERTNORMAL; + if (nanoply::NanoPlyWrapper::SaveModel( + plyFilename.c_str(), *this, mask, false) != 0) { + return false; + } + return true; +} void VCGEdgeMesh::GeneratedRegularSquaredPattern( const double angleDeg, std::vector> &pattern, @@ -90,7 +100,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t squareGridDimension) { } bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth, - const size_t desiredHeight) { + const size_t desiredHeight) { std::cout << "Grid of dimensions:" << desiredWidth << "," << desiredHeight << std::endl; const VCGEdgeMesh::CoordType n(0, 0, 1); diff --git a/elementalmesh.hpp b/elementalmesh.hpp index 463aa98..e60ecdf 100644 --- a/elementalmesh.hpp +++ b/elementalmesh.hpp @@ -7,8 +7,8 @@ struct Element; struct Node; -using CrossSectionType = RectangularBeamDimensions; -// using CrossSectionType = CylindricalBeamDimensions; +// using CrossSectionType = RectangularBeamDimensions; +using CrossSectionType = CylindricalBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: diff --git a/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin new file mode 100644 index 0000000000000000000000000000000000000000..66491894beba32304f33a475f7a9332d3605ae1d GIT binary patch literal 6736 zcmb`Ic{JDA|Hms4MM2qD4ry(u%SqlnB|wh{;ZP;ZqD+ z42sYvH9`wzDMjh~cJFQes59qx&hLzW?s2`}kMlmyjn6RpBgD@C`lzqSbo-ef4Q4MQ z45ZCN;pp#a`g67X;pe;uYTe>zV0qYnw^Q5fA$>tcymhq@=;l_AJQ6h#8?3jeUv;xb zk!gRWESYr%wfAQ_su=m>>`c>|yEdWtqJ8dQe$5Sl#H|3hX(j%w;?8`nvGws^_qSGVBPKBy-o z5)by*(9&T37B%w}yeD&R?(WcBs|hX=4(l#L z(%7mVpA1_V;r?X2h5Se~4cj;{&+|cHqMV1k`ZHTRlDGEre65RUayq}IeN#MU2#mR< zzBC8(d_4}V5G=*(!y_(xKe!L?N8BHnT37)7(uo1{&n80WA+LaYodIBa#yG~&^)NWk zu$fl-rx32HcE2dFEL&CB`bcm2wuvP7ouwI0Cj_)bVFUCjNGQq=H8Lm zWZi2SEohImF$p6K?gnCiK-Q*~DM{D~S4JlL7JTD=J0XVg;~8JZnes2MqI~msKf>R; zi}H>3QT|BA_Yb7}7vJ#666t*~#=-7HduPIalXK^f9oa`WQ2PBtB4mq?kSe z*;Jnn=I>PI?+)g#7xOnxoceo1h5B39zmE7@HsOu)?U#Nw?N?7L?bm<32VN!%We1z1Ve6#X zCGsJ;0M=XYc5wZYoH;5|t?E8}%j>sOwe0n)HQaq;hPiJ#_uM>V1Noy@a-qHK#jP(- zqv6~4McRp7y>Ar?(jVLB1e>Nn;<;_LYnSA}g6qF+vcv-T9P7D1(z^f*b?sM1bLaA0 z0h96Ez4Czj<3z64mYIKCe8BBTg^ll*_C2BhaKG^vf10+;#=swhtL!>wAG`*&q6d71 z)Z@V`NiwF#GZ7}o{gi%TViL^BI+pFjouRTCA3V7GO($NQkXDNaL0yBbYi+JU_d7GE zySn~A@HNINE;gKg2D(;XJmBl<168$o)(1CU1PAk<_9XWPfTa6w-!*pvL8D^+$F`e@}m5u&iuf)(CeJzW@``o!*{xOyV`<9N~v;d zr7bwm3(MUSa2P~eN4U@Lu!BSTqcgT}r{WF8O|snmVcGbw@!WIaYD=45I@^MlpXSE7 z2kd{~zj(InbCHlZe5xI8sx94DXmQoZKx=h(;gi}3(Yfw@g+2w#8v_NpRRf$}R+b5O ztApqAIO8wkjxH+n(aN*>>!W9L*pA<)k|G?ijgYfapm9HioUo}!?=MZ1j zcnaaG-m6?a#4l_-*EYm2tZaTZ#8-_C6C2{It`Djn;;T-L=o#YUH($SVtV+M4un{b3}>1D6CBU#m$e1zqo<= zi$&C5RIMcbVhi;b+oltL(c=U07sp>G{{Fx}-%s+1Ee}XOahn*)CvKqm#EMvwPrOa@ zi77On=s@#{qBNiALG$?oU&fyF3l(U;Q2#vX7dp~@A&>S8m1)1wf%XfJ(SD%}?H76< zB>nn_?_^5$0JG^HVBJ04n zUn7-W0xGIkON>RT!EE8LqDM07;ji(V*SziE!=zP{F7nSjhs<64@k8frt6GgJ<{6!7#=D-?trIcU)F z3LsnI@9wG&&}a`@(!-tU;Ss_+x%)Rp?E;P5bE_YB$7FHurTQe;DRG?d`AZ6e_^s*n zaQf{vflTiCBDMAnj-?$?@K|2I%(e>*8a-8q^sSg(|c+BZw)|%>~$9>o7ymx&S-1{d}e;+A4vllFDW!j!y zdk>)5=+Vrbb}k0upJjU?;unc?*SY6L`4pMdzIg@Wy&sBd-td3m|1G}Ccm3Jtu<=z) z;Of6R;N#cP9+BL4+C@!le97^6VNR+~=k~#65_p8OijK`WJoFnZIt#-wft&KJ!I z$MTuY^7#wP=W3SEvhVTz6O<*R_A? z*Y|$^e}DgKtGirDf`?I`lLyjD@-eFQ*^+95Oq7%#`>0+y2Yr8-pb zWY{xr!fg@d9N=4D`MAU|6GoI;B^}`9!=}k1&xcRo!B)o`@+~gK7+`qVQ(zhoolm~) z`V^gy?OCSjYgS~U#G(y7rQ>rjaBbVFINfAio4Q_2P(K-r?yoh@7teviae{>>R%F5) z-wh3sG5MgKzVlk;G#*Sh+t#wrwHUTNndzcrU5(|+O7_NAtMFr_Rk-2Ra*SE0z!Mkc z;hm&JNR`_GDLPuTl1 z0d|;FdKy*~fs&AU?9r}r@Ii@t^D3)hZ;e>n$+iYqdslz!u0ID z;0C-Wwsl3!vufO-KIIP$n+o(Z?%46is$#TVdSYc|Mk0EA?RsFCnFy*Os@u(07enH+ zG3|c)D_~v!#*J#vt075Avfnni0p936UZeRJACee<5#tMurF=`qpULJ!2Axx)0ZWcs)b|?DeIRG>sJQb0~@vn&tK3zXjn}5;HwnfgVCGm9;mZD$Y*;{cZcqQ8QX&=JvOm9 z+#a0v3U%Mc^;Ry+Vq*{2uUEqJLM^!-Ox@9#H>H}77E9%f#JIEmc!$(c?*5{Klj;bs z$Efy1u-}$@@2E&LA6$#^bBMj=S8lHc>pqCw;r698X5$;7uqs@B@9g1I(U0($!+R&E zwg&Y6^m-7K!u-oqRgAzHY z*OD?LKqUtcsk&%+nIxm)qU=7KZOOQPoyoKAR|%L^zIWi%&IEkww(yX+FWokK_~HKPCBuS2UjxPxA>)G@o#p z<`d#*KEa;m6CTriLSQS&Cm7OvLIdp=nAVeifimqErvvF1tfc*deA+KqN&5wkv|n&l zlk^MJ1xUZ3qKNbh4$yvq?t8KaVBkUa09rH29)KU+14!FW_5k`6$R0ot-2?EoBYOb6 zX0iv6CQSAKZfve^E%6Az>|)5QocFzD%@&P5YQMM3CpvF8S+e` zK|=8SLND29c&#%kq{Sl|;wFY{dB%O`?YUBr^(F*bk2tF9-iSxvYX*FQ@Nf(fJ~cTu zH3+q~7rpmrI)#aM-Q1!(Tp`-rsV>Xn5^M<33(_wPhv_yNLhChRfY-3o<(+37C_AWa zwtklg?z`{2-uo^D4L+3AkE;rXOWkY7d)$md9nm*N1%j9G1k{T7+W25a?OR#xVa{kb z_R8c>o(`a~_|o)xzAq%4x32L&a0%YEHmi?X5DxQWQhs0lBnoz2xMl5FdmU`{%#C%Z z4Mxe6RXP5}K_I3xUCce}GKPi3#1IDM; zH7<`k4Qc(VJ+hKMz#ohsc<3Amw`XcUx!)cF`oEV%m{rEY9N8|1-IYNYr9Z{yrH?xt zxAu`J=sJn-tztC-Ut{K0z*o za_96Mi*li@A_%V*;69pVeK5?*-&uL)!iaJM2 zF1?p!&#M)34L^@2NJNl(Wq&^B(I}zGoU!`LL}+P_avJ+}j~m2@%FQd zR`Ie=AV-F#lBz@?cxy3kc!C2A`ELehr3jd$3UdN$nBWIx2n&9I!QO~=a%o9v#xYQpy`umKVh8uesmch&+kO!q!nYnj`)%gJMLSw-5;HQi9eU= zE*`~v){?A*a zb+Ot4H`8kN-R!p`Zt37<7@qg8(1MpzliGI}`dpC&eKuPq#Jk^=4u)plj zI<(pWw;TAygx3i!QgmZ7-t)xN(?|n$G!XrRh%2<94BZ^J2xCg?t+mb1=i2P6y*Nkg z&`68Nphsy!Zr2dT`R>#(PfH?f=t?$zUj2^m{>D`W37<{IKzdI}h^~MQVTHf%&BMCk z^fOYWQ4jPUMPwbuJBE3zfAJmi4Kj8=_T52YUAtOifjh_A!e8&RgX!K^5rI=qSik*) z`t}etj6JEjBsMGH$+%(#_u(sWKLDBh)=SaC(QOlZ6#eQ%b|q%@&F7)gb4dnoMw zUY=rRHbH$5hl?J%G4a(%mPVQYdrOa78199Y3`*-fgCaafHu1s xQTHUv`p8mbZT`5oJVFhXi?vgkj3y+R<}sDmv}oaiynTF50-u)FMg`N==s!xi>^uMf literal 0 HcmV?d00001 diff --git a/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin b/formFinder_unitTestFiles/SimpleBeam_displacements.eigenBin new file mode 100644 index 0000000000000000000000000000000000000000..61a9a2344b97dada5f43dc7dd08f1a6f2d7ef6bd GIT binary patch literal 256 zcmZQ&fB-fq4WkmTtv;i0Yx918fjnlX@*VrX+3!=Sn7VuateGHiaz9L6-QIkyeU%sO ze_9k=s$#!o@2;r1Ytx@g_Aq%Ez3x$2+|`PU_J7YWFt=s9W$!cP@1ncEG1V3RlM_5D z%CMhFvVHlUwCnqn_Z^o|-N|S#-FoE^_jGQ1`7h$rj?aI#|K4P){VRd=wPQXEH9-1K zjP~^$Abr;5$0-nP+G#PTqsetEnabled(true); } } - static std::unordered_map> - constructFixedVerticesSpanGrid(const size_t &spanGridSize, - const size_t &gridVertices) { - std::unordered_map> fixedVertices; - for (size_t vi = 0; vi < spanGridSize - 1; vi++) { - fixedVertices[vi] = {0, 1, 2}; - } - for (size_t vi = gridVertices - 1 - (spanGridSize - 2); vi < gridVertices; - vi++) { - fixedVertices[vi] = {0, 1, 2}; - } - for (size_t vi = spanGridSize - 1; - vi < gridVertices - 1 - (spanGridSize - 2) - spanGridSize + 1; - vi += spanGridSize + 1) { - fixedVertices[vi] = {0, 1, 2}; - fixedVertices[vi + spanGridSize] = {0, 1, 2}; - } - - return fixedVertices; - } }; +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 + struct SimulationResults { SimulationHistory history; std::vector displacements; double executionTime{0}; - std::string simulationLabel{"NoLabel"}; + std::string label{"NoLabel"}; void registerForDrawing(const SimulationJob &job) const { polyscope::options::groundPlaneEnabled = false; @@ -136,12 +140,12 @@ struct SimulationResults { } /* else { polyscope::removeAllStructures(); }*/ - const std::string undeformedMeshName = "Undeformed_" + simulationLabel; + const std::string undeformedMeshName = "Undeformed_" + label; polyscope::registerCurveNetwork(undeformedMeshName, job.mesh->getEigenVertices(), job.mesh->getEigenEdges()); - const std::string deformedMeshName = "Deformed_" + simulationLabel; + const std::string deformedMeshName = "Deformed_" + label; polyscope::registerCurveNetwork(deformedMeshName, job.mesh->getEigenVertices(), job.mesh->getEigenEdges()); @@ -214,6 +218,24 @@ struct SimulationResults { // polyscope::show(); } + + void save() { + const std::string filename(label + "_displacements.eigenBin"); + + Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); + Eigen::write_binary(filename, m); + } + + // The comparison of the results happens comparing the 6-dof nodal + // displacements + bool isEqual(const Eigen::MatrixXd &nodalDisplacements) { + assert(nodalDisplacements.cols() == 6); + Eigen::MatrixXd eigenDisplacements = + Utilities::toEigenMatrix(this->displacements); + const double errorNorm = (eigenDisplacements - nodalDisplacements).norm(); + return errorNorm < 1e-10; + // return eigenDisplacements.isApprox(nodalDisplacements); + } }; #endif // SIMULATIONHISTORY_HPP diff --git a/utilities.hpp b/utilities.hpp index 5571974..b1d640d 100644 --- a/utilities.hpp +++ b/utilities.hpp @@ -6,53 +6,6 @@ #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++) { @@ -145,6 +98,61 @@ struct Vector6d : public std::array { } }; +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]); }); +} + +inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { + Eigen::MatrixXd m(v.size(), 6); + + for (size_t vi = 0; vi < v.size(); vi++) { + const Vector6d &vec = v[vi]; + for (size_t i = 0; i < 6; i++) { + m(vi, i) = vec[i]; + } + } + + return m; +} + +// 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 + // namespace ConfigurationFile { // inline void getPlyFilename(const std::string jsonFilepath,