Changeset 306


Ignore:
Timestamp:
05/07/09 15:06:35 (16 years ago)
Author:
Eric.Larour
Message:

temporary commit of new diagnostic hutter model processor

Location:
issm/trunk/src/c/ModelProcessorx
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp

    r300 r306  
    142142        int el1,el2;
    143143         
    144         /*Penalty partitioning: */
    145         int num_grids3d_collapsed;
    146         double* double_penalties_grids3d_collapsed=NULL;
    147         double* double_penalties_grids3d_noncollapsed=NULL;
    148         int     grid_ids[6];
    149         int     num_grid_ids;
    150         int     grid_id;
    151 
    152 
    153144        /*First create the elements, nodes and material properties: */
    154145        elements  = new DataSet(ElementsEnum());
    … …  
    336327       
    337328                /*Free data : */
    338                 /*xfree((void**)&model->elements);
     329                xfree((void**)&model->elements);
    339330                xfree((void**)&model->thickness);
    340331                xfree((void**)&model->surface);
    … …  
    345336                xfree((void**)&model->elementoniceshelf);
    346337                xfree((void**)&model->B);
    347                 xfree((void**)&model->n);*/
     338                xfree((void**)&model->n);
    348339        }
    349340        else{ //        if (strcmp(meshtype,"2d")==0)
    … …  
    668659        #ifdef _PARALLEL_
    669660        xfree((void**)&all_numgrids);
     661        xfree((void**)&npart);
    670662        VecFree(&gridborder);
    671663        #endif

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp

    r300 r306  
    1616void    CreateConstraintsDiagnosticHutter(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
    1717
     18        int i;
     19        int count;
    1820
    1921        DataSet* constraints = NULL;
    2022
     23        Spc*    spc  = NULL;
     24
     25        /*spc intermediary data: */
     26        int spc_sid;
     27        int spc_node;
     28        int spc_dof;
     29        double spc_value;
     30       
    2131        /*Create constraints: */
    2232        constraints = new DataSet(ConstraintsEnum());
    … …  
    2434        /*Now, is the flag ishutter on? otherwise, do nothing: */
    2535        if (!model->ishutter)goto cleanup_and_return;
     36
     37        count=0;
    2638       
    27         cleanup_and_return:     
     39        /*Fetch data: */
     40        ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
     41
     42        /*vx and vy are spc'd if we are not on gridonhutter: */
     43        for (i=0;i<model->numberofnodes;i++){
     44        #ifdef _PARALLEL_
     45        /*keep only this partition's nodes:*/
     46        if((model->npart[i]==1)){
     47        #endif
     48
     49                if ((int)model->gridonhutter[i]){
     50       
     51                        spc_sid=count;
     52                        spc_node=i+1;
     53                        spc_dof=1; //we enforce first x translation degree of freedom
     54                        spc_value=0;
     55
     56                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
     57                        constraints->AddObject(spc);
     58                        count++;
     59
     60                        spc_sid=count;
     61                        spc_node=i+1;
     62                        spc_dof=2; //we enforce first y translation degree of freedom
     63                        spc_value=0;
     64                       
     65                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
     66                        constraints->AddObject(spc);
     67                        count++;
     68                }
     69
     70        #ifdef _PARALLEL_
     71        } //if((npart[i]==1))
     72        #endif
     73        }
     74
     75        /*Free data: */
     76        xfree((void**)&model->gridonhutter);
     77
     78        //deal with mpcs for 2d-3d mesh transitions
     79       
     80        /*Fetch data: */
     81        ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
     82
     83        if (strcmp(model->meshtype,"3d")==0){
     84       
     85                if (model->numpenalties){
     86                        for (i=0;i<model->numpenalties;i++){
     87                                for (j=0;j<model->numlayers-1;j++){
     88
     89                                        //constrain first dof
     90                                        rgb_id=count;
     91                                        rgb_dof=1;
     92                                        rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j);
     93                                        rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j+1);
     94                                       
     95                                        rgb = new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
     96                                        constraints->AddObject(rgb);
     97                                        count++;
     98
     99                                        //constrain second dof
     100                                        rgb_id=count;
     101                                        rgb_dof=2;
     102                                        rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j);
     103                                        rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j+1);
     104                                       
     105                                        rgb = new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
     106                                        constraints->AddObject(rgb);
     107                                        count++;
     108       
     109                                }
     110                        }
     111                }
     112        }
     113
     114        /*Free data: */
     115        xfree((void**)&model->penalties);
     116
     117        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     118         * datasets, it will not be redone: */
     119        constraints->Presort();
     120
     121        cleanup_and_return:
    28122       
    29123        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp

    r300 r306  
    2828        materials = new DataSet(MaterialsEnum());
    2929
     30        /*Objects: */
     31        Node*       node   = NULL;
     32        Matice*     matice = NULL;
     33        Matpar*     matpar = NULL;
     34        Beam*       beam   = NULL;
     35        Sing*       sing   = NULL;
     36
     37        int         analysis_type;
     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        /*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        /*sing constructor input: */
     55        int   sing_id;
     56        int   sing_mid;
     57        int   sing_mparid;
     58        int   sing_g;
     59        double sing_h,sing_k;
     60
     61        /*beam constructor input: */
     62        int   beam_id;
     63        int   beam_mid;
     64        int   beam_mparid;
     65        int   beam_g[2];
     66        double beam_h[2];
     67        double beam_s[2];
     68        double beam_b[2];
     69        double beam_k[2];
     70                                       
     71        /*matpar constructor input: */
     72        int     matpar_mid;
     73        double matpar_rho_ice;
     74        double matpar_rho_water;
     75        double matpar_heatcapacity;
     76        double matpar_thermalconductivity;
     77        double matpar_latentheat;
     78        double matpar_beta;
     79        double matpar_meltingpoint;
     80        double matpar_mixed_layer_capacity;
     81        double matpar_thermal_exchange_velocity;
     82        double matpar_g;
     83
     84        /* node constructor input: */
     85        int node_id;
     86        int node_partitionborder=0;
     87        double node_x[3];
     88        int node_onbed;
     89        int node_onsurface;
     90        int node_upper_node_id;
     91        int node_numdofs;
     92
     93        /*First create the elements, nodes and material properties: */
     94        elements  = new DataSet(ElementsEnum());
     95        nodes     = new DataSet(NodesEnum());
     96        materials = new DataSet(MaterialsEnum());
     97
    3098        /*Now, is the flag ishutter on? otherwise, do nothing: */
    3199        if (!model->ishutter)goto cleanup_and_return;
    32100
     101        /*Get analysis_type: */
     102        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     103
     104        /*Width of elements: */
     105        if(strcmp(model->meshtype,"2d")==0){
     106                elements_width=3; //tria elements
     107        }
     108        else{
     109                elements_width=6; //penta elements
     110        }
     111
     112
     113        #ifdef _PARALLEL_
     114        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
     115        if(strcmp(model->meshtype,"2d")==0){
     116                /*load elements: */
     117                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     118        }
     119        else{
     120                /*load elements2d: */
     121                ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
     122        }
     123
     124        MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
     125
     126        /*Free elements and elements2d: */
     127        xfree((void**)&model->elements);
     128        xfree((void**)&model->elements2d);
     129
     130                       
     131        #endif
     132
     133        /*Hutter elements can be partitioned using npart, instead of epart for more classic tria and penta elements. The reason is
     134         * that each hutter elements either lies on a node (in 2d), or a pair of vertically juxtaposed nodes (in 3d): */
     135
     136        /*Fetch data temporarily needed: */
     137        ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
     138        ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
     139        ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
     140        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     141        ModelFetchData((void**)&model->uppergrids,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
     142        ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
     143        ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
     144        ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
     145
     146        /*2d mesh: */
     147        if (strcmp(model->meshtype,"2d")==0){
     148
     149                for (i=0;i<model->numberofnodes;i++){
     150                #ifdef _PARALLEL_
     151                /*keep only this partition's nodes:*/
     152                if(my_rank==npart[i]){
     153                #endif
     154
     155                        if(model->gridonhutter[i]){
     156                               
     157                                /*Deal with sing element: */
     158                                sing_id=i+1;
     159                                sing_mid=i+1; //refers to the corresponding material property card
     160                                sing_mparid=model->numberofnodes+1;//refers to the corresponding matpar property card
     161                                sing_g=i+1;
     162                                sing_h=model->thickness[i];
     163                                sing_k=model->drag[i];
     164
     165                                /*Create sing element using its constructor:*/
     166                                sing=new Sing(sing_id, sing_mid, sing_mparid, sing_g, sing_h, sing_k);
     167
     168                                /*Add tria element to elements dataset: */
     169                                elements->AddObject(sing);
     170
     171                                /*Deal with material property card: */
     172                                matice_mid=i+1; //same as the material id from the geom2 elements.
     173                                matice_B=model->B[i];   
     174                                matice_n=(double)model->n[i];
     175                       
     176                                /*Create matice using its constructor:*/
     177                                matice= new Matice(matice_mid,matice_B,matice_n);
     178       
     179                                /*Add matice element to materials dataset: */
     180                                materials->AddObject(matice);
     181                        }
     182
     183                #ifdef _PARALLEL_
     184                } //if(my_rank==npart[i])
     185                #endif
     186
     187                } //for (i=0;i<model->numberofnodes;i++)
     188        } //if (strcmp(model->meshtype,"2d")==0)
     189        else{
     190
     191                for (i=0;i<model->numberofnodes;i++){
     192                #ifdef _PARALLEL_
     193                /*keep only this partition's nodes:*/
     194                if(my_rank==npart[i]){
     195                #endif
     196
     197                        if(model->gridonhutter[i]){
     198
     199                                if(!model->gridonsurface[i]){
     200                                       
     201                                        /*Deal with sing element: */
     202                                        beam_id=i+1;
     203                                        beam_mid=i+1; //refers to the corresponding material property card
     204                                        beam_mparid=model->numberofnodes+1;//refers to the corresponding matpar property card
     205                                        beam_g[0]=i+1;
     206                                        beam_g[1]=model->uppergrids[i]; //grid that lays right on top
     207                                        beam_h[0]=model->thickness[i];
     208                                        beam_h[1]=model->thickness[model->uppergrids[i]-1];
     209                                        beam_s[0]=model->surface[i];
     210                                        beam_s[1]=model->surface[model->uppergrids[i]-1];
     211                                        beam_b[0]=model->bed[i];
     212                                        beam_b[1]=model->bed[model->uppergrids[i]-1];
     213                                        beam_k[0]=model->drag[i];
     214                                        beam_k[1]=model->drag[model->uppergrids[i]-1];
     215
     216                                        /*Create beam element ubeam its constructor:*/
     217                                        beam=new Beam(beam_id, beam_mid, beam_mparid, beam_g, beam_h, beam_s,beam_b,beam_k);
     218
     219                                        /*Add tria element to elements dataset: */
     220                                        elements->AddObject(beam);
     221
     222                                        /*Deal with material property card: */
     223                                        matice_mid=i+1; //same as the material id from the geom2 elements.
     224                                        matice_B=model->B[i];   
     225                                        matice_n=(double)model->n[i];
     226                               
     227                                        /*Create matice ubeam its constructor:*/
     228                                        matice= new Matice(matice_mid,matice_B,matice_n);
     229               
     230                                        /*Add matice element to materials dataset: */
     231                                        materials->AddObject(matice);
     232                                }
     233                        }
     234
     235                #ifdef _PARALLEL_
     236                } //if(my_rank==npart[i])
     237                #endif
     238
     239                } //for (i=0;i<model->numberofnodes;i++)
     240
     241        } //if (strcmp(model->meshtype,"2d")==0)
     242       
     243
     244        /*Free data: */
     245        xfree((void**)&model->gridonhutter);
     246        xfree((void**)&model->thickness);
     247        xfree((void**)&model->surface);
     248        xfree((void**)&model->gridonsurface);
     249        xfree((void**)&model->uppergrids);
     250        xfree((void**)&model->drag);
     251        xfree((void**)&model->B);
     252        xfree((void**)&model->n);
     253       
     254
     255        /*Add one constant material property to materials: */
     256        matpar_mid=model->numberofnodes+1; //put it at the end of the materials
     257        matpar_g=model->g;
     258        matpar_rho_ice=model->rho_ice;
     259        matpar_rho_water=model->rho_water;
     260        matpar_thermalconductivity=model->thermalconductivity;
     261        matpar_heatcapacity=model->heatcapacity;
     262        matpar_latentheat=model->latentheat;
     263        matpar_beta=model->beta;
     264        matpar_meltingpoint=model->meltingpoint;
     265        matpar_mixed_layer_capacity=model->mixed_layer_capacity;
     266        matpar_thermal_exchange_velocity=model->thermal_exchange_velocity;
     267
     268       
     269        /*Create nodes from x,y,z, as well as the spc values on those grids: */
     270               
     271        /*First fetch data: */
     272        if (strcmp(model->meshtype,"3d")==0){
     273                ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
     274                ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
     275        }
     276        ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
     277        ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
     278        ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
     279        ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
     280        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     281
     282       
     283        /*Get number of dofs per node: */
     284        DistributeNumDofs(&node_numdofs,analysis_type);
     285
     286        for (i=0;i<model->numberofnodes;i++){
     287        #ifdef _PARALLEL_
     288        /*keep only this partition's nodes:*/
     289        if(npart[i]==my_rank){
     290        #endif
     291
     292                node_id=i+1; //matlab indexing
     293                       
     294                node_partitionborder=0;
     295
     296                node_x[0]=model->x[i];
     297                node_x[1]=model->y[i];
     298                node_x[2]=model->z[i];
     299
     300               
     301                node_onbed=(int)model->gridonbed[i];
     302                node_onsurface=(int)model->gridonsurface[i];   
     303                if (strcmp(model->meshtype,"3d")==0){
     304                        if (isnan(model->uppernodes[i])){
     305                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     306                        }
     307                        else{
     308                                node_upper_node_id=(int)model->uppernodes[i];
     309                        }
     310                }
     311                else{
     312                        /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
     313                        node_upper_node_id=node_id;
     314                }
     315
     316                /*Create node using its constructor: */
     317                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
     318
     319                /*set single point constraints.: */
     320                if (strcmp(model->meshtype,"3d")==0){
     321                        /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
     322                        if (model->deadgrids[i]){
     323                                for(k=1;k<=node_numdofs;k++){
     324                                        node->FreezeDof(k);
     325                                }
     326                        }
     327                }
     328                /*Add node to nodes dataset: */
     329                nodes->AddObject(node);
     330
     331        #ifdef _PARALLEL_
     332        } //if(npart[i]==my_rank)
     333        #endif
     334        }
     335
     336        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     337         * datasets, it will not be redone: */
     338        elements->Presort();
     339        nodes->Presort();
     340        materials->Presort();
     341
     342        /*Clean fetched data: */
     343        xfree((void**)&model->deadgrids);
     344        xfree((void**)&model->x);
     345        xfree((void**)&model->y);
     346        xfree((void**)&model->z);
     347        xfree((void**)&model->gridonbed);
     348        xfree((void**)&model->gridonsurface);
     349        xfree((void**)&model->uppernodes);
     350               
    33351        cleanup_and_return:
    34352
    … …  
    37355        *pnodes=nodes;
    38356        *pmaterials=materials;
     357
     358        /*Keep partitioning information into model*/
     359        xfree((void**)&epart);
     360        model->npart=npart;
    39361}

  • issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp

    r300 r306  
    432432        #ifdef _PARALLEL_
    433433        xfree((void**)&all_numgrids);
     434        xfree((void**)&npart);
    434435        VecFree(&gridborder);
    435436        #endif

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r300 r306  
    147147        /*exterior data: */
    148148        int* epart; /*!element partitioning.*/
     149        int* npart; /*!node partitioning.*/
    149150        int* my_grids; /*! grids that belong to this cpu*/
    150151        double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/
Note: See TracChangeset for help on using the changeset viewer.