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

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

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