source: issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp@ 12520

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

cosmetics
File size: 3.4 KB
Line 
1/*!\file: solver_thermal_nonlinear.cpp
2 * \brief: core of the thermal solution
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
11void solver_thermal_nonlinear(FemModel* femmodel){
12
13 /*solution : */
14 Vector* tg=NULL;
15 Vector* tf=NULL;
16 Vector* tf_old=NULL;
17 Vector* ys=NULL;
18 IssmDouble melting_offset;
19
20 /*intermediary: */
21 Matrix* Kff=NULL;
22 Matrix* Kfs=NULL;
23 Vector* pf=NULL;
24 Vector* df=NULL;
25
26 bool converged;
27 int constraints_converged;
28 int num_unstable_constraints;
29 int count;
30 int thermal_penalty_threshold;
31 int thermal_maxiter;
32
33 /*parameters:*/
34 int kflag,pflag;
35 bool lowmem=0;
36 int configuration_type;
37
38
39 /*Recover parameters: */
40 kflag=1; pflag=1;
41 femmodel->parameters->FindParam(&lowmem,SettingsLowmemEnum);
42 femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
43 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
44 femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
45
46 count=1;
47 converged=false;
48
49 if(VerboseSolution()) _pprintLine_("starting direct shooting method");
50 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
51 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum);
52 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
53
54 for(;;){
55
56 SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
57 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
58 Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf);
59 Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
60 xdelete(&tf_old); tf_old=tf->Duplicate();
61 xdelete(&Kff);xdelete(&pf);xdelete(&tg); xdelete(&df);
62 Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys);
63 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg);
64
65 ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
66
67 if (!converged){
68 if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints);
69 if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
70 if (count>=thermal_maxiter){
71 converged=true;
72 _pprintLine_(" maximum number of iterations (" << thermal_maxiter << ") exceeded");
73 }
74 }
75 count++;
76
77 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
78
79 if(converged)break;
80 }
81
82 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg);
83 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,melting_offset,MeltingOffsetEnum);
84
85 /*Free ressources: */
86 xdelete(&tg);
87 xdelete(&tf);
88 xdelete(&tf_old);
89 xdelete(&ys);
90}
Note: See TracBrowser for help on using the repository browser.