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