diff --git a/beamformfinder.cpp b/beamformfinder.cpp index d9e1060..f0e22ea 100644 --- a/beamformfinder.cpp +++ b/beamformfinder.cpp @@ -752,21 +752,19 @@ 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.youngsModulus * - element.properties.A / (2 * element.initialLength); + const double axialPE = pow(e_k, 2) * element.material.youngsModulus * + element.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.G * element.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.material.youngsModulus * - element.properties.I2 / + element.material.youngsModulus * element.I2 / (2 * element.initialLength); const double &theta3_j = element.rotationalDisplacements_j.theta3; const double &theta3_jplus1 = element.rotationalDisplacements_jplus1.theta3; @@ -774,8 +772,7 @@ double FormFinder::computeTotalPotentialEnergy() const { 4 * theta3_j * theta3_jplus1 + 4 * pow(theta3_jplus1, 2); const double secondBendingPE = secondBendingPE_inBracketsTerm * 2 * - element.properties.material.youngsModulus * - element.properties.I3 / + element.material.youngsModulus * element.I3 / element.initialLength; totalInternalPotentialEnergy += @@ -1145,18 +1142,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.youngsModulus * - elem.properties.A / elem.length; + const double SkTranslational = + elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); const long double SkRotational_I2 = - elem.properties.material.youngsModulus * elem.properties.I2 / + elem.material.youngsModulus * elem.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const long double SkRotational_I3 = - elem.properties.material.youngsModulus * elem.properties.I3 / + elem.material.youngsModulus * elem.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const long double SkRotational_J = - elem.properties.material.youngsModulus * elem.properties.J / + elem.material.youngsModulus * elem.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1453,18 +1450,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.youngsModulus * - elem.properties.A / elem.length; + const double SkTranslational = + elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); const double SkRotational_I2 = - elem.properties.material.youngsModulus * elem.properties.I2 / + elem.material.youngsModulus * elem.I2 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_I3 = - elem.properties.material.youngsModulus * elem.properties.I3 / + elem.material.youngsModulus * elem.I3 / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia const double SkRotational_J = - elem.properties.material.youngsModulus * elem.properties.J / + elem.material.youngsModulus * elem.J / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; @@ -1695,10 +1692,10 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) { std::vector I3(mesh->EN()); for (const EdgeType &e : mesh->edge) { const size_t ei = mesh->getIndex(e); - A[ei] = mesh->elements[e].properties.A; - J[ei] = mesh->elements[e].properties.J; - I2[ei] = mesh->elements[e].properties.I2; - I3[ei] = mesh->elements[e].properties.I3; + A[ei] = mesh->elements[e].A; + J[ei] = mesh->elements[e].J; + I2[ei] = mesh->elements[e].I2; + I3[ei] = mesh->elements[e].I3; } polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); diff --git a/elementalmesh.cpp b/elementalmesh.cpp index 79a1458..e8da8e1 100644 --- a/elementalmesh.cpp +++ b/elementalmesh.cpp @@ -80,8 +80,8 @@ void SimulationMesh::computeElementalProperties() { for (const EdgeType &e : edge) { const EdgeIndex ei = getIndex(e); - elements[e].properties = - Element::Properties{elementalDimensions[ei], elementalMaterials[ei]}; + elements[e].setDimensions(elementalDimensions[ei]); + elements[e].setMaterial(elementalMaterials[ei]); } } @@ -184,9 +184,9 @@ void SimulationMesh::initializeElements() { Element &element = elements[e]; element.ei = getIndex(e); // Initialize dimensions - element.properties.dimensions = CrossSectionType(); + element.dimensions = CrossSectionType(); // Initialize material - element.properties.material = ElementMaterial(); + element.material = ElementMaterial(); // Initialize lengths const VCGEdgeMesh::CoordType p0 = e.cP(0); const VCGEdgeMesh::CoordType p1 = e.cP(1); @@ -231,15 +231,15 @@ void SimulationMesh::updateElementalLengths() { void SimulationMesh::setBeamCrossSection( const CrossSectionType &beamDimensions) { for (size_t ei = 0; ei < EN(); ei++) { - elements[ei].properties.dimensions = beamDimensions; - elements[ei].properties.computeDimensionsProperties(beamDimensions); + elements[ei].dimensions = beamDimensions; + elements[ei].computeDimensionsProperties(beamDimensions); elements[ei].updateRigidity(); } } 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].setMaterial(ElementMaterial{pr, ym}); elements[ei].updateRigidity(); } } @@ -247,7 +247,7 @@ void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { std::vector SimulationMesh::getBeamDimensions() { std::vector beamDimensions(EN()); for (size_t ei = 0; ei < EN(); ei++) { - beamDimensions[ei] = elements[ei].properties.dimensions; + beamDimensions[ei] = elements[ei].dimensions; } return beamDimensions; } @@ -255,7 +255,7 @@ std::vector SimulationMesh::getBeamDimensions() { std::vector SimulationMesh::getBeamMaterial() { std::vector beamMaterial(EN()); for (size_t ei = 0; ei < EN(); ei++) { - beamMaterial[ei] = elements[ei].properties.material; + beamMaterial[ei] = elements[ei].material; } return beamMaterial; } @@ -315,14 +315,14 @@ bool SimulationMesh::loadPly(const string &plyFilename) { // if (!handleBeamProperties._handle->data.empty()) { // for (size_t ei = 0; ei < EN(); ei++) { - // elements[ei].properties = + // elements[ei] = // 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].setDimensions(handleBeamDimensions[ei]); + elements[ei].setMaterial(handleBeamMaterial[ei]); elements[ei].updateRigidity(); } @@ -343,7 +343,7 @@ bool SimulationMesh::savePly(const std::string &plyFilename) { material.data()); // std::vector> beamProperties(EN()); // for (size_t ei = 0; ei < EN(); ei++) { - // auto props = elements[ei].properties.toArray(); + // auto props = elements[ei].toArray(); // for (auto p : props) { // std::cout << p << " "; // } @@ -422,12 +422,11 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1, return angle; } -void Element::Properties::computeMaterialProperties( - const ElementMaterial &material) { +void Element::computeMaterialProperties(const ElementMaterial &material) { this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); } -void Element::Properties::computeDimensionsProperties( +void Element::computeDimensionsProperties( const RectangularBeamDimensions &dimensions) { assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); A = (dimensions.b * dimensions.h); @@ -436,7 +435,7 @@ void Element::Properties::computeDimensionsProperties( J = I2 + I3; } -void Element::Properties::computeDimensionsProperties( +void Element::computeDimensionsProperties( const CylindricalBeamDimensions &dimensions) { assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); @@ -445,44 +444,21 @@ void Element::Properties::computeDimensionsProperties( J = I2 + I3; } -void Element::Properties::setDimensions(const CrossSectionType &dimensions) { +void Element::setDimensions(const CrossSectionType &dimensions) { this->dimensions = dimensions; computeDimensionsProperties(dimensions); + updateRigidity(); } -void Element::Properties::setMaterial(const ElementMaterial &material) { +void Element::setMaterial(const ElementMaterial &material) { this->material = material; computeMaterialProperties(material); -} - -Element::Properties::Properties(const CrossSectionType &dimensions, - const ElementMaterial &material) { - setDimensions(dimensions); - setMaterial(material); -} - -Element::Properties::Properties(const std::array &arr) { - assert(arr.size() == 6); - material.youngsModulus = arr[0]; - G = arr[1]; - A = arr[2]; - I2 = arr[3]; - I3 = arr[4]; - J = arr[5]; -} - -std::array Element::Properties::toArray() const { - return std::array({material.youngsModulus, G, A, - static_cast(I2), - static_cast(I3), J}); + updateRigidity(); } void Element::updateRigidity() { - rigidity.axial = - properties.material.youngsModulus * properties.A / initialLength; - rigidity.torsional = properties.G * properties.J / initialLength; - rigidity.firstBending = - 2 * properties.material.youngsModulus * properties.I2 / initialLength; - rigidity.secondBending = - 2 * properties.material.youngsModulus * properties.I3 / initialLength; + rigidity.axial = material.youngsModulus * A / initialLength; + rigidity.torsional = G * J / initialLength; + rigidity.firstBending = 2 * material.youngsModulus * I2 / initialLength; + rigidity.secondBending = 2 * material.youngsModulus * I3 / initialLength; } diff --git a/elementalmesh.hpp b/elementalmesh.hpp index 1436b2c..25cc456 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: @@ -53,29 +53,19 @@ public: struct Element { - struct Properties { - CrossSectionType dimensions; - ElementMaterial material; - double G{0}; - 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); - void - computeDimensionsProperties(const RectangularBeamDimensions &dimensions); - void - computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); - void setDimensions(const CrossSectionType &dimensions); - void setMaterial(const ElementMaterial &material); - Properties(const CrossSectionType &dimensions, - const ElementMaterial &material); - Properties(const std::array &arr); - Properties() {} + CrossSectionType dimensions; + ElementMaterial material; + double G{0}; + 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) - public: - std::array toArray() const; - }; + void computeMaterialProperties(const ElementMaterial &material); + void computeDimensionsProperties(const RectangularBeamDimensions &dimensions); + void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions); + void setDimensions(const CrossSectionType &dimensions); + void setMaterial(const ElementMaterial &material); struct LocalFrame { VectorType t1; @@ -85,7 +75,6 @@ struct Element { EdgeIndex ei; double length{0}; - Properties properties; double initialLength; LocalFrame frame;