Merge remote-tracking branch 'origin/master' into Pierre-Yves

README.md

 # Multiphysics 0471 
 
-Multiphysics Integrated Computational Project
+Multiphysics Integrated Computational Project
\ No newline at end of file
+
+FOR WINDOWS: 
+
+To build the FEM solver:
+- go to the srcs/FEM folder (cd srcs/FEM)
+- create a build folder (mkdir build)
+- go to this build folder (cd build)
+- cmake ..
+- make
+
+To execute the FEM solver:
+- stay in the build folder and solver.exe ..\my_geo.geo
+- for simple tension configuration: solver.exe ..\simple_tension.geo
+- for second validation configuration: solver.exe ..\uniform_charge.geo
\ No newline at end of file

examples/gmsh_api/my_code.cpp

@@ -128,6 +128,15 @@ int main(int argc, char **argv)
             for (int j = 0; j < elementTags[i].size(); ++j)
                 std::cout << elementTags[i][j] << " ";
             std::cout << '\n';
+
+            // getting element properties
+            std::string elementName;
+            int dim, order, numNodes, numPrimaryNodes;
+            std::vector<double> localNodeCoord;
+            gmsh::model::mesh::getElementProperties(elementTypes[i],
+                                            elementName, dim, order,
+                                            numNodes, localNodeCoord, 
+                                            numPrimaryNodes);
             
             // get Gauss points coordinates and weights
             std::vector<double> localCoords, weights;
@@ -148,7 +157,7 @@ int main(int argc, char **argv)
 
             // convert to Eigen format
             Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>
-                localCoords_E(&localCoords[0], localCoords.size()/3, 3);
+                localCoords_E(&localCoords[0], 3, localCoords.size()/3);
             std::cout << "\tlocalCoords (Eigen):\n";
             Eigen::IOFormat fmt(4, 0, ", ", "\n", "\t\t[", "]");
             std::cout << localCoords_E.format(fmt) << '\n';
@@ -184,39 +193,43 @@ int main(int argc, char **argv)
                 std::cout << basisFunctions[j] << " ";
             std::cout << '\n';
 
-            /*---------computing the B matrix of the first element's (i.e elementTags[i][0]) first GP:------*/
-            // ne serait-ça pas plutôt utile de directement convertir les Jacobiens en 2D ? travailler en full 2D
-            // et faire pareil pour l'évaluation des gradients des basis functions ci-dessus ?
-            std::vector<double> b_matrix;
-
-            // convert first jacobian to Eigen format
-            Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
-                jacob(&jacobians[0], 3, 3);
-            // passage en 2D xD
-            Eigen::Matrix<double, 2, 2> jacob2D;
-            jacob2D(0,0) = jacob(0,0);
-            jacob2D(1,0) = jacob(1,0);
-            jacob2D(0,1) = jacob(0,1);
-            jacob2D(1,1) = jacob(1,1);
-            std::cout << "\tfirst jacobian (Eigen):\n" 
-                      << jacob2D.format(fmt) << '\n';
-                    
-            // compute the inverse with Eigen
-            Eigen::Matrix<double, 2, 2> jacobinv = jacob2D.inverse();
-            std::cout << "\tinverse of the first jacobian (Eigen):\n" 
-                      << jacobinv.format(fmt) << "\n";
-
-            // compute the transpose of the inverse with Eigen
-            Eigen::Matrix<double, 2, 2> jacobinvtrans = jacobinv.transpose();
-            std::cout << "\ttranspose of the inverse of the first jacobian (Eigen):\n" 
-                      << jacobinvtrans.format(fmt) << "\n";
-            Eigen::Matrix<double, 3, Eigen::Dynamic> B_matrix;
-            B_matrix(0,0) = 0.;
-            std::cout << "\ttranspose of the inverse of the first jacobian (Eigen):\n" 
-                      << B_matrix.format(fmt) << "\n";
-            for (int j = 0; j < weights.size(); ++j){ // looping over the number of gauss points ( = number of shape functions (A VERIF))
+            /*---------computing the K matrix of the first element------*/
+            Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> K_matrix(2*numNodes,2*numNodes);
+            for (int j = 0; j < 2*numNodes; ++j){
+                for (int l = 0; l < 2*numNodes; ++l){
+                    K_matrix(j,l)=0.;
+                }
+            }
+            for (int j = 0; j < weights.size(); ++j){
+                // convert first jacobian to Eigen format
+                Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
+                    jacob(&jacobians[9*j], 3, 3);
+                // passage en 2D xD
+                Eigen::Matrix<double, 2, 2> jacob2D;
+                jacob2D(0,0) = jacob(0,0);
+                jacob2D(1,0) = jacob(1,0);
+                jacob2D(0,1) = jacob(0,1);
+                jacob2D(1,1) = jacob(1,1);
+                        
+                // compute the inverse with Eigen
+                Eigen::Matrix<double, 2, 2> jacobinv = jacob2D.inverse();
+
+                // compute the transpose of the inverse with Eigen
+                Eigen::Matrix<double, 2, 2> jacobinvtrans = jacobinv.transpose();
+
+                Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B_matrix(3,2*numNodes);
+                for (int l = 0; l < numNodes; ++l){ // looping over the number of shape functions
+                    B_matrix(0,2*l) = jacobinvtrans(0,0)*basisFunctions[3*l+numNodes*3*j]+jacobinvtrans(0,1)*basisFunctions[3*l+numNodes*3*j+1];
+                    B_matrix(1,2*l) = 0.;
+                    B_matrix(2,2*l) = jacobinvtrans(1,0)*basisFunctions[3*l+numNodes*3*j]+jacobinvtrans(1,1)*basisFunctions[3*l+numNodes*3*j+1];
+                    B_matrix(0,2*l+1) = 0.;
+                    B_matrix(1,2*l+1) = B_matrix(2,2*l);
+                    B_matrix(2,2*l+1) = B_matrix(0,2*l);
+                }
                 
+                K_matrix = K_matrix + B_matrix.transpose()*H_matrix*B_matrix*determinants[j]*weights[j];
             }
+            std::cout << "\tfirst elemental K matrix:\n" << K_matrix << "\n";
         }
     }
     

