![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3424
- Timestamp:
- 04/07/10 15:51:50 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp
r3417 r3424 2 2 * CreateElementsNodesAndMaterialsPrognostic.c: 3 3 */ 4 5 4 6 5 #include "../../DataSet/DataSet.h" … … 13 12 #include "../IoModel.h" 14 13 15 16 14 void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 17 18 15 19 16 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 20 17 21 22 18 int i,j,k,n; 23 extern int my_rank;24 extern int num_procs;25 19 26 20 /*DataSets: */ … … 30 24 DataSet* materials = NULL; 31 25 32 /*Objects: */33 Node* node = NULL;34 Vertex* vertex = NULL;35 Tria* tria = NULL;36 Penta* penta = NULL;37 Matice* matice = NULL;38 Matpar* matpar = NULL;39 ElementProperties* tria_properties=NULL;40 ElementProperties* penta_properties=NULL;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 elements51 double B_avg;52 53 /*tria constructor input: */54 int tria_id;55 int tria_matice_id;56 int tria_matpar_id;57 int tria_numpar_id;58 int tria_node_ids[3];59 double tria_h[3];60 double tria_s[3];61 double tria_b[3];62 int tria_shelf;63 bool tria_onwater;64 65 /*matice constructor input: */66 int matice_mid;67 double matice_B;68 double matice_n;69 70 /*penta constructor input: */71 72 int penta_id;73 int penta_matice_id;74 int penta_matpar_id;75 int penta_numpar_id;76 int penta_node_ids[6];77 double penta_h[6];78 double penta_s[6];79 double penta_b[6];80 int penta_shelf;81 int penta_onbed;82 int penta_onsurface;83 int penta_collapse;84 bool penta_onwater;85 86 /*matpar constructor input: */87 int matpar_mid;88 double matpar_rho_ice;89 double matpar_rho_water;90 double matpar_heatcapacity;91 double matpar_thermalconductivity;92 double matpar_latentheat;93 double matpar_beta;94 double matpar_meltingpoint;95 double matpar_mixed_layer_capacity;96 double matpar_thermal_exchange_velocity;97 double matpar_g;98 99 /* node constructor input: */100 int node_id;101 int vertex_id;102 int node_partitionborder=0;103 int node_onbed;104 int node_onsurface;105 int node_onshelf;106 int node_onsheet;107 int node_upper_node_id;108 int node_numdofs;109 110 111 #ifdef _PARALLEL_112 /*Metis partitioning: */113 int range;114 Vec gridborder=NULL;115 int my_numgrids;116 int* all_numgrids=NULL;117 int gridcount;118 int count;119 #endif120 int first_grid_index;121 122 26 /*First create the elements, nodes and material properties: */ 123 27 elements = new DataSet(ElementsEnum()); … … 126 30 materials = new DataSet(MaterialsEnum()); 127 31 128 /*Width of elements: */ 129 if(strcmp(iomodel->meshtype,"2d")==0){ 130 elements_width=3; //tria elements 131 } 132 else{ 133 elements_width=6; //penta elements 134 } 135 136 #ifdef _PARALLEL_ 137 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 138 if(strcmp(iomodel->meshtype,"2d")==0){ 139 /*load elements: */ 140 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 141 } 142 else{ 143 /*load elements2d: */ 144 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 145 } 146 147 148 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 149 150 /*Free elements and elements2d: */ 151 xfree((void**)&iomodel->elements); 152 xfree((void**)&iomodel->elements2d); 153 154 /*Used later on: */ 155 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); 156 #endif 157 158 159 160 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ 32 /*Partition elements and vertices and nodes: */ 33 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 161 34 162 35 /*2d mesh: */ … … 173 46 for (i=0;i<iomodel->numberofelements;i++){ 174 47 175 #ifdef _PARALLEL_ 176 /*!All elements have been partitioned above, only create elements for this CPU: */ 177 if(my_rank==epart[i]){ 178 #endif 179 180 181 /*ids: */ 182 tria_id=i+1; //matlab indexing. 183 tria_matice_id=-1; //no need for materials 184 tria_matpar_id=-1; //no need for materials 185 tria_numpar_id=1; 48 if(iomodel->my_elements[i]){ 186 49 187 /*vertices offsets: */ 188 tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0); 189 tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1); 190 tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2); 50 /*Create and add tria element to elements dataset: */ 51 elements->AddObject(new Tria(i,iomodel)); 191 52 192 /*thickness,surface and bed:*/ 193 tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. 194 tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 195 tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ; 196 197 tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 198 tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 199 tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 200 201 tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 202 tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 203 tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 204 205 /*element on iceshelf?:*/ 206 tria_shelf=(int)*(iomodel->elementoniceshelf+i); 207 tria_onwater=(bool)*(iomodel->elementonwater+i); 208 209 /*Create properties: */ 210 tria_properties=new ElementProperties(3,tria_h, tria_s, tria_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, tria_shelf, UNDEF, tria_onwater, UNDEF, UNDEF, UNDEF); 211 212 /*Create tria element using its constructor:*/ 213 tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties); 214 215 /*Delete properties: */ 216 delete tria_properties; 217 218 /*Add tria element to elements dataset: */ 219 elements->AddObject(tria); 220 221 #ifdef _PARALLEL_ 222 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 223 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 224 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 225 will hold which grids belong to this partition*/ 226 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; 227 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; 228 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; 229 #endif 230 231 #ifdef _PARALLEL_ 232 }//if(my_rank==epart[i]) 233 #endif 234 53 /*Create and add material property to materials dataset: */ 54 materials->AddObject(new Matice(i,iomodel,3)); 55 } 235 56 }//for (i=0;i<numberofelements;i++) 236 57 … … 258 79 259 80 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 81 if(iomodel->my_elements[i]){ 82 /*Create and add penta element to elements dataset: */ 83 elements->AddObject(new Penta(i,iomodel)); 264 84 265 266 /*name and id: */ 267 penta_id=i+1; //matlab indexing. 268 penta_matice_id=-1; 269 penta_matpar_id=-1; //no need for materials 270 penta_numpar_id=1; 271 272 /*vertices,thickness,surface,bed and drag: */ 273 for(j=0;j<6;j++){ 274 penta_node_ids[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)); 85 /*Create and add material property to materials dataset: */ 86 materials->AddObject(new Matice(i,iomodel,6)); 278 87 } 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 /*Create properties: */288 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);289 290 /*Create Penta using its constructor:*/291 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);292 293 /*Delete properties: */294 delete penta_properties;295 296 /*Add penta element to elements dataset: */297 elements->AddObject(penta);298 299 #ifdef _PARALLEL_300 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use301 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)302 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids303 will hold which grids belong to this partition*/304 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;305 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;306 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;307 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;308 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;309 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;310 #endif311 312 #ifdef _PARALLEL_313 }//if(my_rank==epart[i])314 #endif315 316 88 }//for (i=0;i<numberofelements;i++) 317 89 … … 328 100 } //if (strcmp(meshtype,"2d")==0) 329 101 330 #ifdef _PARALLEL_ 331 /*From the element partitioning, we can determine which grids are on the inside of this cpu's 332 *element partition, and which are on its border with other nodes:*/ 333 gridborder=NewVec(iomodel->numberofnodes); 334 335 for (i=0;i<iomodel->numberofnodes;i++){ 336 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); 337 } 338 VecAssemblyBegin(gridborder); 339 VecAssemblyEnd(gridborder); 340 341 #ifdef _ISSM_DEBUG_ 342 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); 343 #endif 344 345 VecToMPISerial(&my_bordergrids,gridborder); 346 347 #ifdef _ISSM_DEBUG_ 348 if(my_rank==0){ 349 for (i=0;i<iomodel->numberofnodes;i++){ 350 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); 351 } 352 } 353 #endif 354 #endif 355 356 /*Add one constant material property to materials: */ 357 matpar_mid=1; //put it at the end of the materials 358 matpar_g=iomodel->g; 359 matpar_rho_ice=iomodel->rho_ice; 360 matpar_rho_water=iomodel->rho_water; 361 matpar_thermalconductivity=iomodel->thermalconductivity; 362 matpar_heatcapacity=iomodel->heatcapacity; 363 matpar_latentheat=iomodel->latentheat; 364 matpar_beta=iomodel->beta; 365 matpar_meltingpoint=iomodel->meltingpoint; 366 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 367 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 368 369 /*Create matpar object using its constructor: */ 370 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, 371 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, 372 matpar_thermal_exchange_velocity,matpar_g); 373 374 /*Add to materials datset: */ 375 materials->AddObject(matpar); 376 377 /*Partition penalties in 3d: */ 378 if(strcmp(iomodel->meshtype,"3d")==0){ 379 380 /*Get penalties: */ 381 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); 382 383 if(iomodel->numpenalties){ 384 385 iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); 386 #ifdef _SERIAL_ 387 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1; 388 #else 389 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1; 390 391 for(i=0;i<iomodel->numpenalties;i++){ 392 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1); 393 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids 394 /*All grids that are being penalised belong to this node's internal grid partition.:*/ 395 iomodel->penaltypartitioning[i]=1; 396 } 397 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border 398 iomodel->penaltypartitioning[i]=0; 399 } 400 } 401 #endif 402 } 403 404 /*Free penalties: */ 405 xfree((void**)&iomodel->penalties); 406 } 407 408 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 409 We can therefore determine which grids are internal to this node's partition 410 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 411 that, go and create the grids*/ 412 413 /*Create nodes from x,y,z, as well as the spc values on those grids: */ 102 /*Add new constrant material property to materials, at the end: */ 103 materials->AddObject(new Matpar(iomodel)); 414 104 415 105 /*First fetch data: */ … … 428 118 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 429 119 120 for (i=0;i<iomodel->numberofvertices;i++){ 430 121 431 /*Get number of dofs per node: */432 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);122 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 123 if(iomodel->my_vertices[i]){ 433 124 434 for (i=0;i<iomodel->numberofnodes;i++){ 435 #ifdef _PARALLEL_ 436 /*keep only this partition's nodes:*/ 437 if((my_grids[i]==1)){ 438 #endif 125 /*Add vertex to vertices dataset: */ 126 vertices->AddObject(new Vertex(i,iomodel)); 439 127 440 /*create vertex: */ 441 vertex_id=i+1; 442 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); 443 vertices->AddObject(vertex); 128 /*Add node to nodes dataset: */ 129 nodes->AddObject(new Node(i,iomodel)); 444 130 445 446 node_id=i+1; //matlab indexing447 448 #ifdef _PARALLEL_449 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border450 node_partitionborder=1;451 131 } 452 else{453 node_partitionborder=0;454 }455 #else456 node_partitionborder=0;457 #endif458 459 node_properties.SetProperties(460 (int)iomodel->gridonbed[i],461 (int)iomodel->gridonsurface[i],462 (int)iomodel->gridoniceshelf[i],463 (int)iomodel->gridonicesheet[i]);464 465 466 if (strcmp(iomodel->meshtype,"3d")==0){467 if (isnan(iomodel->uppernodes[i])){468 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.469 }470 else{471 node_upper_node_id=(int)iomodel->uppernodes[i];472 }473 }474 else{475 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/476 node_upper_node_id=node_id;477 }478 479 /*Create node using its constructor: */480 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);481 482 /*set single point constraints.: */483 if (strcmp(iomodel->meshtype,"3d")==0){484 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */485 if (!iomodel->gridonbed[i]){486 for(k=1;k<=node_numdofs;k++){487 node->FreezeDof(k);488 }489 }490 }491 /*Add node to nodes dataset: */492 nodes->AddObject(node);493 494 #ifdef _PARALLEL_495 } //if((my_grids[i]==1))496 #endif497 132 } 498 499 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these500 * datasets, it will not be redone: */501 elements->Presort();502 nodes->Presort();503 vertices->Presort();504 materials->Presort();505 133 506 134 /*Clean fetched data: */ … … 517 145 xfree((void**)&iomodel->gridoniceshelf); 518 146 519 520 /*Keep partitioning information into iomodel*/ 521 iomodel->epart=epart; 522 iomodel->my_grids=my_grids; 523 iomodel->my_bordergrids=my_bordergrids; 524 525 /*Free ressources:*/ 526 #ifdef _PARALLEL_ 527 xfree((void**)&all_numgrids); 528 xfree((void**)&npart); 529 VecFree(&gridborder); 530 #endif 147 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 148 * datasets, it will not be redone: */ 149 elements->Presort(); 150 nodes->Presort(); 151 vertices->Presort(); 152 materials->Presort(); 531 153 532 154 cleanup_and_return:
Note:
See TracChangeset
for help on using the changeset viewer.
