diff --git a/include/options.h b/include/options.h index 56b2fc0..afbbadc 100644 --- a/include/options.h +++ b/include/options.h @@ -29,104 +29,124 @@ #include namespace fea { - /** - * @brief Provides a method for customizing the finite element analysis. - */ - struct Options { - /** - * @brief Default constructor - * @details This tries to set up reasonable defaults for analysis. It is - * recommended that the user read and overwrite the defaults based on the information - * desired from the analysis. - */ - Options() { - epsilon = 1e-14; +/** + * @brief Provides a method for customizing the finite element analysis. + */ +struct Options { + /** + * @brief Default constructor + * @details This tries to set up reasonable defaults for analysis. It is + * recommended that the user read and overwrite the defaults based on the + * information desired from the analysis. + */ + Options() { + epsilon = 1e-14; - csv_precision = 14; - csv_delimiter = ","; + csv_precision = 14; + csv_delimiter = ","; - save_nodal_displacements = false; - save_nodal_forces = false; - save_tie_forces = false; - verbose = false; - save_report = false; + save_nodal_displacements = false; + save_nodal_forces = false; + save_tie_forces = false; + verbose = false; + save_report = false; - nodal_displacements_filename = "nodal_displacements.csv"; - nodal_forces_filename = "nodal_forces.csv"; - tie_forces_filename = "tie_forces.csv"; - report_filename = "report.txt"; - } + nodal_displacements_filename = "nodal_displacements.csv"; + nodal_forces_filename = "nodal_forces.csv"; + tie_forces_filename = "tie_forces.csv"; + report_filename = "report.txt"; + } - /** - * Values of forces and nodal displacements which have a magnitude less than `epsilon` will be rounded to 0.0. - * Default = `1e-14`. This is a simple way to deal with machine precision when doing calculations. - */ - double epsilon; + /** + * Values of forces and nodal displacements which have a magnitude less than + * `epsilon` will be rounded to 0.0. Default = `1e-14`. This is a simple way + * to deal with machine precision when doing calculations. + */ + double epsilon; - /** - * Number of decimal places to use when saving nodal forces and displacements if either - * `save_nodal_displacements` or `save_nodal_forces` is set to `true`. Default = 14. - */ - unsigned int csv_precision; + /** + * Number of decimal places to use when saving nodal forces and displacements + * if either `save_nodal_displacements` or `save_nodal_forces` is set to + * `true`. Default = 14. + */ + unsigned int csv_precision; - /** - * Delimiter to use when saving nodal forces and displacements if either - * `save_nodal_displacements` or `save_nodal_forces` is set to `true`. Default = ",". - */ - std::string csv_delimiter; + /** + * Delimiter to use when saving nodal forces and displacements if either + * `save_nodal_displacements` or `save_nodal_forces` is set to `true`. Default + * = ",". + */ + std::string csv_delimiter; - /** - * Specifies if the nodal displacements should be saved to a file. Default = `false`. - * If `true` the nodal displacements will be saved to the file indicated by `nodal_displacements_filename`. - */ - bool save_nodal_displacements; + /** + * Specifies if the nodal displacements should be saved to a file. Default = + * `false`. If `true` the nodal displacements will be saved to the file + * indicated by `nodal_displacements_filename`. + */ + bool save_nodal_displacements; - /** - * Specifies if the nodal forces should be saved to a file. Default = `false`. - * If `true` the nodal forces will be saved to the file indicated by `nodal_forces_filename`. - */ - bool save_nodal_forces; + /** + * Specifies if the nodal forces should be saved to a file. Default = `false`. + * If `true` the nodal forces will be saved to the file indicated by + * `nodal_forces_filename`. + */ + bool save_nodal_forces; - /** - * Specifies if the forces associated with tie elements should be saved to a file. Default = `false`. - * If `true` the tie forces will be saved to the file indicated by `tie_forces_filename`. - */ - bool save_tie_forces; + /** + * Specifies if the elemental forces should be saved to a file. Default = + * `false`. If `true` the elemental forces will be saved to the file indicated + * by `elemental_forces_filename`. + */ + bool save_elemental_forces; - /** - * Specifies if progress of the analysis should be written to std::cout. Default = `false`. - * If `true` information on current step and time taken per step will be reported. - */ - bool verbose; + /** + * Specifies if the forces associated with tie elements should be saved to a + * file. Default = `false`. If `true` the tie forces will be saved to the file + * indicated by `tie_forces_filename`. + */ + bool save_tie_forces; - /** - * Specifies if a text file should be written detailing information on the analysis. Default = `false`. - * If `true` the a report will be saved to the file indicated by `report_filename`. - */ - bool save_report; + /** + * Specifies if progress of the analysis should be written to std::cout. + * Default = `false`. If `true` information on current step and time taken per + * step will be reported. + */ + bool verbose; - /** - * File name to save the nodal displacements to when `save_nodal_displacements == true`. - */ - std::string nodal_displacements_filename; + /** + * Specifies if a text file should be written detailing information on the + * analysis. Default = `false`. If `true` the a report will be saved to the + * file indicated by `report_filename`. + */ + bool save_report; - /** - * File name to save the nodal forces to when `save_nodal_forces == true`. - */ - std::string nodal_forces_filename; + /** + * File name to save the nodal displacements to when `save_nodal_displacements + * == true`. + */ + std::string nodal_displacements_filename; - /** - * File name to save the nodal forces to when `save_tie_forces == true`. - */ - std::string tie_forces_filename; + /** + * File name to save the nodal forces to when `save_nodal_forces == true`. + */ + std::string nodal_forces_filename; - /** - * File name to save the nodal forces to when `save_report == true`. - */ - std::string report_filename; + /** + * File name to save the nodal forces to when `save_nodal_forces == true`. + */ + std::string elemental_forces_filename; - }; + /** + * File name to save the nodal forces to when `save_tie_forces == true`. + */ + std::string tie_forces_filename; + + /** + * File name to save the nodal forces to when `save_report == true`. + */ + std::string report_filename; +}; } // namespace fea -#endif //THREEDBEAMFEA_OPTIONS_H +#endif // THREEDBEAMFEA_OPTIONS_H diff --git a/include/threed_beam_fea.h b/include/threed_beam_fea.h index ff63662..f71c47a 100644 --- a/include/threed_beam_fea.h +++ b/include/threed_beam_fea.h @@ -75,6 +75,11 @@ typedef Eigen::Matrix ForceVector; */ typedef Eigen::SparseMatrix SparseMat; +/** + * A 3x3 rotation matrix. + */ +using RotationMatrix = Eigen::Matrix; + /** * @brief Calculates the distance between 2 nodes. * @details Calculates the original Euclidean distance between 2 nodes in the @@ -143,7 +148,7 @@ public: * @param[in] ny `Eigen::Matrix3d`. Unit normal vector in global space * parallel to the beam element's local y-direction. */ - void calcAelem(const Eigen::Vector3d &nx, const Eigen::Vector3d &nz); + void calcAelem(const RotationMatrix &r); /** * @brief Returns the currently stored elemental stiffness matrix. diff --git a/src/threed_beam_fea.cpp b/src/threed_beam_fea.cpp index 13bea6d..5d35ef5 100644 --- a/src/threed_beam_fea.cpp +++ b/src/threed_beam_fea.cpp @@ -35,6 +35,8 @@ #include "threed_beam_fea.h" +#define DEBUG + namespace fea { namespace { @@ -126,14 +128,17 @@ void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) { Klocal(11, 11) = 4.0 * tmp1z; // calculate unit normal vector along local x-direction - Eigen::Vector3d nx = job.nodes[nn2] - job.nodes[nn1]; - nx.normalize(); - + const Eigen::Vector3d nx = (job.nodes[nn2] - job.nodes[nn1]).normalized(); // calculate unit normal vector along y-direction const Eigen::Vector3d ny = job.props[i].normal_vec.normalized(); - + // calculate the unit normal vector in local z direction + const Eigen::Vector3d nz = nx.cross(ny).normalized(); + RotationMatrix r; + r.row(0) = nx; + r.row(1) = ny; + r.row(2) = nz; // update rotation matrices - calcAelem(nx, ny); + calcAelem(r); AelemT = Aelem.transpose(); std::ofstream KlocalFile("Klocal.csv"); @@ -159,102 +164,12 @@ void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) { perElemKlocalAelem[i] = KlocalAelem; }; -void GlobalStiffAssembler::calcAelem(const Eigen::Vector3d &nx, - const Eigen::Vector3d &ny) { - // calculate the unit normal vector in local z direction - Eigen::Vector3d nz; - nz = nx.cross(ny).normalized(); - +void GlobalStiffAssembler::calcAelem(const RotationMatrix &r) { // update rotation matrix - Aelem(0, 0) = nx(0); - Aelem(0, 1) = nx(1); - Aelem(0, 2) = nx(2); - Aelem(1, 0) = ny(0); - Aelem(1, 1) = ny(1); - Aelem(1, 2) = ny(2); - Aelem(2, 0) = nz(0); - Aelem(2, 1) = nz(1); - Aelem(2, 2) = nz(2); - - Aelem(3, 3) = nx(0); - Aelem(3, 4) = nx(1); - Aelem(3, 5) = nx(2); - Aelem(4, 3) = ny(0); - Aelem(4, 4) = ny(1); - Aelem(4, 5) = ny(2); - Aelem(5, 3) = nz(0); - Aelem(5, 4) = nz(1); - Aelem(5, 5) = nz(2); - - Aelem(6, 6) = nx(0); - Aelem(6, 7) = nx(1); - Aelem(6, 8) = nx(2); - Aelem(7, 6) = ny(0); - Aelem(7, 7) = ny(1); - Aelem(7, 8) = ny(2); - Aelem(8, 6) = nz(0); - Aelem(8, 7) = nz(1); - Aelem(8, 8) = nz(2); - - Aelem(9, 9) = nx(0); - Aelem(9, 10) = nx(1); - Aelem(9, 11) = nx(2); - Aelem(10, 9) = ny(0); - Aelem(10, 10) = ny(1); - Aelem(10, 11) = ny(2); - Aelem(11, 9) = nz(0); - Aelem(11, 10) = nz(1); - Aelem(11, 11) = nz(2); - - // std::ofstream AelemFile("Aelem.csv"); - // if (AelemFile.is_open()) { - // const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, - // Eigen::DontAlignCols, ", ", - // "\n"); - // AelemFile << Aelem.format(CSVFormat) << '\n'; - // AelemFile.close(); - // } - - // update transposed rotation matrix - AelemT(0, 0) = nx(0); - AelemT(0, 1) = ny(0); - AelemT(0, 2) = nz(0); - AelemT(1, 0) = nx(1); - AelemT(1, 1) = ny(1); - AelemT(1, 2) = nz(1); - AelemT(2, 0) = nx(2); - AelemT(2, 1) = ny(2); - AelemT(2, 2) = nz(2); - - AelemT(3, 3) = nx(0); - AelemT(3, 4) = ny(0); - AelemT(3, 5) = nz(0); - AelemT(4, 3) = nx(1); - AelemT(4, 4) = ny(1); - AelemT(4, 5) = nz(1); - AelemT(5, 3) = nx(2); - AelemT(5, 4) = ny(2); - AelemT(5, 5) = nz(2); - - AelemT(6, 6) = nx(0); - AelemT(6, 7) = ny(0); - AelemT(6, 8) = nz(0); - AelemT(7, 6) = nx(1); - AelemT(7, 7) = ny(1); - AelemT(7, 8) = nz(1); - AelemT(8, 6) = nx(2); - AelemT(8, 7) = ny(2); - AelemT(8, 8) = nz(2); - - AelemT(9, 9) = nx(0); - AelemT(9, 10) = ny(0); - AelemT(9, 11) = nz(0); - AelemT(10, 9) = nx(1); - AelemT(10, 10) = ny(1); - AelemT(10, 11) = nz(1); - AelemT(11, 9) = nx(2); - AelemT(11, 10) = ny(2); - AelemT(11, 11) = nz(2); + Aelem.block<3, 3>(0, 0) = r; + Aelem.block<3, 3>(3, 3) = r; + Aelem.block<3, 3>(6, 6) = r; + Aelem.block<3, 3>(9, 9) = r; } std::vector GlobalStiffAssembler::getPerElemKlocalAelem() const { @@ -475,6 +390,7 @@ Summary solve(const Job &job, const std::vector &BCs, << " ms.\nNow preprocessing factorization..." << std::endl; unsigned int numberOfDoF = DOF::NUM_DOFS * job.nodes.size(); +#ifdef DEBUG Eigen::MatrixXd KgNoBCDense(Kg.block(0, 0, numberOfDoF, numberOfDoF)); std::ofstream kgNoBCFile("KgNoBC.csv"); if (kgNoBCFile.is_open()) { @@ -483,8 +399,8 @@ Summary solve(const Job &job, const std::vector &BCs, kgNoBCFile << KgNoBCDense.format(CSVFormat) << '\n'; kgNoBCFile.close(); } - std::cout << KgNoBCDense << std::endl; - +// std::cout << KgNoBCDense << std::endl; +#endif // load prescribed boundary conditions into stiffness matrix and force vector loadBCs(Kg, force_vec, BCs, job.nodes.size()); @@ -496,14 +412,16 @@ Summary solve(const Job &job, const std::vector &BCs, if (forces.size() > 0) { loadForces(force_vec, forces); } +#ifdef DEBUG Eigen::MatrixXd KgDense(Kg); - std::cout << KgDense << std::endl; + // std::cout << KgDense << std::endl; Eigen::VectorXd forcesVectorDense(force_vec); - std::cout << forcesVectorDense << std::endl; - +// std::cout << forcesVectorDense << std::endl; +#endif // compress global stiffness matrix since all non-zero values have been added. Kg.prune(1.e-14); Kg.makeCompressed(); +#ifdef DEBUG std::ofstream kgFile("Kg.csv"); if (kgFile.is_open()) { const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, @@ -518,9 +436,7 @@ Summary solve(const Job &job, const std::vector &BCs, forcesFile << forcesVectorDense.format(CSVFormat) << '\n'; forcesFile.close(); } - Eigen::FullPivLU lu(Kg); - // assert(lu.isInvertible()); - +#endif // initialize solver based on whether MKL should be used #ifdef EIGEN_USE_MKL_ALL Eigen::PardisoLU solver; @@ -636,6 +552,8 @@ Summary solve(const Job &job, const std::vector &BCs, CSVParser csv; start_time = std::chrono::high_resolution_clock::now(); if (options.save_nodal_displacements) { + std::cout << "Writing to:" + options.nodal_displacements_filename + << std::endl; csv.write(options.nodal_displacements_filename, disp_vec, options.csv_precision, options.csv_delimiter); } @@ -701,17 +619,23 @@ Summary solve(const Job &job, const std::vector &BCs, perElemKlocalAelem[elemIndex] * elemDisps; for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) { summary.element_forces[elemIndex][static_cast(dofIndex)] = - elemForces(dofIndex); + -elemForces(dofIndex); } // meaning of sign = reference to first node: // + = compression for axial, - traction for (int dofIndex = DOF::NUM_DOFS; dofIndex < 2 * DOF::NUM_DOFS; dofIndex++) { summary.element_forces[elemIndex][static_cast(dofIndex)] = - //-elemForces(dofIndex); +elemForces(dofIndex); } } + + if (options.save_elemental_forces) { + std::cout << "Writing to:" + options.elemental_forces_filename << std::endl; + csv.write(options.elemental_forces_filename, summary.element_forces, + options.csv_precision, options.csv_delimiter); + } + return summary; };