source: issm/trunk/src/c/parallel/ProcessResults.cpp@ 675

Last change on this file since 675 was 675, checked in by Mathieu Morlighem, 16 years ago

fixed output from diagnostic (2,3 or 4 dofs)
File size: 8.6 KB
Line 
1/*!\file: ProcessResults.cpp
2 * \brief: go through results dataset, and for each result, process it for easier retrieval
3 * by the Matlab side. This usually means splitting the velocities from the g-size nodeset
4 * to the grid set (ug->vx,vy,vz), same for pressure (p_g->pressure), etc ... It also implies
5 * departitioning of the results.
6 * This phase is necessary prior to outputting the results on disk.
7 */
8
9#ifdef HAVE_CONFIG_H
10 #include "config.h"
11#else
12#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
13#endif
14
15#undef __FUNCT__
16#define __FUNCT__ "ProcessResults"
17
18#include "../DataSet/DataSet.h"
19#include "../objects/objects.h"
20#include "../EnumDefinitions/EnumDefinitions.h"
21#include "../shared/shared.h"
22
23void ProcessResults(DataSet** presults,FemModel* fems,int analysis_type){
24
25 int i,n;
26 Result* result=NULL;
27 Result* newresult=NULL;
28
29 /*input: */
30 DataSet* results=NULL;
31
32 /*output: */
33 DataSet* newresults=NULL;
34
35 /*fem diagnostic models: */
36 FemModel* fem_dh=NULL;
37 FemModel* fem_dv=NULL;
38 FemModel* fem_dhu=NULL;
39 FemModel* fem_ds=NULL;
40 FemModel* fem_sl=NULL;
41
42 /*fem prognostic models: */
43 FemModel* fem_p=NULL;
44
45 int ishutter;
46 int ismacayealpattyn;
47 int isstokes;
48 int dim;
49
50 /*intermediary: */
51 Vec u_g=NULL;
52 double* u_g_serial=NULL;
53 double* vx=NULL;
54 double* vy=NULL;
55 double* vz=NULL;
56 double* vel=NULL;
57 Vec p_g=NULL;
58 double* p_g_serial=NULL;
59 double* pressure=NULL;
60 double* parameter=NULL;
61 double* partition=NULL;
62 double yts;
63
64 Vec h_g=NULL;
65 double* h_g_serial=NULL;
66 double* thickness=NULL;
67
68 int numberofnodes;
69
70 /*recover input results: */
71 results=*presults;
72
73 /*Initialize new results: */
74 newresults=new DataSet(ResultsEnum());
75
76 /*Recover femmodels first: */
77 if(analysis_type==DiagnosticAnalysisEnum()){
78
79 fem_dh=fems+0;
80 fem_dv=fems+1;
81 fem_ds=fems+2;
82 fem_dhu=fems+3;
83 fem_sl=fems+4;
84
85 /*some flags needed: */
86 fem_dh->parameters->FindParam((void*)&dim,"dim");
87 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
88 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
89 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
90 }
91
92 if(analysis_type==PrognosticAnalysisEnum()){
93 fem_p=fems+0;
94 }
95
96 for(n=0;n<results->Size();n++){
97
98 result=(Result*)results->GetObjectByOffset(n);
99
100 if(strcmp(result->GetFieldName(),"u_g")==0){
101 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
102 *Stokes on 4 dofs: */
103 result->GetField(&u_g);
104 VecToMPISerial(&u_g_serial,u_g);
105
106 //2d results -> 2 dofs per node
107 if (dim==2){
108 /*ok, 2 dofs, on number of nodes: */
109 if(ismacayealpattyn){
110 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
111 VecToMPISerial(&partition,fem_dh->partition);
112 fem_dh->parameters->FindParam((void*)&yts,"yts");
113 }
114 else{
115 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
116 VecToMPISerial(&partition,fem_dhu->partition);
117 fem_dhu->parameters->FindParam((void*)&yts,"yts");
118 }
119 vx=(double*)xmalloc(numberofnodes*sizeof(double));
120 vy=(double*)xmalloc(numberofnodes*sizeof(double));
121 vz=(double*)xmalloc(numberofnodes*sizeof(double));
122 vel=(double*)xmalloc(numberofnodes*sizeof(double));
123
124 for(i=0;i<numberofnodes;i++){
125 vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
126 vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
127 vz[i]=0;
128 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
129 }
130 }
131
132 //3d results -> 3 or 4 (stokes) dofs per node
133 else{
134
135 if(!isstokes){
136 /*ok, 3 dofs, on number of nodes: */
137 if(ismacayealpattyn){
138 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
139 VecToMPISerial(&partition,fem_dh->partition);
140 fem_dh->parameters->FindParam((void*)&yts,"yts");
141 }
142 else{
143 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
144 VecToMPISerial(&partition,fem_dhu->partition);
145 fem_dhu->parameters->FindParam((void*)&yts,"yts");
146 }
147 vx=(double*)xmalloc(numberofnodes*sizeof(double));
148 vy=(double*)xmalloc(numberofnodes*sizeof(double));
149 vz=(double*)xmalloc(numberofnodes*sizeof(double));
150 vel=(double*)xmalloc(numberofnodes*sizeof(double));
151
152 for(i=0;i<numberofnodes;i++){
153 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
154 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
155 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
156 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
157 }
158 }
159
160 else{
161 /* 4 dofs on number of nodes. discard pressure: */
162 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
163 VecToMPISerial(&partition,fem_ds->partition);
164 fem_ds->parameters->FindParam((void*)&yts,"yts");
165 vx=(double*)xmalloc(numberofnodes*sizeof(double));
166 vy=(double*)xmalloc(numberofnodes*sizeof(double));
167 vz=(double*)xmalloc(numberofnodes*sizeof(double));
168 vel=(double*)xmalloc(numberofnodes*sizeof(double));
169 for(i=0;i<numberofnodes;i++){
170 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
171 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
172 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
173 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
174 }
175 }
176 }
177
178 /*Ok, add vx,vy and vz to newresults: */
179 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
180 newresults->AddObject(newresult);
181
182 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
183 newresults->AddObject(newresult);
184
185 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
186 newresults->AddObject(newresult);
187
188 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
189 newresults->AddObject(newresult);
190
191 /*do some cleanup: */
192 xfree((void**)&u_g_serial);
193 xfree((void**)&partition);
194 }
195
196 else if(strcmp(result->GetFieldName(),"p_g")==0){
197 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
198 result->GetField(&p_g);
199 VecToMPISerial(&p_g_serial,p_g);
200
201 if(!isstokes){
202 if(ismacayealpattyn){
203 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
204 VecToMPISerial(&partition,fem_dh->partition);
205 }
206 else{
207 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
208 VecToMPISerial(&partition,fem_dhu->partition);
209 }
210 }
211 else{
212 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
213 VecToMPISerial(&partition,fem_ds->partition);
214 }
215
216 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
217
218 for(i=0;i<numberofnodes;i++){
219 pressure[i]=p_g_serial[(int)partition[i]];
220 }
221
222 /*Ok, add pressure,vy and vz to newresults: */
223 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
224 newresults->AddObject(newresult);
225
226 /*do some cleanup: */
227 xfree((void**)&p_g_serial);
228 xfree((void**)&partition);
229 }
230
231 else if(strcmp(result->GetFieldName(),"h_g")==0){
232 /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
233 result->GetField(&h_g);
234 VecToMPISerial(&h_g_serial,h_g);
235 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
236 VecToMPISerial(&partition,fem_p->partition);
237
238 thickness=(double*)xmalloc(numberofnodes*sizeof(double));
239
240 for(i=0;i<numberofnodes;i++){
241 thickness[i]=h_g_serial[(int)partition[i]];
242 }
243
244 /*Ok, add pressure to newresults: */
245 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
246 newresults->AddObject(newresult);
247
248 /*do some cleanup: */
249 xfree((void**)&h_g_serial);
250 xfree((void**)&partition);
251 }
252
253 else if(strcmp(result->GetFieldName(),"param_g")==0){
254 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
255 result->GetField(&p_g);
256 VecToMPISerial(&p_g_serial,p_g);
257 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
258 VecToMPISerial(&partition,fem_p->partition);
259
260 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
261
262 for(i=0;i<numberofnodes;i++){
263 parameter[i]=p_g_serial[(int)partition[i]];
264 }
265
266 /*Ok, add parameter to newresults: */
267 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
268 newresults->AddObject(newresult);
269
270 /*do some cleanup: */
271 xfree((void**)&p_g_serial);
272 xfree((void**)&partition);
273 }
274
275 else{
276 /*Just copy the result into the new results dataset: */
277 newresults->AddObject(result);
278 }
279 }
280
281 /*Delete results: */
282 delete results;
283
284
285 /*Assign output pointers:*/
286 *presults=newresults;
287}
288
Note: See TracBrowser for help on using the repository browser.