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

Changed Element properties

This commit is contained in:
Iason 2021-01-14 15:39:19 +02:00
parent d933f92eeb
commit 2e7aa6e1e0
5 changed files with 82 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
beam.hpp
View File

@ -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

86
beamformfinder.cpp
View File

@ -8,8 +8,6 @@
#include <execution>
#include <omp.h>
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<SimulationJob>(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<SimulationJob> &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<SimulationJob> &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();
}
}

3
beamformfinder.hpp
View File

@ -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<VertexIndex, VectorType> nodalForcedRotations);
void printCurrentState();
void printCurrentState() const;
void printDebugInfo() const;

47
elementalmesh.cpp
View File

@ -335,18 +335,18 @@ bool SimulationMesh::savePly(const std::string &plyFilename) {
customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
plyPropertyBeamMaterialID, nanoply::NNP_LIST_INT8_FLOAT32,
material.data());
std::vector<std::array<double, 6>> 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<std::array<double, 6>, double, 6>(
plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64,
beamProperties.data());
// std::vector<std::array<double, 6>> 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<std::array<double, 6>, double, 6>(
// plyPropertyBeamProperties, nanoply::NNP_LIST_INT8_FLOAT64,
// beamProperties.data());
// Load the ply file
unsigned int mask = 0;
mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::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<double, 6> &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<double, 6> &arr) {
}
std::array<double, 6> Element::Properties::toArray() const {
return std::array<double, 6>({E, G, A, I2, I3, J});
return std::array<double, 6>({material.E, material.G, A,
static_cast<double>(I2),
static_cast<double>(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;
}

6
elementalmesh.hpp
View File

@ -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