diff --git a/beam.hpp b/beam.hpp index 5380171..f629662 100644 --- a/beam.hpp +++ b/beam.hpp @@ -31,13 +31,14 @@ struct CylindricalBeamDimensions { }; struct ElementMaterial { - float poissonsRatio; // poisson's ratio - float E; // ym in pascal - ElementMaterial(const float &poissonsRatio, const float &youngsModulus) - : poissonsRatio(poissonsRatio), E(youngsModulus) { + float poissonsRatio; // NOTE: if I make this double the unit + // tests produced different results and thus fail + double youngsModulus; + ElementMaterial(const float &poissonsRatio, const double &youngsModulus) + : poissonsRatio(poissonsRatio), youngsModulus(youngsModulus) { assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); } - ElementMaterial() : poissonsRatio(0.3), E(200 * 1e9) {} + ElementMaterial() : poissonsRatio(0.3), youngsModulus(200) {} }; #endif // BEAM_HPP diff --git a/beamformfinder.cpp b/beamformfinder.cpp index 826817c..d9e1060 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -43,7 +43,7 @@ void FormFinder::runUnitTests() { settings.xi = 0.9969; settings.maxDRMIterations = 0; settings.totalResidualForcesNormThreshold = 1; - settings.shouldDraw = true; + // settings.shouldDraw = true; settings.beVerbose = true; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation( std::make_shared(beamSimulationJob), settings); @@ -752,7 +752,8 @@ 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.material.E * + const double axialPE = pow(e_k, 2) * + element.properties.material.youngsModulus * element.properties.A / (2 * element.initialLength); const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 - element.rotationalDisplacements_j.theta1; @@ -763,17 +764,19 @@ double FormFinder::computeTotalPotentialEnergy() const { 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.material.E * - element.properties.I2 / (2 * element.initialLength); + const double firstBendingPE = firstBendingPE_inBracketsTerm * + element.properties.material.youngsModulus * + 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.material.E * - element.properties.I3 / element.initialLength; + const double secondBendingPE = secondBendingPE_inBracketsTerm * 2 * + element.properties.material.youngsModulus * + element.properties.I3 / + element.initialLength; totalInternalPotentialEnergy += axialPE + torsionalPE + firstBendingPE + secondBendingPE; @@ -1142,18 +1145,18 @@ void FormFinder::updateNodalMasses() { for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { const size_t ei = mesh->getIndex(ep); const Element &elem = mesh->elements[ep]; - const double SkTranslational = - elem.properties.material.E * elem.properties.A / elem.length; + const double SkTranslational = elem.properties.material.youngsModulus * + elem.properties.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); const long double SkRotational_I2 = - elem.properties.material.E * elem.properties.I2 / + elem.properties.material.youngsModulus * elem.properties.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const long double SkRotational_I3 = - elem.properties.material.E * elem.properties.I3 / + elem.properties.material.youngsModulus * elem.properties.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const long double SkRotational_J = - elem.properties.material.E * elem.properties.J / + elem.properties.material.youngsModulus * elem.properties.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1450,18 +1453,18 @@ void FormFinder::updatePositionsOnTheFly( double rotationalSumSk_J = 0; for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { const Element &elem = mesh->elements[ep]; - const double SkTranslational = - elem.properties.material.E * elem.properties.A / elem.length; + const double SkTranslational = elem.properties.material.youngsModulus * + elem.properties.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); const double SkRotational_I2 = - elem.properties.material.E * elem.properties.I2 / + elem.properties.material.youngsModulus * elem.properties.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_I3 = - elem.properties.material.E * elem.properties.I3 / + elem.properties.material.youngsModulus * elem.properties.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_J = - elem.properties.material.E * elem.properties.J / + elem.properties.material.youngsModulus * elem.properties.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; diff --git a/elementalmesh.cpp b/elementalmesh.cpp index 6bed494..79a1458 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -282,12 +282,12 @@ bool SimulationMesh::loadPly(const string &plyFilename) { customAttrib.AddEdgeAttribDescriptor( plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32, nullptr); - VCGEdgeMesh::PerEdgeAttributeHandle> - handleBeamProperties = - vcg::tri::Allocator::AddPerEdgeAttribute< - std::array>(*this, plyPropertyBeamProperties); - customAttrib.AddEdgeAttribDescriptor, double, 6>( - plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); + // VCGEdgeMesh::PerEdgeAttributeHandle> + // handleBeamProperties = + // vcg::tri::Allocator::AddPerEdgeAttribute< + // std::array>(*this, plyPropertyBeamProperties); + // customAttrib.AddEdgeAttribDescriptor, double, 6>( + // plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); // VCGEdgeMesh::PerEdgeAttributeHandle // handleBeamRigidityContants; @@ -313,11 +313,17 @@ bool SimulationMesh::loadPly(const string &plyFilename) { initializeElements(); updateEigenEdgeAndVertices(); - if (!handleBeamProperties._handle->data.empty()) { - for (size_t ei = 0; ei < EN(); ei++) { - elements[ei].properties = Element::Properties(handleBeamProperties[ei]); - elements[ei].updateRigidity(); - } + // if (!handleBeamProperties._handle->data.empty()) { + // for (size_t ei = 0; ei < EN(); ei++) { + // elements[ei].properties = + // Element::Properties(handleBeamProperties[ei]); + // elements[ei].updateRigidity(); + // } + // } + for (size_t ei = 0; ei < EN(); ei++) { + elements[ei].properties = + Element::Properties(handleBeamDimensions[ei], handleBeamMaterial[ei]); + elements[ei].updateRigidity(); } return true; @@ -418,8 +424,7 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, void Element::Properties::computeMaterialProperties( const ElementMaterial &material) { - this->material = material; - this->G = material.E / (2 * (1 + material.poissonsRatio)); + this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); } void Element::Properties::computeDimensionsProperties( @@ -451,15 +456,14 @@ void Element::Properties::setMaterial(const ElementMaterial &material) { } Element::Properties::Properties(const CrossSectionType &dimensions, - const ElementMaterial &material) - : dimensions(dimensions), material(material) { - computeDimensionsProperties(dimensions); - computeMaterialProperties(material); + const ElementMaterial &material) { + setDimensions(dimensions); + setMaterial(material); } Element::Properties::Properties(const std::array &arr) { assert(arr.size() == 6); - material.E = arr[0]; + material.youngsModulus = arr[0]; G = arr[1]; A = arr[2]; I2 = arr[3]; @@ -468,15 +472,17 @@ Element::Properties::Properties(const std::array &arr) { } std::array Element::Properties::toArray() const { - return std::array( - {material.E, G, A, static_cast(I2), static_cast(I3), J}); + return std::array({material.youngsModulus, G, A, + static_cast(I2), + static_cast(I3), J}); } void Element::updateRigidity() { - rigidity.axial = properties.material.E * properties.A / initialLength; + rigidity.axial = + properties.material.youngsModulus * properties.A / initialLength; rigidity.torsional = properties.G * properties.J / initialLength; rigidity.firstBending = - 2 * properties.material.E * properties.I2 / initialLength; + 2 * properties.material.youngsModulus * properties.I2 / initialLength; rigidity.secondBending = - 2 * properties.material.E * properties.I3 / initialLength; + 2 * properties.material.youngsModulus * properties.I3 / initialLength; }