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

Last change on this file since 1025 was 1025, checked in by Mathieu Morlighem, 16 years ago

Use sigma to update node position
File size: 18.9 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 double node_sigma;
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->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
466 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
467 ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
468 ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
469 ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
470 ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
471
472
473 /*Get number of dofs per node: */
474 DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
475
476 for (i=0;i<model->numberofnodes;i++){
477 #ifdef _PARALLEL_
478 /*keep only this partition's nodes:*/
479 if((my_grids[i]==1)){
480 #endif
481
482 node_id=i+1; //matlab indexing
483
484
485
486 #ifdef _PARALLEL_
487 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
488 node_partitionborder=1;
489 }
490 else{
491 node_partitionborder=0;
492 }
493 #else
494 node_partitionborder=0;
495 #endif
496
497 node_x[0]=model->x[i];
498 node_x[1]=model->y[i];
499 node_x[2]=model->z[i];
500 node_sigma=(model->z[i]-model->bed[i])/(model->thickness[i]);
501
502 node_onbed=(int)model->gridonbed[i];
503 node_onsurface=(int)model->gridonsurface[i];
504 node_onshelf=(int)model->gridoniceshelf[i];
505 node_onsheet=(int)model->gridonicesheet[i];
506
507 if (strcmp(model->meshtype,"3d")==0){
508 if (isnan(model->uppernodes[i])){
509 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.
510 }
511 else{
512 node_upper_node_id=(int)model->uppernodes[i];
513 }
514 }
515 else{
516 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
517 node_upper_node_id=node_id;
518 }
519
520 /*Create node using its constructor: */
521 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);
522
523 /*set single point constraints.: */
524 if (strcmp(model->meshtype,"3d")==0){
525 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
526 if (!model->gridonbed[i]){
527 for(k=1;k<=node_numdofs;k++){
528 node->FreezeDof(k);
529 }
530 }
531 }
532 /*Add node to nodes dataset: */
533 nodes->AddObject(node);
534
535 #ifdef _PARALLEL_
536 } //if((my_grids[i]==1))
537 #endif
538 }
539
540 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
541 * datasets, it will not be redone: */
542 elements->Presort();
543 nodes->Presort();
544 materials->Presort();
545
546 /*Clean fetched data: */
547 xfree((void**)&model->deadgrids);
548 xfree((void**)&model->x);
549 xfree((void**)&model->y);
550 xfree((void**)&model->z);
551 xfree((void**)&model->thickness);
552 xfree((void**)&model->bed);
553 xfree((void**)&model->gridonbed);
554 xfree((void**)&model->gridonsurface);
555 xfree((void**)&model->uppernodes);
556 xfree((void**)&model->gridonicesheet);
557 xfree((void**)&model->gridoniceshelf);
558
559
560 /*Keep partitioning information into model*/
561 model->epart=epart;
562 model->my_grids=my_grids;
563 model->my_bordergrids=my_bordergrids;
564
565 /*Free ressources:*/
566 #ifdef _PARALLEL_
567 xfree((void**)&all_numgrids);
568 xfree((void**)&npart);
569 VecFree(&gridborder);
570 #endif
571
572 cleanup_and_return:
573
574 /*Assign output pointer: */
575 *pelements=elements;
576 *pnodes=nodes;
577 *pmaterials=materials;
578
579}
Note: See TracBrowser for help on using the repository browser.