From 60f360c8c34683007b4aa7a70af135b1f2d89d52 Mon Sep 17 00:00:00 2001
From: Louis Denis <louis.denis@student.uliege.be>
Date: Fri, 20 May 2022 23:28:03 +0200
Subject: [PATCH] final update .geo files

---
 srcs/FEM/complex_validation.geo               |  6 ++---
 srcs/FEM/hybrid_complex.geo                   | 18 ++++++-------
 srcs/FEM/hybrid_geo.geo                       | 14 ++++++-----
 srcs/FEM/large_rotation_validation.geo        |  6 ++---
 ...rge_rotation_validation_single_surface.geo |  2 ++
 srcs/FEM/self_weight.geo                      |  4 +--
 srcs/FEM/simple_tension.geo                   | 14 +++++------
 srcs/FEM/solverFEM.cpp                        |  5 +++-
 srcs/FEM/uniform_charge.geo                   | 25 +++++++++----------
 srcs/FoldedFlexureBeam.geo                    |  2 +-
 srcs/main.cpp                                 | 10 ++++----
 11 files changed, 53 insertions(+), 53 deletions(-)

diff --git a/srcs/FEM/complex_validation.geo b/srcs/FEM/complex_validation.geo
index c6497d7..9d0cd28 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 97260c6..43fec41 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 4c6fb60..ef30163 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 113db4e..2386bc5 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 d15291b..6db9cbc 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 c04a6a8..601c9eb 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 f1452e5..81de6b7 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 7000141..2449aa9 100644
--- a/srcs/FEM/solverFEM.cpp
+++ b/srcs/FEM/solverFEM.cpp
@@ -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 930e441..9baa462 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 551c004..87a395b 100644
--- a/srcs/FoldedFlexureBeam.geo
+++ b/srcs/FoldedFlexureBeam.geo
@@ -1,6 +1,6 @@
 scale = 1e-6;
 
-nFEM = 4; // FEM element density
+nFEM = 2; // FEM element density
 
 nBEM = 5/nFEM; // BEM element density (mettre ça dans quatrième coordonnée)
 
diff --git a/srcs/main.cpp b/srcs/main.cpp
index 1bacc6c..9441044 100644
--- a/srcs/main.cpp
+++ b/srcs/main.cpp
@@ -44,10 +44,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 +56,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 +89,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;
-- 
GitLab