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

Refactoring

This commit is contained in:
Iason 2020-12-09 17:59:18 +02:00
parent a67003cdd2
commit 5d8445fbca
8 changed files with 147 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
beam.hpp
View File

@ -12,21 +12,21 @@ struct RectangularBeamDimensions {
: b(width), h(height) { : b(width), h(height) {
assert(width > 0 && height > 0); assert(width > 0 && height > 0);
} }
RectangularBeamDimensions() : b(0.01), h(0.01) {} RectangularBeamDimensions() : b(1), h(1) {}
}; };
struct CylindricalElementDimensions { struct CylindricalBeamDimensions {
float od; // Cylinder outside diameter float od; // Cylinder outside diameter
float float
id; // Cylinder inside diameter id; // Cylinder inside diameter
// https://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html // https://www.engineeringtoolbox.com/area-moment-inertia-d_1328.html
CylindricalElementDimensions(const float &outsideDiameter, CylindricalBeamDimensions(const float &outsideDiameter,
const float &insideDiameter) const float &insideDiameter)
: od(outsideDiameter), id(insideDiameter) { : od(outsideDiameter), id(insideDiameter) {
assert(outsideDiameter > 0 && insideDiameter > 0 && assert(outsideDiameter > 0 && insideDiameter > 0 &&
outsideDiameter > insideDiameter); outsideDiameter > insideDiameter);
} }
CylindricalElementDimensions() : od(0.03), id(0.026) {} CylindricalBeamDimensions() : od(0.03), id(0.026) {}
}; };
struct ElementMaterial { struct ElementMaterial {
@ -37,21 +37,7 @@ struct ElementMaterial {
youngsModulusGPascal(youngsModulusGPascal) { youngsModulusGPascal(youngsModulusGPascal) {
assert(poissonsRatio <= 0.5 && poissonsRatio >= -1); assert(poissonsRatio <= 0.5 && poissonsRatio >= -1);
} }
ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(7.5) {} ElementMaterial() : poissonsRatio(0.3), youngsModulusGPascal(200) {}
};
struct BeamProperties {
float crossArea;
float I2;
float I3;
float polarInertia;
BeamProperties(const RectangularBeamDimensions &dimensions,
const ElementMaterial &material) {
crossArea = (dimensions.b * dimensions.h);
I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
polarInertia = (I2 + I3);
}
}; };
#endif // BEAM_HPP #endif // BEAM_HPP

262
beamformfinder.cpp
View File

@ -51,7 +51,7 @@ VectorType FormFinder::computeDerivativeOfNormal(
const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2);
if (dui.dofi == 3) { if (dui.dofi == 3) {
if (nxnyMagnitude >= 1) { if (nxnyMagnitude + 1e-5 >= 1) {
const double normalDerivativeX = const double normalDerivativeX =
1 / sqrt(nxnyMagnitude) - 1 / sqrt(nxnyMagnitude) -
std::pow(nx, 2) / std::pow(nxnyMagnitude, 1.5); std::pow(nx, 2) / std::pow(nxnyMagnitude, 1.5);
@ -67,7 +67,7 @@ VectorType FormFinder::computeDerivativeOfNormal(
VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ); VectorType(normalDerivativeX, normalDerivativeY, normalDerivativeZ);
} }
} else if (dui.dofi == 4) { } else if (dui.dofi == 4) {
if (nxnyMagnitude >= 1) { if (nxnyMagnitude + 1e-5 >= 1) {
const double normalDerivativeX = -nx * ny / std::pow(nxnyMagnitude, 1.5); const double normalDerivativeX = -nx * ny / std::pow(nxnyMagnitude, 1.5);
const double normalDerivativeY = const double normalDerivativeY =
1 / sqrt(nxnyMagnitude) - 1 / sqrt(nxnyMagnitude) -
@ -88,6 +88,8 @@ VectorType FormFinder::computeDerivativeOfNormal(
std::isfinite(normalDerivative[1]) && std::isfinite(normalDerivative[1]) &&
std::isfinite(normalDerivative[2]); std::isfinite(normalDerivative[2]);
assert(normalDerivativeIsFinite); assert(normalDerivativeIsFinite);
const bool shouldBreak =
mCurrentSimulationStep == 118 && vi == 1 && dui.dofi == 3;
return normalDerivative; return normalDerivative;
} }
@ -223,6 +225,8 @@ FormFinder::computeDerivativeT2(const EdgeType &e,
const double t2DerivNorm = t2Deriv.Norm(); const double t2DerivNorm = t2Deriv.Norm();
assert(std::isfinite(t2DerivNorm)); assert(std::isfinite(t2DerivNorm));
const bool shouldBreak = mCurrentSimulationStep == 118 &&
(vi_jplus1 == 1 || vi_j == 1) && dofi == 3;
return t2Deriv; return t2Deriv;
} }
@ -297,6 +301,8 @@ double FormFinder::computeDerivativeTheta1(const EdgeType &e,
const double theta1Derivative = const double theta1Derivative =
derivativeT1 * Cross(t3, n) + derivativeT1 * Cross(t3, n) +
t1 * (Cross(derivativeT3, n) + Cross(t3, nDerivative)); t1 * (Cross(derivativeT3, n) + Cross(t3, nDerivative));
const bool shouldBreak =
mCurrentSimulationStep == 118 && vi == 1 && dwrt_dofi == 3;
return theta1Derivative; return theta1Derivative;
} }
@ -703,6 +709,8 @@ void FormFinder::updateResidualForcesOnTheFly(
secondBendingForce_inBracketsTerm_3; secondBendingForce_inBracketsTerm_3;
double secondBendingForce_dofi = secondBendingForce_inBracketsTerm * double secondBendingForce_dofi = secondBendingForce_inBracketsTerm *
element.secondBendingConstFactor; element.secondBendingConstFactor;
const bool shouldBreak =
mCurrentSimulationStep == 118 && edgeNode.vi == 1 && dofi == 3;
internalForcesContributionFromThisEdge[evi].second[dofi] = internalForcesContributionFromThisEdge[evi].second[dofi] =
axialForce_dofi + firstBendingForce_dofi + axialForce_dofi + firstBendingForce_dofi +
secondBendingForce_dofi + torsionalForce_dofi; secondBendingForce_dofi + torsionalForce_dofi;
@ -1024,7 +1032,7 @@ void FormFinder::updateNodePosition(
node.previousLocation = previousLocation; node.previousLocation = previousLocation;
v.P() = node.initialLocation + displacementVector; v.P() = node.initialLocation + displacementVector;
if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) { if (shouldApplyInitialDistortion && mCurrentSimulationStep < 40) {
VectorType desiredInitialDisplacement(0, 0, 0.1); VectorType desiredInitialDisplacement(0, 0, 0.001);
v.P() = v.P() + desiredInitialDisplacement; v.P() = v.P() + desiredInitialDisplacement;
} }
} }
@ -1033,8 +1041,14 @@ void FormFinder::updateNodeNormal(
VertexType &v, VertexType &v,
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices) { &fixedVertices) {
Node &node = mesh->nodes[v]; Node &node = mesh->nodes[v];
const VertexIndex &vi = node.vi; const VertexIndex &vi = node.vi;
// if (vi == 1) {
// std::cout << "PRE:" << mesh->vert[1].N()[0] << " " <<
// mesh->vert[1].N()[1]
// << " " << mesh->vert[1].N()[2] << std::endl;
// }
VectorType normalDisplacementVector(0, 0, 0); VectorType normalDisplacementVector(0, 0, 0);
if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(3)) { if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(3)) {
normalDisplacementVector += VectorType(node.displacements[3], 0, 0); normalDisplacementVector += VectorType(node.displacements[3], 0, 0);
@ -1042,9 +1056,12 @@ void FormFinder::updateNodeNormal(
if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(4)) { if (!fixedVertices.contains(vi) || !fixedVertices.at(vi).contains(4)) {
normalDisplacementVector += VectorType(0, node.displacements[4], 0); normalDisplacementVector += VectorType(0, node.displacements[4], 0);
} }
const double nxnyMagnitudePre = std::pow(v.N()[0], 2) + std::pow(v.N()[1], 2);
v.N() = node.initialNormal + normalDisplacementVector; v.N() = node.initialNormal + normalDisplacementVector;
const double &nx = v.N()[0], ny = v.N()[1]; const double &nx = v.N()[0], ny = v.N()[1], nz = v.N()[2];
const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2); const double nxnyMagnitude = std::pow(nx, 2) + std::pow(ny, 2);
VectorType n = v.N();
const bool shouldBreak = mCurrentSimulationStep == 118 && vi == 1;
if (nxnyMagnitude > 1) { if (nxnyMagnitude > 1) {
VectorType newNormal(nx / std::sqrt(nxnyMagnitude), VectorType newNormal(nx / std::sqrt(nxnyMagnitude),
ny / std::sqrt(nxnyMagnitude), 0); ny / std::sqrt(nxnyMagnitude), 0);
@ -1067,6 +1084,11 @@ void FormFinder::updateNodeNormal(
VectorType g1 = Cross(refT1, t1Initial); VectorType g1 = Cross(refT1, t1Initial);
node.nR = g1.dot(v.cN()); node.nR = g1.dot(v.cN());
} }
// if (vi == 1) {
// std::cout << "AFTER:" << mesh->vert[1].N()[0] << " " <<
// mesh->vert[1].N()[1]
// << " " << mesh->vert[1].N()[2] << std::endl;
// }
} }
void FormFinder::applyDisplacements( void FormFinder::applyDisplacements(
@ -1103,6 +1125,25 @@ void FormFinder::applyForcedDisplacements(
VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1), VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1),
vertexDisplacement(2)); vertexDisplacement(2));
mesh->vert[vi].P() = node.initialLocation + displacementVector; 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]});
}
}
void FormFinder::applyForcedNormals(
const std::unordered_map<VertexIndex, VectorType> nodalForcedRotations) {
for (const std::pair<VertexIndex, VectorType> vertexIndexDesiredNormalPair :
nodalForcedRotations) {
const VertexIndex vi = vertexIndexDesiredNormalPair.first;
Node &node = mesh->nodes[vi];
mesh->vert[vi].N() = vertexIndexDesiredNormalPair.second;
node.displacements = Vector6d(
{node.displacements[0], node.displacements[1], node.displacements[2],
vertexIndexDesiredNormalPair.second[0] - node.initialNormal[0],
vertexIndexDesiredNormalPair.second[1] - node.initialNormal[1],
node.displacements[5]});
} }
} }
@ -1547,12 +1588,14 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
constrainedVertices[vi].insert({0, 1, 2}); constrainedVertices[vi].insert({0, 1, 2});
} }
} }
if (!job.nodalForcedNormals.empty()) {
for (std::pair<VertexIndex, Eigen::Vector3d> viNormalPair :
job.nodalForcedDisplacements) {
assert(viNormalPair.second[2] >= 0);
}
}
// VCGEdgeMesh *jobEdgeMesh = job.mesh.get();
mesh = std::make_unique<SimulationMesh>(*job.mesh); mesh = std::make_unique<SimulationMesh>(*job.mesh);
// ElementalMesh *jobElementalMesh = job.mesh.get();
// vcg::tri::Append<ElementalMesh, ConstElementalMesh>::MeshCopy(
// *(this->mesh), *jobElementalMesh);
if (beVerbose) { if (beVerbose) {
std::cout << "Executing simulation for mesh with:" << mesh->VN() std::cout << "Executing simulation for mesh with:" << mesh->VN()
<< " nodes and " << mesh->EN() << " elements." << std::endl; << " nodes and " << mesh->EN() << " elements." << std::endl;
@ -1564,75 +1607,21 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
if (!polyscope::state::initialized) { if (!polyscope::state::initialized) {
initPolyscope(); initPolyscope();
} }
if (!polyscope::hasCurveNetwork(deformedMeshName)) { polyscope::registerCurveNetwork(deformedMeshName, mesh->getEigenVertices(),
polyscope::registerCurveNetwork( mesh->getEigenEdges());
deformedMeshName, mesh->getEigenVertices(), mesh->getEigenEdges());
}
// registerWorldAxes(); // registerWorldAxes();
} }
// const double beamRadius = mesh.getBeamDimensions()[0].od / 2;
// std::cout << "BeamRadius:" << beamRadius << std::endl;
// polyscope::getCurveNetwork("Deformed edge mesh")
// //
// // ->setRadius(beamRadius);
// ->setRadius(0.01);
// const bool hasExternalLoads =
// !job.nodalExternalForces.empty() ||
// !job.nodalForcedDisplacements.empty();
// assert(hasExternalLoads);
for (auto fixedVertex : job.fixedVertices) { for (auto fixedVertex : job.fixedVertices) {
assert(fixedVertex.first < mesh->VN()); assert(fixedVertex.first < mesh->VN());
} }
updateElementalFrames(); updateElementalFrames();
updateNodalMasses(); updateNodalMasses();
if (job.nodalExternalForces.empty()) { if (!job.nodalForcedDisplacements.empty() &&
job.nodalExternalForces.empty()) {
shouldApplyInitialDistortion = true; shouldApplyInitialDistortion = true;
} }
// std::unordered_map<VertexIndex, Eigen::Vector3d> temporaryForces{
// // // {2, Eigen::Vector3d(0, 0, 1000)},
// // // {12, Eigen::Vector3d(0, 0, 1000)},
// // // {18, Eigen::Vector3d(0, 0, 1000)},
// // // {8, Eigen::Vector3d(0, 0, 1000)}};
// {10, Eigen::Vector3d(0, 0, 10000)}};
// std::unordered_map<VertexIndex, Eigen::Vector3d> temporaryForces;
// for (VertexIndex vi = 0; vi < mesh.VN(); vi++) {
// for (VertexType &v : mesh.vert) {
// const VertexIndex vi = mesh.getIndex(v);
// VectorType allowedDoFType(1, 1, 1);
// if (constrainedVertices.contains(vi)) {
// std::unordered_set<gsl::index> constrainedDof =
// constrainedVertices.at(vi);
// if (constrainedDof.contains(0)) {
// allowedDoFType[0] = 0;
// } else if (constrainedDof.contains(1)) {
// allowedDoFType[1] = 0;
// } else if (constrainedDof.contains(2)) {
// allowedDoFType[2] = 0;
// }
// }
// VectorType desiredDisplacement = VectorType(0, 0, 0.2);
// VectorType displacement(allowedDoFType[0] * desiredDisplacement[0],
// allowedDoFType[1] * desiredDisplacement[1],
// allowedDoFType[2] * desiredDisplacement[2]);
// v.P() = v.P() + displacement;
// temporaryForces.insert({vi, Eigen::Vector3d(0, 0, 100000)});
// }
// updateNodalExternalForces(temporaryForces, job.fixedVertices);
// appliedTemporaryForce = true;
// } else {
// std::unordered_map<VertexIndex, Eigen::Vector3d> appliedLoad =
// job.nodalExternalForces;
// const size_t numberOfStepsForApplyingExternalLoads = 10;
// int externalLoadStep = 1;
// std::for_each(appliedLoad.begin(), appliedLoad.end(), [](auto &pair) {
// pair.second /= numberOfStepsForApplyingExternalLoads;
// });
// updateNodalExternalForces(appliedLoad, constrainedVertices);
updateNodalExternalForces(job.nodalExternalForces, constrainedVertices); updateNodalExternalForces(job.nodalExternalForces, constrainedVertices);
// }
while (mCurrentSimulationStep < maxDRMIterations) { while (mCurrentSimulationStep < maxDRMIterations) {
// while (true) { // while (true) {
@ -1641,36 +1630,7 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
updateRDerivatives(); updateRDerivatives();
updateT2Derivatives(); updateT2Derivatives();
updateT3Derivatives(); updateT3Derivatives();
// updateRotationalDisplacements();
// if (currentSimulationStep > lastPulseStepIndex) {
// const std::vector<Vector6d> internalForces =
updateResidualForcesOnTheFly(constrainedVertices); updateResidualForcesOnTheFly(constrainedVertices);
//#pragma omp parallel for schedule(static) num_threads(8)
// double totalResidualForcesNorm = 0;
// for (size_t vi = 0; vi < mesh.VN(); vi++) {
// Node::Forces &force = mesh.nodes[vi].force;
// // const Vector6d &parallelForce = internalForcesParallel[vi];
// // const Vector6d &serialForce = internalForces[vi];
// // const double normOfDifference = (serialForce +
// (parallelForce
// *
// // -1)).norm(); assert(normOfDifference < 0.000001);
// force.residual = force.external - internalForces[vi];
// const double residualForceNorm = (force.residual).norm();
// totalResidualForcesNorm += residualForceNorm;
// // assert(residualForceNorm < std::pow(10, 20));
// }
// mesh.totalResidualForcesNorm = totalResidualForcesNorm;
// } else {
// if (currentSimulationStep == lastPulseStepIndex &&
// appliedTemporaryForce) {
// for (const VertexType &v : mesh.vert) {
// Node &node = mesh.nodes[v];
// node.force.external = std::optional<Vector6d>();
// }
// }
// updateNodalInternalForces(job.fixedVertices);
// }
// TODO: write parallel function for updating positions // TODO: write parallel function for updating positions
// TODO: Make a single function out of updateResidualForcesOnTheFly // TODO: Make a single function out of updateResidualForcesOnTheFly
@ -1686,16 +1646,11 @@ FormFinder::executeSimulation(const SimulationJob &job, const bool &beVerbose,
job.nodalForcedDisplacements); job.nodalForcedDisplacements);
} }
if (!job.nodalForcedNormals.empty()) {
applyForcedNormals(job.nodalForcedNormals);
}
updateElementalLengths(); updateElementalLengths();
// if (!std::isfinite(mesh.currentTotalKineticEnergy)) {
// std::cerr << "Infinite kinetic energy. Breaking simulation.."
// << std::endl;
// break;
// }
// assert(std::isfinite(mesh.currentTotalKineticEnergy));
// mesh.printVertexCoordinates(mesh.VN() / 2);
// draw();
if (mShouldDraw && if (mShouldDraw &&
mCurrentSimulationStep % mDrawingStep == 0 /* && mCurrentSimulationStep % mDrawingStep == 0 /* &&
currentSimulationStep > maxDRMIterations*/) { currentSimulationStep > maxDRMIterations*/) {
@ -1704,22 +1659,8 @@ currentSimulationStep > maxDRMIterations*/) {
// .append("Screenshots") // .append("Screenshots")
// .string(); // .string();
// draw(saveTo); // draw(saveTo);
// SimulationResultsReporter::createPlot( std::cout << "Residual forces norm: " << mesh->totalResidualForcesNorm
// "Number of iterations", "Log of Kinetic Energy", << std::endl;
// history.kineticEnergy,
// std::filesystem::path(saveTo).append(
// std::to_string(currentSimulationStep) + "_kin.png"));
// draw();
// auto t2 = std::chrono::high_resolution_clock::now();
// auto timePerNodePerIteration =
// std::chrono::duration_cast<std::chrono::microseconds>(t2 -
// t1)
// .count() /
// (mesh.VN() * (currentSimulationStep + 1));
// std::cout << "Time/(node*iteration) "
// << timePerNodePerIteration * std::pow(10, -6) << "s"
// << std::endl;
draw(); draw();
} }
if (mCurrentSimulationStep != 0) { if (mCurrentSimulationStep != 0) {
@ -1745,45 +1686,11 @@ currentSimulationStep > maxDRMIterations*/) {
std::cout << "Total potential:" << mesh->totalPotentialEnergykN std::cout << "Total potential:" << mesh->totalPotentialEnergykN
<< " kNm" << std::endl; << " kNm" << std::endl;
} }
// if (externalLoadStep !=
// numberOfStepsForApplyingExternalLoads)
// {
// std::for_each(
// appliedLoad.begin(), appliedLoad.end(), [&](auto
// &pair)
// {
// pair.second +=
// job.nodalExternalForces.at(pair.first)
// /
// numberOfStepsForApplyingExternalLoads;
// });
// updateNodalExternalForces(appliedLoad,
// constrainedVertices);
// } else {
break; break;
// } // }
} }
// history.markRed(currentSimulationStep);
// std::cout << "Total potential:" << mesh.totalPotentialEnergykN
// << " kNm"
// << std::endl;
// reset displacements to previous position where the peak occured
// for (VertexType &v : mesh->vert) {
// Node &node = mesh->nodes[v];
// node.displacements = node.displacements - node.velocity * Dt;
// }
// applyDisplacements(constrainedVertices, shouldDraw);
// if (!job.nodalForcedDisplacements.empty()) {
// applyForcedDisplacements(
// job.nodalForcedDisplacements);
// }
// updateElementalLengths();
resetVelocities(); resetVelocities();
Dt = Dt * xi; Dt = Dt * xi;
// std::cout << "Residual forces norm:" <<
// mesh.totalResidualForcesNorm
// << std::endl;
} }
} }
if (mCurrentSimulationStep == maxDRMIterations) { if (mCurrentSimulationStep == maxDRMIterations) {
@ -1809,19 +1716,42 @@ currentSimulationStep > maxDRMIterations*/) {
polyscope::removeCurveNetwork(deformedMeshName); polyscope::removeCurveNetwork(deformedMeshName);
// } // }
} }
// debugger.createPlots();
SimulationResults results = computeResults(*job.mesh); SimulationResults results = computeResults(*job.mesh);
// Eigen::write_binary("optimizedResult.eigenBin",
// results.nodalDisplacements);
// const std::string groundOfTruthBinaryFilename =
// "../../groundOfTruths/grid6_groundOfTruth.eigenBin";
// // "../../groundOfTruths/longSpanGridshell_groundOfTruth.eigenBin";
// assert(std::filesystem::exists(
// std::filesystem::path(groundOfTruthBinaryFilename)));
// Eigen::MatrixX3d groundOfTruthMatrix;
// Eigen::read_binary(groundOfTruthBinaryFilename, groundOfTruthMatrix);
// assert(results.nodalDisplacements.isApprox(groundOfTruthMatrix));
return results; return results;
// return history; }
void FormFinder::runUnitTests() {
FormFinder formFinder;
VCGEdgeMesh mesh;
// const size_t spanGridSize = 11;
// mesh.createSpanGrid(spanGridSize);
mesh.loadFromPly("/home/iason/Models/simple_beam_paper_example.ply");
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3};
fixedVertices[mesh.VN() - 1] = std::unordered_set<DoFType>{1, 2};
std::unordered_map<VertexIndex, Vector6d> nodalForces{
{2, Vector6d({0, 1000, 1000, 0, 0, 0})}};
// Forced displacements
std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
nodalForcedDisplacements.insert({mesh.VN() - 1, Eigen::Vector3d(-0.2, 0, 0)});
SimulationJob beamSimulationJob{
std::make_shared<SimulationMesh>(mesh),
// SimulationJob::constructFixedVerticesSpanGrid(spanGridSize,
// mesh.VN()),
fixedVertices, nodalForces, nodalForcedDisplacements};
beamSimulationJob.mesh->setBeamMaterial(0.3, 200);
assert(typeid(CrossSectionType) == typeid(CylindricalBeamDimensions));
beamSimulationJob.mesh->setBeamCrossSection(CrossSectionType{0.03, 0.026});
SimulationResults beamSimulationResults =
formFinder.executeSimulation(beamSimulationJob, true, false);
const bool testSuccessful =
beamSimulationResults.displacements[2][0] - (-0.09556) < 1e-5 &&
beamSimulationResults.displacements[2][1] - 0.40666 < 1e-5 &&
beamSimulationResults.displacements[2][2] - 0.40666 < 1e-5;
assert(testSuccessful);
beamSimulationResults.simulationLabel = "Beam";
beamSimulationResults.registerForDrawing(beamSimulationJob);
polyscope::show();
} }

