delete

README.md

@@ -18,4 +18,69 @@ To build the independent FEM solver:
 - similar procedure but to be done in the srcs/FEM folder
 
 To build the independent BEM solver:
-- similar procedure but to be done in the srcs/BEM folder
\ No newline at end of file
+- similar procedure but to be done in the srcs/BEM folder
+
+## Creation of a .geo file 
+
+To create a .geo file compatible with the program, some rules must be respected:  
+- For the elastic FEM part:    
+  1. The FEM domain should be defined as a `Physical Surface` with the precise name `FEM_domain`.
+        First, a `Curve Loop` gathering the lines representing the boundary of the FEM domain must be created.
+        Then, a `Plane Surface` containing the corresponding curve loop must be created.
+        Finally, if the plane surface is the first of the .geo file, the physical surface must be created as: 
+        `Physical Surface("FEM_domain", x) = {1};`,
+        where x is the tag of the physical surface.
+        Multiple sub domains can be included in the FEM domain. For example, if plane surfaces have been created for
+        five sub domains, the physical surface must be created as:
+        <pre><code>Physical Surface("FEM_domain", x) = {1, 2, 3, 4, 5};.<code><pre> 
+    
+  2. The mechanical boundary conditions, material properties and volumic forces must be specified for the FEM domain:  
+    - **BOUNDARY CONDITIONS**:
+      1. The edges on which the user wants to impose a boundary condition must be defined as `Physical Curve`.
+                As an example, if the user want to impose a condition on the bottom edge related to the `Line(1)`: 
+                <pre><code>`Physical Curve("bottom_edge", x) = {1};`<code><pre>
+      2. All boundary conditions must be written as `SetNumber("Boundary Conditions/name_of_the_physical_curve/...")` 
+                with `name_of_the_physical_curve` the physical curve on which the condition is imposed.
+      3. The Dirichlet boundary conditions, i.e the imposed horizontal and vertical displacement `u_x` and `u_y` of an edge,    
+                expressed in [m], must be written as `Setnumber("Boundary Conditions/name_of_the_physical_curve/ux", desired_value);`, same applies for `u_y`. Note that the horizontal displacement can be imposed on a physical curve without imposing the vertical displacement, and vice-versa.
+                As an example, if the user wants the left edge to be clamped: 
+                    <pre><code>SetNumber("Boundary Conditions/left_edge/ux", 0.);
+                        SetNumber("Boundary Conditions/left_edge/uy", 0.);<code><pre>
+      4. The Neumann boundary conditions, i.e surface traction in the horizontal (`t_x`) or vertical (`t_y`) direction imposed on
+                an edge, expressed in [Pa], must be written as `SetNumber("Boundary Conditions/name_of_the_physical_curve/tx", desired_value)`, 
+                same applies for `t_y`. Note that on a specific edge, both horizontal and vertical surface tractions must be 
+                prescribed. If the user only wants to prescribe `t_x`, `t_y` must also be set to zero.
+      5. Dirichlet boundary conditions can also be imposed on a specific point of the domain. In this case, a `Physical Point`
+                must be created.
+    - **MATERIAL PROPERTIES**:
+                All the material properties must be written as `SetNumber("Materials/FEM_domain/name_of_property", value_of_property);`.
+                The different properties that must be specified are:
+      - Young's modulus, expressed in [Pa], named `Young`.
+      - Poisson's ratio, without physical units, named `Poisson`.
+      - Mass density, expressed in [kg/m3], named `rho`.
+    - **VOLUMIC FORCES**:
+                Volumic forces in the horizontal (`b_x`) or vertical (`b_y`) direction, expressed in [m/s2], must be written as `SetNumber("Volumic Forces/FEM_domain/b_x",0);`, same applies for `b_y`.
+- For the electrostatic BEM part:
+  1. The lines corresponding to the boundary of the BEM surfaces must be created in such a way that the area of the BEM surface is located on the left of the Curve Loop.
+  2. Multiple BEM domains can be defined. One BEM domain should be defined as a `Physical Surface` with the precise name `BEM_domain_X`, with `X` the identifier of the BEM domain (`X=1` if it is the first domain, `X=2` if it is the second one, ...).
+  3. Boundary conditions and material properties must then be specified for the BEM domain: 
+    - **BOUNDARY CONDITIONS**:
+      1. Dirichlet boundary conditions, corresponding to imposing an electric potential expressed in [V], can be imposed on a physical curve as:
+       `SetNumber("Boundary Conditions/name_of_the_physical_curve/BEM_domain_X/dirichlet", value_of_potential);`.
+      2. Neumann boundary conditions, corresponding to imposing the normal component of the electric field expressed in [V/m], can be imposed on a physical curve as: `SetNumber("Boundary Conditions/name_of_the_physical_curve/BEM_domain_X/neumann", value_of_field);`.
+    - **MATERIAL PROPERTIES**: The dielectric permittivity of the material inside the BEM domain must be specified. It is done writing `SetNumber("Materials/BEM_domain_X/Epsilon", value_of_Epsilon);`.
+
+  **Note**: These operations must be done for each BEM domain the user wants to create, replacing `X` by the given identifier of the BEM domain. 
+- For the coupling between the FEM and BEM domains:
+  1. All the lines belonging to the interface between the FEM domain and the different BEM domains must be specified as a single `Physical Curve("BEM_FEM_boundary, x)={List_of_the_interface_lines}`. As an example, if `Line(1)`, `Line(4)` and `Line(5)` belong to the intersection of the BEM domains and the FEM domain: 
+        <pre><code> Physical Curve("BEM_FEM_boundary", x) = {1, 4, 5};<code><pre>
+- For the mechanical or the coupled solver, the user can choose between the linear or non-linear iterative solver, by setting 
+  `SetNumber("Non_linear_solver",0);` for the linear solver, or `SetNumber("Non_linear_solver",1);` for the non-linear solver.
+
+
+
+
+
+
+
+

