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

Last change on this file since 928 was 928, checked in by seroussi, 16 years ago

missing fem_t in transient for 3d
File size: 12.5 KB
RevLine 
[643]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
[659]23void ProcessResults(DataSet** presults,FemModel* fems,int analysis_type){
[643]24
25 int i,n;
26 Result* result=NULL;
27 Result* newresult=NULL;
[659]28
29 /*input: */
30 DataSet* results=NULL;
[643]31
32 /*output: */
33 DataSet* newresults=NULL;
34
[667]35 /*fem diagnostic models: */
[643]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
[694]42 /*fem thermal models: */
43 FemModel* fem_t=NULL;
44
[667]45 /*fem prognostic models: */
46 FemModel* fem_p=NULL;
47
[758]48 /*fem control models: */
49 FemModel* fem_c=NULL;
50
[694]51 /*some parameters*/
[643]52 int ishutter;
53 int ismacayealpattyn;
54 int isstokes;
[675]55 int dim;
[643]56
57 /*intermediary: */
58 Vec u_g=NULL;
59 double* u_g_serial=NULL;
60 double* vx=NULL;
61 double* vy=NULL;
62 double* vz=NULL;
[659]63 double* vel=NULL;
[643]64 Vec p_g=NULL;
65 double* p_g_serial=NULL;
66 double* pressure=NULL;
67 double* partition=NULL;
68 double yts;
69
[758]70 double* param_g=NULL;
71 double* parameter=NULL;
72
[853]73 Vec s_g=NULL;
74 double* s_g_serial=NULL;
75 double* surface=NULL;
76
77 Vec b_g=NULL;
78 double* b_g_serial=NULL;
79 double* bed=NULL;
80
[667]81 Vec h_g=NULL;
82 double* h_g_serial=NULL;
83 double* thickness=NULL;
84
[694]85 Vec t_g=NULL;
86 double* t_g_serial=NULL;
87 double* temperature=NULL;
88
89 Vec m_g=NULL;
90 double* m_g_serial=NULL;
91 double* melting=NULL;
92
[643]93 int numberofnodes;
94
[659]95 /*recover input results: */
96 results=*presults;
97
[643]98 /*Initialize new results: */
99 newresults=new DataSet(ResultsEnum());
100
101 /*Recover femmodels first: */
102 if(analysis_type==DiagnosticAnalysisEnum()){
103
104 fem_dh=fems+0;
105 fem_dv=fems+1;
106 fem_ds=fems+2;
107 fem_dhu=fems+3;
108 fem_sl=fems+4;
109
110 /*some flags needed: */
[675]111 fem_dh->parameters->FindParam((void*)&dim,"dim");
[643]112 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
113 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
114 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
115 }
116
[853]117 if(analysis_type==TransientAnalysisEnum()){
118
119 fem_dh=fems+0;
120 fem_dv=fems+1;
121 fem_ds=fems+2;
122 fem_dhu=fems+3;
123 fem_sl=fems+4;
124 fem_p=fems+5;
125
126 /*some flags needed: */
127 fem_dh->parameters->FindParam((void*)&dim,"dim");
128 fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
129 fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
130 fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
[928]131
132 if (dim==3){
133 fem_t=fems+6;
134 }
[853]135 }
136
[667]137 if(analysis_type==PrognosticAnalysisEnum()){
138 fem_p=fems+0;
139 }
[643]140
[694]141 if(analysis_type==ThermalAnalysisEnum()){
142 fem_t=fems+0;
143 }
144
[758]145 if(analysis_type==ControlAnalysisEnum()){
146 fem_dh=fems+0; //load u_g
147 fem_c=fems+0; //load param_g
148 }
149
[643]150 for(n=0;n<results->Size();n++){
151 result=(Result*)results->GetObjectByOffset(n);
[853]152
[643]153 if(strcmp(result->GetFieldName(),"u_g")==0){
[853]154
[675]155 /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or
[643]156 *Stokes on 4 dofs: */
157 result->GetField(&u_g);
158 VecToMPISerial(&u_g_serial,u_g);
159
[675]160 //2d results -> 2 dofs per node
161 if (dim==2){
[643]162 /*ok, 2 dofs, on number of nodes: */
163 if(ismacayealpattyn){
164 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
165 VecToMPISerial(&partition,fem_dh->partition);
166 fem_dh->parameters->FindParam((void*)&yts,"yts");
167 }
168 else{
169 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
170 VecToMPISerial(&partition,fem_dhu->partition);
171 fem_dhu->parameters->FindParam((void*)&yts,"yts");
172 }
173 vx=(double*)xmalloc(numberofnodes*sizeof(double));
174 vy=(double*)xmalloc(numberofnodes*sizeof(double));
175 vz=(double*)xmalloc(numberofnodes*sizeof(double));
[659]176 vel=(double*)xmalloc(numberofnodes*sizeof(double));
[643]177
178 for(i=0;i<numberofnodes;i++){
179 vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
180 vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
181 vz[i]=0;
[659]182 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
[643]183 }
[675]184 }
185 //3d results -> 3 or 4 (stokes) dofs per node
[643]186 else{
[675]187 if(!isstokes){
188 /*ok, 3 dofs, on number of nodes: */
189 if(ismacayealpattyn){
190 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
191 VecToMPISerial(&partition,fem_dh->partition);
192 fem_dh->parameters->FindParam((void*)&yts,"yts");
193 }
194 else{
195 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
196 VecToMPISerial(&partition,fem_dhu->partition);
197 fem_dhu->parameters->FindParam((void*)&yts,"yts");
198 }
199 vx=(double*)xmalloc(numberofnodes*sizeof(double));
200 vy=(double*)xmalloc(numberofnodes*sizeof(double));
201 vz=(double*)xmalloc(numberofnodes*sizeof(double));
202 vel=(double*)xmalloc(numberofnodes*sizeof(double));
203
204 for(i=0;i<numberofnodes;i++){
205 vx[i]=u_g_serial[3*(int)partition[i]+0]*yts;
206 vy[i]=u_g_serial[3*(int)partition[i]+1]*yts;
207 vz[i]=u_g_serial[3*(int)partition[i]+2]*yts;
208 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
209 }
[643]210 }
[675]211 else{
212 /* 4 dofs on number of nodes. discard pressure: */
213 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
214 VecToMPISerial(&partition,fem_ds->partition);
215 fem_ds->parameters->FindParam((void*)&yts,"yts");
216 vx=(double*)xmalloc(numberofnodes*sizeof(double));
217 vy=(double*)xmalloc(numberofnodes*sizeof(double));
218 vz=(double*)xmalloc(numberofnodes*sizeof(double));
219 vel=(double*)xmalloc(numberofnodes*sizeof(double));
220 for(i=0;i<numberofnodes;i++){
221 vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
222 vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
223 vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
224 vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
225 }
226 }
[643]227 }
228
229 /*Ok, add vx,vy and vz to newresults: */
230 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
231 newresults->AddObject(newresult);
232
233 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
234 newresults->AddObject(newresult);
235
236 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
237 newresults->AddObject(newresult);
238
[659]239 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
240 newresults->AddObject(newresult);
241
[643]242 /*do some cleanup: */
243 xfree((void**)&u_g_serial);
244 xfree((void**)&partition);
245 }
246 else if(strcmp(result->GetFieldName(),"p_g")==0){
247 /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
248 result->GetField(&p_g);
249 VecToMPISerial(&p_g_serial,p_g);
250
251 if(!isstokes){
252 if(ismacayealpattyn){
253 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
254 VecToMPISerial(&partition,fem_dh->partition);
255 }
256 else{
257 fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
258 VecToMPISerial(&partition,fem_dhu->partition);
259 }
260 }
261 else{
262 fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
263 VecToMPISerial(&partition,fem_ds->partition);
264 }
265
266 pressure=(double*)xmalloc(numberofnodes*sizeof(double));
267
268 for(i=0;i<numberofnodes;i++){
269 pressure[i]=p_g_serial[(int)partition[i]];
270 }
271
272 /*Ok, add pressure,vy and vz to newresults: */
273 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
274 newresults->AddObject(newresult);
275
276 /*do some cleanup: */
277 xfree((void**)&p_g_serial);
278 xfree((void**)&partition);
279 }
[694]280 else if(strcmp(result->GetFieldName(),"t_g")==0){
281 /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
282 result->GetField(&t_g);
283 VecToMPISerial(&t_g_serial,t_g);
284 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
285 VecToMPISerial(&partition,fem_t->partition);
286
287 temperature=(double*)xmalloc(numberofnodes*sizeof(double));
288
289 for(i=0;i<numberofnodes;i++){
290 temperature[i]=t_g_serial[(int)partition[i]];
291 }
292
293 /*Ok, add pressure to newresults: */
294 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes);
295 newresults->AddObject(newresult);
296
297 /*do some cleanup: */
298 xfree((void**)&t_g_serial);
299 xfree((void**)&partition);
300 }
301 else if(strcmp(result->GetFieldName(),"m_g")==0){
302 /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
303 result->GetField(&m_g);
304 VecToMPISerial(&m_g_serial,m_g);
305 fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
[695]306 fem_t->parameters->FindParam((void*)&yts,"yts");
[694]307 VecToMPISerial(&partition,fem_t->partition);
308
309 melting=(double*)xmalloc(numberofnodes*sizeof(double));
310
311 for(i=0;i<numberofnodes;i++){
[695]312 melting[i]=m_g_serial[(int)partition[i]]*yts;
[694]313 }
314
315 /*Ok, add pressure to newresults: */
316 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes);
317 newresults->AddObject(newresult);
318
319 /*do some cleanup: */
320 xfree((void**)&m_g_serial);
321 xfree((void**)&partition);
322 }
[667]323 else if(strcmp(result->GetFieldName(),"h_g")==0){
324 /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
325 result->GetField(&h_g);
326 VecToMPISerial(&h_g_serial,h_g);
327 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
328 VecToMPISerial(&partition,fem_p->partition);
329
330 thickness=(double*)xmalloc(numberofnodes*sizeof(double));
331
332 for(i=0;i<numberofnodes;i++){
333 thickness[i]=h_g_serial[(int)partition[i]];
334 }
335
336 /*Ok, add pressure to newresults: */
337 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofnodes);
338 newresults->AddObject(newresult);
339
340 /*do some cleanup: */
341 xfree((void**)&h_g_serial);
342 xfree((void**)&partition);
343 }
[853]344 else if(strcmp(result->GetFieldName(),"s_g")==0){
345 /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */
346 result->GetField(&s_g);
347 VecToMPISerial(&s_g_serial,s_g);
348 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
349 VecToMPISerial(&partition,fem_p->partition);
350
351 surface=(double*)xmalloc(numberofnodes*sizeof(double));
352
353 for(i=0;i<numberofnodes;i++){
354 surface[i]=s_g_serial[(int)partition[i]];
355 }
356
357 /*Ok, add pressure to newresults: */
358 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes);
359 newresults->AddObject(newresult);
360
361 /*do some cleanup: */
362 xfree((void**)&s_g_serial);
363 xfree((void**)&partition);
364 }
365 else if(strcmp(result->GetFieldName(),"b_g")==0){
366 /*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */
367 result->GetField(&b_g);
368 VecToMPISerial(&b_g_serial,b_g);
369 fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
370 VecToMPISerial(&partition,fem_p->partition);
371
372 bed=(double*)xmalloc(numberofnodes*sizeof(double));
373
374 for(i=0;i<numberofnodes;i++){
375 bed[i]=b_g_serial[(int)partition[i]];
376 }
377
378 /*Ok, add pressure to newresults: */
379 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes);
380 newresults->AddObject(newresult);
381
382 /*do some cleanup: */
383 xfree((void**)&b_g_serial);
384 xfree((void**)&partition);
385 }
[669]386 else if(strcmp(result->GetFieldName(),"param_g")==0){
387 /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
[758]388 result->GetField(&param_g);
389 fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
390 VecToMPISerial(&partition,fem_dh->partition);
[669]391
392 parameter=(double*)xmalloc(numberofnodes*sizeof(double));
393
394 for(i=0;i<numberofnodes;i++){
[758]395 parameter[i]=param_g[2*(int)partition[i]];
[669]396 }
397
398 /*Ok, add parameter to newresults: */
399 newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes);
400 newresults->AddObject(newresult);
401
402 /*do some cleanup: */
403 xfree((void**)&partition);
404 }
405 else{
406 /*Just copy the result into the new results dataset: */
[758]407 newresults->AddObject(result->copy());
[669]408 }
[643]409 }
410
[659]411 /*Delete results: */
412 delete results;
413
414
[643]415 /*Assign output pointers:*/
[659]416 *presults=newresults;
[643]417}
Note: See TracBrowser for help on using the repository browser.