diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp index 6fa578ec232001e86458493880a74d7cafe4cf04..afe69531d6cce4a57ba613d755bfa44ad9cf7f04 100644 --- a/srcs/BEM/solverBEM.cpp +++ b/srcs/BEM/solverBEM.cpp @@ -86,7 +86,7 @@ void solverBEM(std::map<int, double> &electrostaticPressure, std::map<int,std::p exit(0); } - #pragma omp parallel for + //#pragma omp parallel for for(int i = 1; i <= nbDomains; i++) { singleDomainSolverBEM(electrostaticPressure, boundaryDisplacementMap, nbViews, postProcessing, iteration, i); diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp index 4f69cc031e1fd0bbfa33a3c14b875264433cc55d..4680d1d8f3513273a5fc98d4d20f43b295f78973 100644 --- a/srcs/FEM/functionsFEM.cpp +++ b/srcs/FEM/functionsFEM.cpp @@ -50,7 +50,7 @@ bool getNonLinearParameter() // corresponding to each nodeTag, looking at entity tags "tags" of dimension "dim" corresponding to the FEM domain. void RetrieveFEMNodeTagsAndCoords(std::vector<std::size_t> & FEM_nodeTags, std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, - int & dim, std::vector<int> & tags) + const int & dim, const std::vector<int> & tags) { //loop over the entities of the FEM domain to retrieve all nodeTags by taking care not to take some nodes twice //using sort and set_union methods + retrieving the node coordinates and storing them in a map @@ -90,7 +90,8 @@ void RetrieveFEMNodeTagsAndCoords(std::vector<std::size_t> & FEM_nodeTags, // associates an initialized elementData structure to each element tag through the "FEM_kinematicsMap" map, // using the "FEM_nodeCoordsMap" map gathering the coordinates of each node in the FEM domain. int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags) + std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, + const int & dim, const std::vector<int> & tags) { std::vector<int> elementTypes; std::vector<std::vector<std::size_t>> elementTags; @@ -179,7 +180,7 @@ int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, ele // updates the kinematics matrices of one element "element" based on its updated nodal positions and rotation angle. // if "computeFl" is set as true, also computes the local nodal forces vector "fl = Kl*pl", only OK if // the local stiffness matrix Kl has already been initialized, using "fillElementalKandF". -void updateKinMatrices(elementData* element, bool computeFl) +void updateKinMatrices(elementData* element, const bool computeFl) { // for the different formulas, please refer to the report size_t numNodes = element->nodes.size(); @@ -285,11 +286,11 @@ Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu) // "nodeIndexMap" mapping each nodeTag to its global index. // filling the local stiffness matrices stored in the elementData structure contained in "FEM_kinematicsMap". // numerical integration performed using Gauss integration based on "integration_rule". -void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, - const std::string & integration_rule, - int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, - std::vector<double> & full_f_vector, - std::map<int, elementData> & FEM_kinematicsMap) +void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double & rho, const double & bx, + const double & by, const std::string & integration_rule, + const int & dim, const std::vector<int> & tags, std::map<int, int> & nodeIndexMap, + std::vector<double> & full_f_vector, + std::map<int, elementData> & FEM_kinematicsMap) { std::vector<int> elementTypes; std::vector<std::vector<std::size_t>> elementTags; @@ -419,11 +420,8 @@ void computeBmatrix(const int &gaussIndex, const int &numNodes, const std::vecto // assembles the local "f_vector" into the full f vector, "elIndex" is the element index with respect to "nodeTags" // having "numNodes" nodes, "nodeTags" gathering all the nodeTags of the current element type in the current entity. // "nodeIndexMap" associates a global index to each nodetag. -void assembleFlocal(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int &numNodes, std::map<int, int> & nodeIndexMap, - std::vector<double> & full_f_vector, - std::vector<double> & f_vector) +void assembleFlocal(const int &elIndex, const std::vector<std::size_t> & nodeTags, const int &numNodes, + std::map<int, int> & nodeIndexMap, std::vector<double> & full_f_vector, const std::vector<double> & f_vector) { int first_node_index; int first_node; @@ -444,8 +442,7 @@ void assembleFlocal(const int &elIndex, // "BCkeys" gathers the keys associated to boundary conditions in the .geo file (using SetNumber(...)). // "groups" associates a dimension and a tag to each physical group. void includeNeumannBC(const std::vector<std::string> & BCkeys, const std::string & integration_rule, - std::map<int, int> & nodeIndexMap, - std::map<std::string, std::pair<int, int>> & groups, + std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::vector<double> & full_f_vector) { double tx = 0; // surface traction along x axis @@ -670,11 +667,9 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p // corresponding to dirichlet boundary conditions have been removed. A removed index will correspond to "new_DOFindices"= -1. // dir_BC[index] will contain the dirichlet boundary condition prescribed to the node whose global index is retrieved using // "nodeIndexMap". -void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys, - std::map<int, int> & nodeIndexMap, - std::map<std::string, std::pair<int, int>> & groups, - std::vector<int> & new_DOFindices, - std::vector<double> & dir_BC, int & number_removed_lines) +void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys, std::map<int, int> & nodeIndexMap, + std::map<std::string, std::pair<int, int>> & groups, + std::vector<int> & new_DOFindices, std::vector<double> & dir_BC, int & number_removed_lines) { for (std::size_t j = 0; j < new_DOFindices.size(); ++j){ @@ -779,8 +774,8 @@ void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys, // updating the kinematics of the "FEM_elementTags" contained in the "FEM_kinematicsMap", // based on the "full_d_vector" gathering the nodal displacements associated to the global // node indices contained in "nodeIndexMap". -void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector, - std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags) +void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, const std::vector<double> & full_d_vector, + std::map<int, int> & nodeIndexMap, const std::vector<std::size_t> & FEM_elementTags) { size_t numNodes; size_t tmp_nodeTag; @@ -871,7 +866,7 @@ void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vecto // "FEM_elementTags" contained in the "FEM_kinematicsMap", // each node Tag corresponds to a global index given by "nodeIndexMap". void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags) + std::map<int, int> & nodeIndexMap, const std::vector<std::size_t> & FEM_elementTags) { for(size_t i = 0; i < final_nodal_forces.size(); i++) { @@ -901,8 +896,8 @@ void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map // the nodes are renumbered using "new_DOFindices". // each node Tag corresponds to a global index given by "nodeIndexMap". void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, - std::vector<std::size_t> & FEM_elementTags) + std::map<int, int> & nodeIndexMap, const std::vector<int> & new_DOFindices, + const std::vector<std::size_t> & FEM_elementTags) { #pragma omp parallel for for(size_t i = 0; i < FEM_elementTags.size(); i++) @@ -934,8 +929,8 @@ void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elemen // "new_DOFindices" maps the global initial node index (from "nodeIndexMap") to the index // obtained without considering the fixed DOFs. void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, - std::vector<std::size_t> & FEM_elementTags) + std::map<int, int> & nodeIndexMap, const std::vector<int> & new_DOFindices, + const std::vector<std::size_t> & FEM_elementTags) { // to allow efficient parallelization, one triplet list is created for every thread, then the different lists are assembled // to further increase efficiency, space is reserved for each list @@ -1015,8 +1010,8 @@ void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std:: // updating the "full_d_vector" based on the "global_p_vector", using the global and reduced index mapping: // "new_DOFindices". -void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector, - std::vector<int> & new_DOFindices) +void updateNodalDispVector(std::vector<double> & full_d_vector, const std::vector<double> & global_p_vector, + const std::vector<int> & new_DOFindices) { for(size_t i = 0; i < full_d_vector.size(); i++) { @@ -1083,7 +1078,7 @@ void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boun // based on the "final_nodal_forces" and the "nodeIndexMap" to retrieve the global index associated // to one given node tag. void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups, - std::map<int, int> & nodeIndexMap, std::vector<double> & final_nodal_forces, + std::map<int, int> & nodeIndexMap, const std::vector<double> & final_nodal_forces, const std::vector<std::string> & BCkeys) { double Rx; @@ -1154,7 +1149,7 @@ void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups, // plotting the nodal displacement values associated to "FEM_nodeTags" in the "names[0]" model, // based on the "full_d_vector". void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector, - std::vector<std::string> & names) + const std::vector<std::string> & names) { int viewtagUx = gmsh::view::add("ux"); int viewtagUy = gmsh::view::add("uy"); @@ -1178,9 +1173,9 @@ void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<d // based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", looping over the // "nbelems" elements contained in "FEM_kinematicsMap". // returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag". -void computeStressStrainElNodal(int & dim, std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix, +void computeStressStrainElNodal(const int & dim, const std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap, - int & nbelems, double & max_VM_stress, int & max_VM_nodeTag) + const int & nbelems, double & max_VM_stress, int & max_VM_nodeTag) { std::vector<std::string> names; gmsh::model::list(names); @@ -1369,8 +1364,8 @@ void computeStressStrainElNodal(int & dim, std::vector<int> & tags, const Eigen: // computing the work done by external forces, as the nodal forces are energy consistent, it // is sufficient to compute the "externalWork" as the dot product of the nodal displacements // "full_d_vector" and the nodal forces "final_nodal_forces". -void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector, - std::vector<double> & final_nodal_forces) +void computeWorkDoneByExternalForces(double &externalWork, const std::vector<double> & full_d_vector, + const std::vector<double> & final_nodal_forces) { for(size_t i = 0; i < full_d_vector.size(); i++) { @@ -1382,7 +1377,7 @@ void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & // entity tags "tags"), using the elementData stored in "FEM_kinematicsMap" gathering the informations for all // FEM elements. The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains. // Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule". -void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & tags, +void computeStrainEnergy(double & strainEnergy, const int & dim, const std::vector<int> & tags, std::map<int, elementData> & FEM_kinematicsMap, const Eigen::Matrix<double, 3, 3> &H_matrix, const std::string & integration_rule) { @@ -1484,7 +1479,7 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta } //displays the time of the desired code section "function_name", if "display_time" is set as true. -void displayTime(double &start, double &end, std::string &function_name, bool &display_time) +void displayTime(const double &start, const double &end, const std::string &function_name, const bool &display_time) { if(display_time) { @@ -1513,17 +1508,26 @@ void displayTime(double &start, double &end, std::string &function_name, bool &d // node indexation defined by "nodeIndexMap". "H_matrix" is the plane stress Hooke's matrix, "rho" the density of the // material, "bx" the volumic force along x, "by" the volumic force along y. // Numerical integration using gauss integration "integration_rule". -void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, - const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, - Eigen::SparseMatrix<double> & full_K_matrix, std::vector<double> & full_f_vector) +void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> & H_matrix, const double & rho, const double & bx, const double & by, + const std::string & integration_rule, const int & dim, const std::vector<int> & tags, + std::map<int, int> & nodeIndexMap, Eigen::SparseMatrix<double> & full_K_matrix, + std::vector<double> & full_f_vector) { - // the full K matrix will first be assembled as a triplet list (vector) (to increase efficiency), - // before initializing the full sparse matrix in one single time; + // to allow efficient parallelization, one triplet list (as a vector) is created for every thread, then the different lists are assembled + // to further increase efficiency, space is reserved for each list + int estimated_number_of_entries = (int) full_f_vector.size()/2 *8*8; // accurate in the case of Q4 elements, otherwise it provides a good estimation + std::vector<Eigen::Triplet <double>> K_tripletList; + K_tripletList.reserve(estimated_number_of_entries); + int max_threads = omp_get_max_threads(); - // one triplet vector declared for every thread (for the K triplet list) (will be gathered at the end of the method) - std::vector<std::vector<Eigen::Triplet <double>>> my_k_list(max_threads); - // one pair vector declared for every thread (for the f vector) (will be gathered at the end of the method) - std::vector<std::vector<std::pair<int,double>>> thread_f_vector(max_threads); + int estimated_number_of_entries_per_thread = (int) estimated_number_of_entries / max_threads; + std::vector<std::vector<Eigen::Triplet <double>>> thread_K_list(max_threads); + for(int i = 0; i < max_threads; i++) + { + thread_K_list[i].reserve(estimated_number_of_entries_per_thread); + } + + std::vector<Eigen::Triplet <double>>* my_K_list; std::vector<int> elementTypes; std::vector<std::vector<std::size_t>> elementTags; @@ -1581,7 +1585,7 @@ void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double #pragma omp parallel for for (int el = 0; el < elementTags[i].size(); el++){ #else - #pragma omp parallel for + #pragma omp parallel for private (my_K_list) for (std::size_t el = 0; el < elementTags[i].size(); el++) { #endif @@ -1591,6 +1595,9 @@ void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double K_matrix.setZero(); for (int j = 0; j < 2*numNodes; ++j)f_vector[j] = 0.; + int threadNum = omp_get_thread_num(); + my_K_list = & thread_K_list[threadNum]; + // looping over the Gauss points for (std::size_t gaussIndex = 0; gaussIndex < weights.size(); ++gaussIndex){ // converting jacobian to Eigen format @@ -1623,40 +1630,30 @@ void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double } /*------------- ASSEMBLY PROCESS ------------*/ // USE LIST OF TRIPLETS TO PUSH ONLY ONCE -> more rapid and easier to parallelize - assembleKtriplet(el, nodeTags[i], numNodes, nodeIndexMap, my_k_list[omp_get_thread_num()], K_matrix); - assembleFthread(el, nodeTags[i], numNodes, nodeIndexMap, thread_f_vector[omp_get_thread_num()], f_vector); + assembleKtriplet(el, nodeTags[i], numNodes, nodeIndexMap, *my_K_list, K_matrix); + assembleFlocal(el, nodeTags[i], numNodes, nodeIndexMap, full_f_vector, f_vector); } } } // gathering all triplet vectors in one single vector - std::vector<Eigen::Triplet <double>> k_tripletList = my_k_list[0]; + K_tripletList = thread_K_list[0]; if(max_threads > 1) { for (int i = 1; i < max_threads; i++) { - k_tripletList.insert(k_tripletList.end(), my_k_list[i].begin(), my_k_list[i].end()); + K_tripletList.insert(K_tripletList.end(), thread_K_list[i].begin(), thread_K_list[i].end()); } } // assembling the full sparse K matrix - full_K_matrix.setFromTriplets(k_tripletList.begin(), k_tripletList.end()); - - // re-assembling the full f vector based on the separated thread vectors - for(int i = 0; i < max_threads; i++){ - std::vector<std::pair<int,double>> &tmp_vector = thread_f_vector[i]; - for(auto &p : tmp_vector) { - full_f_vector[p.first] += p.second; - } - } + full_K_matrix.setFromTriplets(K_tripletList.begin(), K_tripletList.end()); } // assembles the elemental "K_matrix" into a vector of triplets "k_tripletList" using // the node tags provided by "nodeTags" vector (index of the element is "elIndex"). // The element has "numNodes" nodes. "nodeIndexMap" provides the global nodal index corresponding to a given node tag. -void assembleKtriplet(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int numNodes, std::map<int, int> & nodeIndexMap, - std::vector<Eigen::Triplet <double>> & k_tripletList, - Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix) +void assembleKtriplet(const int &elIndex, const std::vector<std::size_t> & nodeTags, const int numNodes, + std::map<int, int> & nodeIndexMap, std::vector<Eigen::Triplet <double>> & k_tripletList, + const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix) { //we go through the lower half of the elemental K_matrix before assembling in the triplet vector for (int j = 0; j < 2*numNodes; ++j){ @@ -1676,35 +1673,12 @@ void assembleKtriplet(const int &elIndex, } } -// assembles the elemental "f_vector" in a vector of pairs called "thread_f_vector" particular to the parallel thread. -// the element index in "nodeTags" is "elIndex". The element has "numNodes" nodes and their global nodal indices are given -// by the "nodeIndexMap". -void assembleFthread(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int &numNodes, std::map<int, int> & nodeIndexMap, - std::vector<std::pair<int,double>> & thread_f_vector, - std::vector<double> & f_vector) -{ - // we simply go through the elemental f vector - int first_node_index; - int first_node; - int my_row; - for (int j = 0; j < 2*numNodes; ++j) - { - first_node_index = (int) j/2; - first_node = nodeTags[numNodes*elIndex+first_node_index]; - my_row = nodeIndexMap[first_node] + j % 2; - thread_f_vector.push_back(std::pair<int, double>(my_row, f_vector[j])); - } -} - // the "full_f_vector" and "full_K_matrix" are reduced in "reduced_f_vector" and "reduced_K_matrix" respectively, // taking into account the dirichlet boundary conditions contained in "dir_BC". // "new_DOFindices" gathers the new indices of the DOFs when taking the boundary conditions into account. void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::SparseMatrix<double> & full_K_matrix, - const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC, - Eigen::SparseMatrix<double> & reduced_K_matrix, - std::vector<double> & reduced_f_vector) + const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC, + Eigen::SparseMatrix<double> & reduced_K_matrix, std::vector<double> & reduced_f_vector) { //the correction term on the f vector comes from the "Finite Element Method" theoretical course (JP Ponthot, ULiege), Chap2 Slide 38/50 std::list<Eigen::Triplet <double>> reduced_k_tripletList; @@ -1751,9 +1725,9 @@ void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::Spars // displacement vector "full_d_vector" and the global nodal index map "nodeIndexMap". // "nbelems" is the number of 2D elements contained in the FEM domain. // returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag". -void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector, - const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, - std::map<int, int> & nodeIndexMap, int & nbelems, double & max_VM_stress, int & max_VM_nodeTag) +void computeStressStrainLinearWay(const int & dim, const std::vector<int> & tags, const std::vector<double> &full_d_vector, + const Eigen::Matrix<double, 3, 3> &H_matrix, const double &nu, const double &E, + std::map<int, int> & nodeIndexMap, const int &nbelems, double &max_VM_stress, int &max_VM_nodeTag) { std::vector<std::string> names; gmsh::model::list(names); @@ -1926,7 +1900,7 @@ void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std: // entity tags "tags"), using the nodal displacements in "full_d_vector". The nodes are indexed using the "nodeIndexMap". // The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains. // Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule". -void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<int> & tags, +void computeStrainEnergyLinearWay(double & strainEnergy, const int & dim, const std::vector<int> & tags, const std::vector<double> &full_d_vector, const Eigen::Matrix<double, 3, 3> &H_matrix, std::map<int, int> & nodeIndexMap, const std::string & integration_rule) { diff --git a/srcs/FEM/functionsFEM.hpp b/srcs/FEM/functionsFEM.hpp index 1e9df823933dcb29e35ac886c1e97ba80444201f..6f5d81ad887984f481b12384741282f927b2c3f9 100644 --- a/srcs/FEM/functionsFEM.hpp +++ b/srcs/FEM/functionsFEM.hpp @@ -73,19 +73,20 @@ bool getNonLinearParameter(); // corresponding to each nodeTag, looking at entity tags "tags" of dimension "dim" corresponding to the FEM domain. void RetrieveFEMNodeTagsAndCoords(std::vector<std::size_t> & FEM_nodeTags, std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, - int & dim, std::vector<int> & tags); + const int & dim, const std::vector<int> & tags); // fills the "FEM_elementTags" vector with all element tags of the elements of dimension "dim" (in practice: 2) // in the entities associated to the entity tags "tags" (corresponding to the FEM domain). // associates an initialized elementData structure to each element tag through the "FEM_kinematicsMap" map, // using the "FEM_nodeCoordsMap" map gathering the coordinates of each node in the FEM domain. int initKinematics(std::vector<std::size_t> & FEM_elementTags, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, int & dim, std::vector<int> & tags); + std::map<int,std::pair<double,double>> & FEM_nodeCoordsMap, + const int & dim, const std::vector<int> & tags); // updates the kinematics matrices of one element "element" based on its updated nodal positions and rotation angle. // if "computeFl" is set as true, also computes the local nodal forces vector "fl = Kl*pl", only OK if // the local stiffness matrix Kl has already been initialized, using "fillElementalKandF". -void updateKinMatrices(elementData* element, bool computeFl); +void updateKinMatrices(elementData* element, const bool computeFl); // gets physical properties of FEM_domain: Young's modulus "E", Poisson's ratio "nu", // density "rho", volumic force along x "bx", volumic force along y "by". @@ -99,11 +100,11 @@ Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu); // "nodeIndexMap" mapping each nodeTag to its global index. // filling the local stiffness matrices stored in the elementData structure contained in "FEM_kinematicsMap". // numerical integration performed using Gauss integration based on "integration_rule". -void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, - const std::string & integration_rule, - int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, - std::vector<double> & full_f_vector, - std::map<int, elementData> & FEM_kinematicsMap); +void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double & rho, const double & bx, + const double & by, const std::string & integration_rule, + const int & dim, const std::vector<int> & tags, std::map<int, int> & nodeIndexMap, + std::vector<double> & full_f_vector, + std::map<int, elementData> & FEM_kinematicsMap); // computing the strain displacement "B_matrix" at the "gaussIndex"-th gaussian point of the element // consisting of "numNodes" nodes, based on the "basisFunctionsGradient", "determinants" and the inverse of the jacobian @@ -116,11 +117,8 @@ void computeBmatrix(const int &gaussIndex, const int &numNodes, const std::vecto // assembles the local "f_vector" into the full f vector, "elIndex" is the element index with respect to "nodeTags" // having "numNodes" nodes, "nodeTags" gathering all the nodeTags of the current element type in the current entity. // "nodeIndexMap" associates a global index to each nodetag. -void assembleFlocal(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int &numNodes, std::map<int, int> & nodeIndexMap, - std::vector<double> & full_f_vector, - std::vector<double> & f_vector); +void assembleFlocal(const int &elIndex, const std::vector<std::size_t> & nodeTags, const int &numNodes, + std::map<int, int> & nodeIndexMap, std::vector<double> & full_f_vector, const std::vector<double> & f_vector); // fills the "full_f_vector" with the surfacic contribution of the neumann boundary conditions based on the // "nodeIndexMap" associating a global index to each node tag, the local integration of the surfacic forces is @@ -128,8 +126,7 @@ void assembleFlocal(const int &elIndex, // "BCkeys" gathers the keys associated to boundary conditions in the .geo file (using SetNumber(...)). // "groups" associates a dimension and a tag to each physical group. void includeNeumannBC(const std::vector<std::string> & BCkeys, const std::string & integration_rule, - std::map<int, int> & nodeIndexMap, - std::map<std::string, std::pair<int, int>> & groups, + std::map<int, int> & nodeIndexMap, std::map<std::string, std::pair<int, int>> & groups, std::vector<double> & full_f_vector); // for the coupled solver, fills the "full_f_vector" with the electrostatic pressure contained in the @@ -146,23 +143,21 @@ void applyElecPressure(const std::string & integration_rule, std::map<int,std::p // corresponding to dirichlet boundary conditions have been removed. A removed index will correspond to "new_DOFindices"= -1. // dir_BC[index] will contain the dirichlet boundary condition prescribed to the node whose global index is retrieved using // "nodeIndexMap". -void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys, - std::map<int, int> & nodeIndexMap, - std::map<std::string, std::pair<int, int>> & groups, - std::vector<int> & new_DOFindices, - std::vector<double> & dir_BC, int & number_removed_lines); +void fillDOFindicesDirBC(const std::vector<std::string> & BCkeys, std::map<int, int> & nodeIndexMap, + std::map<std::string, std::pair<int, int>> & groups, + std::vector<int> & new_DOFindices, std::vector<double> & dir_BC, int & number_removed_lines); // updating the kinematics of the "FEM_elementTags" contained in the "FEM_kinematicsMap", // based on the "full_d_vector" gathering the nodal displacements associated to the global // node indices contained in "nodeIndexMap". -void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, std::vector<double> & full_d_vector, - std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags); +void updateKinematics(std::map<int, elementData> & FEM_kinematicsMap, const std::vector<double> & full_d_vector, + std::map<int, int> & nodeIndexMap, const std::vector<std::size_t> & FEM_elementTags); // retrieving the global nodal forces based on local nodal forces of the // "FEM_elementTags" contained in the "FEM_kinematicsMap", // each node Tag corresponds to a global index given by "nodeIndexMap". void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<std::size_t> & FEM_elementTags); + std::map<int, int> & nodeIndexMap, const std::vector<std::size_t> & FEM_elementTags); // assembles the "global_f_vector" based on the local nodal forces vector of the // "FEM_elementTags" contained in the "FEM_kinematicsMap", @@ -170,21 +165,21 @@ void retrieveTotalNodalForces(std::vector<double> & final_nodal_forces, std::map // the nodes are renumbered using "new_DOFindices". // each node Tag corresponds to a global index given by "nodeIndexMap". void assembleGlobalF(std::vector<double> & global_f_vector, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, - std::vector<std::size_t> & FEM_elementTags); + std::map<int, int> & nodeIndexMap, const std::vector<int> & new_DOFindices, + const std::vector<std::size_t> & FEM_elementTags); // assembling the global tangent stiffness matrix based on the local elemental stiffness matrices contained // in the "FEM_kinematicsMap" associated to each "FEM_elementTags". // "new_DOFindices" maps the global initial node index (from "nodeIndexMap") to the index // obtained without considering the fixed DOFs. void assembleGlobalKg(Eigen::SparseMatrix<double> & global_tangent_matrix, std::map<int, elementData> & FEM_kinematicsMap, - std::map<int, int> & nodeIndexMap, std::vector<int> & new_DOFindices, - std::vector<std::size_t> & FEM_elementTags); + std::map<int, int> & nodeIndexMap, const std::vector<int> & new_DOFindices, + const std::vector<std::size_t> & FEM_elementTags); // updating the "full_d_vector" based on the "global_p_vector", using the global and reduced index mapping: // "new_DOFindices". -void updateNodalDispVector(std::vector<double> & full_d_vector, std::vector<double> & global_p_vector, - std::vector<int> & new_DOFindices); +void updateNodalDispVector(std::vector<double> & full_d_vector, const std::vector<double> & global_p_vector, + const std::vector<int> & new_DOFindices); // retrieving the nodal displacements on the BEM_FEM_boundary in order to communicate them to the // BEM solver in the two-way coupled solver. Using the "groups" to find the physical group associated @@ -199,38 +194,38 @@ void retrieveBoundaryDisplacements(std::map<int,std::pair<double,double>> & boun // based on the "final_nodal_forces" and the "nodeIndexMap" to retrieve the global index associated // to one given node tag. void computeReactionForces(std::map<std::string, std::pair<int, int>> & groups, - std::map<int, int> & nodeIndexMap, std::vector<double> & final_nodal_forces, + std::map<int, int> & nodeIndexMap, const std::vector<double> & final_nodal_forces, const std::vector<std::string> & BCkeys); // plotting the nodal displacement values associated to "FEM_nodeTags" in the "names[0]" model, // based on the "full_d_vector". void plotUxUy(const std::vector<std::size_t> & FEM_nodeTags, const std::vector<double> & full_d_vector, - std::vector<std::string> & names); + const std::vector<std::string> & names); // computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"), // based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", looping over the // "nbelems" elements contained in "FEM_kinematicsMap". // returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag". -void computeStressStrainElNodal(int & dim, std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix, +void computeStressStrainElNodal(const int & dim, const std::vector<int> & tags, const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, std::map<int, elementData> & FEM_kinematicsMap, - int & nbelems, double & max_VM_stress, int & max_VM_nodeTag); + const int & nbelems, double & max_VM_stress, int & max_VM_nodeTag); // computing the work done by external forces, as the nodal forces are energy consistent, it // is sufficient to compute the "externalWork" as the dot product of the nodal displacements // "full_d_vector" and the nodal forces "final_nodal_forces". -void computeWorkDoneByExternalForces(double &externalWork, std::vector<double> & full_d_vector, - std::vector<double> & final_nodal_forces); +void computeWorkDoneByExternalForces(double &externalWork, const std::vector<double> & full_d_vector, + const std::vector<double> & final_nodal_forces); // computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and // entity tags "tags"), using the elementData stored in "FEM_kinematicsMap" gathering the informations for all // FEM elements. The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains. // Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule". -void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & tags, +void computeStrainEnergy(double & strainEnergy, const int & dim, const std::vector<int> & tags, std::map<int, elementData> & FEM_kinematicsMap, const Eigen::Matrix<double, 3, 3> &H_matrix, const std::string & integration_rule); //displays the time of the desired code section "function_name", if "display_time" is set as true. -void displayTime(double &start, double &end, std::string &function_name, bool &display_time); +void displayTime(const double &start, const double &end, const std::string &function_name, const bool &display_time); /*-------------- FUNCTIONS SPECIFIC TO LINEAR FEM ------------------*/ @@ -239,49 +234,38 @@ void displayTime(double &start, double &end, std::string &function_name, bool &d // node indexation defined by "nodeIndexMap". "H_matrix" is the plane stress Hooke's matrix, "rho" the density of the // material, "bx" the volumic force along x, "by" the volumic force along y. // Numerical integration using gauss integration "integration_rule". -void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, - const std::string & integration_rule, int & dim, std::vector<int> & tags, std::map<int, int> & nodeIndexMap, - Eigen::SparseMatrix<double> & full_K_matrix, std::vector<double> & full_f_vector); +void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> & H_matrix, const double & rho, const double & bx, const double & by, + const std::string & integration_rule, const int & dim, const std::vector<int> & tags, + std::map<int, int> & nodeIndexMap, Eigen::SparseMatrix<double> & full_K_matrix, + std::vector<double> & full_f_vector); // assembles the elemental "K_matrix" into a vector of triplets "k_tripletList" using // the node tags provided by "nodeTags" vector (index of the element is "elIndex"). // The element has "numNodes" nodes. "nodeIndexMap" provides the global nodal index corresponding to a given node tag. -void assembleKtriplet(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int numNodes, std::map<int, int> & nodeIndexMap, - std::vector<Eigen::Triplet <double>> & k_tripletList, - Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix); - -// assembles the elemental "f_vector" in a vector of pairs called "thread_f_vector" particular to the parallel thread. -// the element index in "nodeTags" is "elIndex". The element has "numNodes" nodes and their global nodal indices are given -// by the "nodeIndexMap". -void assembleFthread(const int &elIndex, - const std::vector<std::size_t> & nodeTags, - const int &numNodes, std::map<int, int> & nodeIndexMap, - std::vector<std::pair<int,double>> & thread_f_vector, - std::vector<double> & f_vector); +void assembleKtriplet(const int &elIndex, const std::vector<std::size_t> & nodeTags, const int numNodes, + std::map<int, int> & nodeIndexMap, std::vector<Eigen::Triplet <double>> & k_tripletList, + const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & K_matrix); // the "full_f_vector" and "full_K_matrix" are reduced in "reduced_f_vector" and "reduced_K_matrix" respectively, // taking into account the dirichlet boundary conditions contained in "dir_BC". // "new_DOFindices" gathers the new indices of the DOFs when taking the boundary conditions into account. void fillReducedKf(const std::vector<double> & full_f_vector, const Eigen::SparseMatrix<double> & full_K_matrix, - const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC, - Eigen::SparseMatrix<double> & reduced_K_matrix, - std::vector<double> & reduced_f_vector); + const std::vector<int> & new_DOFindices, const std::vector<double> & dir_BC, + Eigen::SparseMatrix<double> & reduced_K_matrix, std::vector<double> & reduced_f_vector); // computing the elemental-nodal strains and stresses in the FEM domain (dimension "dim" and entity tags "tags"), // based on Hooke's matrix in plane stress "H_matrix", Young's modulus "E", Poisson's ratio "nu", and the full nodal // displacement vector "full_d_vector" and the global nodal index map "nodeIndexMap". // "nbelems" is the number of 2D elements contained in the FEM domain. // returning the "max_VM_stress" and the corresponding nodetag "max_VM_nodeTag". -void computeStressStrainLinearWay(int & dim, std::vector<int> & tags, const std::vector<double> &full_d_vector, - const Eigen::Matrix<double, 3, 3> &H_matrix, const double nu, const double E, - std::map<int, int> & nodeIndexMap, int & nbelems, double & max_VM_stress, int & max_VM_nodeTag); +void computeStressStrainLinearWay(const int & dim, const std::vector<int> & tags, const std::vector<double> &full_d_vector, + const Eigen::Matrix<double, 3, 3> &H_matrix, const double &nu, const double &E, + std::map<int, int> & nodeIndexMap, const int &nbelems, double &max_VM_stress, int &max_VM_nodeTag); // computing the total (internal) strain energy "strain_energy" in the whole FEM domain (dimension "dim" and // entity tags "tags"), using the nodal displacements in "full_d_vector". The nodes are indexed using the "nodeIndexMap". // The Hooke's matrix in plane stress "H_matrix" is used to compute the stresses from the strains. // Computed as a sum of integrals over the finite elements, using Gauss integration "integration_rule". -void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<int> & tags, +void computeStrainEnergyLinearWay(double & strainEnergy, const int & dim, const std::vector<int> & tags, const std::vector<double> &full_d_vector, const Eigen::Matrix<double, 3, 3> &H_matrix, std::map<int, int> & nodeIndexMap, const std::string & integration_rule); \ No newline at end of file diff --git a/srcs/FEM/solverFEM.cpp b/srcs/FEM/solverFEM.cpp index b27c1c8dd556c12bbd274d059ed2939915012ac4..02b5caa8ea32f7a60429190d44b414c34e365277 100644 --- a/srcs/FEM/solverFEM.cpp +++ b/srcs/FEM/solverFEM.cpp @@ -478,6 +478,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) /*--------Integration rule used in the code----------*/ std::string integration_rule = "CompositeGauss4"; + bool display_time = true; //displays the time performance of the different parts of the code + + double start_init = omp_get_wtime(); + /*--------get physical groups--------------*/ std::map<std::string, std::pair<int, int>> groups; gmsh::vectorpair dimTags; @@ -524,6 +528,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) /*-----------------H matrix involved in stress computation------------------*/ Eigen::Matrix<double, 3, 3> H_matrix = computeHmatrix(E, nu); + double end_init = omp_get_wtime(); + std::string step_name = "initialization"; + displayTime(start_init, end_init, step_name, display_time); + /*------------computing K matrix and f vector------------------------*/ // we first consider the K matrix and the volumic part of the f vector // K matrix declared as sparse matrix for efficiency. @@ -537,6 +545,11 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) // Looping over entities and compute full K and volumic F assembleKVolumicF(H_matrix, rho, bx, by, integration_rule, dim, tags, nodeIndexMap, full_K_matrix, full_f_vector); + double end_assembleK = omp_get_wtime(); + step_name = "K assembly"; + displayTime(end_init, end_assembleK, step_name, display_time); + + /*-----------------BOUNDARY CONDITIONS------------*/ /*-----------get group_names related to BC's----------------------*/ std::vector<std::string> BCkeys; @@ -550,6 +563,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) applyElecPressure(integration_rule, FEM_nodeCoordsMap, nodeIndexMap, groups, electrostaticPressure, full_f_vector); } + double end_neumannBC = omp_get_wtime(); + step_name = "neumann BC and elec pressure"; + displayTime(end_assembleK, end_neumannBC, step_name, display_time); + /*------------FOCUS ON DIRICHLET BC'S---------------*/ // new_DOFindices vector will contain the new indices of the K matrix after the rows/columns corresponding // to dirichlet boundary conditions have been removed @@ -565,6 +582,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) int number_removed_lines = 0; fillDOFindicesDirBC(BCkeys, nodeIndexMap, groups, new_DOFindices, dir_BC, number_removed_lines); + double end_dirichletBC = omp_get_wtime(); + step_name = "dirichlet BC"; + displayTime(end_neumannBC, end_dirichletBC, step_name, display_time); + /*------------FILLING THE REDUCED K MATRIX AND F VECTOR--------------*/ //the correction term on the f vector comes from the "Finite Element Method" theoretical course (JP Ponthot, ULiege), Chap2 Slide 38/50 Eigen::SparseMatrix<double> reduced_K_matrix(2*FEM_nodeTags.size() - number_removed_lines,2*FEM_nodeTags.size() - number_removed_lines); @@ -572,6 +593,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) fillReducedKf(full_f_vector, full_K_matrix, new_DOFindices, dir_BC, reduced_K_matrix, reduced_f_vector); + double end_reducedFilling = omp_get_wtime(); + step_name = "filling reduced K and f"; + displayTime(end_dirichletBC, end_reducedFilling, step_name, display_time); + /*--------------SOLVING THE REDUCED SYSTEM--------------------*/ Eigen::VectorXd final_reduced_d(2*FEM_nodeTags.size() - number_removed_lines); @@ -586,6 +611,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) final_reduced_d = solver.solve(final_reduced_f); + double end_solveSystem = omp_get_wtime(); + step_name = "solving the linear system"; + displayTime(end_reducedFilling, end_solveSystem, step_name, display_time); + /*-----------RECONSTRUCTING FULL NODAL VALUES--------------*/ std::vector<double> full_d_vector(2*FEM_nodeTags.size()); for (std::size_t j = 0; j < full_d_vector.size(); ++j){ @@ -611,6 +640,10 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) nodal_forces[i] = nodalForces(i); } + double end_reconstructing = omp_get_wtime(); + step_name = "reconstructing full nodal disp and forces"; + displayTime(end_solveSystem, end_reconstructing, step_name, display_time); + /*--------------PLOTTING AND SAVING NODAL VALUES------------------*/ std::vector<std::string> names; gmsh::model::list(names); @@ -716,4 +749,8 @@ void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) std::cout << "-----------------------------------------\n"; std::cout << "Total potential energy: " << strainEnergy - externalWork << " [J/m].\n"; + + double end_postPro = omp_get_wtime(); + step_name = "post processing"; + displayTime(end_reconstructing, end_postPro, step_name, display_time); } \ No newline at end of file