srcs/FEM/CMakeLists.txt

@@ -45,6 +45,8 @@ IF(NOT EIGEN_INCLUDE_DIRS)
 ENDIF()
 INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
 
+SET(CMAKE_CXX_FLAGS "${CXX_FLAGS} -Wall")
+
 # ------------------------------------------------------------------------------
 
 FILE(GLOB SRCS *.cpp)

srcs/FEM/complex_validation.geo

+//THIS CODE DEFINES THE GEOMETRY OF A CLAMPED BAR WITH TWO TYPES OF MESHES, TRIANGULAR AND QUADRANGULAR
+
+Lx = 5;
+Ly = 2;
+nx = 50;
+ny = 20;
+
+Point(1) = {0, 0, 0, 0.1};
+Point(2) = {Lx, 0, 0, 0.2};
+Point(3) = {Lx, Ly, 0, 0.2};
+Point(4) = {0, Ly, 0, 0.1};
+Point(5) = {Lx/2, 0, 0, 0.15}; 
+Point(6) = {Lx/2, Ly, 0, 0.15}; 
+Line(1) = {1, 5};
+Line(2) = {2, 5};
+Line(3) = {3, 2};
+Line(4) = {6, 3};
+Line(5) = {6, 4};
+Line(6) = {4, 1};
+Line(7) = {5, 6};
+Curve Loop(1) = {1, 7, 5, 6};
+Curve Loop(2) = {7, 4, 3, 2};
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+//Transfinite Curve {1, 5} = nx+1 Using Progression 1;
+//Transfinite Curve {7, 6} = ny+1 Using Progression 1;
+//Transfinite Curve {2, 4} = nx+1 Using Progression 1;
+//Transfinite Curve {3, 7} = ny+1 Using Progression 1;
+//Transfinite Surface {1};
+//Transfinite Surface {2};
+
+
+Mesh.ElementOrder = 2;
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Curve("left_edge", 5) = {6};
+Physical Surface("domain", 6) = {1, 2}; // the trick is to include both plane surfaces in one single domain
+Physical Curve("top_edge", 7) = {4, 5};
+Physical Curve("right_edge", 8) = {3};
+Physical Point("fixed_node", 9) = {1};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.); // ALWAYS NEED TO IMPOSE BOTH ux AND uy ON A GIVEN EDGE !! (pas très réaliste, faut y réfléchir)
+//SetNumber("Boundary Conditions/left_edge/uy", 0.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Materials/domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", 0.); // ALWAYS NEED TO IMPOSE BOTH tx AND ty ON A GIVEN EDGE (realiste, OK) !
+SetNumber("Boundary Conditions/top_edge/ty", 0.); //set to some other value for vertical deflection
+SetNumber("Boundary Conditions/right_edge/tx", 21e3); //for simple tension conditions
+SetNumber("Boundary Conditions/right_edge/ty", 0.);
+SetNumber("Volumic Forces/domain/bx",0.);
+SetNumber("Volumic Forces/domain/by",0.); //set to -9.81 for gravity
+SetNumber("Boundary Conditions/fixed_node/uy",2.);

