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

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

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