From e7c70f635baaf16846fe48f46226266030f06981 Mon Sep 17 00:00:00 2001
From: Louis Denis <louis.denis@student.uliege.be>
Date: Fri, 13 May 2022 20:48:36 +0200
Subject: [PATCH] final cleaning FEM code + debug solver BEM
---
srcs/BEM/solverBEM.cpp | 2 +-
srcs/FEM/functionsFEM.cpp | 162 ++++++++++++++++----------------------
srcs/FEM/functionsFEM.hpp | 106 +++++++++++--------------
srcs/FEM/solverFEM.cpp | 37 +++++++++
4 files changed, 151 insertions(+), 156 deletions(-)
diff --git a/srcs/BEM/solverBEM.cpp b/srcs/BEM/solverBEM.cpp
index 6fa578e..afe6953 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 4f69cc0..4680d1d 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 1e9df82..6f5d81a 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 b27c1c8..02b5caa 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
--
GitLab