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 7000141bdb6d86985fb81d1639fedbe0c5c8674d..2449aa9d805320566ed2f217d85e529db82683a8 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 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 551c0048db565e2db38ea9148df45ac1dfe5254c..87a395b08a9a518b072e5afa992ce5fa1ff85957 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 1bacc6c05934e5e4f838becfe2b4f1a0f656bf6f..9441044f46fb72de62c7539748e282f8a59afcb3 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;