Changeset 3426


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

added ModelProcessor for melting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp

    r3417 r3426  
    1212#include "../IoModel.h"
    1313
    14 
    1514void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1615
    17 
    1816        /*output: int* epart, int* my_grids, double* my_bordergrids*/
    19 
    20 
    2117        int i,j,k,n;
    22         extern int my_rank;
    23         extern int num_procs;
    2418
    2519        /*DataSets: */
     
    2923        DataSet*    materials = NULL;
    3024       
    31         /*Objects: */
    32         Node*       node   = NULL;
    33         Vertex*     vertex = NULL;
    34         Penta*      penta = NULL;
    35         Matice*     matice  = NULL;
    36         Matpar*     matpar  = NULL;
    37         ElementProperties* penta_properties=NULL;
    38        
    39         /*output: */
    40         int* epart=NULL; //element partitioning.
    41         int* npart=NULL; //node partitioning.
    42         int* my_grids=NULL;
    43         double* my_bordergrids=NULL;
    44 
    45 
    46         /*intermediary: */
    47         int elements_width; //size of elements
    48         double B_avg;
    49 
    50         /*matice constructor input: */
    51         int    matice_mid;
    52         double matice_B;
    53         double matice_n;
    54        
    55         /*penta constructor input: */
    56 
    57         int penta_id;
    58         int penta_matice_id;
    59         int penta_matpar_id;
    60         int penta_numpar_id;
    61         int penta_node_ids[6];
    62         double penta_h[6];
    63         double penta_s[6];
    64         double penta_b[6];
    65         double penta_k[6];
    66         int penta_friction_type;
    67         double penta_p;
    68         double penta_q;
    69         int penta_shelf;
    70         int penta_onbed;
    71         int penta_onsurface;
    72         int penta_collapse;
    73         double penta_melting[6];
    74         double penta_accumulation[6];
    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 
    11325        /*First create the elements, nodes and material properties: */
    11426        elements  = new DataSet(ElementsEnum());
     
    11729        materials = new DataSet(MaterialsEnum());
    11830
    119         /*Width of elements: */
    120         if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
    121         elements_width=6; //penta elements
    122 
    123 
    124         #ifdef _PARALLEL_
    125         /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
    126         IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
    127 
    128 
    129         MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
    130 
    131         /*Free elements and elements2d: */
    132         xfree((void**)&iomodel->elements2d);
    133 
    134         /*Used later on: */
    135         my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
    136         #endif
    137 
    138 
    139 
    140         /*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);
    14133
    14234        /*Fetch data needed: */
     
    15951       
    16052        for (i=0;i<iomodel->numberofelements;i++){
    161         #ifdef _PARALLEL_
    162         /*We are using our element partition to decide which elements will be created on this node: */
    163         if(my_rank==epart[i]){
    164         #endif
     53               
     54                        if(iomodel->my_elements[i]){
    16555
    166                
    167                 /*name and id: */
    168                 penta_id=i+1; //matlab indexing.
    169                 penta_matice_id=i+1; //refers to the corresponding material property card
    170                 penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card
    171                 penta_numpar_id=1;
     56                                /*Create and add tria element to elements dataset: */
     57                                elements->AddObject(new Tria(i,iomodel));
    17258
    173                 /*vertices,thickness,surface,bed and drag: */
    174                 for(j=0;j<6;j++){
    175                         penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);
    176                         penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    177                         penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    178                         penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    179                         penta_k[j]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+j)-1));
    180                         penta_melting[j]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+j)-1));
    181                         penta_accumulation[j]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+j)-1));
    182                 }
    183 
    184                 /*basal drag:*/
    185                 penta_friction_type=(int)iomodel->drag_type;
    186 
    187                 penta_p=iomodel->p[i];
    188                 penta_q=iomodel->q[i];
    189 
    190                 /*diverse: */
    191                 penta_shelf=(int)*(iomodel->elementoniceshelf+i);
    192                 penta_onbed=(int)*(iomodel->elementonbed+i);
    193                 penta_onsurface=(int)*(iomodel->elementonsurface+i);
    194                 penta_onwater=(bool)*(iomodel->elementonwater+i);
    195                 penta_collapse=1;
    196 
    197                 /*Create properties: */
    198                 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
    199 
    200                 /*Create Penta using its constructor:*/
    201                 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);
    202 
    203                 /*Delete properties: */
    204                 delete penta_properties;
    205 
    206                 /*Add penta element to elements dataset: */
    207                 elements->AddObject(penta);
    208 
    209 
    210                 /*Deal with material:*/
    211                 matice_mid=i+1; //same as the material id from the geom2 elements.
    212                 /*Average B over 6 element grids: */
    213                 B_avg=0;
    214                 for(j=0;j<6;j++){
    215                         B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
    216                 }
    217                 B_avg=B_avg/6;
    218                 matice_B= B_avg;
    219                 matice_n=(double)*(iomodel->n+i);
    220 
    221                 /*Create matice using its constructor:*/
    222                 matice= new Matice(matice_mid,matice_B,matice_n);
    223 
    224                 /*Add matice element to materials dataset: */
    225                 materials->AddObject(matice);
    226                
    227                 #ifdef _PARALLEL_
    228                 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    229                  *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    230                  into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    231                  will hold which grids belong to this partition*/
    232                 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
    233                 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
    234                 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
    235                 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
    236                 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
    237                 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
    238                 #endif
    239 
    240         #ifdef _PARALLEL_
    241         }//if(my_rank==epart[i])
    242         #endif
     59                                /*Create and add material property to materials dataset: */
     60                                materials->AddObject(new Matice(i,iomodel,3));
     61                        }
    24362
    24463        }//for (i=0;i<numberofelements;i++)
     
    26281        xfree((void**)&iomodel->elementonwater);
    26382
     83        /*Add new constrant material property tgo materials, at the end: */
     84        materials->AddObject(new Matpar(iomodel));
    26485
    265         /*Add one constant material property to materials: */
    266         matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
    267         matpar_g=iomodel->g;
    268         matpar_rho_ice=iomodel->rho_ice;
    269         matpar_rho_water=iomodel->rho_water;
    270         matpar_thermalconductivity=iomodel->thermalconductivity;
    271         matpar_heatcapacity=iomodel->heatcapacity;
    272         matpar_latentheat=iomodel->latentheat;
    273         matpar_beta=iomodel->beta;
    274         matpar_meltingpoint=iomodel->meltingpoint;
    275         matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
    276         matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
    277 
    278         /*Create matpar object using its constructor: */
    279         matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
    280                         matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
    281                         matpar_thermal_exchange_velocity,matpar_g);
    282                
    283         /*Add to materials datset: */
    284         materials->AddObject(matpar);
    285        
    286         #ifdef _PARALLEL_
    287                 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
    288                  *element partition, and which are on its border with other nodes:*/
    289                 gridborder=NewVec(iomodel->numberofnodes);
    290 
    291                 for (i=0;i<iomodel->numberofnodes;i++){
    292                         if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
    293                 }
    294                 VecAssemblyBegin(gridborder);
    295                 VecAssemblyEnd(gridborder);
    296 
    297                 #ifdef _ISSM_DEBUG_
    298                 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
    299                 #endif
    300                
    301                 VecToMPISerial(&my_bordergrids,gridborder);
    302 
    303                 #ifdef _ISSM_DEBUG_
    304                 if(my_rank==0){
    305                         for (i=0;i<iomodel->numberofnodes;i++){
    306                                 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
    307                         }
    308                 }
    309                 #endif
    310         #endif
    311 
    312         /*Partition penalties in 3d: */
    313         if(strcmp(iomodel->meshtype,"3d")==0){
    314        
    315                 /*Get penalties: */
    316                 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
    317 
    318                 if(iomodel->numpenalties){
    319 
    320                         iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
    321                         #ifdef _SERIAL_
    322                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
    323                         #else
    324                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
    325 
    326                         for(i=0;i<iomodel->numpenalties;i++){
    327                                 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
    328                                 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
    329                                         /*All grids that are being penalised belong to this node's internal grid partition.:*/
    330                                         iomodel->penaltypartitioning[i]=1;
    331                                 }
    332                                 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
    333                                         iomodel->penaltypartitioning[i]=0;
    334                                 }
    335                         }
    336                         #endif
    337                 }
    338 
    339                 /*Free penalties: */
    340                 xfree((void**)&iomodel->penalties);
    341         }
    342 
    343         /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
    344          We can therefore determine  which grids are internal to this node's partition
    345          and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
    346          that, go and create the grids*/
    347 
    348         /*Create nodes from x,y,z, as well as the spc values on those grids: */
    349                
    35086        /*First fetch data: */
    35187        if (strcmp(iomodel->meshtype,"3d")==0){
     
    36399        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
    364100
     101        for (i=0;i<iomodel->numberofvertices;i++){
    365102
    366         /*Get number of dofs per node: */
    367         DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
     103                /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
     104                if(iomodel->my_vertices[i]){
     105                       
     106                        /*Add vertex to vertices dataset: */
     107                        vertices->AddObject(new Vertex(i,iomodel));
    368108
    369         for (i=0;i<iomodel->numberofnodes;i++){
    370         #ifdef _PARALLEL_
    371         /*keep only this partition's nodes:*/
    372         if((my_grids[i]==1)){
    373         #endif
     109                        /*Add node to nodes dataset: */
     110                        nodes->AddObject(new Node(i,iomodel));
    374111
    375                 /*create vertex: */
    376                 vertex_id=i+1;
    377                 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
    378                 vertices->AddObject(vertex);
    379 
    380 
    381                 node_id=i+1; //matlab indexing
    382                
    383                 #ifdef _PARALLEL_
    384                 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
    385                         node_partitionborder=1;
    386112                }
    387                 else{
    388                         node_partitionborder=0;
    389                 }
    390                 #else
    391                         node_partitionborder=0;
    392                 #endif
    393 
    394                 node_properties.SetProperties(
    395                         (int)iomodel->gridonbed[i],
    396                         (int)iomodel->gridonsurface[i],
    397                         (int)iomodel->gridoniceshelf[i],
    398                         (int)iomodel->gridonicesheet[i]);
    399 
    400                 if (strcmp(iomodel->meshtype,"3d")==0){
    401                         if (isnan(iomodel->uppernodes[i])){
    402                                 node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
    403                         }
    404                         else{
    405                                 node_upper_node_id=(int)iomodel->uppernodes[i];
    406                         }
    407                 }
    408                 else{
    409                         /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
    410                         node_upper_node_id=node_id;
    411                 }
    412 
    413                 /*Create node using its constructor: */
    414                 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
    415 
    416                 /*set single point constraints.: */
    417                 if (!iomodel->gridonbed[i]){
    418                         for(k=1;k<=node_numdofs;k++){
    419                                 node->FreezeDof(k);
    420                         }
    421                 }
    422                 /*Add node to nodes dataset: */
    423                 nodes->AddObject(node);
    424 
    425         #ifdef _PARALLEL_
    426         } //if((my_grids[i]==1))
    427         #endif
    428113        }
    429 
    430         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    431          * datasets, it will not be redone: */
    432         elements->Presort();
    433         nodes->Presort();
    434         vertices->Presort();
    435         materials->Presort();
    436114
    437115        /*Clean fetched data: */
     
    448126        xfree((void**)&iomodel->gridoniceshelf);
    449127
    450         /*Keep partitioning information into iomodel*/
    451         iomodel->epart=epart;
    452         iomodel->my_grids=my_grids;
    453         iomodel->my_bordergrids=my_bordergrids;
    454 
    455         /*Free ressources:*/
    456         #ifdef _PARALLEL_
    457         xfree((void**)&all_numgrids);
    458         xfree((void**)&npart);
    459         VecFree(&gridborder);
    460         #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();
    461134
    462135        cleanup_and_return:
Note: See TracChangeset for help on using the changeset viewer.