![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3423
- Timestamp:
- 04/07/10 15:50:00 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 4 edited
-
ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (modified) (1 diff)
-
ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (modified) (9 diffs)
-
ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (modified) (1 diff)
-
objects/Node.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
r3421 r3423 23 23 DataSet* materials = NULL; 24 24 25 /*Objects: */26 Node* node = NULL;27 Vertex* vertex = NULL;28 Tria* tria = NULL;29 Penta* penta = NULL;30 Matice* matice = NULL;31 Matpar* matpar = NULL;32 33 25 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 34 26 if (!iomodel->ismacayealpattyn)goto cleanup_and_return; -
issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
r3408 r3423 2 2 * CreateElementsNodesAndMaterialsDiagnosticHutter.c: 3 3 */ 4 5 4 6 5 #include "../../DataSet/DataSet.h" … … 13 12 #include "../IoModel.h" 14 13 14 void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 15 15 16 void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 17 18 19 int i,k; 20 extern int my_rank; 21 extern int num_procs; 16 /*Intermediary*/ 17 int i,j,k,n; 22 18 23 19 /*DataSets: */ 24 20 DataSet* elements = NULL; 25 21 DataSet* nodes = NULL; 22 DataSet* vertices = NULL; 26 23 DataSet* materials = NULL; 27 28 /*Objects: */29 Node* node = NULL;30 Matice* matice = NULL;31 Matpar* matpar = NULL;32 Beam* beam = NULL;33 Sing* sing = NULL;34 ElementProperties* sing_properties=NULL;35 ElementProperties* beam_properties=NULL;36 37 /*output: */38 int* epart=NULL; //element partitioning.39 int* npart=NULL; //node partitioning.40 int* my_grids=NULL;41 double* my_bordergrids=NULL;42 43 /*intermediary: */44 int elements_width; //size of elements45 double B_avg;46 47 /*matice constructor input: */48 int matice_mid;49 double matice_B;50 double matice_n;51 52 /*sing constructor input: */53 int sing_id;54 int sing_matice_id;55 int sing_matpar_id;56 int sing_numpar_id;57 int sing_node_ids[1];58 double sing_h[1],sing_k[1];59 60 /*beam constructor input: */61 int beam_id;62 int beam_matice_id;63 int beam_matpar_id;64 int beam_numpar_id;65 int beam_node_ids[2];66 double beam_h[2];67 double beam_s[2];68 double beam_b[2];69 double beam_k[2];70 bool beam_onbed;71 72 /*matpar constructor input: */73 int matpar_mid;74 double matpar_rho_ice;75 double matpar_rho_water;76 double matpar_heatcapacity;77 double matpar_thermalconductivity;78 double matpar_latentheat;79 double matpar_beta;80 double matpar_meltingpoint;81 double matpar_mixed_layer_capacity;82 double matpar_thermal_exchange_velocity;83 double matpar_g;84 85 /* node constructor input: */86 int node_id;87 int node_partitionborder=0;88 double node_x[3];89 double node_sigma;90 int node_onbed;91 int node_onsurface;92 int node_onshelf;93 int node_onsheet;94 int node_upper_node_id;95 int node_numdofs;96 97 #ifdef _PARALLEL_98 /*Metis partitioning: */99 int range;100 Vec gridborder=NULL;101 int my_numgrids;102 int* all_numgrids=NULL;103 int gridcount;104 int count;105 #endif106 107 /*First create the elements, nodes and material properties: */108 elements = new DataSet(ElementsEnum());109 nodes = new DataSet(NodesEnum());110 materials = new DataSet(MaterialsEnum());111 24 112 25 /*Now, is the flag ishutter on? otherwise, do nothing: */ 113 26 if (!iomodel->ishutter)goto cleanup_and_return; 114 115 /*Width of elements: */116 if(strcmp(iomodel->meshtype,"2d")==0){117 elements_width=3; //sing elements118 }119 else{120 elements_width=6; //beam elements121 }122 123 #ifdef _PARALLEL_124 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/125 if(strcmp(iomodel->meshtype,"2d")==0){126 /*load elements: */127 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");128 }129 else{130 /*load elements2d: */131 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");132 }133 134 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);135 136 /*Free elements and elements2d: */137 xfree((void**)&iomodel->elements);138 xfree((void**)&iomodel->elements2d);139 140 /*Used later on: */141 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));142 #endif143 144 #ifdef _PARALLEL_145 if(strcmp(iomodel->meshtype,"2d")==0){146 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");147 for (i=0;i<iomodel->numberofelements;i++){148 if(my_rank==epart[i]){149 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use150 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)151 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids152 will hold which grids belong to this partition*/153 154 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;155 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;156 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;157 }158 }159 }160 else{161 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");162 for (i=0;i<iomodel->numberofelements;i++){163 if(my_rank==epart[i]){164 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;165 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;166 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;167 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;168 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;169 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;170 }171 }172 }173 174 /*From the element partitioning, we can determine which grids are on the inside of this cpu's175 *element partition, and which are on its border with other nodes:*/176 gridborder=NewVec(iomodel->numberofnodes);177 178 for (i=0;i<iomodel->numberofnodes;i++){179 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);180 }181 VecAssemblyBegin(gridborder);182 VecAssemblyEnd(gridborder);183 184 #ifdef _ISSM_DEBUG_185 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);186 #endif187 188 VecToMPISerial(&my_bordergrids,gridborder);189 190 #ifdef _ISSM_DEBUG_191 if(my_rank==0){192 for (i=0;i<iomodel->numberofnodes;i++){193 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);194 }195 }196 #endif197 #endif198 27 199 28 /*Hutter elements can be partitioned using epart, even if … … 216 45 217 46 for (i=0;i<iomodel->numberofnodes;i++){ 218 #ifdef _PARALLEL_219 /*keep only this partition's nodes:*/220 if(my_grids[i]){221 #endif222 47 223 if(iomodel->gridonhutter[i]){ 224 225 /*Deal with sing element: */ 226 sing_id=i+1; 227 sing_matice_id=i+1; //refers to the corresponding material property card 228 sing_matpar_id=iomodel->numberofnodes+1;//refers to the corresponding matpar property card 229 sing_numpar_id=1; 230 sing_node_ids[0]=i+1; 231 sing_h[0]=iomodel->thickness[i]; 232 sing_k[0]=iomodel->drag[i]; 48 if(iomodel->my_nodes[i]){ 233 49 234 /*Create properties: */235 sing_properties=new ElementProperties(1,sing_h,NULL,NULL, sing_k,NULL,NULL, NULL, UNDEF, (double)UNDEF, (double)UNDEF, UNDEF, UNDEF,true, UNDEF,UNDEF,UNDEF);50 /*Create and add penta element to elements dataset: */ 51 elements->AddObject(new Sing(i,iomodel)); 236 52 237 /*Create sing element using its constructor:*/238 sing=new Sing(sing_id, sing_node_ids, sing_matice_id, sing_matpar_id, sing_numpar_id, sing_properties);53 /*Create and add material property to materials dataset: */ 54 materials->AddObject(new Matice(i,iomodel,1)); 239 55 240 /*delete properties: */241 delete sing_properties;242 243 /*Add sing element to elements dataset: */244 elements->AddObject(sing);245 246 /*Deal with material property card: */247 matice_mid=i+1; //same as the material id from the geom2 elements.248 matice_B=iomodel->B[i];249 matice_n=(double)iomodel->n[1]; //n defined on elements not grids, so take the first value everywhere250 251 /*Create matice using its constructor:*/252 matice= new Matice(matice_mid,matice_B,matice_n);253 254 /*Add matice element to materials dataset: */255 materials->AddObject(matice);256 56 } 257 258 #ifdef _PARALLEL_259 } //if(my_rank==npart[i])260 #endif261 57 262 58 } //for (i=0;i<iomodel->numberofnodes;i++) … … 265 61 266 62 for (i=0;i<iomodel->numberofnodes;i++){ 267 #ifdef _PARALLEL_268 /*keep only this partition's nodes:*/269 if(my_grids[i]){270 #endif271 63 272 if(iomodel->gridonhutter[i]){ 64 if(iomodel->my_nodes[i]){ 65 if(iomodel->gridonhutter[i]){ 66 if(!iomodel->gridonsurface[i]){ 273 67 274 if(!iomodel->gridonsurface[i]){ 275 276 /*Deal with sing element: */ 277 beam_id=i+1; 278 beam_matice_id=i+1; //refers to the corresponding material property card 279 beam_matpar_id=iomodel->numberofnodes-iomodel->numberofnodes2d+1;//refers to the corresponding matpar property card 280 beam_numpar_id=1; 281 beam_node_ids[0]=i+1; 282 beam_node_ids[1]=(int)iomodel->uppernodes[i]; //grid that lays right on top 283 beam_h[0]=iomodel->thickness[i]; 284 beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)]; 285 beam_s[0]=iomodel->surface[i]; 286 beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)]; 287 beam_b[0]=iomodel->bed[i]; 288 beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)]; 289 beam_k[0]=iomodel->drag[i]; 290 beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)]; 68 /*Create and add penta element to elements dataset: */ 69 elements->AddObject(new Beam(i,iomodel)); 291 70 292 beam_onbed=(bool)iomodel->gridonbed[i]; 71 /*Create and add material property to materials dataset: */ 72 materials->AddObject(new Matice(i,iomodel,2)); 293 73 294 /*Create element properties: */ 295 beam_properties=new ElementProperties(2,beam_h, beam_s, beam_b, beam_k, NULL,NULL, NULL, UNDEF, NULL, NULL, UNDEF, beam_onbed, (bool)UNDEF, UNDEF, UNDEF, UNDEF); 296 297 /*Create Beam using its constructor:*/ 298 beam= new Beam(beam_id,beam_node_ids, beam_matice_id, beam_matpar_id, beam_numpar_id, beam_properties); 299 300 /*delete properties: */ 301 delete beam_properties; 302 303 /*Add tria element to elements dataset: */ 304 elements->AddObject(beam); 305 306 /*Deal with material property card: */ 307 matice_mid=i+1; //same as the material id from the geom2 elements. 308 matice_B=iomodel->B[i]; 309 matice_n=(double)iomodel->n[1]; //n defined on elements not grids, so take the first value everywhere 310 311 /*Create matice ubeam its constructor:*/ 312 matice= new Matice(matice_mid,matice_B,matice_n); 313 314 /*Add matice element to materials dataset: */ 315 materials->AddObject(matice); 74 } 316 75 } 317 76 } 318 319 #ifdef _PARALLEL_320 } //if(my_rank==npart[i])321 #endif322 77 323 78 } //for (i=0;i<iomodel->numberofnodes;i++) … … 325 80 } //if (strcmp(iomodel->meshtype,"2d")==0) 326 81 327 328 82 /*Free data: */ 329 83 xfree((void**)&iomodel->elements); … … 337 91 xfree((void**)&iomodel->B); 338 92 xfree((void**)&iomodel->n); 339 340 93 341 /*Add one constant material property to materials: */ 342 matpar_mid=iomodel->numberofnodes-iomodel->numberofnodes2d+1; //put it at the end of the materials 343 matpar_g=iomodel->g; 344 matpar_rho_ice=iomodel->rho_ice; 345 matpar_rho_water=iomodel->rho_water; 346 matpar_thermalconductivity=iomodel->thermalconductivity; 347 matpar_heatcapacity=iomodel->heatcapacity; 348 matpar_latentheat=iomodel->latentheat; 349 matpar_beta=iomodel->beta; 350 matpar_meltingpoint=iomodel->meltingpoint; 351 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 352 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 353 354 /*Create matpar object using its constructor: */ 355 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, 356 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, 357 matpar_thermal_exchange_velocity,matpar_g); 358 359 /*Add to materials datset: */ 360 materials->AddObject(matpar); 361 362 /*Create nodes from x,y,z, as well as the spc values on those grids: */ 94 /*Add new constrant material property to materials, at the end: */ 95 materials->AddObject(new Matpar(iomodel)); 363 96 364 97 /*First fetch data: */ … … 377 110 IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); 378 111 IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); 379 380 /*Get number of dofs per node: */381 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);382 112 383 113 for (i=0;i<iomodel->numberofnodes;i++){ 384 #ifdef _PARALLEL_385 /*keep only this partition's nodes:*/386 if(my_grids[i]){387 #endif388 114 389 node_id=i+1; //matlab indexing390 115 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 116 if(iomodel->my_vertices[i]){ 391 117 392 #ifdef _PARALLEL_ 393 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border 394 node_partitionborder=1; 118 /*Add vertex to vertices dataset: */ 119 vertices->AddObject(new Vertex(i,iomodel)); 120 121 /*Add node to nodes dataset: */ 122 nodes->AddObject(new Node(i,iomodel)); 123 395 124 } 396 else{397 node_partitionborder=0;398 }399 #else400 node_partitionborder=0;401 #endif402 403 node_x[0]=iomodel->x[i];404 node_x[1]=iomodel->y[i];405 node_x[2]=iomodel->z[i];406 node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);407 408 node_onbed=(int)iomodel->gridonbed[i];409 node_onsurface=(int)iomodel->gridonsurface[i];410 node_onshelf=(int)iomodel->gridoniceshelf[i];411 node_onsheet=(int)iomodel->gridonicesheet[i];412 413 if (strcmp(iomodel->meshtype,"3d")==0){414 if (isnan(iomodel->uppernodes[i])){415 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.416 }417 else{418 node_upper_node_id=(int)iomodel->uppernodes[i];419 }420 }421 else{422 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/423 node_upper_node_id=node_id;424 }425 426 /*Create node using its constructor: */427 node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);428 429 /*set single point constraints.: */430 if (!iomodel->gridonhutter[i]){431 for(k=1;k<=node_numdofs;k++){432 node->FreezeDof(k);433 }434 }435 436 /*Add node to nodes dataset: */437 nodes->AddObject(node);438 439 #ifdef _PARALLEL_440 } //if(my_grids[i])441 #endif442 125 } 443 444 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these445 * datasets, it will not be redone: */446 elements->Presort();447 nodes->Presort();448 materials->Presort();449 126 450 127 /*Clean fetched data: */ … … 461 138 xfree((void**)&iomodel->gridonicesheet); 462 139 xfree((void**)&iomodel->gridoniceshelf); 463 464 140 465 /*Keep partitioning information into iomodel*/ 466 iomodel->epart=epart; 467 iomodel->my_grids=my_grids; 468 iomodel->my_bordergrids=my_bordergrids; 469 470 /*Keep partitioning information into iomodel*/ 471 #ifdef _PARALLEL_ 472 xfree((void**)&all_numgrids); 473 xfree((void**)&npart); 474 VecFree(&gridborder); 475 #endif 476 iomodel->npart=npart; 141 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 142 * datasets, it will not be redone: */ 143 elements->Presort(); 144 nodes->Presort(); 145 vertices->Presort(); 146 materials->Presort(); 477 147 478 148 cleanup_and_return: … … 481 151 *pelements=elements; 482 152 *pnodes=nodes; 153 *pvertices=vertices; 483 154 *pmaterials=materials; 484 155 -
issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
r3422 r3423 22 22 DataSet* vertices = NULL; 23 23 DataSet* materials = NULL; 24 25 /*Objects: */26 Node* node = NULL;27 Vertex* vertex = NULL;28 Penta* penta = NULL;29 Matice* matice = NULL;30 Matpar* matpar = NULL;31 24 32 25 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ -
issm/trunk/src/c/objects/Node.cpp
r3422 r3423 98 98 99 99 /*set single point constraints.: */ 100 /*FROM DIAGNOSTICHORIZ*/ 100 101 if (strcmp(iomodel->meshtype,"3d")==0){ 101 102 /*We have a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ … … 111 112 } 112 113 } 114 /*FROM DIAGNOSTICSTOKES*/ 113 115 /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ 114 116 if (iomodel->borderstokes[i]){ … … 121 123 for(k=1;k<=numdofs;k++){ 122 124 this->FreezeDof(k); 125 } 126 } 127 /*FROM DIAGNOSTICHUTTER*/ 128 if (!iomodel->gridonhutter[i]){ 129 for(k=1;k<=numdofs;k++){ 130 node->FreezeDof(k); 123 131 } 124 132 }
Note:
See TracChangeset
for help on using the changeset viewer.
