From 03f3fe2e523874e18760925c6223fccb9a933f1d Mon Sep 17 00:00:00 2001 From: Iason Date: Tue, 21 Jan 2020 18:22:19 +0100 Subject: [PATCH] Computation of the element forces and addition of those to the result struct. --- include/containers.h | 702 ++++++++++++----------- include/summary.h | 201 +++---- include/threed_beam_fea.h | 372 +++++++------ src/threed_beam_fea.cpp | 1113 +++++++++++++++++++++---------------- 4 files changed, 1307 insertions(+), 1081 deletions(-) diff --git a/include/containers.h b/include/containers.h index bddca50..7f5d68a 100644 --- a/include/containers.h +++ b/include/containers.h @@ -1,11 +1,11 @@ /*! -* \file containers.h -* -* \author Ryan Latture -* \date 8-12-15 -* -* Contains the structs used to organize the job, BCs, and ties for 3D beam FEA. -*/ + * \file containers.h + * + * \author Ryan Latture + * \date 8-12-15 + * + * Contains the structs used to organize the job, BCs, and ties for 3D beam FEA. + */ // Copyright 2015. All rights reserved. // @@ -35,346 +35,392 @@ #ifndef FEA_CONTAINERS_H #define FEA_CONTAINERS_H -#include #include +#include namespace fea { +/** + * @brief Convenience enumerator for specifying the active degree of freedom in + * a constraint. + */ +enum DOF { + /** + * Displacement along the global x-axis. + */ + DISPLACEMENT_X, + + /** + * Displacement along the global y-axis. + */ + DISPLACEMENT_Y, + + /** + * Displacement along the global z-axis. + */ + DISPLACEMENT_Z, + + /** + * Rotation about the global x-axis. + */ + ROTATION_X, + + /** + * Rotation about the global y-axis. + */ + ROTATION_Y, + + /** + * Rotation about the global z-axis. + */ + ROTATION_Z, + /** + * Number of degrees of freedom per node. + */ + NUM_DOFS +}; + +/** + * @brief A node that describes a mesh. Uses Eigen's predefined Vector class for + * added functionality. + * @details See the Eigen documentation on the Vector3d class for more options + * of what can be done with `Nodes`. \n Examples of constucting a `Node` at + * \f$(x, y, z)=(0,1,2)\f$: + * @code + * // specify values on constuction + * fea::Node n1(1.0, 2.0, 3.0); + * + * // construct a Node then insert values + * fea::Node n2; + * n2 << 0.0, 1.0, 2.0; + * @endcode + */ +typedef Eigen::Vector3d Node; + +/** + * @brief A boundary condition to enforce. + * @details Set by specifying the node to constrain, which degree of freedom, + * and the value to hold the node at. + * + * @code + * // define the node number to constrain + * unsigned int nn1 = 0; + * // define the value to hold the nodal DOF at + * double value = 0.0; + * fea::BC bc(nn1, fea::DOF::DISPLACEMENT_X, value); + * @endcode + */ +struct BC { + unsigned int node; /** normal_vec = {0.0, 0.0, 1.0}; + * fea::Props props(EA, EIz, EIy, GJ, normal_vec); + * @endcode + */ +struct Props { + double EA; /**. Vector normal to element + * (`normal_vec.size()==3`). Direction should be parallel to the beam + * element's local y-axis. + */ + Props(double _EA, double _EIz, double _EIy, double _GJ, + const std::vector &_normal_vec) + : EA(_EA), EIz(_EIz), EIy(_EIy), GJ(_GJ) { + normal_vec << _normal_vec[0], _normal_vec[1], _normal_vec[2]; + }; +}; + +/** + * @brief Places linear springs between all degrees of freedom of 2 nodes. + * @details To form a tie specify the 2 nodes that will be linked as well as the + * spring constants for translational and rotational degrees of freedom. All + * translational degrees of freedom will be assigned the same spring constant. + * The same is true for rotational degrees of freedom, although the spring + * constant does not have to be the same as that used for the translational + * DOFs. + * @code + * // create the tie between node2 and node3 + * unsigned int nn1 = 1; // i.e. the second node in the node list + * unsigned int nn2 = 2; // i.e. the third node in the node list + * + * // define the spring constant for x, y, and z translational DOFs + * double lmult = 100.0; + * + * // define the spring constant for x, y, and z rotational DOFs + * double rmult = 100.0; + * + * // form the tie + * fea::Tie tie1(nn1, nn2, lmult, rmult); + * @endcode + */ +struct Tie { + unsigned int + node_number_1; /** terms; /**`. A list of terms that sum to zero.k + */ + Equation(const std::vector &_terms) : terms(_terms){}; +}; - /** - * @brief Constructor - * @param[in] node `unsigned int`. The index of the node. - * @param[in] dof `unsigned int`. Degree of freedom to constrain (See fea::DOF). - * @param[in] value `double`. The prescribed value for the boundary condition. - */ - BC(unsigned int _node, unsigned int _dof, double _value) : node(_node), dof(_dof), value(_value) { }; - }; +/** + * @brief An element of the mesh. Contains the indices of the two `fea::Node`'s + * that form the element as well as the properties of the element given by the + * `fea::Props` struct. + */ +struct Elem { + Eigen::Vector2i + node_numbers; /** nodes; /** elems; /** + props; /**. The node list that defines the mesh. + * @param[in] elems std::vector. The elements that define the mesh. + * An element is defined by the connectivity list and the + * associated properties. + */ + Job(const std::vector &_nodes, const std::vector _elems) + : nodes(_nodes) { + unsigned int num_elems = _elems.size(); + elems.reserve(num_elems); + props.reserve(num_elems); - /** - * @brief The set of properties associated with an element. - * @details The properties must define the extensional stiffness, \f$EA\f$, bending stiffness parallel to the local z-axis \f$EI_{z}\f$, - * bending stiffness parallel to the local y-axis\f$EI_{y}\f$, the torsional stiffness, \f$GJ\f$, and a vector pointing along the beam elements local y-axis. - * - * @code - * double EA = 1000.0; - * double EIz = 100.0; - * double EIy = 100.0; - * double GJ = 200.0; - * std::vector normal_vec = {0.0, 0.0, 1.0}; - * fea::Props props(EA, EIz, EIy, GJ, normal_vec); - * @endcode - */ - struct Props { - double EA;/**. Vector normal to element (`normal_vec.size()==3`). Direction should be parallel to the beam element's local y-axis. - */ - Props(double _EA, double _EIz, double _EIy, double _GJ, const std::vector &_normal_vec) - : EA(_EA), EIz(_EIz), EIy(_EIy), GJ(_GJ) { - normal_vec << _normal_vec[0], _normal_vec[1], _normal_vec[2]; - }; - }; - - /** - * @brief Places linear springs between all degrees of freedom of 2 nodes. - * @details To form a tie specify the 2 nodes that will be linked as well as the spring constants for translational and rotational degrees of freedom. - * All translational degrees of freedom will be assigned the same spring constant. - * The same is true for rotational degrees of freedom, although the spring constant does not have to be the same as that used for the translational DOFs. - * @code - * // create the tie between node2 and node3 - * unsigned int nn1 = 1; // i.e. the second node in the node list - * unsigned int nn2 = 2; // i.e. the third node in the node list - * - * // define the spring constant for x, y, and z translational DOFs - * double lmult = 100.0; - * - * // define the spring constant for x, y, and z rotational DOFs - * double rmult = 100.0; - * - * // form the tie - * fea::Tie tie1(nn1, nn2, lmult, rmult); - * @endcode - */ - struct Tie { - unsigned int node_number_1;/** terms;/**`. A list of terms that sum to zero.k - */ - Equation(const std::vector &_terms) : terms(_terms) {}; - }; - - /** - * @brief An element of the mesh. Contains the indices of the two `fea::Node`'s that form the element as well - * as the properties of the element given by the `fea::Props` struct. - */ - struct Elem { - Eigen::Vector2i node_numbers;/** nodes;/** elems;/** props;/**. The node list that defines the mesh. - * @param[in] elems std::vector. The elements that define the mesh. - * An element is defined by the connectivity list and the associated properties. - */ - Job(const std::vector &_nodes, const std::vector _elems) : nodes(_nodes) { - unsigned int num_elems = _elems.size(); - elems.reserve(num_elems); - props.reserve(num_elems); - - for (unsigned int i = 0; i < num_elems; i++) { - elems.push_back(_elems[i].node_numbers); - props.push_back(_elems[i].props); - } - }; - }; - - /** - * @brief Convenience enumerator for specifying the active degree of freedom in a constraint. - */ - enum DOF { - /** - * Displacement along the global x-axis. - */ - DISPLACEMENT_X, - - /** - * Displacement along the global y-axis. - */ - DISPLACEMENT_Y, - - /** - * Displacement along the global z-axis. - */ - DISPLACEMENT_Z, - - /** - * Rotation about the global x-axis. - */ - ROTATION_X, - - /** - * Rotation about the global y-axis. - */ - ROTATION_Y, - - /** - * Rotation about the global z-axis. - */ - ROTATION_Z, - /** - * Number of degrees of freedom per node. - */ - NUM_DOFS - }; + for (unsigned int i = 0; i < num_elems; i++) { + elems.push_back(_elems[i].node_numbers); + props.push_back(_elems[i].props); + } + }; +}; } // namespace fea diff --git a/include/summary.h b/include/summary.h index 7956880..faa448a 100644 --- a/include/summary.h +++ b/include/summary.h @@ -31,122 +31,129 @@ namespace fea { - /** - * @brief Contains the results of an analysis after calling `fea::solve`. - */ - struct Summary { +/** + * @brief Contains the results of an analysis after calling `fea::solve`. + */ +struct Summary { - /** - * Default constructor - */ - Summary(); + /** + * Default constructor + */ + Summary(); - /** - * @brief Returns a message containing the results of the analysis. - */ - std::string FullReport() const; + /** + * @brief Returns a message containing the results of the analysis. + */ + std::string FullReport() const; - /** - * The total time of the FE analysis. - */ - long long total_time_in_ms; + /** + * The total time of the FE analysis. + */ + long long total_time_in_ms; - /** - * The time it took to assemble the global stiffness matrix. - */ - long long assembly_time_in_ms; + /** + * The time it took to assemble the global stiffness matrix. + */ + long long assembly_time_in_ms; - /** - * The time to reorder the nonzero elements of the the global stiffness matrix, - * such that the factorization step creates less fill-in. - */ - long long preprocessing_time_in_ms; + /** + * The time to reorder the nonzero elements of the the global stiffness + * matrix, such that the factorization step creates less fill-in. + */ + long long preprocessing_time_in_ms; - /** - * The time to compute the factors of the coefficient matrix. - */ - long long factorization_time_in_ms; + /** + * The time to compute the factors of the coefficient matrix. + */ + long long factorization_time_in_ms; - /** - * The time to compute the solution of the linear system. - */ - long long solve_time_in_ms; + /** + * The time to compute the solution of the linear system. + */ + long long solve_time_in_ms; - /** - * The time to compute the nodal forces. - */ - long long nodal_forces_solve_time_in_ms; + /** + * The time to compute the nodal forces. + */ + long long nodal_forces_solve_time_in_ms; - /** - * The time to compute the nodal forces. - */ - long long tie_forces_solve_time_in_ms; + /** + * The time to compute the nodal forces. + */ + long long tie_forces_solve_time_in_ms; - /** - * The time to compute the nodal forces. - * Does not include the time to save the summary itself if `Options::save_report == true` since all data must be - * computed prior to saving the report. - */ - long long file_save_time_in_ms; + /** + * The time to compute the nodal forces. + * Does not include the time to save the summary itself if + * `Options::save_report == true` since all data must be computed prior to + * saving the report. + */ + long long file_save_time_in_ms; - /** - * The number of nodes in the analysis. - */ - unsigned long num_nodes; + /** + * The number of nodes in the analysis. + */ + unsigned long num_nodes; - /** - * The number of elements in the analysis. - */ - unsigned long num_elems; + /** + * The number of elements in the analysis. + */ + unsigned long num_elems; - /** - * The number of boundary conditions in the analysis. - */ - unsigned long num_bcs; + /** + * The number of boundary conditions in the analysis. + */ + unsigned long num_bcs; - /** - * The number of prescribed forces in the analysis. - */ - unsigned long num_forces; + /** + * The number of prescribed forces in the analysis. + */ + unsigned long num_forces; - /** - * The number of tie constraints in the analysis. - */ - unsigned long num_ties; + /** + * The number of tie constraints in the analysis. + */ + unsigned long num_ties; - /** - * The number of equation constraints in the analysis. - */ - unsigned long num_eqns; + /** + * The number of equation constraints in the analysis. + */ + unsigned long num_eqns; - /** - * The resultant nodal displacement from the FE analysis. - * `nodal_displacements` is a 2D vector where each row - * correspond to a node, and the columns correspond to - * `[d_x, d_y, d_z, theta_x, theta_y, theta_z]` - */ - std::vector > nodal_displacements; + /** + * The resultant nodal displacement from the FE analysis. + * `nodal_displacements` is a 2D vector where each row + * correspond to a node, and the columns correspond to + * `[d_x, d_y, d_z, theta_x, theta_y, theta_z]` + */ + std::vector> nodal_displacements; - /** - * The resultant nodal forces from the FE analysis. - * `nodal_displacements` is a 2D vector where each row - * correspond to a node, and the columns correspond to - * `[f_x, f_y, f_z, m_x, m_y, m_z]` - */ - std::vector > nodal_forces; + /** + * The resultant nodal forces from the FE analysis. + * `nodal_displacements` is a 2D vector where each row + * correspond to a node, and the columns correspond to + * `[f_x, f_y, f_z, m_x, m_y, m_z]` + */ + std::vector> nodal_forces; - /** - * The resultant forces associated with ties between nodes. - * `tie_forces` is a 2D vector where each row - * correspond to a tie, and the columns correspond to - * `[f_x, f_y, f_z, f_rot_x, f_rot_y, f_rot_z]` - */ - std::vector > tie_forces; + /** + * The resultant forces associated with ties between nodes. + * `tie_forces` is a 2D vector where each row + * correspond to a tie, and the columns correspond to + * `[f_x, f_y, f_z, f_rot_x, f_rot_y, f_rot_z]` + */ + std::vector> tie_forces; - }; + /** + * The resultant forces associated each element. + * `element_forces` is a 2D vector where each row + * correspond to an element, and the columns correspond to + * `[f1_x, f1_y, f1_z, f1_rot_x, f1_rot_y, f1_rot_z, + * f2_x, f2_y, f2_z, f2_rot_x, f2_rot_y, f2_rot_z]` + */ + std::vector> element_forces; +}; -} //namespace fea +} // namespace fea - - -#endif //THREEDBEAMFEA_SUMMARY_H +#endif // THREEDBEAMFEA_SUMMARY_H diff --git a/include/threed_beam_fea.h b/include/threed_beam_fea.h index bab7aa6..ff63662 100644 --- a/include/threed_beam_fea.h +++ b/include/threed_beam_fea.h @@ -1,11 +1,11 @@ /*! -* \file threed_beam_fea.h -* -* \author Ryan Latture -* \date 8-12-15 -* -* Contains declarations for 3D beam FEA functions. -*/ + * \file threed_beam_fea.h + * + * \author Ryan Latture + * \date 8-12-15 + * + * Contains declarations for 3D beam FEA functions. + */ // Copyright 2015. All rights reserved. // @@ -46,192 +46,228 @@ #include #include "containers.h" +#include "csv_parser.h" #include "options.h" #include "summary.h" -#include "csv_parser.h" namespace fea { - /** - * Dense global stiffness matrix - */ - typedef Eigen::Matrix GlobalStiffMatrix; +/** + * Dense global stiffness matrix + */ +typedef Eigen::Matrix GlobalStiffMatrix; - /** - * An elemental matrix in local coordinates. Will either be the elemental stiffness matrix or the global-to-local rotation matrix - */ - typedef Eigen::Matrix LocalMatrix; +/** + * An elemental matrix in local coordinates. Will either be the elemental + * stiffness matrix or the global-to-local rotation matrix + */ +typedef Eigen::Matrix LocalMatrix; - /** - * Vector that stores the nodal forces, i.e. the variable \f$[F]\f$ in \f$[K][Q]=[F]\f$, - * where \f$[K]\f$ is the global stiffness matrix and \f$[Q]\f$ contains the nodal displacements - */ - typedef Eigen::Matrix ForceVector; +/** + * Vector that stores the nodal forces, i.e. the variable \f$[F]\f$ in + * \f$[K][Q]=[F]\f$, where \f$[K]\f$ is the global stiffness matrix and + * \f$[Q]\f$ contains the nodal displacements + */ +typedef Eigen::Matrix ForceVector; - /** - * Sparse matrix that is used internally to hold the global stiffness matrix - */ - typedef Eigen::SparseMatrix SparseMat; +/** + * Sparse matrix that is used internally to hold the global stiffness matrix + */ +typedef Eigen::SparseMatrix SparseMat; - /** - * @brief Calculates the distance between 2 nodes. - * @details Calculates the original Euclidean distance between 2 nodes in the x-y plane. - * - * @param[in] n1 `fea::Node`. Nodal coordinates of first point. - * @param[in] n2 `fea::Node`. Nodal coordinates of second point. - * - * @return Distance `double`. The distance between the nodes. - */ - inline double norm(const Node &n1, const Node &n2); +/** + * @brief Calculates the distance between 2 nodes. + * @details Calculates the original Euclidean distance between 2 nodes in the + * x-y plane. + * + * @param[in] n1 `fea::Node`. Nodal coordinates of first point. + * @param[in] n2 `fea::Node`. Nodal coordinates of second point. + * + * @return Distance `double`. The distance between the nodes. + */ +inline double norm(const Node &n1, const Node &n2); - /** - * @brief Assembles the global stiffness matrix. - */ - class GlobalStiffAssembler { +/** + * @brief Assembles the global stiffness matrix. + */ +class GlobalStiffAssembler { - public: +public: + /** + * @brief Default constructor. + * @details Initializes all entries in member matrices to 0.0. + */ + GlobalStiffAssembler() { + Kelem.setZero(); + Klocal.setZero(); + Aelem.setZero(); + AelemT.setZero(); + SparseKelem.resize(12, 12); + SparseKelem.reserve(40); + }; - /** - * @brief Default constructor. - * @details Initializes all entries in member matrices to 0.0. - */ - GlobalStiffAssembler() { - Kelem.setZero(); - Klocal.setZero(); - Aelem.setZero(); - AelemT.setZero(); - SparseKelem.resize(12, 12); - SparseKelem.reserve(40); - }; + /** + * @brief Assembles the global stiffness matrix. + * @details The input stiffness matrix is modified in place to contain the + * correct values for the given job. Assumes that the input stiffness matrix + * has the correct dimensions and all values are initially set to zero. + * + * @param Kg `fea::GlobalStiffMatrix`. Modified in place. After evaluation, Kg + * contains the correct values for the global stiffness matrix due to the + * input Job. + * @param[in] job `fea::Job`. Current Job to analyze contains node, element, + * and property lists. + * @param[in] ties `std::vector`. Vector of ties that apply to + * attach springs of specified stiffness to all nodal degrees of freedom + * between each set of nodes indicated. + */ + void operator()(SparseMat &Kg, const Job &job, const std::vector &ties); - /** - * @brief Assembles the global stiffness matrix. - * @details The input stiffness matrix is modified in place to contain the correct values for the given job. - * Assumes that the input stiffness matrix has the correct dimensions and all values are initially set to zero. - * - * @param Kg `fea::GlobalStiffMatrix`. Modified in place. After evaluation, Kg contains the correct values for - * the global stiffness matrix due to the input Job. - * @param[in] job `fea::Job`. Current Job to analyze contains node, element, and property lists. - * @param[in] ties `std::vector`. Vector of ties that apply to attach springs of specified stiffness to - * all nodal degrees of freedom between each set of nodes indicated. - */ - void operator()(SparseMat &Kg, const Job &job, const std::vector &ties); + /** + * @brief Updates the elemental stiffness matrix for the `ith` element. + * + * @param[in] i `unsigned int`. Specifies the ith element for which the + * elemental stiffness matrix is calculated. + * @param[in] job `Job`. Current `fea::Job` to analyze contains node, element, + * and property lists. + */ + void calcKelem(unsigned int i, const Job &job); - /** - * @brief Updates the elemental stiffness matrix for the `ith` element. - * - * @param[in] i `unsigned int`. Specifies the ith element for which the elemental stiffness matrix is calculated. - * @param[in] job `Job`. Current `fea::Job` to analyze contains node, element, and property lists. - */ - void calcKelem(unsigned int i, const Job &job); + /** + * @brief Updates the rotation and transposed rotation matrices. + * @details The rotation matrices `Aelem` and `AelemT` are updated based on + * the 2 specified unit normal vectors along the local x and y directions. + * + * @param[in] nx `Eigen::Matrix3d`. Unit normal vector in global space + * parallel to the beam element's local x-direction. + * @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); - /** - * @brief Updates the rotation and transposed rotation matrices. - * @details The rotation matrices `Aelem` and `AelemT` are updated based on the 2 specified unit - * normal vectors along the local x and y directions. - * - * @param[in] nx `Eigen::Matrix3d`. Unit normal vector in global space parallel to the beam element's local x-direction. - * @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); + /** + * @brief Returns the currently stored elemental stiffness matrix. + * @return Elemental stiffness matrix `fea::LocalMatrix`. + */ + LocalMatrix getKelem() { return Kelem; } - /** - * @brief Returns the currently stored elemental stiffness matrix. - * @return Elemental stiffness matrix `fea::LocalMatrix`. - */ - LocalMatrix getKelem() { - return Kelem; - } + /** + * @brief Returns the currently stored rotation matrix. + * @return Rotation matrix `fea::LocalMatrix`. + */ + LocalMatrix getAelem() { return Aelem; } - /** - * @brief Returns the currently stored rotation matrix. - * @return Rotation matrix `fea::LocalMatrix`. - */ - LocalMatrix getAelem() { - return Aelem; - } + std::vector getPerElemKlocalAelem() const; - private: - LocalMatrix Kelem; - /**`. Vector of `BC`'s to apply to the current analysis. - * @param[in] num_nodes `unsigned int`. The number of nodes in the current job being analyzed. - * Used to calculate the position to insert border coefficients associated - * with enforcing boundary conditions via Langrange multipliers. - */ - void loadBCs(SparseMat &Kg, ForceVector &force_vec, const std::vector &BCs, unsigned int num_nodes); + std::vector + perElemKlocalAelem; //[i] holds the local stiffness matrix of the ith beam + // element times the transformation matrix from local + // to global. This is used after the simulation in + // order to find the per element forces +}; - void loadEquations(SparseMat &Kg, const std::vector &equations, unsigned int num_nodes, unsigned int num_bcs); +/** + * @brief Loads the boundary conditions into the global stiffness matrix and + * force vector. + * @details Boundary conditions are enforced via Lagrange multipliers. The + * reaction force due to imposing the boundary condition will be appended + * directly onto the returned nodal displacements in the order the boundary + * conditions were specified. + * + * @param Kg `fea::GlobalStiffnessMatrix`. Coefficients are modified in place to + * reflect Lagrange multipliers. Assumes `Kg` has the correct dimensions. + * @param force_vec `fea::ForceVector`. Right hand side of the \f$[K][Q]=[F]\f$ + * equation of the FE analysis. + * @param[in] BCs `std::vector`. Vector of `BC`'s to apply to the + * current analysis. + * @param[in] num_nodes `unsigned int`. The number of nodes in the current job + * being analyzed. Used to calculate the position to insert border coefficients + * associated with enforcing boundary conditions via Langrange multipliers. + */ +void loadBCs(SparseMat &Kg, ForceVector &force_vec, const std::vector &BCs, + unsigned int num_nodes); - /** - * @brief Loads any tie constraints into the set of triplets that will become the global stiffness matrix. - * @details Tie constraints are enforced via linear springs between the 2 specified - * nodes. The `lmult` member variable is used as the spring constant for displacement degrees - * of freedom, e.g. 0, 1, and 2. `rmult` is used for rotational degrees of freedom, e.g. 3, 4, and 5. - * - * @param triplets `std::vector< Eigen::Triplet< double > >`. A vector of triplets that store data in the - * form (i, j, value) that will be become the sparse global stiffness matrix. - * @param[in] ties `std::vector`. Vector of `Tie`'s to apply to the current analysis. - */ - void loadTies(std::vector > &triplets, const std::vector &ties); +void loadEquations(SparseMat &Kg, const std::vector &equations, + unsigned int num_nodes, unsigned int num_bcs); +/** + * @brief Loads any tie constraints into the set of triplets that will become + * the global stiffness matrix. + * @details Tie constraints are enforced via linear springs between the 2 + * specified nodes. The `lmult` member variable is used as the spring constant + * for displacement degrees of freedom, e.g. 0, 1, and 2. `rmult` is used for + * rotational degrees of freedom, e.g. 3, 4, and 5. + * + * @param triplets `std::vector< Eigen::Triplet< double > >`. A vector of + * triplets that store data in the form (i, j, value) that will be become the + * sparse global stiffness matrix. + * @param[in] ties `std::vector`. Vector of `Tie`'s to apply to the + * current analysis. + */ +void loadTies(std::vector> &triplets, + const std::vector &ties); - /** - * @brief Computes the forces in the tie elements based on the nodal displacements of the FE - * analysis and the spring constants provided in `ties`. - * - * @param[in] ties `std::vector`. Vector of `fea::Tie`'s to applied to the current analysis. - * @param[in] nodal_displacements `std::vector < std::vector < double > >`. The resultant nodal displacements of the analysis. - * @return Tie forces. `std::vector < std::vector < double > >` - */ - std::vector > computeTieForces(const std::vector &ties, - const std::vector > &nodal_displacements); +/** + * @brief Computes the forces in the tie elements based on the nodal + * displacements of the FE analysis and the spring constants provided in `ties`. + * + * @param[in] ties `std::vector`. Vector of `fea::Tie`'s to applied to the + * current analysis. + * @param[in] nodal_displacements `std::vector < std::vector < double > >`. The + * resultant nodal displacements of the analysis. + * @return Tie forces. `std::vector < std::vector < double > >` + */ +std::vector> +computeTieForces(const std::vector &ties, + const std::vector> &nodal_displacements); - /** - * @brief Loads the prescribed forces into the force vector. - * - * @param force_vec `ForceVector`. Right hand side of the \f$[K][Q]=[F]\f$ equation of the FE analysis. - * @param[in] forces std::vector. Vector of prescribed forces to apply to the current analysis. - */ - void loadForces(SparseMat &force_vec, const std::vector &forces); +/** + * @brief Loads the prescribed forces into the force vector. + * + * @param force_vec `ForceVector`. Right hand side of the \f$[K][Q]=[F]\f$ + * equation of the FE analysis. + * @param[in] forces std::vector. Vector of prescribed forces to apply to + * the current analysis. + */ +void loadForces(SparseMat &force_vec, const std::vector &forces); - /** - * @brief Solves the finite element analysis defined by the input Job, boundary conditions, and prescribed nodal forces. - * @details Solves \f$[K][Q]=[F]\f$ for \f$[Q]\f$, where \f$[K]\f$ is the global stiffness matrix, - * \f$[Q]\f$ contains the nodal displacements, and \f$[Q]\f$ contains the nodal forces. - * - * @param[in] job `fea::Job`. Contains the node, element, and property lists for the mesh. - * @param[in] BCs `std::vector`. Vector of boundary conditions to apply to the nodal degrees of freedom contained in the job. - * @param[in] forces `std::vector`. Vector of prescribed forces to apply to the nodal degrees of freedom contained in the job. - * @param[in] ties `std::vector`. Vector of ties that apply to attach springs of specified stiffness to - * all nodal degrees of freedom between each set of nodes indicated. - * - * @return Summary `fea::Summary`. Summary containing the results of the analysis. - */ - Summary solve(const Job &job, - const std::vector &BCs, - const std::vector &forces, - const std::vector &ties, - const std::vector &equations, - const Options &options); +/** + * @brief Solves the finite element analysis defined by the input Job, boundary + * conditions, and prescribed nodal forces. + * @details Solves \f$[K][Q]=[F]\f$ for \f$[Q]\f$, where \f$[K]\f$ is the global + * stiffness matrix, \f$[Q]\f$ contains the nodal displacements, and \f$[Q]\f$ + * contains the nodal forces. + * + * @param[in] job `fea::Job`. Contains the node, element, and property lists for + * the mesh. + * @param[in] BCs `std::vector`. Vector of boundary conditions to apply + * to the nodal degrees of freedom contained in the job. + * @param[in] forces `std::vector`. Vector of prescribed forces to + * apply to the nodal degrees of freedom contained in the job. + * @param[in] ties `std::vector`. Vector of ties that apply to attach + * springs of specified stiffness to all nodal degrees of freedom between each + * set of nodes indicated. + * + * @return Summary `fea::Summary`. Summary containing the results of the + * analysis. + */ +Summary solve(const Job &job, const std::vector &BCs, + const std::vector &forces, const std::vector &ties, + const std::vector &equations, const Options &options); } // namespace fea #endif // THREED_BEAM_FEA_H diff --git a/src/threed_beam_fea.cpp b/src/threed_beam_fea.cpp index a080683..13bea6d 100644 --- a/src/threed_beam_fea.cpp +++ b/src/threed_beam_fea.cpp @@ -23,6 +23,7 @@ // // Author: ryan.latture@gmail.com (Ryan Latture) +#include #include #include #include @@ -36,546 +37,682 @@ namespace fea { - namespace { - void writeStringToTxt(std::string filename, std::string data) { - std::ofstream output_file; - output_file.open(filename); +namespace { +void writeStringToTxt(std::string filename, std::string data) { + std::ofstream output_file; + output_file.open(filename); - if (!output_file.is_open()) { - throw std::runtime_error( - (boost::format("Error opening file %s.") % filename).str() - ); - } - output_file << data; - output_file.close(); + if (!output_file.is_open()) { + throw std::runtime_error( + (boost::format("Error opening file %s.") % filename).str()); + } + output_file << data; + output_file.close(); +} +} // namespace + +inline double norm(const Node &n1, const Node &n2) { + const Node dn = n2 - n1; + return dn.norm(); +} + +void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) { + // extract element properties + const double EA = job.props[i].EA; // Young's modulus * cross area + const double EIz = job.props[i].EIz; // Young's modulus* I3 + const double EIy = job.props[i].EIy; // Young's modulus* I2 + const double GJ = job.props[i].GJ; + + // store node indices of current element + const int nn1 = job.elems[i][0]; + const int nn2 = job.elems[i][1]; + + // calculate the length of the element + const double length = norm(job.nodes[nn1], job.nodes[nn2]); + + // store the entries in the (local) elemental stiffness matrix as temporary + // values to avoid recalculation + const double tmpEA = EA / length; + const double tmpGJ = GJ / length; + + const double tmp12z = 12.0 * EIz / (length * length * length); // ==k3 + const double tmp6z = 6.0 * EIz / (length * length); + const double tmp1z = EIz / length; + + const double tmp12y = 12.0 * EIy / (length * length * length); + const double tmp6y = 6.0 * EIy / (length * length); + const double tmp1y = EIy / length; + + // update local elemental stiffness matrix + Klocal(0, 0) = tmpEA; + Klocal(0, 6) = -tmpEA; + Klocal(1, 1) = tmp12z; + Klocal(1, 5) = tmp6z; + Klocal(1, 7) = -tmp12z; + Klocal(1, 11) = tmp6z; + Klocal(2, 2) = tmp12y; + Klocal(2, 4) = -tmp6y; + Klocal(2, 8) = -tmp12y; + Klocal(2, 10) = -tmp6y; + Klocal(3, 3) = tmpGJ; + Klocal(3, 9) = -tmpGJ; + Klocal(4, 2) = -tmp6y; + Klocal(4, 4) = 4.0 * tmp1y; + Klocal(4, 8) = tmp6y; + Klocal(4, 10) = 2.0 * tmp1y; + Klocal(5, 1) = tmp6z; + Klocal(5, 5) = 4.0 * tmp1z; + Klocal(5, 7) = -tmp6z; + Klocal(5, 11) = 2.0 * tmp1z; + Klocal(6, 0) = -tmpEA; + Klocal(6, 6) = tmpEA; + Klocal(7, 1) = -tmp12z; + Klocal(7, 5) = -tmp6z; + Klocal(7, 7) = tmp12z; + Klocal(7, 11) = -tmp6z; + Klocal(8, 2) = -tmp12y; + Klocal(8, 4) = tmp6y; + Klocal(8, 8) = tmp12y; + Klocal(8, 10) = tmp6y; + Klocal(9, 3) = -tmpGJ; + Klocal(9, 9) = tmpGJ; + Klocal(10, 2) = -tmp6y; + Klocal(10, 4) = 2.0 * tmp1y; + Klocal(10, 8) = tmp6y; + Klocal(10, 10) = 4.0 * tmp1y; + Klocal(11, 1) = tmp6z; + Klocal(11, 5) = 2.0 * tmp1z; + Klocal(11, 7) = -tmp6z; + 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(); + + // calculate unit normal vector along y-direction + const Eigen::Vector3d ny = job.props[i].normal_vec.normalized(); + + // update rotation matrices + calcAelem(nx, ny); + AelemT = Aelem.transpose(); + + std::ofstream KlocalFile("Klocal.csv"); + if (KlocalFile.is_open()) { + const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, + Eigen::DontAlignCols, ", ", "\n"); + KlocalFile << Klocal.format(CSVFormat) << '\n'; + KlocalFile.close(); + } + // update Kelem + Kelem = AelemT * Klocal * Aelem; + // std::ofstream KelemFile("Kelem.csv"); + // if (KelemFile.is_open()) { + // const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, + // Eigen::DontAlignCols, ", ", + // "\n"); + // KelemFile << Kelem.format(CSVFormat) << '\n'; + // KelemFile.close(); + // } + + LocalMatrix KlocalAelem; + KlocalAelem = Klocal * Aelem; + 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(); + + // 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); +} + +std::vector GlobalStiffAssembler::getPerElemKlocalAelem() const { + return perElemKlocalAelem; +}; + +void GlobalStiffAssembler::operator()(SparseMat &Kg, const Job &job, + const std::vector &ties) { + int nn1, nn2; + unsigned int row, col; + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + + // form vector to hold triplets that will be used to assemble global stiffness + // matrix + std::vector> triplets; + triplets.reserve(40 * job.elems.size() + 4 * dofs_per_elem * ties.size()); + + perElemKlocalAelem.resize(job.elems.size()); + for (unsigned int i = 0; i < job.elems.size(); ++i) { + // update Kelem with current elemental stiffness matrix + calcKelem(i, job); // 12x12 matrix + + // get sparse representation of the current elemental stiffness matrix + SparseKelem = Kelem.sparseView(); + + // pull out current node indices + nn1 = job.elems[i][0]; + nn2 = job.elems[i][1]; + + for (unsigned int j = 0; j < SparseKelem.outerSize(); ++j) { + for (SparseMat::InnerIterator it(SparseKelem, j); it; ++it) { + row = it.row(); + col = it.col(); + + // check position in local matrix and update corresponding global + // position + if (row < 6) { + // top left + if (col < 6) { + triplets.push_back(Eigen::Triplet(dofs_per_elem * nn1 + row, + dofs_per_elem * nn1 + col, + it.value())); + } + // top right + else { + triplets.push_back(Eigen::Triplet( + dofs_per_elem * nn1 + row, dofs_per_elem * (nn2 - 1) + col, + it.value())); // I: Giati nn2-1 kai oxi nn2? + } + } else { + // bottom left + if (col < 6) { + triplets.push_back(Eigen::Triplet( + dofs_per_elem * (nn2 - 1) + row, dofs_per_elem * nn1 + col, + it.value())); // I: Giati nn2-1 kai oxi nn2? Epeidi ta nn2 ston + // Kelem exoun idi offset 6 + } + // bottom right + else { + triplets.push_back(Eigen::Triplet( + dofs_per_elem * (nn2 - 1) + row, + dofs_per_elem * (nn2 - 1) + col, it.value())); + } } + } } + } - inline double norm(const Node &n1, const Node &n2) { - const Node dn = n2 - n1; - return dn.norm(); + loadTies(triplets, ties); + + Kg.setFromTriplets(triplets.begin(), triplets.end()); +}; + +void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector &BCs, + unsigned int num_nodes) { + unsigned int bc_idx; + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + // calculate the index that marks beginning of Lagrange multiplier + // coefficients + const unsigned int global_add_idx = dofs_per_elem * num_nodes; + + for (size_t i = 0; i < BCs.size(); ++i) { + bc_idx = dofs_per_elem * BCs[i].node + BCs[i].dof; + + // update global stiffness matrix + Kg.insert(bc_idx, global_add_idx + i) = 1; + Kg.insert(global_add_idx + i, bc_idx) = 1; + + // update force vector. All values are already zero. Only update if BC if + // non-zero. + if (std::abs(BCs[i].value) > std::numeric_limits::epsilon()) { + force_vec.insert(global_add_idx + i, 0) = BCs[i].value; } + } +}; - void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) { - // extract element properties - const double EA = job.props[i].EA; - const double EIz = job.props[i].EIz; - const double EIy = job.props[i].EIy; - const double GJ = job.props[i].GJ; +void loadEquations(SparseMat &Kg, const std::vector &equations, + unsigned int num_nodes, unsigned int num_bcs) { + size_t row_idx, col_idx; + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + const unsigned int global_add_idx = dofs_per_elem * num_nodes + num_bcs; - // store node indices of current element - const int nn1 = job.elems[i][0]; - const int nn2 = job.elems[i][1]; - - // calculate the length of the element - const double length = norm(job.nodes[nn1], job.nodes[nn2]); - - // store the entries in the (local) elemental stiffness matrix as temporary values to avoid recalculation - const double tmpEA = EA / length; - const double tmpGJ = GJ / length; - - const double tmp12z = 12.0 * EIz / (length * length * length); - const double tmp6z = 6.0 * EIz / (length * length); - const double tmp1z = EIz / length; - - const double tmp12y = 12.0 * EIy / (length * length * length); - const double tmp6y = 6.0 * EIy / (length * length); - const double tmp1y = EIy / length; - - // update local elemental stiffness matrix - Klocal(0, 0) = tmpEA; - Klocal(0, 6) = -tmpEA; - Klocal(1, 1) = tmp12z; - Klocal(1, 5) = tmp6z; - Klocal(1, 7) = -tmp12z; - Klocal(1, 11) = tmp6z; - Klocal(2, 2) = tmp12y; - Klocal(2, 4) = -tmp6y; - Klocal(2, 8) = -tmp12y; - Klocal(2, 10) = -tmp6y; - Klocal(3, 3) = tmpGJ; - Klocal(3, 9) = -tmpGJ; - Klocal(4, 2) = -tmp6y; - Klocal(4, 4) = 4.0 * tmp1y; - Klocal(4, 8) = tmp6y; - Klocal(4, 10) = 2.0 * tmp1y; - Klocal(5, 1) = tmp6z; - Klocal(5, 5) = 4.0 * tmp1z; - Klocal(5, 7) = -tmp6z; - Klocal(5, 11) = 2.0 * tmp1z; - Klocal(6, 0) = -tmpEA; - Klocal(6, 6) = tmpEA; - Klocal(7, 1) = -tmp12z; - Klocal(7, 5) = -tmp6z; - Klocal(7, 7) = tmp12z; - Klocal(7, 11) = -tmp6z; - Klocal(8, 2) = -tmp12y; - Klocal(8, 4) = tmp6y; - Klocal(8, 8) = tmp12y; - Klocal(8, 10) = tmp6y; - Klocal(9, 3) = -tmpGJ; - Klocal(9, 9) = tmpGJ; - Klocal(10, 2) = -tmp6y; - Klocal(10, 4) = 2.0 * tmp1y; - Klocal(10, 8) = tmp6y; - Klocal(10, 10) = 4.0 * tmp1y; - Klocal(11, 1) = tmp6z; - Klocal(11, 5) = 2.0 * tmp1z; - Klocal(11, 7) = -tmp6z; - 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(); - - // calculate unit normal vector along y-direction - const Eigen::Vector3d ny = job.props[i].normal_vec.normalized(); - - // update rotation matrices - calcAelem(nx, ny); - - // update Kelem - Kelem = AelemT * Klocal * Aelem; - }; - - 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); - const double dlz = nz.squaredNorm(); - nz /= dlz; - - // 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); - - // 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); - }; - - void GlobalStiffAssembler::operator()(SparseMat &Kg, const Job &job, const std::vector &ties) { - int nn1, nn2; - unsigned int row, col; - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - - // form vector to hold triplets that will be used to assemble global stiffness matrix - std::vector > triplets; - triplets.reserve(40 * job.elems.size() + 4 * dofs_per_elem * ties.size()); - - for (unsigned int i = 0; i < job.elems.size(); ++i) { - // update Kelem with current elemental stiffness matrix - calcKelem(i, job); - - // get sparse representation of the current elemental stiffness matrix - SparseKelem = Kelem.sparseView(); - - // pull out current node indices - nn1 = job.elems[i][0]; - nn2 = job.elems[i][1]; - - for (unsigned int j = 0; j < SparseKelem.outerSize(); ++j) { - for (SparseMat::InnerIterator it(SparseKelem, j); it; ++it) { - row = it.row(); - col = it.col(); - - // check position in local matrix and update corresponding global position - if (row < 6) { - // top left - if (col < 6) { - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn1 + row, - dofs_per_elem * nn1 + col, - it.value())); - } - // top right - else { - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn1 + row, - dofs_per_elem * (nn2 - 1) + col, - it.value())); - } - } - else { - // bottom left - if (col < 6) { - triplets.push_back(Eigen::Triplet(dofs_per_elem * (nn2 - 1) + row, - dofs_per_elem * nn1 + col, - it.value())); - } - // bottom right - else { - triplets.push_back(Eigen::Triplet(dofs_per_elem * (nn2 - 1) + row, - dofs_per_elem * (nn2 - 1) + col, - it.value())); - } - } - } - } - } - - loadTies(triplets, ties); - - Kg.setFromTriplets(triplets.begin(), triplets.end()); - }; - - void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector &BCs, unsigned int num_nodes) { - unsigned int bc_idx; - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - // calculate the index that marks beginning of Lagrange multiplier coefficients - const unsigned int global_add_idx = dofs_per_elem * num_nodes; - - for (size_t i = 0; i < BCs.size(); ++i) { - bc_idx = dofs_per_elem * BCs[i].node + BCs[i].dof; - - // update global stiffness matrix - Kg.insert(bc_idx, global_add_idx + i) = 1; - Kg.insert(global_add_idx + i, bc_idx) = 1; - - // update force vector. All values are already zero. Only update if BC if non-zero. - if (std::abs(BCs[i].value) > std::numeric_limits::epsilon()) { - force_vec.insert(global_add_idx + i, 0) = BCs[i].value; - } - } - }; - - void loadEquations(SparseMat &Kg, const std::vector &equations, unsigned int num_nodes, unsigned int num_bcs) { - size_t row_idx, col_idx; - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - const unsigned int global_add_idx = dofs_per_elem * num_nodes + num_bcs; - - for (size_t i = 0; i < equations.size(); ++i) { - row_idx = global_add_idx + i; - for (size_t j = 0; j < equations[i].terms.size(); ++j) { - col_idx = dofs_per_elem * equations[i].terms[j].node_number + equations[i].terms[j].dof; - Kg.insert(row_idx, col_idx) = equations[i].terms[j].coefficient; - Kg.insert(col_idx, row_idx) = equations[i].terms[j].coefficient; - } - } - }; - - void loadTies(std::vector > &triplets, const std::vector &ties) { - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - unsigned int nn1, nn2; - double lmult, rmult, spring_constant; - - for (size_t i = 0; i < ties.size(); ++i) { - nn1 = ties[i].node_number_1; - nn2 = ties[i].node_number_2; - lmult = ties[i].lmult; - rmult = ties[i].rmult; - - for (unsigned int j = 0; j < dofs_per_elem; ++j) { - // first 3 DOFs are linear DOFs, second 2 are rotational, last is torsional - spring_constant = j < 3 ? lmult : rmult; - - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn1 + j, - dofs_per_elem * nn1 + j, - spring_constant)); - - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn2 + j, - dofs_per_elem * nn2 + j, - spring_constant)); - - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn1 + j, - dofs_per_elem * nn2 + j, - -spring_constant)); - - triplets.push_back(Eigen::Triplet(dofs_per_elem * nn2 + j, - dofs_per_elem * nn1 + j, - -spring_constant)); - } - } - }; - - std::vector > computeTieForces(const std::vector &ties, - const std::vector > &nodal_displacements) { - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - unsigned int nn1, nn2; - double lmult, rmult, spring_constant, delta1, delta2; - - std::vector > tie_forces(ties.size(), std::vector(dofs_per_elem)); - - for (size_t i = 0; i < ties.size(); ++i) { - nn1 = ties[i].node_number_1; - nn2 = ties[i].node_number_2; - lmult = ties[i].lmult; - rmult = ties[i].rmult; - - for (unsigned int j = 0; j < dofs_per_elem; ++j) { - // first 3 DOFs are linear DOFs, second 2 are rotational, last is torsional - spring_constant = j < 3 ? lmult : rmult; - delta1 = nodal_displacements[nn1][j]; - delta2 = nodal_displacements[nn2][j]; - tie_forces[i][j] = spring_constant * (delta2 - delta1); - } - } - return tie_forces; + for (size_t i = 0; i < equations.size(); ++i) { + row_idx = global_add_idx + i; + for (size_t j = 0; j < equations[i].terms.size(); ++j) { + col_idx = dofs_per_elem * equations[i].terms[j].node_number + + equations[i].terms[j].dof; + Kg.insert(row_idx, col_idx) = equations[i].terms[j].coefficient; + Kg.insert(col_idx, row_idx) = equations[i].terms[j].coefficient; } + } +}; - void loadForces(SparseMat &force_vec, const std::vector &forces) { - const unsigned int dofs_per_elem = DOF::NUM_DOFS; - unsigned int idx; +void loadTies(std::vector> &triplets, + const std::vector &ties) { + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + unsigned int nn1, nn2; + double lmult, rmult, spring_constant; - for (size_t i = 0; i < forces.size(); ++i) { - idx = dofs_per_elem * forces[i].node + forces[i].dof; - force_vec.insert(idx, 0) = forces[i].value; - } - }; + for (size_t i = 0; i < ties.size(); ++i) { + nn1 = ties[i].node_number_1; + nn2 = ties[i].node_number_2; + lmult = ties[i].lmult; + rmult = ties[i].rmult; - Summary solve(const Job &job, - const std::vector &BCs, - const std::vector &forces, - const std::vector &ties, - const std::vector &equations, - const Options &options) { - auto initial_start_time = std::chrono::high_resolution_clock::now(); + for (unsigned int j = 0; j < dofs_per_elem; ++j) { + // first 3 DOFs are linear DOFs, second 2 are rotational, last is + // torsional + spring_constant = j < 3 ? lmult : rmult; - Summary summary; - summary.num_nodes = job.nodes.size(); - summary.num_elems = job.elems.size(); - summary.num_bcs = BCs.size(); - summary.num_ties = ties.size(); + triplets.push_back(Eigen::Triplet( + dofs_per_elem * nn1 + j, dofs_per_elem * nn1 + j, spring_constant)); - const unsigned int dofs_per_elem = DOF::NUM_DOFS; + triplets.push_back(Eigen::Triplet( + dofs_per_elem * nn2 + j, dofs_per_elem * nn2 + j, spring_constant)); - // calculate size of global stiffness matrix and force vector - const unsigned long size = dofs_per_elem * job.nodes.size() + BCs.size() + equations.size(); + triplets.push_back(Eigen::Triplet( + dofs_per_elem * nn1 + j, dofs_per_elem * nn2 + j, -spring_constant)); - // construct global stiffness matrix and force vector - SparseMat Kg(size, size); - SparseMat force_vec(size, 1); - force_vec.reserve(forces.size() + BCs.size()); + triplets.push_back(Eigen::Triplet( + dofs_per_elem * nn2 + j, dofs_per_elem * nn1 + j, -spring_constant)); + } + } +}; - // construct global assembler object and assemble global stiffness matrix - auto start_time = std::chrono::high_resolution_clock::now(); - GlobalStiffAssembler assembleK3D = GlobalStiffAssembler(); - assembleK3D(Kg, job, ties); - auto end_time = std::chrono::high_resolution_clock::now(); - auto delta_time = std::chrono::duration_cast(end_time - start_time).count(); - summary.assembly_time_in_ms = delta_time; +std::vector> +computeTieForces(const std::vector &ties, + const std::vector> &nodal_displacements) { + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + unsigned int nn1, nn2; + double lmult, rmult, spring_constant, delta1, delta2; - if (options.verbose) - std::cout << "Global stiffness matrix assembled in " - << delta_time - << " ms.\nNow preprocessing factorization..." << std::endl; + std::vector> tie_forces( + ties.size(), std::vector(dofs_per_elem)); - // load prescribed boundary conditions into stiffness matrix and force vector - loadBCs(Kg, force_vec, BCs, job.nodes.size()); + for (size_t i = 0; i < ties.size(); ++i) { + nn1 = ties[i].node_number_1; + nn2 = ties[i].node_number_2; + lmult = ties[i].lmult; + rmult = ties[i].rmult; - if (equations.size() > 0) { - loadEquations(Kg, equations, job.nodes.size(), BCs.size()); - } + for (unsigned int j = 0; j < dofs_per_elem; ++j) { + // first 3 DOFs are linear DOFs, second 2 are rotational, last is + // torsional + spring_constant = j < 3 ? lmult : rmult; + delta1 = nodal_displacements[nn1][j]; + delta2 = nodal_displacements[nn2][j]; + tie_forces[i][j] = spring_constant * (delta2 - delta1); + } + } + return tie_forces; +} - // load prescribed forces into force vector - if (forces.size() > 0) { - loadForces(force_vec, forces); - } +void loadForces(SparseMat &force_vec, const std::vector &forces) { + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + unsigned int idx; - // compress global stiffness matrix since all non-zero values have been added. - Kg.prune(1.e-14); - Kg.makeCompressed(); + for (size_t i = 0; i < forces.size(); ++i) { + idx = dofs_per_elem * forces[i].node + forces[i].dof; + force_vec.insert(idx, 0) = forces[i].value; + } +}; - // initialize solver based on whether MKL should be used +Summary solve(const Job &job, const std::vector &BCs, + const std::vector &forces, const std::vector &ties, + const std::vector &equations, const Options &options) { + auto initial_start_time = std::chrono::high_resolution_clock::now(); + + Summary summary; + summary.num_nodes = job.nodes.size(); + summary.num_elems = job.elems.size(); + summary.num_bcs = BCs.size(); + summary.num_ties = ties.size(); + + const unsigned int dofs_per_elem = DOF::NUM_DOFS; + + // calculate size of global stiffness matrix and force vector + const unsigned long size = + dofs_per_elem * job.nodes.size() + BCs.size() + equations.size(); + + // construct global stiffness matrix and force vector + SparseMat Kg(size, size); + SparseMat force_vec(size, 1); + force_vec.reserve(forces.size() + BCs.size()); + + // construct global assembler object and assemble global stiffness matrix + auto start_time = std::chrono::high_resolution_clock::now(); + GlobalStiffAssembler assembleK3D = GlobalStiffAssembler(); + assembleK3D(Kg, job, ties); + auto end_time = std::chrono::high_resolution_clock::now(); + auto delta_time = std::chrono::duration_cast( + end_time - start_time) + .count(); + summary.assembly_time_in_ms = delta_time; + + if (options.verbose) + std::cout << "Global stiffness matrix assembled in " << delta_time + << " ms.\nNow preprocessing factorization..." << std::endl; + + unsigned int numberOfDoF = DOF::NUM_DOFS * job.nodes.size(); + Eigen::MatrixXd KgNoBCDense(Kg.block(0, 0, numberOfDoF, numberOfDoF)); + std::ofstream kgNoBCFile("KgNoBC.csv"); + if (kgNoBCFile.is_open()) { + const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, + Eigen::DontAlignCols, ", ", "\n"); + kgNoBCFile << KgNoBCDense.format(CSVFormat) << '\n'; + kgNoBCFile.close(); + } + std::cout << KgNoBCDense << std::endl; + + // load prescribed boundary conditions into stiffness matrix and force vector + loadBCs(Kg, force_vec, BCs, job.nodes.size()); + + if (equations.size() > 0) { + loadEquations(Kg, equations, job.nodes.size(), BCs.size()); + } + + // load prescribed forces into force vector + if (forces.size() > 0) { + loadForces(force_vec, forces); + } + Eigen::MatrixXd KgDense(Kg); + std::cout << KgDense << std::endl; + Eigen::VectorXd forcesVectorDense(force_vec); + std::cout << forcesVectorDense << std::endl; + + // compress global stiffness matrix since all non-zero values have been added. + Kg.prune(1.e-14); + Kg.makeCompressed(); + std::ofstream kgFile("Kg.csv"); + if (kgFile.is_open()) { + const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, + Eigen::DontAlignCols, ", ", "\n"); + kgFile << KgDense.format(CSVFormat) << '\n'; + kgFile.close(); + } + std::ofstream forcesFile("forces.csv"); + if (forcesFile.is_open()) { + const static Eigen::IOFormat CSVFormat(Eigen::StreamPrecision, + Eigen::DontAlignCols, ", ", "\n"); + forcesFile << forcesVectorDense.format(CSVFormat) << '\n'; + forcesFile.close(); + } + Eigen::FullPivLU lu(Kg); + // assert(lu.isInvertible()); + + // initialize solver based on whether MKL should be used #ifdef EIGEN_USE_MKL_ALL - Eigen::PardisoLU solver; + Eigen::PardisoLU solver; #else - Eigen::SparseLU solver; + Eigen::SparseLU solver; + #endif - //Compute the ordering permutation vector from the structural pattern of Kg - start_time = std::chrono::high_resolution_clock::now(); - solver.analyzePattern(Kg); - end_time = std::chrono::high_resolution_clock::now(); - delta_time = std::chrono::duration_cast(end_time - start_time).count(); + // Compute the ordering permutation vector from the structural pattern of Kg + start_time = std::chrono::high_resolution_clock::now(); + solver.analyzePattern(Kg); + end_time = std::chrono::high_resolution_clock::now(); + delta_time = std::chrono::duration_cast(end_time - + start_time) + .count(); - summary.preprocessing_time_in_ms = delta_time; + summary.preprocessing_time_in_ms = delta_time; - if (options.verbose) - std::cout << "Preprocessing step of factorization completed in " - << delta_time - << " ms.\nNow factorizing global stiffness matrix..." << std::endl; + if (options.verbose) + std::cout << "Preprocessing step of factorization completed in " + << delta_time + << " ms.\nNow factorizing global stiffness matrix..." + << std::endl; - // Compute the numerical factorization - start_time = std::chrono::high_resolution_clock::now(); - solver.factorize(Kg); - end_time = std::chrono::high_resolution_clock::now(); + // Compute the numerical factorization + start_time = std::chrono::high_resolution_clock::now(); + solver.factorize(Kg); + std::cout << solver.lastErrorMessage() << std::endl; + assert(solver.info() == Eigen::Success); + end_time = std::chrono::high_resolution_clock::now(); - delta_time = std::chrono::duration_cast(end_time - start_time).count(); - summary.factorization_time_in_ms = delta_time; + delta_time = std::chrono::duration_cast(end_time - + start_time) + .count(); + summary.factorization_time_in_ms = delta_time; - if (options.verbose) - std::cout << "Factorization completed in " - << delta_time - << " ms. Now solving system..." << std::endl; + if (options.verbose) + std::cout << "Factorization completed in " << delta_time + << " ms. Now solving system..." << std::endl; - //Use the factors to solve the linear system - start_time = std::chrono::high_resolution_clock::now(); - SparseMat dispSparse = solver.solve(force_vec); - end_time = std::chrono::high_resolution_clock::now(); - delta_time = std::chrono::duration_cast(end_time - start_time).count(); + // Use the factors to solve the linear system + start_time = std::chrono::high_resolution_clock::now(); + SparseMat dispSparse = solver.solve(force_vec); + end_time = std::chrono::high_resolution_clock::now(); + delta_time = std::chrono::duration_cast(end_time - + start_time) + .count(); - summary.solve_time_in_ms = delta_time; + summary.solve_time_in_ms = delta_time; - if (options.verbose) - std::cout << "System was solved in " - << delta_time - << " ms.\n" << std::endl; + if (options.verbose) + std::cout << "System was solved in " << delta_time << " ms.\n" << std::endl; - // convert to dense matrix - Eigen::VectorXd disp(dispSparse); + // convert to dense matrix + Eigen::VectorXd disp(dispSparse); - // convert from Eigen vector to std vector - std::vector > disp_vec(job.nodes.size(), std::vector(dofs_per_elem)); - for (size_t i = 0; i < disp_vec.size(); ++i) { - for (unsigned int j = 0; j < dofs_per_elem; ++j) - // round all values close to 0.0 - disp_vec[i][j] = - std::abs(disp(dofs_per_elem * i + j)) < options.epsilon ? 0.0 : disp(dofs_per_elem * i + j); - } - summary.nodal_displacements = disp_vec; + // convert from Eigen vector to std vector + std::vector> disp_vec(job.nodes.size(), + std::vector(dofs_per_elem)); + for (size_t i = 0; i < disp_vec.size(); ++i) { + for (unsigned int j = 0; j < dofs_per_elem; ++j) + // round all values close to 0.0 + disp_vec[i][j] = std::abs(disp(dofs_per_elem * i + j)) < options.epsilon + ? 0.0 + : disp(dofs_per_elem * i + j); + } + summary.nodal_displacements = disp_vec; - // [calculate nodal forces - start_time = std::chrono::high_resolution_clock::now(); + // [calculate nodal forces + start_time = std::chrono::high_resolution_clock::now(); - Kg = Kg.topLeftCorner(dofs_per_elem * job.nodes.size(), dofs_per_elem * job.nodes.size()); - dispSparse = dispSparse.topRows(dofs_per_elem * job.nodes.size()); + Kg = Kg.topLeftCorner(dofs_per_elem * job.nodes.size(), + dofs_per_elem * job.nodes.size()); + dispSparse = dispSparse.topRows(dofs_per_elem * job.nodes.size()); - SparseMat nodal_forces_sparse = Kg * dispSparse; + SparseMat nodal_forces_sparse = Kg * dispSparse; - Eigen::VectorXd nodal_forces_dense(nodal_forces_sparse); + Eigen::VectorXd nodal_forces_dense(nodal_forces_sparse); - std::vector > nodal_forces_vec(job.nodes.size(), std::vector(dofs_per_elem)); - for (size_t i = 0; i < nodal_forces_vec.size(); ++i) { - for (unsigned int j = 0; j < dofs_per_elem; ++j) - // round all values close to 0.0 - nodal_forces_vec[i][j] = std::abs(nodal_forces_dense(dofs_per_elem * i + j)) < options.epsilon ? 0.0 - : nodal_forces_dense( - dofs_per_elem * i + j); - } - summary.nodal_forces = nodal_forces_vec; + std::vector> nodal_forces_vec( + job.nodes.size(), std::vector(dofs_per_elem)); + for (size_t i = 0; i < nodal_forces_vec.size(); ++i) { + for (unsigned int j = 0; j < dofs_per_elem; ++j) + // round all values close to 0.0 + nodal_forces_vec[i][j] = + std::abs(nodal_forces_dense(dofs_per_elem * i + j)) < options.epsilon + ? 0.0 + : nodal_forces_dense(dofs_per_elem * i + j); + } + summary.nodal_forces = nodal_forces_vec; - end_time = std::chrono::high_resolution_clock::now(); + end_time = std::chrono::high_resolution_clock::now(); - summary.nodal_forces_solve_time_in_ms = std::chrono::duration_cast( - end_time - start_time).count(); - //] + summary.nodal_forces_solve_time_in_ms = + std::chrono::duration_cast(end_time - + start_time) + .count(); + //] - // [ calculate forces associated with ties - if (ties.size() > 0) { - start_time = std::chrono::high_resolution_clock::now(); - summary.tie_forces = computeTieForces(ties, disp_vec); - end_time = std::chrono::high_resolution_clock::now(); - summary.tie_forces_solve_time_in_ms = std::chrono::duration_cast( - end_time - start_time).count(); - } - // ] + // [ calculate forces associated with ties + if (ties.size() > 0) { + start_time = std::chrono::high_resolution_clock::now(); + summary.tie_forces = computeTieForces(ties, disp_vec); + end_time = std::chrono::high_resolution_clock::now(); + summary.tie_forces_solve_time_in_ms = + std::chrono::duration_cast(end_time - + start_time) + .count(); + } + // ] - // [save files specified in options - CSVParser csv; - start_time = std::chrono::high_resolution_clock::now(); - if (options.save_nodal_displacements) { - csv.write(options.nodal_displacements_filename, disp_vec, options.csv_precision, options.csv_delimiter); - } + // [save files specified in options + CSVParser csv; + start_time = std::chrono::high_resolution_clock::now(); + if (options.save_nodal_displacements) { + csv.write(options.nodal_displacements_filename, disp_vec, + options.csv_precision, options.csv_delimiter); + } - if (options.save_nodal_forces) { - csv.write(options.nodal_forces_filename, nodal_forces_vec, options.csv_precision, options.csv_delimiter); - } + if (options.save_nodal_forces) { + csv.write(options.nodal_forces_filename, nodal_forces_vec, + options.csv_precision, options.csv_delimiter); + } - if (options.save_tie_forces) { - csv.write(options.tie_forces_filename, summary.tie_forces, options.csv_precision, options.csv_delimiter); - } + if (options.save_tie_forces) { + csv.write(options.tie_forces_filename, summary.tie_forces, + options.csv_precision, options.csv_delimiter); + } - end_time = std::chrono::high_resolution_clock::now(); - summary.file_save_time_in_ms = std::chrono::duration_cast( - end_time - start_time).count(); - // ] + end_time = std::chrono::high_resolution_clock::now(); + summary.file_save_time_in_ms = + std::chrono::duration_cast(end_time - + start_time) + .count(); + // ] - auto final_end_time = std::chrono::high_resolution_clock::now(); + auto final_end_time = std::chrono::high_resolution_clock::now(); - delta_time = std::chrono::duration_cast(final_end_time - initial_start_time).count(); - summary.total_time_in_ms = delta_time; + delta_time = std::chrono::duration_cast( + final_end_time - initial_start_time) + .count(); + summary.total_time_in_ms = delta_time; - if (options.save_report) { - writeStringToTxt(options.report_filename, summary.FullReport()); - } + if (options.save_report) { + writeStringToTxt(options.report_filename, summary.FullReport()); + } - if (options.verbose) - std::cout << summary.FullReport(); + if (options.verbose) + std::cout << summary.FullReport(); - return summary; - }; + // Compute per element forces + summary.element_forces.resize(job.elems.size(), + std::vector(2 * DOF::NUM_DOFS)); + std::vector perElemKlocalAelem = + assembleK3D.getPerElemKlocalAelem(); + + for (size_t elemIndex = 0; elemIndex < job.elems.size(); elemIndex++) { + // Assemble the displacements of the nodes of the beam + Eigen::Matrix elemDisps; + // Displacements of the first node of the beam + const size_t nodeIndex0 = job.elems[elemIndex](0); + elemDisps(0) = summary.nodal_displacements[nodeIndex0][0]; + elemDisps(1) = summary.nodal_displacements[nodeIndex0][1]; + elemDisps(2) = summary.nodal_displacements[nodeIndex0][2]; + elemDisps(3) = summary.nodal_displacements[nodeIndex0][3]; + elemDisps(4) = summary.nodal_displacements[nodeIndex0][4]; + elemDisps(5) = summary.nodal_displacements[nodeIndex0][5]; + // Displacemen ts of the second node of the beam + const size_t nodeIndex1 = job.elems[elemIndex](1); + elemDisps(6) = summary.nodal_displacements[nodeIndex1][0]; + elemDisps(7) = summary.nodal_displacements[nodeIndex1][1]; + elemDisps(8) = summary.nodal_displacements[nodeIndex1][2]; + elemDisps(9) = summary.nodal_displacements[nodeIndex1][3]; + elemDisps(10) = summary.nodal_displacements[nodeIndex1][4]; + elemDisps(11) = summary.nodal_displacements[nodeIndex1][5]; + + Eigen::Matrix elemForces = + perElemKlocalAelem[elemIndex] * elemDisps; + for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) { + summary.element_forces[elemIndex][static_cast(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); + } + } + return summary; +}; } // namespace fea