diff --git a/beam.hpp b/beam.hpp index c50f12c..2410c81 100644 --- a/beam.hpp +++ b/beam.hpp @@ -31,14 +31,13 @@ struct CylindricalBeamDimensions { }; struct ElementMaterial { - float poissonsRatio; - float youngsModulusGPascal; - ElementMaterial(const float &poissonsRatio, const float &youngsModulusGPascal) - : poissonsRatio(poissonsRatio), - youngsModulusGPascal(youngsModulusGPascal) { + float G; // poisson's ratio + float E; // ym in pascal + ElementMaterial(const float &poissonsRatio, const float &youngsModulus) + : G(poissonsRatio), E(youngsModulus) { assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); } - ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(200) {} + ElementMaterial() : G(0.3), E(200 * 1e10) {} }; #endif // BEAM_HPP diff --git a/beamformfinder.cpp b/beamformfinder.cpp index a485a96..1abd61b 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -8,8 +8,6 @@ #include #include -const bool debug = true; - void FormFinder::runUnitTests() { const std::filesystem::path groundOfTruthFolder{ "/home/iason/Coding/Libraries/MySources/formFinder_unitTestFiles"}; @@ -36,7 +34,7 @@ void FormFinder::runUnitTests() { // SimulationJob::constructFixedVerticesSpanGrid(spanGridSize, // mesh.VN()), fixedVertices, nodalForces, nodalForcedDisplacements}; - beamSimulationJob.pMesh->setBeamMaterial(0.3, 200); + beamSimulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e10); assert(CrossSectionType::name == CylindricalBeamDimensions::name); beamSimulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); @@ -45,6 +43,8 @@ void FormFinder::runUnitTests() { settings.xi = 0.9969; settings.maxDRMIterations = 0; settings.totalResidualForcesNormThreshold = 1; + settings.shouldDraw = true; + settings.beVerbose = true; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation( std::make_shared(beamSimulationJob), settings); simpleBeam_simulationResults.save(); @@ -59,7 +59,7 @@ void FormFinder::runUnitTests() { if (!simpleBeam_simulationResults.isEqual( simpleBeam_groundOfTruthDisplacements)) { std::cerr << "Error:Beam simulation produces wrong results." << std::endl; - return; + // return; } // Second example of the paper @@ -91,7 +91,7 @@ void FormFinder::runUnitTests() { fixedVertices, {}, nodalForcedDisplacements}; - shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200); + shortSpanGridshellSimulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e10); assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); shortSpanGridshellSimulationJob.pMesh->setBeamCrossSection( CrossSectionType{0.03, 0.026}); @@ -169,7 +169,7 @@ void FormFinder::runUnitTests() { fixedVertices, {}, nodalForcedDisplacements}; - longSpanGridshell_simulationJob.pMesh->setBeamMaterial(0.3, 200); + longSpanGridshell_simulationJob.pMesh->setBeamMaterial(0.3, 200 * 1e10); if (typeid(CrossSectionType) != typeid(CylindricalBeamDimensions)) { std::cerr << "A cylindrical cross section is expected for running the " "paper examples." @@ -752,27 +752,28 @@ double FormFinder::computeTotalPotentialEnergy() const { const Element &element = mesh->elements[e]; const EdgeIndex ei = mesh->getIndex(e); const double e_k = element.length - element.initialLength; - const double axialPE = pow(e_k, 2) * element.properties.E * + const double axialPE = pow(e_k, 2) * element.properties.material.E * element.properties.A / (2 * element.initialLength); const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 - element.rotationalDisplacements_j.theta1; - const double torsionalPE = element.properties.G * element.properties.J * - pow(theta1Diff, 2) / (2 * element.initialLength); + const double torsionalPE = element.properties.material.G * + element.properties.J * pow(theta1Diff, 2) / + (2 * element.initialLength); const double &theta2_j = element.rotationalDisplacements_j.theta2; const double &theta2_jplus1 = element.rotationalDisplacements_jplus1.theta2; const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) + 4 * theta2_j * theta2_jplus1 + 4 * pow(theta2_jplus1, 2); - const double firstBendingPE = firstBendingPE_inBracketsTerm * - element.properties.E * element.properties.I2 / - (2 * element.initialLength); + const double firstBendingPE = + firstBendingPE_inBracketsTerm * element.properties.material.E * + element.properties.I2 / (2 * element.initialLength); const double &theta3_j = element.rotationalDisplacements_j.theta3; const double &theta3_jplus1 = element.rotationalDisplacements_jplus1.theta3; const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) + 4 * theta3_j * theta3_jplus1 + 4 * pow(theta3_jplus1, 2); const double secondBendingPE = - secondBendingPE_inBracketsTerm * 2 * element.properties.E * + secondBendingPE_inBracketsTerm * 2 * element.properties.material.E * element.properties.I3 / element.initialLength; totalInternalPotentialEnergy += @@ -1143,17 +1144,17 @@ void FormFinder::updateNodalMasses() { const size_t ei = mesh->getIndex(ep); const Element &elem = mesh->elements[ep]; const double SkTranslational = - elem.properties.E * elem.properties.A / elem.length; + elem.properties.material.E * elem.properties.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const double SkRotational_I2 = - elem.properties.E * elem.properties.I2 / + const long double SkRotational_I2 = + elem.properties.material.E * elem.properties.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_I3 = - elem.properties.E * elem.properties.I3 / + const long double SkRotational_I3 = + elem.properties.material.E * elem.properties.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const double SkRotational_J = - elem.properties.E * elem.properties.J / + const long double SkRotational_J = + elem.properties.material.E * elem.properties.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1344,7 +1345,7 @@ void FormFinder::applyDisplacements( updateNodeNormal(v, fixedVertices); } updateElementalFrames(); - if (mSettings.shouldDraw || true) { + if (mSettings.shouldDraw) { mesh->updateEigenEdgeAndVertices(); } } @@ -1370,9 +1371,15 @@ 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]}); + // node.displacements = Vector6d( + // {vertexDisplacement(0), vertexDisplacement(1), + // vertexDisplacement(2), + // node.displacements[3], node.displacements[4], + // node.displacements[5]}); + } + + if (mSettings.shouldDraw) { + mesh->updateEigenEdgeAndVertices(); } } @@ -1447,17 +1454,17 @@ void FormFinder::updatePositionsOnTheFly( for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { const Element &elem = mesh->elements[ep]; const double SkTranslational = - elem.properties.E * elem.properties.A / elem.length; + elem.properties.material.E * elem.properties.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); const double SkRotational_I2 = - elem.properties.E * elem.properties.I2 / + elem.properties.material.E * elem.properties.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_I3 = - elem.properties.E * elem.properties.I3 / + elem.properties.material.E * elem.properties.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_J = - elem.properties.E * elem.properties.J / + elem.properties.material.E * elem.properties.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1738,7 +1745,7 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { } } -void FormFinder::printCurrentState() { +void FormFinder::printCurrentState() const { std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl; std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm @@ -1750,6 +1757,7 @@ void FormFinder::printCurrentState() { void FormFinder::printDebugInfo() const { std::cout << mesh->elements[0].rigidity.toString() << std::endl; std::cout << "Number of dampings:" << numOfDampings << endl; + printCurrentState(); } SimulationResults @@ -1760,13 +1768,13 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, reset(); mSettings = settings; - if (!pJob->nodalExternalForces.empty()) { - double externalForcesNorm = 0; - for (const auto &externalForce : pJob->nodalExternalForces) { - externalForcesNorm += externalForce.second.norm(); - } - mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-3; - } + // if (!pJob->nodalExternalForces.empty()) { + // double externalForcesNorm = 0; + // for (const auto &externalForce : pJob->nodalExternalForces) { + // externalForcesNorm += externalForce.second.norm(); + // } + // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; + // } constrainedVertices = pJob->constrainedVertices; if (!pJob->nodalForcedDisplacements.empty()) { @@ -1796,7 +1804,7 @@ FormFinder::executeSimulation(const std::shared_ptr &pJob, } polyscope::registerCurveNetwork( meshPolyscopeLabel, mesh->getEigenVertices(), mesh->getEigenEdges()); - registerWorldAxes(); + // registerWorldAxes(); } for (auto fixedVertex : pJob->constrainedVertices) { assert(fixedVertex.first < mesh->VN()); @@ -1912,7 +1920,7 @@ currentSimulationStep > maxDRMIterations*/ break; } // Residual forces norm convergence - if (mesh->previousTotalKineticEnergy > mesh->currentTotalKineticEnergy + if (mesh->previousTotalKineticEnergy >= mesh->currentTotalKineticEnergy /*|| mesh->previousTotalPotentialEnergykN > mesh->currentTotalPotentialEnergykN*/ @@ -1944,7 +1952,7 @@ mesh->currentTotalPotentialEnergykN*/ ++numOfDampings; } - if (debug) { + if (mSettings.debugMessages) { printDebugInfo(); } } diff --git a/beamformfinder.hpp b/beamformfinder.hpp index 70ae184..8ccee34 100644 --- a/beamformfinder.hpp +++ b/beamformfinder.hpp @@ -25,6 +25,7 @@ struct DifferentiateWithRespectTo { class FormFinder { public: struct Settings { + bool debugMessages{false}; bool shouldDraw{false}; bool beVerbose{false}; bool shouldCreatePlots{false}; @@ -194,7 +195,7 @@ private: void applyForcedNormals( const std::unordered_map nodalForcedRotations); - void printCurrentState(); + void printCurrentState() const; void printDebugInfo() const; diff --git a/elementalmesh.cpp b/elementalmesh.cpp index beacd15..3d0b839 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -335,18 +335,18 @@ bool SimulationMesh::savePly(const std::string &plyFilename) { customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32, material.data()); - std::vector> beamProperties(EN()); - for (size_t ei = 0; ei < EN(); ei++) { - auto props = elements[ei].properties.toArray(); - for (auto p : props) { - std::cout << p << " "; - } - std::cout << std::endl; - beamProperties[ei] = props; - } - customAttrib.AddEdgeAttribDescriptor, double, 6>( - plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, - beamProperties.data()); + // std::vector> beamProperties(EN()); + // for (size_t ei = 0; ei < EN(); ei++) { + // auto props = elements[ei].properties.toArray(); + // for (auto p : props) { + // std::cout << p << " "; + // } + // std::cout << std::endl; + // beamProperties[ei] = props; + // } + // customAttrib.AddEdgeAttribDescriptor, double, 6>( + // plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, + // beamProperties.data()); // Load the ply file unsigned int mask = 0; mask |= nanoply::NanoPlyWrapper::IO_VERTCOORD; @@ -418,12 +418,12 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, void Element::Properties::computeMaterialProperties( const ElementMaterial &material) { - E = material.youngsModulusGPascal * std::pow(10, 9); - G = E / (2 * (1 + material.poissonsRatio)); + this->material = material; } void Element::Properties::computeDimensionsProperties( const RectangularBeamDimensions &dimensions) { + assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); A = (dimensions.b * dimensions.h); I2 = dimensions.b * std::pow(dimensions.h, 3) / 12; I3 = dimensions.h * std::pow(dimensions.b, 3) / 12; @@ -432,6 +432,7 @@ void Element::Properties::computeDimensionsProperties( void Element::Properties::computeDimensionsProperties( const CylindricalBeamDimensions &dimensions) { + assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); 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; @@ -457,8 +458,8 @@ Element::Properties::Properties(const CrossSectionType &dimensions, Element::Properties::Properties(const std::array &arr) { assert(arr.size() == 6); - E = arr[0]; - G = arr[1]; + material.E = arr[0]; + material.G = arr[1]; A = arr[2]; I2 = arr[3]; I3 = arr[4]; @@ -466,12 +467,16 @@ Element::Properties::Properties(const std::array &arr) { } std::array Element::Properties::toArray() const { - return std::array({E, G, A, I2, I3, J}); + return std::array({material.E, material.G, A, + static_cast(I2), + static_cast(I3), J}); } void Element::updateRigidity() { - rigidity.axial = properties.E * properties.A / initialLength; - rigidity.torsional = properties.G * properties.J / initialLength; - rigidity.firstBending = 2 * properties.E * properties.I2 / initialLength; - rigidity.secondBending = 2 * properties.E * properties.I3 / initialLength; + rigidity.axial = properties.material.E * properties.A / initialLength; + rigidity.torsional = properties.material.G * properties.J / initialLength; + rigidity.firstBending = + 2 * properties.material.E * properties.I2 / initialLength; + rigidity.secondBending = + 2 * properties.material.E * properties.I3 / initialLength; } diff --git a/elementalmesh.hpp b/elementalmesh.hpp index b125642..4c8a572 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: @@ -56,8 +56,6 @@ 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