srcs/BEM/functionsBEM.hpp

@@ -206,4 +206,4 @@ double computeAnalyticalGradHx(const double x, const double y, const elementStru
 double computeAnalyticalGradHy(const double x, const double y, const elementStruct &element);
 double computeAnalyticalGradGx(const double x, const double y, const elementStruct &element);
 double computeAnalyticalGradGy(const double x, const double y, const elementStruct &element);
-double computeAnalyticalG(const double x, const double y, const elementStruct &element);
\ No newline at end of file
+double computeAnalyticalG(const double x, const double y, const elementStruct &element);

srcs/BEM/mainBEM.cpp

@@ -22,7 +22,7 @@ int main(int argc, char **argv)
     std::map<int,std::pair<double,double>> nodalDisplacements;
     const bool postProcessing = true;
     const bool meshUntangler = false;
-    
+
     solverBEM(electrostaticPressure, nbViews, nodalDisplacements, postProcessing, 1, meshUntangler); // iteration number randomly set to 1
 
     if(postProcessing)

srcs/FEM/complex_validation.geo

@@ -2,15 +2,15 @@
 
 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}; 
+nx = 250;
+ny = 100;
+
+Point(1) = {0, 0, 0, 0.05}; //4th number is mesh density
+Point(2) = {Lx, 0, 0, 0.1};
+Point(3) = {Lx, Ly, 0, 0.1};
+Point(4) = {0, Ly, 0, 0.05};
+Point(5) = {Lx/2, 0, 0, 0.075}; 
+Point(6) = {Lx/2, Ly, 0, 0.075}; 
 Line(1) = {1, 5};
 Line(2) = {2, 5};
 Line(3) = {3, 2};
@@ -52,6 +52,8 @@ SetNumber("Boundary Conditions/right_edge/tx", 21e3); //for simple tension condi
 SetNumber("Boundary Conditions/right_edge/ty", 0.);
 SetNumber("Volumic Forces/FEM_domain/bx",0.);
 SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
-SetNumber("Boundary Conditions/fixed_node/uy",0.);
+SetNumber("Boundary Conditions/fixed_node/uy",2);
 
 Physical Curve("BEM_FEM_boundary", 10) = {1,2,3,4,5,6};
+
+SetNumber("Non_linear_solver", 0);
\ No newline at end of file

srcs/FEM/large_rotation_validation.geo

