k matrix in FEM srcs

srcs/FEM/my_code.cpp

+// This program takes a .geo file as argument, generates the mesh of the 
+// surfaces and displays the result in the Gmsh window.
+// 
+// This example introduces:
+// - the gmsh.h header
+// - initialize()/finalize() functions
+// - C++ namespaces / standard library (std::cout)
+//
+// How to build/run? (windows) 
+//   [in examples/gmsh_api]
+//   mkdir build
+//   cd build
+//   cmake .. && make && code1_loadgeo.exe ..\rectangle.geo
+
+#include <gmsh.h>       // gmsh main header (READ THIS FILE)
+#include <iostream>     // for std::cout
+#include <map>
+#include <sstream>  // for std::stringstream
+#include <Eigen/Dense> // Eigen library
+
+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(); // first thing to do before using any gmsh SDK function
+
+    gmsh::open(argv[1]); // similar to "File / Open" in the GUI
+    gmsh::model::mesh::generate(2);  // similar to "Mesh / 2D" in the GUI
+
+    /*--------get physical groups (more particularly: the domain)--------------*/
+    std::map<std::string, std::pair<int, int>> groups;
+    gmsh::vectorpair dimTags;
+    gmsh::model::getPhysicalGroups(dimTags);
+    for (int 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 physical properties of domain----------------------*/
+    std::vector<std::string> keys;
+    gmsh::onelab::getNames(keys, "(Materials).+");
+    double E;
+    double nu;
+    for (auto &key : keys)
+    {
+        // get corresponding value
+        std::vector<double> value;
+        gmsh::onelab::getNumber(key, value);
+        // expected key structure is "type/group_name/field"
+        // => split the key string into components
+        std::stringstream ss(key);
+        std::vector<std::string> words;
+        std::string word;
+        while(std::getline(ss, word, '/')) // read string until '/'
+            words.push_back(word);
+        if(words.size()==3) 
+        {
+            if(words[2].compare("Young") == 0){
+                E = value[0];
+            }
+            if(words[2].compare("Poisson") == 0){
+                nu = value[0];
+            }
+        }
+    }
+    // std::cout << "E = " << E << "\n";
+    // std::cout << "nu = " << nu << "\n";
+
+    /*-----------------H matrix------------------*/
+    Eigen::Matrix<double, 3, 3> H_matrix;
+    H_matrix(0,0) = 1.;
+    H_matrix(1,0) = nu;
+    H_matrix(2,0) = 0.;
+    H_matrix(0,1) = nu;
+    H_matrix(1,1) = 1.;
+    H_matrix(2,1) = 0.;
+    H_matrix(0,2) = 0.;
+    H_matrix(1,2) = 0.;
+    H_matrix(2,2) = (1-nu)/2;
+    H_matrix = E/(1-nu*nu)*H_matrix;
+    // std::cout << H_matrix << '\n'; 
+
+    /*------------loop over elements of domain------------------------*/
+    std::string groupname = "domain";
+
+    int dim = groups[groupname].first;
+    int tag = groups[groupname].second;
+    if (tag == 0)
+    {
+        std::cerr << "Group '" << groupname << "' does not exist!\n";
+        return 1;
+    }
+
+    // get entities of the domain
+    std::vector<int> tags;
+    gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+    std::cout << "Entities in group named '" << groupname << "': ";
+    for (int i = 0; i < tags.size(); ++i)
+        std::cout << tags[i] << " ";
+    std::cout << '\n';
+
+    // loop over entities
+    for (int k = 0; k < tags.size(); ++k)
+    {
+        std::cout << "Entity (" << dim << "," << tags[k] << "):\n";
+        std::vector<int> elementTypes;
+        std::vector<std::vector<std::size_t>> elementTags;
+        std::vector<std::vector<std::size_t>> nodeTags;
+        gmsh::model::mesh::getElements(elementTypes, elementTags,
+                                       nodeTags, dim, tags[k]);
+        
+        // loop over element types
+        for (int i = 0; i < elementTypes.size(); ++i)
+        {
+            std::cout << "\telement type " << elementTypes[i] << ":\n";
+            std::cout << "\t\telements: ";
+            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;
+            gmsh::model::mesh::getIntegrationPoints(elementTypes[i], 
+                                                    "CompositeGauss2", // "Gauss2"
+                                                    localCoords, weights);
+            std::cout << "\t" << weights.size() << " integration points\n";
+
+            std::cout << "\tweights: ";
+            for (int j = 0; j < weights.size(); ++j)
+                std::cout << weights[j] << " ";
+            std::cout << '\n';
+
+            std::cout << "\tlocalCoords (gmsh): ";
+            for (int j = 0; j < localCoords.size(); ++j)
+                std::cout << localCoords[j] << " ";
+            std::cout << '\n';
+
+            // convert to Eigen format
+            Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>
+                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';
+
+            // get determinants and Jacobians
+            std::vector<double> jacobians, determinants, coords;
+            gmsh::model::mesh::getJacobians(elementTypes[i], 
+                                            localCoords, jacobians, 
+                                            determinants, coords);
+            std::cout << "\tGot " << determinants.size() 
+                    << " Jacobians for type " << elementTypes[i] << "\n";
+
+            std::cout << "\tdeterminants (gmsh): ";
+            for (int j = 0; j < determinants.size(); ++j)
+                std::cout << determinants[j] << " ";
+            std::cout << '\n';
+
+            std::cout << "\tjacobians (gmsh): ";
+            for (int j = 0; j < jacobians.size(); ++j)
+                std::cout << jacobians[j] << " ";
+            std::cout << '\n';
+
+            // 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
+            int numComponents, numOrientations;
+            std::vector<double> basisFunctions;
+            gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords, 
+                                                 "GradLagrange", numComponents, 
+                                                 basisFunctions, 
+                                                 numOrientations);
+            std::cout << "\tbasis functions (gmsh): ";
+            for (int j = 0; j < basisFunctions.size(); ++j)
+                std::cout << basisFunctions[j] << " ";
+            std::cout << '\n';
+
+            /*---------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";
+        }
+    }
+    
+    gmsh::finalize(); // clean-up things before ending the program
+    return 0;
+}

srcs/FEM/my_geo.geo

+Lx = 5;
+Ly = 2;
+nx = 5;
+ny = 1;
+
+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};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.);
+SetNumber("Boundary Conditions/left_edge/uy", 0.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", -100.);
\ No newline at end of file