source: issm/oecreview/Archive/16554-17801/ISSM-16775-16776.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 5.6 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
2===================================================================
3--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16775)
4+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 16776)
5@@ -26,6 +26,7 @@
6 void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
7 void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
8 void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
9+ void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
10 void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
11 };
12 #endif
13Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
14===================================================================
15--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16775)
16+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16776)
17@@ -952,7 +952,7 @@
18 InputUpdateFromSolutionHOFS(solution,element);
19 return;
20 case SSAFSApproximationEnum:
21- _error_("not implemented yet");
22+ InputUpdateFromSolutionSSAFS(solution,element);
23 return;
24 default:
25 _error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
26@@ -1423,6 +1423,112 @@
27 xDelete<int>(doflist);
28 if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
29 }/*}}}*/
30+void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
31+
32+ int i;
33+ IssmDouble rho_ice,g,FSreconditioning;
34+ int* doflistSSA = NULL;
35+ int* doflistFSv = NULL;
36+ int* doflistFSp = NULL;
37+
38+ /*we have to add results of this element for FS and results from the element
39+ * at base for SSA, so we need to have the pointer toward the basal element*/
40+ Element* basalelement=element->GetBasalElement();
41+ if(basalelement->ObjectEnum()!=PentaEnum){
42+ _error_("Coupling not supported for "<<EnumToStringx(basalelement->ObjectEnum()));
43+ }
44+
45+ /*Fetch number of nodes and dof for this finite element*/
46+ int numnodes = 6;
47+ int numdof2d = numnodes;
48+ int numdofSSA = 6*2;
49+ int numdofFSv = 6*3;
50+ int numdofFSp = 6;
51+
52+ /*Fetch dof list and allocate solution vectors*/
53+ element->GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
54+ element->GetDofListPressure(&doflistFSp,GsetEnum);
55+ basalelement->GetDofList(&doflistSSA, SSAApproximationEnum, GsetEnum);
56+ IssmDouble* SSAvalues = xNew<IssmDouble>(numdofSSA);
57+ IssmDouble* FSvalues = xNew<IssmDouble>(numdofFSv+numdofFSp);
58+ IssmDouble* vx = xNew<IssmDouble>(numnodes);
59+ IssmDouble* vy = xNew<IssmDouble>(numnodes);
60+ IssmDouble* vz = xNew<IssmDouble>(numnodes);
61+ IssmDouble* vzSSA = xNew<IssmDouble>(numnodes);
62+ IssmDouble* vzFS = xNew<IssmDouble>(numnodes);
63+ IssmDouble* vel = xNew<IssmDouble>(numnodes);
64+ IssmDouble* pressure = xNew<IssmDouble>(numnodes);
65+
66+ /*Prepare coordinate system list*/
67+ int* cs_list = xNew<int>(2*numnodes);
68+ for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
69+ for(i=0;i<numnodes;i++) cs_list[numnodes+i] = PressureEnum;
70+
71+ /*Use the dof list to index into the solution vector: */
72+ element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
73+ for(i=0;i<numdof2d;i++){
74+ SSAvalues[i] = solution[doflistSSA[i]];
75+ SSAvalues[i+numdof2d] = solution[doflistSSA[i]];
76+ }
77+ for(i=0;i<numdofFSv;i++) FSvalues[i]=solution[doflistFSv[i]];
78+ for(i=0;i<numdofFSp;i++) FSvalues[numdofFSv+i]=solution[doflistFSp[i]];
79+
80+ /*Transform solution in Cartesian Space*/
81+ element->TransformSolutionCoord(FSvalues,2*numnodes,cs_list);
82+ element->TransformSolutionCoord(SSAvalues,numnodes,XYEnum);
83+
84+ /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
85+
86+ for(i=0;i<numnodes;i++){
87+ vx[i] = FSvalues[i*3+0]+SSAvalues[i*2+0];
88+ vy[i] = FSvalues[i*3+1]+SSAvalues[i*2+1];
89+ vzFS[i] = FSvalues[i*3+2];
90+ pressure[i] = FSvalues[numnodes*3+i]*FSreconditioning;
91+
92+ /*Check solution*/
93+ if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
94+ if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
95+ if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
96+ if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
97+ }
98+
99+ /*Get Vz and compute vel*/
100+ element->GetInputListOnVertices(vzSSA,VzSSAEnum);
101+ for(i=0;i<numnodes;i++){
102+ vz[i] = vzSSA[i]+vzFS[i];
103+ vel[i]= sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
104+ }
105+
106+ /*Now, we have to move the previous Vx and Vy inputs to old
107+ * status, otherwise, we'll wipe them off: */
108+ element->InputChangeName(VxEnum,VxPicardEnum);
109+ element->InputChangeName(VyEnum,VyPicardEnum);
110+ element->InputChangeName(VzEnum,VzPicardEnum);
111+ element->InputChangeName(PressureEnum,PressurePicardEnum);
112+
113+ /*Add vx and vy as inputs to element: */
114+ element->AddInput(VxEnum,vx,P1Enum);
115+ element->AddInput(VyEnum,vy,P1Enum);
116+ element->AddInput(VzEnum,vz,P1Enum);
117+ element->AddInput(VzFSEnum,vzFS,P1Enum);
118+ element->AddInput(VelEnum,vel,P1Enum);
119+ element->AddInput(PressureEnum,pressure,P1Enum);
120+
121+ /*Free ressources:*/
122+ xDelete<IssmDouble>(pressure);
123+ xDelete<IssmDouble>(vel);
124+ xDelete<IssmDouble>(vz);
125+ xDelete<IssmDouble>(vzSSA);
126+ xDelete<IssmDouble>(vzFS);
127+ xDelete<IssmDouble>(vy);
128+ xDelete<IssmDouble>(vx);
129+ xDelete<IssmDouble>(FSvalues);
130+ xDelete<IssmDouble>(SSAvalues);
131+ xDelete<int>(doflistFSp);
132+ xDelete<int>(doflistFSv);
133+ xDelete<int>(doflistSSA);
134+ xDelete<int>(cs_list);
135+}/*}}}*/
136 void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
137
138 int i,meshtype;
Note: See TracBrowser for help on using the repository browser.