From a67003cdd2ef621a065c370fb03ea681d7cf0a62 Mon Sep 17 00:00:00 2001 From: Iason Date: Thu, 3 Dec 2020 20:56:03 +0200 Subject: [PATCH] Refactoring --- beam.hpp | 2 +- beamformfinder.cpp | 183 +++++++++++++++++++++++++------------------ beamformfinder.hpp | 23 ++++-- edgemesh.cpp | 93 ++++++---------------- edgemesh.hpp | 11 +-- elementalmesh.cpp | 179 +++++++++++++++++++++++++++++++++++++----- elementalmesh.hpp | 52 ++++++------ flatpattern.cpp | 4 +- flatpattern.hpp | 2 +- simulationresult.hpp | 43 +++++++--- 10 files changed, 371 insertions(+), 221 deletions(-) diff --git a/beam.hpp b/beam.hpp index 3657070..463261e 100644 --- a/beam.hpp +++ b/beam.hpp @@ -12,7 +12,7 @@ struct RectangularBeamDimensions { : b(width), h(height) { assert(width > 0 && height > 0); } - RectangularBeamDimensions() : b(1), h(1) {} + RectangularBeamDimensions() : b(0.01), h(0.01) {} }; struct CylindricalElementDimensions { diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 157837a..794d43a 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -10,13 +10,13 @@ void FormFinder::reset() { Dt = Dtini; - t = 0; - currentSimulationStep = 0; + mCurrentSimulationStep = 0; history.clear(); - // theta3Derivatives = - // Eigen::Tensor(DoF::NumDoFType, mesh->VN(), mesh->EN(), - // mesh->VN()); - // theta3Derivatives.setZero(); + constrainedVertices.clear(); + rigidSupports.clear(); + mesh.release(); + plotYValues.clear(); + plotHandle.reset(); } VectorType FormFinder::computeDisplacementDifferenceDerivative( @@ -40,6 +40,7 @@ VectorType FormFinder::computeDisplacementDifferenceDerivative( VectorType FormFinder::computeDerivativeOfNormal( const VertexType &v, const DifferentiateWithRespectTo &dui) const { + const size_t vi = mesh->getIndex(v); VectorType normalDerivative(0, 0, 0); if (&dui.v != &v || (dui.dofi == 0 || dui.dofi == 1 || dui.dofi == 2 || dui.dofi == 5)) { @@ -184,7 +185,9 @@ FormFinder::computeDerivativeT2(const EdgeType &e, const DoFType dofi = dui.dofi; const VertexType &v_j = *e.cV(0); + const size_t vi_j = mesh->getIndex(v_j); const VertexType &v_jplus1 = *e.cV(1); + const size_t vi_jplus1 = mesh->getIndex(v_jplus1); const VectorType r = (v_j.cN() + v_jplus1.cN()).Normalize(); const VectorType derivativeR_j = @@ -231,7 +234,9 @@ FormFinder::computeDerivativeT3(const EdgeType &e, const VectorType &t2 = element.frame.t2; const VectorType t1CrossT2 = Cross(t1, t2); const VertexType &v_j = *e.cV(0); + const size_t vi_j = mesh->getIndex(v_j); const VertexType &v_jplus1 = *e.cV(1); + const size_t vi_jplus1 = mesh->getIndex(v_jplus1); const VectorType derivativeT1_j = dui.dofi < 3 && &v_j == &dui.v ? mesh->elements[e].derivativeT1_j[dui.dofi] @@ -279,6 +284,7 @@ double FormFinder::computeDerivativeTheta1(const EdgeType &e, const VertexIndex &dwrt_evi, const DoFType &dwrt_dofi) const { const VertexType &v = *e.cV(evi); + const size_t vi = mesh->getIndex(v); const Element &element = mesh->elements[e]; const VectorType derivativeT1 = element.derivativeT1[dwrt_evi][dwrt_dofi]; const VectorType derivativeT3 = element.derivativeT3[dwrt_evi][dwrt_dofi]; @@ -617,6 +623,7 @@ void FormFinder::updateResidualForcesOnTheFly( internalForcesContributionFromThisEdge[evi + 2].first = refElemOtherVertexNode.vi; } + const size_t vi = edgeNode.vi; // #pragma omp parallel for schedule(static) num_threads(6) for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { const bool isDofConstrainedFor_ev = @@ -765,11 +772,16 @@ void FormFinder::updateResidualForcesOnTheFly( } } for (size_t vi = 0; vi < mesh->VN(); vi++) { - const double residualForceNorm = (mesh->nodes[vi].force.residual).norm(); + const Vector6d &nodeResidualForce = mesh->nodes[vi].force.residual; + const double residualForceNorm = nodeResidualForce.norm(); totalResidualForcesNorm += residualForceNorm; } mesh->previousTotalResidualForcesNorm = mesh->totalResidualForcesNorm; mesh->totalResidualForcesNorm = totalResidualForcesNorm; + + // plotYValues.push_back(totalResidualForcesNorm); + // auto xPlot = matplot::linspace(0, plotYValues.size(), plotYValues.size()); + // plotHandle = matplot::scatter(xPlot, plotYValues); } void FormFinder::updateNodalExternalForces( @@ -901,6 +913,7 @@ FormFinder::FormFinder() {} void FormFinder::updateNodalMasses() { const double gamma = 0.8; for (VertexType &v : mesh->vert) { + const size_t vi = mesh->getIndex(v); double translationalSumSk = 0; double rotationalSumSk_I2 = 0; double rotationalSumSk_I3 = 0; @@ -955,6 +968,7 @@ void FormFinder::updateNodalMasses() { void FormFinder::updateNodalAccelerations() { for (VertexType &v : mesh->vert) { Node &node = mesh->nodes[v]; + const VertexIndex vi = mesh->getIndex(v); for (DoFType dofi = DoF::Ux; dofi < DoF::NumDoF; dofi++) { if (dofi == DoF::Ux || dofi == DoF::Uy || dofi == DoF::Uz) { node.acceleration[dofi] = @@ -964,10 +978,10 @@ void FormFinder::updateNodalAccelerations() { node.force.residual[dofi] / node.rotationalMass_J; } else if (dofi == DoF::Ny) { node.acceleration[dofi] = - node.force.residual[dofi] / node.rotationalMass_I2; + node.force.residual[dofi] / node.rotationalMass_I3; } else if (dofi == DoF::Nr) { node.acceleration[dofi] = - node.force.residual[dofi] / node.rotationalMass_I3; + node.force.residual[dofi] / node.rotationalMass_I2; } } } @@ -975,6 +989,7 @@ void FormFinder::updateNodalAccelerations() { void FormFinder::updateNodalVelocities() { for (VertexType &v : mesh->vert) { + const VertexIndex vi = mesh->getIndex(v); Node &node = mesh->nodes[v]; node.velocity = node.velocity + node.acceleration * Dt; } @@ -1008,52 +1023,63 @@ void FormFinder::updateNodePosition( } node.previousLocation = previousLocation; v.P() = node.initialLocation + displacementVector; - if (shouldApplyInitialDistortion && currentSimulationStep < 40) { + if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) { VectorType desiredInitialDisplacement(0, 0, 0.1); v.P() = v.P() + desiredInitialDisplacement; } } +void FormFinder::updateNodeNormal( + VertexType &v, + const std::unordered_map> + &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) { + VectorType newNormal(nx / std::sqrt(nxnyMagnitude), + ny / std::sqrt(nxnyMagnitude), 0); + v.N() = newNormal; + } 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()); + } +} + 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()); - } + updateNodeNormal(v, fixedVertices); } updateElementalFrames(); + if (mShouldDraw) { + mesh->updateEigenEdgeAndVertices(); + } } void FormFinder::updateElementalFrames() { @@ -1098,10 +1124,10 @@ void FormFinder::updateKineticEnergy() { // 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); + // 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; } @@ -1244,7 +1270,7 @@ FormFinder::computeResults(const SimulationMesh &initialMesh) { displacements[vi] = mesh->nodes[vi].displacements; } - history.numberOfSteps = currentSimulationStep; + history.numberOfSteps = mCurrentSimulationStep; return SimulationResults{history, displacements}; } @@ -1478,6 +1504,9 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { ImGui::InputInt("Simulation drawing step", &mDrawingStep); // set a int variable + ImGui::Checkbox("Enable drawing", + &mShouldDraw); // set a int variable + ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); ImGui::PopItemWidth(); }; @@ -1493,7 +1522,7 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { } polyscope::screenshot( std::filesystem::path(screenshotsFolder) - .append(std::to_string(currentSimulationStep) + ".png") + .append(std::to_string(mCurrentSimulationStep) + ".png") .string(), false); } else { @@ -1507,6 +1536,8 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, const size_t &maxDRMIterations) { assert(job.mesh.operator bool()); auto t1 = std::chrono::high_resolution_clock::now(); + reset(); + mShouldDraw = shouldDraw; constrainedVertices = job.fixedVertices; if (!job.nodalForcedDisplacements.empty()) { @@ -1526,19 +1557,19 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, std::cout << "Executing simulation for mesh with:" << mesh->VN() << " nodes and " << mesh->EN() << " elements." << std::endl; } - reset(); computeRigidSupports(); - if (shouldDraw) { + + const std::string deformedMeshName = "Deformed edge mesh"; + if (mShouldDraw) { if (!polyscope::state::initialized) { initPolyscope(); } - const std::string deformedMeshName = "Deformed edge mesh"; if (!polyscope::hasCurveNetwork(deformedMeshName)) { polyscope::registerCurveNetwork( deformedMeshName, mesh->getEigenVertices(), mesh->getEigenEdges()); } - registerWorldAxes(); + // registerWorldAxes(); } // const double beamRadius = mesh.getBeamDimensions()[0].od / 2; // std::cout << "BeamRadius:" << beamRadius << std::endl; @@ -1603,7 +1634,7 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, updateNodalExternalForces(job.nodalExternalForces, constrainedVertices); // } - while (currentSimulationStep < maxDRMIterations) { + while (mCurrentSimulationStep < maxDRMIterations) { // while (true) { updateNormalDerivatives(); updateT1Derivatives(); @@ -1665,8 +1696,8 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, // assert(std::isfinite(mesh.currentTotalKineticEnergy)); // mesh.printVertexCoordinates(mesh.VN() / 2); // draw(); - if (shouldDraw && - currentSimulationStep % mDrawingStep == 0 /* && + if (mShouldDraw && + mCurrentSimulationStep % mDrawingStep == 0 /* && currentSimulationStep > maxDRMIterations*/) { // std::string saveTo = std::filesystem::current_path() // .append("Debugging_files") @@ -1691,11 +1722,11 @@ currentSimulationStep > maxDRMIterations*/) { // << std::endl; draw(); } - if (currentSimulationStep != 0) { + if (mCurrentSimulationStep != 0) { history.stepPulse(*mesh); } // t = t + Dt; - currentSimulationStep++; + mCurrentSimulationStep++; // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy // << std::endl; // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm @@ -1704,8 +1735,8 @@ currentSimulationStep > maxDRMIterations*/) { if (/*mesh.totalPotentialEnergykN < 10 ||*/ mesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { if (beVerbose) { - std::cout << "Convergence after " << currentSimulationStep << " steps" - << std::endl; + std::cout << "Convergence after " << mCurrentSimulationStep + << " steps" << std::endl; std::cout << "Residual force of magnitude:" << mesh->previousTotalResidualForcesNorm << std::endl; std::cout << "Kinetic energy:" << mesh->currentTotalKineticEnergy @@ -1737,17 +1768,17 @@ currentSimulationStep > maxDRMIterations*/) { // << " 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( + // for (VertexType &v : mesh->vert) { + // Node &node = mesh->nodes[v]; + // node.displacements = node.displacements - node.velocity * Dt; + // } + // applyDisplacements(constrainedVertices, shouldDraw); + // if (!job.nodalForcedDisplacements.empty()) { + // applyForcedDisplacements( - job.nodalForcedDisplacements); - } - updateElementalLengths(); + // job.nodalForcedDisplacements); + // } + // updateElementalLengths(); resetVelocities(); Dt = Dt * xi; // std::cout << "Residual forces norm:" << @@ -1755,7 +1786,7 @@ currentSimulationStep > maxDRMIterations*/) { // << std::endl; } } - if (currentSimulationStep == maxDRMIterations) { + if (mCurrentSimulationStep == maxDRMIterations) { std::cout << "Did not reach equilibrium before reaching the maximum number " "of DRM steps (" << maxDRMIterations << "). Breaking simulation" << std::endl; @@ -1768,16 +1799,16 @@ currentSimulationStep > maxDRMIterations*/) { << std::endl; std::cout << "Time/(node*iteration) " << simulationDuration / - (static_cast(currentSimulationStep) * mesh->VN()) + (static_cast(mCurrentSimulationStep) * mesh->VN()) << "s" << std::endl; } // mesh.printVertexCoordinates(mesh.VN() / 2); - if (shouldDraw) { + if (mShouldDraw) { draw(); + // if (!polyscope::hasCurveNetwork(deformedMeshName)) { + polyscope::removeCurveNetwork(deformedMeshName); + // } } - // if (!polyscope::hasCurveNetwork(deformedMeshName)) { - // polyscope::removeCurveNetwork(deformedMeshName); - // } // debugger.createPlots(); SimulationResults results = computeResults(*job.mesh); // Eigen::write_binary("optimizedResult.eigenBin", diff --git a/beamformfinder.hpp b/beamformfinder.hpp index bd268da..30cffc0 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -2,6 +2,7 @@ #define BEAMFORMFINDER_HPP #include "elementalmesh.hpp" +#include "matplot/matplot.h" #include "polyscope/curve_network.h" #include "polyscope/polyscope.h" #include "simulationresult.hpp" @@ -25,11 +26,14 @@ 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}; + const double totalResidualForcesNormThreshold{10}; + size_t mCurrentSimulationStep{0}; + int mDrawingStep{5}; + bool mShouldDraw{false}; + matplot::line_handle plotHandle; + std::vector plotYValues; + std::unique_ptr mesh; std::unordered_map> constrainedVertices; @@ -63,7 +67,7 @@ private: const std::unordered_map nodalForcedDisplacements); - ::Vector6d computeElementTorsionalForce( + Vector6d computeElementTorsionalForce( const Element &element, const Vector6d &displacementDifference, const std::unordered_set &constrainedDof); @@ -97,12 +101,12 @@ private: const Vector6d &elementInternalForce, const std::unordered_set &nodalFixedDof); - ::Vector6d computeElementInternalForce( + 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; + Vector6d computeElementAxialForce(const ::EdgeType &e) const; VectorType computeDisplacementDifferenceDerivative( const EdgeType &e, const DifferentiateWithRespectTo &dui) const; double @@ -167,6 +171,11 @@ private: const std::unordered_map> &fixedVertices); + void updateNodeNormal( + VertexType &v, + const std::unordered_map> + &fixedVertices); + public: FormFinder(); SimulationResults executeSimulation( diff --git a/edgemesh.cpp b/edgemesh.cpp index 9c3d5b6..4ca26fd 100644 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -12,36 +12,7 @@ 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; -} +bool VCGEdgeMesh::savePly(const std::string plyFilename) {} void VCGEdgeMesh::GeneratedRegularSquaredPattern( const double angleDeg, std::vector> &pattern, @@ -114,11 +85,11 @@ void VCGEdgeMesh::createSpiral(const float °reesOfArm, // setDefaultAttributes(); } -bool VCGEdgeMesh::createGrid(const size_t squareGridDimension) { - return createGrid(squareGridDimension, squareGridDimension); +bool VCGEdgeMesh::createSpanGrid(const size_t squareGridDimension) { + return createSpanGrid(squareGridDimension, squareGridDimension); } -bool VCGEdgeMesh::createGrid(const size_t desiredWidth, +bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth, const size_t desiredHeight) { std::cout << "Grid of dimensions:" << desiredWidth << "," << desiredHeight << std::endl; @@ -229,44 +200,38 @@ 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) { + if (nanoply::NanoPlyWrapper::LoadModel(plyFilename.c_str(), + *this, mask) != 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; - } +// 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); -} +// const std::vector &edgePropertyVector = +// edgeElemIt->propVec; +// return hasPlyEdgeProperty(plyFilename, edgePropertyVector, +// plyPropertyBeamDimensionsID) && +// hasPlyEdgeProperty(plyFilename, edgePropertyVector, +// plyPropertyBeamMaterialID); +//} bool VCGEdgeMesh::hasPlyEdgeProperty( const std::string &plyFilename, @@ -316,13 +281,7 @@ void VCGEdgeMesh::getEdges(Eigen::MatrixX3d &edgeStartingPoints, } } -VCGEdgeMesh::VCGEdgeMesh() { - handleBeamDimensions = vcg::tri::Allocator::AddPerEdgeAttribute< - CylindricalElementDimensions>(*this, plyPropertyBeamDimensionsID); - handleBeamMaterial = - vcg::tri::Allocator::AddPerEdgeAttribute( - *this, plyPropertyBeamMaterialID); -} +VCGEdgeMesh::VCGEdgeMesh() {} void VCGEdgeMesh::updateEigenEdgeAndVertices() { getEdges(eigenEdges); diff --git a/edgemesh.hpp b/edgemesh.hpp index af2fadd..9f8ab96 100644 --- a/edgemesh.hpp +++ b/edgemesh.hpp @@ -29,11 +29,6 @@ class VCGEdgeMeshEdgeType 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; @@ -75,15 +70,13 @@ public: bool savePly(const std::string plyFilename); - bool createGrid(const size_t squareGridDimension); - bool createGrid(const size_t desiredWidth, const size_t desiredHeight); + bool createSpanGrid(const size_t squareGridDimension); + bool createSpanGrid(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; diff --git a/elementalmesh.cpp b/elementalmesh.cpp index cb91033..752aae7 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -1,5 +1,21 @@ #include "elementalmesh.hpp" +SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) { + vcg::tri::MeshAssert::VertexNormalNormalized(mesh); + + vcg::tri::Append::MeshCopy(*this, mesh); + 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 = mesh.getLabel(); + eigenEdges = mesh.getEigenEdges(); + eigenVertices = mesh.getEigenVertices(); +} + SimulationMesh::SimulationMesh(FlatPattern &pattern) { vcg::tri::MeshAssert::VertexNormalNormalized(pattern); @@ -16,9 +32,8 @@ SimulationMesh::SimulationMesh(FlatPattern &pattern) { eigenVertices = pattern.getEigenVertices(); } -SimulationMesh::SimulationMesh(SimulationMesh &elementalMesh) { - vcg::tri::Append::MeshCopy(*this, - elementalMesh); +SimulationMesh::SimulationMesh(SimulationMesh &mesh) { + vcg::tri::Append::MeshCopy(*this, mesh); elements = vcg::tri::Allocator::GetPerEdgeAttribute( *this, std::string("Elements")); nodes = vcg::tri::Allocator::GetPerVertexAttribute( @@ -27,16 +42,16 @@ SimulationMesh::SimulationMesh(SimulationMesh &elementalMesh) { initializeNodes(); for (size_t ei = 0; ei < EN(); ei++) { - elements[ei] = elementalMesh.elements[ei]; + elements[ei] = mesh.elements[ei]; } - label = label; - eigenEdges = elementalMesh.getEigenEdges(); - eigenVertices = elementalMesh.getEigenVertices(); + label = mesh.label; + eigenEdges = mesh.getEigenEdges(); + eigenVertices = mesh.getEigenVertices(); } void SimulationMesh::computeElementalProperties() { - const std::vector elementalDimensions = + const std::vector elementalDimensions = getBeamDimensions(); const std::vector elementalMaterials = getBeamMaterial(); assert(EN() == elementalDimensions.size() && @@ -102,6 +117,10 @@ void SimulationMesh::initializeElements() { for (const EdgeType &e : edge) { Element &element = elements[e]; element.ei = getIndex(e); + // Initialize dimensions + element.properties.dimensions = CrossSectionType(1, 1); + // Initialize material + element.properties.material = ElementMaterial(); // Initialize lengths const VCGEdgeMesh::CoordType p0 = e.cP(0); const VCGEdgeMesh::CoordType p1 = e.cP(1); @@ -109,16 +128,7 @@ void SimulationMesh::initializeElements() { 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.updateConstFactors(); element.derivativeT1.resize( 2, std::vector(6, VectorType(0, 0, 0))); element.derivativeT2.resize( @@ -152,6 +162,91 @@ void SimulationMesh::updateElementalLengths() { } } +void SimulationMesh::setBeamCrossSection(const double &firstDimension, + const double &secondDimension) { + for (size_t ei = 0; ei < EN(); ei++) { + elements[ei].properties.setDimensions( + CrossSectionType{firstDimension, secondDimension}); + elements[ei].updateConstFactors(); + } +} + +void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { + for (size_t ei = 0; ei < EN(); ei++) { + elements[ei].properties.setMaterial(ElementMaterial{pr, ym}); + elements[ei].updateConstFactors(); + } +} + +std::vector SimulationMesh::getBeamDimensions() { + std::vector beamDimensions(EN()); + for (size_t ei = 0; ei < EN(); ei++) { + beamDimensions[ei] = elements[ei].properties.dimensions; + } + return beamDimensions; +} + +std::vector SimulationMesh::getBeamMaterial() { + std::vector beamMaterial(EN()); + for (size_t ei = 0; ei < EN(); ei++) { + beamMaterial[ei] = elements[ei].properties.material; + } + return beamMaterial; +} + +bool SimulationMesh::loadFromPly(const string &plyFilename) { + this->Clear(); + // assert(plyFileHasAllRequiredFields(plyFilename)); + // Load the ply file + VCGEdgeMesh::PerEdgeAttributeHandle + handleBeamDimensions; + VCGEdgeMesh::PerEdgeAttributeHandle handleBeamMaterial; + handleBeamDimensions = + vcg::tri::Allocator::AddPerEdgeAttribute( + *this, plyPropertyBeamDimensionsID); + handleBeamMaterial = + vcg::tri::Allocator::AddPerEdgeAttribute( + *this, plyPropertyBeamMaterialID); + 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); + 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) != 0) { + return false; + } + return true; +} + +bool SimulationMesh::savePly(const std::string &plyFilename) { + nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; + customAttrib.GetMeshAttrib(plyFilename); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, + getBeamDimensions().data()); + customAttrib.AddEdgeAttribDescriptor( + plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, + getBeamMaterial().data()); + // 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; +} + SimulationMesh::EdgePointer SimulationMesh::getReferenceElement(const VCGEdgeMesh::VertexType &v) { const VertexIndex vi = getIndex(v); @@ -207,3 +302,51 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, } return angle; } + +void Element::Properties::computeMaterialProperties( + const ElementMaterial &material) { + E = material.youngsModulusGPascal * std::pow(10, 9); + G = E / (2 * (1 + material.poissonsRatio)); +} + +void Element::Properties::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; + J = I2 + I3; +} + +void Element::Properties::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; + J = I2 + I3; +} + +void Element::Properties::setDimensions(const CrossSectionType &dimensions) { + this->dimensions = dimensions; + computeDimensionsProperties(dimensions); +} + +void Element::Properties::setMaterial(const ElementMaterial &material) { + this->material = material; + computeMaterialProperties(material); +} + +Element::Properties::Properties(const CrossSectionType &dimensions, + const ElementMaterial &material) + : dimensions(dimensions), material(material) { + computeDimensionsProperties(dimensions); + computeMaterialProperties(material); +} +void Element::updateConstFactors() { + + axialConstFactor = properties.E * properties.A / initialLength; + torsionConstFactor = properties.G * properties.J / initialLength; + firstBendingConstFactor = 2 * properties.E * properties.I2 / initialLength; + secondBendingConstFactor = 2 * properties.E * properties.I3 / initialLength; + int i = 0; + i++; +} diff --git a/elementalmesh.hpp b/elementalmesh.hpp index 682277b..7357309 100644 --- a/elementalmesh.hpp +++ b/elementalmesh.hpp @@ -7,6 +7,7 @@ struct Element; struct Node; +using CrossSectionType = RectangularBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: @@ -15,6 +16,9 @@ private: void initializeElements(); EdgePointer getReferenceElement(const VertexType &v); + const std::string plyPropertyBeamDimensionsID{"beam_dimensions"}; + const std::string plyPropertyBeamMaterialID{"beam_material"}; + public: PerEdgeAttributeHandle elements; PerVertexAttributeHandle nodes; @@ -22,55 +26,43 @@ public: SimulationMesh(ConstVCGEdgeMesh &edgeMesh); SimulationMesh(SimulationMesh &elementalMesh); void updateElementalLengths(); + void setBeamCrossSection(const double &firstDimension, + const double &secondDimension); + void setBeamMaterial(const double &pr, const double &ym); std::vector getIncidentElements(const VCGEdgeMesh::VertexType &v); + bool loadFromPly(const string &plyFilename); + std::vector getBeamDimensions(); + std::vector getBeamMaterial(); double previousTotalKineticEnergy{0}; double previousTotalResidualForcesNorm{0}; double currentTotalKineticEnergy{0}; double totalResidualForcesNorm{0}; double totalPotentialEnergykN{0}; + bool savePly(const std::string &plyFilename); }; struct Element { struct Properties { + CrossSectionType dimensions; + ElementMaterial material; 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 computeMaterialProperties(const ElementMaterial &material); 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; - } + computeDimensionsProperties(const RectangularBeamDimensions &dimensions); + void + computeDimensionsProperties(const CylindricalElementDimensions &dimensions); + void setDimensions(const CrossSectionType &dimensions); + void setMaterial(const ElementMaterial &material); + Properties(const CrossSectionType &dimensions, + const ElementMaterial &material); Properties() {} }; @@ -80,6 +72,8 @@ struct Element { VectorType t3; }; + void updateConstFactors(); + EdgeIndex ei; double length{0}; Properties properties; diff --git a/flatpattern.cpp b/flatpattern.cpp index 381bbbf..acb2629 100644 --- a/flatpattern.cpp +++ b/flatpattern.cpp @@ -17,8 +17,9 @@ FlatPattern::FlatPattern(const string &filename, bool addNormalsIfAbsent) { } vcg::tri::UpdateTopology::VertexEdge(*this); + scale(); - // scale(); + updateEigenEdgeAndVertices(); } FlatPattern::FlatPattern(const std::vector &numberOfNodesPerSlot, @@ -31,6 +32,7 @@ FlatPattern::FlatPattern(const std::vector &numberOfNodesPerSlot, v.N() = CoordType(0, 0, 1); } } + scale(); updateEigenEdgeAndVertices(); } diff --git a/flatpattern.hpp b/flatpattern.hpp index 47b42d5..4878eef 100644 --- a/flatpattern.hpp +++ b/flatpattern.hpp @@ -27,7 +27,7 @@ private: void removeDuplicateVertices(); void scale(); - const double desiredBaseTriangleCentralEdgeSize{0.25 / std::tan(M_PI / 6)}; + const double desiredBaseTriangleCentralEdgeSize{0.025}; }; #endif // FLATPATTERN_HPP diff --git a/simulationresult.hpp b/simulationresult.hpp index 33c990c..9d2c578 100644 --- a/simulationresult.hpp +++ b/simulationresult.hpp @@ -97,6 +97,26 @@ struct SimulationJob { ->setEnabled(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; + } }; struct SimulationResults { @@ -163,33 +183,32 @@ struct SimulationResults { }); // per node external forces std::vector> externalForces(job.mesh->VN()); + std::vector> externalMoments(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]}; + externalMoments[index] = {force[3], force[4], 0}; } polyscope::getCurveNetwork(undeformedMeshName) - ->addNodeColorQuantity("Boundary conditions", nodeColors); - polyscope::getCurveNetwork(undeformedMeshName) - ->getQuantity("Boundary conditions") + ->addNodeColorQuantity("Boundary conditions", nodeColors) ->setEnabled(true); polyscope::getCurveNetwork(undeformedMeshName) - ->addNodeVectorQuantity("External force", externalForces); - polyscope::getCurveNetwork(undeformedMeshName) - ->getQuantity("External force") + ->addNodeVectorQuantity("External force", externalForces) ->setEnabled(true); + + polyscope::getCurveNetwork(undeformedMeshName) + ->addNodeVectorQuantity("External moments", externalMoments) + ->setEnabled(true); + polyscope::getCurveNetwork(deformedMeshName) - ->addNodeColorQuantity("Boundary conditions", nodeColors); - polyscope::getCurveNetwork(deformedMeshName) - ->getQuantity("Boundary conditions") + ->addNodeColorQuantity("Boundary conditions", nodeColors) ->setEnabled(false); polyscope::getCurveNetwork(deformedMeshName) - ->addNodeVectorQuantity("External force", externalForces); - polyscope::getCurveNetwork(deformedMeshName) - ->getQuantity("External force") + ->addNodeVectorQuantity("External force", externalForces) ->setEnabled(true); // polyscope::show();