From 5d8445fbca10b54ce908ccc5fa43a7126be5e362 Mon Sep 17 00:00:00 2001 From: Iason Date: Wed, 9 Dec 2020 17:59:18 +0200 Subject: [PATCH] Refactoring --- beam.hpp | 26 +---- beamformfinder.cpp | 262 ++++++++++++++++--------------------------- beamformfinder.hpp | 9 +- elementalmesh.cpp | 41 ++++--- elementalmesh.hpp | 8 +- flatpattern.hpp | 2 +- simulationresult.hpp | 1 + utilities.hpp | 8 +- 8 files changed, 147 insertions(+), 210 deletions(-) diff --git a/beam.hpp b/beam.hpp index 463261e..4851276 100644 --- a/beam.hpp +++ b/beam.hpp @@ -12,21 +12,21 @@ struct RectangularBeamDimensions { : b(width), h(height) { assert(width > 0 && height > 0); } - RectangularBeamDimensions() : b(0.01), h(0.01) {} + RectangularBeamDimensions() : b(1), h(1) {} }; -struct CylindricalElementDimensions { +struct CylindricalBeamDimensions { float od; // Cylinder outside diameter float id; // Cylinder inside diameter // https://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html - CylindricalElementDimensions(const float &outsideDiameter, - const float &insideDiameter) + CylindricalBeamDimensions(const float &outsideDiameter, + const float &insideDiameter) : od(outsideDiameter), id(insideDiameter) { assert(outsideDiameter > 0 && insideDiameter > 0 && outsideDiameter > insideDiameter); } - CylindricalElementDimensions() : od(0.03), id(0.026) {} + CylindricalBeamDimensions() : od(0.03), id(0.026) {} }; struct ElementMaterial { @@ -37,21 +37,7 @@ struct ElementMaterial { youngsModulusGPascal(youngsModulusGPascal) { assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); } - ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(7.5) {} -}; - -struct BeamProperties { - float crossArea; - float I2; - float I3; - float polarInertia; - BeamProperties(const RectangularBeamDimensions &dimensions, - const ElementMaterial &material) { - crossArea = (dimensions.b * dimensions.h); - I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; - I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; - polarInertia = (I2 + I3); - } + ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(200) {} }; #endif // BEAM_HPP diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 794d43a..1d39e20 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -51,7 +51,7 @@ VectorType FormFinder::computeDerivativeOfNormal( const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); if (dui.dofi == 3) { - if (nxnyMagnitude >= 1) { + if (nxnyMagnitude + 1e-5 >= 1) { const double normalDerivativeX = 1 / sqrt(nxnyMagnitude) - std::pow(nx, 2) / std::pow(nxnyMagnitude, 1.5); @@ -67,7 +67,7 @@ VectorType FormFinder::computeDerivativeOfNormal( VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); } } else if (dui.dofi == 4) { - if (nxnyMagnitude >= 1) { + if (nxnyMagnitude + 1e-5 >= 1) { const double normalDerivativeX = -nx * ny / std::pow(nxnyMagnitude, 1.5); const double normalDerivativeY = 1 / sqrt(nxnyMagnitude) - @@ -88,6 +88,8 @@ VectorType FormFinder::computeDerivativeOfNormal( std::isfinite(normalDerivative[1]) && std::isfinite(normalDerivative[2]); assert(normalDerivativeIsFinite); + const bool shouldBreak = + mCurrentSimulationStep == 118 && vi == 1 && dui.dofi == 3; return normalDerivative; } @@ -223,6 +225,8 @@ FormFinder::computeDerivativeT2(const EdgeType &e, const double t2DerivNorm = t2Deriv.Norm(); assert(std::isfinite(t2DerivNorm)); + const bool shouldBreak = mCurrentSimulationStep == 118 && + (vi_jplus1 == 1 || vi_j == 1) && dofi == 3; return t2Deriv; } @@ -297,6 +301,8 @@ double FormFinder::computeDerivativeTheta1(const EdgeType &e, const double theta1Derivative = derivativeT1 * Cross(t3, n) + t1 * (Cross(derivativeT3, n) + Cross(t3, nDerivative)); + const bool shouldBreak = + mCurrentSimulationStep == 118 && vi == 1 && dwrt_dofi == 3; return theta1Derivative; } @@ -703,6 +709,8 @@ void FormFinder::updateResidualForcesOnTheFly( secondBendingForce_inBracketsTerm_3; double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * element.secondBendingConstFactor; + const bool shouldBreak = + mCurrentSimulationStep == 118 && edgeNode.vi == 1 && dofi == 3; internalForcesContributionFromThisEdge[evi].second[dofi] = axialForce_dofi + firstBendingForce_dofi + secondBendingForce_dofi + torsionalForce_dofi; @@ -1024,7 +1032,7 @@ void FormFinder::updateNodePosition( node.previousLocation = previousLocation; v.P() = node.initialLocation + displacementVector; if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) { - VectorType desiredInitialDisplacement(0, 0, 0.1); + VectorType desiredInitialDisplacement(0, 0, 0.001); v.P() = v.P() + desiredInitialDisplacement; } } @@ -1033,8 +1041,14 @@ void FormFinder::updateNodeNormal( VertexType &v, const std::unordered_map> &fixedVertices) { + Node &node = mesh->nodes[v]; const VertexIndex &vi = node.vi; + // if (vi == 1) { + // std::cout << "PRE:" << mesh->vert[1].N()[0] << " " << + // mesh->vert[1].N()[1] + // << " " << mesh->vert[1].N()[2] << std::endl; + // } VectorType normalDisplacementVector(0, 0, 0); if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(3)) { normalDisplacementVector += VectorType(node.displacements[3], 0, 0); @@ -1042,9 +1056,12 @@ void FormFinder::updateNodeNormal( if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(4)) { normalDisplacementVector += VectorType(0, node.displacements[4], 0); } + const double nxnyMagnitudePre = std::pow(v.N()[0], 2) + std::pow(v.N()[1], 2); v.N() = node.initialNormal + normalDisplacementVector; - const double &nx = v.N()[0], ny = v.N()[1]; + 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; if (nxnyMagnitude > 1) { VectorType newNormal(nx / std::sqrt(nxnyMagnitude), ny / std::sqrt(nxnyMagnitude), 0); @@ -1067,6 +1084,11 @@ void FormFinder::updateNodeNormal( VectorType g1 = Cross(refT1, t1Initial); node.nR = g1.dot(v.cN()); } + // if (vi == 1) { + // std::cout << "AFTER:" << mesh->vert[1].N()[0] << " " << + // mesh->vert[1].N()[1] + // << " " << mesh->vert[1].N()[2] << std::endl; + // } } void FormFinder::applyDisplacements( @@ -1103,6 +1125,25 @@ void FormFinder::applyForcedDisplacements( VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1), vertexDisplacement(2)); mesh->vert[vi].P() = node.initialLocation + displacementVector; + node.displacements = Vector6d( + {vertexDisplacement(0), vertexDisplacement(1), vertexDisplacement(2), + node.displacements[3], node.displacements[4], node.displacements[5]}); + } +} + +void FormFinder::applyForcedNormals( + const std::unordered_map nodalForcedRotations) { + for (const std::pair vertexIndexDesiredNormalPair : + nodalForcedRotations) { + const VertexIndex vi = vertexIndexDesiredNormalPair.first; + + Node &node = mesh->nodes[vi]; + mesh->vert[vi].N() = vertexIndexDesiredNormalPair.second; + node.displacements = Vector6d( + {node.displacements[0], node.displacements[1], node.displacements[2], + vertexIndexDesiredNormalPair.second[0] - node.initialNormal[0], + vertexIndexDesiredNormalPair.second[1] - node.initialNormal[1], + node.displacements[5]}); } } @@ -1547,12 +1588,14 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, constrainedVertices[vi].insert({0, 1, 2}); } } + if (!job.nodalForcedNormals.empty()) { + for (std::pair viNormalPair : + job.nodalForcedDisplacements) { + assert(viNormalPair.second[2] >= 0); + } + } - // VCGEdgeMesh *jobEdgeMesh = job.mesh.get(); mesh = std::make_unique(*job.mesh); - // ElementalMesh *jobElementalMesh = job.mesh.get(); - // vcg::tri::Append::MeshCopy( - // *(this->mesh), *jobElementalMesh); if (beVerbose) { std::cout << "Executing simulation for mesh with:" << mesh->VN() << " nodes and " << mesh->EN() << " elements." << std::endl; @@ -1564,75 +1607,21 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, if (!polyscope::state::initialized) { initPolyscope(); } - if (!polyscope::hasCurveNetwork(deformedMeshName)) { - polyscope::registerCurveNetwork( - deformedMeshName, mesh->getEigenVertices(), mesh->getEigenEdges()); - } - + polyscope::registerCurveNetwork(deformedMeshName, mesh->getEigenVertices(), + mesh->getEigenEdges()); // registerWorldAxes(); } - // const double beamRadius = mesh.getBeamDimensions()[0].od / 2; - // std::cout << "BeamRadius:" << beamRadius << std::endl; - // polyscope::getCurveNetwork("Deformed edge mesh") - // // - // // ->setRadius(beamRadius); - // ->setRadius(0.01); - - // const bool hasExternalLoads = - // !job.nodalExternalForces.empty() || - // !job.nodalForcedDisplacements.empty(); - // assert(hasExternalLoads); for (auto fixedVertex : job.fixedVertices) { assert(fixedVertex.first < mesh->VN()); } updateElementalFrames(); updateNodalMasses(); - if (job.nodalExternalForces.empty()) { + if (!job.nodalForcedDisplacements.empty() && + job.nodalExternalForces.empty()) { shouldApplyInitialDistortion = true; } - // std::unordered_map temporaryForces{ - // // // {2, Eigen::Vector3d(0, 0, 1000)}, - // // // {12, Eigen::Vector3d(0, 0, 1000)}, - // // // {18, Eigen::Vector3d(0, 0, 1000)}, - // // // {8, Eigen::Vector3d(0, 0, 1000)}}; - // {10, Eigen::Vector3d(0, 0, 10000)}}; - // std::unordered_map temporaryForces; - // for (VertexIndex vi = 0; vi < mesh.VN(); vi++) { - // for (VertexType &v : mesh.vert) { - // const VertexIndex vi = mesh.getIndex(v); - // VectorType allowedDoFType(1, 1, 1); - // if (constrainedVertices.contains(vi)) { - // std::unordered_set constrainedDof = - // constrainedVertices.at(vi); - // if (constrainedDof.contains(0)) { - // allowedDoFType[0] = 0; - // } else if (constrainedDof.contains(1)) { - // allowedDoFType[1] = 0; - // } else if (constrainedDof.contains(2)) { - // allowedDoFType[2] = 0; - // } - // } - // VectorType desiredDisplacement = VectorType(0, 0, 0.2); - // VectorType displacement(allowedDoFType[0] * desiredDisplacement[0], - // allowedDoFType[1] * desiredDisplacement[1], - // allowedDoFType[2] * desiredDisplacement[2]); - // v.P() = v.P() + displacement; - // temporaryForces.insert({vi, Eigen::Vector3d(0, 0, 100000)}); - // } - // updateNodalExternalForces(temporaryForces, job.fixedVertices); - // appliedTemporaryForce = true; - // } else { - // std::unordered_map appliedLoad = - // job.nodalExternalForces; - // const size_t numberOfStepsForApplyingExternalLoads = 10; - // int externalLoadStep = 1; - // std::for_each(appliedLoad.begin(), appliedLoad.end(), [](auto &pair) { - // pair.second /= numberOfStepsForApplyingExternalLoads; - // }); - // updateNodalExternalForces(appliedLoad, constrainedVertices); updateNodalExternalForces(job.nodalExternalForces, constrainedVertices); - // } while (mCurrentSimulationStep < maxDRMIterations) { // while (true) { @@ -1641,36 +1630,7 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, updateRDerivatives(); updateT2Derivatives(); updateT3Derivatives(); - // updateRotationalDisplacements(); - // if (currentSimulationStep > lastPulseStepIndex) { - // const std::vector internalForces = updateResidualForcesOnTheFly(constrainedVertices); - //#pragma omp parallel for schedule(static) num_threads(8) - // double totalResidualForcesNorm = 0; - // for (size_t vi = 0; vi < mesh.VN(); vi++) { - // Node::Forces &force = mesh.nodes[vi].force; - // // const Vector6d ¶llelForce = internalForcesParallel[vi]; - // // const Vector6d &serialForce = internalForces[vi]; - // // const double normOfDifference = (serialForce + - // (parallelForce - // * - // // -1)).norm(); assert(normOfDifference < 0.000001); - // force.residual = force.external - internalForces[vi]; - // const double residualForceNorm = (force.residual).norm(); - // totalResidualForcesNorm += residualForceNorm; - // // assert(residualForceNorm < std::pow(10, 20)); - // } - // mesh.totalResidualForcesNorm = totalResidualForcesNorm; - // } else { - // if (currentSimulationStep == lastPulseStepIndex && - // appliedTemporaryForce) { - // for (const VertexType &v : mesh.vert) { - // Node &node = mesh.nodes[v]; - // node.force.external = std::optional(); - // } - // } - // updateNodalInternalForces(job.fixedVertices); - // } // TODO: write parallel function for updating positions // TODO: Make a single function out of updateResidualForcesOnTheFly @@ -1686,16 +1646,11 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose, job.nodalForcedDisplacements); } + if (!job.nodalForcedNormals.empty()) { + applyForcedNormals(job.nodalForcedNormals); + } updateElementalLengths(); - // if (!std::isfinite(mesh.currentTotalKineticEnergy)) { - // std::cerr << "Infinite kinetic energy. Breaking simulation.." - // << std::endl; - // break; - // } - // assert(std::isfinite(mesh.currentTotalKineticEnergy)); - // mesh.printVertexCoordinates(mesh.VN() / 2); - // draw(); if (mShouldDraw && mCurrentSimulationStep % mDrawingStep == 0 /* && currentSimulationStep > maxDRMIterations*/) { @@ -1704,22 +1659,8 @@ currentSimulationStep > maxDRMIterations*/) { // .append("Screenshots") // .string(); // draw(saveTo); - // SimulationResultsReporter::createPlot( - // "Number of iterations", "Log of Kinetic Energy", - // history.kineticEnergy, - // std::filesystem::path(saveTo).append( - // std::to_string(currentSimulationStep) + "_kin.png")); - - // draw(); - // auto t2 = std::chrono::high_resolution_clock::now(); - // auto timePerNodePerIteration = - // std::chrono::duration_cast(t2 - - // t1) - // .count() / - // (mesh.VN() * (currentSimulationStep + 1)); - // std::cout << "Time/(node*iteration) " - // << timePerNodePerIteration * std::pow(10, -6) << "s" - // << std::endl; + std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm + << std::endl; draw(); } if (mCurrentSimulationStep != 0) { @@ -1745,45 +1686,11 @@ currentSimulationStep > maxDRMIterations*/) { std::cout << "Total potential:" << mesh->totalPotentialEnergykN << " kNm" << std::endl; } - // if (externalLoadStep != - // numberOfStepsForApplyingExternalLoads) - // { - // std::for_each( - // appliedLoad.begin(), appliedLoad.end(), [&](auto - // &pair) - // { - // pair.second += - // job.nodalExternalForces.at(pair.first) - // / - // numberOfStepsForApplyingExternalLoads; - // }); - // updateNodalExternalForces(appliedLoad, - // constrainedVertices); - // } else { break; // } } - // history.markRed(currentSimulationStep); - // std::cout << "Total potential:" << mesh.totalPotentialEnergykN - // << " kNm" - // << std::endl; - // reset displacements to previous position where the peak occured - // for (VertexType &v : mesh->vert) { - // Node &node = mesh->nodes[v]; - // node.displacements = node.displacements - node.velocity * Dt; - // } - // applyDisplacements(constrainedVertices, shouldDraw); - // if (!job.nodalForcedDisplacements.empty()) { - // applyForcedDisplacements( - - // job.nodalForcedDisplacements); - // } - // updateElementalLengths(); resetVelocities(); Dt = Dt * xi; - // std::cout << "Residual forces norm:" << - // mesh.totalResidualForcesNorm - // << std::endl; } } if (mCurrentSimulationStep == maxDRMIterations) { @@ -1809,19 +1716,42 @@ currentSimulationStep > maxDRMIterations*/) { polyscope::removeCurveNetwork(deformedMeshName); // } } - // debugger.createPlots(); SimulationResults results = computeResults(*job.mesh); - // Eigen::write_binary("optimizedResult.eigenBin", - // results.nodalDisplacements); - - // const std::string groundOfTruthBinaryFilename = - // "../../groundOfTruths/grid6_groundOfTruth.eigenBin"; - // // "../../groundOfTruths/longSpanGridshell_groundOfTruth.eigenBin"; - // assert(std::filesystem::exists( - // std::filesystem::path(groundOfTruthBinaryFilename))); - // Eigen::MatrixX3d groundOfTruthMatrix; - // Eigen::read_binary(groundOfTruthBinaryFilename, groundOfTruthMatrix); - // assert(results.nodalDisplacements.isApprox(groundOfTruthMatrix)); return results; - // return history; +} + +void FormFinder::runUnitTests() { + FormFinder formFinder; + VCGEdgeMesh mesh; + // const size_t spanGridSize = 11; + // mesh.createSpanGrid(spanGridSize); + mesh.loadFromPly("/home/iason/Models/simple_beam_paper_example.ply"); + std::unordered_map> fixedVertices; + fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; + fixedVertices[mesh.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)}); + + SimulationJob beamSimulationJob{ + std::make_shared(mesh), + // 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); + + 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(); } diff --git a/beamformfinder.hpp b/beamformfinder.hpp index 30cffc0..2e6c08f 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -27,7 +27,7 @@ private: const double Dtini{0.1}; double Dt{Dtini}; const double xi{0.9969}; - const double totalResidualForcesNormThreshold{10}; + const double totalResidualForcesNormThreshold{0.01}; size_t mCurrentSimulationStep{0}; int mDrawingStep{5}; bool mShouldDraw{false}; @@ -176,13 +176,18 @@ private: const std::unordered_map> &fixedVertices); + void applyForcedNormals( + const std::unordered_map nodalForcedRotations); + public: FormFinder(); SimulationResults executeSimulation( const SimulationJob &job, const bool &beVerbose = false, const bool &shouldDraw = false, const size_t &maxDRMIterations = FormFinder::maxDRMIterations); - inline static const size_t maxDRMIterations{50000}; + inline static const size_t maxDRMIterations{500000}; + + static void runUnitTests(); }; template PointType Cross(PointType p1, PointType p2) { diff --git a/elementalmesh.cpp b/elementalmesh.cpp index 752aae7..c049b96 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -2,6 +2,21 @@ SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) { vcg::tri::MeshAssert::VertexNormalNormalized(mesh); + // bool containsNormals = true; + // for (VertexIterator vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi) + // if (!vi->IsD()) { + // if (fabs(vi->cN().Norm() - 1.0) > 0.000001) { + // containsNormals = false; + // break; + // } + // } + + // if (!containsNormals) { + // for (VertexIterator vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi) + // if (!vi->IsD()) { + // vi->N() = CoordType(1, 0, 0); + // } + // } vcg::tri::Append::MeshCopy(*this, mesh); elements = vcg::tri::Allocator::GetPerEdgeAttribute( @@ -51,8 +66,7 @@ SimulationMesh::SimulationMesh(SimulationMesh &mesh) { } void SimulationMesh::computeElementalProperties() { - const std::vector elementalDimensions = - getBeamDimensions(); + const std::vector elementalDimensions = getBeamDimensions(); const std::vector elementalMaterials = getBeamMaterial(); assert(EN() == elementalDimensions.size() && elementalDimensions.size() == elementalMaterials.size()); @@ -118,7 +132,7 @@ void SimulationMesh::initializeElements() { Element &element = elements[e]; element.ei = getIndex(e); // Initialize dimensions - element.properties.dimensions = CrossSectionType(1, 1); + element.properties.dimensions = CrossSectionType(); // Initialize material element.properties.material = ElementMaterial(); // Initialize lengths @@ -162,11 +176,11 @@ void SimulationMesh::updateElementalLengths() { } } -void SimulationMesh::setBeamCrossSection(const double &firstDimension, - const double &secondDimension) { +void SimulationMesh::setBeamCrossSection( + const CrossSectionType &beamDimensions) { for (size_t ei = 0; ei < EN(); ei++) { - elements[ei].properties.setDimensions( - CrossSectionType{firstDimension, secondDimension}); + elements[ei].properties.dimensions = beamDimensions; + elements[ei].properties.computeDimensionsProperties(beamDimensions); elements[ei].updateConstFactors(); } } @@ -178,8 +192,8 @@ void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { } } -std::vector SimulationMesh::getBeamDimensions() { - std::vector beamDimensions(EN()); +std::vector SimulationMesh::getBeamDimensions() { + std::vector beamDimensions(EN()); for (size_t ei = 0; ei < EN(); ei++) { beamDimensions[ei] = elements[ei].properties.dimensions; } @@ -198,8 +212,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) { this->Clear(); // assert(plyFileHasAllRequiredFields(plyFilename)); // Load the ply file - VCGEdgeMesh::PerEdgeAttributeHandle - handleBeamDimensions; + VCGEdgeMesh::PerEdgeAttributeHandle handleBeamDimensions; VCGEdgeMesh::PerEdgeAttributeHandle handleBeamMaterial; handleBeamDimensions = vcg::tri::Allocator::AddPerEdgeAttribute( @@ -209,7 +222,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) { *this, plyPropertyBeamMaterialID); nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; customAttrib.GetMeshAttrib(plyFilename); - customAttrib.AddEdgeAttribDescriptor( + customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); @@ -228,7 +241,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) { bool SimulationMesh::savePly(const std::string &plyFilename) { nanoply::NanoPlyWrapper::CustomAttributeDescriptor customAttrib; customAttrib.GetMeshAttrib(plyFilename); - customAttrib.AddEdgeAttribDescriptor( + customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, getBeamDimensions().data()); customAttrib.AddEdgeAttribDescriptor( @@ -318,7 +331,7 @@ void Element::Properties::computeDimensionsProperties( } void Element::Properties::computeDimensionsProperties( - const CylindricalElementDimensions &dimensions) { + const CylindricalBeamDimensions &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; diff --git a/elementalmesh.hpp b/elementalmesh.hpp index 7357309..463aa98 100644 --- a/elementalmesh.hpp +++ b/elementalmesh.hpp @@ -8,6 +8,7 @@ struct Element; struct Node; using CrossSectionType = RectangularBeamDimensions; +// using CrossSectionType = CylindricalBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: @@ -26,9 +27,6 @@ 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); @@ -42,6 +40,8 @@ public: double totalResidualForcesNorm{0}; double totalPotentialEnergykN{0}; bool savePly(const std::string &plyFilename); + void setBeamCrossSection(const CrossSectionType &beamDimensions); + void setBeamMaterial(const double &pr, const double &ym); }; struct Element { @@ -58,7 +58,7 @@ struct Element { void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); void - computeDimensionsProperties(const CylindricalElementDimensions &dimensions); + computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); void setDimensions(const CrossSectionType &dimensions); void setMaterial(const ElementMaterial &material); Properties(const CrossSectionType &dimensions, diff --git a/flatpattern.hpp b/flatpattern.hpp index 4878eef..d6b4008 100644 --- a/flatpattern.hpp +++ b/flatpattern.hpp @@ -27,7 +27,7 @@ private: void removeDuplicateVertices(); void scale(); - const double desiredBaseTriangleCentralEdgeSize{0.025}; + const double desiredBaseTriangleCentralEdgeSize{0.03}; }; #endif // FLATPATTERN_HPP diff --git a/simulationresult.hpp b/simulationresult.hpp index 9d2c578..da4492a 100644 --- a/simulationresult.hpp +++ b/simulationresult.hpp @@ -36,6 +36,7 @@ struct SimulationJob { std::unordered_map> fixedVertices; std::unordered_map nodalExternalForces; std::unordered_map nodalForcedDisplacements; + std::unordered_map nodalForcedNormals; void registerForDrawing() const { initPolyscope(); diff --git a/utilities.hpp b/utilities.hpp index 73d5107..5571974 100644 --- a/utilities.hpp +++ b/utilities.hpp @@ -237,11 +237,13 @@ inline void registerWorldAxes() { if (!polyscope::state::initialized) { initPolyscope(); } + Eigen::MatrixX3d axesPositions(4, 3); axesPositions.row(0) = Eigen::Vector3d(0, 0, 0); - axesPositions.row(1) = Eigen::Vector3d(1, 0, 0); - axesPositions.row(2) = Eigen::Vector3d(0, 1, 0); - axesPositions.row(3) = Eigen::Vector3d(0, 0, 1); + axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0); + axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0); + axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale); + Eigen::MatrixX2i axesEdges(3, 2); axesEdges.row(0) = Eigen::Vector2i(0, 1); axesEdges.row(1) = Eigen::Vector2i(0, 2);