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

Refactoring. Moved the calculation of the stiffness matrix to a function for avoid recomputation when needed.

This commit is contained in:
iasonmanolas 2021-12-20 19:04:08 +02:00
parent 13096dbd6c
commit dd894b7685
6 changed files with 472 additions and 407 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists.txt
View File

@ -1,7 +1,8 @@
cmake_minimum_required(VERSION 2.8.11)
cmake_minimum_required(VERSION 3.0)
project(ThreedBeamFEA)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++11")
#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++11")
#include(FindOpenMP)
#if(OPENMP_FOUND)
@ -27,3 +28,5 @@ target_include_directories(${PROJECT_NAME}
PUBLIC ${EXT_RAPIDJSON_ROOT}
PUBLIC ${EXT_TCLAP_ROOT}
)
#set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE TRUE)

36
include/containers.h
View File

@ -386,19 +386,19 @@ struct Elem {
* @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
struct BeamStructure
{
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.*/
std::vector<Props> props; /**<A vector that contains the properties of each element.*/
/**
/**
* @brief Default constructor
*/
Job() : nodes(0), elems(0), props(0){};
BeamStructure() : 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
@ -409,17 +409,17 @@ struct Job {
* 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);
BeamStructure(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);
}
};
for (unsigned int i = 0; i < num_elems; i++) {
elems.push_back(_elems[i].node_numbers);
props.push_back(_elems[i].props);
}
};
};
} // namespace fea

2
include/setup.h
View File

