diff --git a/README.md b/README.md
index ae94909629b33953ad17c02389d8c63533c68bf6..032b488f310bece60979e8553ec36ac2bd4fd760 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,7 @@
# Multiphysics 0471
+The code has been written by Louis Denis, Pierre-Yves Goffin, Kévin Maltez and Valentin Vanraes in the context of the Multiphysics Integrated Computational Project (Spring 2022), as part of the Master of Science in Engineering Physics of the Faculty of Applied Sciences at University of Liège, Belgium.
-Multiphysics Integrated Computational Project
+## Building the solver
FOR WINDOWS:
@@ -67,7 +68,7 @@ To create a .geo file compatible with the program, some rules must be respected:
- **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);`.
+ 2. Neumann boundary conditions, corresponding to imposing the opposite of the exterior 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.
diff --git a/srcs/BEM/mainBEM.cpp b/srcs/BEM/mainBEM.cpp
index cd639bbb7d0f038b5fcbc83f760d2f1b87ee232f..0a4e31a8bb9272f5de47b6282f7617e6888b8155 100644
--- a/srcs/BEM/mainBEM.cpp
+++ b/srcs/BEM/mainBEM.cpp
@@ -15,6 +15,7 @@ int main(int argc, char **argv)
omp_set_num_threads(nthreads);
gmsh::initialize();
+ omp_set_num_threads(nthreads);
gmsh::open(argv[1]);
gmsh::model::mesh::generate(2);
diff --git a/srcs/BEM/postProcessing.cpp b/srcs/BEM/postProcessing.cpp
index 691daac77f047786ed3cf5140f35429bf0bba996..45919f256dcb71344900469bc0dd09e5b3b704c6 100644
--- a/srcs/BEM/postProcessing.cpp
+++ b/srcs/BEM/postProcessing.cpp
@@ -324,8 +324,11 @@ void boundaryElementTags(std::vector<std::size_t> &domainBoundaryElementTags,
// #pragma omp for schedule(guided)
for(std::size_t i = 0; i < elementVector.size(); i++)
{
- std::vector<size_t> elementTags;
- gmsh::model::mesh::getElementsByCoordinates(elementVector[i].x1, elementVector[i].y1, 0, elementTags, 2);
+ int elementType;
+ std::vector<size_t> nodeTags;
+ int entityDim;
+ int entityTag;
+ std::vector<int> physicalTags;
// Adds the elements tags to the local vector.
for(std::size_t j = 0; j < elementTags.size(); j++)
diff --git a/srcs/FEM/complex_validation.geo b/srcs/FEM/complex_validation.geo
index c6497d704d992ca7c8e18d45c13bd1a1477d9b89..9d0cd28d57118abc3035c2ade7597c6903f2f4f7 100644
--- a/srcs/FEM/complex_validation.geo
+++ b/srcs/FEM/complex_validation.geo
@@ -30,7 +30,7 @@ Plane Surface(2) = {2};
//Transfinite Surface {2};
-Mesh.ElementOrder = 2;
+Mesh.ElementOrder = 1;
Recombine Surface {1}; // quads instead of triangles
@@ -54,6 +54,4 @@ 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",2);
-Physical Curve("BEM_FEM_boundary", 10) = {1,2,3,4,5,6};
-
-SetNumber("Non_linear_solver", 0);
\ No newline at end of file
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/FEM/hybrid_complex.geo b/srcs/FEM/hybrid_complex.geo
index 97260c6e00c6fcd4f876e3402157891d1611c5b9..43fec41702f91afebca03d8b50122f7ebf6603a3 100644
--- a/srcs/FEM/hybrid_complex.geo
+++ b/srcs/FEM/hybrid_complex.geo
@@ -3,9 +3,11 @@ Ly_poutre = 2;
h_poutre = 2;
h_tot = 10;
width = 10;
-nx = 10; // prend beaucoup de temps à pd de 200x40
+nx = 10;
ny = 4;
+// THIS .geo FILE IS NOT USED IN THE REPORT.
+
Point(1) = {0, h_poutre, 0, 2};
Point(2) = {Lx_poutre, h_poutre, 0, 2};
Point(3) = {Lx_poutre, h_poutre + Ly_poutre, 0, 2};
@@ -18,7 +20,7 @@ Point(8) = {0, h_tot, 0, 2};
Point(9) = {Lx_poutre/2, h_poutre, 0, 2};
Point(10) = {Lx_poutre/2, h_poutre + Ly_poutre, 0, 2};
-// A IMPLEMENTER LA SUITE
+
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
@@ -60,17 +62,15 @@ Physical Point("fixed_node", 9) = {1};
Physical Surface("BEM_domain", 10) = {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/ux", 0.);
//SetNumber("Boundary Conditions/left_edge/uy", 2.);
SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_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("Materials/FEM_domain/rho",7800);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", 0.);
+SetNumber("Boundary Conditions/right_edge/tx", 21e3);
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",2.);
-
-Physical Curve("BEM_FEM_boundary", 11) = {1,2,3}; // useless but necessary to make it work
diff --git a/srcs/FEM/hybrid_geo.geo b/srcs/FEM/hybrid_geo.geo
index 4c6fb60c284587664a2c48d52c0ac044ff742854..ef301636dcbafdc37568b5c226f5b7f27225c0be 100644
--- a/srcs/FEM/hybrid_geo.geo
+++ b/srcs/FEM/hybrid_geo.geo
@@ -6,6 +6,8 @@ width = 10;
nx = 20; // prend beaucoup de temps à pd de 200x40
ny = 8;
+// THIS .geo FILE IS NOT USED IN THE REPORT.
+
Point(1) = {0, h_poutre, 0, 2};
Point(2) = {Lx_poutre, h_poutre, 0, 2};
Point(3) = {Lx_poutre, h_poutre + Ly_poutre, 0, 2};
@@ -58,17 +60,17 @@ Physical Point("fixed_node", 9) = {1};
Physical Surface("BEM_domain", 10) = {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/ux", 0.);
//SetNumber("Boundary Conditions/left_edge/uy", 2.);
SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_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("Materials/FEM_domain/rho",7800);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", 0.);
+SetNumber("Boundary Conditions/right_edge/tx", 21e3);
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("Volumic Forces/FEM_domain/by",0.);
SetNumber("Boundary Conditions/fixed_node/uy",2.);
// BEM geometry
diff --git a/srcs/FEM/large_rotation_validation.geo b/srcs/FEM/large_rotation_validation.geo
index 113db4e8a52d2345b5b8075dc49a0b23ff2dab77..2386bc5e882b643b509b866431ef7c3edd19b1a1 100644
--- a/srcs/FEM/large_rotation_validation.geo
+++ b/srcs/FEM/large_rotation_validation.geo
@@ -2,7 +2,7 @@
h = 1;
H = 10;
-n = 16;
+n = 4;
Point(1) = {0, 0, 0, 0.5};
Point(2) = {0.9*H, 0, 0, 0.5};
@@ -83,6 +83,4 @@ SetNumber("Boundary Conditions/top_edge/ty", 0); //set to some other value for v
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
+SetNumber("Non_linear_solver",1); // activate non linear solver
\ No newline at end of file
diff --git a/srcs/FEM/large_rotation_validation_single_surface.geo b/srcs/FEM/large_rotation_validation_single_surface.geo
index d15291b19f3553b05de6fcccd5bd78e967e6b1a7..6db9cbc05ed3dcf33cc393fd15a5794c509cb3fd 100644
--- a/srcs/FEM/large_rotation_validation_single_surface.geo
+++ b/srcs/FEM/large_rotation_validation_single_surface.geo
@@ -4,6 +4,8 @@ H = 2;
n = 10;
+// THIS .geo FILE IS NOT USED IN THE REPORT.
+
Point(1) = {0, 0, 0, 0.5};
Point(2) = {(H-h), 0, 0, 0.5};
Point(3) = {H-h, h, 0, 0.5};
diff --git a/srcs/FEM/self_weight.geo b/srcs/FEM/self_weight.geo
index c04a6a8a44f4dd0dc62137b00f36ea44a94cd666..601c9eb4549efded7daab671b919b87f66fc85c1 100644
--- a/srcs/FEM/self_weight.geo
+++ b/srcs/FEM/self_weight.geo
@@ -3,6 +3,8 @@ Ly = 2;
nx = 50;
ny = 20;
+// THIS .geo FILE IS NOT USED IN THE REPORT
+
Point(1) = {0, 0, 0, 1.0};
Point(2) = {Lx, 0, 0, 1.0};
Point(3) = {Lx, Ly, 0, 1.0};
@@ -37,6 +39,4 @@ 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
-
SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/FEM/simple_tension.geo b/srcs/FEM/simple_tension.geo
index f1452e518608936d6c4b17101e93a32c0f296abd..81de6b72f9d36524324c9f78553b5eda0093e814 100644
--- a/srcs/FEM/simple_tension.geo
+++ b/srcs/FEM/simple_tension.geo
@@ -28,19 +28,17 @@ 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/ux", 0.);
//SetNumber("Boundary Conditions/left_edge/uy", 0.);
SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_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("Materials/FEM_domain/rho",7800);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", 0.);
SetNumber("Boundary Conditions/right_edge/tx", 21e3); // for simple tension conditions
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("Volumic Forces/FEM_domain/by",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
+SetNumber("Non_linear_solver",0);
\ No newline at end of file
diff --git a/srcs/FEM/solverFEM.cpp b/srcs/FEM/solverFEM.cpp
index 97bcda367a82fab7ef26824712b9e056fac5069f..088dc6aa5bec43af8f5c496b7e69c91e4378d368 100644
--- a/srcs/FEM/solverFEM.cpp
+++ b/srcs/FEM/solverFEM.cpp
@@ -190,7 +190,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure,
previousTotalNodalForces[i] = 0;
currentTotalNodalForces[i] = 0;
}
- double residualTolerance = 1e-12; // empirical
+ double residualTolerance = 1e-10; // empirical
double start5 = omp_get_wtime();
step_name = "further initialization of global variables for the iterative algorithm";
@@ -250,7 +250,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure,
relativeForces = sqrt(relativeForces/(2*FEM_nodeTags.size())/currentMaxNodalForce); // norm made relative
else //to avoid division by zero
relativeForces = 1e17;
- std::cout << "step: " << step << ", relativeForces: " << relativeForces << "\n";
+ //std::cout << "step: " << step << ", relativeForces: " << relativeForces << "\n";
if(previousRelativeForces < relativeForces) //the current step induces an increase of nodal difference
{
@@ -434,7 +434,10 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure,
/*--------------RETRIEVING NODAL DISPLACEMENTS ON BEM-FEM BOUNDARY------------*/
// those displacements will be sent to the BEM code in the non-linear iterative solver
// mapping nodeTag (of node on the boundary) into (u,v) displacements of the node
- retrieveBoundaryDisplacements(boundaryDisplacementMap, full_d_vector, nodeIndexMap, groups);
+ if(electrostaticPressure.size())
+ {
+ retrieveBoundaryDisplacements(boundaryDisplacementMap, full_d_vector, nodeIndexMap, groups);
+ }
// if mesh untangler is activated, all nodal coordinates are updated according to their final
// nodal displacements. It is useful when coupled to the BEM solver in order to untangle
diff --git a/srcs/FEM/uniform_charge.geo b/srcs/FEM/uniform_charge.geo
index 930e4410c83958dd777d5c1f19e42af13acb4e13..9baa462b391fb24b45fc7018edc0523c63df3a3f 100644
--- a/srcs/FEM/uniform_charge.geo
+++ b/srcs/FEM/uniform_charge.geo
@@ -1,7 +1,6 @@
Lx = 10;
Ly = 1;
-nx = 200;
-ny = 20;
+n = 64;
Point(1) = {0, 0, 0, 1.0};
Point(2) = {Lx, 0, 0, 1.0};
@@ -13,8 +12,8 @@ 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 Curve {3, 1} = Lx*n+1 Using Progression 1;
+Transfinite Curve {4, 2} = Ly*n+1 Using Progression 1;
Transfinite Surface {1};
Recombine Surface {1}; // quads instead of triangles
@@ -24,19 +23,19 @@ Physical Surface("FEM_domain", 6) = {1};
Physical Curve("top_edge", 7) = {3};
Physical Curve("right_edge", 8) = {2};
+//Mesh.ElementOrder = 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/ux", 0.);
SetNumber("Boundary Conditions/left_edge/uy", 0.);
SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_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("Materials/FEM_domain/rho",7800);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", -1);
+SetNumber("Boundary Conditions/right_edge/tx", 0.);
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
-
-Physical Curve("BEM_FEM_boundary", 9) = {1,2,3};
+SetNumber("Volumic Forces/FEM_domain/by",0.);
-SetNumber("Non_linear_solver", 1);
\ No newline at end of file
+SetNumber("Non_linear_solver", 0);
\ No newline at end of file
diff --git a/srcs/FoldedFlexureBeam.geo b/srcs/FoldedFlexureBeam.geo
index 1d71824d2c4a7997b4f1c1ae6c3182f5c3a7636b..87a395b08a9a518b072e5afa992ce5fa1ff85957 100644
--- a/srcs/FoldedFlexureBeam.geo
+++ b/srcs/FoldedFlexureBeam.geo
@@ -1,10 +1,10 @@
scale = 1e-6;
-nBEM = 10; // BEM element density (mettre ça dans quatrième coordonnée)
+nFEM = 2; // FEM element density
-nFEM = 0.5; // FEM element density
+nBEM = 5/nFEM; // BEM element density (mettre ça dans quatrième coordonnée)
-phi = 100;
+phi = 110;
SetNumber("Boundary Conditions/BEM_FEM_boundary/BEM_domain_1/dirichlet", 0);
SetNumber("Boundary Conditions/top_electrode/BEM_domain_1/dirichlet", phi);
SetNumber("Boundary Conditions/rest_of_outside/BEM_domain_1/neumann", 0);
@@ -16,7 +16,7 @@ SetNumber("Boundary Conditions/anchor/uy", 0);
SetNumber("Boundary Conditions/left_edge/ux", 0);
SetNumber("Materials/FEM_domain/Young", 160e9); // A DETERMINER PRECISEMENT
SetNumber("Materials/FEM_domain/Poisson", 0.22);
-SetNumber("Materials/FEM_domain/rho",7800); //MODIF
+SetNumber("Materials/FEM_domain/rho",2300);
SetNumber("Volumic Forces/FEM_domain/bx",0);
SetNumber("Volumic Forces/FEM_domain/by",0);
SetNumber("Non_linear_solver",1);
@@ -53,12 +53,12 @@ Ncombs = 4; // ne pas changer, pas généralisé
// fixed bottom anchor
y = -(Lg-2-4)*scale;
x = (5*g + 4.5*Wc)*scale;
-Point(1) = {x, y-Wss-2*scale, 0, nBEM*scale}; // limite avec suite
-Point(2) = {x, y-Wts, 0, nBEM*scale};
-Point(3) = {x-Wts, y-Wts, 0, nBEM*scale};
-Point(4) = {x-Wts, y, 0, nBEM*scale};
-Point(5) = {x, y, 0, nBEM*scale};
-Point(6) = {x, y-2*scale, 0, nBEM*scale}; // limite avec suite
+Point(1) = {x, y-Wss-2*scale, 0, 10*scale}; // limite avec suite
+Point(2) = {x, y-Wts, 0, 10*scale};
+Point(3) = {x-Wts, y-Wts, 0, 10*scale};
+Point(4) = {x-Wts, y, 0, 10*scale};
+Point(5) = {x, y, 0, 10*scale};
+Point(6) = {x, y-2*scale, 0, 10*scale}; // limite avec suite
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
@@ -75,8 +75,8 @@ Transfinite Surface{1} = {2:5};
Recombine Surface{1};
// bottom long beam
-Point(7) = {x+Lss, y-2*scale, 0, nBEM*scale}; // limite avec suite
-Point(8) = {x+Lss, y-Wss-2*scale, 0, nBEM*scale}; // limite avec suite
+Point(7) = {x+Lss, y-2*scale, 0, 10*scale}; // limite avec suite
+Point(8) = {x+Lss, y-Wss-2*scale, 0, 10*scale}; // limite avec suite
Line(7) = {6, 7};
Line(8) = {7, 8}; // limite avec suite
Line(9) = {8, 1};
@@ -90,12 +90,12 @@ Recombine Surface{2};
//truss
x = x + Lss;
y = y - 10*scale - Wss;
-Point(9) = {x, y + 8*scale + Wss + Lgs, 0 , nBEM*scale}; // limite avec suite
-Point(10) = {x, y + 8*scale + 2*Wss + Lgs, 0 , nBEM*scale}; // limite avec suite
-Point(11) = {x, y + 16*scale + 2*Wss + Lgs, 0 , nBEM*scale};
-Point(12) = {x + Wts, y + 16*scale + 2*Wss + Lgs, 0 , nBEM*scale};
-Point(13) = {x + Wts, y, 0 , nBEM*scale};
-Point(14) = {x, y, 0 , nBEM*scale};
+Point(9) = {x, y + 8*scale + Wss + Lgs, 0 , 10*scale}; // limite avec suite
+Point(10) = {x, y + 8*scale + 2*Wss + Lgs, 0 , 10*scale}; // limite avec suite
+Point(11) = {x, y + 16*scale + 2*Wss + Lgs, 0 , 10*scale};
+Point(12) = {x + Wts, y + 16*scale + 2*Wss + Lgs, 0 , 10*scale};
+Point(13) = {x + Wts, y, 0 , 10*scale};
+Point(14) = {x, y, 0 , 10*scale};
Line(10) = {7, 9};
Line(11) = {9, 10}; // limite avec suite
Line(12) = {10, 11};
@@ -115,8 +115,8 @@ Recombine Surface{3};
//top long beam
x = (5*g + 4.5*Wc)*scale;
-Point(15) = {x, 4*scale, 0 , nBEM*scale}; // limite avec suite
-Point(16) = {x, 4*scale + Wss, 0 , nBEM*scale}; // limite avec suite
+Point(15) = {x, 4*scale, 0 , 10*scale}; // limite avec suite
+Point(16) = {x, 4*scale + Wss, 0 , 10*scale}; // limite avec suite
Line(17) = {9, 15};
Line(18) = {15, 16}; // limite avec suite
Line(19) = {16, 10};
@@ -128,10 +128,10 @@ Transfinite Surface{4};
Recombine Surface{4};
//bottom base moveable combs
-Point(17) = {x, 0, 0 , nBEM*scale};
-Point(18) = {0, 0, 0 , nBEM*scale}; // roulements
-Point(19) = {0, 14*scale + Wss, 0 , nBEM*scale}; // roulements + limite avec suite
-Point(20) = {x, 14*scale + Wss, 0 , nBEM*scale}; // limite avec suite
+Point(17) = {x, 0, 0 , 10*scale};
+Point(18) = {0, 0, 0 , 10*scale}; // roulements
+Point(19) = {0, 14*scale + Wss, 0 , 10*scale}; // roulements + limite avec suite
+Point(20) = {x, 14*scale + Wss, 0 , 10*scale}; // limite avec suite
Line(20) = {15, 17};
Line(21) = {17, 18};
Line(22) = {18, 19}; // roulements
@@ -245,14 +245,14 @@ Recombine Surface{7:10};
Physical Curve("anchor", 1) = {1:6};
Physical Surface("FEM_domain", 2) = {1:10};
-Physical Curve("BEM_FEM_boundary", 3) = {1:5, 7, 10, 17, 20, 21, 26:nbLines2, 24, 19, 12:16, 9};
+Physical Curve("BEM_FEM_boundary", 3) = {26:nbLines2, 24, 19, 12:13};
Physical Curve("left_edge", 4) = {22,25};
//top electrode
y = y + Lcs + (Lcs-y0s) + 12*scale;
xtot = Ncombs*(2*gs+2*Wcs) + Wcs/2;
-Point(nbPts2 + 1) = {0, y, 0, nBEM*scale};
-Point(nbPts2 + 2) = {xtot, y, 0, nBEM*scale};
+Point(nbPts2 + 1) = {0, y, 0, 10*scale};
+Point(nbPts2 + 2) = {xtot, y, 0, 10*scale};
For i In {1:Ncombs}
Point(nbPts2 + 2 + 4*i - 3) = {xtot-(i-1)*(2*Wcs+2*gs), y-12*scale-Lcs, 0, 0.2*nBEM*scale};
Point(nbPts2 + 2 + 4*i - 2) = {xtot-(i-1)*(2*Wcs+2*gs) - Wcs, y-12*scale - Lcs, 0, 0.2*nBEM*scale};
@@ -260,8 +260,8 @@ For i In {1:Ncombs}
Point(nbPts2 + 2 + 4*i) = {xtot-i*(2*Wcs+2*gs), y-12*scale, 0, 0.2*nBEM*scale};
EndFor
nbPts3 = nbPts2 + 2 + 4*Ncombs;
-Point(nbPts3 + 1) = {Wcs/2, y-12*scale-Lcs, 0, nBEM*scale};
-Point(nbPts3 + 2) = {0, y-12*scale-Lcs, 0, nBEM*scale};
+Point(nbPts3 + 1) = {Wcs/2, y-12*scale-Lcs, 0, 0.2*nBEM*scale};
+Point(nbPts3 + 2) = {0, y-12*scale-Lcs, 0, 0.2*nBEM*scale};
nbPts4 = nbPts3 + 2;
For i In {nbPts2+1:nbPts4-1}
@@ -274,20 +274,18 @@ Physical Curve("top_electrode", 5) = {nbLines3+1:nbLines4-1};
// rest of outside
Line(nbLines4 + 1) = {nbPts4, 21};
+y = -(Lg-2-4)*scale - 10*scale - Wss + 16*scale + 2*Wss + Lgs;
+Point(nbPts4 + 1) = {Wtots, y, 0, 20*scale};
+Point(nbPts4 + 2) = {Wtots, Htops, 0, 20*scale};
+Point(nbPts4 + 3) = {0, Htops, 0, 20*scale};
-Point(nbPts4 + 1) = {0, -Hbottoms, 0, 2*nBEM*scale};
-Point(nbPts4 + 2) = {Wtots, -Hbottoms, 0, 2*nBEM*scale};
-Point(nbPts4 + 3) = {Wtots, Htops, 0, 2*nBEM*scale};
-Point(nbPts4 + 4) = {0, Htops, 0, 2*nBEM*scale};
-
-Line(nbLines4 + 2) = {18, nbPts4 + 1};
+Line(nbLines4 + 2) = {12, nbPts4 + 1};
Line(nbLines4 + 3) = {nbPts4 + 1, nbPts4 + 2};
Line(nbLines4 + 4) = {nbPts4 + 2, nbPts4 + 3};
-Line(nbLines4 + 5) = {nbPts4 + 3, nbPts4 + 4};
-Line(nbLines4 + 6) = {nbPts4 + 4, nbPts2 + 1};
+Line(nbLines4 + 5) = {nbPts4 + 3, nbPts2 + 1};
-Curve Loop(11) = {nbLines4 + 2: nbLines4 + 6, nbLines3 + 1: nbLines4 - 1, nbLines4 + 1, 26:nbLines2, 24, 19, 12:16, 9, 1:5, 7, 10, 17, 20, 21};
+Curve Loop(11) = {nbLines4 + 2: nbLines4 + 5, nbLines3 + 1: nbLines4 - 1, nbLines4 + 1, 26:nbLines2, 24, 19, 12:13};
Plane Surface(11) = {11};
Physical Surface("BEM_domain_1", 6) = {11};
-Physical Curve("rest_of_outside", 7) = {nbLines4+1: nbLines4+6};
\ No newline at end of file
+Physical Curve("rest_of_outside", 7) = {nbLines4+1: nbLines4+5};
\ No newline at end of file
diff --git a/srcs/MEMS.geo b/srcs/MEMS.geo
index a918f61c3068a787a1b17244848dfd096e19f79e..1caef0d5d166b5ad1b61ec8fcf77a4c51d5462cb 100644
--- a/srcs/MEMS.geo
+++ b/srcs/MEMS.geo
@@ -9,7 +9,7 @@ Lx_diaph = 2*scale;
nx = 20;
ny = 50;
-// monter à 500+ steps en FEM non linéaire car converge très très lentement
+// THIS .geo FILE is mentioned once IN THE REPORT, in section 5.3.
Point(1) = {0, 0, 0, 2};
Point(2) = {Lx_foot, 0, 0, 2};
@@ -90,7 +90,7 @@ SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
SetNumber("Materials/FEM_domain/rho",7800);
-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/tx", 0.);
SetNumber("Boundary Conditions/top_edge/ty", 0.);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
SetNumber("Volumic Forces/FEM_domain/by",0.);
diff --git a/srcs/beamActuation.geo b/srcs/beamActuation.geo
index a37ee6d35de807819cba90a988a91c0491262946..fbec4f3386b33cc611be8ea64f4faf12c9e2671e 100644
--- a/srcs/beamActuation.geo
+++ b/srcs/beamActuation.geo
@@ -1,5 +1,7 @@
scale = 1e-6;
+// THIS .geo FILE IS NOT USED IN THE REPORT
+
bt = 1; // Beam tickness
bh = 1; // Beam height (position)
diff --git a/srcs/beamActuationGeneral.geo b/srcs/beamActuationGeneral.geo
deleted file mode 100644
index 3ad6cf8a80d55447284b1467115cedc983abd0e0..0000000000000000000000000000000000000000
--- a/srcs/beamActuationGeneral.geo
+++ /dev/null
@@ -1,229 +0,0 @@
-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
diff --git a/srcs/clampedMicroBeam.geo b/srcs/clampedMicroBeam.geo
index 5e94f6d3570ef5bd1d8605fdd2833d07b4ee2f59..8b19c2193394b977c8f093c4bd56cf441d51a9cc 100644
--- a/srcs/clampedMicroBeam.geo
+++ b/srcs/clampedMicroBeam.geo
@@ -10,13 +10,13 @@ n = 6; // FEM MESH DENSITY
// additional parameters given to the solver
SetNumber("Boundary Conditions/left_edge/ux", 0.);
SetNumber("Boundary Conditions/left_edge/uy", 0);
-SetNumber("Materials/domain/Young", 150e9);
-SetNumber("Materials/domain/Poisson", 0.27);
-SetNumber("Materials/domain/rho",2300);
+SetNumber("Materials/FEM_domain/Young", 150e9);
+SetNumber("Materials/FEM_domain/Poisson", 0.27);
+SetNumber("Materials/FEM_domain/rho",2300);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
-SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
+SetNumber("Volumic Forces/FEM_domain/by",0.);
-phi_top = 112; // Ã modifier valou
+phi_top = 111.2;
SetNumber("Boundary Conditions/mass/BEM_domain_1/dirichlet", 0);
SetNumber("Boundary Conditions/BEM_FEM_boundary/BEM_domain_1/dirichlet", phi_top);
SetNumber("Boundary Conditions/rest_of_outside/BEM_domain_1/neumann", 0);
diff --git a/srcs/coupling_validation.geo b/srcs/coupling_validation.geo
index e5510d1e2ee35de5eebe1b186d0f0be9839fefe1..78440018e8001ebf4a75ebb44cf8abcf605323eb 100644
--- a/srcs/coupling_validation.geo
+++ b/srcs/coupling_validation.geo
@@ -44,19 +44,19 @@ Physical Surface("FEM_domain", 2) = {1};
Physical Point("fixed_node", 3) = {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/ux", 0.);
//SetNumber("Boundary Conditions/left_edge/uy", 0.);
SetNumber("Boundary Conditions/fixed_node/uy", 0.);
SetNumber("Materials/FEM_domain/Young", 210e3);
SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_domain/rho",7800); //volumic mass of acier
+SetNumber("Materials/FEM_domain/rho",7800);
Physical Curve("BEM_FEM_boundary", 4) = {2};
Physical Surface("BEM_domain_1", 5) = {2};
Physical Curve("right_BEM", 6) = {6};
Physical Curve("homogeneous_field", 7) = {5, 7};
-phi = 30;
+phi = 90;
SetNumber("Boundary Conditions/right_BEM/BEM_domain_1/dirichlet", 0);
SetNumber("Boundary Conditions/BEM_FEM_boundary/BEM_domain_1/dirichlet", phi);
SetNumber("Boundary Conditions/homogeneous_field/BEM_domain_1/neumann", 0);
diff --git a/srcs/hybrid_geo.geo b/srcs/hybrid_geo.geo
index 53ba5b7907d95934a50e77a778638bc84538074a..ef858c7410c7b7db6318730c3a5832abe20d0f95 100644
--- a/srcs/hybrid_geo.geo
+++ b/srcs/hybrid_geo.geo
@@ -4,9 +4,11 @@ Ly_poutre = 2*scale;
h_poutre = 2*scale;
h_tot = 10*scale;
width = 10*scale;
-nx = 20; // prend beaucoup de temps à pd de 200x40
+nx = 20;
ny = 8;
+// THIS .geo FILE IS NOT USED IN THE REPORT
+
Point(1) = {0, h_poutre, 0, 2};
Point(2) = {Lx_poutre, h_poutre, 0, 2};
Point(3) = {Lx_poutre, h_poutre + Ly_poutre, 0, 2};
@@ -59,17 +61,17 @@ Physical Point("fixed_node", 9) = {1};
Physical Surface("BEM_domain_1", 10) = {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/ux", 0.);
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("Materials/domain/rho",7800);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", 0.);
+//SetNumber("Boundary Conditions/right_edge/tx", 21e3);
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("Volumic Forces/FEM_domain/by",0.);
//SetNumber("Boundary Conditions/fixed_node/uy",2.);
// BEM geometry
@@ -81,8 +83,5 @@ SetNumber("Boundary Conditions/mass/BEM_domain_1/dirichlet", 0);
SetNumber("Boundary Conditions/BEM_FEM_boundary/BEM_domain_1/dirichlet", phi_top);
SetNumber("Boundary Conditions/rest_of_outside/BEM_domain_1/neumann", 0);
SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
-//Physical Curve("BEM_boundary", 14) = {1,3,5,6,7,8,9,10,2};
-Physical Curve("Electrode", 15) = {1};
-//SetNumber("Boundary Conditions/bottom_edge/tx", 0);
-//SetNumber("Boundary Conditions/bottom_edge/ty", -25000);
\ No newline at end of file
+Physical Curve("Electrode", 15) = {1};
\ No newline at end of file
diff --git a/srcs/longitudinalCombDevice.geo b/srcs/longitudinalCombDevice.geo
index 098922a6252baea88986052ff3e46fad2b5af584..8398239152c737c0473d347a3969cd5b81471d3b 100644
--- a/srcs/longitudinalCombDevice.geo
+++ b/srcs/longitudinalCombDevice.geo
@@ -1,13 +1,15 @@
scale = 2e-6;
// USE WITH MINIMUM 2 FINS, else use longitudinal_comb.geo
-N_fins = 3; // number of fins on one side of the comb // MAX 12 si on change ps la longueur
+N_fins = 4; // number of fins on one side of the comb // MAX 12 si on change ps la longueur
// WARNING: when using more fins the pull-in voltage decreases
-n = 1; // FEM elements density
+n = 6; // FEM elements density
nBEM = 1; // BEM elements density
+a = 0; // acceleration verticale
+
// mechanical properties and boundary conditions
SetNumber("Boundary Conditions/left/ux", 0.); // encastrement
SetNumber("Boundary Conditions/left/uy", 0);
@@ -17,10 +19,12 @@ SetNumber("Materials/FEM_domain/Young", 150e9);
SetNumber("Materials/FEM_domain/Poisson", 0.27);
SetNumber("Materials/FEM_domain/rho",2300);
SetNumber("Volumic Forces/FEM_domain/bx",0);
-SetNumber("Volumic Forces/FEM_domain/by",00); // acceleration of accelerometer
+SetNumber("Volumic Forces/FEM_domain/by",a); // acceleration of accelerometer
+
+SetNumber("Non_linear_solver", 1);
// BEM properties in bottom domain
-phi_1 = 50;
+phi_1 = 20;
SetNumber("Boundary Conditions/BEM_FEM_boundary_1/BEM_domain_1/dirichlet", 0);
SetNumber("Boundary Conditions/electrode_1/BEM_domain_1/dirichlet", phi_1);
SetNumber("Boundary Conditions/outside_1/BEM_domain_1/neumann", 0);
@@ -42,7 +46,7 @@ h_fin_elec = 2.4*scale; // length of the fins of the electrode, can be longer th
//l_bord = 10*scale;
//l_tot = 34.4*scale;
-l_tot = 40*scale;
+l_tot = 60*scale;
l_fin = 0.8*scale;
l_space = 0.8*scale;
t_electrode = 1*scale; // width of one electrode
diff --git a/srcs/longitudinal_comb.geo b/srcs/longitudinal_comb.geo
deleted file mode 100644
index 6c89eec081168abdcadd2fbd9403a00c051875fd..0000000000000000000000000000000000000000
--- a/srcs/longitudinal_comb.geo
+++ /dev/null
@@ -1,147 +0,0 @@
-scale = 1e-5;
-
-h_tot = 20*scale;
-h_base = 1*scale;
-h_pin = 2.5*scale;
-h_space = 1.2*scale;
-
-l_tot = 20*scale;
-l_pin = 0.8*scale;
-l_space = 0.5*scale;
-
-t_electrode = 1*scale; // width of fixed electrodes
-// FEM domain
-Point(1) = {-l_tot/2, h_base/2, 0, 0.2*scale};
-Point(2) = {-l_pin/2, h_base/2, 0, 0.2*scale};
-Point(3) = {-l_pin/2, h_base/2 + h_pin, 0, 0.2*scale};
-Point(4) = {l_pin/2, h_base/2 + h_pin, 0, 0.2*scale};
-Point(5) = {l_pin/2, h_base/2, 0, 0.2*scale};
-Point(6) = {l_tot/2, h_base/2, 0, 0.2*scale};
-Point(7) = {l_tot/2, -h_base/2, 0, 0.2*scale};
-Point(8) = {l_pin/2, -h_base/2, 0, 0.2*scale};
-Point(9) = {l_pin/2, -(h_base/2 + h_pin), 0, 0.2*scale};
-Point(10) = {-l_pin/2, -(h_base/2 + h_pin), 0, 0.2*scale};
-Point(11) = {-l_pin/2, -h_base/2, 0, 0.2*scale};
-Point(12) = {-l_tot/2, -h_base/2, 0, 0.2*scale};
-
-// FEM domain
-Line(1) = {1, 2};
-Line(2) = {2, 3};
-Line(3) = {3, 4};
-Line(4) = {4, 5};
-Line(5) = {5, 6};
-Line(6) = {6, 7};
-Line(7) = {7, 8};
-Line(8) = {8, 9};
-Line(9) = {9, 10};
-Line(10) = {10, 11};
-Line(11) = {11, 12};
-Line(12) = {12, 1};
-
-Curve Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
-Plane Surface(1) = {1};
-
-Recombine Surface {1}; // quads instead of triangles
-
-Physical Curve("left", 1) = {12};
-Physical Curve("right", 2) = {6};
-
-Physical Surface("FEM_domain", 3) = {1};
-Physical Curve("BEM_FEM_boundary", 4) = {1, 2, 3, 4, 5, 7, 8, 9, 10, 11};
-
-// mechanical properties and boundary conditions
-SetNumber("Boundary Conditions/left/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/left/uy", 0);
-SetNumber("Boundary Conditions/right/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/right/uy", 0);
-SetNumber("Materials/domain/Young", 210e3); // A DETERMINER PRECISEMENT
-SetNumber("Materials/domain/Poisson", 0.3);
-SetNumber("Materials/domain/rho",7800); //volumic mass of acier
-SetNumber("Volumic Forces/FEM_domain/bx",0); // acceleration of accelerometer
-SetNumber("Volumic Forces/FEM_domain/by",-200.);
-
-// BEM domain 1
-Point(13) = {-(l_pin/2 + l_space + t_electrode), -(h_base/2 + h_pin + h_space + t_electrode), 0, 0.2*scale}; // first fixed electrode
-Point(14) = {-(l_pin/2 + l_space + t_electrode), -(h_base/2 + h_space), 0, 0.2*scale};
-Point(15) = {-(l_pin/2 + l_space), -(h_base/2 + h_space), 0, 0.2*scale};
-Point(16) = {-(l_pin/2 + l_space), -(h_base/2 + h_pin + h_space), 0, 0.2*scale};
-Point(17) = {l_pin/2 + l_space, -(h_base/2 + h_pin + h_space), 0, 0.2*scale};
-Point(18) = {l_pin/2 + l_space, -(h_base/2 + h_space), 0, 0.2*scale};
-Point(19) = {l_pin/2 + l_space + t_electrode, -(h_base/2 + h_space), 0, 0.2*scale};
-Point(20) = {l_pin/2 + l_space + t_electrode, -(h_base/2 + h_pin + h_space + t_electrode), 0, 0.2*scale};
-
-Line(13) = {13, 14};
-Line(14) = {14, 15};
-Line(15) = {15, 16};
-Line(16) = {16, 17};
-Line(17) = {17, 18};
-Line(18) = {18, 19};
-Line(19) = {19, 20};
-Line(20) = {20, 13};
-
-Point(29) = {-l_tot/2, -h_tot/2, 0, 0.2*scale};
-Point(30) = {l_tot/2, -h_tot/2, 0, 0.2*scale};
-
-Line(29) = {12, 29};
-Line(30) = {29, 30};
-Line(31) = {30, 7};
-
-Curve Loop(2) = {29, 30, 31, 7, 8, 9, 10, 11};
-Curve Loop(3) = {13, 14, 15, 16, 17, 18, 19, 20};
-Plane Surface(2) = {2, 3};
-
-Physical Surface("BEM_domain_1", 5) = {2};
-Physical Curve("BEM_FEM_boundary_1", 6) = {7, 8, 9, 10, 11};
-Physical Curve("electrode_1", 7) = {13, 14, 15, 16, 17, 18, 19, 20};
-Physical Curve("outside_1", 8) = {29, 30, 31};
-
-phi_1 = 0.1;
-SetNumber("Boundary Conditions/BEM_FEM_boundary_1/BEM_domain_1/dirichlet", 0);
-SetNumber("Boundary Conditions/electrode_1/BEM_domain_1/dirichlet", phi_1);
-SetNumber("Boundary Conditions/outside_1/BEM_domain_1/neumann", 0);
-SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
-
-// BEM domain 2
-Point(21) = {-(l_pin/2 + l_space + t_electrode), h_base/2 + h_space, 0, 0.2*scale}; // second fixed electrode
-Point(22) = {-(l_pin/2 + l_space + t_electrode), h_base/2 + h_pin + h_space + t_electrode, 0, 0.2*scale};
-Point(23) = {l_pin/2 + l_space + t_electrode, h_base/2 + h_pin + h_space + t_electrode, 0, 0.2*scale};
-Point(24) = {l_pin/2 + l_space + t_electrode, h_base/2 + h_space, 0, 0.2*scale};
-Point(25) = {l_pin/2 + l_space, h_base/2 + h_space, 0, 0.2*scale};
-Point(26) = {l_pin/2 + l_space, h_base/2 + h_pin + h_space, 0, 0.2*scale};
-Point(27) = {-(l_pin/2 + l_space), h_base/2 + h_pin + h_space, 0, 0.2*scale};
-Point(28) = {-(l_pin/2 + l_space), h_base/2 + h_space, 0, 0.2*scale};
-
-
-Line(21) = {21, 22};
-Line(22) = {22, 23};
-Line(23) = {23, 24};
-Line(24) = {24, 25};
-Line(25) = {25, 26};
-Line(26) = {26, 27};
-Line(27) = {27, 28};
-Line(28) = {28, 21};
-
-Point(31) = {l_tot/2, h_tot/2, 0, 0.2*scale};
-Point(32) = {-l_tot/2, h_tot/2, 0, 0.2*scale};
-
-Line(32) = {6, 31};
-Line(33) = {31, 32};
-Line(34) = {32, 1};
-
-Curve Loop(4) = {1, 2, 3, 4, 5, 32, 33, 34};
-Curve Loop(5) = {21, 22, 23, 24, 25, 26, 27, 28};
-Plane Surface(3) = {4, 5};
-
-Physical Surface("BEM_domain_2", 9) = {3};
-Physical Curve("BEM_FEM_boundary_2", 10) = {1, 2, 3, 4, 5};
-Physical Curve("electrode_2", 11) = {21, 22, 23, 24, 25, 26, 27, 28};
-Physical Curve("outside_2", 12) = {32, 33, 34};
-
-phi_2 = 25;
-SetNumber("Boundary Conditions/BEM_FEM_boundary_2/BEM_domain_2/dirichlet", 0);
-SetNumber("Boundary Conditions/electrode_2/BEM_domain_2/dirichlet", phi_2);
-SetNumber("Boundary Conditions/outside_2/BEM_domain_2/neumann", 0);
-SetNumber("Materials/BEM_domain_2/Epsilon", 8.8541878128e-12); // dielectric permittivity
-
-
-Mesh.Algorithm = 8;
\ No newline at end of file
diff --git a/srcs/main.cpp b/srcs/main.cpp
index 1bacc6c05934e5e4f838becfe2b4f1a0f656bf6f..817334806b0ad1d63d3f7b4069e6c25a9b3f25f5 100644
--- a/srcs/main.cpp
+++ b/srcs/main.cpp
@@ -33,6 +33,9 @@ int main(int argc, char **argv)
if(nonLinearSolver)
{
+ std::cout << "-----------------------------------------\n";
+ std::cout << "NON-LINEAR TWO-WAY COUPLED ITERATIVE SOLVER\n";
+
bool postProcessing = 0; // pass it as an argument to the solver, only compute post processing at last iteration
std::map<int, double> electrostaticPressure;
@@ -44,10 +47,9 @@ int main(int argc, char **argv)
int iteration;
- int viewTagU = gmsh::view::add("u"); // Ã modifier plus tard
+ int viewTagU = gmsh::view::add("u"); // passed as an argument to the FEM solver
int viewTagF = gmsh::view::add("f");
-
solverBEM(electrostaticPressure, nbViews, boundaryDisplacementMap, postProcessing, 1, false);
solverFEMnonLinear(electrostaticPressure, boundaryDisplacementMap, nbViews, postProcessing, 1, viewTagU, viewTagF, false);
@@ -57,9 +59,10 @@ int main(int argc, char **argv)
std::vector<double> relativeDifference(boundaryDisplacementMap.size());
double relativeDisplacement;
- const double criticalNorm = 0.0000001;
- const int maxNbIteration = 50; // Ã modifier
+ const double criticalNorm = 1e-7;
+ const int maxNbIteration = 50;
+ // see report section 5.2 "stopping criterion" for the theoretical explanation
int i = 0;
for(auto p : boundaryDisplacementMap)
{
@@ -89,7 +92,7 @@ int main(int argc, char **argv)
i++;
}
relativeDisplacement = sqrt(relativeDisplacement/i);
- //std::cout << iteration << ", " << relativeDisplacement << ";\n";
+
std::cout << "Iteration: " << iteration << "\t-> Averaged relative displacement = " << relativeDisplacement << "\n";
if(relativeDisplacement < criticalNorm)
break;
@@ -114,6 +117,9 @@ int main(int argc, char **argv)
else
{
+ std::cout << "-----------------------------------------\n";
+ std::cout << "LINEAR ONE-WAY COUPLED SOLVER\n";
+
std::map<int, double> electrostaticPressure;
std::map<int,std::pair<double,double>> boundaryDisplacementMap;
diff --git a/srcs/squeeze.geo b/srcs/squeeze.geo
index ad1a4a9bd13f24d997c7c2e461eb74453df44709..fc23efd5c4281e22c98503bab48d0d68d5ad7be6 100644
--- a/srcs/squeeze.geo
+++ b/srcs/squeeze.geo
@@ -1,25 +1,36 @@
scale = 4e-6;
-h_tot = 5*scale;
-w_tot = 10*scale;
-h_foot = 4*scale;
-h_diaph = 4.4*scale;
-h_end_diaph = 4.6*scale;
-Lx_diaph = 1*scale;
-
-nx = 10; // prend beaucoup de temps à pd de 200x40
-ny = 4;
-
-Point(1) = {0, 0, 0, 2};
-Point(2) = {w_tot, 0, 0, 2};
-Point(3) = {w_tot, h_foot, 0, 2};
+h_tot0 = 5;
+w_tot0 = 10;
+h_foot0 = 4;
+h_diaph0 = 4.4;
+h_end_diaph0 = 4.6;
+Lx_diaph0 = 1;
+
+unit_l = 0.2; // serves as a reference
+
+h_tot = h_tot0*scale;
+w_tot = w_tot0*scale;
+h_foot = h_foot0*scale;
+h_diaph = h_diaph0*scale;
+h_end_diaph = h_end_diaph0*scale;
+Lx_diaph = Lx_diaph0*scale;
+
+n = 8; // FEM elements density
+nBEM = 0.25; // BEM elements density
+
+// THIS .geo FILE IS NOT USED IN THE REPORT
+
+Point(1) = {0, 0, 0, 0.8*nBEM*scale};
+Point(2) = {w_tot, 0, 0, 0.8*nBEM*scale};
+Point(3) = {w_tot, h_foot, 0, 0.4*nBEM*scale};
Point(4) = {w_tot, h_tot, 0, 2};
Point(5) = {w_tot-Lx_diaph, h_end_diaph, 0, 2};
Point(6) = {Lx_diaph, h_end_diaph, 0, 2};
Point(7) = {0, h_tot, 0, 2};
-Point(8) = {0, h_foot, 0, 2};
-Point(9) = {Lx_diaph, h_diaph, 0, 2};
-Point(10) = {w_tot-Lx_diaph, h_diaph, 0, 2};
+Point(8) = {0, h_foot, 0, 0.4*nBEM*scale};
+Point(9) = {Lx_diaph, h_diaph, 0, 0.1*nBEM*scale};
+Point(10) = {w_tot-Lx_diaph, h_diaph, 0, 0.1*nBEM*scale};
Line(1) = {9, 8};
Line(2) = {10, 9};
@@ -34,22 +45,31 @@ Line(9) = {8, 1};
Line(10) = {1, 2};
Line(11) = {2 ,3};
-Curve Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8};
+Line(13) = {9, 6};
+Line(14) = {5, 10};
+
+Curve Loop(1) = {-13, 1, 8, 7};
Plane Surface(1) = {1};
-Curve Loop(2) = {9, 10, 11, 3, 2, 1};
-Plane Surface(2) = {2};
+Curve Loop(2) = {2, 13, 6, 14};
+Plane Surface(2) = {2};
+
+Curve Loop(3) = {-14, 5, 4, 3};
+Plane Surface(3) = {3};
+
+Curve Loop(4) = {9, 10, 11, 3, 2, 1};
+Plane Surface(4) = {4};
-Transfinite Curve {10} = 10*nx+1 Using Progression 1;
-Transfinite Curve {1, 2, 3} = 3*nx+1 Using Progression 1;
-Transfinite Curve {7, 6, 5} = 3*nx+1 Using Progression 1;
-Transfinite Curve {8, 4} = 2*ny+1 Using Progression 1;
-Transfinite Curve {7, 1, 5, 3} = 1*ny+1 Using Progression 1;
-Transfinite Curve {9, 11} = 4*ny+1 Using Progression 1;
+Transfinite Curve {13, 14} = (h_end_diaph0 - h_diaph0)/unit_l*n+1 Using Progression 1;
+Transfinite Curve {2, 6} = (w_tot0 - 2*Lx_diaph0)/unit_l*n+1 Using Progression 1;
+Transfinite Curve {4, 8} = (h_tot0-h_foot0)/unit_l*n+1 Using Progression 1;
+Transfinite Curve {1, 3, 5, 7} = (h_tot0-h_foot0)/unit_l*n+1 Using Progression 1;
-//Transfinite Surface {1};
+Transfinite Surface {2};
Recombine Surface {1};
+Recombine Surface {2};
+Recombine Surface {3};
Mesh.ElementOrder = 1;
@@ -62,8 +82,8 @@ Physical Curve ("right_up_foot", 10) = {5};
Physical Curve ("up_diaph", 11) = {6};
Physical Curve ("low_diaph", 12) = {2};
-Physical Surface("FEM_domain", 13) = {1};
-Physical Surface("BEM_domain_1", 14) = {2};
+Physical Surface("FEM_domain", 13) = {1,2,3};
+Physical Surface("BEM_domain_1", 14) = {4};
SetNumber("Boundary Conditions/left_edge_foot/ux", 0.);
@@ -79,9 +99,9 @@ SetNumber("Boundary Conditions/right_low_foot/uy", 0);
SetNumber("Boundary Conditions/right_up_foot/ux", 0.);
SetNumber("Boundary Conditions/right_up_foot/uy", 0);
-SetNumber("Materials/FEM_domain/Young", 210e3);
-SetNumber("Materials/FEM_domain/Poisson", 0.3);
-SetNumber("Materials/FEM_domain/rho",7800);
+SetNumber("Materials/FEM_domain/Young", 150e9);
+SetNumber("Materials/FEM_domain/Poisson", 0.27);
+SetNumber("Materials/FEM_domain/rho",2300);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
SetNumber("Volumic Forces/FEM_domain/by",0.);
diff --git a/srcs/transverse_comb.geo b/srcs/transverse_comb.geo
deleted file mode 100644
index 21122d65ca96ae9a77466cb6eb8c97cf2f302e65..0000000000000000000000000000000000000000
--- a/srcs/transverse_comb.geo
+++ /dev/null
@@ -1,157 +0,0 @@
-scale = 1e-5;
-
-h_tot = 20*scale;
-h_base = 3*scale;
-h_pin = 2.5*scale;
-h_space = 1*scale;
-
-l_tot = 20*scale;
-l_base = 1*scale;
-l_pin = 0.8*scale;
-l_space = 1.2*scale;
-
-t_electrode = 1*scale; // width of fixed electrodes
-// FEM domain
-Point(1) = {-l_tot/2, -h_tot/2, 0, 0.2*scale};
-Point(2) = {-l_tot/2, h_tot/2, 0, 0.2*scale};
-Point(3) = {-(l_tot/2 - l_base), h_tot/2, 0, 0.2*scale};
-Point(4) = {-(l_tot/2 - l_base), h_base/2, 0, 0.2*scale};
-Point(5) = {-l_pin/2, h_base/2, 0, 0.2*scale};
-Point(6) = {-l_pin/2, h_base/2 + h_pin, 0, 0.2*scale};
-Point(7) = {l_pin/2, h_base/2 + h_pin, 0, 0.2*scale};
-Point(8) = {l_pin/2, h_base/2, 0, 0.2*scale};
-Point(9) = {(l_tot/2 - l_base), h_base/2, 0, 0.2*scale};
-Point(10) = {(l_tot/2 - l_base), h_tot/2, 0, 0.2*scale};
-Point(11) = {l_tot/2, h_tot/2, 0, 0.2*scale};
-Point(12) = {l_tot/2, -h_tot/2, 0, 0.2*scale};
-Point(13) = {(l_tot/2 - l_base), -h_tot/2, 0, 0.2*scale};
-Point(14) = {(l_tot/2 - l_base), -h_base/2, 0, 0.2*scale};
-Point(15) = {l_pin/2, -h_base/2, 0, 0.2*scale};
-Point(16) = {l_pin/2, -(h_base/2 + h_pin), 0, 0.2*scale};
-Point(17) = {-l_pin/2, -(h_base/2 + h_pin), 0, 0.2*scale};
-Point(18) = {-l_pin/2, -h_base/2, 0, 0.2*scale};
-Point(19) = {-(l_tot/2 - l_base), -h_base/2, 0, 0.2*scale};
-Point(20) = {-(l_tot/2 - l_base), -h_tot/2, 0, 0.2*scale};
-
-// FEM domain
-Line(1) = {1, 2};
-Line(2) = {2, 3};
-Line(3) = {3, 4};
-Line(4) = {4, 5};
-Line(5) = {5, 6};
-Line(6) = {6, 7};
-Line(7) = {7, 8};
-Line(8) = {8, 9};
-Line(9) = {9, 10};
-Line(10) = {10, 11};
-Line(11) = {11, 12};
-Line(12) = {12, 13};
-Line(13) = {13, 14};
-Line(14) = {14, 15};
-Line(15) = {15, 16};
-Line(16) = {16, 17};
-Line(17) = {17, 18};
-Line(18) = {18, 19};
-Line(19) = {19, 20};
-Line(20) = {20, 1};
-
-Curve Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
-Plane Surface(1) = {1};
-
-Recombine Surface {1}; // quads instead of triangles
-
-Physical Curve("bottom_left", 1) = {20};
-Physical Curve("top_left", 2) = {2};
-Physical Curve("top_right", 3) = {10};
-Physical Curve("bottom_right", 4) = {12};
-
-Physical Surface("FEM_domain", 5) = {1};
-Physical Curve("BEM_FEM_boundary", 6) = {3, 4, 5, 6, 7, 8, 9, 13, 14, 15, 16, 17, 18, 19};
-
-// mechanical properties and boundary conditions
-SetNumber("Boundary Conditions/bottom_left/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/bottom_left/uy", 0);
-SetNumber("Boundary Conditions/top_left/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/top_left/uy", 0);
-SetNumber("Boundary Conditions/top_right/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/top_right/uy", 0);
-SetNumber("Boundary Conditions/bottom_right/ux", 0.); // encastrement
-SetNumber("Boundary Conditions/bottom_right/uy", 0);
-SetNumber("Materials/domain/Young", 210e3); // A DETERMINER PRECISEMENT
-SetNumber("Materials/domain/Poisson", 0.3);
-SetNumber("Materials/domain/rho",7800); //volumic mass of acier
-SetNumber("Volumic Forces/FEM_domain/bx",50); // acceleration of accelerometer
-SetNumber("Volumic Forces/FEM_domain/by",0.);
-
-// BEM domain 1
-Point(21) = {-(l_pin/2 + l_space + t_electrode), -(h_base/2 + h_pin + h_space + t_electrode), 0, 0.2*scale}; // first fixed electrode
-Point(22) = {-(l_pin/2 + l_space + t_electrode), -(h_base/2 + h_space), 0, 0.2*scale};
-Point(23) = {-(l_pin/2 + l_space), -(h_base/2 + h_space), 0, 0.2*scale};
-Point(24) = {-(l_pin/2 + l_space), -(h_base/2 + h_pin + h_space), 0, 0.2*scale};
-Point(25) = {l_pin/2 + l_space, -(h_base/2 + h_pin + h_space), 0, 0.2*scale};
-Point(26) = {l_pin/2 + l_space, -(h_base/2 + h_space), 0, 0.2*scale};
-Point(27) = {l_pin/2 + l_space + t_electrode, -(h_base/2 + h_space), 0, 0.2*scale};
-Point(28) = {l_pin/2 + l_space + t_electrode, -(h_base/2 + h_pin + h_space + t_electrode), 0, 0.2*scale};
-
-Line(21) = {21, 22};
-Line(22) = {22, 23};
-Line(23) = {23, 24};
-Line(24) = {24, 25};
-Line(25) = {25, 26};
-Line(26) = {26, 27};
-Line(27) = {27, 28};
-Line(28) = {28, 21};
-
-Line(37) = {20, 13};
-
-Curve Loop(2) = {37, 13, 14, 15, 16, 17, 18, 19};
-Curve Loop(3) = {21, 22, 23, 24, 25, 26, 27, 28};
-Plane Surface(2) = {2, 3};
-
-Physical Surface("BEM_domain_1", 7) = {2};
-Physical Curve("BEM_FEM_boundary_1", 8) = {13, 14, 15, 16, 17, 18, 19};
-Physical Curve("electrode_1", 9) = {21, 22, 23, 24, 25, 26, 27, 28};
-Physical Curve("outside_1", 10) = {37};
-
-phi_1 = 30;
-SetNumber("Boundary Conditions/BEM_FEM_boundary_1/BEM_domain_1/dirichlet", 0);
-SetNumber("Boundary Conditions/electrode_1/BEM_domain_1/dirichlet", phi_1);
-SetNumber("Boundary Conditions/outside_1/BEM_domain_1/neumann", 0);
-SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
-
-// BEM domain 2
-Point(29) = {-(l_pin/2 + l_space + t_electrode), h_base/2 + h_space, 0, 0.2*scale}; // second fixed electrode
-Point(30) = {-(l_pin/2 + l_space + t_electrode), h_base/2 + h_pin + h_space + t_electrode, 0, 0.2*scale};
-Point(31) = {l_pin/2 + l_space + t_electrode, h_base/2 + h_pin + h_space + t_electrode, 0, 0.2*scale};
-Point(32) = {l_pin/2 + l_space + t_electrode, h_base/2 + h_space, 0, 0.2*scale};
-Point(33) = {l_pin/2 + l_space, h_base/2 + h_space, 0, 0.2*scale};
-Point(34) = {l_pin/2 + l_space, h_base/2 + h_pin + h_space, 0, 0.2*scale};
-Point(35) = {-(l_pin/2 + l_space), h_base/2 + h_pin + h_space, 0, 0.2*scale};
-Point(36) = {-(l_pin/2 + l_space), h_base/2 + h_space, 0, 0.2*scale};
-
-
-Line(29) = {29, 30};
-Line(30) = {30, 31};
-Line(31) = {31, 32};
-Line(32) = {32, 33};
-Line(33) = {33, 34};
-Line(34) = {34, 35};
-Line(35) = {35, 36};
-Line(36) = {36, 29};
-
-Line(38) = {10, 3};
-
-Curve Loop(4) = {38, 3, 4, 5, 6, 7, 8, 9};
-Curve Loop(5) = {29, 30, 31, 32, 33, 34, 35, 36};
-Plane Surface(3) = {4, 5};
-
-Physical Surface("BEM_domain_2", 11) = {3};
-Physical Curve("BEM_FEM_boundary_2", 12) = {3, 4, 5, 6, 7, 8, 9};
-Physical Curve("electrode_2", 13) = {29, 30, 31, 32, 33, 34, 35, 36};
-Physical Curve("outside_2", 14) = {38};
-
-phi_2 = 30;
-SetNumber("Boundary Conditions/BEM_FEM_boundary_2/BEM_domain_2/dirichlet", 0);
-SetNumber("Boundary Conditions/electrode_2/BEM_domain_2/dirichlet", phi_2);
-SetNumber("Boundary Conditions/outside_2/BEM_domain_2/neumann", 0);
-SetNumber("Materials/BEM_domain_2/Epsilon", 8.8541878128e-12); // dielectric permittivity
\ No newline at end of file
diff --git a/srcs/two_BEM_test.geo b/srcs/two_BEM_test.geo
index 89926c27ade94d3510c858b0f78b71bc8ac2cc14..449a3efd30ec9a30eccfceae47ea7e4505dbf785 100644
--- a/srcs/two_BEM_test.geo
+++ b/srcs/two_BEM_test.geo
@@ -7,6 +7,8 @@ h_FEM = 4.5*scale;
nx = 40;
ny = 4;
+// THIS .geo FILE IS NOT USED IN THE REPORT
+
Point(1) = {0, 0, 0, 2};
Point(2) = {w_tot, 0, 0, 2};
Point(3) = {w_tot, h_FEM, 0, 2};