Changeset 3427
- Timestamp:
- 04/07/10 16:07:57 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp ¶
r3417 r3427 15 15 void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 16 16 17 /*Intermediary*/ 17 18 int i,j,k,n; 18 extern int my_rank;19 extern int num_procs;20 19 21 20 /*DataSets: */ 22 21 DataSet* elements = NULL; 23 22 DataSet* nodes = NULL; 24 DataSet* 23 DataSet* vertices = NULL; 25 24 DataSet* materials = NULL; 26 27 /*Objects: */28 Node* node = NULL;29 Vertex* vertex = NULL;30 Tria* tria = NULL;31 Penta* penta = NULL;32 Matice* matice = NULL;33 Matpar* matpar = NULL;34 ElementProperties* tria_properties=NULL;35 36 /*output: */37 int* epart=NULL; //element partitioning.38 int* npart=NULL; //node partitioning.39 int* my_grids=NULL;40 double* my_bordergrids=NULL;41 42 /*intermediary: */43 int elements_width; //size of elements44 double B_avg;45 46 /*tria constructor input: */47 int tria_id;48 int tria_matice_id;49 int tria_matpar_id;50 int tria_numpar_id;51 int tria_node_ids[3];52 double tria_h[3];53 double tria_s[3];54 double tria_b[3];55 int tria_shelf;56 bool tria_onwater;57 58 /*matice constructor input: */59 int matice_mid;60 double matice_B;61 double matice_n;62 63 /*penta constructor input: */64 65 int penta_id;66 int penta_mid;67 int penta_mparid;68 int penta_numparid;69 int penta_g[6];70 double penta_h[6];71 double penta_s[6];72 double penta_b[6];73 double penta_k[6];74 double penta_p;75 double penta_q;76 int penta_shelf;77 int penta_onbed;78 int penta_onsurface;79 int penta_collapse;80 double penta_melting[6];81 double penta_accumulation[6];82 int penta_thermal_steadystate;83 bool penta_onwater;84 85 /*matpar constructor input: */86 int matpar_mid;87 double matpar_rho_ice;88 double matpar_rho_water;89 double matpar_heatcapacity;90 double matpar_thermalconductivity;91 double matpar_latentheat;92 double matpar_beta;93 double matpar_meltingpoint;94 double matpar_mixed_layer_capacity;95 double matpar_thermal_exchange_velocity;96 double matpar_g;97 98 /* node constructor input: */99 int node_id;100 int vertex_id;101 int vertex_partitionborder=0;102 int node_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 /*Nodes cloning*/120 double e1,e2;121 int i1,i2;122 int pos;123 #endif124 int first_grid_index;125 25 126 26 /*First create the elements, nodes and material properties: */ … … 130 30 materials = new DataSet(MaterialsEnum()); 131 31 132 /*Width of elements: */ 133 if(strcmp(iomodel->meshtype,"2d")==0){ 134 elements_width=3; //tria elements 135 } 136 else{ 137 ISSMERROR("not implemented yet!"); 138 elements_width=6; //penta elements 139 } 140 141 #ifdef _PARALLEL_ 142 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 143 if(strcmp(iomodel->meshtype,"2d")==0){ 144 /*load elements: */ 145 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 146 } 147 else{ 148 /*load elements2d: */ 149 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 150 } 151 152 153 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 154 155 /*Free elements and elements2d: */ 156 xfree((void**)&iomodel->elements); 157 xfree((void**)&iomodel->elements2d); 158 159 /*Used later on: */ 160 my_grids=(int*)xcalloc(3*iomodel->numberofelements,sizeof(int)); 161 #endif 32 /*Partition elements and vertices and nodes: */ 33 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 162 34 163 35 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ 164 165 36 /*2d mesh: */ 166 37 if (strcmp(iomodel->meshtype,"2d")==0){ … … 176 47 for (i=0;i<iomodel->numberofelements;i++){ 177 48 178 #ifdef _PARALLEL_ 179 /*!All elements have been partitioned above, only create elements for this CPU: */ 180 if(my_rank==epart[i]){ 181 #endif 182 183 /*ids: */ 184 tria_id=i+1; //matlab indexing. 185 tria_matice_id=-1; //no need for materials 186 tria_matpar_id=-1; //no need for materials 187 tria_numpar_id=1; 49 if(iomodel->my_elements[i]){ 188 50 189 /*vertices offsets: */ 190 tria_node_ids[0]=3*i+1; 191 tria_node_ids[1]=3*i+2; 192 tria_node_ids[2]=3*i+3; 51 /*Create and add tria element to elements dataset: */ 52 elements->AddObject(new Tria(i,iomodel)); 193 53 194 /*thickness,surface and bed:*/ 195 tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. 196 tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 197 tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ; 198 199 tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 200 tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 201 tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 202 203 tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); 204 tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); 205 tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); 206 207 /*element on iceshelf?:*/ 208 tria_shelf=(int)*(iomodel->elementoniceshelf+i); 209 tria_onwater=(bool)*(iomodel->elementonwater+i); 210 211 /*Create properties: */ 212 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); 213 214 /*Create tria element using its constructor:*/ 215 tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties); 216 217 /*Delete properties: */ 218 delete tria_properties; 219 220 /*Add tria element to elements dataset: */ 221 elements->AddObject(tria); 222 223 #ifdef _PARALLEL_ 224 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 225 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 226 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 227 will hold which grids belong to this partition*/ 228 my_grids[3*i+0]=1; 229 my_grids[3*i+1]=1; 230 my_grids[3*i+2]=1; 231 #endif 232 233 #ifdef _PARALLEL_ 234 }//if(my_rank==epart[i]) 235 #endif 236 54 /*Create and add material property to materials dataset: */ 55 materials->AddObject(new Matice(i,iomodel,3)); 56 } 237 57 }//for (i=0;i<numberofelements;i++) 238 58 … … 249 69 } //if (strcmp(meshtype,"2d")==0) 250 70 251 #ifdef _PARALLEL_ 252 /*If we are in parallel, we must add the nodes that are not in this partition 253 * but share edges of elements that are in this partition. This is needed for 254 * the loads as numerical fluxes involve the dofs of all nodes on a given edge*/ 71 /*Add new constrant material property tgo materials, at the end: */ 72 materials->AddObject(new Matpar(iomodel)); 255 73 256 /*Get edges and elements*/ 257 IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges"); 258 259 /*!All elements have been partitioned above, only create elements for this CPU: */ 260 for (i=0;i<iomodel->numberofedges;i++){ 261 262 /*Get left and right elements*/ 263 e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] 264 e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] 265 266 /* 1) If the element e1 is in the current partition 267 * 2) and if the edge of the element is shared by another element 268 * 3) and if this element is not in the same partition: 269 * we must clone the nodes on this partition so that the loads (Numericalflux) 270 * will have access to their properties (dofs,...)*/ 271 if(my_rank==epart[(int)e1] && !isnan(e2) && my_rank!=epart[(int)e2]){ 272 273 /*1: Get vertices ids*/ 274 i1=(int)iomodel->edges[4*i+0]; 275 i2=(int)iomodel->edges[4*i+1]; 276 277 /*2: Get the column where these ids are located in the index*/ 278 pos==UNDEF; 279 for(j=0;j<3;j++){ 280 if (iomodel->elements[3*(int)e2+j]==i1) pos=j+1; 281 } 282 ISSMASSERT(pos!=UNDEF); 283 284 /*3: We have the id of the elements and the position of the vertices in the index 285 * we can now create the corresponding nodes:*/ 286 my_grids[(int)e2*3+pos-1]=1; 287 my_grids[(int)e2*3+((pos+1)%3)]=1; 288 } 289 } 290 291 /*Free data: */ 292 xfree((void**)&iomodel->edges); 293 294 /*From the element partitioning, we can determine which nodes are on the inside of this cpu's 295 *element partition, and which are on its border with other nodes:*/ 296 gridborder=NewVec(3*iomodel->numberofelements); 297 298 for (i=0;i<3*iomodel->numberofelements;i++){ 299 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); 300 } 301 VecAssemblyBegin(gridborder); 302 VecAssemblyEnd(gridborder); 303 304 VecToMPISerial(&my_bordergrids,gridborder); 305 306 #endif 307 308 /*Add one constant material property to materials: */ 309 matpar_mid=1; //put it at the end of the materials 310 matpar_g=iomodel->g; 311 matpar_rho_ice=iomodel->rho_ice; 312 matpar_rho_water=iomodel->rho_water; 313 matpar_thermalconductivity=iomodel->thermalconductivity; 314 matpar_heatcapacity=iomodel->heatcapacity; 315 matpar_latentheat=iomodel->latentheat; 316 matpar_beta=iomodel->beta; 317 matpar_meltingpoint=iomodel->meltingpoint; 318 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 319 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 320 321 /*Create matpar object using its constructor: */ 322 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, 323 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, 324 matpar_thermal_exchange_velocity,matpar_g); 325 326 /*Add to materials datset: */ 327 materials->AddObject(matpar); 328 329 /*Ok, let's summarise. Now, every CPU has the following array: my_grids 330 We can therefore determine which grids are internal to this node's partition 331 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 332 that, go and create the grids*/ 333 334 /*Create nodes from x,y,z, as well as the spc values on those grids: */ 335 336 /*First fetch data: */ 337 if (strcmp(iomodel->meshtype,"3d")==0){ 338 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 339 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 340 } 74 /*Create nodes and vertices: */ 341 75 IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); 342 76 IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); … … 346 80 IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); 347 81 IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); 82 IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); 348 83 IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); 349 84 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 350 351 /*Get number of dofs per node: */ 352 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); 353 354 /*Create vertices: */ 355 for (i=0;i<iomodel->numberofnodes;i++){ 356 #ifdef _PARALLEL_ 357 /*keep only this partition's nodes:*/ 358 if((my_grids[i]==1)){ 359 #endif 360 361 /*create vertex: */ 362 vertex_id=i+1; 363 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); 364 365 vertices->AddObject(vertex); 366 367 #ifdef _PARALLEL_ 368 } //if((my_grids[i]==1)) 369 #endif 85 if (strcmp(iomodel->meshtype,"3d")==0){ 86 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 87 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 370 88 } 371 89 372 /*Build Nodes dataset -> 3 for each element*/ 90 /*Build Vertices dataset*/ 91 for (i=0;i<iomodel->numberofvertices;i++){ 92 93 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 94 if(iomodel->my_vertices[i]){ 95 96 /*Add vertex to vertices dataset: */ 97 vertices->AddObject(new Vertex(i,iomodel)); 98 99 } 100 } 101 102 /*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/ 373 103 for (i=0;i<iomodel->numberofelements;i++){ 374 104 for (j=0;j<3;j++){ 375 105 376 #ifdef _PARALLEL_ 377 /*!All elements have been partitioned above, only create elements for this CPU: */ 378 if(my_grids[3*i+j]){ 379 #endif 106 if(my_nodes[3*i+j]){ 380 107 381 108 //Get id of the vertex on which the current node is located 382 109 k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing) 383 ISSMASSERT(k>=0 && k<iomodel->numberofnodes); 384 385 //Get node properties 386 node_id=i*3+j+1; 387 #ifdef _PARALLEL_ 388 if(my_bordergrids[node_id-1]>1.0) { //this grid belongs to a partition border 389 node_partitionborder=1; 390 } 391 else{ 392 node_partitionborder=0; 393 } 394 #else 395 node_partitionborder=0; 396 #endif 397 398 node_properties.SetProperties( 399 (int)iomodel->gridonbed[k], 400 (int)iomodel->gridonsurface[k], 401 (int)iomodel->gridoniceshelf[k], 402 (int)iomodel->gridonicesheet[k]); 403 404 if (strcmp(iomodel->meshtype,"3d")==0){ 405 if (isnan(iomodel->uppernodes[k])){ 406 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. 407 } 408 else{ 409 node_upper_node_id=(int)iomodel->uppernodes[k]; 410 } 411 } 412 else{ 413 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ 414 node_upper_node_id=node_id; 415 } 416 417 /*Create node using its constructor: */ 418 node=new Node(node_id, k, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); 110 ISSMASSERT(k>=0 && k<iomodel->numberofvertices); 419 111 420 112 /*Add node to nodes dataset: */ 421 nodes->AddObject(n ode);422 #ifdef _PARALLEL_ 113 nodes->AddObject(new Node(k,i,iomodel)); 114 423 115 } 424 #endif425 116 } 426 117 } 118 119 /*Clean fetched data: */ 120 xfree((void**)&iomodel->deadgrids); 121 xfree((void**)&iomodel->x); 122 xfree((void**)&iomodel->y); 123 xfree((void**)&iomodel->z); 124 xfree((void**)&iomodel->thickness); 125 xfree((void**)&iomodel->bed); 126 xfree((void**)&iomodel->gridonbed); 127 xfree((void**)&iomodel->gridonsurface); 128 xfree((void**)&iomodel->gridonhutter); 129 xfree((void**)&iomodel->uppernodes); 130 xfree((void**)&iomodel->gridonicesheet); 131 xfree((void**)&iomodel->gridoniceshelf); 427 132 428 133 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these … … 432 137 vertices->Presort(); 433 138 materials->Presort(); 434 435 /*Clean fetched data: */436 xfree((void**)&iomodel->deadgrids);437 xfree((void**)&iomodel->elements);438 xfree((void**)&iomodel->x);439 xfree((void**)&iomodel->y);440 xfree((void**)&iomodel->z);441 xfree((void**)&iomodel->thickness);442 xfree((void**)&iomodel->bed);443 xfree((void**)&iomodel->gridonbed);444 xfree((void**)&iomodel->gridonsurface);445 xfree((void**)&iomodel->uppernodes);446 xfree((void**)&iomodel->gridonicesheet);447 xfree((void**)&iomodel->gridoniceshelf);448 449 /*Keep partitioning information into iomodel*/450 iomodel->epart=epart;451 iomodel->my_grids=my_grids;452 iomodel->my_bordergrids=my_bordergrids;453 454 /*Free ressources:*/455 #ifdef _PARALLEL_456 xfree((void**)&all_numgrids);457 xfree((void**)&npart);458 VecFree(&gridborder);459 #endif460 139 461 140 cleanup_and_return: -
TabularUnified issm/trunk/src/c/objects/Node.cpp ¶
r3423 r3427 52 52 } 53 53 /*}}}*/ 54 /*FUNCTION Node constructor from iomodel{{{2*/54 /*FUNCTION Node constructor from iomodel continuous Galerkin{{{2*/ 55 55 Node::Node(int i, IoModel* iomodel){ //i is the node index 56 56 … … 131 131 } 132 132 } 133 } 134 /*}}}*/ 135 /*FUNCTION Node constructor from iomodel discontinuous Galerkin{{{2*/ 136 Node::Node(int i,int j,IoModel* iomodel){ 137 /* i -> index of the vertex in C indexing 138 * j -> index of the node in C indexing*/ 139 140 int numdofs; 141 int partitionborder; 142 int vertex_id; 143 int upper_node_id; 144 145 /*id: */ 146 this->id=j+1; //matlab indexing 147 148 /*indexing:*/ 149 DistributeNumDofs(&numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); //number of dofs per node 150 if(iomodel->my_bordervertices[i])partitionborder=1; else partitionborder=0;//is this node on a partition border? 151 152 this->indexing.Init(numdofs,partitionborder); 153 154 /*properties (come from vertex number i): */ 155 this->properties.Init( 156 (int)iomodel->gridonbed[i], 157 (int)iomodel->gridonsurface[i], 158 (int)iomodel->gridoniceshelf[i], 159 (int)iomodel->gridonicesheet[i]); 160 161 /*hooks: */ 162 vertex_id=i+1; //matlab indexing 163 164 if (strcmp(iomodel->meshtype,"3d")==0){ 165 if (isnan(iomodel->uppernodes[i])){ 166 upper_node_id=this->id; //nodes on surface do not have upper nodes, only themselves. 167 } 168 else{ 169 upper_node_id=(int)iomodel->uppernodes[i]; 170 } 171 } 172 else{ 173 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ 174 upper_node_id=this->id; 175 } 176 177 this->hvertex.Init(&vertex_id,1); 178 this->hupper_node.Init(&upper_node_id,1); 179 133 180 } 134 181 /*}}}*/ -
TabularUnified issm/trunk/src/c/objects/Node.h ¶
r3420 r3427 44 44 Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node); 45 45 Node(int i, IoModel* iomodel); 46 Node(int i,int j,IoModel* iomodel); 46 47 ~Node(); 47 48 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.