debug multiple FEM domains

srcs/FEM/functionsFEM.cpp

@@ -129,7 +129,7 @@ void fillElementalKandF(const Eigen::Matrix<double, 3, 3> H_matrix, const double
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[i], 
                                             localCoords, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[k]);
 
             // values and derivatives of the shape functions at the Gauss points of the reference element
             // we are here working with derivatives in the reference space, same for all elements
@@ -334,7 +334,7 @@ void includeNeumannBC(const std::vector<std::string> & keys, const std::string &
                         std::vector<double> jacobians, determinants, coords;
                         gmsh::model::mesh::getJacobians(elementTypes[i], 
                                                         localCoords, jacobians, 
-                                                        determinants, coords,tags[k]); // ai dû rajouter mon tags[k] pr focus on this particular entity
+                                                        determinants, coords, tags[k]); // ai dû rajouter mon tags[k] pr focus on this particular entity
 
                         gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords, 
                                                             "Lagrange", numComponents, 
@@ -600,7 +600,7 @@ void computeStressStrainElNodal(int & dim, std::vector<int> & tags,
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[ityp], 
                                             paramCoord, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[entityTag]);
 
             // derivatives of the shape functions at the nodes of the reference element
             // we are here working with derivatives in the reference space, same for all elements
@@ -1275,7 +1275,7 @@ void computeStrainEnergy(double & strainEnergy, int & dim, std::vector<int> & ta
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[i], 
                                             localCoords, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[k]);
 
             // values and derivatives of the shape functions at the Gauss points of the reference element
             // we are here working with derivatives in the reference space, same for all elements
@@ -1424,7 +1424,9 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[i], 
                                             localCoords, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[k]); //big modif
+
+            //std::cout << "entity " << k << ", element type " << i << ",jacobians size: " << jacobians.size() << ", weights size: " << weights.size() << "\n";
 
             // values and derivatives of the shape functions at the Gauss points of the reference element
             // we are here working with derivatives in the reference space, same for all elements
@@ -1438,6 +1440,8 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
                                                  basisFunctionsGradient, 
                                                  numOrientations);
 
+            //std::cout << "entity " << k << ", element tags " << elementTags[i].size() << ",jacobians size: " << jacobians.size() << ",determinants size: " << determinants.size() << ", weights size: " << weights.size() << ", basis functions size: " << basisFunctions.size() << ", basis functions gradient size: " << basisFunctionsGradient.size() << "\n";
+
             /*---------computing the K matrix (and f vector, only volumic part for now)------*/
             // looping over the elements (computing elemental K matrices)
             //double test1 = omp_get_wtime();
@@ -1468,6 +1472,8 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
                     jacob2D(1,0) = jacob(1,0);
                     jacob2D(0,1) = jacob(0,1);
                     jacob2D(1,1) = jacob(1,1);
+
+                    //std::cout << "element tag: " << elementTags[i][el] << "jacob2d: " << jacob2D << "\n";
                             
                     // computing the inverse
                     Eigen::Matrix<double, 2, 2> jacobinv = jacob2D.inverse();
@@ -1480,6 +1486,8 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
 
 
                     computeBmatrix(j, numNodes, basisFunctionsGradient, determinants, jacobinvtrans, B_matrix);
+
+                    //std::cout << "element tag: " << elementTags[i][el] << "B matrix: " << B_matrix << "\n";
                     
                     for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
                         f_vector[2*l] = f_vector[2*l] + rho*basisFunctions[l+numNodes*j]*bx*determinants[j+weights.size()*el]*weights[j]; //pas oublier le dét !
@@ -1487,11 +1495,14 @@ void assembleKlistVolumicF(const Eigen::Matrix<double, 3, 3> H_matrix, const dou
                     }
 
                     K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[j+weights.size()*el]*weights[j];
+                    //std::cout << "determinant: " << determinants[j+weights.size()*el] << ", weight: " << weights[j] << "\n";
                 }
                 /*------------- ASSEMBLY PROCESS ------------*/
                 // USE LIST OF TRIPLETS TO PUSH ONLY ONCE -> more rapid and easier to parallelize
                 assembleKtriplet(i, el, nodeTags, numNodes, nodeIndexMap, my_k_list[omp_get_thread_num()], K_matrix, f_vector);
                 assembleFthread(i, el, nodeTags, numNodes, nodeIndexMap, thread_f_vector[omp_get_thread_num()], f_vector);
+
+                //std::cout << "element tag: " << elementTags[i][el] << "local K matrix: " << K_matrix << "\n";
             }
             //double test2 = omp_get_wtime() - test1;
             //std::cout << "loop time: " << test2 << "\n";
@@ -1593,7 +1604,7 @@ void computeStressStrainLinearWay(int & dim, std::vector<int> & tags,
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[ityp], 
                                             paramCoord, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[entityTag]);
 
             // derivatives of the shape functions at the nodes of the reference element
             // we are here working with derivatives in the reference space, same for all elements
@@ -1705,7 +1716,7 @@ void computeStrainEnergyLinearWay(double & strainEnergy, int & dim, std::vector<
             std::vector<double> jacobians, determinants, coords;
             gmsh::model::mesh::getJacobians(elementTypes[i], 
                                             localCoords, jacobians, 
-                                            determinants, coords);
+                                            determinants, coords, tags[k]);
 
             // values and derivatives of the shape functions at the Gauss points of the reference element
             // we are here working with derivatives in the reference space, same for all elements