Changeset 1832


Ignore:
Timestamp:
08/24/09 17:53:21 (15 years ago)
Author:
Eric.Larour
Message:

Brachning back from issm.controlstatic, by hand

Location:
issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r1762 r1832  
    1111#include "../../objects/objects.h"
    1212#include "../../shared/shared.h"
    13 #include "../Model.h"
     13#include "../IoModel.h"
    1414
    1515
    16 void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
     16void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1717
    1818
    … …  
    4444
    4545        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    46         if (!model->ismacayealpattyn)goto cleanup_and_return;
     46        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
    4747
    4848        /*Fetch data: */
    49         ModelFetchData((void**)&spcvelocity,NULL,NULL,model_handle,"spcvelocity","Matrix","Mat");
    50         ModelFetchData((void**)&gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
     49        IoModelFetchData((void**)&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity","Matrix","Mat");
     50        IoModelFetchData((void**)&gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter","Matrix","Mat");
    5151
    5252        count=0;
    5353
    5454        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    55         for (i=0;i<model->numberofnodes;i++){
     55        for (i=0;i<iomodel->numberofnodes;i++){
    5656        #ifdef _PARALLEL_
    5757        /*keep only this partition's nodes:*/
    58         if((model->my_grids[i]==1)){
     58        if((iomodel->my_grids[i]==1)){
    5959        #endif
    6060
    … …  
    6666                        spc_node=i+1;
    6767                        spc_dof=1; //we enforce first x translation degree of freedom
    68                         spc_value=*(spcvelocity+6*i+3)/model->yts;
     68                        spc_value=*(spcvelocity+6*i+3)/iomodel->yts;
    6969
    7070                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
    … …  
    8080                        spc_node=i+1;
    8181                        spc_dof=2; //we enforce first y translation degree of freedom
    82                         spc_value=*(spcvelocity+6*i+4)/model->yts;
     82                        spc_value=*(spcvelocity+6*i+4)/iomodel->yts;
    8383                       
    8484                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
    … …  
    9595
    9696        /*First fetch penalties: */
    97         if (strcmp(model->meshtype,"3d")==0){
    98                 ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
     97        if (strcmp(iomodel->meshtype,"3d")==0){
     98                IoModelFetchData((void**)&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties","Matrix","Mat");
    9999        }
    100100
    101101        count=0;
    102         if(model->numpenalties){
     102        if(iomodel->numpenalties){
    103103
    104104                /*First deal with internal grids: */
    105                 for (i=0;i<model->numpenalties;i++){
    106                         if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
     105                for (i=0;i<iomodel->numpenalties;i++){
     106                        if (iomodel->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
    107107
    108                                 for(j=1;j<model->numlayers;j++){
     108                                for(j=1;j<iomodel->numlayers;j++){
    109109                                        /*We are pairing grids along a vertical profile.*/
    110                                         grid1=(int)*(model->penalties+model->numlayers*i);
    111                                         grid2=(int)*(model->penalties+model->numlayers*i+j);
     110                                        grid1=(int)*(iomodel->penalties+iomodel->numlayers*i);
     111                                        grid2=(int)*(iomodel->penalties+iomodel->numlayers*i+j);
    112112
    113113                                        rgb_id=count+1; //matlab indexing
    … …  
    132132
    133133        /*Free data: */
    134         xfree((void**)&model->penalties);
     134        xfree((void**)&iomodel->penalties);
    135135
    136136        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these

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

    r1784 r1832  
    1212#include "../../shared/shared.h"
    1313#include "../../MeshPartitionx/MeshPartitionx.h"
    14 #include "../Model.h"
    15 
    16 
    17 void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
     14#include "../IoModel.h"
     15
     16
     17void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1818
    1919
    … …  
    149149
    150150        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    151         if (!model->ismacayealpattyn)goto cleanup_and_return;
     151        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
    152152
    153153        /*Width of elements: */
    154         if(strcmp(model->meshtype,"2d")==0){
     154        if(strcmp(iomodel->meshtype,"2d")==0){
    155155                elements_width=3; //tria elements
    156156        }
    … …  
    163163        #ifdef _PARALLEL_
    164164        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
    165         if(strcmp(model->meshtype,"2d")==0){
     165        if(strcmp(iomodel->meshtype,"2d")==0){
    166166                /*load elements: */
    167                 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     167                IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
    168168        }
    169169        else{
    170170                /*load elements2d: */
    171                 ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
    172         }
    173 
    174         MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
     171                IoModelFetchData((void**)&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d","Matrix","Mat");
     172        }
     173
     174        MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
    175175
    176176        /*Free elements and elements2d: */
    177         xfree((void**)&model->elements);
    178         xfree((void**)&model->elements2d);
     177        xfree((void**)&iomodel->elements);
     178        xfree((void**)&iomodel->elements2d);
    179179
    180180
    181181        /*Deal with rifts, they have to be included into one partition only, not several: */
    182         ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
    183 
    184         for(i=0;i<model->numrifts;i++){
    185                 el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
    186                 el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
     182        IoModelFetchData((void**)&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo","Matrix","Mat");
     183
     184        for(i=0;i<iomodel->numrifts;i++){
     185                el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
     186                el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
    187187                epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
    188188        }
    189189        /*Free rifts: */
    190         xfree((void**)&model->riftinfo);
     190        xfree((void**)&iomodel->riftinfo);
    191191
    192192        /*Used later on: */
    193         my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
     193        my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
    194194        #endif
    195195
    … …  
    199199
    200200        /*2d mesh: */
    201         if (strcmp(model->meshtype,"2d")==0){
     201        if (strcmp(iomodel->meshtype,"2d")==0){
    202202
    203203                /*Fetch data needed: */
    204                 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
    205                 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
    206                 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
    207                 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
    208                 ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
    209                 ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
    210                 ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
    211                 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
    212                 ModelFetchData((void**)&model->elementonwater,NULL,NULL,model_handle,"elementonwater","Matrix","Mat");
    213                 ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
    214                 ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
    215                 ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
    216                 ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
    217                 ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
    218                
    219                 for (i=0;i<model->numberofelements;i++){
     204                IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
     205                IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
     206                IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
     207                IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
     208                IoModelFetchData((void**)&iomodel->drag,NULL,NULL,iomodel_handle,"drag","Matrix","Mat");
     209                IoModelFetchData((void**)&iomodel->p,NULL,NULL,iomodel_handle,"p","Matrix","Mat");
     210                IoModelFetchData((void**)&iomodel->q,NULL,NULL,iomodel_handle,"q","Matrix","Mat");
     211                IoModelFetchData((void**)&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf","Matrix","Mat");
     212                IoModelFetchData((void**)&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater","Matrix","Mat");
     213                IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
     214                IoModelFetchData((void**)&iomodel->B,NULL,NULL,iomodel_handle,"B","Matrix","Mat");
     215                IoModelFetchData((void**)&iomodel->n,NULL,NULL,iomodel_handle,"n","Matrix","Mat");
     216                IoModelFetchData((void**)&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation","Matrix","Mat");
     217                IoModelFetchData((void**)&iomodel->melting,NULL,NULL,iomodel_handle,"melting","Matrix","Mat");
     218               
     219                for (i=0;i<iomodel->numberofelements;i++){
    220220
    221221                #ifdef _PARALLEL_
    … …  
    224224                #endif
    225225                       
    226                         if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
     226                        if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
    227227                               
    228228                                /*ids: */
    229229                                tria_id=i+1; //matlab indexing.
    230230                                tria_mid=i+1; //refers to the corresponding material property card
    231                                 tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
     231                                tria_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
    232232
    233233                                /*vertices offsets: */
    234                                 tria_g[0]=(int)*(model->elements+elements_width*i+0);
    235                                 tria_g[1]=(int)*(model->elements+elements_width*i+1);
    236                                 tria_g[2]=(int)*(model->elements+elements_width*i+2);
     234                                tria_g[0]=(int)*(iomodel->elements+elements_width*i+0);
     235                                tria_g[1]=(int)*(iomodel->elements+elements_width*i+1);
     236                                tria_g[2]=(int)*(iomodel->elements+elements_width*i+2);
    237237
    238238                                /*thickness,surface and bed:*/
    239                                 tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
    240                                 tria_h[1]=*(model->thickness+  ((int)*(model->elements+elements_width*i+1)-1));
    241                                 tria_h[2]=*(model->thickness+  ((int)*(model->elements+elements_width*i+2)-1)) ;
    242 
    243                                 tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1));
    244                                 tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1));
    245                                 tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1));
    246 
    247                                 tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1));
    248                                 tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1));
    249                                 tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1));
     239                                tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
     240                                tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1));
     241                                tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
     242
     243                                tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1));
     244                                tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1));
     245                                tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1));
     246
     247                                tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1));
     248                                tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1));
     249                                tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1));
    250250
    251251                                /*basal drag:*/
    252                                 tria_friction_type=(int)model->drag_type;
    253 
    254                                 tria_k[0]=*(model->drag+        ((int)*(model->elements+elements_width*i+0)-1));
    255                                 tria_k[1]=*(model->drag+        ((int)*(model->elements+elements_width*i+1)-1));
    256                                 tria_k[2]=*(model->drag+        ((int)*(model->elements+elements_width*i+2)-1));
     252                                tria_friction_type=(int)iomodel->drag_type;
     253
     254                                tria_k[0]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+0)-1));
     255                                tria_k[1]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+1)-1));
     256                                tria_k[2]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+2)-1));
    257257                               
    258                                 tria_p=model->p[i];
    259                                 tria_q=model->q[i];
     258                                tria_p=iomodel->p[i];
     259                                tria_q=iomodel->q[i];
    260260
    261261                                /*meling and accumulation*/
    262                                 tria_melting[0]=*(model->melting+        ((int)*(model->elements+elements_width*i+0)-1));
    263                                 tria_melting[1]=*(model->melting+        ((int)*(model->elements+elements_width*i+1)-1));
    264                                 tria_melting[2]=*(model->melting+        ((int)*(model->elements+elements_width*i+2)-1));
    265 
    266                                 tria_accumulation[0]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+0)-1));
    267                                 tria_accumulation[1]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+1)-1));
    268                                 tria_accumulation[2]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+2)-1));
     262                                tria_melting[0]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+0)-1));
     263                                tria_melting[1]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+1)-1));
     264                                tria_melting[2]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+2)-1));
     265
     266                                tria_accumulation[0]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+0)-1));
     267                                tria_accumulation[1]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+1)-1));
     268                                tria_accumulation[2]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+2)-1));
    269269
    270270                                /*element on iceshelf, water?:*/
    271                                 tria_shelf=(int)*(model->elementoniceshelf+i);
    272                                 tria_onwater=(bool)*(model->elementonwater+i);
    273 
    274                                 tria_meanvel=model->meanvel;
    275                                 tria_epsvel=model->epsvel;
     271                                tria_shelf=(int)*(iomodel->elementoniceshelf+i);
     272                                tria_onwater=(bool)*(iomodel->elementonwater+i);
     273
     274                                tria_meanvel=iomodel->meanvel;
     275                                tria_epsvel=iomodel->epsvel;
    276276
    277277                                /*viscosity_overshoot*/
    278                                 tria_viscosity_overshoot=model->viscosity_overshoot;
     278                                tria_viscosity_overshoot=iomodel->viscosity_overshoot;
    279279
    280280                                /*Create tria element using its constructor:*/
    … …  
    290290                                B_avg=0;
    291291                                for(j=0;j<3;j++){
    292                                         B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
     292                                        B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
    293293                                }
    294294                                B_avg=B_avg/3;
    295295                                matice_B=B_avg;
    296                                 matice_n=(double)*(model->n+i);
     296                                matice_n=(double)*(iomodel->n+i);
    297297                               
    298298                                /*Create matice using its constructor:*/
    … …  
    309309                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    310310                         will hold which grids belong to this partition*/
    311                         my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
    312                         my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
    313                         my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     311                        my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
     312                        my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
     313                        my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
    314314                        #endif
    315315
    … …  
    322322       
    323323                /*Free data : */
    324                 xfree((void**)&model->elements);
    325                 xfree((void**)&model->thickness);
    326                 xfree((void**)&model->surface);
    327                 xfree((void**)&model->bed);
    328                 xfree((void**)&model->drag);
    329                 xfree((void**)&model->p);
    330                 xfree((void**)&model->q);
    331                 xfree((void**)&model->elementoniceshelf);
    332                 xfree((void**)&model->elementonwater);
    333                 xfree((void**)&model->B);
    334                 xfree((void**)&model->n);
    335                 xfree((void**)&model->elements_type);
     324                xfree((void**)&iomodel->elements);
     325                xfree((void**)&iomodel->thickness);
     326                xfree((void**)&iomodel->surface);
     327                xfree((void**)&iomodel->bed);
     328                xfree((void**)&iomodel->drag);
     329                xfree((void**)&iomodel->p);
     330                xfree((void**)&iomodel->q);
     331                xfree((void**)&iomodel->elementoniceshelf);
     332                xfree((void**)&iomodel->elementonwater);
     333                xfree((void**)&iomodel->B);
     334                xfree((void**)&iomodel->n);
     335                xfree((void**)&iomodel->elements_type);
    336336
    337337        }
    … …  
    339339
    340340                /*Fetch data needed: */
    341                 ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
    342                 ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
    343                 ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
    344                 ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
    345                 ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
    346                 ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
    347                 ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
    348                 ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
    349                 ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
    350                 ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
    351                 ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
    352                 ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
    353                 ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
    354                 ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
    355                 ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
    356                 ModelFetchData((void**)&model->elementonwater,NULL,NULL,model_handle,"elementonwater","Matrix","Mat");
    357                
    358                 for (i=0;i<model->numberofelements;i++){
     341                IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
     342                IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
     343                IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
     344                IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
     345                IoModelFetchData((void**)&iomodel->drag,NULL,NULL,iomodel_handle,"drag","Matrix","Mat");
     346                IoModelFetchData((void**)&iomodel->p,NULL,NULL,iomodel_handle,"p","Matrix","Mat");
     347                IoModelFetchData((void**)&iomodel->q,NULL,NULL,iomodel_handle,"q","Matrix","Mat");
     348                IoModelFetchData((void**)&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf","Matrix","Mat");
     349                IoModelFetchData((void**)&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed","Matrix","Mat");
     350                IoModelFetchData((void**)&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface","Matrix","Mat");
     351                IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
     352                IoModelFetchData((void**)&iomodel->B,NULL,NULL,iomodel_handle,"B","Matrix","Mat");
     353                IoModelFetchData((void**)&iomodel->n,NULL,NULL,iomodel_handle,"n","Matrix","Mat");
     354                IoModelFetchData((void**)&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation","Matrix","Mat");
     355                IoModelFetchData((void**)&iomodel->melting,NULL,NULL,iomodel_handle,"melting","Matrix","Mat");
     356                IoModelFetchData((void**)&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater","Matrix","Mat");
     357               
     358                for (i=0;i<iomodel->numberofelements;i++){
    359359                #ifdef _PARALLEL_
    360360                /*We are using our element partition to decide which elements will be created on this node: */
    … …  
    362362                #endif
    363363
    364                         if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum() | *(model->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
     364                        if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
    365365                       
    366366                                /*name and id: */
    367367                                penta_id=i+1; //matlab indexing.
    368368                                penta_mid=i+1; //refers to the corresponding material property card
    369                                 penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
     369                                penta_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
    370370
    371371                                /*vertices,thickness,surface,bed and drag: */
    372372                                for(j=0;j<6;j++){
    373                                         penta_g[j]=(int)*(model->elements+elements_width*i+j);
    374                                         penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1));
    375                                         penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1));
    376                                         penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1));
    377                                         penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1));
    378                                         penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
    379                                         penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
     373                                        penta_g[j]=(int)*(iomodel->elements+elements_width*i+j);
     374                                        penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1));
     375                                        penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1));
     376                                        penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1));
     377                                        penta_k[j]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+j)-1));
     378                                        penta_melting[j]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+j)-1));
     379                                        penta_accumulation[j]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+j)-1));
    380380                                }
    381381
    382382                                /*basal drag:*/
    383                                 penta_friction_type=(int)model->drag_type;
    384                
    385                                 penta_p=model->p[i];
    386                                 penta_q=model->q[i];
     383                                penta_friction_type=(int)iomodel->drag_type;
     384               
     385                                penta_p=iomodel->p[i];
     386                                penta_q=iomodel->q[i];
    387387
    388388                                /*diverse: */
    389                                 penta_shelf=(int)*(model->elementoniceshelf+i);
    390                                 penta_onbed=(int)*(model->elementonbed+i);
    391                                 penta_onsurface=(int)*(model->elementonsurface+i);
    392                                 penta_meanvel=model->meanvel;
    393                                 penta_epsvel=model->epsvel;
    394                                 penta_onwater=(bool)*(model->elementonwater+i);
     389                                penta_shelf=(int)*(iomodel->elementoniceshelf+i);
     390                                penta_onbed=(int)*(iomodel->elementonbed+i);
     391                                penta_onsurface=(int)*(iomodel->elementonsurface+i);
     392                                penta_meanvel=iomodel->meanvel;
     393                                penta_epsvel=iomodel->epsvel;
     394                                penta_onwater=(bool)*(iomodel->elementonwater+i);
    395395                               
    396396                                /*viscosity_overshoot*/
    397                                 penta_viscosity_overshoot=model->viscosity_overshoot;
    398 
    399                                 if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
     397                                penta_viscosity_overshoot=iomodel->viscosity_overshoot;
     398
     399                                if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
    400400                                        penta_collapse=1;
    401401                                }
    … …  
    420420                                B_avg=0;
    421421                                for(j=0;j<6;j++){
    422                                         B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
     422                                        B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
    423423                                }
    424424                                B_avg=B_avg/6;
    425425                                matice_B= B_avg;
    426                                 matice_n=(double)*(model->n+i);
     426                                matice_n=(double)*(iomodel->n+i);
    427427               
    428428                                /*Create matice using its constructor:*/
    … …  
    439439                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    440440                         will hold which grids belong to this partition*/
    441                         my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
    442                         my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
    443                         my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
    444                         my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
    445                         my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
    446                         my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
     441                        my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
     442                        my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
     443                        my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
     444                        my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
     445                        my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
     446                        my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
    447447                        #endif
    448448
    … …  
    454454
    455455                /*Free data: */
    456                 xfree((void**)&model->elements);
    457                 xfree((void**)&model->thickness);
    458                 xfree((void**)&model->surface);
    459                 xfree((void**)&model->bed);
    460                 xfree((void**)&model->drag);
    461                 xfree((void**)&model->p);
    462                 xfree((void**)&model->q);
    463                 xfree((void**)&model->elementoniceshelf);
    464                 xfree((void**)&model->elementonbed);
    465                 xfree((void**)&model->elementonsurface);
    466                 xfree((void**)&model->elements_type);
    467                 xfree((void**)&model->n);
    468                 xfree((void**)&model->B);
    469                 xfree((void**)&model->accumulation);
    470                 xfree((void**)&model->melting);
    471                 xfree((void**)&model->elementonwater);
     456                xfree((void**)&iomodel->elements);
     457                xfree((void**)&iomodel->thickness);
     458                xfree((void**)&iomodel->surface);
     459                xfree((void**)&iomodel->bed);
     460                xfree((void**)&iomodel->drag);
     461                xfree((void**)&iomodel->p);
     462                xfree((void**)&iomodel->q);
     463                xfree((void**)&iomodel->elementoniceshelf);
     464                xfree((void**)&iomodel->elementonbed);
     465                xfree((void**)&iomodel->elementonsurface);
     466                xfree((void**)&iomodel->elements_type);
     467                xfree((void**)&iomodel->n);
     468                xfree((void**)&iomodel->B);
     469                xfree((void**)&iomodel->accumulation);
     470                xfree((void**)&iomodel->melting);
     471                xfree((void**)&iomodel->elementonwater);
    472472
    473473        } //if (strcmp(meshtype,"2d")==0)
    474474
    475475        /*Add one constant material property to materials: */
    476         matpar_mid=model->numberofelements+1; //put it at the end of the materials
    477         matpar_g=model->g;
    478         matpar_rho_ice=model->rho_ice;
    479         matpar_rho_water=model->rho_water;
    480         matpar_thermalconductivity=model->thermalconductivity;
    481         matpar_heatcapacity=model->heatcapacity;
    482         matpar_latentheat=model->latentheat;
    483         matpar_beta=model->beta;
    484         matpar_meltingpoint=model->meltingpoint;
    485         matpar_mixed_layer_capacity=model->mixed_layer_capacity;
    486         matpar_thermal_exchange_velocity=model->thermal_exchange_velocity;
     476        matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
     477        matpar_g=iomodel->g;
     478        matpar_rho_ice=iomodel->rho_ice;
     479        matpar_rho_water=iomodel->rho_water;
     480        matpar_thermalconductivity=iomodel->thermalconductivity;
     481        matpar_heatcapacity=iomodel->heatcapacity;
     482        matpar_latentheat=iomodel->latentheat;
     483        matpar_beta=iomodel->beta;
     484        matpar_meltingpoint=iomodel->meltingpoint;
     485        matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
     486        matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
    487487
    488488        /*Create matpar object using its constructor: */
    … …  
    497497                /*From the element partitioning, we can determine which grids are on the inside of this cpu's
    498498                 *element partition, and which are on its border with other nodes:*/
    499                 gridborder=NewVec(model->numberofnodes);
    500 
    501                 for (i=0;i<model->numberofnodes;i++){
     499                gridborder=NewVec(iomodel->numberofnodes);
     500
     501                for (i=0;i<iomodel->numberofnodes;i++){
    502502                        if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
    503503                }
    … …  
    513513                #ifdef _DEBUG_
    514514                if(my_rank==0){
    515                         for (i=0;i<model->numberofnodes;i++){
     515                        for (i=0;i<iomodel->numberofnodes;i++){
    516516                                printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
    517517                        }
    … …  
    521521
    522522        /*Partition penalties in 3d: */
    523         if(strcmp(model->meshtype,"3d")==0){
     523        if(strcmp(iomodel->meshtype,"3d")==0){
    524524       
    525525                /*Get penalties: */
    526                 ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
    527 
    528                 if(model->numpenalties){
    529 
    530                         model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
     526                IoModelFetchData((void**)&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties","Matrix","Mat");
     527
     528                if(iomodel->numpenalties){
     529
     530                        iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
    531531                        #ifdef _SERIAL_
    532                         for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
     532                        for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
    533533                        #else
    534                         for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
    535 
    536                         for(i=0;i<model->numpenalties;i++){
    537                                 first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
     534                        for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
     535
     536                        for(i=0;i<iomodel->numpenalties;i++){
     537                                first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
    538538                                if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
    539539                                        /*All grids that are being penalised belong to this node's internal grid partition.:*/
    540                                         model->penaltypartitioning[i]=1;
     540                                        iomodel->penaltypartitioning[i]=1;
    541541                                }
    542542                                if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
    543                                         model->penaltypartitioning[i]=0;
     543                                        iomodel->penaltypartitioning[i]=0;
    544544                                }
    545545                        }
    … …  
    548548
    549549                /*Free penalties: */
    550                 xfree((void**)&model->penalties);
     550                xfree((void**)&iomodel->penalties);
    551551        }
    552552
    … …  
    559559               
    560560        /*First fetch data: */
    561         if (strcmp(model->meshtype,"3d")==0){
    562                 ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
    563                 ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
    564         }
    565         ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
    566         ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
    567         ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
    568         ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
    569         ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
    570         ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
    571         ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
    572         ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
    573         ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
    574         ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
     561        if (strcmp(iomodel->meshtype,"3d")==0){
     562                IoModelFetchData((void**)&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids","Matrix","Mat");
     563                IoModelFetchData((void**)&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids","Matrix","Mat");
     564        }
     565        IoModelFetchData((void**)&iomodel->x,NULL,NULL,iomodel_handle,"x","Matrix","Mat");
     566        IoModelFetchData((void**)&iomodel->y,NULL,NULL,iomodel_handle,"y","Matrix","Mat");
     567        IoModelFetchData((void**)&iomodel->z,NULL,NULL,iomodel_handle,"z","Matrix","Mat");
     568        IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
     569        IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
     570        IoModelFetchData((void**)&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed","Matrix","Mat");
     571        IoModelFetchData((void**)&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface","Matrix","Mat");
     572        IoModelFetchData((void**)&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter","Matrix","Mat");
     573        IoModelFetchData((void**)&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet","Matrix","Mat");
     574        IoModelFetchData((void**)&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf","Matrix","Mat");
    575575
    576576        /*Get number of dofs per node: */
    577         DistributeNumDofs(&node_numdofs,model->analysis_type,model->sub_analysis_type);
    578 
    579         for (i=0;i<model->numberofnodes;i++){
     577        DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
     578
     579        for (i=0;i<iomodel->numberofnodes;i++){
    580580        #ifdef _PARALLEL_
    581581
    … …  
    598598                #endif
    599599
    600                 node_x[0]=model->x[i];
    601                 node_x[1]=model->y[i];
    602                 node_x[2]=model->z[i];
    603                 node_sigma=(model->z[i]-model->bed[i])/(model->thickness[i]);
    604                
    605                 node_onbed=(int)model->gridonbed[i];
    606                 node_onsurface=(int)model->gridonsurface[i];   
    607                 node_onshelf=(int)model->gridoniceshelf[i];     
    608                 node_onsheet=(int)model->gridonicesheet[i];     
    609 
    610                 if (strcmp(model->meshtype,"3d")==0){
    611                         if (isnan(model->uppernodes[i])){
     600                node_x[0]=iomodel->x[i];
     601                node_x[1]=iomodel->y[i];
     602                node_x[2]=iomodel->z[i];
     603                node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
     604               
     605                node_onbed=(int)iomodel->gridonbed[i];
     606                node_onsurface=(int)iomodel->gridonsurface[i]; 
     607                node_onshelf=(int)iomodel->gridoniceshelf[i];   
     608                node_onsheet=(int)iomodel->gridonicesheet[i];   
     609
     610                if (strcmp(iomodel->meshtype,"3d")==0){
     611                        if (isnan(iomodel->uppernodes[i])){
    612612                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
    613613                        }
    614614                        else{
    615                                 node_upper_node_id=(int)model->uppernodes[i];
     615                                node_upper_node_id=(int)iomodel->uppernodes[i];
    616616                        }
    617617                }
    … …  
    625625
    626626                /*set single point constraints.: */
    627                 if (strcmp(model->meshtype,"3d")==0){
     627                if (strcmp(iomodel->meshtype,"3d")==0){
    628628                        /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
    629                         if (model->deadgrids[i]){
     629                        if (iomodel->deadgrids[i]){
    630630                                for(k=1;k<=node_numdofs;k++){
    631631                                        node->FreezeDof(k);
    … …  
    633633                        }
    634634                }
    635                 if (model->gridonhutter[i]){
     635                if (iomodel->gridonhutter[i]){
    636636                        for(k=1;k<=node_numdofs;k++){
    637637                                node->FreezeDof(k);
    … …  
    653653
    654654        /*Clean fetched data: */
    655         xfree((void**)&model->deadgrids);
    656         xfree((void**)&model->x);
    657         xfree((void**)&model->y);
    658         xfree((void**)&model->z);
    659         xfree((void**)&model->thickness);
    660         xfree((void**)&model->bed);
    661         xfree((void**)&model->gridonbed);
    662         xfree((void**)&model->gridonsurface);
    663         xfree((void**)&model->gridonhutter);
    664         xfree((void**)&model->uppernodes);
    665         xfree((void**)&model->gridonicesheet);
    666         xfree((void**)&model->gridoniceshelf);
     655        xfree((void**)&iomodel->deadgrids);
     656        xfree((void**)&iomodel->x);
     657        xfree((void**)&iomodel->y);
     658        xfree((void**)&iomodel->z);
     659        xfree((void**)&iomodel->thickness);
     660        xfree((void**)&iomodel->bed);
     661        xfree((void**)&iomodel->gridonbed);
     662        xfree((void**)&iomodel->gridonsurface);
     663        xfree((void**)&iomodel->gridonhutter);
     664        xfree((void**)&iomodel->uppernodes);
     665        xfree((void**)&iomodel->gridonicesheet);
     666        xfree((void**)&iomodel->gridoniceshelf);
    667667       
    668668
    669         /*Keep partitioning information into model*/
    670         model->epart=epart;
    671         model->my_grids=my_grids;
    672         model->my_bordergrids=my_bordergrids;
     669        /*Keep partitioning information into iomodel*/
     670        iomodel->epart=epart;
     671        iomodel->my_grids=my_grids;
     672        iomodel->my_bordergrids=my_bordergrids;
    673673
    674674        /*Free ressources:*/

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

    r1805 r1832  
    1111#include "../../shared/shared.h"
    1212#include "../../include/macros.h"
    13 #include "../Model.h"
    14 
    15 
    16 void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){
     13#include "../IoModel.h"
     14
     15
     16void    CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1717
    1818        int i;
    … …  
    5555        double riftfront_friction;
    5656        double riftfront_fraction;
    57         double riftfront_fractionincrement;
    5857        bool   riftfront_shelf;
    5958        double riftfront_penalty_offset;
    … …  
    6968        double friction;
    7069        double fraction;
    71         double fractionincrement;
    7270       
    7371        /*Create loads: */
    … …  
    7573
    7674        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    77         if (!model->ismacayealpattyn)goto cleanup_and_return;
     75        if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
    7876       
    7977        /*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
    8078         * referenced by a certain load must belong to the cluster node): */
    81         ModelFetchData((void**)&model->pressureload,&model->numberofpressureloads,NULL,model_handle,"pressureload","Matrix","Mat");
    82         ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
    83         ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
    84         ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
     79        IoModelFetchData((void**)&iomodel->pressureload,&iomodel->numberofpressureloads,NULL,iomodel_handle,"pressureload","Matrix","Mat");
     80        IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
     81        IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
     82        IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
    8583
    8684        /*First load data:*/
    87         for (i=0;i<model->numberofpressureloads;i++){
    88                
    89                 if (strcmp(model->meshtype,"2d")==0){
     85        for (i=0;i<iomodel->numberofpressureloads;i++){
     86               
     87                if (strcmp(iomodel->meshtype,"2d")==0){
    9088                        segment_width=3;
    9189                        element_type=TriaEnum();
    … …  
    9795
    9896
    99                 element=(int)(*(model->pressureload+segment_width*i+segment_width-1)-1); //element is in the last column
     97                element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-1)-1); //element is in the last column
    10098
    10199                #ifdef _PARALLEL_
    102                 if (model->epart[element]!=my_rank){
     100                if (iomodel->epart[element]!=my_rank){
    103101                        /*This load does not belong to this cluster node, as it references an element
    104102                         *that does not belong to this node's partition. Drop this 'i':*/
    … …  
    107105                #endif
    108106       
    109                 icefront_mparid=model->numberofelements+1; //matlab indexing
     107                icefront_mparid=iomodel->numberofelements+1; //matlab indexing
    110108                icefront_sid=i+1; //matlab indexing
    111                 icefront_eid=(int)*(model->pressureload+segment_width*i+segment_width-1); //matlab indexing
     109                icefront_eid=(int)*(iomodel->pressureload+segment_width*i+segment_width-1); //matlab indexing
    112110                icefront_element_type=element_type;
    113111
    114                 i1=(int)*(model->pressureload+segment_width*i+0);
    115                 i2=(int)*(model->pressureload+segment_width*i+1);
     112                i1=(int)*(iomodel->pressureload+segment_width*i+0);
     113                i2=(int)*(iomodel->pressureload+segment_width*i+1);
    116114                       
    117115                icefront_node_ids[0]=i1;
    118116                icefront_node_ids[1]=i2;
    119117               
    120                 icefront_h[0]=model->thickness[i1-1];
    121                 icefront_h[1]=model->thickness[i2-1];
    122 
    123                 icefront_b[0]=model->bed[i1-1];
    124                 icefront_b[1]=model->bed[i2-1];
    125                
    126                 if ((int)*(model->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d
     118                icefront_h[0]=iomodel->thickness[i1-1];
     119                icefront_h[1]=iomodel->thickness[i2-1];
     120
     121                icefront_b[0]=iomodel->bed[i1-1];
     122                icefront_b[1]=iomodel->bed[i2-1];
     123               
     124                if ((int)*(iomodel->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d
    127125                        strcpy(icefront_type,"segment");
    128126                }
    129                 else if ((int)*(model->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d.
     127                else if ((int)*(iomodel->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d.
    130128                        strcpy(icefront_type,"quad");
    131                         i3=(int)*(model->pressureload+segment_width*i+2);
    132                         i4=(int)*(model->pressureload+segment_width*i+3);
     129                        i3=(int)*(iomodel->pressureload+segment_width*i+2);
     130                        i4=(int)*(iomodel->pressureload+segment_width*i+3);
    133131                        icefront_node_ids[2]=i3;
    134132                        icefront_node_ids[3]=i4;
    135133                       
    136                         icefront_h[2]=model->thickness[i3-1];
    137                         icefront_h[3]=model->thickness[i4-1];
    138 
    139                         icefront_b[2]=model->bed[i3-1];
    140                         icefront_b[3]=model->bed[i4-1];
     134                        icefront_h[2]=iomodel->thickness[i3-1];
     135                        icefront_h[3]=iomodel->thickness[i4-1];
     136
     137                        icefront_b[2]=iomodel->bed[i3-1];
     138                        icefront_b[3]=iomodel->bed[i4-1];
    141139                }
    142140                else{
    143                         throw ErrorException(__FUNCT__,exprintf(" element type %i not supported yet",(int)*(model->elements_type+2*element+0)));
     141                        throw ErrorException(__FUNCT__,exprintf(" element type %i not supported yet",(int)*(iomodel->elements_type+2*element+0)));
    144142                }
    145143
    … …  
    150148        }
    151149        /*Free data: */
    152         xfree((void**)&model->pressureload);
    153         xfree((void**)&model->elements_type);
    154         xfree((void**)&model->thickness);
    155         xfree((void**)&model->bed);
     150        xfree((void**)&iomodel->pressureload);
     151        xfree((void**)&iomodel->elements_type);
     152        xfree((void**)&iomodel->thickness);
     153        xfree((void**)&iomodel->bed);
    156154
    157155        /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
    158156        /*First fetch rifts: */
    159         ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
    160         ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
    161         ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
    162         ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
    163         ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
    164        
    165         if(model->numrifts){
    166                 for(i=0;i<model->numrifts;i++){
     157        IoModelFetchData((void**)&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo","Matrix","Mat");
     158        IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
     159        IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
     160        IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
     161        IoModelFetchData((void**)&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf","Matrix","Mat");
     162       
     163        if(iomodel->numrifts){
     164                for(i=0;i<iomodel->numrifts;i++){
    167165                               
    168                         el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2);
     166                        el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2);
    169167                        #ifdef _PARALLEL_
    170                         if (model->epart[el1-1]!=my_rank){
     168                        if (iomodel->epart[el1-1]!=my_rank){
    171169                                /*This load does not belong to this cluster node, as it references an element
    172170                                 *that does not belong to this node's partition. Drop this 'j':*/
    … …  
    176174
    177175                        /*Ok, retrieve all the data needed to add a penalty between the two grids: */
    178                         el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3);
    179 
    180                         grid1=(int)*(model->riftinfo+RIFTINFOSIZE*i+0);
    181                         grid2=(int)*(model->riftinfo+RIFTINFOSIZE*i+1);
    182                        
    183                         normal[0]=*(model->riftinfo+RIFTINFOSIZE*i+4);
    184                         normal[1]=*(model->riftinfo+RIFTINFOSIZE*i+5);
    185                         length=*(model->riftinfo+RIFTINFOSIZE*i+6);
    186                        
    187                         fill = (int)*(model->riftinfo+RIFTINFOSIZE*i+7);
    188                         friction=*(model->riftinfo+RIFTINFOSIZE*i+8);
    189                         fraction=*(model->riftinfo+RIFTINFOSIZE*i+9);
     176                        el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3);
     177
     178                        grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0);
     179                        grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1);
     180                       
     181                        normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4);
     182                        normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5);
     183                        length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
     184                       
     185                        fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7);
     186                        friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8);
     187                        fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
    190188                        fractionincrement=*(model->riftinfo+RIFTINFOSIZE*i+10);
    191189       
    … …  
    194192                        riftfront_node_ids[0]=grid1;
    195193                        riftfront_node_ids[1]=grid2;
    196                         riftfront_mparid=model->numberofelements+1; //matlab indexing
    197 
    198                         riftfront_h[0]=model->thickness[grid1-1];
    199                         riftfront_h[1]=model->thickness[grid2-1];
    200 
    201                         riftfront_b[0]=model->bed[grid1-1];
    202                         riftfront_b[1]=model->bed[grid2-1];
    203 
    204                         riftfront_s[0]=model->surface[grid1-1];
    205                         riftfront_s[1]=model->surface[grid2-1];
     194                        riftfront_mparid=iomodel->numberofelements+1; //matlab indexing
     195
     196                        riftfront_h[0]=iomodel->thickness[grid1-1];
     197                        riftfront_h[1]=iomodel->thickness[grid2-1];
     198
     199                        riftfront_b[0]=iomodel->bed[grid1-1];
     200                        riftfront_b[1]=iomodel->bed[grid2-1];
     201
     202                        riftfront_s[0]=iomodel->surface[grid1-1];
     203                        riftfront_s[1]=iomodel->surface[grid2-1];
    206204
    207205                        riftfront_normal[0]=normal[0];
    … …  
    215213                        riftfront_shelf=(bool)model->gridoniceshelf[grid1-1];
    216214
    217                         riftfront_penalty_offset=model->penalty_offset;
    218                         riftfront_penalty_lock=model->penalty_lock;
     215                        riftfront_penalty_offset=iomodel->penalty_offset;
     216                        riftfront_penalty_lock=iomodel->penalty_lock;
    219217
    220218                        riftfront_active=0;
    … …  
    222220                        riftfront_prestable=0;
    223221                       
    224                         riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
     222                        riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
    225223
    226224                        loads->AddObject(riftfront);
    … …  
    230228
    231229        /*free ressources: */
    232         xfree((void**)&model->riftinfo);
    233         xfree((void**)&model->thickness);
    234         xfree((void**)&model->bed);
    235         xfree((void**)&model->surface);
    236         xfree((void**)&model->gridoniceshelf);
     230        xfree((void**)&iomodel->riftinfo);
     231        xfree((void**)&iomodel->thickness);
     232        xfree((void**)&iomodel->bed);
     233        xfree((void**)&iomodel->surface);
     234        xfree((void**)&iomodel->gridoniceshelf);
    237235
    238236

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

    r800 r1832  
    1111#include "../../objects/objects.h"
    1212#include "../../shared/shared.h"
    13 #include "../Model.h"
     13#include "../IoModel.h"
    1414
    15 void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
     15void CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
    1616       
    1717        Param*   param = NULL;
    … …  
    3131
    3232        /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
    33 //      if (!model->ismacayealpattyn)return;
     33//      if (!iomodel->ismacayealpattyn)return;
    3434
    3535        /*Get vx and vy: */
    36         ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Matrix","Mat");
    37         ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Matrix","Mat");
    38         ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat");
     36        IoModelFetchData((void**)&vx,NULL,NULL,iomodel_handle,"vx","Matrix","Mat");
     37        IoModelFetchData((void**)&vy,NULL,NULL,iomodel_handle,"vy","Matrix","Mat");
     38        IoModelFetchData((void**)&vz,NULL,NULL,iomodel_handle,"vz","Matrix","Mat");
    3939
    40         ug=(double*)xcalloc(model->numberofnodes*3,sizeof(double));
     40        ug=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
    4141
    42         if(vx) for(i=0;i<model->numberofnodes;i++)ug[3*i+0]=vx[i]/model->yts;
    43         if(vy) for(i=0;i<model->numberofnodes;i++)ug[3*i+1]=vy[i]/model->yts;
    44         if(vz) for(i=0;i<model->numberofnodes;i++)ug[3*i+2]=vz[i]/model->yts;
     42        if(vx) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+0]=vx[i]/iomodel->yts;
     43        if(vy) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+1]=vy[i]/iomodel->yts;
     44        if(vz) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+2]=vz[i]/iomodel->yts;
    4545
    4646        count++;
    4747        param= new Param(count,"u_g",DOUBLEVEC);
    48         param->SetDoubleVec(ug,3*model->numberofnodes,3);
     48        param->SetDoubleVec(ug,3*iomodel->numberofnodes,3);
    4949        parameters->AddObject(param);
    5050
Note: See TracChangeset for help on using the changeset viewer.