manolas
/
threed-beam-fea
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Added converged bool variable to the summary

This commit is contained in:
iasonmanolas 2021-05-01 21:33:10 +03:00
parent e49c911524
commit 13096dbd6c
2 changed files with 35 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
include/summary.h
View File

@ -120,6 +120,11 @@ struct Summary {
*/ */
unsigned long num_eqns; unsigned long num_eqns;
/**
* @brief true if the solution was aquired error-less false otherwise
*/
bool converged;
/** /**
* The resultant nodal displacement from the FE analysis. * The resultant nodal displacement from the FE analysis.
* `nodal_displacements` is a 2D vector where each row * `nodal_displacements` is a 2D vector where each row

43
src/threed_beam_fea.cpp
View File

@ -130,8 +130,10 @@ void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) {
// calculate unit normal vector along local x-direction // calculate unit normal vector along local x-direction
const Eigen::Vector3d nx = (job.nodes[nn2] - job.nodes[nn1]).normalized(); const Eigen::Vector3d nx = (job.nodes[nn2] - job.nodes[nn1]).normalized();
// calculate unit normal vector along y-direction // calculate unit normal vector along y-direction
const Eigen::Vector3d ny = job.props[i].normal_vec.normalized(); // const Eigen::Vector3d ny = job.props[i].normal_vec.normalized();
const Eigen::Vector3d ny = job.props[i].normal_vec.normalized().cross(nx).normalized();
// calculate the unit normal vector in local z direction // calculate the unit normal vector in local z direction
// const Eigen::Vector3d nz = nx.cross(ny).normalized();
const Eigen::Vector3d nz = nx.cross(ny).normalized(); const Eigen::Vector3d nz = nx.cross(ny).normalized();
RotationMatrix r; RotationMatrix r;
r.row(0) = nx; r.row(0) = nx;
@ -415,14 +417,14 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
loadForces(force_vec, forces); loadForces(force_vec, forces);
} }
#ifdef DEBUG #ifdef DEBUG
// Eigen::MatrixXd KgDense(Kg); Eigen::MatrixXd KgDense(Kg);
// std::cout << KgDense << std::endl; // std::cout << KgDense << std::endl;
// Eigen::VectorXd forcesVectorDense(force_vec); Eigen::VectorXd forcesVectorDense(force_vec);
// std::cout << forcesVectorDense << std::endl; // std::cout << forcesVectorDense << std::endl;
#endif #endif
// compress global stiffness matrix since all non-zero values have been added. // compress global stiffness matrix since all non-zero values have been added.
Kg.prune(1.e-14); Kg.prune(1.e-14);
Kg.makeCompressed(); Kg.makeCompressed();
#ifdef DEBUG #ifdef DEBUG
std::ofstream kgFile("Kg.csv"); std::ofstream kgFile("Kg.csv");
if (kgFile.is_open()) { if (kgFile.is_open()) {
@ -466,10 +468,11 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
// Compute the numerical factorization // Compute the numerical factorization
start_time = std::chrono::high_resolution_clock::now(); start_time = std::chrono::high_resolution_clock::now();
solver.factorize(Kg); solver.factorize(Kg);
if(options.verbose){ #ifdef DEBUG
std::cout << solver.lastErrorMessage() << std::endl; Eigen::VectorXd kgCol(Kg.col(7));
} Eigen::FullPivLU<Eigen::MatrixXd> kg_invCheck(Kg);
assert(solver.info() == Eigen::Success); assert(kg_invCheck.isInvertible());
#endif
end_time = std::chrono::high_resolution_clock::now(); end_time = std::chrono::high_resolution_clock::now();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
@ -481,6 +484,13 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
std::cout << "Factorization completed in " << delta_time std::cout << "Factorization completed in " << delta_time
<< " ms. Now solving system..." << std::endl; << " ms. Now solving system..." << std::endl;
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
assert(solver.info() == Eigen::Success);
return summary;
}
// Use the factors to solve the linear system // Use the factors to solve the linear system
start_time = std::chrono::high_resolution_clock::now(); start_time = std::chrono::high_resolution_clock::now();
SparseMat dispSparse = solver.solve(force_vec); SparseMat dispSparse = solver.solve(force_vec);
@ -493,6 +503,12 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
if (options.verbose) if (options.verbose)
std::cout << "System was solved in " << delta_time << " ms.\n" << std::endl; std::cout << "System was solved in " << delta_time << " ms.\n" << std::endl;
assert(solver.info() == Eigen::Success);
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
return summary;
}
// convert to dense matrix // convert to dense matrix
Eigen::VectorXd disp(dispSparse); Eigen::VectorXd disp(dispSparse);
@ -622,8 +638,7 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
Eigen::Matrix<double, 12, 1> elemForces = Eigen::Matrix<double, 12, 1> elemForces =
perElemKlocalAelem[elemIndex] * elemDisps; perElemKlocalAelem[elemIndex] * elemDisps;
for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) { for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) {
summary.element_forces[elemIndex][static_cast<size_t>(dofIndex)] = summary.element_forces[elemIndex][static_cast<size_t>(dofIndex)] = -elemForces(dofIndex);
-elemForces(dofIndex);
} }
// meaning of sign = reference to first node: // meaning of sign = reference to first node:
// + = compression for axial, - traction // + = compression for axial, - traction
@ -640,6 +655,8 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
options.csv_precision, options.csv_delimiter); options.csv_precision, options.csv_delimiter);
} }
summary.converged = true;
return summary; return summary;
}; };