diff --git a/srcs/beamActuation.geo b/srcs/beamActuation.geo
new file mode 100644
index 0000000000000000000000000000000000000000..a37ee6d35de807819cba90a988a91c0491262946
--- /dev/null
+++ b/srcs/beamActuation.geo
@@ -0,0 +1,177 @@
+scale = 1e-6;
+
+bt = 1; // Beam tickness
+bh = 1; // Beam height (position)
+
+eL = 3; // Electrode base length
+el = 2.5; // Electrode tip length
+et = 1; // Electrode tickness
+es = 1.5; // Electrode spacing
+
+W = 3*es + 3*eL; // Width of the space
+H = bh + 3*(bt+et); // Height of the space
+
+nEs = 9;
+nEL = 12;
+nBt = 6;
+nEt = 3;
+nBh = 3;
+
+bt = bt * scale;
+bh = bh * scale;
+eL = eL * scale;
+el = el * scale;
+et = et * scale;
+es = es * scale;
+W = W * scale;
+H = H * scale;
+
+e = (eL - el)/2;    // Electrode constant
+
+/*
+ *  Points
+ */
+
+// Space
+Point(1) = {0, 0, 0, 1};
+Point(2) = {W, 0, 0, 1};
+Point(3) = {W, H, 0, 1};
+Point(4) = {0, H, 0, 1};
+
+// Beam
+Point(5) = {0, bh+bt, 0, 1};
+Point(6) = {es, bh+bt, 0, 1};
+Point(7) = {es+eL, bh+bt, 0, 1};
+Point(8) = {2*es+eL, bh+bt, 0, 1};
+Point(9) = {2*es+2*eL, bh+bt, 0, 1};
+Point(10) = {3*es+2*eL, bh+bt, 0, 1};
+
+Point(11) = {3*es+2*eL, bh, 0, 1};
+Point(12) = {2*es+2*eL, bh, 0, 1};
+Point(13) = {2*es+eL, bh, 0, 1};
+Point(14) = {es+eL, bh, 0, 1};
+Point(15) = {es, bh, 0, 1};
+Point(16) = {0, bh, 0, 1};
+
+// Electrodes
+Point(17) = {es+e, bh+bt+et, 0, 1};
+Point(18) = {es+e+el, bh+bt+et, 0, 1};
+Point(19) = {2*es+eL+e, bh+bt+et, 0, 1};
+Point(20) = {2*es+eL+e+el, bh+bt+et, 0, 1};
+
+/*
+ *  Lines
+ */
+
+ // Space
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 5};
+Line(5) = {16, 1};
+
+// Beam
+Line(6) = {5, 6};
+Line(7) = {6, 7};
+Line(8) = {7, 8};
+Line(9) = {8, 9};
+Line(10) = {9, 10};
+Line(11) = {10, 11};
+Line(12) = {11, 12};
+Line(13) = {12, 13};
+Line(14) = {13, 14};
+Line(15) = {14, 15};
+Line(16) = {15, 16};
+Line(17) = {16, 5};
+Line(18) = {15, 6};
+Line(19) = {14, 7};
+Line(20) = {13, 8};
+Line(21) = {12, 9};
+
+// Electrode
+Line(22) = {6, 17};
+Line(23) = {17, 18};
+Line(24) = {18, 7};
+Line(25) = {8, 19};
+Line(26) = {19, 20};
+Line(27) = {20, 9};
+
+/*
+ *  Transfinite Curve
+ */
+Transfinite Curve{6, 8, 10, 12, 14, 16} = nEs + 1 Using Progression 1;
+Transfinite Curve{7, 9, 13, 15, 23, 26} = nEL + 1 Using Progression 1;
+Transfinite Curve{17, 18, 19, 20, 21, 11} = nBt + 1 Using Progression 1;
+Transfinite Curve{22, 24, 25, 27} = nEt + 1 Using Progression 1;
+
+Transfinite Curve{1, 3} = 3*nEs + 3*nEL + 1 Using Progression 1;
+Transfinite Curve(2) = 4 * (nBt + nEt) + 1 Using Progression 1;
+Transfinite Curve(4) = 3*nEt + 2*nBt + 1 Using Progression 1;
+Transfinite Curve(5) = nBh + 1 Using Progression 1;
+
+/*
+ *  Curve Loop
+ */
+Curve Loop(1) = {6, -18, 16, 17};
+Curve Loop(2) = {7, -19, 15, 18};
+Curve Loop(3) = {8, -20, 14, 19};
+Curve Loop(4) = {9, -21, 13, 20};
+Curve Loop(5) = {10, 11, 12, 21};
+Curve Loop(6) = {23, 24, -7, 22};
+Curve Loop(7) = {26, 27, -9, 25};
+Curve Loop(8) = {1, 2, 3, 4, 6, 22, 23, 24, 8, 25, 26, 27, 10, 11, 12, 13, 14, 15, 16, 5};
+
+/*
+ *  Plane Surface
+ */
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {3};
+Plane Surface(4) = {4};
+Plane Surface(5) = {5};
+Plane Surface(6) = {6};
+Plane Surface(7) = {7};
+Plane Surface(8) = {8};
+
+/*
+ *  Surface Manipulation
+ */
+Transfinite Surface{1, 2, 3, 4, 5, 6, 7};
+Recombine Surface{1, 2, 3, 4, 5, 6, 7};
+
+/*
+ *  Physical Curve
+ */
+Physical Curve("clamp", 1) = {17};
+Physical Curve("outside", 2) = {1, 2, 3, 4, 5};
+Physical Curve("left_electrode", 3) = {24};
+Physical Curve("right_electrode", 4) = {25};
+Physical Curve("beam", 5) = {6, 22, 23, 8, 26, 10, 11, 12, 13, 14, 15, 16};
+Physical Curve("BEM_FEM_boundary", 6) = {6, 22, 23, 24, 8, 26, 26, 27, 10, 11, 12, 13, 14, 15, 16};
+
+/*
+ *  Physical Surface
+ */
+Physical Surface("FEM_domain", 1) = {1, 2, 3, 4, 5, 6, 7};
+Physical Surface("BEM_domain_1", 2) = {8};
+
+/*
+ *  Parameters
+ */
+SetNumber("Boundary Conditions/clamp/ux", 0.0);
+SetNumber("Boundary Conditions/clamp/uy", 0.0);
+
+SetNumber("Materials/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Poisson", 0.3);
+SetNumber("Materials/FEM_domain/rho",7800);
+
+SetNumber("Volumic Forces/FEM_domain/bx", 0.0);
+SetNumber("Volumic Forces/FEM_domain/by", 0.0);
+
+SetNumber("Boundary Conditions/outside/BEM_domain_1/neumann", 0);
+SetNumber("Boundary Conditions/beam/BEM_domain_1/neumann", 0);
+SetNumber("Boundary Conditions/left_electrode/BEM_domain_1/dirichlet", 200);
+SetNumber("Boundary Conditions/right_electrode/BEM_domain_1/dirichlet", 0);
+SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
+
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/beamActuationGeneral.geo b/srcs/beamActuationGeneral.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3ad6cf8a80d55447284b1467115cedc983abd0e0
--- /dev/null
+++ b/srcs/beamActuationGeneral.geo
@@ -0,0 +1,229 @@
+scale = 1e-6;
+
+bt = 0.2; // Beam tickness
+bh = 1; // Beam height (position)
+
+eL = 3; // Electrode base length
+el = 1.5; // Electrode tip length
+et = 1; // Electrode tickness
+es = 1.5; // Electrode spacing
+
+nbElectrodes = 10;
+
+W = nbElectrodes*(es + eL) + 2*es;    // Width of the space
+H = W;  // Height of the space
+
+nEs = 12;
+nEL = 15;
+nBt = 8;
+nEt = 3;
+nBh = 3;
+
+bt = bt * scale;
+bh = bh * scale;
+eL = eL * scale;
+el = el * scale;
+et = et * scale;
+es = es * scale;
+W = W * scale;
+H = H * scale;
+
+e = (eL - el)/2;    // Electrode constant
+
+/*
+ *  Points
+ */
+
+For i In {0:nbElectrodes-1}
+    Point(6*i + 1) = {(es+eL)*i, bh + bt, 0, 1};
+    Point(6*i + 2) = {(es+eL)*i + es, bh + bt, 0, 1};
+    Point(6*i + 3) = {(es+eL)*i + es+e, bh + bt + et, 0, 1};
+    Point(6*i + 4) = {(es+eL)*i + es+e+el, bh + bt + et, 0, 1};
+    Point(6*i + 5) = {(es+eL)*i + es, bh, 0, 1};
+    Point(6*i + 6) = {(es+eL)*i, bh, 0, 1};
+EndFor
+Point(6*nbElectrodes + 1) = {nbElectrodes*(es+eL), bh + bt, 0, 1};
+Point(6*nbElectrodes + 2) = {nbElectrodes*(es+eL) + es, bh + bt, 0, 1};
+Point(6*nbElectrodes + 3) = {nbElectrodes*(es+eL) + es, bh, 0, 1};
+Point(6*nbElectrodes + 6) = {nbElectrodes*(es+eL), bh, 0, 1};
+
+Point(6*nbElectrodes + 4) = {0, 0, 0, 1};
+Point(6*nbElectrodes + 5) = {W, 0, 0, 1};
+Point(6*nbElectrodes + 7) = {W, H, 0, 1};
+Point(6*nbElectrodes + 8) = {0, H, 0, 1};
+
+For i In {0:nbElectrodes-1}
+    Point(6*nbElectrodes + 8 + i + 1) = {i*(es+eL) + es + e + el/2, bh + bt + et - ((el*e)/(2*et)), 0, 1};
+EndFor
+
+/*
+ *  Lines
+ */
+Line(9*nbElectrodes + 1) = {6*nbElectrodes + 4, 6*nbElectrodes + 5};
+Line(9*nbElectrodes + 2) = {6*nbElectrodes + 5, 6*nbElectrodes + 7};
+Line(9*nbElectrodes + 3) = {6*nbElectrodes + 7, 6*nbElectrodes + 8};
+Line(9*nbElectrodes + 4) = {6*nbElectrodes + 8, 1};
+Line(9*nbElectrodes + 5) = {6, 6*nbElectrodes + 4};
+
+Line(9*nbElectrodes + 9) = {6*nbElectrodes + 1, 6*nbElectrodes + 2};
+Line(9*nbElectrodes + 6) = {6*nbElectrodes + 2, 6*nbElectrodes + 3};
+Line(9*nbElectrodes + 7) = {6*nbElectrodes + 3, 6*nbElectrodes + 6};
+Line(9*nbElectrodes + 8) = {6*nbElectrodes + 6, 6*nbElectrodes + 1};
+For i In {0:nbElectrodes-1}
+    Line(9*i + 1) = {6*i + 1, 6*i + 2};
+    Line(9*i + 2) = {6*i + 2, 6*i + 3};
+    // Line(9*i + 3) = {6*i + 3, 6*i + 4};
+    Circle(9*i + 3) = {6*i + 3, 6*nbElectrodes + 8 + i + 1, 6*i + 4};
+    Line(9*i + 4) = {6*i + 4, 6*(i+1) + 1};
+    Line(9*i + 5) = {6*(i+1) + 1, 6*i + 2};
+    Line(9*i + 6) = {6*i + 5, 6*i + 2};
+    Line(9*i + 7) = {6*i + 5, 6*i + 6};
+    Line(9*i + 8) = {6*i + 6, 6*i + 1};
+    Line(9*i + 9) = {6*(i+1) + 6, 6*i + 5};
+EndFor
+
+/*
+ *  Transfinite Curve
+ */
+For i In {0:nbElectrodes-1}
+    Transfinite Curve{9*i + 1, 9*i + 7} = nEs + 1 Using Progression 1;
+    Transfinite Curve{9*i + 2, 9*i + 4} = nEt + 1 Using Progression 1;
+    Transfinite Curve{9*i + 3, 9*i + 5, 9*i + 9} = nEL + 1 Using Progression 1;
+    Transfinite Curve{9*i + 8, 9*i + 6} = nBt + 1 Using Progression 1;
+EndFor
+Transfinite Curve{9*nbElectrodes + 9, 9*nbElectrodes + 7} = nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 8, 9*nbElectrodes + 6} = nBt + 1 Using Progression 1;
+
+Transfinite Curve{9*nbElectrodes + 1, 9*nbElectrodes + 3} = nbElectrodes*(nEs + nEL) + 2*nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 2} = nbElectrodes*(nEs + nEL) + 2*nEs + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 4} = nbElectrodes*(nEs + nEL) + 2*nEs - nBt - nBh + 1 Using Progression 1;
+Transfinite Curve{9*nbElectrodes + 5} = nBh + 1 Using Progression 1;
+
+/*
+ *  Curve Loop
+ */
+
+For i In {0:nbElectrodes-1}
+    Curve Loop(3*i + 1) = {9*i + 1, -(9*i + 6), 9*i + 7, 9*i + 8};
+    Curve Loop(3*i + 2) = {9*i + 3, 9*i + 4, 9*i + 5, 9*i + 2};
+    Curve Loop(3*i + 3) = {-(9*i + 5), -(9*(i+1) + 8), 9*i + 9, 9*i + 6};
+EndFor
+Curve Loop(3*nbElectrodes + 1) = {9*nbElectrodes + 9, 9*nbElectrodes + 6, 9*nbElectrodes + 7, 9*nbElectrodes + 8};
+// How to fill a vector in a loop ?
+
+bemBoundary = {1:6*nbElectrodes + 8};
+bemBoundaryWithoutE = {1:4*nbElectrodes + 8 + 2};
+bemFemBoundary = {1:6*nbElectrodes + 3};
+leftElectrodesFull = {1:nbElectrodes};
+rightElectrodesFull = {1:nbElectrodes};
+For i In {0:nbElectrodes-1}
+    bemBoundary[4*i+1 - 1] = 9*i + 1;
+    bemBoundary[4*i+2 - 1] = 9*i + 2;
+    bemBoundary[4*i+3 - 1] = 9*i + 3;
+    bemBoundary[4*i+4 - 1] = 9*i + 4;
+
+    bemBoundaryWithoutE[2*i+1 - 1] = 9*i + 1;
+    bemBoundaryWithoutE[2*i+2 - 1] = 9*i + 3;
+
+    bemFemBoundary[4*i+1 - 1] = 9*i + 1;
+    bemFemBoundary[4*i+2 - 1] = 9*i + 2;
+    bemFemBoundary[4*i+3 - 1] = 9*i + 3;
+    bemFemBoundary[4*i+4 - 1] = 9*i + 4;
+
+    leftElectrodesFull[i] = 9*i + 2;
+    rightElectrodesFull[i] = 9*i + 4;
+EndFor
+bemBoundary[4*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemBoundary[4*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemBoundary[4*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+
+bemBoundaryWithoutE[2*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemBoundaryWithoutE[2*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemBoundaryWithoutE[2*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+
+bemFemBoundary[4*nbElectrodes + 1 - 1] = 9*nbElectrodes + 9;
+bemFemBoundary[4*nbElectrodes + 2 - 1] = 9*nbElectrodes + 6;
+bemFemBoundary[4*nbElectrodes + 3 - 1] = 9*nbElectrodes + 7;
+For i In {0:nbElectrodes-1}
+    bemBoundary[4*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+    bemBoundary[4*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+
+    bemBoundaryWithoutE[2*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+    bemBoundaryWithoutE[2*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+
+    bemFemBoundary[4*nbElectrodes + 3 + 2*i+1 - 1] = 9*i + 9;
+    bemFemBoundary[4*nbElectrodes + 3 + 2*i+2 - 1] = 9*i + 7;
+EndFor
+bemBoundary[6*nbElectrodes + 4 - 1] = 9*nbElectrodes + 5;
+bemBoundary[6*nbElectrodes + 5 - 1] = 9*nbElectrodes + 1;
+bemBoundary[6*nbElectrodes + 6 - 1] = 9*nbElectrodes + 2;
+bemBoundary[6*nbElectrodes + 7 - 1] = 9*nbElectrodes + 3;
+bemBoundary[6*nbElectrodes + 8 - 1] = 9*nbElectrodes + 4;
+
+bemBoundaryWithoutE[4*nbElectrodes + 4 - 1] = 9*nbElectrodes + 5;
+bemBoundaryWithoutE[4*nbElectrodes + 5 - 1] = 9*nbElectrodes + 1;
+bemBoundaryWithoutE[4*nbElectrodes + 6 - 1] = 9*nbElectrodes + 2;
+bemBoundaryWithoutE[4*nbElectrodes + 7 - 1] = 9*nbElectrodes + 3;
+bemBoundaryWithoutE[4*nbElectrodes + 8 - 1] = 9*nbElectrodes + 4;
+
+leftElectrodes = {1:nbElectrodes-1};
+rightElectrodes = {1:nbElectrodes-1};
+For i In {0:nbElectrodes-2}
+    leftElectrodes[i] = leftElectrodesFull[i+1];
+    rightElectrodes[i] = rightElectrodesFull[i];
+EndFor
+
+bemBoundaryWithoutE[4*nbElectrodes + 9 - 1] = 2;
+bemBoundaryWithoutE[4*nbElectrodes + 10 - 1] = 9*nbElectrodes + 4;
+
+Curve Loop(3*nbElectrodes + 2) = bemBoundary[];
+
+/*
+ *  Plane Surface
+ */
+For i In {0:nbElectrodes-1}
+    Plane Surface(3*i + 1) = {3*i + 1};
+    Plane Surface(3*i + 2) = {3*i + 2};
+    Plane Surface(3*i + 3) = {3*i + 3};
+EndFor
+Plane Surface(3*nbElectrodes + 1) = {3*nbElectrodes + 1};
+Plane Surface(3*nbElectrodes + 2) = {3*nbElectrodes + 2};
+
+/*
+ *  Surface operations
+ */
+Recombine Surface{1:3*nbElectrodes + 1};
+Transfinite Surface{1:3*nbElectrodes + 1};
+
+/*
+ *  Physical Curve
+ */
+Physical Curve("clamp", 1) = {8};
+Physical Curve("BEM_boundary_without_e", 2) = bemBoundaryWithoutE[];
+Physical Curve("BEM_FEM_boundary", 3) = bemFemBoundary[];
+Physical Curve("left_electrodes", 4) = leftElectrodes[];
+Physical Curve("right_electrodes", 5) = rightElectrodes[];
+
+Physical Surface("FEM_domain") = {1:3*nbElectrodes + 1};
+Physical Surface("BEM_domain_1") = {3*nbElectrodes + 2};
+
+/*
+ *  Parameters
+ */
+ SetNumber("Boundary Conditions/clamp/ux", 0.0);
+ SetNumber("Boundary Conditions/clamp/uy", 0.0);
+ 
+//  SetNumber("Materials/FEM_domain/Young", 3e7);
+SetNumber("Materials/FEM_domain/Young", 1e7);
+ SetNumber("Materials/FEM_domain/Poisson", 0.3);
+ SetNumber("Materials/FEM_domain/rho",7800);
+ 
+ SetNumber("Volumic Forces/FEM_domain/bx", 0.0);
+ SetNumber("Volumic Forces/FEM_domain/by", 0.0);
+ 
+ SetNumber("Boundary Conditions/BEM_boundary_without_e/BEM_domain_1/neumann", 0);
+ SetNumber("Boundary Conditions/left_electrodes/BEM_domain_1/dirichlet", 8);
+ SetNumber("Boundary Conditions/right_electrodes/BEM_domain_1/dirichlet", 0);
+ SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
+ 
+ SetNumber("Non_linear_solver", 1);
\ No newline at end of file