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