![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3428
- Timestamp:
- 04/07/10 16:10:03 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp
r3417 r3428 14 14 void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 15 15 16 17 16 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 18 19 20 17 int i,j,k,n; 21 extern int my_rank;22 extern int num_procs;23 18 24 19 /*DataSets: */ … … 28 23 DataSet* materials = NULL; 29 24 30 /*Objects: */31 Node* node = NULL;32 Vertex* vertex = NULL;33 Penta* penta = NULL;34 Matice* matice = NULL;35 Matpar* matpar = NULL;36 ElementProperties* penta_properties=NULL;37 38 /*output: */39 int* epart=NULL; //element partitioning.40 int* npart=NULL; //node partitioning.41 int* my_grids=NULL;42 double* my_bordergrids=NULL;43 44 45 /*intermediary: */46 int elements_width; //size of elements47 double B_avg;48 49 /*matice constructor input: */50 int matice_mid;51 double matice_B;52 double matice_n;53 54 /*penta constructor input: */55 56 int penta_id;57 int penta_matice_id;58 int penta_matpar_id;59 int penta_numpar_id;60 int penta_node_ids[6];61 double penta_h[6];62 double penta_s[6];63 double penta_b[6];64 double penta_k[6];65 int penta_friction_type;66 double penta_p;67 double penta_q;68 int penta_shelf;69 int penta_onbed;70 int penta_onsurface;71 int penta_collapse;72 double penta_geothermalflux[6];73 bool penta_onwater;74 75 /*matpar constructor input: */76 int matpar_mid;77 double matpar_rho_ice;78 double matpar_rho_water;79 double matpar_heatcapacity;80 double matpar_thermalconductivity;81 double matpar_latentheat;82 double matpar_beta;83 double matpar_meltingpoint;84 double matpar_mixed_layer_capacity;85 double matpar_thermal_exchange_velocity;86 double matpar_g;87 88 /* node constructor input: */89 int node_id;90 int vertex_id;91 int node_partitionborder=0;92 int node_onbed;93 int node_onsurface;94 int node_onshelf;95 int node_onsheet;96 int node_upper_node_id;97 int node_numdofs;98 99 100 #ifdef _PARALLEL_101 /*Metis partitioning: */102 int range;103 Vec gridborder=NULL;104 int my_numgrids;105 int* all_numgrids=NULL;106 int gridcount;107 int count;108 #endif109 int first_grid_index;110 111 25 /*First create the elements, nodes and material properties: */ 112 26 elements = new DataSet(ElementsEnum()); … … 115 29 materials = new DataSet(MaterialsEnum()); 116 30 117 /*return if 2d mesh*/ 118 if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; 119 120 /*Width of elements: */ 121 elements_width=6; //penta elements 122 123 #ifdef _PARALLEL_ 124 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 125 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 126 127 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 128 129 /*Free elements and elements2d: */ 130 xfree((void**)&iomodel->elements2d); 131 132 /*Used later on: */ 133 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); 134 #endif 135 136 137 /*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); 138 33 139 34 /*Fetch data needed: */ … … 155 50 156 51 for (i=0;i<iomodel->numberofelements;i++){ 157 #ifdef _PARALLEL_ 158 /*We are using our element partition to decide which elements will be created on this node: */ 159 if(my_rank==epart[i]){ 160 #endif 52 if(iomodel->my_elements[i]){ 161 53 162 163 /*name and id: */ 164 penta_id=i+1; //matlab indexing. 165 penta_matice_id=i+1; //refers to the corresponding material property card 166 penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card 167 penta_numpar_id=1; 54 /*Create and add tria element to elements dataset: */ 55 elements->AddObject(new Tria(i,iomodel)); 168 56 169 /*vertices,thickness,surface,bed and drag: */ 170 for(j=0;j<6;j++){ 171 penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); 172 penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 173 penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 174 penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 175 penta_k[j]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 176 penta_geothermalflux[j]=*(iomodel->geothermalflux+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 57 /*Create and add material property to materials dataset: */ 58 materials->AddObject(new Matice(i,iomodel,6)); 177 59 } 178 179 /*basal drag:*/180 penta_friction_type=(int)iomodel->drag_type;181 182 penta_p=iomodel->p[i];183 penta_q=iomodel->q[i];184 185 /*diverse: */186 penta_shelf=(int)*(iomodel->elementoniceshelf+i);187 penta_onbed=(int)*(iomodel->elementonbed+i);188 penta_onsurface=(int)*(iomodel->elementonsurface+i);189 penta_onwater=(bool)*(iomodel->elementonwater+i);190 191 /*We need the field collapse for transient, so that we can use compute B with the average temperature*/192 if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.193 penta_collapse=1;194 }195 else{196 penta_collapse=0;197 }198 199 /*Create properties: */200 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, NULL, NULL, penta_geothermalflux, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);201 202 /*Create Penta using its constructor:*/203 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);204 205 /*Delete properties: */206 delete penta_properties;207 208 /*Add penta element to elements dataset: */209 elements->AddObject(penta);210 211 /*Deal with material:*/212 matice_mid=i+1; //same as the material id from the geom2 elements.213 /*Average B over 6 element grids: */214 B_avg=0;215 for(j=0;j<6;j++){216 B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));217 }218 B_avg=B_avg/6;219 matice_B= B_avg;220 matice_n=(double)*(iomodel->n+i);221 222 /*Create matice using its constructor:*/223 matice= new Matice(matice_mid,matice_B,matice_n);224 225 /*Add matice element to materials dataset: */226 materials->AddObject(matice);227 228 #ifdef _PARALLEL_229 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use230 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)231 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids232 will hold which grids belong to this partition*/233 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;234 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;235 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;236 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;237 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;238 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;239 #endif240 241 #ifdef _PARALLEL_242 }//if(my_rank==epart[i])243 #endif244 245 60 }//for (i=0;i<numberofelements;i++) 246 61 … … 262 77 xfree((void**)&iomodel->elementonwater); 263 78 264 /*Add one constant material property to materials: */ 265 matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials 266 matpar_g=iomodel->g; 267 matpar_rho_ice=iomodel->rho_ice; 268 matpar_rho_water=iomodel->rho_water; 269 matpar_thermalconductivity=iomodel->thermalconductivity; 270 matpar_heatcapacity=iomodel->heatcapacity; 271 matpar_latentheat=iomodel->latentheat; 272 matpar_beta=iomodel->beta; 273 matpar_meltingpoint=iomodel->meltingpoint; 274 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 275 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 276 277 /*Create matpar object using its constructor: */ 278 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, 279 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, 280 matpar_thermal_exchange_velocity,matpar_g); 281 282 /*Add to materials datset: */ 283 materials->AddObject(matpar); 79 /*Add new constrant material property tgo materials, at the end: */ 80 materials->AddObject(new Matpar(iomodel)); 284 81 285 #ifdef _PARALLEL_ 286 /*From the element partitioning, we can determine which grids are on the inside of this cpu's 287 *element partition, and which are on its border with other nodes:*/ 288 gridborder=NewVec(iomodel->numberofnodes); 289 290 for (i=0;i<iomodel->numberofnodes;i++){ 291 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); 292 } 293 VecAssemblyBegin(gridborder); 294 VecAssemblyEnd(gridborder); 295 296 #ifdef _ISSM_DEBUG_ 297 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); 298 #endif 299 300 VecToMPISerial(&my_bordergrids,gridborder); 301 302 #ifdef _ISSM_DEBUG_ 303 if(my_rank==0){ 304 for (i=0;i<iomodel->numberofnodes;i++){ 305 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); 306 } 307 } 308 #endif 309 #endif 310 311 /*Partition penalties in 3d: */ 312 if(strcmp(iomodel->meshtype,"3d")==0){ 313 314 /*Get penalties: */ 315 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); 316 317 if(iomodel->numpenalties){ 318 319 iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); 320 #ifdef _SERIAL_ 321 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1; 322 #else 323 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1; 324 325 for(i=0;i<iomodel->numpenalties;i++){ 326 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1); 327 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids 328 /*All grids that are being penalised belong to this node's internal grid partition.:*/ 329 iomodel->penaltypartitioning[i]=1; 330 } 331 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border 332 iomodel->penaltypartitioning[i]=0; 333 } 334 } 335 #endif 336 } 337 338 /*Free penalties: */ 339 xfree((void**)&iomodel->penalties); 340 } 341 342 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 343 We can therefore determine which grids are internal to this node's partition 344 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 345 that, go and create the grids*/ 346 347 /*Create nodes from x,y,z, as well as the spc values on those grids: */ 348 349 /*First fetch data: */ 82 /*Create nodes and vertices: */ 350 83 if (strcmp(iomodel->meshtype,"3d")==0){ 351 84 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); … … 362 95 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 363 96 364 /*Get number of dofs per node: */ 365 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); 97 for (i=0;i<iomodel->numberofvertices;i++){ 366 98 367 for (i=0;i<iomodel->numberofnodes;i++){368 #ifdef _PARALLEL_369 /*keep only this partition's nodes:*/370 if((my_grids[i]==1)){371 #endif99 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 100 if(iomodel->my_vertices[i]){ 101 102 /*Add vertex to vertices dataset: */ 103 vertices->AddObject(new Vertex(i,iomodel)); 372 104 373 /*create vertex: */ 374 vertex_id=i+1; 375 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); 376 vertices->AddObject(vertex); 105 /*Add node to nodes dataset: */ 106 nodes->AddObject(new Node(i,iomodel)); 377 107 378 379 node_id=i+1; //matlab indexing380 381 #ifdef _PARALLEL_382 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border383 node_partitionborder=1;384 108 } 385 else{386 node_partitionborder=0;387 }388 #else389 node_partitionborder=0;390 #endif391 392 node_properties.SetProperties(393 (int)iomodel->gridonbed[i],394 (int)iomodel->gridonsurface[i],395 (int)iomodel->gridoniceshelf[i],396 (int)iomodel->gridonicesheet[i]);397 398 399 if (strcmp(iomodel->meshtype,"3d")==0){400 if (isnan(iomodel->uppernodes[i])){401 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.402 }403 else{404 node_upper_node_id=(int)iomodel->uppernodes[i];405 }406 }407 else{408 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/409 node_upper_node_id=node_id;410 }411 412 /*Create node using its constructor: */413 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);414 415 /*Add node to nodes dataset: */416 nodes->AddObject(node);417 418 #ifdef _PARALLEL_419 } //if((my_grids[i]==1))420 #endif421 109 } 422 423 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these424 * datasets, it will not be redone: */425 elements->Presort();426 nodes->Presort();427 vertices->Presort();428 materials->Presort();429 110 430 111 /*Clean fetched data: */ … … 441 122 xfree((void**)&iomodel->gridoniceshelf); 442 123 124 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 125 * datasets, it will not be redone: */ 126 elements->Presort(); 127 nodes->Presort(); 128 vertices->Presort(); 129 materials->Presort(); 130 443 131 cleanup_and_return: 444 132 … … 449 137 *pmaterials=materials; 450 138 451 /*Keep partitioning information into iomodel*/452 iomodel->epart=epart;453 iomodel->my_grids=my_grids;454 iomodel->my_bordergrids=my_bordergrids;455 456 /*Free ressources:*/457 #ifdef _PARALLEL_458 xfree((void**)&all_numgrids);459 xfree((void**)&npart);460 VecFree(&gridborder);461 #endif462 463 139 }
Note:
See TracChangeset
for help on using the changeset viewer.
