diff --git a/srcs/FEM/FEM_Functions.cpp b/srcs/FEM/FEM_Functions.cpp index 71dca76d67731b32da25035c815a997ddca703f2..33331a85dbaf5be20228e9c3de21d5f07097aec0 100644 --- a/srcs/FEM/FEM_Functions.cpp +++ b/srcs/FEM/FEM_Functions.cpp @@ -70,6 +70,28 @@ void Compute_B_matrix(const int j, const int numNodes, const std::vector<double> } } +//Compute the normal of one element +void Compute_Normal(std::vector<double> & normalElement, std::vector<int> ez, std::vector<double> all_nodecoord, + std::vector<std::vector<std::size_t>> elementTags, std::vector<std::size_t> all_nodeTags, int i){ + + for (std::size_t el = 0; el < elementTags[i].size(); el++){ + + int it1 = all_nodeTags[2*el]; + int it2 = all_nodeTags[2*el+1]; + double x1 = all_nodecoord[3*(it1 - 1)]; + double y1 = all_nodecoord[3*(it1 - 1) + 1]; + + double x2 = all_nodecoord[3*(it2 - 1)]; + double y2 = all_nodecoord[3*(it2 - 1) + 1]; + + std::vector<double> vectorElement; + vectorElement[0] = x2 - x1; + vectorElement[1] = y2 - y1; + + normalElement[0] = vectorElement[1] * ez[2] ; + normalElement[1] = -vectorElement[1] * ez[2]; + } +} //Get physical properties of domain void Get_physical_Prop_Domain(const std::vector<std::string> keys, double & E, double & nu, double & rho, @@ -111,7 +133,8 @@ void Get_physical_Prop_Domain(const std::vector<std::string> keys, double & E, d // we first consider the K matrix and the volumic part of the f vector void Loop_K_F(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, const double bx, const double by, - const std::string integration_rule, + const std::string integration_rule, const std::vector<double> all_nodecoord, std::vector<double> & normalElement, + const std::vector<int> ez, int & dim, int & tag, std::vector<int> & tags, Eigen::SparseMatrix<double> & full_K_matrix, std::vector<double> & full_f_vector, @@ -157,6 +180,7 @@ void Loop_K_F(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, cons int numComponents, numOrientations; std::vector<double> basisFunctions; std::vector<double> basisFunctionsGradient; + gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords, "Lagrange", numComponents, basisFunctions, @@ -166,7 +190,6 @@ void Loop_K_F(const Eigen::Matrix<double, 3, 3> H_matrix, const double rho, cons "GradLagrange", numComponents, basisFunctionsGradient, numOrientations); - /*---------computing the K matrix (and f vector, only volumic part for now)------*/ // looping over the elements (computing elemental K matrices) for (std::size_t el = 0; el < elementTags[i].size(); el++){ @@ -230,9 +253,11 @@ void Compute_Neumann_BC(const int numNodes, const int i, const std::vector<std:: std::vector<double> & full_f_vector){ for (std::size_t el = 0; el < elementTags[i].size(); el++){ + for (int j = 0; j < 2*numNodes; ++j){ f_vector[j] = 0.; } + // looping over the Gauss points (formula of slide 17) for (std::size_t j = 0; j < weights.size(); ++j){ // looping over the number of shape functions @@ -314,13 +339,14 @@ void Neumann_BC(const std::vector<std::string> keys, const std::string integrati elementName, dim, order, numNodes, localNodeCoord, numPrimaryNodes); - + // get Gauss points coordinates and weights std::vector<double> localCoords, weights; gmsh::model::mesh::getIntegrationPoints(elementTypes[i], integration_rule, //maybe use different integration method on boundary than in domain localCoords, weights); + // get determinants and Jacobians (only determinants will be useful) std::vector<double> jacobians, determinants, coords; gmsh::model::mesh::getJacobians(elementTypes[i],