From 1bb04cf38a2c83c93ee3afbcbb7a94f0bff35306 Mon Sep 17 00:00:00 2001
From: Louis Denis <louis.denis@student.uliege.be>
Date: Tue, 17 May 2022 19:17:46 +0200
Subject: [PATCH] gauss integration
---
srcs/FEM/functionsFEM.cpp | 31 ++++++++++++++++++++-----------
srcs/FEM/functionsFEM.hpp | 13 ++++---------
srcs/FEM/solverFEM.cpp | 12 ++++++------
3 files changed, 30 insertions(+), 26 deletions(-)
diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp
index 4680d1d..4e39d98 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 6f5d81a..cadea10 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 a69c7e9..97bcda3 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";
--
GitLab