Changeset 3427


Ignore:
Timestamp:
04/07/10 16:07:57 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added prognostic2

Location:
issm/trunk/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp

    r3417 r3427  
    1515void    CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1616
     17        /*Intermediary*/
    1718        int i,j,k,n;
    18         extern int my_rank;
    19         extern int num_procs;
    2019
    2120        /*DataSets: */
    2221        DataSet* elements  = NULL;
    2322        DataSet* nodes = NULL;
    24         DataSet*    vertices = NULL;
     23        DataSet* vertices = NULL;
    2524        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 elements
    44         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         #endif
    124         int  first_grid_index;
    12525
    12626        /*First create the elements, nodes and material properties: */
    … …  
    13030        materials = new DataSet(MaterialsEnum());
    13131
    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);
    16234
    16335        /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
    164 
    16536        /*2d mesh: */
    16637        if (strcmp(iomodel->meshtype,"2d")==0){
    … …  
    17647                for (i=0;i<iomodel->numberofelements;i++){
    17748
    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]){
    18850
    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));
    19353
    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                        }
    23757                }//for (i=0;i<numberofelements;i++)
    23858       
    … …  
    24969        } //if (strcmp(meshtype,"2d")==0)
    25070
    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));
    25573
    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: */
    34175        IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
    34276        IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
    … …  
    34680        IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
    34781        IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
     82        IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
    34883        IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
    34984        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");
    37088        }
    37189
    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)*/
    373103        for (i=0;i<iomodel->numberofelements;i++){
    374104                for (j=0;j<3;j++){
    375105
    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]){
    380107
    381108                                //Get id of the vertex on which the current node is located
    382109                                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);
    419111
    420112                                /*Add node to nodes dataset: */
    421                                 nodes->AddObject(node);
    422                         #ifdef _PARALLEL_
     113                                nodes->AddObject(new Node(k,i,iomodel));
     114
    423115                        }
    424                         #endif
    425116                }
    426117        }
     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);
    427132
    428133        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    … …  
    432137        vertices->Presort();
    433138        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         #endif
    460139
    461140        cleanup_and_return:

  • issm/trunk/src/c/objects/Node.cpp

    r3423 r3427  
    5252}
    5353/*}}}*/
    54 /*FUNCTION Node constructor  from iomodel{{{2*/
     54/*FUNCTION Node constructor from iomodel continuous Galerkin{{{2*/
    5555Node::Node(int i, IoModel* iomodel){ //i is the node index
    5656
    … …  
    131131                }
    132132        }
     133}
     134/*}}}*/
     135/*FUNCTION Node constructor from iomodel discontinuous Galerkin{{{2*/
     136Node::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
    133180}
    134181/*}}}*/

  • issm/trunk/src/c/objects/Node.h

    r3420 r3427  
    4444                Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node);
    4545                Node(int i, IoModel* iomodel);
     46                Node(int i,int j,IoModel* iomodel);
    4647                ~Node();
    4748                /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.