+
+h = 1;
+H = 10;
+
+n = 16;
+
+Point(1) = {0, 0, 0, 0.5};
+Point(2) = {0.9*H, 0, 0, 0.5};
+Point(3) = {0.9*H, h, 0, 0.5};
+Point(4) = {0, h, 0, 0.5};
+
+Point(5) = {0.95*H, 0, 0, 0.5};
+Point(6) = {0.95*H, h, 0, 0.5};
+
+Point(7) = {0.95*H, H, 0, 0.5};
+Point(8) = {0.9*H, H, 0, 0.5};
+
+Point(9) = {H, 0, 0, 0.5};
+Point(10) = {H, h, 0, 0.5};
+Point(11) = {H, H, 0, 0.5};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {1, 2, 3, 4};
+Plane Surface(1) = {1};
+
+Line(5) = {2, 5};
+Line(6) = {5, 6};
+Line(7) = {6, 3};
+Curve Loop(2) = {5, 6, 7, -2};
+Plane Surface(2) = {2};
+
+Line(8) = {6, 7};
+Line(9) = {7, 8};
+Line(10) = {8, 3};
+Curve Loop(3) = {-7, 8, 9, 10};
+Plane Surface(3) = {3};
+
+Line(11) = {5, 9};
+Line(12) = {9, 10};
+Line(13) = {10, 6};
+Curve Loop(4) = {11, 12, 13, -6};
+Plane Surface(4) = {4};
+
+Line(14) = {10, 11};
+Line(15) = {11, 7};
+Curve Loop(5) = {14, 15, -8, -13};
+Plane Surface(5) = {5};
+
+Transfinite Curve {1, 3, 8, 10, 14} = 18*n+1 Using Progression 1;
+Transfinite Curve {2, 4, 6, 12} = 2*n+1 Using Progression 1;
+Transfinite Curve {5, 7, 9, 11, 13, 15} = n+1 Using Progression 1;
+Transfinite Surface {1};
+Transfinite Surface {2};
+Transfinite Surface {3};
+Transfinite Surface {4};
+Transfinite Surface {5};
+
+Recombine Surface {1};
+Recombine Surface {2};
+Recombine Surface {3};
+Recombine Surface {4};
+Recombine Surface {5};
+
+Physical Curve("left_edge", 1) = {4};
+Physical Surface("FEM_domain", 2) = {1, 2, 3, 4, 5}; // the trick is to include both plane surfaces in one single domain
+Physical Curve("top_edge", 3) = {9, 15};
+
+Physical Point("point_A", 5) = {7};
+
+F = 40000;
+
+// 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/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Poisson", 0.3);
+SetNumber("Materials/FEM_domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", F); // 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("Volumic Forces/FEM_domain/bx",0.);
+SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
+
+SetNumber("Non_linear_solver",1); // activate non linear solver
+
+Physical Curve("BEM_FEM_boundary", 4) = {4};
\ No newline at end of file

srcs/FEM/large_rotation_validation_single_surface.geo

+
+h = 1;
+H = 2;
+
+n = 10;
+
+Point(1) = {0, 0, 0, 0.5};
+Point(2) = {(H-h), 0, 0, 0.5};
+Point(3) = {H-h, h, 0, 0.5};
+Point(4) = {0, h, 0, 0.5};
+
+Point(5) = {H, 0, 0, 0.5};
+Point(6) = {H, h, 0, 0.5};
+
+Point(7) = {H, H, 0, 0.5};
+Point(8) = {H-h, H, 0, 0.5};
+
+Line(1) = {1, 2};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line(5) = {2, 5};
+Line(6) = {5, 6};
+
+Line(8) = {6, 7};
+Line(9) = {7, 8};
+Line(10) = {8, 3};
+
+Curve Loop(1) = {1, 5, 6, 8, 9, 10, 3, 4};
+Plane Surface(1) = {1};
+
+Transfinite Curve {1, 3, 8, 10} = (H-h)*n+1 Using Progression 1;
+Transfinite Curve {4, 5, 6, 9} = n+1 Using Progression 1;
+
+Recombine Surface {1};
+
+Physical Curve("left_edge", 1) = {4};
+Physical Surface("FEM_domain", 2) = {1}; // the trick is to include both plane surfaces in one single domain
+Physical Curve("top_edge", 3) = {9};
+
+F = 40000;
+
+// 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/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Poisson", 0.3);
+SetNumber("Materials/FEM_domain/rho",7800); //volumic mass of acier
+SetNumber("Boundary Conditions/top_edge/tx", F); // 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("Volumic Forces/FEM_domain/bx",0.);
+SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
+
+Physical Curve("BEM_FEM_boundary", 4) = {4};//+
+
+Mesh.Algorithm = 8;

