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

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

Large commit, to strip out Petsc from the code. Created new Matrix and Vector
wrapper objects to the Petsc Mat and Vec objects. These wrappers make it possible
to hook up a different package to handle matrix and vector operations.

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 double 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 _printf_(VerboseSolution(),"%s\n","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 _printf_(VerboseConvergence(),"%s%i\n"," #unstable constraints = ",num_unstable_constraints);
69 if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
70 if (count>=thermal_maxiter){
71 converged=true;
72 _printf_(true," maximum number of iterations (%i) exceeded\n",thermal_maxiter);
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.