first attempt normal computation

srcs/FEM/solver_FEM.cpp

+#include <gmsh.h>       // gmsh main header
+#include <iostream>     // for std::cout
+#include <map>
+#include <sstream>  // for std::stringstream
+#include <Eigen/Dense> // Eigen library
+#include <Eigen/Sparse> // Eigen library for sparse matrices
+#include <Eigen/SparseCholesky> // solving sparse linear systems
+#include <cmath>
+#include "FEM.hpp"
+
+
+int main(int argc, char **argv){
+
+    /*--------------INIT----------------*/
+    // the program requires 1 argument: a .geo file.
+    if (argc < 2)
+    {
+        std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+        return 0;
+    }
+
+    gmsh::initialize();
+
+    gmsh::open(argv[1]);
+    gmsh::model::mesh::generate(2);
+
+    /*--------Integration rule used in the code----------*/
+    // could be useful to pass it as an argument in the .geo file
+    std::string integration_rule = "CompositeGauss4"; // tested with other methods, OK too.
+
+    /*--------get physical groups--------------*/
+    std::map<std::string, std::pair<int, int>> groups;
+    gmsh::vectorpair dimTags;
+    gmsh::model::getPhysicalGroups(dimTags);
+    for (std::size_t i = 0; i < dimTags.size(); ++i)
+    {
+        int dim = dimTags[i].first;
+        int tag = dimTags[i].second;
+        std::string name;
+        gmsh::model::getPhysicalName(dim, tag, name);
+        groups[name] = {dim, tag};
+    }
+
+    /*--------get number of nodes-----------------*/
+    std::vector<std::size_t> all_nodeTags;
+    std::vector<double> all_nodecoord;
+    std::vector<double> all_nodeparametricCoord;
+
+    gmsh::model::mesh::getNodes(all_nodeTags, all_nodecoord, all_nodeparametricCoord);
+    /*
+    for (int i = 0; i < all_nodecoord.size(); ++i){
+        std::cout << "NodeCoord '" << all_nodecoord[i] << "\n";
+    }
+    for (int i = 0; i < all_nodeTags.size(); ++i){
+        std::cout << "NodeTags '" << all_nodeTags[i] << "\n";
+    
+    }
+    /*
+
+    /*-------map between nodeTag and index of first nodal displacement associated to the node (second one is the same +1)------*/
+    std::map<int, int> node_index;
+    for (std::size_t i = 0; i < all_nodeTags.size(); ++i)
+    {
+        int nodeTag = all_nodeTags[i];// use i à la place de nodeTag pour traiter cas ou le nodeTag n'est pas une suite continue
+        int node_first_index = 2*i; 
+        node_index[nodeTag] = node_first_index;
+        //std::cout << node_index[nodeTag] << " ";
+    }
+
+    /*-----------get physical properties of domain----------------------*/
+    std::vector<std::string> keys;
+    gmsh::onelab::getNames(keys, "(Volumic Forces|Materials).+");
+    double E = 0; // Young modulus
+    double nu = 0; // Poisson ratio
+    double rho = 0; // volumic mass
+    double bx = 0; // volumic force along x axis
+    double by = 0; // volumic force along y axis
+
+    Get_physical_Prop_Domain(keys, E, nu, rho, bx, by);
+
+    /*-----------------H matrix involved in stress computation------------------*/
+    Eigen::Matrix<double, 3, 3> H_matrix = Compute_H_matrix(E, nu);
+
+    /*------------computing K matrix and f vector------------------------*/
+    // we first consider the K matrix and the volumic part of the f vector
+    Eigen::SparseMatrix<double> full_K_matrix(2*all_nodeTags.size(),2*all_nodeTags.size());
+    std::vector<double> full_f_vector(2*all_nodeTags.size());//Encore besoin de ce vecteur ?
+
+    std::vector<Eigen::Triplet <double>> k_tripletList;
+
+    std::string groupname = "domain"; // we consider the domain as containing the whole 2d model
+
+    int dim = groups[groupname].first;
+    int tag = groups[groupname].second;
+    if (tag == 0)
+    {
+        std::cerr << "Group '" << groupname << "' does not exist!\n";
+        return 1;
+    }
+
+    // Getting entities of the domain
+    std::vector<int> tags;
+    gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+    std::vector<double> normalElement(2);
+    std::vector<int> ez = {0, 0, 1};
+
+    // Looping over entities and compute K and F
+    Loop_K_F(H_matrix, rho, bx, by, integration_rule, all_nodecoord, normalElement, ez, dim, tag, tags, full_K_matrix, 
+              full_f_vector, k_tripletList);
+
+    full_K_matrix.setFromTriplets(k_tripletList.begin(), k_tripletList.end());
+
+
+    /*-----------------BOUNDARY CONDITIONS------------*/
+    /*-----------get group_names related to BC's----------------------*/
+    gmsh::onelab::getNames(keys, "(Boundary Conditions).+");
+    /*-----------first focus on neumann BC's (surface traction) into full_f_vector------------*/
+    Neumann_BC(keys, integration_rule, dim, tag, tags, node_index, groups, groupname, full_f_vector);
+
+    /*------------FOCUS ON DIRICHLET BC'S---------------*/
+    // new_indices vector will contain the new indices of the K matrix after the rows/columns corresponding
+    // to dirichlet boundary conditions have been removed
+    // IN THIS IMPLEMENTATION WE DECIDE TO REMOVE EVERY ROW/COL ASSOCIATED TO DIR. BC's
+    // non-homo. dir BC will induce a correction in the RHS term f of the other rows, while homogeneous not
+    // if the given index has been removed, it will be set to -1 in new_indices
+    std::vector<int> new_indices(2*all_nodeTags.size());
+    for (std::size_t j = 0; j < new_indices.size(); ++j){
+        new_indices[j] = (int) j;
+    } 
+
+    // dir_BC[index] will contain the dirichlet boundary condition prescribed to this node
+    std::vector<double> dir_BC(2*all_nodeTags.size());
+    for (std::size_t j = 0; j < dir_BC.size(); ++j){
+        dir_BC[j] = 0;
+    }
+
+    // Fill the two vectors defined just above
+    Fill_New_indices_Dir_BC(keys, dim, tag, tags, node_index, groups, groupname, new_indices, dir_BC);
+
+    int number_removed_lines = 0;
+    for (std::size_t j = 0; j < new_indices.size(); ++j){
+        if(new_indices[j] < 0){ // line has to be removed
+            number_removed_lines++;
+        }
+        else{
+            new_indices[j] -= number_removed_lines;
+        }
+    }
+
+
+    /*------------FILLING THE REDUCED K MATRIX AND F VECTOR--------------*/
+    //the correction term on the f vector comes from the "FEM" theoretical course (JP Ponthot), Chap2 Slide 38/50
+    Eigen::SparseMatrix<double> reduced_K_matrix(2*all_nodeTags.size() - number_removed_lines,2*all_nodeTags.size() - number_removed_lines);
+    std::vector<double> reduced_f_vector(2*all_nodeTags.size() - number_removed_lines);
+
+    std::vector<Eigen::Triplet <double>> reduced_k_tripletList;
+
+    Fill_K_F_Reduced(full_f_vector, full_K_matrix, new_indices, dir_BC, reduced_K_matrix, reduced_f_vector, reduced_k_tripletList);
+
+    /*--------------SOLVING THE REDUCED SYSTEM--------------------*/
+    Eigen::VectorXd final_reduced_d(2*all_nodeTags.size() - number_removed_lines);
+
+    Eigen::VectorXd final_reduced_f(2*all_nodeTags.size() - number_removed_lines);
+    for (std::size_t j = 0; j < reduced_f_vector.size(); ++j){
+        final_reduced_f(j) = reduced_f_vector[j];
+    }
+
+    // solve Kd=f
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+    solver.compute(reduced_K_matrix);
+
+    final_reduced_d = solver.solve(final_reduced_f);
+
+
+    /*-----------RECONSTRUCTING FULL NODAL VALUES--------------*/
+    std::vector<double> full_d_vector(2*all_nodeTags.size());
+    for (std::size_t j = 0; j < full_d_vector.size(); ++j){
+        if(new_indices[j] < 0){ // line was previously removed due to dirichlet boundary conditions
+            full_d_vector[j] = dir_BC[j];
+        }
+        else{
+            full_d_vector[j] = final_reduced_d(new_indices[j]);
+        }
+    }
+
+    /*--------------COMPUTING NODAL FORCES-------------*/
+    Eigen::VectorXd eigen_d(2*all_nodeTags.size());
+    for(size_t i = 0 ; i < 2*all_nodeTags.size(); ++i){
+        eigen_d(i) = full_d_vector[i];
+    }
+    Eigen::VectorXd nodal_forces(2*all_nodeTags.size());
+    nodal_forces = full_K_matrix*eigen_d;
+
+    /*--------------REACTION FORCES COMPUTATION------------------*/
+
+    double Rx = 0.;
+    double Ry = 0.;
+    Reaction_Forces(Rx, Ry, dim, tag, tags, groupname, groups, node_index, nodal_forces, keys);
+
+    /*--------------PLOTTING AND SAVING NODAL VALUES------------------*/
+    // loop over the nodes
+    std::vector<std::string> names = Plot_Save_ux_uy(all_nodeTags, full_d_vector);
+
+    /*-------------POST PROCESSING---------------*/
+    /*-------------Computing and saving strains/stresses-----------------*/
+
+    // get all the elements for dim = 2
+    std::vector<int> elementTypes;
+    std::vector<std::vector<std::size_t>> elementTags;
+    std::vector<std::vector<std::size_t>> elnodeTags;
+    gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, 2);
+    
+    // compute total number of 2D elements
+    int nbelems=0;
+    for(std::size_t i=0; i< elementTags.size(); ++i)
+        nbelems+=elementTags[i].size();
+
+
+    /*-------------------ELEMENTAL-NODAL VALUES-----------------------*/
+    // éventuellement faire un tableau ici qui regroupe tout le bazar pour ne pas avoir la même ligne 7 fois ...
+    int viewtag3 = gmsh::view::add("nodal_eps_x");
+    int viewtag4 = gmsh::view::add("nodal_eps_y");
+    int viewtag5 = gmsh::view::add("nodal_2eps_xy");
+    int viewtag6 = gmsh::view::add("nodal_eps_z");
+    int viewtag7 = gmsh::view::add("nodal_sig_x");
+    int viewtag8 = gmsh::view::add("nodal_sig_y");
+    int viewtag9 = gmsh::view::add("nodal_sig_xy");
+    int viewtag10 = gmsh::view::add("nodal_sig_VM");
+    std::vector<std::vector<double>> data3(nbelems);
+    std::vector<std::vector<double>> data4(nbelems);
+    std::vector<std::vector<double>> data5(nbelems);
+    std::vector<std::vector<double>> data6(nbelems);
+    std::vector<std::vector<double>> data7(nbelems);
+    std::vector<std::vector<double>> data8(nbelems);
+    std::vector<std::vector<double>> data9(nbelems);
+    std::vector<std::vector<double>> data10(nbelems);
+    std::vector<std::size_t> elems2D(nbelems);
+
+
+    std::vector<int> nb_adjacent_elements(all_nodeTags.size()); // number of 2d elements adjacent to node (i+1)
+    std::vector<double> nodal_sigma_x(all_nodeTags.size()); // will store nodal values (may be helpful later)
+    std::vector<double> nodal_sigma_y(all_nodeTags.size());
+    std::vector<double> nodal_sigma_xy(all_nodeTags.size());
+
+    Compute_Stress_Strain_Elemental_Nodal(elementTypes, elementTags, elnodeTags, full_d_vector, elems2D, H_matrix, 
+                                    nu, E, nb_adjacent_elements, nodal_sigma_x, nodal_sigma_y, nodal_sigma_xy, node_index,
+                                    data3, data4, data5, data6, data7, data8, data9, data10);
+
+    gmsh::view::addModelData(viewtag3, 0, names[0], "ElementNodeData",
+                                elems2D, data3, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag4, 0, names[0], "ElementNodeData",
+                                elems2D, data4, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag5, 0, names[0], "ElementNodeData",
+                                elems2D, data5, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag6, 0, names[0], "ElementNodeData",
+                                elems2D, data6, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag7, 0, names[0], "ElementNodeData",
+                                elems2D, data7, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag8, 0, names[0], "ElementNodeData",
+                                elems2D, data8, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag9, 0, names[0], "ElementNodeData",
+                                elems2D, data9, 1);  // the last ,1 is important!
+    gmsh::view::addModelData(viewtag10, 0, names[0], "ElementNodeData",
+                                elems2D, data10, 1);  // the last ,1 is important!
+
+    // representing vector displacement field and nodal forces
+    int viewtag11 = gmsh::view::add("u");
+    int viewtag12 = gmsh::view::add("f");
+    std::vector<std::vector<double>> data11(all_nodeTags.size());
+    std::vector<std::vector<double>> data12(all_nodeTags.size());
+    for (std::size_t i = 0; i < all_nodeTags.size(); ++i)
+    {
+        data11[i].resize(3);
+        data12[i].resize(3);
+        data11[i][0] = full_d_vector[2*i]; //ux
+        data11[i][1] = full_d_vector[2*i+1]; //uy
+        data11[i][2] = 0.; //uz
+        data12[i][0] = nodal_forces(2*i); //fx
+        data12[i][1] = nodal_forces(2*i+1); //fy
+        data12[i][2] = 0.; //fz
+    }
+
+    gmsh::view::addModelData(viewtag11, 0, names[0], "NodeData",
+                             all_nodeTags, data11); //independent of time here
+    gmsh::view::addModelData(viewtag12, 0, names[0], "NodeData",
+                             all_nodeTags, data12); //independent of time here
+
+
+
+    // empeche que tout se plot l'un sur l'autre
+    int nb_views = 12; //hardcoded
+    for( int i = 0; i < nb_views; ++i){
+        std::string my_string = "View[" + std::to_string(i) + "].Visible";
+        gmsh::option::setNumber(my_string, 0);
+    }
+
+
+    
+    gmsh::fltk::run(); // opens the gmsh window
+    
+    gmsh::finalize();
+    return 0;
+}
\ No newline at end of file