source: issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp@ 15002

Last change on this file since 15002 was 15002, checked in by Eric.Larour, 12 years ago

CHG: solutionsequences should include shared/shared.h and classes/classes.h only
File size: 4.6 KB
Line 
1/*!\file: solutionsequence_nonlinear.cpp
2 * \brief: core of a non-linear solution, using fixed-point method
3 */
4
5#include "../toolkits/toolkits.h"
6#include "../classes/classes.h"
7#include "../shared/shared.h"
8#include "../modules/modules.h"
9#include "../analyses/analyses.h"
10
11void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads){
12
13 /*intermediary: */
14 Matrix<IssmDouble>* Kff = NULL;
15 Matrix<IssmDouble>* Kfs = NULL;
16 Vector<IssmDouble>* ug = NULL;
17 Vector<IssmDouble>* old_ug = NULL;
18 Vector<IssmDouble>* uf = NULL;
19 Vector<IssmDouble>* old_uf = NULL;
20 Vector<IssmDouble>* pf = NULL;
21 Vector<IssmDouble>* df = NULL;
22 Vector<IssmDouble>* ys = NULL;
23
24 Loads* savedloads=NULL;
25 bool converged;
26 int constraints_converged;
27 int num_unstable_constraints;
28 int count;
29
30 /*parameters:*/
31 int min_mechanical_constraints;
32 int max_nonlinear_iterations;
33 int configuration_type;
34
35 /*Recover parameters: */
36 femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum);
37 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
38 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
39 femmodel->UpdateConstraintsx();
40
41 /*Were loads requested as output? : */
42 if(conserve_loads){
43 savedloads = static_cast<Loads*>(femmodel->loads->Copy());
44 }
45
46 count=1;
47 converged=false;
48
49 /*Start non-linear iteration using input velocity: */
50 GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
51 Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
52
53 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
54 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
55 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
56
57 for(;;){
58
59 //save pointer to old velocity
60 delete old_ug;old_ug=ug;
61 delete old_uf;old_uf=uf;
62
63 femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
64 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
65 Reduceloadx(pf, Kfs, ys); delete Kfs;
66 Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
67 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
68
69 convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); delete Kff; delete pf; delete df;
70 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
71 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
72
73 ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
74 if(VerboseConvergence()) _pprintLine_(" number of unstable constraints: " << num_unstable_constraints);
75
76 //rift convergence
77 if (!constraints_converged) {
78 if (converged){
79 if (num_unstable_constraints <= min_mechanical_constraints) converged=true;
80 else converged=false;
81 }
82 }
83
84 /*Increase count: */
85 count++;
86 if(converged==true){
87 bool max_iteration_state=false;
88 int tempStep=1;
89 IssmDouble tempTime=1.0;
90 femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
91 break;
92 }
93 if(count>=max_nonlinear_iterations){
94 _pprintLine_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded");
95 converged=true;
96 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
97 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
98 bool max_iteration_state=true;
99 int tempStep=1;
100 IssmDouble tempTime=1.0;
101 femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime));
102 break;
103 }
104 }
105
106 if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1);
107
108 /*clean-up*/
109 if(conserve_loads){
110 delete femmodel->loads;
111 femmodel->loads=savedloads;
112 }
113 delete uf;
114 delete ug;
115 delete old_ug;
116 delete old_uf;
117}
Note: See TracBrowser for help on using the repository browser.