srcs/FEM/my_geo.geo

+Lx = 5;
+Ly = 2;
+nx = 50; // prend beaucoup de temps àpd de 200x40
+ny = 20;
+
+Point(1) = {0, 0, 0, 0.1};
+Point(2) = {Lx, 0, 0, 0.1};
+Point(3) = {Lx, Ly, 0, 0.1};
+Point(4) = {0, Ly, 0, 0.2};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {2, 4} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Mesh.ElementOrder = 1;
+
+Physical Curve("left_edge", 5) = {4,5};
+Physical Surface("domain", 6) = {1};
+Physical Curve("top_edge", 7) = {3};
+Physical Curve("right_edge", 8) = {2};
+Physical Point("fixed_node", 9) = {1};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.); // HERE YOU DO NOT HAVE TO IMPOSE BOTH ux and uy simultaneously ! (permet aussi de simuler appuis à roulettes)
+//SetNumber("Boundary Conditions/left_edge/uy", 2.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Materials/domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", 0.); // ALWAYS NEED TO IMPOSE BOTH tx AND ty ON A GIVEN EDGE (realiste, OK) !
+SetNumber("Boundary Conditions/top_edge/ty", 0.); //set to some non-zero value to induce vertical deflection
+SetNumber("Boundary Conditions/right_edge/tx", 21e3); // for simple tension conditions
+SetNumber("Boundary Conditions/right_edge/ty", 0.);
+SetNumber("Volumic Forces/domain/bx",0.);
+SetNumber("Volumic Forces/domain/by",0.); //set to -9.81 for gravity
+SetNumber("Boundary Conditions/fixed_node/uy",2.);

srcs/FEM/simple_tension.geo

+Lx = 5;
+Ly = 2;
+nx = 50;
+ny = 20;
+
+Point(1) = {0, 0, 0, 1.0};
+Point(2) = {Lx, 0, 0, 1.0};
+Point(3) = {Lx, Ly, 0, 1.0};
+Point(4) = {0, Ly, 0, 1.0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {4, 2} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Curve("left_edge", 5) = {4};
+Physical Surface("domain", 6) = {1};
+Physical Curve("top_edge", 7) = {3};
+Physical Curve("right_edge", 8) = {2};
+Physical Point("fixed_node", 9) = {1};
+
+Mesh.ElementOrder = 1;
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.); // HERE YOU DO NOT HAVE TO IMPOSE BOTH ux and uy simultaneously ! (permet aussi de simuler appuis à roulettes)
+//SetNumber("Boundary Conditions/left_edge/uy", 0.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Materials/domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", 0.); // ALWAYS NEED TO IMPOSE BOTH tx AND ty ON A GIVEN EDGE (realiste, OK) !
+SetNumber("Boundary Conditions/top_edge/ty", 0.); //set to some non-zero value to induce vertical deflection
+SetNumber("Boundary Conditions/right_edge/tx", 21e3); // for simple tension conditions
+SetNumber("Boundary Conditions/right_edge/ty", 0.);
+SetNumber("Volumic Forces/domain/bx",0.);
+SetNumber("Volumic Forces/domain/by",0.); //set to -9.81 for gravity
+SetNumber("Boundary Conditions/fixed_node/uy",2.);

srcs/FEM/uniform_charge.geo

+Lx = 10;
+Ly = 1;
+nx = 200;
+ny = 20;
+
+Point(1) = {0, 0, 0, 1.0};
+Point(2) = {Lx, 0, 0, 1.0};
+Point(3) = {Lx, Ly, 0, 1.0};
+Point(4) = {0, Ly, 0, 1.0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {4, 2} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Curve("left_edge", 5) = {4};
+Physical Surface("domain", 6) = {1};
+Physical Curve("top_edge", 7) = {3};
+Physical Curve("right_edge", 8) = {2};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.); // HERE YOU DO NOT HAVE TO IMPOSE BOTH ux and uy simultaneously ! (permet aussi de simuler appuis à roulettes)
+SetNumber("Boundary Conditions/left_edge/uy", 0.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Materials/domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", 0.); // ALWAYS NEED TO IMPOSE BOTH tx AND ty ON A GIVEN EDGE (realiste, OK) !
+SetNumber("Boundary Conditions/top_edge/ty", -1); //set to some non-zero value to induce vertical deflection
+SetNumber("Boundary Conditions/right_edge/tx", 0.); // for simple tension conditions
+SetNumber("Boundary Conditions/right_edge/ty", 0.);
+SetNumber("Volumic Forces/domain/bx",0.);
+SetNumber("Volumic Forces/domain/by",0.); //set to -9.81 for gravity