Changeset 3428


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

added ModelProcessor for Thermal

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp

    r3417 r3428  
    1414void    CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1515
    16 
    1716        /*output: int* epart, int* my_grids, double* my_bordergrids*/
    18 
    19 
    2017        int i,j,k,n;
    21         extern int my_rank;
    22         extern int num_procs;
    2318
    2419        /*DataSets: */
     
    2823        DataSet*    materials = NULL;
    2924       
    30         /*Objects: */
    31         Node*       node   = NULL;
    32         Vertex*     vertex = NULL;
    33         Penta*      penta = NULL;
    34         Matice*     matice  = NULL;
    35         Matpar*     matpar  = NULL;
    36         ElementProperties* penta_properties=NULL;
    37 
    38         /*output: */
    39         int* epart=NULL; //element partitioning.
    40         int* npart=NULL; //node partitioning.
    41         int* my_grids=NULL;
    42         double* my_bordergrids=NULL;
    43 
    44 
    45         /*intermediary: */
    46         int elements_width; //size of elements
    47         double B_avg;
    48                        
    49         /*matice constructor input: */
    50         int    matice_mid;
    51         double matice_B;
    52         double matice_n;
    53        
    54         /*penta constructor input: */
    55 
    56         int penta_id;
    57         int penta_matice_id;
    58         int penta_matpar_id;
    59         int penta_numpar_id;
    60         int penta_node_ids[6];
    61         double penta_h[6];
    62         double penta_s[6];
    63         double penta_b[6];
    64         double penta_k[6];
    65         int penta_friction_type;
    66         double penta_p;
    67         double penta_q;
    68         int penta_shelf;
    69         int penta_onbed;
    70         int penta_onsurface;
    71         int penta_collapse;
    72         double penta_geothermalflux[6];
    73         bool   penta_onwater;
    74 
    75         /*matpar constructor input: */
    76         int     matpar_mid;
    77         double matpar_rho_ice;
    78         double matpar_rho_water;
    79         double matpar_heatcapacity;
    80         double matpar_thermalconductivity;
    81         double matpar_latentheat;
    82         double matpar_beta;
    83         double matpar_meltingpoint;
    84         double matpar_mixed_layer_capacity;
    85         double matpar_thermal_exchange_velocity;
    86         double matpar_g;
    87 
    88         /* node constructor input: */
    89         int node_id;
    90         int vertex_id;
    91         int node_partitionborder=0;
    92         int node_onbed;
    93         int node_onsurface;
    94         int node_onshelf;
    95         int node_onsheet;
    96         int node_upper_node_id;
    97         int node_numdofs;
    98 
    99                
    100         #ifdef _PARALLEL_
    101         /*Metis partitioning: */
    102         int  range;
    103         Vec  gridborder=NULL;
    104         int  my_numgrids;
    105         int* all_numgrids=NULL;
    106         int  gridcount;
    107         int  count;
    108         #endif
    109         int  first_grid_index;
    110 
    11125        /*First create the elements, nodes and material properties: */
    11226        elements  = new DataSet(ElementsEnum());
     
    11529        materials = new DataSet(MaterialsEnum());
    11630
    117         /*return if 2d mesh*/
    118         if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
    119 
    120         /*Width of elements: */
    121         elements_width=6; //penta elements
    122 
    123         #ifdef _PARALLEL_
    124         /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
    125         IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
    126 
    127         MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
    128 
    129         /*Free elements and elements2d: */
    130         xfree((void**)&iomodel->elements2d);
    131 
    132         /*Used later on: */
    133         my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
    134         #endif
    135 
    136 
    137         /*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);
    13833
    13934        /*Fetch data needed: */
     
    15550       
    15651        for (i=0;i<iomodel->numberofelements;i++){
    157         #ifdef _PARALLEL_
    158         /*We are using our element partition to decide which elements will be created on this node: */
    159         if(my_rank==epart[i]){
    160         #endif
     52                if(iomodel->my_elements[i]){
    16153
    162                
    163                 /*name and id: */
    164                 penta_id=i+1; //matlab indexing.
    165                 penta_matice_id=i+1; //refers to the corresponding material property card
    166                 penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card
    167                 penta_numpar_id=1;
     54                        /*Create and add tria element to elements dataset: */
     55                        elements->AddObject(new Tria(i,iomodel));
    16856
    169                 /*vertices,thickness,surface,bed and drag: */
    170                 for(j=0;j<6;j++){
    171                         penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);
    172                         penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    173                         penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    174                         penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1));
    175                         penta_k[j]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+j)-1));
    176                         penta_geothermalflux[j]=*(iomodel->geothermalflux+        ((int)*(iomodel->elements+elements_width*i+j)-1));
     57                        /*Create and add material property to materials dataset: */
     58                        materials->AddObject(new Matice(i,iomodel,6));
    17759                }
    178 
    179                 /*basal drag:*/
    180                 penta_friction_type=(int)iomodel->drag_type;
    181 
    182                 penta_p=iomodel->p[i];
    183                 penta_q=iomodel->q[i];
    184 
    185                 /*diverse: */
    186                 penta_shelf=(int)*(iomodel->elementoniceshelf+i);
    187                 penta_onbed=(int)*(iomodel->elementonbed+i);
    188                 penta_onsurface=(int)*(iomodel->elementonsurface+i);
    189                 penta_onwater=(bool)*(iomodel->elementonwater+i);
    190 
    191                 /*We need the field collapse for transient, so that we can use compute B with the average temperature*/
    192                 if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
    193                         penta_collapse=1;
    194                 }
    195                 else{
    196                         penta_collapse=0;
    197                 }
    198 
    199                 /*Create properties: */
    200                 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, NULL, NULL, penta_geothermalflux, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
    201 
    202                 /*Create Penta using its constructor:*/
    203                 penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);
    204 
    205                 /*Delete properties: */
    206                 delete penta_properties;
    207 
    208                 /*Add penta element to elements dataset: */
    209                 elements->AddObject(penta);
    210 
    211                 /*Deal with material:*/
    212                 matice_mid=i+1; //same as the material id from the geom2 elements.
    213                 /*Average B over 6 element grids: */
    214                 B_avg=0;
    215                 for(j=0;j<6;j++){
    216                         B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
    217                 }
    218                 B_avg=B_avg/6;
    219                 matice_B= B_avg;
    220                 matice_n=(double)*(iomodel->n+i);
    221 
    222                 /*Create matice using its constructor:*/
    223                 matice= new Matice(matice_mid,matice_B,matice_n);
    224 
    225                 /*Add matice element to materials dataset: */
    226                 materials->AddObject(matice);
    227                
    228                 #ifdef _PARALLEL_
    229                 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    230                  *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    231                  into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    232                  will hold which grids belong to this partition*/
    233                 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
    234                 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
    235                 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
    236                 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
    237                 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
    238                 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
    239                 #endif
    240 
    241         #ifdef _PARALLEL_
    242         }//if(my_rank==epart[i])
    243         #endif
    244 
    24560        }//for (i=0;i<numberofelements;i++)
    24661
     
    26277        xfree((void**)&iomodel->elementonwater);
    26378
    264         /*Add one constant material property to materials: */
    265         matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
    266         matpar_g=iomodel->g;
    267         matpar_rho_ice=iomodel->rho_ice;
    268         matpar_rho_water=iomodel->rho_water;
    269         matpar_thermalconductivity=iomodel->thermalconductivity;
    270         matpar_heatcapacity=iomodel->heatcapacity;
    271         matpar_latentheat=iomodel->latentheat;
    272         matpar_beta=iomodel->beta;
    273         matpar_meltingpoint=iomodel->meltingpoint;
    274         matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
    275         matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
    276 
    277         /*Create matpar object using its constructor: */
    278         matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
    279                         matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
    280                         matpar_thermal_exchange_velocity,matpar_g);
    281                
    282         /*Add to materials datset: */
    283         materials->AddObject(matpar);
     79        /*Add new constrant material property tgo materials, at the end: */
     80        materials->AddObject(new Matpar(iomodel));
    28481       
    285         #ifdef _PARALLEL_
    286                 /*From the element partitioning, we can determine which grids are on the inside of this cpu's
    287                  *element partition, and which are on its border with other nodes:*/
    288                 gridborder=NewVec(iomodel->numberofnodes);
    289 
    290                 for (i=0;i<iomodel->numberofnodes;i++){
    291                         if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
    292                 }
    293                 VecAssemblyBegin(gridborder);
    294                 VecAssemblyEnd(gridborder);
    295 
    296                 #ifdef _ISSM_DEBUG_
    297                 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
    298                 #endif
    299                
    300                 VecToMPISerial(&my_bordergrids,gridborder);
    301 
    302                 #ifdef _ISSM_DEBUG_
    303                 if(my_rank==0){
    304                         for (i=0;i<iomodel->numberofnodes;i++){
    305                                 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
    306                         }
    307                 }
    308                 #endif
    309         #endif
    310 
    311         /*Partition penalties in 3d: */
    312         if(strcmp(iomodel->meshtype,"3d")==0){
    313        
    314                 /*Get penalties: */
    315                 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
    316 
    317                 if(iomodel->numpenalties){
    318 
    319                         iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
    320                         #ifdef _SERIAL_
    321                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
    322                         #else
    323                         for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
    324 
    325                         for(i=0;i<iomodel->numpenalties;i++){
    326                                 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
    327                                 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
    328                                         /*All grids that are being penalised belong to this node's internal grid partition.:*/
    329                                         iomodel->penaltypartitioning[i]=1;
    330                                 }
    331                                 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
    332                                         iomodel->penaltypartitioning[i]=0;
    333                                 }
    334                         }
    335                         #endif
    336                 }
    337 
    338                 /*Free penalties: */
    339                 xfree((void**)&iomodel->penalties);
    340         }
    341 
    342         /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
    343          We can therefore determine  which grids are internal to this node's partition
    344          and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
    345          that, go and create the grids*/
    346 
    347         /*Create nodes from x,y,z, as well as the spc values on those grids: */
    348                
    349         /*First fetch data: */
     82        /*Create nodes and vertices: */
    35083        if (strcmp(iomodel->meshtype,"3d")==0){
    35184                IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
     
    36295        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
    36396
    364         /*Get number of dofs per node: */
    365         DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
     97        for (i=0;i<iomodel->numberofvertices;i++){
    36698
    367         for (i=0;i<iomodel->numberofnodes;i++){
    368         #ifdef _PARALLEL_
    369         /*keep only this partition's nodes:*/
    370         if((my_grids[i]==1)){
    371         #endif
     99                /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
     100                if(iomodel->my_vertices[i]){
     101                       
     102                        /*Add vertex to vertices dataset: */
     103                        vertices->AddObject(new Vertex(i,iomodel));
    372104
    373                 /*create vertex: */
    374                 vertex_id=i+1;
    375                 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]));
    376                 vertices->AddObject(vertex);
     105                        /*Add node to nodes dataset: */
     106                        nodes->AddObject(new Node(i,iomodel));
    377107
    378 
    379                 node_id=i+1; //matlab indexing
    380                
    381                 #ifdef _PARALLEL_
    382                 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
    383                         node_partitionborder=1;
    384108                }
    385                 else{
    386                         node_partitionborder=0;
    387                 }
    388                 #else
    389                         node_partitionborder=0;
    390                 #endif
    391 
    392                 node_properties.SetProperties(
    393                         (int)iomodel->gridonbed[i],
    394                         (int)iomodel->gridonsurface[i],
    395                         (int)iomodel->gridoniceshelf[i],
    396                         (int)iomodel->gridonicesheet[i]);
    397 
    398 
    399                 if (strcmp(iomodel->meshtype,"3d")==0){
    400                         if (isnan(iomodel->uppernodes[i])){
    401                                 node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
    402                         }
    403                         else{
    404                                 node_upper_node_id=(int)iomodel->uppernodes[i];
    405                         }
    406                 }
    407                 else{
    408                         /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
    409                         node_upper_node_id=node_id;
    410                 }
    411 
    412                 /*Create node using its constructor: */
    413                 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
    414 
    415                 /*Add node to nodes dataset: */
    416                 nodes->AddObject(node);
    417 
    418         #ifdef _PARALLEL_
    419         } //if((my_grids[i]==1))
    420         #endif
    421109        }
    422 
    423         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    424          * datasets, it will not be redone: */
    425         elements->Presort();
    426         nodes->Presort();
    427         vertices->Presort();
    428         materials->Presort();
    429110
    430111        /*Clean fetched data: */
     
    441122        xfree((void**)&iomodel->gridoniceshelf);
    442123       
     124        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     125         * datasets, it will not be redone: */
     126        elements->Presort();
     127        nodes->Presort();
     128        vertices->Presort();
     129        materials->Presort();
     130
    443131        cleanup_and_return:
    444132
     
    449137        *pmaterials=materials;
    450138
    451         /*Keep partitioning information into iomodel*/
    452         iomodel->epart=epart;
    453         iomodel->my_grids=my_grids;
    454         iomodel->my_bordergrids=my_bordergrids;
    455 
    456         /*Free ressources:*/
    457         #ifdef _PARALLEL_
    458         xfree((void**)&all_numgrids);
    459         xfree((void**)&npart);
    460         VecFree(&gridborder);
    461         #endif
    462 
    463139}
Note: See TracChangeset for help on using the changeset viewer.