Changeset 3431


Ignore:
Timestamp:
04/07/10 16:16:44 (15 years ago)
Author:
seroussi
Message:

added ModelProcessor Balancedvelocities

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TabularUnified issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp

    r3417 r3431  
    1212#include "../IoModel.h"
    1313
    14 
    1514void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1615
    17 
    18         /*output: int* epart, int* my_grids, double* my_bordergrids*/
    19 
    20 
     16        /*Intermediary*/
    2117        int i,j,k,n;
    22         extern int my_rank;
    23         extern int num_procs;
    2418
    2519        /*DataSets: */
     
    2822        DataSet*    vertices = NULL;
    2923        DataSet*    materials = NULL;
    30         ElementProperties* tria_properties=NULL;
    31         ElementProperties* penta_properties=NULL;
    32        
    33         /*Objects: */
    34         Node*       node   = NULL;
    35         Vertex*     vertex = NULL;
    36         Tria*       tria = NULL;
    37         Penta*      penta = NULL;
    38         Matice*     matice  = NULL;
    39         Matpar*     matpar  = NULL;
    40 
    41         /*output: */
    42         int* epart=NULL; //element partitioning.
    43         int* npart=NULL; //node partitioning.
    44         int* my_grids=NULL;
    45         double* my_bordergrids=NULL;
    46 
    47 
    48         /*intermediary: */
    49         int elements_width; //size of elements
    50         double B_avg;
    51                        
    52         /*tria constructor input: */
    53         int tria_id;
    54         int tria_matice_id;
    55         int tria_matpar_id;
    56         int tria_numpar_id;
    57         int tria_node_ids[3];
    58         double tria_h[3];
    59         double tria_s[3];
    60         double tria_b[3];
    61         int    tria_shelf;
    62         bool   tria_onwater;
    63 
    64         /*matice constructor input: */
    65         int    matice_mid;
    66         double matice_B;
    67         double matice_n;
    68        
    69         /*penta constructor input: */
    70 
    71         int penta_id;
    72         int penta_matice_id;
    73         int penta_matpar_id;
    74         int penta_numpar_id;
    75         int penta_node_ids[6];
    76         double penta_h[6];
    77         double penta_s[6];
    78         double penta_b[6];
    79         int penta_shelf;
    80         int penta_onbed;
    81         int penta_onsurface;
    82         int penta_collapse;
    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 partitionborder=0;
    102         int node_onbed;
    103         int node_onsurface;
    104         int node_onshelf;
    105         int node_onsheet;
    106         int node_upper_node_id;
    107         int node_numdofs;
    108 
    109                
    110         #ifdef _PARALLEL_
    111         /*Metis partitioning: */
    112         int  range;
    113         Vec  gridborder=NULL;
    114         int  my_numgrids;
    115         int* all_numgrids=NULL;
    116         int  gridcount;
    117         int  count;
    118         #endif
    119         int  first_grid_index;
    12024
    12125        /*First create the elements, nodes and material properties: */
     
    12529        materials = new DataSet(MaterialsEnum());
    12630
    127         /*Width of elements: */
    128         if(strcmp(iomodel->meshtype,"2d")==0){
    129                 elements_width=3; //tria elements
    130         }
    131         else{
    132                 elements_width=6; //penta elements
    133         }
    134 
    135         #ifdef _PARALLEL_
    136         /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
    137         if(strcmp(iomodel->meshtype,"2d")==0){
    138                 /*load elements: */
    139                 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    140         }
    141         else{
    142                 /*load elements2d: */
    143                 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
    144         }
    145 
    146 
    147         MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
    148 
    149         /*Free elements and elements2d: */
    150         xfree((void**)&iomodel->elements);
    151         xfree((void**)&iomodel->elements2d);
    152 
    153         /*Used later on: */
    154         my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
    155         #endif
    156 
    157 
    158 
    159         /*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);
    16033
    16134        /*2d mesh: */
     
    17245                for (i=0;i<iomodel->numberofelements;i++){
    17346
    174                 #ifdef _PARALLEL_
    175                 /*!All elements have been partitioned above, only create elements for this CPU: */
    176                 if(my_rank==epart[i]){
    177                 #endif
    178                        
    179                        
    180                         /*ids: */
    181                         tria_id=i+1; //matlab indexing.
    182                         tria_matice_id=-1; //no need for materials
    183                         tria_matpar_id=-1; //no need for materials
    184                         tria_numpar_id=1;
     47                        if(iomodel->my_elements[i]){
    18548
    186                         /*vertices offsets: */
    187                         tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0);
    188                         tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1);
    189                         tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2);
     49                                /*Create and add tria element to elements dataset: */
     50                                elements->AddObject(new Tria(i,iomodel));
    19051
    191                         /*thickness,surface and bed:*/
    192                         tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
    193                         tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1));
    194                         tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
    195 
    196                         tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1));
    197                         tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1));
    198                         tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1));
    199 
    200                         tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1));
    201                         tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1));
    202                         tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1));
    203 
    204                         /*element on iceshelf?:*/
    205                         tria_shelf=(int)*(iomodel->elementoniceshelf+i);
    206                         tria_onwater=(bool)*(iomodel->elementonwater+i);
    207 
    208                         /*Create properties: */
    209                         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);
    210 
    211                         /*Create tria element using its constructor:*/
    212                         tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties);
    213 
    214                         /*Delete properties: */
    215                         delete tria_properties;
    216 
    217                         /*Add tria element to elements dataset: */
    218                         elements->AddObject(tria);
    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 use
    222                          *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_grids
    224                          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                         #endif
    229 
    230                 #ifdef _PARALLEL_
    231                 }//if(my_rank==epart[i])
    232                 #endif
     52                                /*Create and add material property to materials dataset: */
     53                                materials->AddObject(new Matice(i,iomodel,3));
     54                        }
    23355
    23456                }//for (i=0;i<numberofelements;i++)
     
    25577                IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
    25678                IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
    257        
     79
    25880                for (i=0;i<iomodel->numberofelements;i++){
    259                 #ifdef _PARALLEL_
    260                 /*We are using our element partition to decide which elements will be created on this node: */
    261                 if(my_rank==epart[i]){
    262                 #endif
     81                        if(iomodel->my_elements[i]){
     82                                /*Create and add penta element to elements dataset: */
     83                                elements->AddObject(new Penta(i,iomodel));
    26384
    264                        
    265                         /*name and id: */
    266                         penta_id=i+1; //matlab indexing.
    267                         penta_matice_id=-1;
    268                         penta_matpar_id=-1; //no need for materials
    269                         penta_numpar_id=1;
     85                                /*Create and add material property to materials dataset: */
     86                                materials->AddObject(new Matice(i,iomodel,6));
    27087
    271                         /*vertices,thickness,surface,bed and drag: */
    272                         for(j=0;j<6;j++){
    273                                 penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);
    274                                 penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    275                                 penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    276                                 penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    277                         }
    278 
    279                         /*diverse: */
    280                         penta_shelf=(int)*(iomodel->elementoniceshelf+i);
    281                         penta_onbed=(int)*(iomodel->elementonbed+i);
    282                         penta_onsurface=(int)*(iomodel->elementonsurface+i);
    283                         penta_collapse=1;
    284                         penta_onwater=(bool)*(iomodel->elementonwater+i);
    285        
    286                         /*Create properties: */
    287                         penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
    288 
    289                         /*Create Penta using its constructor:*/
    290                         penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);
    291 
    292                         /*Delete properties: */
    293                         delete penta_properties;
    294 
    295                         /*Add penta element to elements dataset: */
    296                         elements->AddObject(penta);
    297        
    298                         #ifdef _PARALLEL_
    299                         /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    300                          *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    301                          into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    302                          will hold which grids belong to this partition*/
    303                         my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
    304                         my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
    305                         my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
    306                         my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
    307                         my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
    308                         my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
    309                         #endif
    310 
    311                 #ifdef _PARALLEL_
    312                 }//if(my_rank==epart[i])
    313                 #endif
     88                        }//if(my_elements[i])
    31489
    31590                }//for (i=0;i<numberofelements;i++)
     
    327102        } //if (strcmp(meshtype,"2d")==0)
    328103
    329         #ifdef _PARALLEL_
    330                 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
    331                  *element partition, and which are on its border with other nodes:*/
    332                 gridborder=NewVec(iomodel->numberofnodes);
    333 
    334                 for (i=0;i<iomodel->numberofnodes;i++){
    335                         if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
    336                 }
    337                 VecAssemblyBegin(gridborder);
    338                 VecAssemblyEnd(gridborder);
    339 
    340                 #ifdef _ISSM_DEBUG_
    341                 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
    342                 #endif
    343                
    344                 VecToMPISerial(&my_bordergrids,gridborder);
    345 
    346                 #ifdef _ISSM_DEBUG_
    347                 if(my_rank==0){
    348                         for (i=0;i<iomodel->numberofnodes;i++){
    349                                 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
    350                         }
    351                 }
    352                 #endif
    353         #endif
    354 
    355         /*Add one constant material property to materials: */
    356         matpar_mid=1; //put it at the end of the materials
    357         matpar_g=iomodel->g;
    358         matpar_rho_ice=iomodel->rho_ice;
    359         matpar_rho_water=iomodel->rho_water;
    360         matpar_thermalconductivity=iomodel->thermalconductivity;
    361         matpar_heatcapacity=iomodel->heatcapacity;
    362         matpar_latentheat=iomodel->latentheat;
    363         matpar_beta=iomodel->beta;
    364         matpar_meltingpoint=iomodel->meltingpoint;
    365         matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
    366         matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
    367 
    368         /*Create matpar object using its constructor: */
    369         matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
    370                         matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
    371                         matpar_thermal_exchange_velocity,matpar_g);
    372                
    373         /*Add to materials datset: */
    374         materials->AddObject(matpar);
    375 
    376         /*Partition penalties in 3d: */
    377         if(strcmp(iomodel->meshtype,"3d")==0){
    378        
    379                 /*Get penalties: */
    380                 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
    381 
    382                 if(iomodel->numpenalties){
    383 
    384                         iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
    385                         #ifdef _SERIAL_
    386                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
    387                         #else
    388                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
    389 
    390                         for(i=0;i<iomodel->numpenalties;i++){
    391                                 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
    392                                 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
    393                                         /*All grids that are being penalised belong to this node's internal grid partition.:*/
    394                                         iomodel->penaltypartitioning[i]=1;
    395                                 }
    396                                 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
    397                                         iomodel->penaltypartitioning[i]=0;
    398                                 }
    399                         }
    400                         #endif
    401                 }
    402 
    403                 /*Free penalties: */
    404                 xfree((void**)&iomodel->penalties);
    405         }
    406 
    407         /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
    408          We can therefore determine  which grids are internal to this node's partition
    409          and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
    410          that, go and create the grids*/
    411 
    412         /*Create nodes from x,y,z, as well as the spc values on those grids: */
     104        /*Add new constrant material property tgo materials, at the end: */
     105        materials->AddObject(new Matpar(iomodel));
    413106               
    414107        /*First fetch data: */
     
    427120        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
    428121
     122        for (i=0;i<iomodel->numberofvertices;i++){
    429123
    430         /*Get number of dofs per node: */
    431         DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
     124                /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
     125                if(iomodel->my_vertices[i]){
     126                       
     127                        /*Add vertex to vertices dataset: */
     128                        vertices->AddObject(new Vertex(i,iomodel));
    432129
    433         for (i=0;i<iomodel->numberofnodes;i++){
    434         #ifdef _PARALLEL_
    435         /*keep only this partition's nodes:*/
    436         if((my_grids[i]==1)){
    437         #endif
     130                        /*Add node to nodes dataset: */
     131                        nodes->AddObject(new Node(i,iomodel));
    438132
    439                 /*create vertex: */
    440                 vertex_id=i+1;
    441                 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
    442 
    443                 vertices->AddObject(vertex);
    444 
    445                 node_id=i+1; //matlab indexing
    446                        
    447                 #ifdef _PARALLEL_
    448                 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
    449                         node_partitionborder=1;
    450133                }
    451                 else{
    452                         node_partitionborder=0;
    453                 }
    454                 #else
    455                         node_partitionborder=0;
    456                 #endif
    457 
    458                 node_properties.SetProperties(
    459                         (int)iomodel->gridonbed[i],
    460                         (int)iomodel->gridonsurface[i],
    461                         (int)iomodel->gridoniceshelf[i],
    462                         (int)iomodel->gridonicesheet[i]);
    463 
    464 
    465                 if (strcmp(iomodel->meshtype,"3d")==0){
    466                         if (isnan(iomodel->uppernodes[i])){
    467                                 node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
    468                         }
    469                         else{
    470                                 node_upper_node_id=(int)iomodel->uppernodes[i];
    471                         }
    472                 }
    473                 else{
    474                         /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
    475                         node_upper_node_id=node_id;
    476                 }
    477 
    478                 /*Create node using its constructor: */
    479                 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
    480 
    481                 /*set single point constraints.: */
    482                 if (strcmp(iomodel->meshtype,"3d")==0){
    483                         /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
    484                         if (!iomodel->gridonbed[i]){
    485                                 for(k=1;k<=node_numdofs;k++){
    486                                         node->FreezeDof(k);
    487                                 }
    488                         }
    489                 }
    490                 /*Add node to nodes dataset: */
    491                 nodes->AddObject(node);
    492 
    493         #ifdef _PARALLEL_
    494         } //if((my_grids[i]==1))
    495         #endif
    496134        }
    497 
    498         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    499          * datasets, it will not be redone: */
    500         elements->Presort();
    501         nodes->Presort();
    502         vertices->Presort();
    503         materials->Presort();
    504135
    505136        /*Clean fetched data: */
     
    516147        xfree((void**)&iomodel->gridoniceshelf);
    517148       
    518 
    519         /*Keep partitioning information into iomodel*/
    520         iomodel->epart=epart;
    521         iomodel->my_grids=my_grids;
    522         iomodel->my_bordergrids=my_bordergrids;
    523 
    524         /*Free ressources:*/
    525         #ifdef _PARALLEL_
    526         xfree((void**)&all_numgrids);
    527         xfree((void**)&npart);
    528         VecFree(&gridborder);
    529         #endif
     149        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     150         * datasets, it will not be redone: */
     151        elements->Presort();
     152        nodes->Presort();
     153        vertices->Presort();
     154        materials->Presort();
    530155
    531156        cleanup_and_return:
Note: See TracChangeset for help on using the changeset viewer.