![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3425
- Timestamp:
- 04/07/10 15:58:33 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
r3417 r3425 2 2 * CreateElementsNodesAndMaterialsSlopeCompute.c: 3 3 */ 4 5 4 6 5 #include "../../DataSet/DataSet.h" … … 15 14 void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 16 15 17 18 16 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 19 17 20 21 18 int i,j,k,n; 22 extern int my_rank;23 extern int num_procs;24 25 19 /*DataSets: */ 26 20 DataSet* elements = NULL; … … 29 23 DataSet* materials = NULL; 30 24 31 /*Objects: */32 Node* node = NULL;33 Vertex* vertex = NULL;34 Tria* tria = NULL;35 Penta* penta = NULL;36 ElementProperties* tria_properties=NULL;37 ElementProperties* penta_properties=NULL;38 39 /*output: */40 int* epart=NULL; //element partitioning.41 int* npart=NULL; //node partitioning.42 int* my_grids=NULL;43 double* my_bordergrids=NULL;44 45 /*intermediary: */46 int elements_width; //size of elements47 48 /*tria constructor input: */49 int tria_id;50 int tria_matice_id;51 int tria_matpar_id;52 int tria_numpar_id;53 int tria_node_ids[3];54 double tria_s[3];55 double tria_b[3];56 bool tria_onwater;57 58 /*penta constructor input: */59 60 int penta_id;61 int penta_matice_id;62 int penta_matpar_id;63 int penta_numpar_id;64 int penta_node_ids[6];65 double penta_s[6];66 double penta_b[6];67 int penta_onbed;68 bool penta_onwater;69 70 /* node constructor input: */71 int node_id;72 int vertex_id;73 int node_partitionborder=0;74 int node_onbed;75 int node_onsurface;76 int node_onshelf;77 int node_onsheet;78 int node_upper_node_id;79 int node_numdofs;80 81 82 #ifdef _PARALLEL_83 /*Metis partitioning: */84 int range;85 Vec gridborder=NULL;86 int my_numgrids;87 int* all_numgrids=NULL;88 int gridcount;89 int count;90 #endif91 int first_grid_index;92 93 25 /*First create the elements, nodes and material properties: */ 94 26 elements = new DataSet(ElementsEnum()); … … 97 29 materials = new DataSet(MaterialsEnum()); 98 30 99 /*Now, is the flag isstokes on? otherwise, do nothing: */ 100 if (!iomodel->isstokes & !iomodel->ishutter)goto cleanup_and_return; 101 102 /*Width of elements: */ 103 if(strcmp(iomodel->meshtype,"2d")==0){ 104 elements_width=3; //tria elements 105 } 106 else{ 107 elements_width=6; //penta elements 108 } 109 110 #ifdef _PARALLEL_ 111 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 112 if(strcmp(iomodel->meshtype,"2d")==0){ 113 /*load elements: */ 114 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 115 } 116 else{ 117 /*load elements2d: */ 118 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 119 } 120 121 122 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 123 124 /*Free elements and elements2d: */ 125 xfree((void**)&iomodel->elements); 126 xfree((void**)&iomodel->elements2d); 127 128 129 /*Used later on: */ 130 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); 131 #endif 132 133 134 135 /*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); 136 33 137 34 /*2d mesh: */ … … 146 43 for (i=0;i<iomodel->numberofelements;i++){ 147 44 148 #ifdef _PARALLEL_ 149 /*!All elements have been partitioned above, only create elements for this CPU: */ 150 if(my_rank==epart[i]){ 151 #endif 152 153 154 /*ids: */ 155 tria_id=i+1; //matlab indexing. 156 tria_matice_id=-1; //no materials 157 tria_matpar_id=-1; //no materials 158 tria_numpar_id=1; //no materials 45 if(iomodel->my_elements[i]){ 46 47 /*Create and add tria element to elements dataset: */ 48 elements->AddObject(new Tria(i,iomodel)); 159 49 160 /*vertices offsets: */ 161 tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0); 162 tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1); 163 tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2); 164 165 /*surface and bed:*/ 166 tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 167 tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 168 tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 169 170 tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 171 tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 172 tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 173 174 /*element on water? : */ 175 tria_onwater=(bool)*(iomodel->elementonwater+i); 176 177 /*Create properties: */ 178 tria_properties=new ElementProperties(3,NULL, tria_s, tria_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, UNDEF, tria_onwater, UNDEF, UNDEF, UNDEF); 179 180 /*Create tria element using its constructor:*/ 181 tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties); 182 183 /*Delete properties: */ 184 delete tria_properties; 185 186 /*Add tria element to elements dataset: */ 187 elements->AddObject(tria); 188 189 #ifdef _PARALLEL_ 190 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 191 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 192 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 193 will hold which grids belong to this partition*/ 194 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; 195 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; 196 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; 197 #endif 198 199 #ifdef _PARALLEL_ 200 }//if(my_rank==epart[i]) 201 #endif 50 /*Create and add material property to materials dataset: */ 51 materials->AddObject(new Matice(i,iomodel,3)); 52 } 202 53 203 54 }//for (i=0;i<numberofelements;i++) … … 221 72 222 73 for (i=0;i<iomodel->numberofelements;i++){ 223 #ifdef _PARALLEL_ 224 /*We are using our element partition to decide which elements will be created on this node: */ 225 if(my_rank==epart[i]){ 226 #endif 74 if(iomodel->my_elements[i]){ 75 /*Create and add penta element to elements dataset: */ 76 elements->AddObject(new Penta(i,iomodel)); 227 77 228 229 /*name and id: */ 230 penta_id=i+1; //matlab indexing. 231 penta_matice_id=-1; //no materials 232 penta_matpar_id=-1; //no materials 233 penta_numpar_id=1; //no materials 234 78 /*Create and add material property to materials dataset: */ 79 materials->AddObject(new Matice(i,iomodel,6)); 235 80 236 /*vertices,thickness,surface,bed and drag: */ 237 for(j=0;j<6;j++){ 238 penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); 239 penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 240 penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 241 } 242 243 /*diverse: */ 244 penta_onbed=(int)*(iomodel->elementonbed+i); 245 penta_onwater=(bool)*(iomodel->elementonwater+i); 246 247 /*Create properties: */ 248 penta_properties=new ElementProperties(6,NULL, penta_s, penta_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF,UNDEF,penta_onbed, penta_onwater, UNDEF, UNDEF, UNDEF); 249 250 /*Create Penta using its constructor:*/ 251 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties); 252 253 /*Delete properties: */ 254 delete penta_properties; 255 256 /*Add penta element to elements dataset: */ 257 elements->AddObject(penta); 258 259 #ifdef _PARALLEL_ 260 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 261 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 262 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 263 will hold which grids belong to this partition*/ 264 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; 265 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; 266 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; 267 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1; 268 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1; 269 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1; 270 #endif 271 272 #ifdef _PARALLEL_ 273 }//if(my_rank==epart[i]) 274 #endif 275 81 }//if(my_elements[i]) 276 82 }//for (i=0;i<numberofelements;i++) 277 83 … … 285 91 } //if (strcmp(meshtype,"2d")==0) 286 92 287 #ifdef _PARALLEL_ 288 /*From the element partitioning, we can determine which grids are on the inside of this cpu's 289 *element partition, and which are on its border with other nodes:*/ 290 gridborder=NewVec(iomodel->numberofnodes); 291 292 for (i=0;i<iomodel->numberofnodes;i++){ 293 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); 294 } 295 VecAssemblyBegin(gridborder); 296 VecAssemblyEnd(gridborder); 297 298 #ifdef _ISSM_DEBUG_ 299 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); 300 #endif 301 302 VecToMPISerial(&my_bordergrids,gridborder); 303 304 #ifdef _ISSM_DEBUG_ 305 if(my_rank==0){ 306 for (i=0;i<iomodel->numberofnodes;i++){ 307 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); 308 } 309 } 310 #endif 311 #endif 312 313 /*Partition penalties in 3d: */ 314 if(strcmp(iomodel->meshtype,"3d")==0){ 93 /*Add new constrant material property tgo materials, at the end: */ 94 materials->AddObject(new Matpar(iomodel)); 315 95 316 /*Get penalties: */317 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");318 319 if(iomodel->numpenalties){320 321 iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));322 #ifdef _SERIAL_323 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;324 #else325 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;326 327 for(i=0;i<iomodel->numpenalties;i++){328 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);329 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids330 /*All grids that are being penalised belong to this node's internal grid partition.:*/331 iomodel->penaltypartitioning[i]=1;332 }333 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border334 iomodel->penaltypartitioning[i]=0;335 }336 }337 #endif338 }339 340 /*Free penalties: */341 xfree((void**)&iomodel->penalties);342 }343 344 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.345 We can therefore determine which grids are internal to this node's partition346 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing347 that, go and create the grids*/348 349 /*Create nodes from x,y,z, as well as the spc values on those grids: */350 351 96 /*First fetch data: */ 352 97 if (strcmp(iomodel->meshtype,"3d")==0){ … … 364 109 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 365 110 111 for (i=0;i<iomodel->numberofvertices;i++){ 366 112 367 /*Get number of dofs per node: */ 368 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); 113 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 114 if(iomodel->my_vertices[i]){ 115 116 /*Add vertex to vertices dataset: */ 117 vertices->AddObject(new Vertex(i,iomodel)); 369 118 370 for (i=0;i<iomodel->numberofnodes;i++){ 371 #ifdef _PARALLEL_ 372 /*keep only this partition's nodes:*/ 373 if((my_grids[i]==1)){ 374 #endif 119 /*Add node to nodes dataset: */ 120 nodes->AddObject(new Node(i,iomodel)); 375 121 376 /*create vertex: */377 vertex_id=i+1;378 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));379 vertices->AddObject(vertex);380 381 382 node_id=i+1; //matlab indexing383 384 #ifdef _PARALLEL_385 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border386 node_partitionborder=1;387 122 } 388 else{389 node_partitionborder=0;390 }391 #else392 node_partitionborder=0;393 #endif394 395 node_properties.SetProperties(396 (int)iomodel->gridonbed[i],397 (int)iomodel->gridonsurface[i],398 (int)iomodel->gridoniceshelf[i],399 (int)iomodel->gridonicesheet[i]);400 401 if (strcmp(iomodel->meshtype,"3d")==0){402 if (isnan(iomodel->uppernodes[i])){403 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.404 }405 else{406 node_upper_node_id=(int)iomodel->uppernodes[i];407 }408 }409 else{410 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/411 node_upper_node_id=node_id;412 }413 414 /*Create node using its constructor: */415 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);416 417 /*set single point constraints.: */418 if (strcmp(iomodel->meshtype,"3d")==0){419 /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */420 if (iomodel->gridonbed[i]==0){421 for(k=1;k<=node_numdofs;k++){422 node->FreezeDof(k);423 }424 }425 }426 /*Add node to nodes dataset: */427 nodes->AddObject(node);428 429 #ifdef _PARALLEL_430 } //if((my_grids[i]==1))431 #endif432 123 } 433 434 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these435 * datasets, it will not be redone: */436 elements->Presort();437 nodes->Presort();438 vertices->Presort();439 materials->Presort();440 124 441 125 /*Clean fetched data: */ … … 452 136 xfree((void**)&iomodel->gridoniceshelf); 453 137 454 455 /*Keep partitioning information into iomodel*/ 456 iomodel->epart=epart; 457 iomodel->my_grids=my_grids; 458 iomodel->my_bordergrids=my_bordergrids; 459 460 /*Free ressources:*/ 461 #ifdef _PARALLEL_ 462 xfree((void**)&all_numgrids); 463 xfree((void**)&npart); 464 VecFree(&gridborder); 465 #endif 138 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 139 * datasets, it will not be redone: */ 140 elements->Presort(); 141 nodes->Presort(); 142 vertices->Presort(); 143 materials->Presort(); 466 144 467 145 cleanup_and_return:
Note:
See TracChangeset
for help on using the changeset viewer.
