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 | //save pointer to old velocity
|
---|
50 | VecFree(&old_ug);old_ug=ug;
|
---|
51 | VecFree(&old_uf);old_uf=uf;
|
---|
52 |
|
---|
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);
|
---|
57 |
|
---|
58 | convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
|
---|
59 |
|
---|
60 | CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
61 |
|
---|
62 | Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);
|
---|
63 | InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
|
---|
64 | SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
|
---|
65 | VecDuplicate(pf,&pJf);
|
---|
66 | MatMultPatch(Kff,uf,pJf);
|
---|
67 | VecAXPY(pJf,-1.,pf);
|
---|
68 | Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
|
---|
69 | MatFree(&Kff);VecFree(&pf); VecFree(&df);
|
---|
70 | VecAXPY(uf,-1.,duf);
|
---|
71 | Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
|
---|
72 |
|
---|
73 | InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
|
---|
74 |
|
---|
75 | /*Increase count: */
|
---|
76 | count++;
|
---|
77 | if(converged==true)break;
|
---|
78 | if(count>=max_nonlinear_iterations){
|
---|
79 | _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
|
---|
80 | break;
|
---|
81 | }
|
---|
82 | }
|
---|
83 |
|
---|
84 | _printf_(VerboseConvergence(),"\n total number of iterations: %i\n",count-1);
|
---|
85 |
|
---|
86 | /*clean-up*/
|
---|
87 | VecFree(&uf);
|
---|
88 | VecFree(&ug);
|
---|
89 | VecFree(&old_ug);
|
---|
90 | VecFree(&old_uf);
|
---|
91 | }
|
---|