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

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

need a matpar in prognostic for transient solution
File size: 18.6 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_onshelf;
126 int node_onsheet;
127 int node_upper_node_id;
128 int node_numdofs;
129
130
131 #ifdef _PARALLEL_
132 /*Metis partitioning: */
133 int range;
134 Vec gridborder=NULL;
135 int my_numgrids;
136 int* all_numgrids=NULL;
137 int gridcount;
138 int count;
139 #endif
140 int first_grid_index;
141
142 /*Rifts:*/
143 int* riftsnumpenaltypairs;
144 double** riftspenaltypairs;
145 int* riftsfill;
146 double* riftsfriction;
147 double* riftpenaltypairs=NULL;
148 int el1,el2;
149
150 /*First create the elements, nodes and material properties: */
151 elements = new DataSet(ElementsEnum());
152 nodes = new DataSet(NodesEnum());
153 materials = new DataSet(MaterialsEnum());
154
155
156 /*Get analysis_type: */
157 analysis_type=AnalysisTypeAsEnum(model->analysis_type);
158 sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
159
160
161 /*Width of elements: */
162 if(strcmp(model->meshtype,"2d")==0){
163 elements_width=3; //tria elements
164 }
165 else{
166 elements_width=6; //penta elements
167 }
168
169 #ifdef _PARALLEL_
170 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
171 if(strcmp(model->meshtype,"2d")==0){
172 /*load elements: */
173 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
174 }
175 else{
176 /*load elements2d: */
177 ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
178 }
179
180
181 MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
182
183 /*Free elements and elements2d: */
184 xfree((void**)&model->elements);
185 xfree((void**)&model->elements2d);
186
187
188 /*Deal with rifts, they have to be included into one partition only, not several: */
189 FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
190
191 for(i=0;i<model->numrifts;i++){
192 riftpenaltypairs=model->riftspenaltypairs[i];
193 for(j=0;j<model->riftsnumpenaltypairs[i];j++){
194 el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
195 el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
196 epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
197 }
198 }
199 /*Free rifts: */
200 xfree((void**)&riftsnumpenaltypairs);
201 for(i=0;i<model->numrifts;i++){
202 double* temp=riftspenaltypairs[i];
203 xfree((void**)&temp);
204 }
205 xfree((void**)&riftspenaltypairs);
206 xfree((void**)&riftsfill);
207 xfree((void**)&riftsfriction);
208
209 /*Used later on: */
210 my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
211 #endif
212
213
214
215 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
216
217 /*2d mesh: */
218 if (strcmp(model->meshtype,"2d")==0){
219
220 /*Fetch data needed: */
221 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
222 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
223 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
224 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
225 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
226
227 for (i=0;i<model->numberofelements;i++){
228
229 #ifdef _PARALLEL_
230 /*!All elements have been partitioned above, only create elements for this CPU: */
231 if(my_rank==epart[i]){
232 #endif
233
234
235 /*ids: */
236 tria_id=i+1; //matlab indexing.
237 tria_mid=-1; //no need for materials
238 tria_mparid=-1; //no need for materials
239
240 /*vertices offsets: */
241 tria_g[0]=(int)*(model->elements+elements_width*i+0);
242 tria_g[1]=(int)*(model->elements+elements_width*i+1);
243 tria_g[2]=(int)*(model->elements+elements_width*i+2);
244
245 /*thickness,surface and bed:*/
246 tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
247 tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1));
248 tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ;
249
250 tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1));
251 tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1));
252 tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1));
253
254 tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1));
255 tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1));
256 tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1));
257
258 /*element on iceshelf?:*/
259 tria_shelf=(int)*(model->elementoniceshelf+i);
260 tria_artdiff=model->artificial_diffusivity;
261
262 /*Create tria element using its constructor:*/
263 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);
264
265 /*Add tria element to elements dataset: */
266 elements->AddObject(tria);
267
268 #ifdef _PARALLEL_
269 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
270 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
271 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
272 will hold which grids belong to this partition*/
273 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
274 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
275 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
276 #endif
277
278 #ifdef _PARALLEL_
279 }//if(my_rank==epart[i])
280 #endif
281
282 }//for (i=0;i<numberofelements;i++)
283
284
285 /*Free data : */
286 xfree((void**)&model->elements);
287 xfree((void**)&model->thickness);
288 xfree((void**)&model->surface);
289 xfree((void**)&model->bed);
290 xfree((void**)&model->elementoniceshelf);
291
292 }
293 else{ // if (strcmp(meshtype,"2d")==0)
294
295 /*Fetch data needed: */
296 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
297 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
298 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
299 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
300 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
301 ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
302 ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
303
304 for (i=0;i<model->numberofelements;i++){
305 #ifdef _PARALLEL_
306 /*We are using our element partition to decide which elements will be created on this node: */
307 if(my_rank==epart[i]){
308 #endif
309
310
311 /*name and id: */
312 penta_id=i+1; //matlab indexing.
313 penta_mid=-1;
314 penta_mparid=-1; //no need for materials
315
316 /*vertices,thickness,surface,bed and drag: */
317 for(j=0;j<6;j++){
318 penta_g[j]=(int)*(model->elements+elements_width*i+j);
319 penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1));
320 penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1));
321 penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1));
322 }
323
324 /*diverse: */
325 penta_shelf=(int)*(model->elementoniceshelf+i);
326 penta_onbed=(int)*(model->elementonbed+i);
327 penta_onsurface=(int)*(model->elementonsurface+i);
328 penta_collapse=1;
329 penta_artdiff=model->artificial_diffusivity;
330
331
332 /*Create Penta using its constructor:*/
333 penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
334 penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
335 penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
336 penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning);
337
338 /*Add penta element to elements dataset: */
339 elements->AddObject(penta);
340
341 #ifdef _PARALLEL_
342 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
343 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
344 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
345 will hold which grids belong to this partition*/
346 my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
347 my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
348 my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
349 my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
350 my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
351 my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
352 #endif
353
354 #ifdef _PARALLEL_
355 }//if(my_rank==epart[i])
356 #endif
357
358 }//for (i=0;i<numberofelements;i++)
359
360 /*Free data: */
361 xfree((void**)&model->elements);
362 xfree((void**)&model->thickness);
363 xfree((void**)&model->surface);
364 xfree((void**)&model->bed);
365 xfree((void**)&model->elementoniceshelf);
366 xfree((void**)&model->elementonbed);
367 xfree((void**)&model->elementonsurface);
368
369 } //if (strcmp(meshtype,"2d")==0)
370
371 #ifdef _PARALLEL_
372 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
373 *element partition, and which are on its border with other nodes:*/
374 gridborder=NewVec(model->numberofnodes);
375
376 for (i=0;i<model->numberofnodes;i++){
377 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
378 }
379 VecAssemblyBegin(gridborder);
380 VecAssemblyEnd(gridborder);
381
382 #ifdef _DEBUG_
383 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
384 #endif
385
386 VecToMPISerial(&my_bordergrids,gridborder);
387
388 #ifdef _DEBUG_
389 if(my_rank==0){
390 for (i=0;i<model->numberofnodes;i++){
391 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
392 }
393 }
394 #endif
395 #endif
396
397 /*Add one constant material property to materials: */
398 matpar_mid=1; //put it at the end of the materials
399 matpar_g=model->g;
400 matpar_rho_ice=model->rho_ice;
401 matpar_rho_water=model->rho_water;
402 matpar_thermalconductivity=model->thermalconductivity;
403 matpar_heatcapacity=model->heatcapacity;
404 matpar_latentheat=model->latentheat;
405 matpar_beta=model->beta;
406 matpar_meltingpoint=model->meltingpoint;
407 matpar_mixed_layer_capacity=model->mixed_layer_capacity;
408 matpar_thermal_exchange_velocity=model->thermal_exchange_velocity;
409
410 /*Create matpar object using its constructor: */
411 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
412 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
413 matpar_thermal_exchange_velocity,matpar_g);
414
415 /*Add to materials datset: */
416 materials->AddObject(matpar);
417
418 /*Partition penalties in 3d: */
419 if(strcmp(model->meshtype,"3d")==0){
420
421 /*Get penalties: */
422 ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
423
424 if(model->numpenalties){
425
426 model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
427 #ifdef _SERIAL_
428 for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
429 #else
430 for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
431
432 for(i=0;i<model->numpenalties;i++){
433 first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
434 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids
435 /*All grids that are being penalised belong to this node's internal grid partition.:*/
436 model->penaltypartitioning[i]=1;
437 }
438 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
439 model->penaltypartitioning[i]=0;
440 }
441 }
442 #endif
443 }
444
445 /*Free penalties: */
446 xfree((void**)&model->penalties);
447 }
448
449 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
450 We can therefore determine which grids are internal to this node's partition
451 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
452 that, go and create the grids*/
453
454 /*Create nodes from x,y,z, as well as the spc values on those grids: */
455
456 /*First fetch data: */
457 if (strcmp(model->meshtype,"3d")==0){
458 ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
459 ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
460 }
461 ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
462 ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
463 ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
464 ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
465 ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
466 ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
467 ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
468
469
470 /*Get number of dofs per node: */
471 DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
472
473 for (i=0;i<model->numberofnodes;i++){
474 #ifdef _PARALLEL_
475 /*keep only this partition's nodes:*/
476 if((my_grids[i]==1)){
477 #endif
478
479 node_id=i+1; //matlab indexing
480
481
482
483 #ifdef _PARALLEL_
484 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
485 node_partitionborder=1;
486 }
487 else{
488 node_partitionborder=0;
489 }
490 #else
491 node_partitionborder=0;
492 #endif
493
494
495 node_x[0]=model->x[i];
496 node_x[1]=model->y[i];
497 node_x[2]=model->z[i];
498
499
500 node_onbed=(int)model->gridonbed[i];
501 node_onsurface=(int)model->gridonsurface[i];
502 node_onshelf=(int)model->gridoniceshelf[i];
503 node_onsheet=(int)model->gridonicesheet[i];
504
505 if (strcmp(model->meshtype,"3d")==0){
506 if (isnan(model->uppernodes[i])){
507 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.
508 }
509 else{
510 node_upper_node_id=(int)model->uppernodes[i];
511 }
512 }
513 else{
514 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
515 node_upper_node_id=node_id;
516 }
517
518 /*Create node using its constructor: */
519 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
520
521 /*set single point constraints.: */
522 if (strcmp(model->meshtype,"3d")==0){
523 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
524 if (!model->gridonbed[i]){
525 for(k=1;k<=node_numdofs;k++){
526 node->FreezeDof(k);
527 }
528 }
529 }
530 /*Add node to nodes dataset: */
531 nodes->AddObject(node);
532
533 #ifdef _PARALLEL_
534 } //if((my_grids[i]==1))
535 #endif
536 }
537
538 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
539 * datasets, it will not be redone: */
540 elements->Presort();
541 nodes->Presort();
542 materials->Presort();
543
544 /*Clean fetched data: */
545 xfree((void**)&model->deadgrids);
546 xfree((void**)&model->x);
547 xfree((void**)&model->y);
548 xfree((void**)&model->z);
549 xfree((void**)&model->gridonbed);
550 xfree((void**)&model->gridonsurface);
551 xfree((void**)&model->uppernodes);
552 xfree((void**)&model->gridonicesheet);
553 xfree((void**)&model->gridoniceshelf);
554
555
556 /*Keep partitioning information into model*/
557 model->epart=epart;
558 model->my_grids=my_grids;
559 model->my_bordergrids=my_bordergrids;
560
561 /*Free ressources:*/
562 #ifdef _PARALLEL_
563 xfree((void**)&all_numgrids);
564 xfree((void**)&npart);
565 VecFree(&gridborder);
566 #endif
567
568 cleanup_and_return:
569
570 /*Assign output pointer: */
571 *pelements=elements;
572 *pnodes=nodes;
573 *pmaterials=materials;
574
575}
Note: See TracBrowser for help on using the repository browser.