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

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

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