1 | /*!\file: solver_nonlinear.cpp
|
---|
2 | * \brief: core of a non-linear solution, using fixed-point method
|
---|
3 | */
|
---|
4 |
|
---|
5 | #include "../toolkits/toolkits.h"
|
---|
6 | #include "../objects/objects.h"
|
---|
7 | #include "../io/io.h"
|
---|
8 | #include "../EnumDefinitions/EnumDefinitions.h"
|
---|
9 | #include "../modules/modules.h"
|
---|
10 | #include "../solutions/solutions.h"
|
---|
11 | #include "./solvers.h"
|
---|
12 |
|
---|
13 | void solver_newton(FemModel* femmodel){
|
---|
14 |
|
---|
15 | /*intermediary: */
|
---|
16 | Mat Kff = NULL, Kfs = NULL, Jff = NULL;
|
---|
17 | Vec ug = NULL, old_ug = NULL;
|
---|
18 | Vec uf = NULL, old_uf = NULL, duf = NULL;
|
---|
19 | Vec pf = NULL, pJf = NULL;
|
---|
20 | Vec df = NULL;
|
---|
21 | Vec ys = NULL;
|
---|
22 |
|
---|
23 | bool converged;
|
---|
24 | int num_unstable_constraints;
|
---|
25 | int count;
|
---|
26 |
|
---|
27 | /*parameters:*/
|
---|
28 | int max_nonlinear_iterations;
|
---|
29 | int configuration_type;
|
---|
30 |
|
---|
31 | /*Recover parameters: */
|
---|
32 | femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
|
---|
33 | femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
|
---|
34 | UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
|
---|
35 |
|
---|
36 | count=1;
|
---|
37 | converged=false;
|
---|
38 |
|
---|
39 | /*Start non-linear iteration using input velocity: */
|
---|
40 | GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
41 | Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);
|
---|
42 |
|
---|
43 | //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
|
---|
44 | InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
|
---|
45 | InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
|
---|
46 |
|
---|
47 | for(;;){
|
---|
48 |
|
---|
49 | VecFree(&old_ug);old_ug=ug;
|
---|
50 | VecFree(&old_uf);old_uf=uf;
|
---|
51 |
|
---|
52 | /*Solver forward model*/
|
---|
53 | SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
54 | CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
|
---|
55 | Reduceloadx(pf,Kfs,ys);MatFree(&Kfs);
|
---|
56 | Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);VecFree(&df);
|
---|
57 | Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
|
---|
58 | InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);VecFree(&ug);
|
---|
59 |
|
---|
60 | /*Check convergence*/
|
---|
61 | convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
|
---|
62 | MatFree(&Kff);VecFree(&pf);
|
---|
63 | count++;
|
---|
64 | if(converged==true)break;
|
---|
65 | if(count>=max_nonlinear_iterations){
|
---|
66 | _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
|
---|
67 | break;
|
---|
68 | }
|
---|
69 |
|
---|
70 | /*Prepare next iteration using Newton's method*/
|
---|
71 | CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
72 | SystemMatricesx(&Kff,NULL,&pf,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
73 | CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
|
---|
74 | VecDuplicate(pf,&pJf);
|
---|
75 | MatMultPatch(Kff,uf,pJf); MatFree(&Kff);
|
---|
76 | VecScale(pJf,-1.);
|
---|
77 | VecAXPY(pJf,+1.,pf); VecFree(&pf);
|
---|
78 | Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
|
---|
79 | MatFree(&Jff); VecFree(&pJf);
|
---|
80 | VecAXPY(uf,+1.,duf); VecFree(&duf);
|
---|
81 | Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
|
---|
82 | InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
|
---|
83 | }
|
---|
84 |
|
---|
85 | _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1);
|
---|
86 |
|
---|
87 | /*clean-up*/
|
---|
88 | VecFree(&uf);
|
---|
89 | VecFree(&ug);
|
---|
90 | VecFree(&old_ug);
|
---|
91 | VecFree(&old_uf);
|
---|
92 | }
|
---|