From ef18646cfd6ad89e3ae8fb19181611683f386d91 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 19 Jul 2022 14:54:39 +0300 Subject: [PATCH] Refactoring --- der_leimer.cpp | 26 +++- der_leimer.hpp | 11 +- edgemesh.cpp | 2 +- formFinder_unitTestFiles/simpleBeam.ply | 3 +- simulation_structs.hpp | 195 +++++++++++++++--------- simulationmesh.cpp | 6 +- 6 files changed, 159 insertions(+), 84 deletions(-) diff --git a/der_leimer.cpp b/der_leimer.cpp index c43d606..e2da1fd 100644 --- a/der_leimer.cpp +++ b/der_leimer.cpp @@ -427,8 +427,8 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptrpMesh->getBeamMaterial()[0].youngsModulus; - const double w = pJob->pMesh->getBeamDimensions()[0].b; - const double h = pJob->pMesh->getBeamDimensions()[0].h; + const double w = pJob->pMesh->getBeamDimensions()[0].getDim1(); + const double h = pJob->pMesh->getBeamDimensions()[0].getDim2(); rodParams.thickness = h; rodParams.kbending = E; rodParams.kstretching = E; @@ -455,6 +455,13 @@ SimulationResults DER_leimer::executeSimulation(const std::shared_ptrnumIntersections() == 0); + + // pRod = std::make_unique(verts); + // const auto cs = CrossSection::construct("rectangle", E, poissonsRatio, {w, h}); + // RodMaterial mat(*cs); + // pRod->setMaterial(mat); + bool done = false; while (!done) { done = simulateOneStepGrid(config); @@ -567,6 +574,19 @@ bool DER_leimer::simulateOneStepGrid(RodConfig *config) forceResidual = lineSearchGrid(*config, updateVec, params, externalForces); std::cout << "Force residual: " << forceResidual << std::endl; + //panetta's rod update + // assert(config->numIntersections() == 0); + std::vector vertexPositions(config->numVertices()); + for (int vi = 0; vi < config->numVertices(); vi++) { + vertexPositions[vi] = config->vertices.row(vi); + } + std::vector thetas(config->numRods()); + for (int ri = 0; ri < config->numRods(); ri++) { + assert(config->rods[ri]->numSegments() == 1); + thetas[ri] = config->rods[ri]->curState.thetas[0]; + } + // pRod->setDeformedConfiguration(vertexPositions, thetas); + rAndJGrid(*config, r, &Jr, linE, gravityForces, params, externalForcesPE, externalForces); iter++; @@ -4346,7 +4366,7 @@ DER_leimer::RodConfig *DER_leimer::readRodGrid( } Eigen::VectorXd segmentWidths_eigen(numSegments_rod); - const double rodWidth = pMesh->elements[rodIndex].dimensions.b; + const double rodWidth = pMesh->elements[rodIndex].dimensions.getDim1(); for (int si_rod = 0; si_rod < numSegments_rod; si_rod++) { segmentWidths_eigen[si_rod] = rodWidth; } diff --git a/der_leimer.hpp b/der_leimer.hpp index 015c383..d983d9c 100644 --- a/der_leimer.hpp +++ b/der_leimer.hpp @@ -1,6 +1,7 @@ #ifndef DER_LEIMER_HPP #define DER_LEIMER_HPP +//#include "ElasticRod.hh" #include "simulationmodel.hpp" #include @@ -8,7 +9,7 @@ class DER_leimer : public SimulationModel { public: DER_leimer(); - virtual SimulationResults executeSimulation(const std::shared_ptr &pJob); + virtual SimulationResults executeSimulation(const std::shared_ptr &pJob) override; struct RodParams { double thickness; @@ -175,15 +176,17 @@ private: const Eigen::VectorXd &update, const SimParams ¶ms, const Eigen::VectorXd &externalForces = Eigen::VectorXd()); + SimulationResults constructSimulationResults(const std::shared_ptr &pJob, + RodConfig *config) const; + Eigen::Vector3d getPerpendicularVector(const Eigen::Vector3d &t) const; double energy; Eigen::VectorXd updateVec; double forceResidual; int iter{0}; - SimulationResults constructSimulationResults(const std::shared_ptr &pJob, - RodConfig *config) const; std::shared_ptr pJob; - Eigen::Vector3d getPerpendicularVector(const Eigen::Vector3d &t) const; + + // std::unique_ptr pRod; }; Eigen::Matrix3d eulerRotation(double e0, double e1, double e2, int deriv = -1, int deriv2 = -1); diff --git a/edgemesh.cpp b/edgemesh.cpp index cb5c36b..9346701 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -202,7 +202,7 @@ bool VCGEdgeMesh::createSpanGrid(const size_t desiredWidth, return true; } -bool VCGEdgeMesh::load(const std::filesystem::__cxx11::path &meshFilePath) +bool VCGEdgeMesh::load(const std::filesystem::path &meshFilePath) { std::string usedPath = meshFilePath; if (std::filesystem::path(meshFilePath).is_relative()) { diff --git a/formFinder_unitTestFiles/simpleBeam.ply b/formFinder_unitTestFiles/simpleBeam.ply index 1d3d728..5681759 100644 --- a/formFinder_unitTestFiles/simpleBeam.ply +++ b/formFinder_unitTestFiles/simpleBeam.ply @@ -7,7 +7,6 @@ property float z { z coordinate, too } property float nx property float ny property float nz - element edge 4 { there are 6 "face" elements in the file } property int vertex1 property int vertex2 @@ -24,4 +23,4 @@ end_header 0 1 2 0.03 0.026 2 0.3 200 1 2 2 0.03 0.026 2 0.3 200 2 3 2 0.03 0.026 2 0.3 200 -3 4 2 0.03 0.026 2 0.3 200 \ No newline at end of file +3 4 2 0.03 0.026 2 0.3 200 diff --git a/simulation_structs.hpp b/simulation_structs.hpp index aa6cbfb..b27bc3e 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -8,6 +8,9 @@ #include #include +#ifdef POLYSCOPE_DEFINED +#include +#endif namespace Eigen { template void write_binary(const std::string &filename, const Matrix &matrix) { @@ -157,28 +160,35 @@ class SimulationJob { // simulation mesh class } jsonLabels; -public: - inline static std::string jsonDefaultFileName = "SimulationJob.json"; - std::shared_ptr pMesh; - std::string label{"empty_job"}; - std::unordered_map> constrainedVertices; - std::unordered_map nodalExternalForces; - std::unordered_map nodalForcedDisplacements; - SimulationJob(const std::shared_ptr &m, - const std::string &label, - const std::unordered_map> &cv, - const std::unordered_map &ef = {}, - const std::unordered_map &fd = {}) - : pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef), - nodalForcedDisplacements(fd) - {} - SimulationJob() {} - SimulationJob(const std::string &jsonFilename) { load(jsonFilename); } +#ifdef POLYSCOPE_DEFINED + const std::string polyscopeLabel_bcAsPointCloud{"BC_spheres"}; +#endif - bool isEmpty() - { - return constrainedVertices.empty() && nodalExternalForces.empty() - && nodalForcedDisplacements.empty() && pMesh == nullptr; + public: + inline static std::string jsonDefaultFileName = "SimulationJob.json"; + std::shared_ptr pMesh; + std::string label{"empty_job"}; + std::unordered_map> constrainedVertices; + std::unordered_map nodalExternalForces; + std::unordered_map nodalForcedDisplacements; + SimulationJob(const std::shared_ptr &m, + const std::string &label, + const std::unordered_map> &cv, + const std::unordered_map &ef = {}, + const std::unordered_map &fd = {}) + : pMesh(m), label(label), constrainedVertices(cv), nodalExternalForces(ef), + nodalForcedDisplacements(fd) + {} + SimulationJob() {} + SimulationJob(const std::string &jsonFilename) + { + load(jsonFilename); + } + + bool isEmpty() + { + return constrainedVertices.empty() && nodalExternalForces.empty() + && nodalForcedDisplacements.empty() && pMesh == nullptr; } void remap(const std::unordered_map &sourceToDestinationViMap, SimulationJob &destination_simulationJob) const @@ -470,6 +480,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m ->addNodeColorQuantity("Boundary conditions", nodeColors) ->setEnabled(shouldEnable); } + drawBcAsSpheres(); // per node external forces std::vector> externalForces(pMesh->VN()); @@ -509,7 +520,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m ->setEnabled(shouldEnable); } } - void unregister(const std::string &meshLabel) + void unregister(const std::string &meshLabel) const { if (polyscope::getCurveNetwork(meshLabel) == nullptr) { return; @@ -519,6 +530,7 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m } if (!constrainedVertices.empty()) { polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions"); + polyscope::removeStructure(polyscopeLabel_bcAsPointCloud, false); } // per node external moments @@ -534,6 +546,38 @@ json[jsonLabels.meshFilename]= std::filesystem::relative(std::filesystem::path(m polyscope::getCurveNetwork(meshLabel)->removeQuantity("External moment"); } } + + void drawBcAsSpheres() const + { + polyscope::removeStructure(polyscopeLabel_bcAsPointCloud, false); + std::vector bcPos; + std::vector bcColors; + for (std::pair> bc : constrainedVertices) { + bcPos.push_back(glm::vec3(pMesh->vert[bc.first].cP()[0], + pMesh->vert[bc.first].cP()[1], + pMesh->vert[bc.first].cP()[2])); + const bool hasRotationalDoFConstrained = bc.second.contains(3) || bc.second.contains(4) + || bc.second.contains(5); + const bool hasTranslationalDoFConstrained = bc.second.contains(0) || bc.second.contains(1) + || bc.second.contains(2); + + std::array color_fixedTranslation{63.0 / 255, 191.0 / 255, 127.0 / 255}; + std::array color_rigid({1, 0, 0}); + if (hasTranslationalDoFConstrained && !hasRotationalDoFConstrained) { + bcColors.push_back(glm::vec3(color_fixedTranslation[0], + color_fixedTranslation[1], + color_fixedTranslation[2])); + } else if (!hasTranslationalDoFConstrained && hasRotationalDoFConstrained) { + bcColors.push_back(glm::vec3(0, 0, 1)); + } else { + bcColors.push_back(glm::vec3(color_rigid[0], color_rigid[1], color_rigid[2])); + } + } + + auto bcPolyscopePointCloud = polyscope::registerPointCloud(polyscopeLabel_bcAsPointCloud, + bcPos); + bcPolyscopePointCloud->addColorQuantity("bc_colors", bcColors)->setEnabled(true); + } #endif // POLYSCOPE_DEFINED }; struct SimulationResults @@ -554,6 +598,7 @@ struct SimulationResults double internalPotentialEnergy{0}; double executionTime{0}; std::vector> perVertexInternalForces; //axial,torsion,bending1,bending2 + std::string simulationModelUsed{"empty label"}; std::string labelPrefix{"deformed"}; inline static char deliminator{' '}; SimulationResults() { pJob = std::make_shared(); } @@ -572,7 +617,8 @@ struct SimulationResults void setLabelPrefix(const std::string &lp) { labelPrefix += deliminator + lp; } std::string getLabel() const { - return labelPrefix + deliminator + pJob->pMesh->getLabel() + deliminator + pJob->getLabel(); + return labelPrefix + deliminator + simulationModelUsed + deliminator + + pJob->pMesh->getLabel() + deliminator + pJob->getLabel(); } bool saveDeformedModel(const std::string &outputFolder = std::string()) @@ -742,46 +788,51 @@ struct SimulationResults pJob->unregister(getLabel()); polyscope::removeCurveNetwork(getLabel()); } + polyscope::CurveNetwork *registerForDrawing( const std::optional> &desiredColor = std::nullopt, - const bool &shouldEnable = true, - const double &desiredRadius = 0.001, - // const double &desiredRadius = 0.0001, - const bool &shouldShowFrames = false) const + const bool &shouldEnable = true /*, + const bool &shouldShowFrames = false*/) const { + if (!converged) { + std::cerr << "Simulation did not converge. Drawing nothing." << std::endl; + return nullptr; + } PolyscopeInterface::init(); - const std::shared_ptr &mesh = pJob->pMesh; + const std::shared_ptr &simulationMesh = pJob->pMesh; polyscope::CurveNetwork *polyscopeHandle_deformedEdmeMesh; - if (!polyscope::hasStructure(getLabel())) { - const auto verts = mesh->getEigenVertices(); - const auto edges = mesh->getEigenEdges(); - polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(getLabel(), + const std::string &meshLabel = getLabel(); + if (!polyscope::hasStructure(meshLabel)) { + const auto &verts = simulationMesh->getEigenVertices(); + const auto &edges = simulationMesh->getEigenEdges(); + polyscopeHandle_deformedEdmeMesh = polyscope::registerCurveNetwork(meshLabel, verts, edges); } else { - polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(getLabel()); + polyscopeHandle_deformedEdmeMesh = polyscope::getCurveNetwork(meshLabel); } polyscopeHandle_deformedEdmeMesh->setEnabled(shouldEnable); - polyscopeHandle_deformedEdmeMesh->setRadius(desiredRadius, true); + const double drawingRadius = simulationMesh->getBeamDimensions()[0].getDrawingRadius(); + polyscopeHandle_deformedEdmeMesh->setRadius(drawingRadius /*0.00025*/, false); if (desiredColor.has_value()) { const glm::vec3 desiredColor_glm(desiredColor.value()[0], desiredColor.value()[1], desiredColor.value()[2]); - polyscopeHandle_deformedEdmeMesh->setColor(/*glm::normalize(*/ desiredColor_glm /*)*/); + polyscopeHandle_deformedEdmeMesh->setColor(desiredColor_glm); } - Eigen::MatrixX3d nodalDisplacements(mesh->VN(), 3); - Eigen::MatrixX3d framesX(mesh->VN(), 3); - Eigen::MatrixX3d framesY(mesh->VN(), 3); - Eigen::MatrixX3d framesZ(mesh->VN(), 3); + Eigen::MatrixX3d nodalDisplacements(simulationMesh->VN(), 3); + Eigen::MatrixX3d framesX(simulationMesh->VN(), 3); + Eigen::MatrixX3d framesY(simulationMesh->VN(), 3); + Eigen::MatrixX3d framesZ(simulationMesh->VN(), 3); - Eigen::MatrixX3d framesX_initial(mesh->VN(), 3); - Eigen::MatrixX3d framesY_initial(mesh->VN(), 3); - Eigen::MatrixX3d framesZ_initial(mesh->VN(), 3); + Eigen::MatrixX3d framesX_initial(simulationMesh->VN(), 3); + Eigen::MatrixX3d framesY_initial(simulationMesh->VN(), 3); + Eigen::MatrixX3d framesZ_initial(simulationMesh->VN(), 3); - // std::unordered_set interfaceNodes{1, 3, 5, 7, 9, 11}; - // std::unordered_set interfaceNodes{3, 7, 11, 15, 19, 23}; + // std::unordered_set interfaceNodes{1, 3, 5, 7, 9, 11}; + std::unordered_set interfaceNodes{3, 7, 11, 15, 19, 23}; // std::unordered_set interfaceNodes{}; - for (VertexIndex vi = 0; vi < mesh->VN(); vi++) { + for (VertexIndex vi = 0; vi < simulationMesh->VN(); vi++) { const Vector6d &nodalDisplacement = displacements[vi]; nodalDisplacements.row(vi) = Eigen::Vector3d(nodalDisplacement[0], nodalDisplacement[1], @@ -790,32 +841,36 @@ struct SimulationResults // Eigen::Quaternion Ry(Eigen::AngleAxis(nodalDisplacement[4],Eigen::Vector3d(0, 1, 0))); // Eigen::Quaternion Rz(Eigen::AngleAxis(nodalDisplacement[5],Eigen::Vector3d(0, 0, 1))); // Eigen::Quaternion R=Rx*Ry*Rz; - // if (interfaceNodes.contains(vi)) { - auto deformedNormal = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 0, 1); - auto deformedFrameY = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(0, 1, 0); - auto deformedFrameX = rotationalDisplacementQuaternion[vi] * Eigen::Vector3d(1, 0, 0); - framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], - deformedFrameX[1], - deformedFrameX[2]); - framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], - deformedFrameY[1], - deformedFrameY[2]); - framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], - deformedNormal[1], - deformedNormal[2]); - framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0); - framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0); - framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1); - // } else { - // framesX.row(vi) = Eigen::Vector3d(0, 0, 0); - // framesY.row(vi) = Eigen::Vector3d(0, 0, 0); - // framesZ.row(vi) = Eigen::Vector3d(0, 0, 0); - // framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - // framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - // framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0); - // } + if (interfaceNodes.contains(vi)) { + assert(!rotationalDisplacementQuaternion.empty()); + auto deformedNormal = rotationalDisplacementQuaternion[vi] + * Eigen::Vector3d(0, 0, 1); + auto deformedFrameY = rotationalDisplacementQuaternion[vi] + * Eigen::Vector3d(0, 1, 0); + auto deformedFrameX = rotationalDisplacementQuaternion[vi] + * Eigen::Vector3d(1, 0, 0); + framesX.row(vi) = Eigen::Vector3d(deformedFrameX[0], + deformedFrameX[1], + deformedFrameX[2]); + framesY.row(vi) = Eigen::Vector3d(deformedFrameY[0], + deformedFrameY[1], + deformedFrameY[2]); + framesZ.row(vi) = Eigen::Vector3d(deformedNormal[0], + deformedNormal[1], + deformedNormal[2]); + framesX_initial.row(vi) = Eigen::Vector3d(1, 0, 0); + framesY_initial.row(vi) = Eigen::Vector3d(0, 1, 0); + framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 1); + } else { + framesX.row(vi) = Eigen::Vector3d(0, 0, 0); + framesY.row(vi) = Eigen::Vector3d(0, 0, 0); + framesZ.row(vi) = Eigen::Vector3d(0, 0, 0); + framesX_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + framesY_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + framesZ_initial.row(vi) = Eigen::Vector3d(0, 0, 0); + } } - polyscopeHandle_deformedEdmeMesh->updateNodePositions(mesh->getEigenVertices() + polyscopeHandle_deformedEdmeMesh->updateNodePositions(simulationMesh->getEigenVertices() + nodalDisplacements); const double frameRadius_default = 0.035; diff --git a/simulationmesh.cpp b/simulationmesh.cpp index b7f95d3..6be0423 100755 --- a/simulationmesh.cpp +++ b/simulationmesh.cpp @@ -257,9 +257,7 @@ void SimulationMesh::updateElementalFrames() { polyscope::CurveNetwork *SimulationMesh::registerForDrawing( const std::optional> &desiredColor, const bool &shouldEnable) { - const double drawingRadius = getBeamDimensions()[0].b - / (std::sqrt( - 2.0)); //polyscope can only draw a circular cross section + const double drawingRadius = getBeamDimensions()[0].getDrawingRadius(); // std::cout << __FUNCTION__ << " revert this" << std::endl; return VCGEdgeMesh::registerForDrawing(desiredColor, /*0.08*/ drawingRadius, shouldEnable); } @@ -500,7 +498,7 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, void Element::setDimensions(const CrossSectionType &dimensions) { this->dimensions = dimensions; assert(this->dimensions.A == dimensions.A); - // computeDimensionsProperties(dimensions); + // computeDimensionsProperties(dimensions); updateRigidity(); }