diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp index 9cd7a082d8a8e7dec9a4c01e9318d2d1032001a7..27ec51c66f0290e3b5486ea01358e326146865f2 100644 --- a/srcs/FEM/functionsFEM.cpp +++ b/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