diff --git a/srcs/FEM/complex_validation.geo b/srcs/FEM/complex_validation.geo
index 47e62fe7faf89f96039b9b0497e4db00896d6c36..2c7851a8f2a460fc47db4ba3bc424a65a27e20be 100644
--- a/srcs/FEM/complex_validation.geo
+++ b/srcs/FEM/complex_validation.geo
@@ -2,10 +2,10 @@
Lx = 5;
Ly = 2;
-nx = 50;
-ny = 20;
+nx = 250;
+ny = 100;
-Point(1) = {0, 0, 0, 0.1};
+Point(1) = {0, 0, 0, 0.1}; //4th number is mesh density
Point(2) = {Lx, 0, 0, 0.2};
Point(3) = {Lx, Ly, 0, 0.2};
Point(4) = {0, Ly, 0, 0.1};
@@ -52,6 +52,8 @@ SetNumber("Boundary Conditions/right_edge/tx", 21e3); //for simple tension condi
SetNumber("Boundary Conditions/right_edge/ty", 0.);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
-SetNumber("Boundary Conditions/fixed_node/uy",0.);
+SetNumber("Boundary Conditions/fixed_node/uy",2);
Physical Curve("BEM_FEM_boundary", 10) = {1,2,3,4,5,6};
+
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/FEM/functionsFEM.cpp b/srcs/FEM/functionsFEM.cpp
index 27ec51c66f0290e3b5486ea01358e326146865f2..e6cb0bde065691b4898abcfb4dcad297221b1d84 100644
--- a/srcs/FEM/functionsFEM.cpp
+++ b/srcs/FEM/functionsFEM.cpp
@@ -20,6 +20,30 @@
#endif
#define PI M_PI
+bool getNonLinearParameter()
+{
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "Non_linear_solver");
+ if(keys.size())
+ {
+ std::vector<double> value;
+ gmsh::onelab::getNumber(keys[0], value);
+ if(value[0])
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+ else
+ {
+ std::cout << "you did not choose the type of solver (linear or not) in the .geo file, 'SetNumber('Non_linear_solver',x);'. Automatically set to non linear type \n";
+ return true;
+ }
+}
+
//Compute the matrix H
Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu)
{
@@ -572,8 +596,6 @@ void computeStressStrainElNodal(int & dim, std::vector<int> & tags,
elementData* element;
double theta;
- //double dN_dx_local;
- //double dN_dy_local;
Eigen::Matrix<double, 2, 2> strain_tensor;
Eigen::Matrix<double, 2, 2> stress_tensor;
diff --git a/srcs/FEM/functionsFEM.hpp b/srcs/FEM/functionsFEM.hpp
index 7c162cfa64d4c089fb811350dd3e09de44de560b..4572228eb2bb786d47ea61f12e49f5f2dd2ed86f 100644
--- a/srcs/FEM/functionsFEM.hpp
+++ b/srcs/FEM/functionsFEM.hpp
@@ -41,6 +41,7 @@ struct elementData{
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> B;
};
+bool getNonLinearParameter();
Eigen::Matrix<double, 3, 3> computeHmatrix(const double E, const double nu);
void assembleFlocal(const int &i, const int &el, const std::vector<std::vector<std::size_t>> & nodeTags, const int &numNodes, std::map<int, int> & nodeIndexMap, std::vector<double> & full_f_vector, std::vector<double> & f_vector);
void computeBmatrix(const int &j, const int &numNodes, const std::vector<double> & basisFunctionsGradient, const std::vector<double> & determinants, const Eigen::Matrix<double, 2, 2> & jacobinvtrans, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & B_matrix);
diff --git a/srcs/FEM/large_rotation_validation.geo b/srcs/FEM/large_rotation_validation.geo
index d10805c5adbfd18f62b5b180e270b426f6e8bfef..f2bcf8057a4af3a6fe7e181587a90ccd4f1aaebc 100644
--- a/srcs/FEM/large_rotation_validation.geo
+++ b/srcs/FEM/large_rotation_validation.geo
@@ -2,19 +2,23 @@
h = 1;
H = 10;
-n = 20;
+n = 48;
Point(1) = {0, 0, 0, 0.5};
Point(2) = {0.9*H, 0, 0, 0.5};
Point(3) = {0.9*H, h, 0, 0.5};
Point(4) = {0, h, 0, 0.5};
-Point(5) = {H, 0, 0, 0.5};
-Point(6) = {H, h, 0, 0.5};
+Point(5) = {0.95*H, 0, 0, 0.5};
+Point(6) = {0.95*H, h, 0, 0.5};
-Point(7) = {H, H, 0, 0.5};
+Point(7) = {0.95*H, H, 0, 0.5};
Point(8) = {0.9*H, H, 0, 0.5};
+Point(9) = {H, 0, 0, 0.5};
+Point(10) = {H, h, 0, 0.5};
+Point(11) = {H, H, 0, 0.5};
+
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
@@ -34,21 +38,39 @@ Line(10) = {8, 3};
Curve Loop(3) = {-7, 8, 9, 10};
Plane Surface(3) = {3};
-Transfinite Curve {1, 3, 8, 10} = 9*n+1 Using Progression 1;
-Transfinite Curve {2, 4, 5, 6, 7, 9} = n+1 Using Progression 1;
+Line(11) = {5, 9};
+Line(12) = {9, 10};
+Line(13) = {10, 6};
+Curve Loop(4) = {11, 12, 13, -6};
+Plane Surface(4) = {4};
+
+Line(14) = {10, 11};
+Line(15) = {11, 7};
+Curve Loop(5) = {14, 15, -8, -13};
+Plane Surface(5) = {5};
+
+Transfinite Curve {1, 3, 8, 10, 14} = 18*n+1 Using Progression 1;
+Transfinite Curve {2, 4, 6, 12} = 2*n+1 Using Progression 1;
+Transfinite Curve {5, 7, 9, 11, 13, 15} = n+1 Using Progression 1;
Transfinite Surface {1};
Transfinite Surface {2};
Transfinite Surface {3};
+Transfinite Surface {4};
+Transfinite Surface {5};
Recombine Surface {1};
Recombine Surface {2};
Recombine Surface {3};
+Recombine Surface {4};
+Recombine Surface {5};
Physical Curve("left_edge", 1) = {4};
-Physical Surface("FEM_domain", 2) = {1, 2, 3}; // the trick is to include both plane surfaces in one single domain
-Physical Curve("top_edge", 3) = {9};
+Physical Surface("FEM_domain", 2) = {1, 2, 3, 4, 5}; // the trick is to include both plane surfaces in one single domain
+Physical Curve("top_edge", 3) = {9, 15};
-F = 500;
+Physical Point("point_A", 5) = {7};
+
+F = 40000;
// additional parameters given to the solver
SetNumber("Boundary Conditions/left_edge/ux", 0.); // ALWAYS NEED TO IMPOSE BOTH ux AND uy ON A GIVEN EDGE !! (pas très réaliste, faut y réfléchir)
@@ -61,4 +83,6 @@ 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
diff --git a/srcs/FEM/mainFEM.cpp b/srcs/FEM/mainFEM.cpp
index 00beeaab76d93e3887b078e0a836bc941ec46e7e..02b644c7b50b97264c33f519de48f169856dcb0d 100644
--- a/srcs/FEM/mainFEM.cpp
+++ b/srcs/FEM/mainFEM.cpp
@@ -6,9 +6,6 @@
int main(int argc, char **argv)
{
-
- bool nonLinearSolver = true; // use non-linear solver or not
-
if (argc < 2)
{
std::cout << "Usage: " << argv[0] << " <geo_file>\n";
@@ -28,6 +25,9 @@ int main(int argc, char **argv)
Eigen::initParallel();
+ // determine if non linear solver must be set
+ bool nonLinearSolver = getNonLinearParameter();
+
std::map<int, double> electrostaticPressure;
int nbViews = 0;
diff --git a/srcs/FEM/my_geo.geo b/srcs/FEM/my_geo.geo
index 5c0005f48a4cc907279bccaf78d0ca08318e49da..fde90809eae403a3dd4574e691689a4ee1ed5da6 100644
--- a/srcs/FEM/my_geo.geo
+++ b/srcs/FEM/my_geo.geo
@@ -1,7 +1,7 @@
Lx = 2;
Ly = 1;
-nx = 2; // prend beaucoup de temps à pd de 200x40
-ny = 1;
+nx = 16; // prend beaucoup de temps à pd de 200x40
+ny = 8;
Point(1) = {0, 0, 0, 0.1};
Point(2) = {Lx, 0, 0, 0.1};
@@ -21,7 +21,7 @@ Recombine Surface {1}; // quads instead of triangles
Mesh.ElementOrder = 1;
-Physical Curve("left_edge", 5) = {4,5};
+Physical Curve("left_edge", 5) = {4};
Physical Surface("FEM_domain", 6) = {1};
Physical Curve("top_edge", 7) = {3};
Physical Curve("right_edge", 8) = {2};
@@ -29,17 +29,17 @@ Physical Point("fixed_node", 9) = {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/uy", 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("Boundary Conditions/right_edge/tx", 40); // for simple tension conditions
+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("Boundary Conditions/fixed_node/ux",0.);
-//SetNumber("Boundary Conditions/fixed_node/uy",0.);
+SetNumber("Boundary Conditions/fixed_node/uy",2.);
Physical Curve("BEM_FEM_boundary", 10) = {1,2,3};
\ No newline at end of file
diff --git a/srcs/FEM/self_weight.geo b/srcs/FEM/self_weight.geo
index d50925376b3a42a31e2508c7fa480d9ccdfba644..c04a6a8a44f4dd0dc62137b00f36ea44a94cd666 100644
--- a/srcs/FEM/self_weight.geo
+++ b/srcs/FEM/self_weight.geo
@@ -37,4 +37,6 @@ SetNumber("Boundary Conditions/right_edge/ty", 0.);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
SetNumber("Volumic Forces/FEM_domain/by",-9.81); //set to -9.81 for gravity
-Physical Curve("BEM_FEM_boundary", 9) = {1,2,3}; // useless but necessary to make it work
\ No newline at end of file
+Physical Curve("BEM_FEM_boundary", 9) = {1,2,3}; // useless but necessary to make it work
+
+SetNumber("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/FEM/simple_tension.geo b/srcs/FEM/simple_tension.geo
index 18783ae0b69b461a559836c84cbd892676cec6a4..f1452e518608936d6c4b17101e93a32c0f296abd 100644
--- a/srcs/FEM/simple_tension.geo
+++ b/srcs/FEM/simple_tension.geo
@@ -1,7 +1,7 @@
Lx = 5;
Ly = 2;
-nx = 5;
-ny = 2;
+nx = 50;
+ny = 20;
Point(1) = {0, 0, 0, 1.0};
Point(2) = {Lx, 0, 0, 1.0};
@@ -39,6 +39,8 @@ SetNumber("Boundary Conditions/right_edge/tx", 21e3); // for simple tension cond
SetNumber("Boundary Conditions/right_edge/ty", 0.);
SetNumber("Volumic Forces/FEM_domain/bx",0.);
SetNumber("Volumic Forces/FEM_domain/by",0.); //set to -9.81 for gravity
-SetNumber("Boundary Conditions/fixed_node/uy",0.);
+SetNumber("Boundary Conditions/fixed_node/uy",2.);
Physical Curve("BEM_FEM_boundary", 10) = {1,2,3}; // useless in this case but necessary to make it work
+
+SetNumber("Non_linear_solver",1);
\ No newline at end of file
diff --git a/srcs/FEM/solverFEM.cpp b/srcs/FEM/solverFEM.cpp
index 770715a40ae1e3bb34fd00708b8ade10787d0626..05eb54d07ea8f3b8876a329b8d072da8e46c4d43 100644
--- a/srcs/FEM/solverFEM.cpp
+++ b/srcs/FEM/solverFEM.cpp
@@ -27,6 +27,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
- important remark: in pg we focus only on the free DOFs (we do not focus on the DOFs fixed by dirichlet BCs)*/
bool display_time = false;
+ bool validation = false;
double start = omp_get_wtime();
//int max_threads = omp_get_max_threads();
@@ -217,13 +218,14 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
Eigen::Matrix<double, Eigen::Dynamic, 1> residual(nb_free_DOFs, 1);
- double residualTolerance = 1e-11; // empirical
+ double residualTolerance = 1e-12; // empirical
// implementation of second stopping criterion
std::vector<double> previousTotalNodalForces(2*FEM_nodeTags.size());
std::vector<double> currentTotalNodalForces(2*FEM_nodeTags.size());
std::vector<double> relativeDifference(2*FEM_nodeTags.size());
double relativeForces = 1e17;
+ double previousRelativeForces = 1e17;
double currentMaxNodalForce;
for(size_t i = 0; i < 2*FEM_nodeTags.size(); i++)
{
@@ -237,8 +239,12 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
/*---------ITERATIVE PART---------------*/
- int max_number_of_steps = 50; // purely arbitrary
- for(int step = 0; step < max_number_of_steps; step++)
+ int max_number_of_steps = 200; // purely arbitrary
+ int step = 0;
+
+ double num_diverged_steps = 0; // stores the number of times the residual has increased between two steps, after 5 times,
+ // we take the assumption it has converged towards an other value
+ for(step = 0; step < max_number_of_steps; step++)
{
/*------------STEP 2 of the algorithm----------------*/
/*
@@ -266,6 +272,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
if(step > 1)
{
+ previousRelativeForces = relativeForces;
relativeForces = 0;
currentMaxNodalForce = 0;
for(size_t i = 0; i < 2*FEM_nodeTags.size(); i++)
@@ -285,6 +292,16 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
else
relativeForces = 1e17;
std::cout << "step: " << step << ", relativeForces: " << relativeForces << "\n";
+
+ if(previousRelativeForces < relativeForces)
+ {
+ num_diverged_steps++;
+ if(num_diverged_steps == 5)
+ {
+ std::cout << "-----------------------------------------\n";
+ std::cout << "Pay attention, the differential residual has not reached the tolerance of " << residualTolerance << ", it has converged towards: " << relativeForces << ", after " << step <<" steps\n";
+ }
+ }
}
@@ -303,7 +320,7 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
displayTime(end_retrieveTotalForces, end_assembleglobalf, step_name, display_time);
/*------------STEP 3: EVALUATING THE RESIDUAL, have we converged ??----------*/
- if (relativeForces < residualTolerance && step > 0)
+ if ((relativeForces < residualTolerance && step > 0) || num_diverged_steps == 5)
break;
/*------------STEP 4: ASSEMBLING THE GLOBAL TANGENT STIFFNESS MATRIX----------*/
@@ -342,6 +359,12 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
displayTime(end_solve, end_updateNodalDisp, step_name, display_time);
}
+ if(step == max_number_of_steps)
+ {
+ std::cout << "-----------------------------------------\n";
+ std::cout << "Pay attention, the differential residual has not reached the tolerance of " << residualTolerance << ", final value: " << relativeForces << "\n";
+ }
+
double start_postPro = omp_get_wtime();
//updating the kinematics according to the results of last iteration
@@ -522,6 +545,22 @@ void solverFEMnonLinear(std::map<int, double> &electrostaticPressure, std::map<i
gmsh::model::mesh::setNode(FEM_nodeTags[i], total_disp, {});
}
}
+
+ if(validation)
+ {
+ // point A
+ std::string groupname = "point_A"; // we consider the domain as containing the FEM 2d model
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+ std::vector<std::size_t> A_nodeTag;
+ std::vector<double> A_nodecoord;
+ std::vector<double> A_nodeparametricCoord;
+ gmsh::model::mesh::getNodes(A_nodeTag, A_nodecoord, A_nodeparametricCoord, dim, tags[0], true);
+ std::cout << "-----------------------------------------\n";
+ std::cout << "point A displacement: u = " << full_d_vector[nodeIndexMap[A_nodeTag[0]]] << ", v = " << full_d_vector[nodeIndexMap[A_nodeTag[0]] + 1] << "\n";
+ }
}
diff --git a/srcs/FEM/uniform_charge.geo b/srcs/FEM/uniform_charge.geo
index 1b19a8ee9cef22b003fa155be8fb51442f9bf7ea..930e4410c83958dd777d5c1f19e42af13acb4e13 100644
--- a/srcs/FEM/uniform_charge.geo
+++ b/srcs/FEM/uniform_charge.geo
@@ -38,3 +38,5 @@ 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("Non_linear_solver", 1);
\ No newline at end of file
diff --git a/srcs/clampedMicroBeam.geo b/srcs/clampedMicroBeam.geo
index da4f0c45e65e7e57f358781ba04573d3d516dc7f..ddbf181966070c65ed1c270d4ec0db1718712afa 100644
--- a/srcs/clampedMicroBeam.geo
+++ b/srcs/clampedMicroBeam.geo
@@ -73,11 +73,13 @@ 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.);
+SetNumber("Non_linear_solver",1);
+
// BEM geometry
Physical Curve("BEM_FEM_boundary", 11) = {1,2,3};
Physical Curve("mass", 12) = {10}; //lower electrode
Physical Curve("rest_of_outside", 13) = {5,6,7,8,9};
-phi_top = 35;
+phi_top = 35.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);
diff --git a/srcs/coupling_validation.geo b/srcs/coupling_validation.geo
index 6a558ca2be560a6d4f8aa23e2baf8e43dad0150b..e5510d1e2ee35de5eebe1b186d0f0be9839fefe1 100644
--- a/srcs/coupling_validation.geo
+++ b/srcs/coupling_validation.geo
@@ -60,4 +60,6 @@ phi = 30;
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);
-SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
\ No newline at end of file
+SetNumber("Materials/BEM_domain_1/Epsilon", 8.8541878128e-12); // dielectric permittivity
+
+SetNumber("Non_linear_solver", 0);
\ No newline at end of file
diff --git a/srcs/longitudinal_comb.geo b/srcs/longitudinal_comb.geo
index 35aec87e3c670ffb3d33f36eee98cd9429209f66..6c89eec081168abdcadd2fbd9403a00c051875fd 100644
--- a/srcs/longitudinal_comb.geo
+++ b/srcs/longitudinal_comb.geo
@@ -141,4 +141,7 @@ 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
\ No newline at end of file
+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 711679cf5aff5fae0f0ea2453b7b057baf4948d9..90b59b98a2b0a7e9c02385c2057e14ed35a6be14 100644
--- a/srcs/main.cpp
+++ b/srcs/main.cpp
@@ -13,7 +13,6 @@ int main(int argc, char **argv)
return 0;
}
- bool nonLinearSolver = true; // use non-linear solver or not -> SET AS AN ARGUMENT OF MAIN, much easier
bool untangleMesh = false; // untangle geometry (deformed geometry, only available for nonLinear)
// If compiled with OpenMP support, gmsh::initialize
@@ -28,6 +27,9 @@ int main(int argc, char **argv)
Eigen::initParallel();
+ // determine if non linear solver must be set
+ bool nonLinearSolver = getNonLinearParameter();
+
if(nonLinearSolver)
{