Changeset 653


Ignore:
Timestamp:
05/29/09 16:17:27 (16 years ago)
Author:
seroussi
Message:

fixed problems with mixed modes: fields needed, elements and constraints problems

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

Legend:

Unmodified
Added
Removed

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

    r586 r653  
    3131        double* dirichletvalues_diag=NULL;
    3232        double* gridondirichlet_diag=NULL;
     33        double* gridonhutter=NULL;
    3334       
    3435        /*Create constraints: */
    … …  
    4142        ModelFetchData((void**)&gridondirichlet_diag,NULL,NULL,model_handle,"gridondirichlet_diag","Matrix","Mat");
    4243        ModelFetchData((void**)&dirichletvalues_diag,NULL,NULL,model_handle,"dirichletvalues_diag","Matrix","Mat");
     44        ModelFetchData((void**)&gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
    4345
    4446        count=0;
    … …  
    5153        #endif
    5254
    53                 if ((int)gridondirichlet_diag[i]){
     55                if ((int)gridondirichlet_diag[i] | (int)gridonhutter[i]){
    5456       
    5557                        /*This grid needs to be spc'd to vx_obs and vy_obs:*/
    … …  
    8789        xfree((void**)&gridondirichlet_diag);
    8890        xfree((void**)&dirichletvalues_diag);
     91        xfree((void**)&gridonhutter);
    8992       
    9093        cleanup_and_return:

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

    r598 r653  
    229229                ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
    230230                ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
     231                ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
    231232                ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
    232233                ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
    … …  
    241242                #endif
    242243                       
    243                        
    244                         /*ids: */
    245                         tria_id=i+1; //matlab indexing.
    246                         tria_mid=i+1; //refers to the corresponding material property card
    247                         tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
    248 
    249                         /*vertices offsets: */
    250                         tria_g[0]=(int)*(model->elements+elements_width*i+0);
    251                         tria_g[1]=(int)*(model->elements+elements_width*i+1);
    252                         tria_g[2]=(int)*(model->elements+elements_width*i+2);
    253 
    254                         /*thickness,surface and bed:*/
    255                         tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
    256                         tria_h[1]=*(model->thickness+  ((int)*(model->elements+elements_width*i+1)-1));
    257                         tria_h[2]=*(model->thickness+  ((int)*(model->elements+elements_width*i+2)-1)) ;
    258 
    259                         tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1));
    260                         tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1));
    261                         tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1));
    262 
    263                         tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1));
    264                         tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1));
    265                         tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1));
    266 
    267                         /*basal drag:*/
    268                         tria_friction_type=(int)model->drag_type;
    269 
    270                         tria_k[0]=*(model->drag+        ((int)*(model->elements+elements_width*i+0)-1));
    271                         tria_k[1]=*(model->drag+        ((int)*(model->elements+elements_width*i+1)-1));
    272                         tria_k[2]=*(model->drag+        ((int)*(model->elements+elements_width*i+2)-1));
    273                        
    274                         tria_p=model->p[i];
    275                         tria_q=model->q[i];
    276 
    277                         /*meling and accumulation*/
    278                         tria_melting[0]=*(model->melting+        ((int)*(model->elements+elements_width*i+0)-1));
    279                         tria_melting[1]=*(model->melting+        ((int)*(model->elements+elements_width*i+1)-1));
    280                         tria_melting[2]=*(model->melting+        ((int)*(model->elements+elements_width*i+2)-1));
    281 
    282                         tria_accumulation[0]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+0)-1));
    283                         tria_accumulation[1]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+1)-1));
    284                         tria_accumulation[2]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+2)-1));
    285 
    286                         /*element on iceshelf?:*/
    287                         tria_shelf=(int)*(model->elementoniceshelf+i);
    288 
    289                         tria_meanvel=model->meanvel;
    290                         tria_epsvel=model->epsvel;
    291 
    292                         /*viscosity_overshoot*/
    293                         tria_viscosity_overshoot=model->viscosity_overshoot;
    294 
    295                         /*Create tria element using its constructor:*/
    296                         tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff);
    297 
    298                         /*Add tria element to elements dataset: */
    299                         elements->AddObject(tria);
    300 
    301                         /*Deal with material property card: */
    302                         matice_mid=i+1; //same as the material id from the geom2 elements.
    303                        
    304                         /*Average B over 3 grid elements: */
    305                         B_avg=0;
    306                         for(j=0;j<3;j++){
    307                                 B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
    308                         }
    309                         B_avg=B_avg/3;
    310                         matice_B=B_avg;
    311                         matice_n=(double)*(model->n+i);
    312                        
    313                         /*Create matice using its constructor:*/
    314                         matice= new Matice(matice_mid,matice_B,matice_n);
    315        
    316                         /*Add matice element to materials dataset: */
    317                         materials->AddObject(matice);
    318        
    319                         #ifdef _PARALLEL_
    320                         /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    321                          *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    322                          into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    323                          will hold which grids belong to this partition*/
    324                         my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
    325                         my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
    326                         my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
    327                         #endif
     244                        if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
     245                               
     246                                /*ids: */
     247                                tria_id=i+1; //matlab indexing.
     248                                tria_mid=i+1; //refers to the corresponding material property card
     249                                tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
     250
     251                                /*vertices offsets: */
     252                                tria_g[0]=(int)*(model->elements+elements_width*i+0);
     253                                tria_g[1]=(int)*(model->elements+elements_width*i+1);
     254                                tria_g[2]=(int)*(model->elements+elements_width*i+2);
     255
     256                                /*thickness,surface and bed:*/
     257                                tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
     258                                tria_h[1]=*(model->thickness+  ((int)*(model->elements+elements_width*i+1)-1));
     259                                tria_h[2]=*(model->thickness+  ((int)*(model->elements+elements_width*i+2)-1)) ;
     260
     261                                tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1));
     262                                tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1));
     263                                tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1));
     264
     265                                tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1));
     266                                tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1));
     267                                tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1));
     268
     269                                /*basal drag:*/
     270                                tria_friction_type=(int)model->drag_type;
     271
     272                                tria_k[0]=*(model->drag+        ((int)*(model->elements+elements_width*i+0)-1));
     273                                tria_k[1]=*(model->drag+        ((int)*(model->elements+elements_width*i+1)-1));
     274                                tria_k[2]=*(model->drag+        ((int)*(model->elements+elements_width*i+2)-1));
     275                               
     276                                tria_p=model->p[i];
     277                                tria_q=model->q[i];
     278
     279                                /*meling and accumulation*/
     280                                tria_melting[0]=*(model->melting+        ((int)*(model->elements+elements_width*i+0)-1));
     281                                tria_melting[1]=*(model->melting+        ((int)*(model->elements+elements_width*i+1)-1));
     282                                tria_melting[2]=*(model->melting+        ((int)*(model->elements+elements_width*i+2)-1));
     283
     284                                tria_accumulation[0]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+0)-1));
     285                                tria_accumulation[1]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+1)-1));
     286                                tria_accumulation[2]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+2)-1));
     287
     288                                /*element on iceshelf?:*/
     289                                tria_shelf=(int)*(model->elementoniceshelf+i);
     290
     291                                tria_meanvel=model->meanvel;
     292                                tria_epsvel=model->epsvel;
     293
     294                                /*viscosity_overshoot*/
     295                                tria_viscosity_overshoot=model->viscosity_overshoot;
     296
     297                                /*Create tria element using its constructor:*/
     298                                tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff);
     299
     300                                /*Add tria element to elements dataset: */
     301                                elements->AddObject(tria);
     302
     303                                /*Deal with material property card: */
     304                                matice_mid=i+1; //same as the material id from the geom2 elements.
     305                               
     306                                /*Average B over 3 grid elements: */
     307                                B_avg=0;
     308                                for(j=0;j<3;j++){
     309                                        B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
     310                                }
     311                                B_avg=B_avg/3;
     312                                matice_B=B_avg;
     313                                matice_n=(double)*(model->n+i);
     314                               
     315                                /*Create matice using its constructor:*/
     316                                matice= new Matice(matice_mid,matice_B,matice_n);
     317               
     318                                /*Add matice element to materials dataset: */
     319                                materials->AddObject(matice);
     320               
     321                                #ifdef _PARALLEL_
     322                                /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
     323                                 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
     324                                 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
     325                                 will hold which grids belong to this partition*/
     326                                my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
     327                                my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
     328                                my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     329                                #endif
     330                        } //if(!HutterEnum)
    328331
    329332                #ifdef _PARALLEL_
    … …  
    345348                xfree((void**)&model->B);
    346349                xfree((void**)&model->n);
     350                xfree((void**)&model->elements_type);
    347351
    348352        }
    … …  
    372376                #endif
    373377
     378                        if (*(model->elements_type+2*i+0)==MacAyealEnum() | *(model->elements_type+2*i+0)==PattynEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
    374379                       
    375                         /*name and id: */
    376                         penta_id=i+1; //matlab indexing.
    377                         penta_mid=i+1; //refers to the corresponding material property card
    378                         penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
    379 
    380                         /*vertices,thickness,surface,bed and drag: */
    381                         for(j=0;j<6;j++){
    382                                 penta_g[j]=(int)*(model->elements+elements_width*i+j);
    383                                 penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1));
    384                                 penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1));
    385                                 penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1));
    386                                 penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1));
    387                                 penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
    388                                 penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
     380                                /*name and id: */
     381                                penta_id=i+1; //matlab indexing.
     382                                penta_mid=i+1; //refers to the corresponding material property card
     383                                penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
     384
     385                                /*vertices,thickness,surface,bed and drag: */
     386                                for(j=0;j<6;j++){
     387                                        penta_g[j]=(int)*(model->elements+elements_width*i+j);
     388                                        penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1));
     389                                        penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1));
     390                                        penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1));
     391                                        penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1));
     392                                        penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
     393                                        penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
     394                                }
     395
     396                                /*basal drag:*/
     397                                penta_friction_type=(int)model->drag_type;
     398               
     399                                penta_p=model->p[i];
     400                                penta_q=model->q[i];
     401
     402                                /*diverse: */
     403                                penta_shelf=(int)*(model->elementoniceshelf+i);
     404                                penta_onbed=(int)*(model->elementonbed+i);
     405                                penta_onsurface=(int)*(model->elementonsurface+i);
     406                                penta_meanvel=model->meanvel;
     407                                penta_epsvel=model->epsvel;
     408                               
     409                                /*viscosity_overshoot*/
     410                                penta_viscosity_overshoot=model->viscosity_overshoot;
     411
     412                                if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
     413                                        penta_collapse=1;
     414                                }
     415                                else{
     416                                        penta_collapse=0;
     417                                }
     418
     419
     420                                /*Create Penta using its constructor:*/
     421                                penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
     422                                                penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
     423                                                penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
     424                                                penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning);
     425
     426                                /*Add penta element to elements dataset: */
     427                                elements->AddObject(penta);
     428               
     429
     430                                /*Deal with material:*/
     431                                matice_mid=i+1; //same as the material id from the geom2 elements.
     432                                /*Average B over 6 element grids: */
     433                                B_avg=0;
     434                                for(j=0;j<6;j++){
     435                                        B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
     436                                }
     437                                B_avg=B_avg/6;
     438                                matice_B= B_avg;
     439                                matice_n=(double)*(model->n+i);
     440               
     441                                /*Create matice using its constructor:*/
     442                                matice= new Matice(matice_mid,matice_B,matice_n);
     443               
     444                                /*Add matice element to materials dataset: */
     445                                materials->AddObject(matice);
     446                               
     447                                #ifdef _PARALLEL_
     448                                /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
     449                                 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
     450                                 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
     451                                 will hold which grids belong to this partition*/
     452                                my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
     453                                my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
     454                                my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     455                                my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
     456                                my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
     457                                my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
     458                                #endif
    389459                        }
    390 
    391                         /*basal drag:*/
    392                         penta_friction_type=(int)model->drag_type;
    393        
    394                         penta_p=model->p[i];
    395                         penta_q=model->q[i];
    396 
    397                         /*diverse: */
    398                         penta_shelf=(int)*(model->elementoniceshelf+i);
    399                         penta_onbed=(int)*(model->elementonbed+i);
    400                         penta_onsurface=(int)*(model->elementonsurface+i);
    401                         penta_meanvel=model->meanvel;
    402                         penta_epsvel=model->epsvel;
    403                        
    404                         /*viscosity_overshoot*/
    405                         penta_viscosity_overshoot=model->viscosity_overshoot;
    406 
    407                         if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
    408                                 penta_collapse=1;
    409                         }
    410                         else{
    411                                 penta_collapse=0;
    412                         }
    413 
    414 
    415                         /*Create Penta using its constructor:*/
    416                         penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
    417                                         penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
    418                                         penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
    419                                         penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning);
    420 
    421                         /*Add penta element to elements dataset: */
    422                         elements->AddObject(penta);
    423        
    424 
    425                         /*Deal with material:*/
    426                         matice_mid=i+1; //same as the material id from the geom2 elements.
    427                         /*Average B over 6 element grids: */
    428                         B_avg=0;
    429                         for(j=0;j<6;j++){
    430                                 B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
    431                         }
    432                         B_avg=B_avg/6;
    433                         matice_B= B_avg;
    434                         matice_n=(double)*(model->n+i);
    435        
    436                         /*Create matice using its constructor:*/
    437                         matice= new Matice(matice_mid,matice_B,matice_n);
    438        
    439                         /*Add matice element to materials dataset: */
    440                         materials->AddObject(matice);
    441                        
    442                         #ifdef _PARALLEL_
    443                         /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
    444                          *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
    445                          into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
    446                          will hold which grids belong to this partition*/
    447                         my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
    448                         my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
    449                         my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
    450                         my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
    451                         my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
    452                         my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
    453                         #endif
    454460
    455461                #ifdef _PARALLEL_
    … …  
    572578        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
    573579
    574 
    575580        /*Get number of dofs per node: */
    576581        DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);

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

    r586 r653  
    5353        model->uppernodes=NULL;
    5454        model->gridonhutter=NULL;
     55        model->gridonmacayeal=NULL;
     56        model->gridonpattyn=NULL;
    5557       
    5658        model->vx_obs=NULL;
    … …  
    203205        xfree((void**)&model->elements_type);
    204206        xfree((void**)&model->gridonhutter);
     207        xfree((void**)&model->gridonmacayeal);
    205208        if (strcmp(model->meshtype,"3d")==0){
    206209                xfree((void**)&model->elements2d);
    207210                xfree((void**)&model->deadgrids);
    208211                xfree((void**)&model->uppernodes);
     212                xfree((void**)&model->gridonpattyn);
    209213        }
    210214        xfree((void**)&model->solverstring);

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

    r586 r653  
    4242        int     isstokes;
    4343        double* gridonhutter;
     44        double* gridonmacayeal;
     45        double* gridonpattyn;
    4446
    4547        /*results: */
Note: See TracChangeset for help on using the changeset viewer.