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

Last change on this file since 1648 was 1648, checked in by Mathieu Morlighem, 15 years ago

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