Correction Compute Normal

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],