![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3434 for issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
- Timestamp:
- 04/07/10 16:24:05 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
r3417 r3434 16 16 17 17 18 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 19 20 18 /*Intermediary*/ 21 19 int i,j,k,n; 22 extern int my_rank;23 extern int num_procs;24 20 25 21 /*DataSets: */ … … 28 24 DataSet* vertices = NULL; 29 25 DataSet* materials = NULL; 30 31 /*Objects: */32 Node* node = NULL;33 NodeProperties node_properties;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 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 26 122 27 /*First create the elements, nodes and material properties: */ … … 126 31 materials = new DataSet(MaterialsEnum()); 127 32 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: */ 33 /*Partition elements and vertices and nodes: */ 34 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 161 35 162 36 /*2d mesh: */ … … 173 47 for (i=0;i<iomodel->numberofelements;i++){ 174 48 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; 49 if(iomodel->my_elements[i]){ 186 50 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); 51 /*Create and add tria element to elements dataset: */ 52 elements->AddObject(new Tria(i,iomodel)); 191 53 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 54 /*Create and add material property to materials dataset: */ 55 materials->AddObject(new Matice(i,iomodel,3)); 56 } 235 57 }//for (i=0;i<numberofelements;i++) 236 58 237 238 59 /*Free data : */ 239 60 xfree((void**)&iomodel->elements); … … 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 288 /*Create properties: */289 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);290 291 /*Create Penta using its constructor:*/292 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);293 294 /*Delete properties: */295 delete penta_properties;296 297 /*Add penta element to elements dataset: */298 elements->AddObject(penta);299 300 #ifdef _PARALLEL_301 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use302 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)303 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids304 will hold which grids belong to this partition*/305 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;306 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;307 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;308 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;309 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;310 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;311 #endif312 313 #ifdef _PARALLEL_314 }//if(my_rank==epart[i])315 #endif316 317 88 }//for (i=0;i<numberofelements;i++) 318 89 … … 329 100 } //if (strcmp(meshtype,"2d")==0) 330 101 331 #ifdef _PARALLEL_ 332 /*From the element partitioning, we can determine which grids are on the inside of this cpu's 333 *element partition, and which are on its border with other nodes:*/ 334 gridborder=NewVec(iomodel->numberofnodes); 102 /*Add new constrant material property to materials, at the end: */ 103 materials->AddObject(new Matpar(iomodel)); 335 104 336 for (i=0;i<iomodel->numberofnodes;i++){337 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);338 }339 VecAssemblyBegin(gridborder);340 VecAssemblyEnd(gridborder);341 342 #ifdef _ISSM_DEBUG_343 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);344 #endif345 346 VecToMPISerial(&my_bordergrids,gridborder);347 348 #ifdef _ISSM_DEBUG_349 if(my_rank==0){350 for (i=0;i<iomodel->numberofnodes;i++){351 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);352 }353 }354 #endif355 #endif356 357 /*Add one constant material property to materials: */358 matpar_mid=1; //put it at the end of the materials359 matpar_g=iomodel->g;360 matpar_rho_ice=iomodel->rho_ice;361 matpar_rho_water=iomodel->rho_water;362 matpar_thermalconductivity=iomodel->thermalconductivity;363 matpar_heatcapacity=iomodel->heatcapacity;364 matpar_latentheat=iomodel->latentheat;365 matpar_beta=iomodel->beta;366 matpar_meltingpoint=iomodel->meltingpoint;367 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;368 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;369 370 /*Create matpar object using its constructor: */371 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,372 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,373 matpar_thermal_exchange_velocity,matpar_g);374 375 /*Add to materials datset: */376 materials->AddObject(matpar);377 378 /*Partition penalties in 3d: */379 if(strcmp(iomodel->meshtype,"3d")==0){380 381 /*Get penalties: */382 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");383 384 if(iomodel->numpenalties){385 386 iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));387 #ifdef _SERIAL_388 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;389 #else390 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;391 392 for(i=0;i<iomodel->numpenalties;i++){393 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);394 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids395 /*All grids that are being penalised belong to this node's internal grid partition.:*/396 iomodel->penaltypartitioning[i]=1;397 }398 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border399 iomodel->penaltypartitioning[i]=0;400 }401 }402 #endif403 }404 405 /*Free penalties: */406 xfree((void**)&iomodel->penalties);407 }408 409 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.410 We can therefore determine which grids are internal to this node's partition411 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing412 that, go and create the grids*/413 414 /*Create nodes from x,y,z, as well as the spc values on those grids: */415 416 105 /*First fetch data: */ 417 106 if (strcmp(iomodel->meshtype,"3d")==0){ … … 429 118 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 430 119 120 for (i=0;i<iomodel->numberofvertices;i++){ 431 121 432 /*Get number of dofs per node: */433 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]){ 434 124 435 for (i=0;i<iomodel->numberofnodes;i++){ 436 #ifdef _PARALLEL_ 437 /*keep only this partition's nodes:*/ 438 if((my_grids[i]==1)){ 439 #endif 125 /*Add vertex to vertices dataset: */ 126 vertices->AddObject(new Vertex(i,iomodel)); 440 127 441 #ifdef _PARALLEL_442 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border443 partitionborder=1; 128 /*Add node to nodes dataset: */ 129 nodes->AddObject(new Node(i,iomodel)); 130 444 131 } 445 else{446 partitionborder=0;447 }448 #else449 partitionborder=0;450 #endif451 452 /*create vertex: */453 vertex_id=i+1;454 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder);455 vertices->AddObject(vertex);456 457 /*create node: */458 node_id=i+1; //matlab indexing459 460 461 node_properties.SetProperties(462 (int)iomodel->gridonbed[i],463 (int)iomodel->gridonsurface[i],464 (int)iomodel->gridoniceshelf[i],465 (int)iomodel->gridonicesheet[i]);466 467 if (strcmp(iomodel->meshtype,"3d")==0){468 if (isnan(iomodel->uppernodes[i])){469 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.470 }471 else{472 node_upper_node_id=(int)iomodel->uppernodes[i];473 }474 }475 else{476 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/477 node_upper_node_id=node_id;478 }479 480 /*Create node using its constructor: */481 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);482 483 484 /*set single point constraints.: */485 if (strcmp(iomodel->meshtype,"3d")==0){486 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */487 if (!iomodel->gridonbed[i]){488 for(k=1;k<=node_numdofs;k++){489 node->FreezeDof(k);490 }491 }492 }493 /*Add node to nodes dataset: */494 nodes->AddObject(node);495 496 #ifdef _PARALLEL_497 } //if((my_grids[i]==1))498 #endif499 132 } 500 501 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these502 * datasets, it will not be redone: */503 elements->Presort();504 vertices->Presort();505 nodes->Presort();506 materials->Presort();507 133 508 134 /*Clean fetched data: */ … … 518 144 xfree((void**)&iomodel->gridonicesheet); 519 145 xfree((void**)&iomodel->gridoniceshelf); 520 521 146 522 /*Keep partitioning information into iomodel*/ 523 iomodel->epart=epart; 524 iomodel->my_grids=my_grids; 525 iomodel->my_bordergrids=my_bordergrids; 526 527 /*Free ressources:*/ 528 #ifdef _PARALLEL_ 529 xfree((void**)&all_numgrids); 530 xfree((void**)&npart); 531 VecFree(&gridborder); 532 #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 vertices->Presort(); 151 nodes->Presort(); 152 materials->Presort(); 533 153 534 154 cleanup_and_return:
Note:
See TracChangeset
for help on using the changeset viewer.
