source: issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp@ 586

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

Added dakota parallel c code capability + Prognostic capability
File size: 17.3 KB
Line 
1/*
2 * CreateElementsNodesAndMaterialsPrognostic.c:
3 */
4
5#undef __FUNCT__
6#define __FUNCT__ "CreateElementsNodesAndMaterialsPrognostic"
7
8#include "../../DataSet/DataSet.h"
9#include "../../toolkits/toolkits.h"
10#include "../../EnumDefinitions/EnumDefinitions.h"
11#include "../../objects/objects.h"
12#include "../../shared/shared.h"
13#include "../../MeshPartitionx/MeshPartitionx.h"
14#include "../Model.h"
15
16
17void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
18
19
20 /*output: int* epart, int* my_grids, double* my_bordergrids*/
21
22
23 int i,j,k,n;
24 extern int my_rank;
25 extern int num_procs;
26
27 /*DataSets: */
28 DataSet* elements = NULL;
29 DataSet* nodes = NULL;
30 DataSet* materials = NULL;
31
32 /*Objects: */
33 Node* node = NULL;
34 Tria* tria = NULL;
35 Penta* penta = NULL;
36 Matice* matice = NULL;
37 Matpar* matpar = NULL;
38
39 int analysis_type;
40 int sub_analysis_type;
41
42 /*output: */
43 int* epart=NULL; //element partitioning.
44 int* npart=NULL; //node partitioning.
45 int* my_grids=NULL;
46 double* my_bordergrids=NULL;
47
48
49 /*intermediary: */
50 int elements_width; //size of elements
51 double B_avg;
52
53 /*tria constructor input: */
54 int tria_id;
55 int tria_mid;
56 int tria_mparid;
57 int tria_g[3];
58 double tria_h[3];
59 double tria_s[3];
60 double tria_b[3];
61 double tria_k[3];
62 double tria_melting[3];
63 double tria_accumulation[3];
64 double tria_geothermalflux[3];
65 int tria_friction_type;
66 double tria_p;
67 double tria_q;
68 int tria_shelf;
69 double tria_meanvel;/*!scaling ratio for velocities*/
70 double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
71 double tria_viscosity_overshoot;
72 int tria_artdiff;
73
74 /*matice constructor input: */
75 int matice_mid;
76 double matice_B;
77 double matice_n;
78
79 /*penta constructor input: */
80
81 int penta_id;
82 int penta_mid;
83 int penta_mparid;
84 int penta_g[6];
85 double penta_h[6];
86 double penta_s[6];
87 double penta_b[6];
88 double penta_k[6];
89 int penta_friction_type;
90 double penta_p;
91 double penta_q;
92 int penta_shelf;
93 int penta_onbed;
94 int penta_onsurface;
95 double penta_meanvel;/*!scaling ratio for velocities*/
96 double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
97 int penta_collapse;
98 double penta_melting[6];
99 double penta_accumulation[6];
100 double penta_geothermalflux[6];
101 int penta_artdiff;
102 int penta_thermal_steadystate;
103 double penta_viscosity_overshoot;
104 double penta_stokesreconditioning;
105
106 /*matpar constructor input: */
107 int matpar_mid;
108 double matpar_rho_ice;
109 double matpar_rho_water;
110 double matpar_heatcapacity;
111 double matpar_thermalconductivity;
112 double matpar_latentheat;
113 double matpar_beta;
114 double matpar_meltingpoint;
115 double matpar_mixed_layer_capacity;
116 double matpar_thermal_exchange_velocity;
117 double matpar_g;
118
119 /* node constructor input: */
120 int node_id;
121 int node_partitionborder=0;
122 double node_x[3];
123 int node_onbed;
124 int node_onsurface;
125 int node_upper_node_id;
126 int node_numdofs;
127
128
129 #ifdef _PARALLEL_
130 /*Metis partitioning: */
131 int range;
132 Vec gridborder=NULL;
133 int my_numgrids;
134 int* all_numgrids=NULL;
135 int gridcount;
136 int count;
137 #endif
138 int first_grid_index;
139
140 /*Rifts:*/
141 int* riftsnumpenaltypairs;
142 double** riftspenaltypairs;
143 int* riftsfill;
144 double* riftsfriction;
145 double* riftpenaltypairs=NULL;
146 int el1,el2;
147
148 /*First create the elements, nodes and material properties: */
149 elements = new DataSet(ElementsEnum());
150 nodes = new DataSet(NodesEnum());
151 materials = new DataSet(MaterialsEnum());
152
153
154 /*Get analysis_type: */
155 analysis_type=AnalysisTypeAsEnum(model->analysis_type);
156 sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
157
158
159 /*Width of elements: */
160 if(strcmp(model->meshtype,"2d")==0){
161 elements_width=3; //tria elements
162 }
163 else{
164 elements_width=6; //penta elements
165 }
166
167 #ifdef _PARALLEL_
168 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
169 if(strcmp(model->meshtype,"2d")==0){
170 /*load elements: */
171 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
172 }
173 else{
174 /*load elements2d: */
175 ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
176 }
177
178
179 MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
180
181 /*Free elements and elements2d: */
182 xfree((void**)&model->elements);
183 xfree((void**)&model->elements2d);
184
185
186 /*Deal with rifts, they have to be included into one partition only, not several: */
187 FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
188
189 for(i=0;i<model->numrifts;i++){
190 riftpenaltypairs=model->riftspenaltypairs[i];
191 for(j=0;j<model->riftsnumpenaltypairs[i];j++){
192 el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
193 el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
194 epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
195 }
196 }
197 /*Free rifts: */
198 xfree((void**)&riftsnumpenaltypairs);
199 for(i=0;i<model->numrifts;i++){
200 double* temp=riftspenaltypairs[i];
201 xfree((void**)&temp);
202 }
203 xfree((void**)&riftspenaltypairs);
204 xfree((void**)&riftsfill);
205 xfree((void**)&riftsfriction);
206
207 /*Used later on: */
208 my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
209 #endif
210
211
212
213 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
214
215 /*2d mesh: */
216 if (strcmp(model->meshtype,"2d")==0){
217
218 /*Fetch data needed: */
219 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
220 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
221 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
222 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
223 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
224
225 for (i=0;i<model->numberofelements;i++){
226
227 #ifdef _PARALLEL_
228 /*!All elements have been partitioned above, only create elements for this CPU: */
229 if(my_rank==epart[i]){
230 #endif
231
232
233 /*ids: */
234 tria_id=i+1; //matlab indexing.
235 tria_mid=-1; //no need for materials
236 tria_mparid=-1; //no need for materials
237
238 /*vertices offsets: */
239 tria_g[0]=(int)*(model->elements+elements_width*i+0);
240 tria_g[1]=(int)*(model->elements+elements_width*i+1);
241 tria_g[2]=(int)*(model->elements+elements_width*i+2);
242
243 /*thickness,surface and bed:*/
244 tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
245 tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1));
246 tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ;
247
248 tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1));
249 tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1));
250 tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1));
251
252 tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1));
253 tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1));
254 tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1));
255
256 /*element on iceshelf?:*/
257 tria_shelf=(int)*(model->elementoniceshelf+i);
258 tria_artdiff=model->artificial_diffusivity;
259
260 /*Create tria element using its constructor:*/
261 tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff);
262
263 /*Add tria element to elements dataset: */
264 elements->AddObject(tria);
265
266 #ifdef _PARALLEL_
267 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
268 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
269 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
270 will hold which grids belong to this partition*/
271 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
272 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
273 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
274 #endif
275
276 #ifdef _PARALLEL_
277 }//if(my_rank==epart[i])
278 #endif
279
280 }//for (i=0;i<numberofelements;i++)
281
282
283 /*Free data : */
284 xfree((void**)&model->elements);
285 xfree((void**)&model->thickness);
286 xfree((void**)&model->surface);
287 xfree((void**)&model->bed);
288 xfree((void**)&model->elementoniceshelf);
289
290 }
291 else{ // if (strcmp(meshtype,"2d")==0)
292
293 /*Fetch data needed: */
294 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
295 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
296 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
297 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
298 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
299 ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
300 ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
301
302 for (i=0;i<model->numberofelements;i++){
303 #ifdef _PARALLEL_
304 /*We are using our element partition to decide which elements will be created on this node: */
305 if(my_rank==epart[i]){
306 #endif
307
308
309 /*name and id: */
310 penta_id=i+1; //matlab indexing.
311 penta_mid=-1;
312 penta_mparid=-1; //no need for materials
313
314 /*vertices,thickness,surface,bed and drag: */
315 for(j=0;j<6;j++){
316 penta_g[j]=(int)*(model->elements+elements_width*i+j);
317 penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1));
318 penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1));
319 penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1));
320 }
321
322 /*diverse: */
323 penta_shelf=(int)*(model->elementoniceshelf+i);
324 penta_onbed=(int)*(model->elementonbed+i);
325 penta_onsurface=(int)*(model->elementonsurface+i);
326 penta_collapse=1;
327 penta_artdiff=model->artificial_diffusivity;
328
329
330 /*Create Penta using its constructor:*/
331 penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
332 penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
333 penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
334 penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning);
335
336 /*Add penta element to elements dataset: */
337 elements->AddObject(penta);
338
339
340 #ifdef _PARALLEL_
341 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
342 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
343 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
344 will hold which grids belong to this partition*/
345 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
346 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
347 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
348 my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
349 my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
350 my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
351 #endif
352
353 #ifdef _PARALLEL_
354 }//if(my_rank==epart[i])
355 #endif
356
357 }//for (i=0;i<numberofelements;i++)
358
359 /*Free data: */
360 xfree((void**)&model->elements);
361 xfree((void**)&model->thickness);
362 xfree((void**)&model->surface);
363 xfree((void**)&model->bed);
364 xfree((void**)&model->elementoniceshelf);
365 xfree((void**)&model->elementonbed);
366 xfree((void**)&model->elementonsurface);
367
368 } //if (strcmp(meshtype,"2d")==0)
369
370 #ifdef _PARALLEL_
371 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
372 *element partition, and which are on its border with other nodes:*/
373 gridborder=NewVec(model->numberofnodes);
374
375 for (i=0;i<model->numberofnodes;i++){
376 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
377 }
378 VecAssemblyBegin(gridborder);
379 VecAssemblyEnd(gridborder);
380
381 #ifdef _DEBUG_
382 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
383 #endif
384
385 VecToMPISerial(&my_bordergrids,gridborder);
386
387 #ifdef _DEBUG_
388 if(my_rank==0){
389 for (i=0;i<model->numberofnodes;i++){
390 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
391 }
392 }
393 #endif
394 #endif
395
396 /*Partition penalties in 3d: */
397 if(strcmp(model->meshtype,"3d")==0){
398
399 /*Get penalties: */
400 ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
401
402 if(model->numpenalties){
403
404 model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
405 #ifdef _SERIAL_
406 for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
407 #else
408 for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
409
410 for(i=0;i<model->numpenalties;i++){
411 first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
412 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids
413 /*All grids that are being penalised belong to this node's internal grid partition.:*/
414 model->penaltypartitioning[i]=1;
415 }
416 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
417 model->penaltypartitioning[i]=0;
418 }
419 }
420 #endif
421 }
422
423 /*Free penalties: */
424 xfree((void**)&model->penalties);
425 }
426
427 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
428 We can therefore determine which grids are internal to this node's partition
429 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
430 that, go and create the grids*/
431
432 /*Create nodes from x,y,z, as well as the spc values on those grids: */
433
434 /*First fetch data: */
435 if (strcmp(model->meshtype,"3d")==0){
436 ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
437 ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
438 }
439 ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
440 ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
441 ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
442 ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
443 ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
444
445
446 /*Get number of dofs per node: */
447 DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
448
449 for (i=0;i<model->numberofnodes;i++){
450 #ifdef _PARALLEL_
451 /*keep only this partition's nodes:*/
452 if((my_grids[i]==1)){
453 #endif
454
455 node_id=i+1; //matlab indexing
456
457
458
459 #ifdef _PARALLEL_
460 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
461 node_partitionborder=1;
462 }
463 else{
464 node_partitionborder=0;
465 }
466 #else
467 node_partitionborder=0;
468 #endif
469
470
471 node_x[0]=model->x[i];
472 node_x[1]=model->y[i];
473 node_x[2]=model->z[i];
474
475
476 node_onbed=(int)model->gridonbed[i];
477 node_onsurface=(int)model->gridonsurface[i];
478 if (strcmp(model->meshtype,"3d")==0){
479 if (isnan(model->uppernodes[i])){
480 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.
481 }
482 else{
483 node_upper_node_id=(int)model->uppernodes[i];
484 }
485 }
486 else{
487 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
488 node_upper_node_id=node_id;
489 }
490
491 /*Create node using its constructor: */
492 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
493
494 /*set single point constraints.: */
495 if (strcmp(model->meshtype,"3d")==0){
496 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
497 if (!model->gridonbed[i]){
498 for(k=1;k<=node_numdofs;k++){
499 node->FreezeDof(k);
500 }
501 }
502 }
503 /*Add node to nodes dataset: */
504 nodes->AddObject(node);
505
506 #ifdef _PARALLEL_
507 } //if((my_grids[i]==1))
508 #endif
509 }
510
511 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
512 * datasets, it will not be redone: */
513 elements->Presort();
514 nodes->Presort();
515 materials->Presort();
516
517 /*Clean fetched data: */
518 xfree((void**)&model->deadgrids);
519 xfree((void**)&model->x);
520 xfree((void**)&model->y);
521 xfree((void**)&model->z);
522 xfree((void**)&model->gridonbed);
523 xfree((void**)&model->gridonsurface);
524 xfree((void**)&model->uppernodes);
525
526
527 /*Keep partitioning information into model*/
528 model->epart=epart;
529 model->my_grids=my_grids;
530 model->my_bordergrids=my_bordergrids;
531
532 /*Free ressources:*/
533 #ifdef _PARALLEL_
534 xfree((void**)&all_numgrids);
535 xfree((void**)&npart);
536 VecFree(&gridborder);
537 #endif
538
539 cleanup_and_return:
540
541 /*Assign output pointer: */
542 *pelements=elements;
543 *pnodes=nodes;
544 *pmaterials=materials;
545
546}
Note: See TracBrowser for help on using the repository browser.