diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp index 4680d1d8f3513273a5fc98d4d20f43b295f78973..4e39d98fc945929feef1a900577489a4f010f6de 100644 --- a/srcs/FEM/functionsFEM.cpp +++ b/srcs/FEM/functionsFEM.cpp @@ -285,12 +285,9 @@ Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu) // with the local volumic forces computed based on the volumic force ("bx","by") and density "rho" with the help of // "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, - 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) + const double & by, 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; @@ -316,10 +313,17 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double elementName, elDim, order, numNodes, localNodeCoord, numPrimaryNodes); + + // adapting the composite integration rule as a function of the order of the elements + std::string adapting_integration_rule = "CompositeGauss14"; + if(order < 8) + { + adapting_integration_rule = "CompositeGauss" + std::to_string(2*order); + } // get Gauss points coordinates and weights gmsh::model::mesh::getIntegrationPoints(elementTypes[i], - integration_rule, + adapting_integration_rule, localCoords, weights); // get determinants and Jacobians @@ -1507,11 +1511,9 @@ void displayTime(const double &start, const double &end, const std::string &func // of the nodal forces "full_f_vector", in the FEM_domain (dimension "dim" and entity tag "tags") using the global // 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, 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) + 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) { // 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 @@ -1553,10 +1555,17 @@ void assembleKVolumicF(const Eigen::Matrix<double, 3, 3> & H_matrix, const doubl elementName, elDim, order, numNodes, localNodeCoord, numPrimaryNodes); + + // adapting the composite integration rule as a function of the order of the elements + std::string adapting_integration_rule = "CompositeGauss14"; + if(order < 8) + { + adapting_integration_rule = "CompositeGauss" + std::to_string(2*order); + } // get Gauss points coordinates and weights gmsh::model::mesh::getIntegrationPoints(elementTypes[i], - integration_rule, + adapting_integration_rule, localCoords, weights); // get determinants and Jacobians diff --git a/srcs/FEM/functionsFEM.hpp b/srcs/FEM/functionsFEM.hpp index 6f5d81ad887984f481b12384741282f927b2c3f9..cadea106899f5c85d6499d90de0a282f1cffb5ec 100644 --- a/srcs/FEM/functionsFEM.hpp +++ b/srcs/FEM/functionsFEM.hpp @@ -99,12 +99,9 @@ Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu); // with the local volumic forces computed based on the volumic force ("bx","by") and density "rho" with the help of // "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, - 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); + const double & by, 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 @@ -233,11 +230,9 @@ void displayTime(const double &start, const double &end, const std::string &func // of the nodal forces "full_f_vector", in the FEM_domain (dimension "dim" and entity tag "tags") using the global // 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, 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); + 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"). diff --git a/srcs/FEM/solverFEM.cpp b/srcs/FEM/solverFEM.cpp index a69c7e95ded2dfc408ba79100709da619345b986..97bcda367a82fab7ef26824712b9e056fac5069f 100644 --- a/srcs/FEM/solverFEM.cpp +++ b/srcs/FEM/solverFEM.cpp @@ -43,8 +43,8 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, double start = omp_get_wtime(); - /*--------Integration rule used in the code----------*/ - std::string integration_rule = "CompositeGauss4"; // tested with other methods, OK too. + /*--------Integration rule used in the code for surface integration (neumann BC)----------*/ + std::string integration_rule = "CompositeGauss4"; /*--------get all physical groups--------------*/ std::map<std::string, std::pair<int, int>> groups; @@ -114,7 +114,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, } // Looping over entities of FEM_dimain and computing local K matrices and volumic f vector - fillElementalKandF(H_matrix, rho, bx, by, integration_rule, dim, tags, nodeIndexMap, full_f_vector, FEM_kinematicsMap); + fillElementalKandF(H_matrix, rho, bx, by, dim, tags, nodeIndexMap, full_f_vector, FEM_kinematicsMap); double start3 = omp_get_wtime(); step_name = "volumic f and elemental Kl"; @@ -475,10 +475,10 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, // BEM_FEM_boundary. "nbViews" is the current number of views already allocated in the gmsh window. void solverFEM(std::map<int, double> &electrostaticPressure, int &nbViews) { - /*--------Integration rule used in the code----------*/ + /*--------Integration rule used in the code for surface integration (neumann BC)----------*/ std::string integration_rule = "CompositeGauss4"; - bool display_time = true; //displays the time performance of the different parts of the code + bool display_time = false; //displays the time performance of the different parts of the code double start_init = omp_get_wtime(); @@ -543,7 +543,7 @@ 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); + assembleKVolumicF(H_matrix, rho, bx, by, dim, tags, nodeIndexMap, full_K_matrix, full_f_vector); double end_assembleK = omp_get_wtime(); step_name = "K assembly";