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

Removed E(youngs modulus) from element properties and left only the material property youngsModulus

This commit is contained in:
Iason 2021-01-15 19:21:17 +02:00
parent ab1206be9b
commit a0833efc3a
3 changed files with 56 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
beam.hpp
View File

@ -31,13 +31,14 @@ struct CylindricalBeamDimensions {
}; };
struct ElementMaterial { struct ElementMaterial {
float poissonsRatio; // poisson's ratio float poissonsRatio; // NOTE: if I make this double the unit
float E; // ym in pascal // tests produced different results and thus fail
ElementMaterial(const float &poissonsRatio, const float &youngsModulus) double youngsModulus;
: poissonsRatio(poissonsRatio), E(youngsModulus) { ElementMaterial(const float &poissonsRatio, const double &youngsModulus)
: poissonsRatio(poissonsRatio), youngsModulus(youngsModulus) {
assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
} }
ElementMaterial() : poissonsRatio(0.3), E(200 * 1e9) {} ElementMaterial() : poissonsRatio(0.3), youngsModulus(200) {}
}; };
#endif // BEAM_HPP #endif // BEAM_HPP

39
beamformfinder.cpp
View File

@ -43,7 +43,7 @@ void FormFinder::runUnitTests() {
settings.xi = 0.9969; settings.xi = 0.9969;
settings.maxDRMIterations = 0; settings.maxDRMIterations = 0;
settings.totalResidualForcesNormThreshold = 1; settings.totalResidualForcesNormThreshold = 1;
settings.shouldDraw = true; // settings.shouldDraw = true;
settings.beVerbose = true; settings.beVerbose = true;
SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation( SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(
std::make_shared<SimulationJob>(beamSimulationJob), settings); std::make_shared<SimulationJob>(beamSimulationJob), settings);
@ -752,7 +752,8 @@ 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) * element.properties.material.E * const double axialPE = pow(e_k, 2) *
element.properties.material.youngsModulus *
element.properties.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;
@ -763,17 +764,19 @@ double FormFinder::computeTotalPotentialEnergy() const {
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 = const double firstBendingPE = firstBendingPE_inBracketsTerm *
firstBendingPE_inBracketsTerm * element.properties.material.E * element.properties.material.youngsModulus *
element.properties.I2 / (2 * element.initialLength); element.properties.I2 /
(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;
const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) + const double secondBendingPE_inBracketsTerm = 4 * pow(theta3_j, 2) +
4 * theta3_j * theta3_jplus1 + 4 * theta3_j * theta3_jplus1 +
4 * pow(theta3_jplus1, 2); 4 * pow(theta3_jplus1, 2);
const double secondBendingPE = const double secondBendingPE = secondBendingPE_inBracketsTerm * 2 *
secondBendingPE_inBracketsTerm * 2 * element.properties.material.E * element.properties.material.youngsModulus *
element.properties.I3 / element.initialLength; element.properties.I3 /
element.initialLength;
totalInternalPotentialEnergy += totalInternalPotentialEnergy +=
axialPE + torsionalPE + firstBendingPE + secondBendingPE; axialPE + torsionalPE + firstBendingPE + secondBendingPE;
@ -1142,18 +1145,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 = const double SkTranslational = elem.properties.material.youngsModulus *
elem.properties.material.E * elem.properties.A / elem.length; elem.properties.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.E * elem.properties.I2 / elem.properties.material.youngsModulus * elem.properties.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.E * elem.properties.I3 / elem.properties.material.youngsModulus * elem.properties.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.E * elem.properties.J / elem.properties.material.youngsModulus * elem.properties.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;
@ -1450,18 +1453,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 = const double SkTranslational = elem.properties.material.youngsModulus *
elem.properties.material.E * elem.properties.A / elem.length; elem.properties.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.E * elem.properties.I2 / elem.properties.material.youngsModulus * elem.properties.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.E * elem.properties.I3 / elem.properties.material.youngsModulus * elem.properties.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.E * elem.properties.J / elem.properties.material.youngsModulus * elem.properties.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;

52
elementalmesh.cpp
View File

@ -282,12 +282,12 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>( customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32, nullptr); plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32, nullptr);
VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>> // VCGEdgeMesh::PerEdgeAttributeHandle<std::array<double, 6>>
handleBeamProperties = // handleBeamProperties =
vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute< // vcg::tri::Allocator<SimulationMesh>::AddPerEdgeAttribute<
std::array<double, 6>>(*this, plyPropertyBeamProperties); // std::array<double, 6>>(*this, plyPropertyBeamProperties);
customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>( // customAttrib.AddEdgeAttribDescriptor<std::array<double, 6>, double, 6>(
plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr); // plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64, nullptr);
// VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> // VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial>
// handleBeamRigidityContants; // handleBeamRigidityContants;
@ -313,11 +313,17 @@ bool SimulationMesh::loadPly(const string &plyFilename) {
initializeElements(); initializeElements();
updateEigenEdgeAndVertices(); updateEigenEdgeAndVertices();
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 = Element::Properties(handleBeamProperties[ei]); // elements[ei].properties =
elements[ei].updateRigidity(); // 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; return true;
@ -418,8 +424,7 @@ double computeAngle(const VectorType &vector0, const VectorType &vector1,
void Element::Properties::computeMaterialProperties( void Element::Properties::computeMaterialProperties(
const ElementMaterial &material) { const ElementMaterial &material) {
this->material = material; this->G = material.youngsModulus / (2 * (1 + material.poissonsRatio));
this->G = material.E / (2 * (1 + material.poissonsRatio));
} }
void Element::Properties::computeDimensionsProperties( void Element::Properties::computeDimensionsProperties(
@ -451,15 +456,14 @@ void Element::Properties::setMaterial(const ElementMaterial &material) {
} }
Element::Properties::Properties(const CrossSectionType &dimensions, Element::Properties::Properties(const CrossSectionType &dimensions,
const ElementMaterial &material) const ElementMaterial &material) {
: dimensions(dimensions), material(material) { setDimensions(dimensions);
computeDimensionsProperties(dimensions); setMaterial(material);
computeMaterialProperties(material);
} }
Element::Properties::Properties(const std::array<double, 6> &arr) { Element::Properties::Properties(const std::array<double, 6> &arr) {
assert(arr.size() == 6); assert(arr.size() == 6);
material.E = arr[0]; material.youngsModulus = arr[0];
G = arr[1]; G = arr[1];
A = arr[2]; A = arr[2];
I2 = arr[3]; I2 = arr[3];
@ -468,15 +472,17 @@ Element::Properties::Properties(const std::array<double, 6> &arr) {
} }
std::array<double, 6> Element::Properties::toArray() const { std::array<double, 6> Element::Properties::toArray() const {
return std::array<double, 6>( return std::array<double, 6>({material.youngsModulus, G, A,
{material.E, G, A, static_cast<double>(I2), static_cast<double>(I3), J}); static_cast<double>(I2),
static_cast<double>(I3), J});
} }
void Element::updateRigidity() { 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.torsional = properties.G * properties.J / initialLength;
rigidity.firstBending = rigidity.firstBending =
2 * properties.material.E * properties.I2 / initialLength; 2 * properties.material.youngsModulus * properties.I2 / initialLength;
rigidity.secondBending = rigidity.secondBending =
2 * properties.material.E * properties.I3 / initialLength; 2 * properties.material.youngsModulus * properties.I3 / initialLength;
} }