@ -101,7 +101,7 @@ namespace fea {
* nodes, elements, and properties.
* @return Job. `fea::Job`.
*/
Job createJobFromJSON(const rapidjson::Document &config_doc);
BeamStructure createJobFromJSON(const rapidjson::Document &config_doc);
/**
* Creates an `fea::Options` object from the configuration document. Any options provided will override the

44
include/threed_beam_fea.h
View File

@ -126,7 +126,7 @@ public:
* 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);
void operator()(SparseMat &Kg, const BeamStructure &job, const std::vector<Tie> &ties);
/**
* @brief Updates the elemental stiffness matrix for the `ith` element.
@ -136,7 +136,10 @@ public:
* @param[in] job `Job`. Current `fea::Job` to analyze contains node, element,
* and property lists.
*/
void calcKelem(unsigned int i, const Job &job);
void calcKelem(const size_t &elementIndex,
const Props &elementProps,
const Node &n1,
const Node &n2);
/**
* @brief Updates the rotation and transposed rotation matrices.
@ -162,9 +165,10 @@ public:
*/
LocalMatrix getAelem() { return Aelem; }
std::vector<LocalMatrix> getPerElemKlocalAelem() const;
//FIX: This allows the client to change the data member. I think there is a proble with eigen..
std::vector<LocalMatrix> &getPerElemKlocalAelem();
private:
private:
LocalMatrix Kelem;
/**<Elemental stiffness matrix in global coordinate system.*/
LocalMatrix Klocal;
@ -185,6 +189,10 @@ private:
// order to find the per element forces
};
void assembleStiffnessMatrix(const BeamStructure &structure,
long long &timeToAssemble,
std::vector<LocalMatrix> &perElemKlocalAelem,
SparseMat &Kg);
/**
* @brief Loads the boundary conditions into the global stiffness matrix and
* force vector.
@ -250,7 +258,9 @@ computeTieForces(const std::vector<Tie> &ties,
*/
void loadForces(SparseMat &force_vec, const std::vector<Force> &forces);
/**
struct ThreedBeamFEA
{
/**
* @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
@ -270,9 +280,27 @@ void loadForces(SparseMat &force_vec, const std::vector<Force> &forces);
* @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 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 BeamStructure &structure,
const std::vector<BC> &BCs,
const std::vector<Force> &forces,
const std::vector<Tie> &ties,
const std::vector<Equation> &equations,
const Options &options);
void initialize(const BeamStructure &structure);
std::vector<LocalMatrix> perElemKlocalAelem;
SparseMat Kg;
bool initialized{false};
BeamStructure structure;
public:
bool wasInitialized() const;
};
} // namespace fea
#endif // THREED_BEAM_FEA_H

4
src/setup.cpp
View File

@ -207,10 +207,10 @@ namespace fea {
return eqns_out;
}
Job createJobFromJSON(const rapidjson::Document &config_doc) {
BeamStructure createJobFromJSON(const rapidjson::Document &config_doc) {
std::vector<Node> nodes = createNodeVecFromJSON(config_doc);
std::vector<Elem> elems = createElemVecFromJSON(config_doc);
return Job(nodes, elems);
return BeamStructure(nodes, elems);
}
Options createOptionsFromJSON(const rapidjson::Document &config_doc) {

786
src/threed_beam_fea.cpp
View File

@ -58,90 +58,90 @@ inline double norm(const Node &n1, const Node &n2) {
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;
void GlobalStiffAssembler::calcKelem(const size_t &elementIndex,
const Props &elementProperties,
const Node &n1,
const Node &n2)
{
// extract element properties
const double EA = elementProperties.EA; // Young's modulus * cross area
const double EIz = elementProperties.EIz; // Young's modulus* I3
const double EIy = elementProperties.EIy; // Young's modulus* I2
const double GJ = elementProperties.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(n1, n2);
// 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;
// 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 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;
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;
// 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
const Eigen::Vector3d nx = (job.nodes[nn2] - job.nodes[nn1]).normalized();
// calculate unit normal vector along y-direction
// const Eigen::Vector3d ny = job.props[i].normal_vec.normalized();
const Eigen::Vector3d ny = job.props[i].normal_vec.normalized().cross(nx).normalized();
// calculate the unit normal vector in local z direction
// const Eigen::Vector3d nz = nx.cross(ny).normalized();
const Eigen::Vector3d nz = nx.cross(ny).normalized();
RotationMatrix r;
r.row(0) = nx;
r.row(1) = ny;
r.row(2) = nz;
// update rotation matrices
calcAelem(r);
AelemT = Aelem.transpose();
// calculate unit normal vector along local x-direction
const Eigen::Vector3d nx = (n2 - n1).normalized();
// calculate unit normal vector along y-direction
// const Eigen::Vector3d ny = job.props[i].normal_vec.normalized();
const Eigen::Vector3d ny = elementProperties.normal_vec.normalized().cross(nx).normalized();
// calculate the unit normal vector in local z direction
// const Eigen::Vector3d nz = nx.cross(ny).normalized();
const Eigen::Vector3d nz = nx.cross(ny).normalized();
RotationMatrix r;
r.row(0) = nx;
r.row(1) = ny;
r.row(2) = nz;
// update rotation matrices
calcAelem(r);
AelemT = Aelem.transpose();
#ifdef DEBUG
std::ofstream KlocalFile("Klocal.csv");
@ -165,8 +165,8 @@ void GlobalStiffAssembler::calcKelem(unsigned int i, const Job &job) {
LocalMatrix KlocalAelem;
KlocalAelem = Klocal * Aelem;
perElemKlocalAelem[i] = KlocalAelem;
};
perElemKlocalAelem[elementIndex] = KlocalAelem;
}
void GlobalStiffAssembler::calcAelem(const RotationMatrix &r) {
// update rotation matrix
@ -176,76 +176,87 @@ void GlobalStiffAssembler::calcAelem(const RotationMatrix &r) {
Aelem.block<3, 3>(9, 9) = r;
}
std::vector<LocalMatrix> GlobalStiffAssembler::getPerElemKlocalAelem() const {
return perElemKlocalAelem;
};
std::vector<LocalMatrix> &GlobalStiffAssembler::getPerElemKlocalAelem()
{
return perElemKlocalAelem;
}
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;
void GlobalStiffAssembler::operator()(SparseMat &Kg,
const BeamStructure &job,
const std::vector<Tie> &ties)
{
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<Eigen::Triplet<double>> triplets;
triplets.reserve(40 * job.elems.size() + 4 * dofs_per_elem * ties.size());
// 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); // 12x12 matrix
perElemKlocalAelem.resize(job.elems.size());
// #ifdef ALL_ELEMENTS_HAVE_THE_SAME_PROPS
for (unsigned int elementIndex = 0; elementIndex < job.elems.size();
++elementIndex) { //NOTE: should be uncommented if each element has a different Kelem
// update Kelem with current elemental stiffness matrix
const Props &elementProperties = job.props[elementIndex];
// get sparse representation of the current elemental stiffness matrix
SparseKelem = Kelem.sparseView();
// store node indices of current element
const int &ni1 = job.elems[elementIndex][0];
const int &ni2 = job.elems[elementIndex][1];
const Node &n1 = job.nodes[ni1];
const Node &n2 = job.nodes[ni2];
calcKelem(elementIndex, elementProperties, n1, n2); // 12x12 matrix
// pull out current node indices
nn1 = job.elems[i][0];
nn2 = job.elems[i][1];
// get sparse representation of the current elemental stiffness matrix
SparseKelem = Kelem.sparseView();
// for (unsigned int i = 0; i < job.elems.size(); ++i) {
// pull out current node indices
for (unsigned int j = 0; j < SparseKelem.outerSize(); ++j) {
for (SparseMat::InnerIterator it(SparseKelem, j); it; ++it) {
row = it.row();
col = it.col();
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<double>(dofs_per_elem * nn1 + row,
dofs_per_elem * nn1 + col,
it.value()));
}
// top right
else {
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 {
// bottom left
if (col < 6) {
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()));
}
// check position in local matrix and update corresponding global
// position
if (row < 6) {
// top left
if (col < 6) {
triplets.push_back(Eigen::Triplet<double>(dofs_per_elem * ni1 + row,
dofs_per_elem * ni1 + col,
it.value()));
}
// top right
else {
triplets.push_back(
Eigen::Triplet<double>(dofs_per_elem * ni1 + row,
dofs_per_elem * (ni2 - 1) + col,
it.value())); // I: Giati nn2-1 kai oxi nn2?
}
} else {
// bottom left
if (col < 6) {
triplets.push_back(Eigen::Triplet<double>(
dofs_per_elem * (ni2 - 1) + row,
dofs_per_elem * ni1 + 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 * (ni2 - 1) + row,
dofs_per_elem * (ni2 - 1) + col,
it.value()));
}
}
}
}
}
}
}
loadTies(triplets, ties);
loadTies(triplets, ties);
Kg.setFromTriplets(triplets.begin(), triplets.end());
};
Kg.setFromTriplets(triplets.begin(), triplets.end());
}
void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector<BC> &BCs,
unsigned int num_nodes) {
@ -268,7 +279,7 @@ void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector<BC> &BCs,
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) {
@ -285,7 +296,7 @@ void loadEquations(SparseMat &Kg, const std::vector<Equation> &equations,
Kg.insert(col_idx, row_idx) = equations[i].terms[j].coefficient;
}
}
};
}
void loadTies(std::vector<Eigen::Triplet<double>> &triplets,
const std::vector<Tie> &ties) {
@ -355,61 +366,69 @@ void loadForces(SparseMat &force_vec, const std::vector<Force> &forces) {
idx = dofs_per_elem * forces[i].node + forces[i].dof;
force_vec.insert(idx, 0) = forces[i].value;
}
};
}
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();
void assembleStiffnessMatrix(const BeamStructure &structure,
long long &timeToAssemble,
std::vector<LocalMatrix> &perElemKlocalAelem,
SparseMat &Kg)
{
auto start_time = std::chrono::high_resolution_clock::now();
GlobalStiffAssembler assembleK3D = GlobalStiffAssembler();
assembleK3D(Kg, structure, {});
perElemKlocalAelem.clear();
perElemKlocalAelem = assembleK3D.getPerElemKlocalAelem();
auto end_time = std::chrono::high_resolution_clock::now();
timeToAssemble = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time)
.count();
}
Summary summary;
summary.num_nodes = job.nodes.size();
summary.num_elems = job.elems.size();
summary.num_bcs = BCs.size();
summary.num_ties = ties.size();
Summary ThreedBeamFEA::solve(const std::vector<BC> &BCs,
const std::vector<Force> &forces,
const std::vector<Tie> &ties,
const std::vector<Equation> &equations,
const Options &options)
{
assert(wasInitialized());
auto initial_start_time = std::chrono::high_resolution_clock::now();
const unsigned int dofs_per_elem = DOF::NUM_DOFS;
Summary summary;
summary.num_nodes = structure.nodes.size();
summary.num_elems = structure.elems.size();
summary.num_bcs = BCs.size();
summary.num_ties = ties.size();
// 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 int dofs_per_elem = DOF::NUM_DOFS;
// construct global stiffness matrix and force vector
SparseMat Kg(size, size);
SparseMat force_vec(size, 1);
force_vec.reserve(forces.size() + BCs.size());
// calculate size of global stiffness matrix and force vector
const unsigned long size = dofs_per_elem * structure.nodes.size() + BCs.size()
+ equations.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<std::chrono::milliseconds>(
end_time - start_time)
.count();
summary.assembly_time_in_ms = delta_time;
// construct global stiffness matrix and force vector
SparseMat force_vec(size, 1);
force_vec.reserve(forces.size() + BCs.size());
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();
SparseMat Kg_local = Kg;
Kg_local.conservativeResize(size, size); //NOTE: Does this remove existing elements?
unsigned int numberOfDoF = DOF::NUM_DOFS * structure.nodes.size();
#ifdef DEBUG
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();
}
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;
#endif
// load prescribed boundary conditions into stiffness matrix and force vector
loadBCs(Kg, force_vec, BCs, job.nodes.size());
loadBCs(Kg_local, force_vec, BCs, structure.nodes.size());
if (equations.size() > 0) {
loadEquations(Kg, equations, job.nodes.size(), BCs.size());
loadEquations(Kg_local, equations, structure.nodes.size(), BCs.size());
}
// load prescribed forces into force vector
@ -423,8 +442,8 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
// std::cout << forcesVectorDense << std::endl;
#endif
// compress global stiffness matrix since all non-zero values have been added.
Kg.prune(1.e-14);
Kg.makeCompressed();
Kg_local.prune(1.e-14);
Kg_local.makeCompressed();
#ifdef DEBUG
std::ofstream kgFile("Kg.csv");
if (kgFile.is_open()) {
@ -450,214 +469,229 @@ Summary solve(const Job &job, const std::vector<BC> &BCs,
#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();
auto start_time = std::chrono::high_resolution_clock::now();
solver.analyzePattern(Kg_local);
auto end_time = std::chrono::high_resolution_clock::now();
auto 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;
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);
solver.factorize(Kg_local);
#ifdef DEBUG
Eigen::VectorXd kgCol(Kg.col(7));
Eigen::FullPivLU<Eigen::MatrixXd> kg_invCheck(Kg);
assert(kg_invCheck.isInvertible());
// Eigen::VectorXd kgCol(Kg.col(7));
// Eigen::FullPivLU<Eigen::MatrixXd> kg_invCheck(Kg);
// assert(kg_invCheck.isInvertible());
#endif
end_time = std::chrono::high_resolution_clock::now();
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
<< " ms. Now solving system..." << std::endl;
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
assert(solver.info() == Eigen::Success);
return summary;
}
// Use the factors to solve the linear system
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();
summary.solve_time_in_ms = delta_time;
if (options.verbose)
std::cout << "System was solved in " << delta_time << " ms.\n" << std::endl;
assert(solver.info() == Eigen::Success);
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
return summary;
}
// convert to dense matrix
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));
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();
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));
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();
summary.nodal_forces_solve_time_in_ms =
std::chrono::duration_cast<std::chrono::milliseconds>(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<std::chrono::milliseconds>(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) {
std::cout << "Writing to:" + options.nodal_displacements_filename
<< std::endl;
csv.write(options.nodal_displacements_filename, disp_vec,
options.csv_precision, options.csv_delimiter);
}
delta_time = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count();
summary.factorization_time_in_ms = delta_time;
if (options.save_nodal_forces) {
csv.write(options.nodal_forces_filename, nodal_forces_vec,
options.csv_precision, options.csv_delimiter);
}
if (options.verbose)
std::cout << "Factorization completed in " << delta_time << " ms. Now solving system..."
<< std::endl;
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<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();
summary.total_time_in_ms = delta_time;
if (options.save_report) {
writeStringToTxt(options.report_filename, summary.FullReport());
}
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);
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
assert(solver.info() == Eigen::Success);
return summary;
}
// 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);
// Use the factors to solve the linear system
start_time = std::chrono::high_resolution_clock::now();
// solver.compute(Kg);
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();
summary.solve_time_in_ms = delta_time;
if (options.verbose)
std::cout << "System was solved in " << delta_time << " ms.\n" << std::endl;
assert(solver.info() == Eigen::Success);
if (solver.info() != Eigen::Success) {
std::cout << solver.lastErrorMessage() << std::endl;
summary.converged = false;
return summary;
}
}
if (options.save_elemental_forces) {
std::cout << "Writing to:" + options.elemental_forces_filename << std::endl;
csv.write(options.elemental_forces_filename, summary.element_forces,
options.csv_precision, options.csv_delimiter);
}
// convert to dense matrix
Eigen::VectorXd disp(dispSparse);
summary.converged = true;
// convert from Eigen vector to std vector
std::vector<std::vector<double>> disp_vec(structure.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);
}
summary.nodal_displacements = disp_vec;
return summary;
};
// [calculate nodal forces
start_time = std::chrono::high_resolution_clock::now();
Kg_local = Kg_local.topLeftCorner(dofs_per_elem * structure.nodes.size(),
dofs_per_elem * structure.nodes.size());
dispSparse = dispSparse.topRows(dofs_per_elem * structure.nodes.size());
SparseMat nodal_forces_sparse = Kg_local * dispSparse;
Eigen::VectorXd nodal_forces_dense(nodal_forces_sparse);
std::vector<std::vector<double>> nodal_forces_vec(structure.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);
}
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();
//]
// [ 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<std::chrono::milliseconds>(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) {
std::cout << "Writing to:" + options.nodal_displacements_filename << std::endl;
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_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<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();
summary.total_time_in_ms = delta_time;
if (options.save_report) {
writeStringToTxt(options.report_filename, summary.FullReport());
}
if (options.verbose)
std::cout << summary.FullReport();
// Compute per element forces
summary.element_forces.resize(structure.elems.size(), std::vector<double>(2 * DOF::NUM_DOFS));
for (size_t elemIndex = 0; elemIndex < structure.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 = structure.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 = structure.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);
}
}
if (options.save_elemental_forces) {
std::cout << "Writing to:" + options.elemental_forces_filename << std::endl;
csv.write(options.elemental_forces_filename,
summary.element_forces,
options.csv_precision,
options.csv_delimiter);
}
summary.converged = true;
return summary;
}
Summary ThreedBeamFEA::solve(const BeamStructure &structure,
const std::vector<BC> &BCs,
const std::vector<Force> &forces,
const std::vector<Tie> &ties,
const std::vector<Equation> &equations,
const Options &options)
{
initialize(structure);
return solve(BCs, forces, ties, equations, options);
}
void ThreedBeamFEA::initialize(const BeamStructure &structure)
{
this->structure = structure;
long long timeToAssemble;
const size_t KgConservativeSize = fea::NUM_DOFS * structure.nodes.size();
Kg.resize(KgConservativeSize, KgConservativeSize);
assembleStiffnessMatrix(structure, timeToAssemble, perElemKlocalAelem, Kg);
initialized = true;
}
bool ThreedBeamFEA::wasInitialized() const
{
return initialized;
}
} // namespace fea