![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3421
- Timestamp:
- 04/07/10 15:14:36 (15 years ago)
- Location:
- issm/trunk/src/c/ModelProcessorx
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
r3420 r3421 12 12 #include "../IoModel.h" 13 13 14 15 14 void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 16 15 17 16 /*Intermediary*/ 18 17 int i,j,k,n; 19 18 … … 146 145 147 146 } //if (strcmp(meshtype,"2d")==0) 148 149 147 150 148 /*Add new constrant material property tgo materials, at the end: */ … … 166 164 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 167 165 } 168 169 166 170 167 for (i=0;i<iomodel->numberofvertices;i++){ -
issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
r3417 r3421 2 2 * CreateElementsNodesAndMaterialsDiagnosticVert.c: 3 3 */ 4 5 4 6 5 #include "../../DataSet/DataSet.h" … … 13 12 #include "../IoModel.h" 14 13 15 16 14 void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 17 15 18 19 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 20 21 16 /*Intermediary*/ 22 17 int i,j,k,n; 23 extern int my_rank;24 extern int num_procs;25 18 26 19 /*DataSets: */ … … 36 29 Matice* matice = NULL; 37 30 Matpar* matpar = NULL; 38 ElementProperties* penta_properties=NULL;39 31 40 /*output: */ 41 int* epart=NULL; //element partitioning. 42 int* npart=NULL; //node partitioning. 43 int* my_grids=NULL; 44 double* my_bordergrids=NULL; 45 46 47 /*intermediary: */ 48 const int elements_width=6; 49 double B_avg; 50 51 /*matice constructor input: */ 52 int matice_mid; 53 double matice_B; 54 double matice_n; 55 56 /*penta constructor input: */ 57 58 int penta_id; 59 int penta_matice_id; 60 int penta_matpar_id; 61 int penta_numpar_id; 62 int penta_node_ids[6]; 63 double penta_h[6]; 64 double penta_s[6]; 65 double penta_b[6]; 66 int penta_shelf; 67 int penta_onbed; 68 int penta_onsurface; 69 int penta_collapse; 70 double penta_melting[6]; 71 double penta_accumulation[6]; 72 int penta_artdiff; 73 double penta_viscosity_overshoot; 74 double penta_stokesreconditioning; 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 #endif 111 int first_grid_index; 112 113 /*Penalty partitioning: */ 114 int num_grids3d_collapsed; 115 double* double_penalties_grids3d_collapsed=NULL; 116 double* double_penalties_grids3d_noncollapsed=NULL; 117 int grid_ids[6]; 118 int num_grid_ids; 119 int grid_id; 120 121 32 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 33 if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; 122 34 123 35 /*First create the elements, nodes and material properties: */ … … 127 39 materials = new DataSet(MaterialsEnum()); 128 40 129 /*Now, is the iomodel running in 3d? : */ 130 if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; 131 132 #ifdef _PARALLEL_ 133 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 134 if(strcmp(iomodel->meshtype,"2d")==0){ 135 /*load elements: */ 136 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 137 } 138 else{ 139 /*load elements2d: */ 140 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 141 } 142 143 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 144 145 /*Free elements and elements2d: */ 146 xfree((void**)&iomodel->elements); 147 xfree((void**)&iomodel->elements2d); 148 149 /*Used later on: */ 150 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); 151 #endif 152 41 /*Partition elements and vertices and nodes: */ 42 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 153 43 154 44 /*Create 3d elements: */ … … 167 57 168 58 for (i=0;i<iomodel->numberofelements;i++){ 169 #ifdef _PARALLEL_170 /*We are using our element partition to decide which elements will be created on this node: */171 if(my_rank==epart[i]){172 #endif173 59 174 175 /*name and id: */ 176 penta_id=i+1; //matlab indexing. 177 penta_matice_id=i+1; //refers to the corresponding material property card 178 penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card 179 penta_numpar_id=1; 60 if(iomodel->my_elements[i]){ 180 61 181 /*vertices,thickness,surface,bed and drag: */ 182 for(j=0;j<6;j++){ 183 penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); 184 penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 185 penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 186 penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); 187 penta_melting[j]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+j)-1))/iomodel->yts; 188 penta_accumulation[j]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+j)-1))/iomodel->yts; 62 /*Create and add penta element to elements dataset: */ 63 elements->AddObject(new Penta(i,iomodel)); 64 65 /*Create and add material property to materials dataset: */ 66 materials->AddObject(new Matice(i,iomodel,6)); 67 189 68 } 190 191 /*diverse: */192 penta_shelf=(int)*(iomodel->elementoniceshelf+i);193 penta_onbed=(int)*(iomodel->elementonbed+i);194 penta_onsurface=(int)*(iomodel->elementonsurface+i);195 penta_collapse=1;196 penta_onwater=(bool)*(iomodel->elementonwater+i);197 198 /*Create properties: */199 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, NULL, penta_melting, penta_accumulation, NULL, UNDEF,UNDEF,UNDEF, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);200 201 /*Create Penta using its constructor:*/202 penta=new Penta(penta_id,penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);203 204 /*delete properties: */205 delete penta_properties;206 207 /*Add penta element to elements dataset: */208 elements->AddObject(penta);209 210 211 /*Deal with material:*/212 matice_mid=i+1; //same as the material id from the geom2 elements.213 214 /*Create matice using its constructor:*/215 matice= new Matice(matice_mid,matice_B,matice_n);216 217 /*Add matice element to materials dataset: */218 materials->AddObject(matice);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 use222 *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_grids224 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 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;229 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;230 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;231 #endif232 233 #ifdef _PARALLEL_234 }//if(my_rank==epart[i])235 #endif236 69 237 70 }//for (i=0;i<numberofelements;i++) … … 250 83 xfree((void**)&iomodel->elementonwater); 251 84 252 253 /*Add one constant material property to materials: */ 254 matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials 255 matpar_g=iomodel->g; 256 matpar_rho_ice=iomodel->rho_ice; 257 matpar_rho_water=iomodel->rho_water; 258 matpar_thermalconductivity=iomodel->thermalconductivity; 259 matpar_heatcapacity=iomodel->heatcapacity; 260 matpar_latentheat=iomodel->latentheat; 261 matpar_beta=iomodel->beta; 262 matpar_meltingpoint=iomodel->meltingpoint; 263 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 264 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 265 266 /*Create matpar object using its constructor: */ 267 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, 268 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, 269 matpar_thermal_exchange_velocity,matpar_g); 270 271 /*Add to materials datset: */ 272 materials->AddObject(matpar); 85 /*Add new constrant material property tgo materials, at the end: */ 86 materials->AddObject(new Matpar(iomodel)); 273 87 274 #ifdef _PARALLEL_275 /*From the element partitioning, we can determine which grids are on the inside of this cpu's276 *element partition, and which are on its border with other nodes:*/277 gridborder=NewVec(iomodel->numberofnodes);278 279 for (i=0;i<iomodel->numberofnodes;i++){280 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);281 }282 VecAssemblyBegin(gridborder);283 VecAssemblyEnd(gridborder);284 285 #ifdef _ISSM_DEBUG_286 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);287 #endif288 289 VecToMPISerial(&my_bordergrids,gridborder);290 291 #ifdef _ISSM_DEBUG_292 if(my_rank==0){293 for (i=0;i<iomodel->numberofnodes;i++){294 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);295 }296 }297 #endif298 #endif299 300 301 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.302 We can therefore determine which grids are internal to this node's partition303 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing304 that, go and create the grids*/305 306 /*Create nodes from x,y,z, as well as the spc values on those grids: */307 308 88 /*First fetch data: */ 309 if (strcmp(iomodel->meshtype,"3d")==0){ 310 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 311 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 312 } 89 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 90 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 313 91 IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); 314 92 IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); … … 321 99 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 322 100 323 324 /*Get number of dofs per node: */ 325 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); 101 for (i=0;i<iomodel->numberofvertices;i++){ 326 102 327 for (i=0;i<iomodel->numberofnodes;i++){ 328 #ifdef _PARALLEL_ 329 /*keep only this partition's nodes:*/ 330 if((my_grids[i]==1)){ 331 #endif 103 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 104 if(iomodel->my_vertices[i]){ 332 105 333 /*create vertex: */ 334 vertex_id=i+1; 335 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); 106 /*Add vertex to vertices dataset: */ 107 vertices->AddObject(new Vertex(i,iomodel)); 336 108 337 vertices->AddObject(vertex); 109 /*Add node to nodes dataset: */ 110 nodes->AddObject(new Node(i,iomodel)); 338 111 339 340 node_id=i+1; //matlab indexing341 342 #ifdef _PARALLEL_343 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border344 node_partitionborder=1;345 112 } 346 else{347 node_partitionborder=0;348 }349 #else350 node_partitionborder=0;351 #endif352 353 node_properties.SetProperties(354 (int)iomodel->gridonbed[i],355 (int)iomodel->gridonsurface[i],356 (int)iomodel->gridoniceshelf[i],357 (int)iomodel->gridonicesheet[i]);358 359 if (isnan(iomodel->uppernodes[i])){360 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.361 }362 else{363 node_upper_node_id=(int)iomodel->uppernodes[i];364 }365 366 /*Create node using its constructor: */367 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);368 369 /*Add node to nodes dataset: */370 nodes->AddObject(node);371 372 #ifdef _PARALLEL_373 } //if((my_grids[i]==1))374 #endif375 113 } 376 377 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these378 * datasets, it will not be redone: */379 elements->Presort();380 nodes->Presort();381 vertices->Presort();382 materials->Presort();383 114 384 115 /*Clean fetched data: */ … … 395 126 xfree((void**)&iomodel->gridoniceshelf); 396 127 397 398 /*Keep partitioning information into iomodel*/ 399 iomodel->epart=epart; 400 iomodel->my_grids=my_grids; 401 iomodel->my_bordergrids=my_bordergrids; 402 403 /*Free ressources:*/ 404 #ifdef _PARALLEL_ 405 xfree((void**)&all_numgrids); 406 xfree((void**)&npart); 407 VecFree(&gridborder); 408 #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(); 409 134 410 135 cleanup_and_return:
Note:
See TracChangeset
for help on using the changeset viewer.
