source: issm/branches/trunk-jpl-damage/src/c/solvers/solver_stokescoupling_nonlinear.cpp@ 11684

Last change on this file since 11684 was 11684, checked in by cborstad, 13 years ago

merged revisions 11428:11680 from trunk-jpl into branches/trunk-jpl-damage
File size: 4.2 KB
Line 
1/*!\file: solver_stokescoupling_nonlinear.cpp
2 * \brief: core of the coupling between stokes and macayealpattyn
3 */
4
5#include "../toolkits/toolkits.h"
6#include "../objects/objects.h"
7#include "../EnumDefinitions/EnumDefinitions.h"
8#include "../io/io.h"
9#include "../modules/modules.h"
10#include "../solutions/solutions.h"
11#include "./solvers.h"
12
13void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){
14
15 /*intermediary: */
16 Matrix* Kff_horiz = NULL;
17 Matrix* Kfs_horiz = NULL;
18 Vector* ug_horiz = NULL;
19 Vector* uf_horiz = NULL;
20 Vector* old_uf_horiz = NULL;
21 Vector* pf_horiz = NULL;
22 Vector* df_horiz = NULL;
23 Matrix* Kff_vert = NULL;
24 Matrix* Kfs_vert = NULL;
25 Vector* ug_vert = NULL;
26 Vector* uf_vert = NULL;
27 Vector* pf_vert = NULL;
28 Vector* df_vert = NULL;
29 Vector* ys = NULL;
30 bool converged;
31 int constraints_converged;
32 int num_unstable_constraints;
33 int count;
34
35 /*parameters:*/
36 int min_mechanical_constraints;
37 int max_nonlinear_iterations;
38 int configuration_type;
39
40 /*Recover parameters: */
41 femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum);
42 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
43 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
44
45 count=1;
46 converged=false;
47
48 /*First get ug_horiz:*/
49 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
50 GetSolutionFromInputsx(&ug_horiz, femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters);
51 Reducevectorgtofx(&uf_horiz, ug_horiz, femmodel->nodes,femmodel->parameters);
52
53 for(;;){
54
55 /*First diagnostic horiz:*/
56 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
57 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
58
59 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
60 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
61 xdelete(&ug_horiz);
62
63 //save pointer to old velocity
64 xdelete(&old_uf_horiz); old_uf_horiz=uf_horiz;
65
66 /*solve: */
67 SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
68 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
69 Reduceloadx(pf_horiz, Kfs_horiz, ys); xdelete(&Kfs_horiz);
70 Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters);
71 Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys);
72 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
73
74 convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); xdelete(&Kff_horiz); xdelete(&pf_horiz); xdelete(&df_horiz);
75
76 /*Second compute vertical velocity: */
77 femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
78 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
79
80 /*solve: */
81 SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
82 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
83 Reduceloadx(pf_vert, Kfs_vert, ys); xdelete(&Kfs_vert);
84 Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); xdelete(&Kff_vert); xdelete(&pf_vert); xdelete(&df_vert);
85 Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters);xdelete(&uf_vert); xdelete(&ys);
86 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert);
87 xdelete(&ug_vert); xdelete(&uf_vert);
88
89 /*Increase count: */
90 count++;
91 if(converged==true)break;
92 if(count>=max_nonlinear_iterations){
93 _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
94 break;
95 }
96 }
97
98 /*clean-up*/
99 xdelete(&old_uf_horiz);
100 xdelete(&uf_horiz);
101 xdelete(&ug_horiz);
102 xdelete(&ys);
103}
Note: See TracBrowser for help on using the repository browser.