source: issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp@ 3354

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

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