srcs/FEM/mainFEM.cpp

@@ -4,17 +4,18 @@
 #include <iostream>
 #include <omp.h>
 
+// this program implements a FEM solver, either linear or non-linear (by taking large displacements into account)
+// to choose the linear solver, please add: SetNumber("Non_linear_solver",0); in the .geo file.
 int main(int argc, char **argv)
 {
-
-    bool nonLinearSolver = true; // use non-linear solver or not
-
     if (argc < 2)
     {
         std::cout << "Usage: " << argv[0] << " <geo_file>\n";
         return 0;
     }
 
+    double start = omp_get_wtime();
+
     // If compiled with OpenMP support, gmsh::initialize 
     // also sets the number of threads to "General.NumThreads".
     int nthreads = omp_get_max_threads();
@@ -26,22 +27,32 @@ int main(int argc, char **argv)
 
     Eigen::initParallel();
 
+    // determine if non linear solver must be set
+    bool nonLinearSolver = getNonLinearParameter(); 
+
+    // maps the elementTags to a corresponding electrostaticPressure (only useful for coupled solver)
     std::map<int, double> electrostaticPressure;
 
+    int nbViews = 0; // number of views passed to the BEM solver (only useful for coupled solver)
+
     if(nonLinearSolver)
     {
-        std::map<int,std::pair<double,double>> boundaryDisplacementMap;
+        std::map<int,std::pair<double,double>> boundaryDisplacementMap; //(only useful for coupled solver)
 
-        int viewTagU = gmsh::view::add("u"); // à modifier plus tard
-        int viewTagF = gmsh::view::add("f");
+        int viewTagU = gmsh::view::add("u"); //viewtag of the displacement view
+        int viewTagF = gmsh::view::add("f"); //viewtag of the nodal forces view
 
-        solverFEMnonLinear(electrostaticPressure, boundaryDisplacementMap, 0, true, 1, viewTagU, viewTagF); //iteration randomly set to 1
+        solverFEMnonLinear(electrostaticPressure, boundaryDisplacementMap, nbViews, true, 1, viewTagU, viewTagF, false);
     }
     else
     {
-        solverFEM(electrostaticPressure, 0);
+        solverFEM(electrostaticPressure, nbViews);
     }
 
+    double total_time = omp_get_wtime() - start;
+    std::cout << "-----------------------------------------\n";
+    std::cout << "total execution time: " << total_time << " [s]. \n";
+
     gmsh::fltk::run();
     
     gmsh::finalize();

srcs/FEM/self_weight.geo

@@ -37,4 +37,6 @@ SetNumber("Boundary Conditions/right_edge/ty", 0.);
 SetNumber("Volumic Forces/FEM_domain/bx",0.);
 SetNumber("Volumic Forces/FEM_domain/by",-9.81); //set to -9.81 for gravity
 
-Physical Curve("BEM_FEM_boundary", 9) = {1,2,3}; // useless but necessary to make it work
\ No newline at end of file
+Physical Curve("BEM_FEM_boundary", 9) = {1,2,3}; // useless but necessary to make it work
+
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file

srcs/FEM/simple_tension.geo

 Lx = 5;
 Ly = 2;
-nx = 5;
-ny = 2;
+nx = 50;
+ny = 20;
 
 Point(1) = {0, 0, 0, 1.0};
 Point(2) = {Lx, 0, 0, 1.0};
@@ -39,6 +39,8 @@ SetNumber("Boundary Conditions/right_edge/tx", 21e3); // for simple tension cond
 SetNumber("Boundary Conditions/right_edge/ty", 0.);
 SetNumber("Volumic Forces/FEM_domain/bx",0.);
 SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
-SetNumber("Boundary Conditions/fixed_node/uy",0.);
+SetNumber("Boundary Conditions/fixed_node/uy",2.);
 
 Physical Curve("BEM_FEM_boundary", 10) = {1,2,3}; // useless in this case but necessary to make it work
+
+SetNumber("Non_linear_solver",1);
\ No newline at end of file