manolas
/
MySources
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Removed the Properties struct from within the Element struct

This commit is contained in:
Iason 2021-01-15 21:12:42 +02:00
parent a0833efc3a
commit bf77c9e64d
3 changed files with 58 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
beamformfinder.cpp
View File

@ -752,21 +752,19 @@ double FormFinder::computeTotalPotentialEnergy() const {
const Element &element = mesh->elements[e]; const Element &element = mesh->elements[e];
const EdgeIndex ei = mesh->getIndex(e); const EdgeIndex ei = mesh->getIndex(e);
const double e_k = element.length - element.initialLength; const double e_k = element.length - element.initialLength;
const double axialPE = pow(e_k, 2) * const double axialPE = pow(e_k, 2) * element.material.youngsModulus *
element.properties.material.youngsModulus * element.A / (2 * element.initialLength);
element.properties.A / (2 * element.initialLength);
const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 - const double theta1Diff = element.rotationalDisplacements_jplus1.theta1 -
element.rotationalDisplacements_j.theta1; element.rotationalDisplacements_j.theta1;
const double torsionalPE = element.properties.G * element.properties.J * const double torsionalPE = element.G * element.J * pow(theta1Diff, 2) /
pow(theta1Diff, 2) / (2 * element.initialLength); (2 * element.initialLength);
const double &theta2_j = element.rotationalDisplacements_j.theta2; const double &theta2_j = element.rotationalDisplacements_j.theta2;
const double &theta2_jplus1 = element.rotationalDisplacements_jplus1.theta2; const double &theta2_jplus1 = element.rotationalDisplacements_jplus1.theta2;
const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) + const double firstBendingPE_inBracketsTerm = 4 * pow(theta2_j, 2) +
4 * theta2_j * theta2_jplus1 + 4 * theta2_j * theta2_jplus1 +
4 * pow(theta2_jplus1, 2); 4 * pow(theta2_jplus1, 2);
const double firstBendingPE = firstBendingPE_inBracketsTerm * const double firstBendingPE = firstBendingPE_inBracketsTerm *
element.properties.material.youngsModulus * element.material.youngsModulus * element.I2 /
element.properties.I2 /
(2 * element.initialLength); (2 * element.initialLength);
const double &theta3_j = element.rotationalDisplacements_j.theta3; const double &theta3_j = element.rotationalDisplacements_j.theta3;
const double &theta3_jplus1 = element.rotationalDisplacements_jplus1.theta3; const double &theta3_jplus1 = element.rotationalDisplacements_jplus1.theta3;
@ -774,8 +772,7 @@ double FormFinder::computeTotalPotentialEnergy() const {
4 * theta3_j * theta3_jplus1 + 4 * theta3_j * theta3_jplus1 +
4 * pow(theta3_jplus1, 2); 4 * pow(theta3_jplus1, 2);
const double secondBendingPE = secondBendingPE_inBracketsTerm * 2 * const double secondBendingPE = secondBendingPE_inBracketsTerm * 2 *
element.properties.material.youngsModulus * element.material.youngsModulus * element.I3 /
element.properties.I3 /
element.initialLength; element.initialLength;
totalInternalPotentialEnergy += totalInternalPotentialEnergy +=
@ -1145,18 +1142,18 @@ void FormFinder::updateNodalMasses() {
for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { for (const EdgePointer &ep : mesh->nodes[v].incidentElements) {
const size_t ei = mesh->getIndex(ep); const size_t ei = mesh->getIndex(ep);
const Element &elem = mesh->elements[ep]; const Element &elem = mesh->elements[ep];
const double SkTranslational = elem.properties.material.youngsModulus * const double SkTranslational =
elem.properties.A / elem.length; elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational; translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3); const double lengthToThe3 = std::pow(elem.length, 3);
const long double SkRotational_I2 = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const long double SkRotational_I3 = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const long double SkRotational_J = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_I3 += SkRotational_I3;
@ -1453,18 +1450,18 @@ void FormFinder::updatePositionsOnTheFly(
double rotationalSumSk_J = 0; double rotationalSumSk_J = 0;
for (const EdgePointer &ep : mesh->nodes[v].incidentElements) { for (const EdgePointer &ep : mesh->nodes[v].incidentElements) {
const Element &elem = mesh->elements[ep]; const Element &elem = mesh->elements[ep];
const double SkTranslational = elem.properties.material.youngsModulus * const double SkTranslational =
elem.properties.A / elem.length; elem.material.youngsModulus * elem.A / elem.length;
translationalSumSk += SkTranslational; translationalSumSk += SkTranslational;
const double lengthToThe3 = std::pow(elem.length, 3); const double lengthToThe3 = std::pow(elem.length, 3);
const double SkRotational_I2 = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_I3 = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
const double SkRotational_J = 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 lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia
rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I2 += SkRotational_I2;
rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_I3 += SkRotational_I3;
@ -1695,10 +1692,10 @@ void FormFinder::draw(const std::string &screenshotsFolder = {}) {
std::vector<double> I3(mesh->EN()); std::vector<double> I3(mesh->EN());
for (const EdgeType &e : mesh->edge) { for (const EdgeType &e : mesh->edge) {
const size_t ei = mesh->getIndex(e); const size_t ei = mesh->getIndex(e);
A[ei] = mesh->elements[e].properties.A; A[ei] = mesh->elements[e].A;
J[ei] = mesh->elements[e].properties.J; J[ei] = mesh->elements[e].J;
I2[ei] = mesh->elements[e].properties.I2; I2[ei] = mesh->elements[e].I2;
I3[ei] = mesh->elements[e].properties.I3; I3[ei] = mesh->elements[e].I3;
} }
polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A); polyscope::getCurveNetwork(meshPolyscopeLabel)->addEdgeScalarQuantity("A", A);

72
elementalmesh.cpp
View File

@ -80,8 +80,8 @@ void SimulationMesh::computeElementalProperties() {
for (const EdgeType &e : edge) { for (const EdgeType &e : edge) {
const EdgeIndex ei = getIndex(e); const EdgeIndex ei = getIndex(e);
elements[e].properties = elements[e].setDimensions(elementalDimensions[ei]);
Element::Properties{elementalDimensions[ei], elementalMaterials[ei]}; elements[e].setMaterial(elementalMaterials[ei]);
} }
} }
@ -184,9 +184,9 @@ void SimulationMesh::initializeElements() {
Element &element = elements[e]; Element &element = elements[e];
element.ei = getIndex(e); element.ei = getIndex(e);
// Initialize dimensions // Initialize dimensions
element.properties.dimensions = CrossSectionType(); element.dimensions = CrossSectionType();
// Initialize material // Initialize material
element.properties.material = ElementMaterial(); element.material = ElementMaterial();
// Initialize lengths // Initialize lengths
const VCGEdgeMesh::CoordType p0 = e.cP(0); const VCGEdgeMesh::CoordType p0 = e.cP(0);
const VCGEdgeMesh::CoordType p1 = e.cP(1); const VCGEdgeMesh::CoordType p1 = e.cP(1);
@ -231,15 +231,15 @@ void SimulationMesh::updateElementalLengths() {
void SimulationMesh::setBeamCrossSection( void SimulationMesh::setBeamCrossSection(
const CrossSectionType &beamDimensions) { const CrossSectionType &beamDimensions) {
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties.dimensions = beamDimensions; elements[ei].dimensions = beamDimensions;
elements[ei].properties.computeDimensionsProperties(beamDimensions); elements[ei].computeDimensionsProperties(beamDimensions);
elements[ei].updateRigidity(); elements[ei].updateRigidity();
} }
} }
void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) { void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties.setMaterial(ElementMaterial{pr, ym}); elements[ei].setMaterial(ElementMaterial{pr, ym});
elements[ei].updateRigidity(); elements[ei].updateRigidity();
} }
} }
@ -247,7 +247,7 @@ void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {
std::vector<CrossSectionType> SimulationMesh::getBeamDimensions() { std::vector<CrossSectionType> SimulationMesh::getBeamDimensions() {
std::vector<CrossSectionType> beamDimensions(EN()); std::vector<CrossSectionType> beamDimensions(EN());
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
beamDimensions[ei] = elements[ei].properties.dimensions; beamDimensions[ei] = elements[ei].dimensions;
} }
return beamDimensions; return beamDimensions;
} }
@ -255,7 +255,7 @@ std::vector<CrossSectionType> SimulationMesh::getBeamDimensions() {
std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() { std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
std::vector<ElementMaterial> beamMaterial(EN()); std::vector<ElementMaterial> beamMaterial(EN());
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
beamMaterial[ei] = elements[ei].properties.material; beamMaterial[ei] = elements[ei].material;
} }
return beamMaterial; return beamMaterial;
} }
@ -315,14 +315,14 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
// if (!handleBeamProperties._handle->data.empty()) { // if (!handleBeamProperties._handle->data.empty()) {
// for (size_t ei = 0; ei < EN(); ei++) { // for (size_t ei = 0; ei < EN(); ei++) {
// elements[ei].properties = // elements[ei] =
// Element::Properties(handleBeamProperties[ei]); // Element::Properties(handleBeamProperties[ei]);
// elements[ei].updateRigidity(); // elements[ei].updateRigidity();
// } // }
// } // }
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties = elements[ei].setDimensions(handleBeamDimensions[ei]);
Element::Properties(handleBeamDimensions[ei], handleBeamMaterial[ei]); elements[ei].setMaterial(handleBeamMaterial[ei]);
elements[ei].updateRigidity(); elements[ei].updateRigidity();
} }
@ -343,7 +343,7 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
material.data()); material.data());
// std::vector<std::array<double, 6>> beamProperties(EN()); // std::vector<std::array<double, 6>> beamProperties(EN());
// for (size_t ei = 0; ei < EN(); ei++) { // for (size_t ei = 0; ei < EN(); ei++) {
// auto props = elements[ei].properties.toArray(); // auto props = elements[ei].toArray();
// for (auto p : props) { // for (auto p : props) {
// std::cout << p << " "; // std::cout << p << " ";
// } // }
@ -422,12 +422,11 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
return angle; return angle;
} }
void Element::Properties::computeMaterialProperties( void Element::computeMaterialProperties(const ElementMaterial &material) {
const ElementMaterial &material) {
this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio)); this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
} }
void Element::Properties::computeDimensionsProperties( void Element::computeDimensionsProperties(
const RectangularBeamDimensions &dimensions) { const RectangularBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions)); assert(typeid(CrossSectionType) == typeid(RectangularBeamDimensions));
A = (dimensions.b * dimensions.h); A = (dimensions.b * dimensions.h);
@ -436,7 +435,7 @@ void Element::Properties::computeDimensionsProperties(
J = I2 + I3; J = I2 + I3;
} }
void Element::Properties::computeDimensionsProperties( void Element::computeDimensionsProperties(
const CylindricalBeamDimensions &dimensions) { const CylindricalBeamDimensions &dimensions) {
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions)); assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2)); 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; J = I2 + I3;
} }
void Element::Properties::setDimensions(const CrossSectionType &dimensions) { void Element::setDimensions(const CrossSectionType &dimensions) {
this->dimensions = dimensions; this->dimensions = dimensions;
computeDimensionsProperties(dimensions); computeDimensionsProperties(dimensions);
updateRigidity();
} }
void Element::Properties::setMaterial(const ElementMaterial &material) { void Element::setMaterial(const ElementMaterial &material) {
this->material = material; this->material = material;
computeMaterialProperties(material); computeMaterialProperties(material);
} updateRigidity();
Element::Properties::Properties(const CrossSectionType &dimensions,
const ElementMaterial &material) {
setDimensions(dimensions);
setMaterial(material);
}
Element::Properties::Properties(const std::array<double, 6> &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<double, 6> Element::Properties::toArray() const {
return std::array<double, 6>({material.youngsModulus, G, A,
static_cast<double>(I2),
static_cast<double>(I3), J});
} }
void Element::updateRigidity() { void Element::updateRigidity() {
rigidity.axial = rigidity.axial = material.youngsModulus * A / initialLength;
properties.material.youngsModulus * properties.A / initialLength; rigidity.torsional = G * J / initialLength;
rigidity.torsional = properties.G * properties.J / initialLength; rigidity.firstBending = 2 * material.youngsModulus * I2 / initialLength;
rigidity.firstBending = rigidity.secondBending = 2 * material.youngsModulus * I3 / initialLength;
2 * properties.material.youngsModulus * properties.I2 / initialLength;
rigidity.secondBending =
2 * properties.material.youngsModulus * properties.I3 / initialLength;
} }

21
elementalmesh.hpp
View File

@ -7,8 +7,8 @@
struct Element; struct Element;
struct Node; struct Node;
// using CrossSectionType = RectangularBeamDimensions; using CrossSectionType = RectangularBeamDimensions;
using CrossSectionType = CylindricalBeamDimensions; // using CrossSectionType = CylindricalBeamDimensions;
class SimulationMesh : public VCGEdgeMesh { class SimulationMesh : public VCGEdgeMesh {
private: private:
@ -53,7 +53,6 @@ public:
struct Element { struct Element {
struct Properties {
CrossSectionType dimensions; CrossSectionType dimensions;
ElementMaterial material; ElementMaterial material;
double G{0}; double G{0};
@ -61,21 +60,12 @@ struct Element {
double I2{0}; // second moment of inertia double I2{0}; // second moment of inertia
double I3{0}; // third moment of inertia double I3{0}; // third moment of inertia
double J{0}; // torsional constant (polar moment of inertia) double J{0}; // torsional constant (polar moment of inertia)
void computeMaterialProperties(const ElementMaterial &material); void computeMaterialProperties(const ElementMaterial &material);
void void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
computeDimensionsProperties(const RectangularBeamDimensions &dimensions); void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void
computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
void setDimensions(const CrossSectionType &dimensions); void setDimensions(const CrossSectionType &dimensions);
void setMaterial(const ElementMaterial &material); void setMaterial(const ElementMaterial &material);
Properties(const CrossSectionType &dimensions,
const ElementMaterial &material);
Properties(const std::array<double, 6> &arr);
Properties() {}
public:
std::array<double, 6> toArray() const;
};
struct LocalFrame { struct LocalFrame {
VectorType t1; VectorType t1;
@ -85,7 +75,6 @@ struct Element {
EdgeIndex ei; EdgeIndex ei;
double length{0}; double length{0};
Properties properties;
double initialLength; double initialLength;
LocalFrame frame; LocalFrame frame;