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

Computation of the element forces and addition of those to the result struct.

This commit is contained in:
Iason 2020-01-21 18:22:19 +01:00
parent a3443cfe66
commit 03f3fe2e52
4 changed files with 1307 additions and 1081 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

644
include/containers.h
View File

@ -35,310 +35,14 @@
#ifndef FEA_CONTAINERS_H
#define FEA_CONTAINERS_H
#include <vector>
#include <Eigen/Core>
#include <vector>
namespace fea {
/**
* @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;/**<The index of the node to constrain*/
/**
* The index of the dof to constrain. The fea::DOF enum can be used for specification or the integer
* values can be used directly `0==d_x`, `1==d_y`, ...
*/
unsigned int dof;
double value;/**<The value to hold the dof at.*/
/**
* @brief Default Constructor
* @details All values are set to zero.
*/
BC() : node(0), dof(0), value(0) { };
/**
* @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 A nodal force to enforce.
* @details Set by specifying the node where the force will be applied, which degree of freedom will be affected,
* and the value to apply.
*
* @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::Force force(nn1, fea::DOF::DISPLACEMENT_X, value);
* @endcode
*/
struct Force {
unsigned int node;/**<The index of the node to apply the force*/
/**
* The index of the dof to constrain. The fea::DOF enum can be used for specification or the integer
* values can be used directly `0==d_x`, `1==d_y`, ...
*/
unsigned int dof;
double value;/**<The value of the force to apply.*/
/**
* @brief Default Constructor
* @details All values are set to zero.
*/
Force() : node(0), dof(0), value(0) { };
/**
* @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 force.
*/
Force(unsigned int _node, unsigned int _dof, double _value) : node(_node), dof(_dof), value(_value) { };
};
/**
* @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<double> normal_vec = {0.0, 0.0, 1.0};
* fea::Props props(EA, EIz, EIy, GJ, normal_vec);
* @endcode
*/
struct Props {
double EA;/**<Extensional stiffness.*/
double EIz;/**<Bending stiffness parallel to local z-axis.*/
double EIy;/**<Bending stiffness parallel to local y-axis.*/
double GJ;/**<Torsional stiffness.*/
Eigen::Vector3d normal_vec;/**<Vector normal to element (size==3). Direction should be parallel to the beam element's local y-axis.*/
Props() : EA(0), EIz(0), EIy(0), GJ(0) { };/**<Default constuctor*/
/**
* @brief Constructor
* @details Allows properties to be set upon initialization.
*
* @param[in] EA double. Extensional stiffness.
* @param[in] EIz double. Bending stiffness parallel to z-axis.
* @param[in] EIy double. Bending stiffness parallel to y-axis.
* @param[in] GJ double. Torsional stiffness.
* @param[in] normal_vec std::vector<double>. 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<double> &_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;/**<The first element's index involved in the constraint.*/
unsigned int node_number_2;/**<The second element's index involved in the constraint.*/
double lmult;/**<multiplier for the linear spring.*/
double rmult; /**<multiplier for the rotational spring.*/
/**
* @brief Default Constructor
* @details All member variables are set to 0.
*/
Tie() : node_number_1(0), node_number_2(0), lmult(0), rmult(0) { };
/**
* @brief Constructor
* @param[in] node_number_1 `unsigned int`. Index of the first node.
* @param[in] node_number_2 `unsigned int`. Index of the second node.
* @param[in] lmult `double`. Spring constant for the translational degrees of freedom.
* @param[in] rmult `double`. Spring constant for the rotational degrees of freedom.
*/
Tie(unsigned int _node_number_1, unsigned int _node_number_2, double _lmult, double _rmult) :
node_number_1(_node_number_1), node_number_2(_node_number_2), lmult(_lmult), rmult(_rmult) { };
};
/**
* @brief A linear multipoint constraint.
* @details Equation constraints are defined by a series of terms that
* sum to zero, e.g. `t1 + t2 + t3 ... = 0`, where `tn` is the `n`th term.
* Each term specifies the node number, degree of freedom and coefficient.
* The node number and degree of freedom specify which nodal variable
* (either nodal displacement or rotation) is involved with the equation
* constraint, and coefficient is multiplied by the specified nodal variable
* when forming the equation. Note, the equation sums to zero, so in order
* to specify that 2 nodal degrees of freedom are equal their coefficients
* should be equal and opposite.
*
* @code
* // Create an empty equation
* fea::Equation eqn;
*
* // Stipulate that the x and y displacement for the first node must be equal
* unsigned int node_number = 0;
* eqn.terms.push_back(fea::Equation::Term(node_number, fea::DOF::DISPLACEMENT_X, 1.0));
* eqn.terms.push_back(fea::Equation::Term(node_number, fea::DOF::DISPLACEMENT_Y, -1.0);
* @endcode
*/
struct Equation {
/**
* @brief A single term in the equation constraint.
* @details Each term defines the node number, degree of freedom and
* coefficient.
*/
struct Term {
unsigned int node_number;/**<Index of the node in the node list.*/
unsigned int dof;/**<Degree of freedom. @sa fea::DOF*/
double coefficient;/**<Coefficient to multiply the nodal variable by.*/
/**
* Default constructor
*/
Term() : node_number(0), dof(0), coefficient(0) {};
/**
* @brief Constructor
* @param node_number. `unsigned int`. Index of the node within the node list.
* @param dof. `unsigned int`. Degree of freedom for the specified node.
* @param coefficient. `double`. coefficient to multiply the corresponding nodal variable by.
*/
Term(unsigned int _node_number, unsigned int _dof, double _coefficient) :
node_number(_node_number), dof(_dof), coefficient(_coefficient) {};
};
std::vector<Term> terms;/**<A list of terms that sum to zero.*/
/**
* Default constructor
*/
Equation() : terms(0) {};
/**
* @brief Constructor
* @param terms. `std::vector<Term>`. A list of terms that sum to zero.k
*/
Equation(const std::vector<Term> &_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;/**<The indices of the node list that define the element.*/
Props props;/**<The set of properties to associate with the element.*/
/**
* @brief Default Constructor
*/
Elem() { };
/**
* @brief Constructor
* @details Constructor if the node indices are passed independently. Assumes 2 nodes per element.
*
* @param[in] node1 unsigned int. The indices of first node associate with the element.
* @param[in] node2 unsigned int. The indices of second node associate with the element.
* @param[in] props Props. The element's properties.
*/
Elem(unsigned int node1, unsigned int node2, const Props &_props) : props(_props) {
node_numbers << node1, node2;
}
};
/**
* @brief Contains a node list, element list, and the properties of each element.
*/
struct Job {
std::vector<Node> nodes;/**<A vector of Node objects that define the mesh.*/
std::vector<Eigen::Vector2i> elems;/**<A 2D vector of ints that defines the connectivity of the node list.*/
std::vector<Props> props;/**<A vector that contains the properties of each element.*/
/**
* @brief Default constructor
*/
Job() : nodes(0), elems(0), props(0) { };
/**
* @brief Constructor
* @details Takes an list of nodes and elements as inputs. The elements are deconstructed into
* an array holding the connectivity list and an array holding the properties.
*
* @param[in] nodes std::vector<Node>. The node list that defines the mesh.
* @param[in] elems std::vector<Elem>. The elements that define the mesh.
* An element is defined by the connectivity list and the associated properties.
*/
Job(const std::vector<Node> &_nodes, const std::vector<Elem> _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.
* @brief Convenience enumerator for specifying the active degree of freedom in
* a constraint.
*/
enum DOF {
/**
@ -376,6 +80,348 @@ namespace fea {
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; /**<The index of the node to constrain*/
/**
* The index of the dof to constrain. The fea::DOF enum can be used for
* specification or the integer values can be used directly `0==d_x`,
* `1==d_y`, ...
*/
unsigned int dof;
double value; /**<The value to hold the dof at.*/
/**
* @brief Default Constructor
* @details All values are set to zero.
*/
BC() : node(0), dof(0), value(0){};
/**
* @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) {
assert(dof < DOF::NUM_DOFS);
};
};
/**
* @brief A nodal force to enforce.
* @details Set by specifying the node where the force will be applied, which
* degree of freedom will be affected, and the value to apply.
*
* @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::Force force(nn1, fea::DOF::DISPLACEMENT_X, value);
* @endcode
*/
struct Force {
unsigned int node; /**<The index of the node to apply the force*/
/**
* The index of the dof to constrain. The fea::DOF enum can be used for
* specification or the integer values can be used directly `0==d_x`,
* `1==d_y`, ...
*/
unsigned int dof;
double value; /**<The value of the force to apply.*/
/**
* @brief Default Constructor
* @details All values are set to zero.
*/
Force() : node(0), dof(0), value(0){};
/**
* @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 force.
*/
Force(unsigned int _node, unsigned int _dof, double _value)
: node(_node), dof(_dof), value(_value) {
assert(dof < DOF::NUM_DOFS);
};
};
/**
* @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<double> normal_vec = {0.0, 0.0, 1.0};
* fea::Props props(EA, EIz, EIy, GJ, normal_vec);
* @endcode
*/
struct Props {
double EA; /**<Extensional stiffness.*/
double EIz; /**<Bending stiffness parallel to local z-axis.*/
double EIy; /**<Bending stiffness parallel to local y-axis.*/
double GJ; /**<Torsional stiffness.*/
Eigen::Vector3d
normal_vec; /**<Vector normal to element (size==3). Direction should be
parallel to the beam element's local y-axis.*/
Props() : EA(0), EIz(0), EIy(0), GJ(0){}; /**<Default constuctor*/
/**
* @brief Constructor
* @details Allows properties to be set upon initialization.
*
* @param[in] EA double. Extensional stiffness.
* @param[in] EIz double. Bending stiffness parallel to z-axis.
* @param[in] EIy double. Bending stiffness parallel to y-axis.
* @param[in] GJ double. Torsional stiffness.
* @param[in] normal_vec std::vector<double>. 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<double> &_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; /**<The first element's index involved in the constraint.*/
unsigned int node_number_2; /**<The second element's index involved in the
constraint.*/
double lmult; /**<multiplier for the linear spring.*/
double rmult; /**<multiplier for the rotational spring.*/
/**
* @brief Default Constructor
* @details All member variables are set to 0.
*/
Tie() : node_number_1(0), node_number_2(0), lmult(0), rmult(0){};
/**
* @brief Constructor
* @param[in] node_number_1 `unsigned int`. Index of the first node.
* @param[in] node_number_2 `unsigned int`. Index of the second node.
* @param[in] lmult `double`. Spring constant for the translational degrees of
* freedom.
* @param[in] rmult `double`. Spring constant for the rotational degrees of
* freedom.
*/
Tie(unsigned int _node_number_1, unsigned int _node_number_2, double _lmult,
double _rmult)
: node_number_1(_node_number_1), node_number_2(_node_number_2),
lmult(_lmult), rmult(_rmult){};
};
/**
* @brief A linear multipoint constraint.
* @details Equation constraints are defined by a series of terms that
* sum to zero, e.g. `t1 + t2 + t3 ... = 0`, where `tn` is the `n`th term.
* Each term specifies the node number, degree of freedom and coefficient.
* The node number and degree of freedom specify which nodal variable
* (either nodal displacement or rotation) is involved with the equation
* constraint, and coefficient is multiplied by the specified nodal variable
* when forming the equation. Note, the equation sums to zero, so in order
* to specify that 2 nodal degrees of freedom are equal their coefficients
* should be equal and opposite.
*
* @code
* // Create an empty equation
* fea::Equation eqn;
*
* // Stipulate that the x and y displacement for the first node must be equal
* unsigned int node_number = 0;
* eqn.terms.push_back(fea::Equation::Term(node_number,
* fea::DOF::DISPLACEMENT_X, 1.0));
* eqn.terms.push_back(fea::Equation::Term(node_number,
* fea::DOF::DISPLACEMENT_Y, -1.0);
* @endcode
*/
struct Equation {
/**
* @brief A single term in the equation constraint.
* @details Each term defines the node number, degree of freedom and
* coefficient.
*/
struct Term {
unsigned int node_number; /**<Index of the node in the node list.*/
unsigned int dof; /**<Degree of freedom. @sa fea::DOF*/
double coefficient; /**<Coefficient to multiply the nodal variable by.*/
/**
* Default constructor
*/
Term() : node_number(0), dof(0), coefficient(0){};
/**
* @brief Constructor
* @param node_number. `unsigned int`. Index of the node within the node
* list.
* @param dof. `unsigned int`. Degree of freedom for the specified node.
* @param coefficient. `double`. coefficient to multiply the corresponding
* nodal variable by.
*/
Term(unsigned int _node_number, unsigned int _dof, double _coefficient)
: node_number(_node_number), dof(_dof), coefficient(_coefficient){};
};
std::vector<Term> terms; /**<A list of terms that sum to zero.*/
/**
* Default constructor
*/
Equation() : terms(0){};
/**
* @brief Constructor
* @param terms. `std::vector<Term>`. A list of terms that sum to zero.k
*/
Equation(const std::vector<Term> &_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; /**<The indices of the node list that define the element.*/
Props props; /**<The set of properties to associate with the element.*/
/**
* @brief Default Constructor
*/
Elem(){};
/**
* @brief Constructor
* @details Constructor if the node indices are passed independently. Assumes
* 2 nodes per element.
*
* @param[in] node1 unsigned int. The indices of first node associate with the
* element.
* @param[in] node2 unsigned int. The indices of second node associate with
* the element.
* @param[in] props Props. The element's properties.
*/
Elem(unsigned int node1, unsigned int node2, const Props &_props)
: props(_props) {
node_numbers << node1, node2;
}
};
/**
* @brief Contains a node list, element list, and the properties of each
* element.
*/
struct Job {
std::vector<Node> nodes; /**<A vector of Node objects that define the mesh.*/
std::vector<Eigen::Vector2i> elems; /**<A 2D vector of ints that defines the
connectivity of the node list.*/
std::vector<Props>
props; /**<A vector that contains the properties of each element.*/
/**
* @brief Default constructor
*/
Job() : nodes(0), elems(0), props(0){};
/**
* @brief Constructor
* @details Takes an list of nodes and elements as inputs. The elements are
* deconstructed into an array holding the connectivity list and an array
* holding the properties.
*
* @param[in] nodes std::vector<Node>. The node list that defines the mesh.
* @param[in] elems std::vector<Elem>. The elements that define the mesh.
* An element is defined by the connectivity list and the
* associated properties.
*/
Job(const std::vector<Node> &_nodes, const std::vector<Elem> _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);
}
};
};
} // namespace fea
#endif // FEA_CONTAINERS_H

19
include/summary.h
View File

@ -57,8 +57,8 @@ namespace fea {
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.
* 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;
@ -84,8 +84,9 @@ namespace fea {
/**
* 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.
* 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;
@ -143,10 +144,16 @@ namespace fea {
*/
std::vector<std::vector<double>> 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<std::vector<double>> element_forces;
};
} // namespace fea
#endif // THREEDBEAMFEA_SUMMARY_H

180
include/threed_beam_fea.h
View File

@ -46,9 +46,9 @@
#include <Eigen/SparseCore>
#include "containers.h"
#include "csv_parser.h"
#include "options.h"
#include "summary.h"
#include "csv_parser.h"
namespace fea {
@ -58,13 +58,15 @@ namespace fea {
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> GlobalStiffMatrix;
/**
* An elemental matrix in local coordinates. Will either be the elemental stiffness matrix or the global-to-local rotation matrix
* An elemental matrix in local coordinates. Will either be the elemental
* stiffness matrix or the global-to-local rotation matrix
*/
typedef Eigen::Matrix<double, 12, 12, Eigen::RowMajor> 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
* 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<double, Eigen::Dynamic, 1> ForceVector;
@ -75,7 +77,8 @@ namespace fea {
/**
* @brief Calculates the distance between 2 nodes.
* @details Calculates the original Euclidean distance between 2 nodes in the x-y plane.
* @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.
@ -90,7 +93,6 @@ namespace fea {
class GlobalStiffAssembler {
public:
/**
* @brief Default constructor.
* @details Initializes all entries in member matrices to 0.0.
@ -106,32 +108,40 @@ namespace fea {
/**
* @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.
* @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<fea::Tie>`. Vector of ties that apply to attach springs of specified stiffness to
* all nodal degrees of freedom between each set of nodes indicated.
* @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<fea::Tie>`. 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<Tie> &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.
* @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.
* @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.
* @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);
@ -139,99 +149,125 @@ namespace fea {
* @brief Returns the currently stored elemental stiffness matrix.
* @return <B>Elemental stiffness matrix</B> `fea::LocalMatrix`.
*/
LocalMatrix getKelem() {
return Kelem;
}
LocalMatrix getKelem() { return Kelem; }
/**
* @brief Returns the currently stored rotation matrix.
* @return <B>Rotation matrix</B> `fea::LocalMatrix`.
*/
LocalMatrix getAelem() {
return Aelem;
}
LocalMatrix getAelem() { return Aelem; }
std::vector<LocalMatrix> getPerElemKlocalAelem() const;
private:
LocalMatrix Kelem;
/**<Elemental stiffness matrix in global coordinate system.*/
LocalMatrix Klocal;
/**<Elemental stiffness matrix in local coordinate system (used as temporary variable).*/
/**<Elemental stiffness matrix in local coordinate system (used as temporary
* variable).*/
LocalMatrix Aelem;
/**<Rotation matrix. Transforms `Klocal` to global coorinate system (`Kelem`).*/
/**<Rotation matrix. Transforms `Klocal` to global coorinate system
* (`Kelem`).*/
LocalMatrix AelemT;
/**<Transposed rotation matrix.*/
SparseMat SparseKelem;/**<Sparse representation of elemental stiffness matrix.*/
SparseMat
SparseKelem; /**<Sparse representation of elemental stiffness matrix.*/
std::vector<LocalMatrix>
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
};
/**
* @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.
* @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<fea::BC>`. 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.
* @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<fea::BC>`. 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<BC> &BCs, unsigned int num_nodes);
void loadBCs(SparseMat &Kg, ForceVector &force_vec, const std::vector<BC> &BCs,
unsigned int num_nodes);
void loadEquations(SparseMat &Kg, const std::vector<Equation> &equations, unsigned int num_nodes, unsigned int num_bcs);
void loadEquations(SparseMat &Kg, const std::vector<Equation> &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.
* @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<fea::Tie>`. Vector of `Tie`'s to apply to the current analysis.
* @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<fea::Tie>`. Vector of `Tie`'s to apply to the
* current analysis.
*/
void loadTies(std::vector<Eigen::Triplet<double> > &triplets, const std::vector<Tie> &ties);
void loadTies(std::vector<Eigen::Triplet<double>> &triplets,
const std::vector<Tie> &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`.
* @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<Tie>`. 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.
* @param[in] ties `std::vector<Tie>`. 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<std::vector<double> > computeTieForces(const std::vector<Tie> &ties,
std::vector<std::vector<double>>
computeTieForces(const std::vector<Tie> &ties,
const std::vector<std::vector<double>> &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<Force>. Vector of prescribed forces to apply to the current analysis.
* @param force_vec `ForceVector`. Right hand side of the \f$[K][Q]=[F]\f$
* equation of the FE analysis.
* @param[in] forces std::vector<Force>. Vector of prescribed forces to apply to
* the current analysis.
*/
void loadForces(SparseMat &force_vec, const std::vector<Force> &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.
* @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<fea::BC>`. Vector of boundary conditions to apply to the nodal degrees of freedom contained in the job.
* @param[in] forces `std::vector<fea::Force>`. Vector of prescribed forces to apply to the nodal degrees of freedom contained in the job.
* @param[in] ties `std::vector<fea::Tie>`. Vector of ties that apply to attach springs of specified stiffness to
* all nodal degrees of freedom between each set of nodes indicated.
* @param[in] job `fea::Job`. Contains the node, element, and property lists for
* the mesh.
* @param[in] BCs `std::vector<fea::BC>`. Vector of boundary conditions to apply
* to the nodal degrees of freedom contained in the job.
* @param[in] forces `std::vector<fea::Force>`. Vector of prescribed forces to
* apply to the nodal degrees of freedom contained in the job.
* @param[in] ties `std::vector<fea::Tie>`. Vector of ties that apply to attach
* springs of specified stiffness to all nodal degrees of freedom between each
* set of nodes indicated.
*
* @return <B>Summary</B> `fea::Summary`. Summary containing the results of the analysis.
* @return <B>Summary</B> `fea::Summary`. Summary containing the results of the
* analysis.
*/
Summary solve(const Job &job,
const std::vector<BC> &BCs,
const std::vector<Force> &forces,
const std::vector<Tie> &ties,
const std::vector<Equation> &equations,
const Options &options);
Summary solve(const Job &job, const std::vector<BC> &BCs,
const std::vector<Force> &forces, const std::vector<Tie> &ties,
const std::vector<Equation> &equations, const Options &options);
} // namespace fea
#endif // THREED_BEAM_FEA_H

309
src/threed_beam_fea.cpp
View File

@ -23,6 +23,7 @@
//
// Author: ryan.latture@gmail.com (Ryan Latture)
#include <Eigen/LU>
#include <boost/format.hpp>
#include <chrono>
#include <cmath>
@ -43,13 +44,12 @@ namespace fea {
if (!output_file.is_open()) {
throw std::runtime_error(
(boost::format("Error opening file %s.") % filename).str()
);
(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;
@ -58,9 +58,9 @@ namespace fea {
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 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
@ -70,11 +70,12 @@ namespace fea {
// 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
// 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 tmp12z = 12.0 * EIz / (length * length * length); // ==k3
const double tmp6z = 6.0 * EIz / (length * length);
const double tmp1z = EIz / length;
@ -133,17 +134,36 @@ namespace fea {
// 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) {
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;
nz = nx.cross(ny).normalized();
// update rotation matrix
Aelem(0, 0) = nx(0);
@ -186,6 +206,15 @@ namespace fea {
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);
@ -226,20 +255,27 @@ namespace fea {
AelemT(11, 9) = nx(2);
AelemT(11, 10) = ny(2);
AelemT(11, 11) = nz(2);
}
std::vector<LocalMatrix> GlobalStiffAssembler::getPerElemKlocalAelem() const {
return perElemKlocalAelem;
};
void GlobalStiffAssembler::operator()(SparseMat &Kg, const Job &job, const std::vector<Tie> &ties) {
void GlobalStiffAssembler::operator()(SparseMat &Kg, const Job &job,
const std::vector<Tie> &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
// form vector to hold triplets that will be used to assemble global stiffness
// matrix
std::vector<Eigen::Triplet<double>> 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);
calcKelem(i, job); // 12x12 matrix
// get sparse representation of the current elemental stiffness matrix
SparseKelem = Kelem.sparseView();
@ -253,7 +289,8 @@ namespace fea {
row = it.row();
col = it.col();
// check position in local matrix and update corresponding global position
// check position in local matrix and update corresponding global
// position
if (row < 6) {
// top left
if (col < 6) {
@ -263,23 +300,23 @@ namespace fea {
}
// top right
else {
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * nn1 + row,
dofs_per_elem * (nn2 - 1) + col,
it.value()));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * nn1 + row, dofs_per_elem * (nn2 - 1) + col,
it.value())); // I: Giati nn2-1 kai oxi nn2?
}
}
else {
} else {
// bottom left
if (col < 6) {
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * (nn2 - 1) + row,
dofs_per_elem * nn1 + col,
it.value()));
triplets.push_back(Eigen::Triplet<double>(
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<double>(dofs_per_elem * (nn2 - 1) + row,
dofs_per_elem * (nn2 - 1) + col,
it.value()));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * (nn2 - 1) + row,
dofs_per_elem * (nn2 - 1) + col, it.value()));
}
}
}
@ -291,10 +328,12 @@ namespace fea {
Kg.setFromTriplets(triplets.begin(), triplets.end());
};
void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector<BC> &BCs, unsigned int num_nodes) {
void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector<BC> &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
// 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) {
@ -304,14 +343,16 @@ namespace fea {
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.
// update force vector. All values are already zero. Only update if BC if
// non-zero.
if (std::abs(BCs[i].value) > std::numeric_limits<double>::epsilon()) {
force_vec.insert(global_add_idx + i, 0) = BCs[i].value;
}
}
};
void loadEquations(SparseMat &Kg, const std::vector<Equation> &equations, unsigned int num_nodes, unsigned int num_bcs) {
void loadEquations(SparseMat &Kg, const std::vector<Equation> &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;
@ -319,14 +360,16 @@ namespace fea {
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;
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<Eigen::Triplet<double> > &triplets, const std::vector<Tie> &ties) {
void loadTies(std::vector<Eigen::Triplet<double>> &triplets,
const std::vector<Tie> &ties) {
const unsigned int dofs_per_elem = DOF::NUM_DOFS;
unsigned int nn1, nn2;
double lmult, rmult, spring_constant;
@ -338,35 +381,34 @@ namespace fea {
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
// first 3 DOFs are linear DOFs, second 2 are rotational, last is
// torsional
spring_constant = j < 3 ? lmult : rmult;
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * nn1 + j,
dofs_per_elem * nn1 + j,
spring_constant));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * nn1 + j, dofs_per_elem * nn1 + j, spring_constant));
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * nn2 + j,
dofs_per_elem * nn2 + j,
spring_constant));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * nn2 + j, dofs_per_elem * nn2 + j, spring_constant));
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * nn1 + j,
dofs_per_elem * nn2 + j,
-spring_constant));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * nn1 + j, dofs_per_elem * nn2 + j, -spring_constant));
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * nn2 + j,
dofs_per_elem * nn1 + j,
-spring_constant));
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * nn2 + j, dofs_per_elem * nn1 + j, -spring_constant));
}
}
};
std::vector<std::vector<double> > computeTieForces(const std::vector<Tie> &ties,
std::vector<std::vector<double>>
computeTieForces(const std::vector<Tie> &ties,
const std::vector<std::vector<double>> &nodal_displacements) {
const unsigned int dofs_per_elem = DOF::NUM_DOFS;
unsigned int nn1, nn2;
double lmult, rmult, spring_constant, delta1, delta2;
std::vector<std::vector<double> > tie_forces(ties.size(), std::vector<double>(dofs_per_elem));
std::vector<std::vector<double>> tie_forces(
ties.size(), std::vector<double>(dofs_per_elem));
for (size_t i = 0; i < ties.size(); ++i) {
nn1 = ties[i].node_number_1;
@ -375,7 +417,8 @@ namespace fea {
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
// 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];
@ -395,12 +438,9 @@ namespace fea {
}
};
Summary solve(const Job &job,
const std::vector<BC> &BCs,
const std::vector<Force> &forces,
const std::vector<Tie> &ties,
const std::vector<Equation> &equations,
const Options &options) {
Summary solve(const Job &job, const std::vector<BC> &BCs,
const std::vector<Force> &forces, const std::vector<Tie> &ties,
const std::vector<Equation> &equations, const Options &options) {
auto initial_start_time = std::chrono::high_resolution_clock::now();
Summary summary;
@ -412,7 +452,8 @@ namespace fea {
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();
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);
@ -424,14 +465,26 @@ namespace fea {
GlobalStiffAssembler assembleK3D = GlobalStiffAssembler();
assembleK3D(Kg, job, ties);
auto end_time = std::chrono::high_resolution_clock::now();
auto delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count();
auto delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(
end_time - start_time)
.count();
summary.assembly_time_in_ms = delta_time;
if (options.verbose)
std::cout << "Global stiffness matrix assembled in "
<< delta_time
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());
@ -443,94 +496,128 @@ namespace fea {
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<Eigen::MatrixXd> lu(Kg);
// assert(lu.isInvertible());
// initialize solver based on whether MKL should be used
#ifdef EIGEN_USE_MKL_ALL
Eigen::PardisoLU<SparseMat> solver;
#else
Eigen::SparseLU<SparseMat> 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<std::chrono::milliseconds>(end_time - start_time).count();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
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;
<< " ms.\nNow factorizing global stiffness matrix..."
<< std::endl;
// 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<std::chrono::milliseconds>(end_time - start_time).count();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
summary.factorization_time_in_ms = delta_time;
if (options.verbose)
std::cout << "Factorization completed in "
<< delta_time
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<std::chrono::milliseconds>(end_time - start_time).count();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
summary.solve_time_in_ms = delta_time;
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;
// convert to dense matrix
Eigen::VectorXd disp(dispSparse);
// convert from Eigen vector to std vector
std::vector<std::vector<double> > disp_vec(job.nodes.size(), std::vector<double>(dofs_per_elem));
std::vector<std::vector<double>> disp_vec(job.nodes.size(),
std::vector<double>(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);
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();
Kg = Kg.topLeftCorner(dofs_per_elem * job.nodes.size(), 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;
Eigen::VectorXd nodal_forces_dense(nodal_forces_sparse);
std::vector<std::vector<double> > nodal_forces_vec(job.nodes.size(), std::vector<double>(dofs_per_elem));
std::vector<std::vector<double>> nodal_forces_vec(
job.nodes.size(), std::vector<double>(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);
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();
summary.nodal_forces_solve_time_in_ms = std::chrono::duration_cast<std::chrono::milliseconds>(
end_time - start_time).count();
summary.nodal_forces_solve_time_in_ms =
std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
//]
// [ calculate forces associated with ties
@ -538,8 +625,10 @@ namespace fea {
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<std::chrono::milliseconds>(
end_time - start_time).count();
summary.tie_forces_solve_time_in_ms =
std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
}
// ]
@ -547,25 +636,32 @@ namespace fea {
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);
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);
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);
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<std::chrono::milliseconds>(
end_time - start_time).count();
summary.file_save_time_in_ms =
std::chrono::duration_cast<std::chrono::milliseconds>(end_time -
start_time)
.count();
// ]
auto final_end_time = std::chrono::high_resolution_clock::now();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(final_end_time - initial_start_time).count();
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(
final_end_time - initial_start_time)
.count();
summary.total_time_in_ms = delta_time;
if (options.save_report) {
@ -575,6 +671,47 @@ namespace fea {
if (options.verbose)
std::cout << summary.FullReport();
// Compute per element forces
summary.element_forces.resize(job.elems.size(),
std::vector<double>(2 * DOF::NUM_DOFS));
std::vector<LocalMatrix> perElemKlocalAelem =
assembleK3D.getPerElemKlocalAelem();
for (size_t elemIndex = 0; elemIndex < job.elems.size(); elemIndex++) {
// Assemble the displacements of the nodes of the beam
Eigen::Matrix<double, 12, 1> 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<double, 12, 1> elemForces =
perElemKlocalAelem[elemIndex] * elemDisps;
for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) {
summary.element_forces[elemIndex][static_cast<size_t>(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<size_t>(dofIndex)] =
//-elemForces(dofIndex);
+elemForces(dofIndex);
}
}
return summary;
};