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

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

Big commit: created Numpar, new object to hold solution parameters necessary
in elements. This lead to creating FetchParams and WriteParams, which now writes
a DataSet* parameters to a matlab workspace structure and vice versa. We now always have
a DataSet* parametes inside the x code. Introduced also a new configuration phase for the paramters
dataset. Also, rewrote the io/ using overloaded functions IoModelFetchData, FetchData and WriteData.
Much cleaner and less error prone, as arguments are consistently checked.

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