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

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

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