diff --git a/CMakeLists.txt b/CMakeLists.txt index a10ac15..e73c561 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/include/containers.h b/include/containers.h index 7f5d68a..a7c671f 100644 --- a/include/containers.h +++ b/include/containers.h @@ -386,19 +386,19 @@ struct Elem { * @brief Contains a node list, element list, and the properties of each * element. */ -struct Job { - std::vector nodes; /** elems; /** nodes; /** elems; /** - props; /** props; /** &_nodes, const std::vector _elems) - : nodes(_nodes) { - unsigned int num_elems = _elems.size(); - elems.reserve(num_elems); - props.reserve(num_elems); + BeamStructure(const std::vector &_nodes, const std::vector _elems) : nodes(_nodes) + { + unsigned int num_elems = _elems.size(); + elems.reserve(num_elems); + props.reserve(num_elems); - for (unsigned int i = 0; i < num_elems; i++) { - elems.push_back(_elems[i].node_numbers); - props.push_back(_elems[i].props); - } - }; + for (unsigned int i = 0; i < num_elems; i++) { + elems.push_back(_elems[i].node_numbers); + props.push_back(_elems[i].props); + } + }; }; } // namespace fea diff --git a/include/setup.h b/include/setup.h index dbd59fb..ef33447 100644 --- a/include/setup.h +++ b/include/setup.h @@ -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 diff --git a/include/threed_beam_fea.h b/include/threed_beam_fea.h index f71c47a..8e542fa 100644 --- a/include/threed_beam_fea.h +++ b/include/threed_beam_fea.h @@ -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 &ties); + void operator()(SparseMat &Kg, const BeamStructure &job, const std::vector &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 getPerElemKlocalAelem() const; + //FIX: This allows the client to change the data member. I think there is a proble with eigen.. + std::vector &getPerElemKlocalAelem(); -private: + private: LocalMatrix Kelem; /** &perElemKlocalAelem, + SparseMat &Kg); /** * @brief Loads the boundary conditions into the global stiffness matrix and * force vector. @@ -250,7 +258,9 @@ computeTieForces(const std::vector &ties, */ void loadForces(SparseMat &force_vec, const std::vector &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 &forces); * @return Summary `fea::Summary`. Summary containing the results of the * analysis. */ -Summary solve(const Job &job, const std::vector &BCs, - const std::vector &forces, const std::vector &ties, - const std::vector &equations, const Options &options); + Summary solve(const std::vector &BCs, + const std::vector &forces, + const std::vector &ties, + const std::vector &equations, + const Options &options); + Summary solve(const BeamStructure &structure, + const std::vector &BCs, + const std::vector &forces, + const std::vector &ties, + const std::vector &equations, + const Options &options); + + void initialize(const BeamStructure &structure); + std::vector perElemKlocalAelem; + SparseMat Kg; + bool initialized{false}; + BeamStructure structure; + +public: + bool wasInitialized() const; +}; } // namespace fea #endif // THREED_BEAM_FEA_H diff --git a/src/setup.cpp b/src/setup.cpp index 8cb65b0..f275968 100644 --- a/src/setup.cpp +++ b/src/setup.cpp @@ -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 nodes = createNodeVecFromJSON(config_doc); std::vector elems = createElemVecFromJSON(config_doc); - return Job(nodes, elems); + return BeamStructure(nodes, elems); } Options createOptionsFromJSON(const rapidjson::Document &config_doc) { diff --git a/src/threed_beam_fea.cpp b/src/threed_beam_fea.cpp index 1888f0b..49b7a5e 100644 --- a/src/threed_beam_fea.cpp +++ b/src/threed_beam_fea.cpp @@ -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 GlobalStiffAssembler::getPerElemKlocalAelem() const { - return perElemKlocalAelem; -}; +std::vector &GlobalStiffAssembler::getPerElemKlocalAelem() +{ + return perElemKlocalAelem; +} -void GlobalStiffAssembler::operator()(SparseMat &Kg, const Job &job, - const std::vector &ties) { - int nn1, nn2; - unsigned int row, col; - const unsigned int dofs_per_elem = DOF::NUM_DOFS; +void GlobalStiffAssembler::operator()(SparseMat &Kg, + const BeamStructure &job, + const std::vector &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> 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> 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(dofs_per_elem * nn1 + row, - dofs_per_elem * nn1 + col, - it.value())); - } - // top right - else { - triplets.push_back(Eigen::Triplet( - dofs_per_elem * nn1 + row, dofs_per_elem * (nn2 - 1) + col, - it.value())); // I: Giati nn2-1 kai oxi nn2? - } - } else { - // bottom left - if (col < 6) { - triplets.push_back(Eigen::Triplet( - dofs_per_elem * (nn2 - 1) + row, dofs_per_elem * nn1 + col, - it.value())); // I: Giati nn2-1 kai oxi nn2? Epeidi ta nn2 ston - // Kelem exoun idi offset 6 - } - // bottom right - else { - triplets.push_back(Eigen::Triplet( - dofs_per_elem * (nn2 - 1) + row, - dofs_per_elem * (nn2 - 1) + col, it.value())); - } + // check position in local matrix and update corresponding global + // position + if (row < 6) { + // top left + if (col < 6) { + triplets.push_back(Eigen::Triplet(dofs_per_elem * ni1 + row, + dofs_per_elem * ni1 + col, + it.value())); + } + // top right + else { + triplets.push_back( + Eigen::Triplet(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( + 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(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 &BCs, unsigned int num_nodes) { @@ -268,7 +279,7 @@ void loadBCs(SparseMat &Kg, SparseMat &force_vec, const std::vector &BCs, force_vec.insert(global_add_idx + i, 0) = BCs[i].value; } } -}; +} void loadEquations(SparseMat &Kg, const std::vector &equations, unsigned int num_nodes, unsigned int num_bcs) { @@ -285,7 +296,7 @@ void loadEquations(SparseMat &Kg, const std::vector &equations, Kg.insert(col_idx, row_idx) = equations[i].terms[j].coefficient; } } -}; +} void loadTies(std::vector> &triplets, const std::vector &ties) { @@ -355,61 +366,69 @@ void loadForces(SparseMat &force_vec, const std::vector &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 &BCs, - const std::vector &forces, const std::vector &ties, - const std::vector &equations, const Options &options) { - auto initial_start_time = std::chrono::high_resolution_clock::now(); +void assembleStiffnessMatrix(const BeamStructure &structure, + long long &timeToAssemble, + std::vector &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(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 &BCs, + const std::vector &forces, + const std::vector &ties, + const std::vector &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( - 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 &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 &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(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(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 kg_invCheck(Kg); - assert(kg_invCheck.isInvertible()); +// Eigen::VectorXd kgCol(Kg.col(7)); +// Eigen::FullPivLU kg_invCheck(Kg); +// assert(kg_invCheck.isInvertible()); #endif - end_time = std::chrono::high_resolution_clock::now(); - - delta_time = std::chrono::duration_cast(end_time - - start_time) - .count(); - summary.factorization_time_in_ms = delta_time; - - 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(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> disp_vec(job.nodes.size(), - std::vector(dofs_per_elem)); - for (size_t i = 0; i < disp_vec.size(); ++i) { - for (unsigned int j = 0; j < dofs_per_elem; ++j) - // round all values close to 0.0 - disp_vec[i][j] = std::abs(disp(dofs_per_elem * i + j)) < options.epsilon - ? 0.0 - : disp(dofs_per_elem * i + j); - } - summary.nodal_displacements = disp_vec; - - // [calculate nodal forces - start_time = std::chrono::high_resolution_clock::now(); - - 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> nodal_forces_vec( - job.nodes.size(), std::vector(dofs_per_elem)); - for (size_t i = 0; i < nodal_forces_vec.size(); ++i) { - for (unsigned int j = 0; j < dofs_per_elem; ++j) - // round all values close to 0.0 - nodal_forces_vec[i][j] = - std::abs(nodal_forces_dense(dofs_per_elem * i + j)) < options.epsilon - ? 0.0 - : nodal_forces_dense(dofs_per_elem * i + j); - } - summary.nodal_forces = nodal_forces_vec; - - end_time = std::chrono::high_resolution_clock::now(); - - summary.nodal_forces_solve_time_in_ms = - std::chrono::duration_cast(end_time - - start_time) - .count(); - //] - - // [ calculate forces associated with ties - if (ties.size() > 0) { - start_time = std::chrono::high_resolution_clock::now(); - summary.tie_forces = computeTieForces(ties, disp_vec); end_time = std::chrono::high_resolution_clock::now(); - summary.tie_forces_solve_time_in_ms = - std::chrono::duration_cast(end_time - - start_time) - .count(); - } - // ] - // [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(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(end_time - - start_time) - .count(); - // ] - - auto final_end_time = std::chrono::high_resolution_clock::now(); - - delta_time = std::chrono::duration_cast( - final_end_time - initial_start_time) - .count(); - summary.total_time_in_ms = delta_time; - - 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(2 * DOF::NUM_DOFS)); - std::vector perElemKlocalAelem = - assembleK3D.getPerElemKlocalAelem(); - - for (size_t elemIndex = 0; elemIndex < job.elems.size(); elemIndex++) { - // Assemble the displacements of the nodes of the beam - Eigen::Matrix elemDisps; - // Displacements of the first node of the beam - const size_t nodeIndex0 = job.elems[elemIndex](0); - elemDisps(0) = summary.nodal_displacements[nodeIndex0][0]; - elemDisps(1) = summary.nodal_displacements[nodeIndex0][1]; - elemDisps(2) = summary.nodal_displacements[nodeIndex0][2]; - elemDisps(3) = summary.nodal_displacements[nodeIndex0][3]; - elemDisps(4) = summary.nodal_displacements[nodeIndex0][4]; - elemDisps(5) = summary.nodal_displacements[nodeIndex0][5]; - // Displacemen ts of the second node of the beam - const size_t nodeIndex1 = job.elems[elemIndex](1); - elemDisps(6) = summary.nodal_displacements[nodeIndex1][0]; - elemDisps(7) = summary.nodal_displacements[nodeIndex1][1]; - elemDisps(8) = summary.nodal_displacements[nodeIndex1][2]; - elemDisps(9) = summary.nodal_displacements[nodeIndex1][3]; - elemDisps(10) = summary.nodal_displacements[nodeIndex1][4]; - elemDisps(11) = summary.nodal_displacements[nodeIndex1][5]; - - Eigen::Matrix elemForces = - perElemKlocalAelem[elemIndex] * elemDisps; - for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) { - summary.element_forces[elemIndex][static_cast(dofIndex)] = -elemForces(dofIndex); + 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(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(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> disp_vec(structure.nodes.size(), + std::vector(dofs_per_elem)); + for (size_t i = 0; i < disp_vec.size(); ++i) { + for (unsigned int j = 0; j < dofs_per_elem; ++j) + // round all values close to 0.0 + disp_vec[i][j] = std::abs(disp(dofs_per_elem * i + j)) < options.epsilon + ? 0.0 + : disp(dofs_per_elem * i + j); + } + summary.nodal_displacements = disp_vec; - 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> nodal_forces_vec(structure.nodes.size(), + std::vector(dofs_per_elem)); + for (size_t i = 0; i < nodal_forces_vec.size(); ++i) { + for (unsigned int j = 0; j < dofs_per_elem; ++j) + // round all values close to 0.0 + nodal_forces_vec[i][j] = std::abs(nodal_forces_dense(dofs_per_elem * i + j)) + < options.epsilon + ? 0.0 + : nodal_forces_dense(dofs_per_elem * i + j); + } + summary.nodal_forces = nodal_forces_vec; + + end_time = std::chrono::high_resolution_clock::now(); + + summary.nodal_forces_solve_time_in_ms + = std::chrono::duration_cast(end_time - start_time).count(); + //] + + // [ calculate forces associated with ties + if (ties.size() > 0) { + start_time = std::chrono::high_resolution_clock::now(); + summary.tie_forces = computeTieForces(ties, disp_vec); + end_time = std::chrono::high_resolution_clock::now(); + summary.tie_forces_solve_time_in_ms + = std::chrono::duration_cast(end_time - start_time).count(); + } + // ] + + // [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(end_time - start_time).count(); + // ] + + auto final_end_time = std::chrono::high_resolution_clock::now(); + + delta_time = std::chrono::duration_cast(final_end_time + - initial_start_time) + .count(); + summary.total_time_in_ms = delta_time; + + 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(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 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 elemForces = perElemKlocalAelem[elemIndex] * elemDisps; + for (int dofIndex = 0; dofIndex < DOF::NUM_DOFS; dofIndex++) { + summary.element_forces[elemIndex][static_cast(dofIndex)] = -elemForces(dofIndex); + } + // meaning of sign = reference to first node: + // + = compression for axial, - traction + for (int dofIndex = DOF::NUM_DOFS; dofIndex < 2 * DOF::NUM_DOFS; dofIndex++) { + summary.element_forces[elemIndex][static_cast(dofIndex)] = +elemForces(dofIndex); + } + } + + 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 &BCs, + const std::vector &forces, + const std::vector &ties, + const std::vector &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