source: issm/trunk-jpl/src/c/solvers/solver_newton.cpp@ 11324

Last change on this file since 11324 was 11324, checked in by Mathieu Morlighem, 13 years ago

Fixed many memory leaks
File size: 3.7 KB
RevLine 
[11322]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
13void 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
[11324]52 /*Solver forward model*/
[11322]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);
[11324]56 Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);VecFree(&df);
[11322]57 Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
[11324]58 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);VecFree(&ug);
[11322]59
[11324]60 /*Check convergence*/
61 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
62 MatFree(&Kff);VecFree(&pf);
[11322]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 }
[11324]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);
[11322]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}
Note: See TracBrowser for help on using the repository browser.