9
beamformfinder.hpp
View File

@ -27,7 +27,7 @@ private:
const double Dtini{0.1}; const double Dtini{0.1};
double Dt{Dtini}; double Dt{Dtini};
const double xi{0.9969}; const double xi{0.9969};
const double totalResidualForcesNormThreshold{10}; const double totalResidualForcesNormThreshold{0.01};
size_t mCurrentSimulationStep{0}; size_t mCurrentSimulationStep{0};
int mDrawingStep{5}; int mDrawingStep{5};
bool mShouldDraw{false}; bool mShouldDraw{false};
@ -176,13 +176,18 @@ private:
const std::unordered_map<VertexIndex, std::unordered_set<DoFType>> const std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
&fixedVertices); &fixedVertices);
void applyForcedNormals(
const std::unordered_map<VertexIndex, VectorType> nodalForcedRotations);
public: public:
FormFinder(); FormFinder();
SimulationResults executeSimulation( SimulationResults executeSimulation(
const SimulationJob &job, const bool &beVerbose = false, const SimulationJob &job, const bool &beVerbose = false,
const bool &shouldDraw = false, const bool &shouldDraw = false,
const size_t &maxDRMIterations = FormFinder::maxDRMIterations); const size_t &maxDRMIterations = FormFinder::maxDRMIterations);
inline static const size_t maxDRMIterations{50000}; inline static const size_t maxDRMIterations{500000};
static void runUnitTests();
}; };
template <typename PointType> PointType Cross(PointType p1, PointType p2) { template <typename PointType> PointType Cross(PointType p1, PointType p2) {

41
elementalmesh.cpp
View File

@ -2,6 +2,21 @@
SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) { SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
vcg::tri::MeshAssert<VCGEdgeMesh>::VertexNormalNormalized(mesh); vcg::tri::MeshAssert<VCGEdgeMesh>::VertexNormalNormalized(mesh);
// bool containsNormals = true;
// for (VertexIterator vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi)
// if (!vi->IsD()) {
// if (fabs(vi->cN().Norm() - 1.0) > 0.000001) {
// containsNormals = false;
// break;
// }
// }
// if (!containsNormals) {
// for (VertexIterator vi = mesh.vert.begin(); vi != mesh.vert.end(); ++vi)
// if (!vi->IsD()) {
// vi->N() = CoordType(1, 0, 0);
// }
// }
vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, mesh); vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, mesh);
elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>( elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
@ -51,8 +66,7 @@ SimulationMesh::SimulationMesh(SimulationMesh &mesh) {
} }
void SimulationMesh::computeElementalProperties() { void SimulationMesh::computeElementalProperties() {
const std::vector<RectangularBeamDimensions> elementalDimensions = const std::vector<CrossSectionType> elementalDimensions = getBeamDimensions();
getBeamDimensions();
const std::vector<ElementMaterial> elementalMaterials = getBeamMaterial(); const std::vector<ElementMaterial> elementalMaterials = getBeamMaterial();
assert(EN() == elementalDimensions.size() && assert(EN() == elementalDimensions.size() &&
elementalDimensions.size() == elementalMaterials.size()); elementalDimensions.size() == elementalMaterials.size());
@ -118,7 +132,7 @@ 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(1, 1); element.properties.dimensions = CrossSectionType();
// Initialize material // Initialize material
element.properties.material = ElementMaterial(); element.properties.material = ElementMaterial();
// Initialize lengths // Initialize lengths
@ -162,11 +176,11 @@ void SimulationMesh::updateElementalLengths() {
} }
} }
void SimulationMesh::setBeamCrossSection(const double &firstDimension, void SimulationMesh::setBeamCrossSection(
const double &secondDimension) { const CrossSectionType &beamDimensions) {
for (size_t ei = 0; ei < EN(); ei++) { for (size_t ei = 0; ei < EN(); ei++) {
elements[ei].properties.setDimensions( elements[ei].properties.dimensions = beamDimensions;
CrossSectionType{firstDimension, secondDimension}); elements[ei].properties.computeDimensionsProperties(beamDimensions);
elements[ei].updateConstFactors(); elements[ei].updateConstFactors();
} }
} }
@ -178,8 +192,8 @@ void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {
} }
} }
std::vector<RectangularBeamDimensions> SimulationMesh::getBeamDimensions() { std::vector<CrossSectionType> SimulationMesh::getBeamDimensions() {
std::vector<RectangularBeamDimensions> 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].properties.dimensions;
} }
@ -198,8 +212,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) {
this->Clear(); this->Clear();
// assert(plyFileHasAllRequiredFields(plyFilename)); // assert(plyFileHasAllRequiredFields(plyFilename));
// Load the ply file // Load the ply file
VCGEdgeMesh::PerEdgeAttributeHandle<RectangularBeamDimensions> VCGEdgeMesh::PerEdgeAttributeHandle<CrossSectionType> handleBeamDimensions;
handleBeamDimensions;
VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial; VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial;
handleBeamDimensions = handleBeamDimensions =
vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<CrossSectionType>( vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<CrossSectionType>(
@ -209,7 +222,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) {
*this, plyPropertyBeamMaterialID); *this, plyPropertyBeamMaterialID);
nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib; nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
customAttrib.GetMeshAttrib(plyFilename); customAttrib.GetMeshAttrib(plyFilename);
customAttrib.AddEdgeAttribDescriptor<RectangularBeamDimensions, float, 2>( customAttrib.AddEdgeAttribDescriptor<CrossSectionType, float, 2>(
plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>( customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr); plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
@ -228,7 +241,7 @@ bool SimulationMesh::loadFromPly(const string &plyFilename) {
bool SimulationMesh::savePly(const std::string &plyFilename) { bool SimulationMesh::savePly(const std::string &plyFilename) {
nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib; nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
customAttrib.GetMeshAttrib(plyFilename); customAttrib.GetMeshAttrib(plyFilename);
customAttrib.AddEdgeAttribDescriptor<RectangularBeamDimensions, float, 2>( customAttrib.AddEdgeAttribDescriptor<CrossSectionType, float, 2>(
plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32,
getBeamDimensions().data()); getBeamDimensions().data());
customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>( customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
@ -318,7 +331,7 @@ void Element::Properties::computeDimensionsProperties(
} }
void Element::Properties::computeDimensionsProperties( void Element::Properties::computeDimensionsProperties(
const CylindricalElementDimensions &dimensions) { const CylindricalBeamDimensions &dimensions) {
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));
I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64; I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
I3 = I2; I3 = I2;

8
elementalmesh.hpp
View File

@ -8,6 +8,7 @@
struct Element; struct Element;
struct Node; struct Node;
using CrossSectionType = RectangularBeamDimensions; using CrossSectionType = RectangularBeamDimensions;
// using CrossSectionType = CylindricalBeamDimensions;
class SimulationMesh : public VCGEdgeMesh { class SimulationMesh : public VCGEdgeMesh {
private: private:
@ -26,9 +27,6 @@ public:
SimulationMesh(ConstVCGEdgeMesh &edgeMesh); SimulationMesh(ConstVCGEdgeMesh &edgeMesh);
SimulationMesh(SimulationMesh &elementalMesh); SimulationMesh(SimulationMesh &elementalMesh);
void updateElementalLengths(); void updateElementalLengths();
void setBeamCrossSection(const double &firstDimension,
const double &secondDimension);
void setBeamMaterial(const double &pr, const double &ym);
std::vector<VCGEdgeMesh::EdgePointer> std::vector<VCGEdgeMesh::EdgePointer>
getIncidentElements(const VCGEdgeMesh::VertexType &v); getIncidentElements(const VCGEdgeMesh::VertexType &v);
@ -42,6 +40,8 @@ public:
double totalResidualForcesNorm{0}; double totalResidualForcesNorm{0};
double totalPotentialEnergykN{0}; double totalPotentialEnergykN{0};
bool savePly(const std::string &plyFilename); bool savePly(const std::string &plyFilename);
void setBeamCrossSection(const CrossSectionType &beamDimensions);
void setBeamMaterial(const double &pr, const double &ym);
}; };
struct Element { struct Element {
@ -58,7 +58,7 @@ struct Element {
void void
computeDimensionsProperties(const RectangularBeamDimensions &dimensions); computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
void void
computeDimensionsProperties(const CylindricalElementDimensions &dimensions); 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, Properties(const CrossSectionType &dimensions,

2
flatpattern.hpp
View File

@ -27,7 +27,7 @@ private:
void removeDuplicateVertices(); void removeDuplicateVertices();
void scale(); void scale();
const double desiredBaseTriangleCentralEdgeSize{0.025}; const double desiredBaseTriangleCentralEdgeSize{0.03};
}; };
#endif // FLATPATTERN_HPP #endif // FLATPATTERN_HPP

1
simulationresult.hpp
View File

@ -36,6 +36,7 @@ struct SimulationJob {
std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices; std::unordered_map<VertexIndex, std::unordered_set<int>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalExternalForces; std::unordered_map<VertexIndex, Vector6d> nodalExternalForces;
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements; std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
std::unordered_map<VertexIndex, VectorType> nodalForcedNormals;
void registerForDrawing() const { void registerForDrawing() const {
initPolyscope(); initPolyscope();

8
utilities.hpp
View File

@ -237,11 +237,13 @@ inline void registerWorldAxes() {
if (!polyscope::state::initialized) { if (!polyscope::state::initialized) {
initPolyscope(); initPolyscope();
} }
Eigen::MatrixX3d axesPositions(4, 3); Eigen::MatrixX3d axesPositions(4, 3);
axesPositions.row(0) = Eigen::Vector3d(0, 0, 0); axesPositions.row(0) = Eigen::Vector3d(0, 0, 0);
axesPositions.row(1) = Eigen::Vector3d(1, 0, 0); axesPositions.row(1) = Eigen::Vector3d(polyscope::state::lengthScale, 0, 0);
axesPositions.row(2) = Eigen::Vector3d(0, 1, 0); axesPositions.row(2) = Eigen::Vector3d(0, polyscope::state::lengthScale, 0);
axesPositions.row(3) = Eigen::Vector3d(0, 0, 1); axesPositions.row(3) = Eigen::Vector3d(0, 0, polyscope::state::lengthScale);
Eigen::MatrixX2i axesEdges(3, 2); Eigen::MatrixX2i axesEdges(3, 2);
axesEdges.row(0) = Eigen::Vector2i(0, 1); axesEdges.row(0) = Eigen::Vector2i(0, 1);
axesEdges.row(1) = Eigen::Vector2i(0, 2); axesEdges.row(1) = Eigen::Vector2i(0, 2);