Changeset 6320


Ignore:
Timestamp:
10/15/10 11:55:16 (15 years ago)
Author:
seroussi
Message:

some cleaning in Penta

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r6271 r6320  
    652652void  Penta::ComputeBasalStress(Vec sigma_b){
    653653
    654         int i,j;
    655         int    doflist[NUMVERTICES];
    656         double xyz_list[NUMVERTICES][3];
    657         double xyz_list_tria[3][3];
    658 
    659         /*Parameters*/
    660         double rho_ice,gravity;
    661         double bed_normal[3];
    662         double bed;
    663         double basalforce[3];
    664         double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    665         double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    666         double stresstensor[6]={0.0};
    667         double viscosity;
    668         int analysis_type;
    669 
    670         int  dofv[3]={0,1,2};
    671         int  dofp[1]={3};
    672         double Jdet2d;
    673 
    674         /*Gauss*/
    675         int     ig;
     654        int         i,j,ig;
     655        int         dofv[3]={0,1,2};
     656        int         dofp[1]={3};
     657        int         analysis_type,approximation;
     658        int         doflist[NUMVERTICES];
     659        double      xyz_list[NUMVERTICES][3];
     660        double      xyz_list_tria[3][3];
     661        double      rho_ice,gravity,stokesreconditioning;
     662        double      pressure,viscosity,bed,Jdet2d;
     663        double      bed_normal[3];
     664        double      basalforce[3];
     665        double      epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     666        double      devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     667        double      stresstensor[6]={0.0};
     668        double      sigma_xx,sigma_yy,sigma_zz;
     669        double      sigma_xy,sigma_xz,sigma_yz;
     670        double      surface=0,value=0;
    676671        GaussPenta* gauss;
    677 
    678         /*Output*/
    679         double pressure;
    680         double sigma_xx,sigma_yy,sigma_zz;
    681         double sigma_xy,sigma_xz,sigma_yz;
    682         double surface=0;
    683         double value=0;
    684        
    685         /*flags: */
    686         int  approximation;
    687 
    688         /*parameters: */
    689         double stokesreconditioning;
    690672
    691673        /*retrive parameters: */
    … …  
    712694        /* Get node coordinates and dof list: */
    713695        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    714         for(i=0;i<3;i++){
    715                 for(j=0;j<3;j++){
    716                         xyz_list_tria[i][j]=xyz_list[i][j];
    717                 }
    718         }
     696        for(i=0;i<3;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    719697
    720698        /*Retrieve all inputs we will be needing: */
    … …  
    11141092
    11151093        int i;
     1094        bool    converged=true;
    11161095        Input** new_inputs=NULL;
    11171096        Input** old_inputs=NULL;
    1118         bool    converged=true;
    11191097
    11201098        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
    … …  
    11401118        /*Return output*/
    11411119        return converged;
    1142 
    11431120}
    11441121/*}}}*/
    … …  
    11891166void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
    11901167
    1191         /*Intermediaries*/
    11921168        int  step,i;
    1193 
     1169        double  xyz_list[NUMVERTICES][3];
     1170        double  Helem_list[NUMVERTICES];
     1171        double  zeros_list[NUMVERTICES]={0.0};
    11941172        Penta* penta=NULL;
    11951173        Input* original_input=NULL;
    … …  
    11991177        Input* total_thickness_input=NULL;
    12001178        Input* depth_averaged_input=NULL;
    1201 
    1202         double  xyz_list[NUMVERTICES][3];
    1203         double  Helem_list[NUMVERTICES];
    1204         double  zeros_list[NUMVERTICES]={0.0};
    12051179
    12061180        /*recover parameters: */
    … …  
    14081382        /*Assign output pointers:*/
    14091383        *pmaxvx=maxvx;
    1410 
    14111384}
    14121385/*}}}*/
    … …  
    14221395        /*Assign output pointers:*/
    14231396        *pmaxvy=maxvy;
    1424 
    14251397}
    14261398/*}}}*/
    … …  
    14361408        /*Assign output pointers:*/
    14371409        *pmaxvz=maxvz;
    1438 
    14391410}
    14401411/*}}}*/
    … …  
    14501421        /*Assign output pointers:*/
    14511422        *pminvel=minvel;
    1452 
    14531423}
    14541424/*}}}*/
    … …  
    14641434        /*Assign output pointers:*/
    14651435        *pminvx=minvx;
    1466 
    14671436}
    14681437/*}}}*/
    … …  
    14781447        /*Assign output pointers:*/
    14791448        *pminvy=minvy;
    1480 
    14811449}
    14821450/*}}}*/
    … …  
    14921460        /*Assign output pointers:*/
    14931461        *pminvz=minvz;
    1494 
    14951462}
    14961463/*}}}*/
    … …  
    14981465double  Penta::TimeAdapt(void){
    14991466
    1500         /*intermediary: */
    1501         double C;
    1502         double maxabsvx;
    1503         double maxabsvy;
    1504         double maxabsvz;
    1505         double maxx,minx;
    1506         double maxy,miny;
    1507         double maxz,minz;
    1508         double dx,dy,dz;
     1467        int    i;
     1468        double C,dx,dy,dz,dt;
     1469        double maxabsvx,maxabsvy,maxabsvz;
     1470        double maxx,minx,maxy,miny,maxz,minz;
    15091471        double xyz_list[NUMVERTICES][3];
    1510         int    i;
    1511 
    1512         /*output: */
    1513         double dt;
    15141472
    15151473        /*get CFL coefficient:*/
    … …  
    15471505
    15481506        return dt;
    1549 
    15501507}
    15511508/*FUNCTION Penta::ThicknessAbsMisfit {{{1*/
    15521509double Penta::ThicknessAbsMisfit(bool process_units){
    15531510
     1511        int    approximation;
    15541512        double J;
    1555         Tria* tria=NULL;
    1556 
    1557         /*inputs: */
    1558         int  approximation;
     1513        Tria*  tria=NULL;
    15591514
    15601515        /*retrieve inputs :*/
    … …  
    15741529double Penta::SurfaceAbsVelMisfit(bool process_units){
    15751530
     1531        int    approximation;
    15761532        double J;
    1577         Tria* tria=NULL;
    1578 
    1579         /*inputs: */
    1580         int  approximation;
     1533        Tria*  tria=NULL;
    15811534
    15821535        /*retrieve inputs :*/
    … …  
    16131566double Penta::SurfaceRelVelMisfit(bool process_units){
    16141567
     1568        int    approximation;
    16151569        double J;
    1616         Tria* tria=NULL;
    1617 
    1618         /*inputs: */
    1619         int  approximation;
     1570        Tria*  tria=NULL;
    16201571
    16211572        /*retrieve inputs :*/
    … …  
    16521603double Penta::SurfaceLogVelMisfit(bool process_units){
    16531604
     1605        int    approximation;
    16541606        double J;
    1655         Tria* tria=NULL;
    1656 
    1657         /*inputs: */
    1658         int  approximation;
     1607        Tria*  tria=NULL;
    16591608
    16601609        /*retrieve inputs :*/
    … …  
    17301679double Penta::SurfaceAverageVelMisfit(bool process_units){
    17311680
     1681        int    approximation;
    17321682        double J;
    1733         Tria* tria=NULL;
    1734 
    1735         /*inputs: */
    1736         int  approximation;
     1683        Tria*  tria=NULL;
    17371684
    17381685        /*retrieve inputs :*/
    … …  
    17691716void  Penta::PatchFill(int* pcount, Patch* patch){
    17701717
    1771         int i;
    1772         int count;
     1718        int i,count;
    17731719        int vertices_ids[6];
    1774 
    17751720
    17761721        /*recover pointer: */
    … …  
    17991744
    18001745        int     i;
    1801        
    1802         /*output: */
    1803         int     numrows     = 0;
    1804         int     numnodes    = 0;
     1746        int     numrows  = 0;
     1747        int     numnodes = 0;
    18051748
    18061749        /*Go through all the results objects, and update the counters: */
    … …  
    18171760        *pnumvertices=NUMVERTICES;
    18181761        *pnumnodes=numnodes;
    1819        
    18201762}
    18211763/*}}}*/
    … …  
    18291771                elementresult->ProcessUnits(this->parameters);
    18301772        }
    1831 
    18321773}
    18331774/*}}}*/
    … …  
    18511792double Penta::SurfaceArea(void){
    18521793
     1794        int    approximation;
    18531795        double S;
    1854         Tria* tria=NULL;
    1855 
    1856         /*inputs: */
    1857         int  approximation;
     1796        Tria*  tria=NULL;
    18581797
    18591798        /*retrieve inputs :*/
    … …  
    19141853
    19151854        /*Recover vertices ids needed to initialize inputs*/
    1916         for(i=0;i<6;i++){
    1917                 penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
    1918         }
     1855        for(i=0;i<6;i++) penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
    19191856
    19201857        /*Recover nodes ids needed to initialize the node hook.*/
    … …  
    20031940                        break;
    20041941        }
    2005 
    20061942}
    20071943/*}}}*/
    … …  
    20101946
    20111947        /*Intermediaries*/
    2012         int i;
     1948        int    i;
    20131949        double rho_ice,rho_water;
    20141950
    … …  
    20391975                this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
    20401976        }
    2041 
    20421977}
    20431978/*}}}*/
    … …  
    20481983
    20491984        int i;
    2050         double v13[3];
    2051         double v23[3];
     1985        double v13[3],v23[3];
    20521986        double normal[3];
    20531987        double normal_norm;
    … …  
    20611995        normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
    20621996        normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
    2063 
    20641997        normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
    20651998
    … …  
    20682001        *(bed_normal+1)=-normal[1]/normal_norm;
    20692002        *(bed_normal+2)=-normal[2]/normal_norm;
    2070 
    20712003}
    20722004/*}}}*/
    … …  
    39393871void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
    39403872
    3941         int i,j;
    3942         int numberofdofs=0;
    3943         int count=0;
    3944 
    3945         /*output: */
     3873        int  i,j,count=0;
     3874        int  numberofdofs=0;
    39463875        int* doflist=NULL;
    39473876
    39483877        /*First, figure out size of doflist: */
    3949         for(i=0;i<6;i++){
    3950                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    3951         }
     3878        for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    39523879
    39533880        /*Allocate: */
    … …  
    39633890        /*Assign output pointers:*/
    39643891        *pdoflist=doflist;
    3965 
    39663892}
    39673893/*}}}*/
    … …  
    39703896       
    39713897        int i;
    3972 
    3973         for(i=0;i<6;i++){
    3974                 doflist[i]=nodes[i]->GetDofList1();
    3975         }
     3898        for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1();
    39763899
    39773900}
    … …  
    39903913        /*return PentaRef field*/
    39913914        return this->element_type;
    3992 
    39933915}
    39943916/*}}}*/
    … …  
    41014023void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
    41024024
    4103         int i;
    4104 
    4105         const int    numdofpervertex=2;
    4106         const int    numdof=numdofpervertex*NUMVERTICES;
     4025        const int    numdof=NDOF2*NUMVERTICES;
     4026
     4027        int          i;
     4028        int          approximation;
     4029        int*         doflist=NULL;
     4030        double       vx,vy;
     4031        double       values[numdof];
    41074032        GaussPenta*  gauss;
    4108         int*         doflist=NULL;
    4109         double       values[numdof];
    4110         double       vx;
    4111         double       vy;
    4112         int          approximation;
    41134033
    41144034        /*Get approximation enum and dof list: */
    … …  
    41304050                vx_input->GetParameterValue(&vx,gauss);
    41314051                vy_input->GetParameterValue(&vy,gauss);
    4132                 values[i*numdofpervertex+0]=vx;
    4133                 values[i*numdofpervertex+1]=vy;
     4052                values[i*NDOF2+0]=vx;
     4053                values[i*NDOF2+1]=vy;
    41344054        }
    41354055
    … …  
    41454065void  Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){
    41464066
    4147         int i;
    4148 
    4149         const int    numdofpervertex=2;
    4150         const int    numdof=numdofpervertex*NUMVERTICES;
     4067        const int    numdof=NDOF2*NUMVERTICES;
     4068
     4069        int          i;
     4070        int*         doflist=NULL;
     4071        double       vx,vy;
     4072        double       values[numdof];
    41514073        GaussPenta*  gauss=NULL;
    4152         int*         doflist=NULL;
    4153         double       values[numdof];
    4154         double       vx;
    4155         double       vy;
    41564074
    41574075        /*Get dof list: */
    … …  
    41644082        gauss=new GaussPenta();
    41654083        for(i=0;i<NUMVERTICES;i++){
    4166 
    41674084                /*Recover vx and vy*/
    41684085                gauss->GaussVertex(i);
    41694086                vx_input->GetParameterValue(&vx,gauss);
    41704087                vy_input->GetParameterValue(&vy,gauss);
    4171                 values[i*numdofpervertex+0]=vx;
    4172                 values[i*numdofpervertex+1]=vy;
     4088                values[i*NDOF2+0]=vx;
     4089                values[i*NDOF2+1]=vy;
    41734090        }
    41744091
    … …  
    41844101void  Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){
    41854102
    4186         int i;
    4187 
    4188         const int    numdofpervertex=1;
    4189         const int    numdof=numdofpervertex*NUMVERTICES;
     4103        const int    numdof=NDOF1*NUMVERTICES;
     4104
     4105        int          i;
     4106        int*         doflist=NULL;
     4107        double       vz;
     4108        double       values[numdof];
    41904109        GaussPenta*  gauss=NULL;
    4191         int*         doflist=NULL;
    4192         double       values[numdof];
    4193         double       vz;
    41944110
    41954111        /*Get dof list: */
    … …  
    42014117        gauss=new GaussPenta();
    42024118        for(i=0;i<NUMVERTICES;i++){
    4203 
    42044119                /*Recover vz */
    42054120                gauss->GaussVertex(i);
    … …  
    42194134void  Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){
    42204135
    4221         int i;
    4222 
    4223         const int    numdofpervertex=4;
    4224         const int    numdof=numdofpervertex*NUMVERTICES;
    4225         GaussPenta  *gauss;
     4136        const int    numdof=NDOF4*NUMVERTICES;
     4137
     4138        int          i;
    42264139        int*         doflist=NULL;
    4227         double       values[numdof];
    42284140        double       vx,vy,vz,p;
    42294141        double       stokesreconditioning;
     4142        double       values[numdof];
     4143        GaussPenta   *gauss;
    42304144
    42314145        /*Get dof list: */
    … …  
    42484162                vz_input->GetParameterValue(&vz,gauss);
    42494163                p_input ->GetParameterValue(&p ,gauss);
    4250                 values[i*numdofpervertex+0]=vx;
    4251                 values[i*numdofpervertex+1]=vy;
    4252                 values[i*numdofpervertex+2]=vz;
    4253                 values[i*numdofpervertex+3]=p/stokesreconditioning;
     4164                values[i*NDOF4+0]=vx;
     4165                values[i*NDOF4+1]=vy;
     4166                values[i*NDOF4+2]=vz;
     4167                values[i*NDOF4+3]=p/stokesreconditioning;
    42544168        }
    42554169
    … …  
    42654179void  Penta::GetSolutionFromInputsThermal(Vec solution){
    42664180
    4267         int i;
    4268 
    4269         const int    numdofpervertex=1;
    4270         const int    numdof=numdofpervertex*NUMVERTICES;
    4271         GaussPenta  *gauss=NULL;
     4181        const int    numdof=NDOF1*NUMVERTICES;
     4182
     4183        int          i;
    42724184        int*         doflist=NULL;
    42734185        double       values[numdof];
    42744186        double       temp;
    4275 
     4187        GaussPenta   *gauss=NULL;
    42764188
    42774189        /*Get dof list: */
    … …  
    42794191        Input* t_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(t_input);
    42804192
    4281         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4282         /*P1 element only for now*/
    42834193        gauss=new GaussPenta();
    42844194        for(i=0;i<NUMVERTICES;i++){
    4285 
    4286                 /*Recover vz */
     4195                /*Recover temperature*/
    42874196                gauss->GaussVertex(i);
    42884197                t_input->GetParameterValue(&temp,gauss);
    … …  
    43274236
    43284237        int i;
    4329 
    43304238        double epsilonvx[5];
    43314239        double epsilonvy[5];
    … …  
    43424250        /*Sum all contributions*/
    43434251        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    4344 
    43454252}
    43464253/*}}}*/
    … …  
    43534260
    43544261        int i;
    4355 
    43564262        double epsilonvx[6];
    43574263        double epsilonvy[6];
    … …  
    43704276        /*Sum all contributions*/
    43714277        for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
    4372 
    43734278}
    43744279/*}}}*/
    … …  
    43824287        penta=this;
    43834288        for(;;){
    4384 
    43854289                /*Stop if we have reached the surface, else, take lower penta*/
    43864290                if (penta->IsOnBed()) break;
    … …  
    43994303
    44004304        Penta* upper_penta=NULL;
     4305
    44014306        upper_penta=(Penta*)neighbors[0]; //first one (0) under, second one (1) above
     4307       
    44024308        return upper_penta;
    4403 
    44044309}
    44054310/*}}}*/
    … …  
    44084313
    44094314        Penta* upper_penta=NULL;
     4315       
    44104316        upper_penta=(Penta*)neighbors[1]; //first one under, second one above
     4317
    44114318        return upper_penta;
    4412 
    44134319}
    44144320/*}}}*/
    … …  
    44174323
    44184324        int    i;
     4325        double z;
    44194326        double xyz_list[NUMVERTICES][3];
    44204327        double z_list[NUMVERTICES];
    4421         double z;
    44224328
    44234329        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4424         for(i=0;i<NUMVERTICES;i++){
    4425                 z_list[i]=xyz_list[i][2];
    4426         }
     4330        for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
    44274331        PentaRef::GetParameterValue(&z,z_list,gauss);
    44284332
    … …  
    44334337void  Penta::GradjB(Vec gradient){
    44344338
    4435         int i;
    4436         Tria* tria=NULL;
     4339        int              i;
     4340        int              approximation;
     4341        Tria*            tria           =NULL;
    44374342        TriaVertexInput* triavertexinput=NULL;
    4438 
    4439         /*retrieve inputs :*/
    4440         int  approximation;
    4441         inputs->GetParameterValue(&approximation,ApproximationEnum);
    44424343
    44434344        /*If on water, skip: */
    44444345        if(IsOnWater())return;
     4346        inputs->GetParameterValue(&approximation,ApproximationEnum);
    44454347
    44464348        if (approximation==MacAyealApproximationEnum){
    … …  
    44764378                this->matice->inputs->DeleteInput(RheologyBbarEnum);
    44774379        }
    4478 
    44794380}
    44804381/*}}}*/
    … …  
    44824383void  Penta::GradjDrag(Vec gradient){
    44834384
    4484         int i;
    4485         Tria* tria=NULL;
     4385        int              i,approximation;
     4386        double           temp_gradient[6]={0,0,0,0,0,0};
     4387        Tria*            tria=NULL;
    44864388        TriaVertexInput* triavertexinput=NULL;
    4487         double temp_gradient[6]={0,0,0,0,0,0};
    44884389
    44894390        /*retrieve inputs :*/
    4490         int approximation;
    44914391        inputs->GetParameterValue(&approximation,ApproximationEnum);
    44924392
    4493         /*If on water, skip: */
    4494         if(IsOnWater())return;
    4495 
    4496         /*If on shelf, skip: */
    4497         if(IsOnShelf())return;
    4498 
    4499         /*Bail out if this element does not touch the bedrock: */
    4500         if (!IsOnBed()) return;
     4393        /*If on water, on shelf or not on bed, skip: */
     4394        if(IsOnWater()|| IsOnShelf() || !IsOnBed())return;
    45014395
    45024396        if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum){
    4503 
    45044397                /*MacAyeal or Pattyn*/
    45054398                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    … …  
    45084401        }
    45094402        else if (approximation==StokesApproximationEnum){
    4510 
    45114403                /*Stokes*/
    45124404                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    … …  
    45184410        }
    45194411        else ISSMERROR("approximation %s not supported yet",EnumToString(approximation));
    4520 
    4521 
    45224412}
    45234413/*}}}*/
    … …  
    45304420        /*Are we on the base, not on the surface?:*/
    45314421        if(IsOnBed()){
    4532 
    45334422                /*OK, we are on bed. we will follow the steps:
    45344423                 * 1: find input and extrude it.
    … …  
    45694458                        /*Stop if we have reached the surface*/
    45704459                        if (penta->IsOnSurface()) break;
    4571 
    45724460                }
    45734461        }
    … …  
    45794467void  Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
    45804468
    4581         /*Inputs*/
    45824469        int  approximation;
    45834470
    … …  
    46124499        }
    46134500}
    4614 
    46154501/*}}}*/
    46164502/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/
    46174503void  Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){
    46184504
    4619         int i;
    4620 
    4621         const int    numdofpervertex=2;
    4622         const int    numdof=numdofpervertex*NUMVERTICES;
    4623         int*         doflist=NULL;
    4624         double       values[numdof];
    4625         double       vx[NUMVERTICES];
    4626         double       vy[NUMVERTICES];
    4627         double       vz[NUMVERTICES];
    4628         double       vel[NUMVERTICES];
    4629         int          dummy;
    4630         double       pressure[NUMVERTICES];
    4631         double       surface[NUMVERTICES];
    4632         double       rho_ice,g;
    4633         double       xyz_list[NUMVERTICES][3];
    4634 
    4635         double *vz_ptr          = NULL;
    4636         Penta  *penta          = NULL;
    4637 
    4638         /*OK, we are on bed. Now recover results*/
     4505        const int    numdof=NDOF2*NUMVERTICES;
     4506
     4507        int     i,dummy;
     4508        double  rho_ice,g;
     4509        double  values[numdof];
     4510        double  vx[NUMVERTICES];
     4511        double  vy[NUMVERTICES];
     4512        double  vz[NUMVERTICES];
     4513        double  vel[NUMVERTICES];
     4514        double  pressure[NUMVERTICES];
     4515        double  surface[NUMVERTICES];
     4516        double  xyz_list[NUMVERTICES][3];
     4517        int    *doflist = NULL;
     4518        double *vz_ptr  = NULL;
     4519        Penta  *penta   = NULL;
    46394520
    46404521        /*Get dof list: */
    … …  
    46424523
    46434524        /*Use the dof list to index into the solution vector: */
    4644         for(i=0;i<numdof;i++){
    4645                 values[i]=solution[doflist[i]];
    4646         }
     4525        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    46474526
    46484527        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    46494528        for(i=0;i<3;i++){
    4650                 vx[i]  =values[i*numdofpervertex+0];
    4651                 vy[i]  =values[i*numdofpervertex+1];
     4529                vx[i]  =values[i*NDOF2+0];
     4530                vy[i]  =values[i*NDOF2+1];
    46524531                vx[i+3]=vx[i];
    46534532                vy[i+3]=vy[i];
    … …  
    47054584void  Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){
    47064585
    4707         int i;
    4708 
    4709         const int    numdofpervertex=2;
    4710         const int    numdof=numdofpervertex*NUMVERTICES;
    4711         const int    numdof2d=numdofpervertex*NUMVERTICES2D;
    4712         int*         doflistp=NULL;
    4713         int*         doflistm=NULL;
    4714         double       macayeal_values[numdof];
    4715         double       pattyn_values[numdof];
    4716         double       vx[NUMVERTICES];
    4717         double       vy[NUMVERTICES];
    4718         double       vz[NUMVERTICES];
    4719         double       vel[NUMVERTICES];
    4720         int          dummy;
    4721         double       pressure[NUMVERTICES];
    4722         double       surface[NUMVERTICES];
    4723         double       rho_ice,g;
    4724         double       xyz_list[NUMVERTICES][3];
    4725 
    4726         double *vz_ptr         = NULL;
    4727         Penta  *penta          = NULL;
     4586        const int    numdof=NDOF2*NUMVERTICES;
     4587        const int    numdof2d=NDOF2*NUMVERTICES2D;
     4588
     4589        int     i,dummy;
     4590        double  rho_ice,g;
     4591        double  macayeal_values[numdof];
     4592        double  pattyn_values[numdof];
     4593        double  vx[NUMVERTICES];
     4594        double  vy[NUMVERTICES];
     4595        double  vz[NUMVERTICES];
     4596        double  vel[NUMVERTICES];
     4597        double  pressure[NUMVERTICES];
     4598        double  surface[NUMVERTICES];
     4599        double  xyz_list[NUMVERTICES][3];
     4600        int*    doflistp = NULL;
     4601        int*    doflistm = NULL;
     4602        double  *vz_ptr  = NULL;
     4603        Penta   *penta   = NULL;
    47284604
    47294605        /*OK, we have to add results of this element for pattyn
    … …  
    47504626        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    47514627        for(i=0;i<NUMVERTICES;i++){
    4752                 vx[i]=macayeal_values[i*numdofpervertex+0]+pattyn_values[i*numdofpervertex+0];
    4753                 vy[i]=macayeal_values[i*numdofpervertex+1]+pattyn_values[i*numdofpervertex+1];
     4628                vx[i]=macayeal_values[i*NDOF2+0]+pattyn_values[i*NDOF2+0];
     4629                vy[i]=macayeal_values[i*NDOF2+1]+pattyn_values[i*NDOF2+1];
    47544630        }
    47554631
    … …  
    47684644
    47694645        /*Now Compute vel*/
    4770         for(i=0;i<NUMVERTICES;i++) {
    4771                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    4772         }
     4646        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    47734647
    47744648        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    … …  
    47994673void  Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(double* solution){
    48004674
    4801         int i;
    4802 
    48034675        const int    numdofm=NDOF2*NUMVERTICES;
    48044676        const int    numdofs=NDOF4*NUMVERTICES;
    48054677        const int    numdof2d=NDOF2*NUMVERTICES2D;
    4806         int*         doflistm=NULL;
    4807         int*         doflists=NULL;
    4808         double       macayeal_values[numdofm];
    4809         double       stokes_values[numdofs];
    4810         double       vx[NUMVERTICES];
    4811         double       vy[NUMVERTICES];
    4812         double       vz[NUMVERTICES];
    4813         double       vzmacayeal[NUMVERTICES];
    4814         double       vzstokes[NUMVERTICES];
    4815         double       vel[NUMVERTICES];
    4816         int          dummy;
    4817         double       pressure[NUMVERTICES];
    4818         double       xyz_list[NUMVERTICES][3];
    4819         double       stokesreconditioning;
    4820 
    4821         double *vzmacayeal_ptr        = NULL;
    4822         Penta  *penta          = NULL;
     4678
     4679        int     i,dummy;
     4680        double  stokesreconditioning;
     4681        double  macayeal_values[numdofm];
     4682        double  stokes_values[numdofs];
     4683        double  vx[NUMVERTICES];
     4684        double  vy[NUMVERTICES];
     4685        double  vz[NUMVERTICES];
     4686        double  vzmacayeal[NUMVERTICES];
     4687        double  vzstokes[NUMVERTICES];
     4688        double  vel[NUMVERTICES];
     4689        double  pressure[NUMVERTICES];
     4690        double  xyz_list[NUMVERTICES][3];
     4691        int*    doflistm        = NULL;
     4692        int*    doflists        = NULL;
     4693        double  *vzmacayeal_ptr = NULL;
     4694        Penta   *penta          = NULL;
    48234695
    48244696        /*OK, we have to add results of this element for macayeal
    … …  
    48934765void  Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){
    48944766       
    4895         int i;
    4896 
    4897         const int    numdofpervertex=2;
    4898         const int    numdof=numdofpervertex*NUMVERTICES;
    4899         int*         doflist=NULL;
    4900         double       values[numdof];
    4901         double       vx[NUMVERTICES];
    4902         double       vy[NUMVERTICES];
    4903         double       vz[NUMVERTICES];
    4904         double       vel[NUMVERTICES];
    4905         int          dummy;
    4906         double       pressure[NUMVERTICES];
    4907         double       surface[NUMVERTICES];
    4908         double       rho_ice,g;
    4909         double       xyz_list[NUMVERTICES][3];
    4910         double *vz_ptr          = NULL;
     4767        const int    numdof=NDOF2*NUMVERTICES;
     4768
     4769        int    i,dummy;
     4770        double rho_ice,g;
     4771        double values[numdof];
     4772        double vx[NUMVERTICES];
     4773        double vy[NUMVERTICES];
     4774        double vz[NUMVERTICES];
     4775        double vel[NUMVERTICES];
     4776        double pressure[NUMVERTICES];
     4777        double surface[NUMVERTICES];
     4778        double xyz_list[NUMVERTICES][3];
     4779        int*   doflist = NULL;
     4780        double *vz_ptr = NULL;
    49114781
    49124782        /*Get dof list: */
    … …  
    49174787
    49184788        /*Use the dof list to index into the solution vector: */
    4919         for(i=0;i<numdof;i++){
    4920                 values[i]=solution[doflist[i]];
    4921         }
     4789        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    49224790
    49234791        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    49244792        for(i=0;i<NUMVERTICES;i++){
    4925                 vx[i]=values[i*numdofpervertex+0];
    4926                 vy[i]=values[i*numdofpervertex+1];
     4793                vx[i]=values[i*NDOF2+0];
     4794                vy[i]=values[i*NDOF2+1];
    49274795        }
    49284796
    … …  
    49654833        xfree((void**)&doflist);
    49664834}
    4967 
    49684835/*}}}*/
    49694836/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{1*/
    49704837void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){
    49714838
    4972         int i;
    4973 
    4974         const int    numdofpervertexp=2;
    4975         const int    numdofpervertexs=4;
    4976         const int    numdofp=numdofpervertexp*NUMVERTICES;
    4977         const int    numdofs=numdofpervertexs*NUMVERTICES;
    4978         int*         doflistp=NULL;
    4979         int*         doflists=NULL;
    4980         double       pattyn_values[numdofp];
    4981         double       stokes_values[numdofs];
    4982         double       vx[NUMVERTICES];
    4983         double       vy[NUMVERTICES];
    4984         double       vz[NUMVERTICES];
    4985         double       vzpattyn[NUMVERTICES];
    4986         double       vzstokes[NUMVERTICES];
    4987         double       vel[NUMVERTICES];
    4988         int          dummy;
    4989         double       pressure[NUMVERTICES];
    4990         double       xyz_list[NUMVERTICES][3];
    4991         double       stokesreconditioning;
    4992 
    4993         double *vzpattyn_ptr         = NULL;
    4994         Penta  *penta          = NULL;
     4839        const int    numdofp=NDOF2*NUMVERTICES;
     4840        const int    numdofs=NDOF4*NUMVERTICES;
     4841
     4842        int    i,dummy;
     4843        double pattyn_values[numdofp];
     4844        double stokes_values[numdofs];
     4845        double vx[NUMVERTICES];
     4846        double vy[NUMVERTICES];
     4847        double vz[NUMVERTICES];
     4848        double vzpattyn[NUMVERTICES];
     4849        double vzstokes[NUMVERTICES];
     4850        double vel[NUMVERTICES];
     4851        double pressure[NUMVERTICES];
     4852        double xyz_list[NUMVERTICES][3];
     4853        double stokesreconditioning;
     4854        int*   doflistp      = NULL;
     4855        int*   doflists      = NULL;
     4856        double *vzpattyn_ptr = NULL;
     4857        Penta  *penta        = NULL;
    49954858
    49964859        /*OK, we have to add results of this element for pattyn
    … …  
    50074870
    50084871        /*Use the dof list to index into the solution vector: */
    5009         for(i=0;i<numdofp;i++){
    5010                 pattyn_values[i]=solution[doflistp[i]];
    5011         }
    5012         for(i=0;i<numdofs;i++){
    5013                 stokes_values[i]=solution[doflists[i]];
    5014         }
     4872        for(i=0;i<numdofp;i++) pattyn_values[i]=solution[doflistp[i]];
     4873        for(i=0;i<numdofs;i++) stokes_values[i]=solution[doflists[i]];
    50154874
    50164875        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    50174876        for(i=0;i<NUMVERTICES;i++){
    5018                 vx[i]=stokes_values[i*numdofpervertexs+0]+pattyn_values[i*numdofpervertexp+0];
    5019                 vy[i]=stokes_values[i*numdofpervertexs+1]+pattyn_values[i*numdofpervertexp+1];
    5020                 vzstokes[i]=stokes_values[i*numdofpervertexs+2];
    5021                 pressure[i]=stokes_values[i*numdofpervertexs+3]*stokesreconditioning;
     4877                vx[i]=stokes_values[i*NDOF4+0]+pattyn_values[i*NDOF2+0];
     4878                vy[i]=stokes_values[i*NDOF4+1]+pattyn_values[i*NDOF2+1];
     4879                vzstokes[i]=stokes_values[i*NDOF4+2];
     4880                pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
    50224881        }
    50234882
    … …  
    50644923void  Penta::InputUpdateFromSolutionDiagnosticHutter(double* solution){
    50654924       
    5066         int i;
    5067 
    5068         const int    numdofpervertex=2;
    5069         const int    numdof=numdofpervertex*NUMVERTICES;
    5070         int*         doflist=NULL;
    5071         double       values[numdof];
    5072         double       vx[NUMVERTICES];
    5073         double       vy[NUMVERTICES];
    5074         double       vz[NUMVERTICES];
    5075         double       vel[NUMVERTICES];
    5076         int          dummy;
    5077         double       pressure[NUMVERTICES];
    5078         double       surface[NUMVERTICES];
    5079         double       rho_ice,g;
    5080         double       xyz_list[NUMVERTICES][3];
    5081         double*      vz_ptr=NULL;
     4925        const int    numdof=NDOF2*NUMVERTICES;
     4926
     4927        int     i,dummy;
     4928        double  rho_ice,g;
     4929        double  values[numdof];
     4930        double  vx[NUMVERTICES];
     4931        double  vy[NUMVERTICES];
     4932        double  vz[NUMVERTICES];
     4933        double  vel[NUMVERTICES];
     4934        double  pressure[NUMVERTICES];
     4935        double  surface[NUMVERTICES];
     4936        double  xyz_list[NUMVERTICES][3];
     4937        int*    doflist = NULL;
     4938        double* vz_ptr  = NULL;
    50824939
    50834940        /*Get dof list: */
    … …  
    50884945
    50894946        /*Use the dof list to index into the solution vector: */
    5090         for(i=0;i<numdof;i++){
    5091                 values[i]=solution[doflist[i]];
    5092         }
     4947        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    50934948
    50944949        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    50954950        for(i=0;i<NUMVERTICES;i++){
    5096                 vx[i]=values[i*numdofpervertex+0];
    5097                 vy[i]=values[i*numdofpervertex+1];
     4951                vx[i]=values[i*NDOF2+0];
     4952                vy[i]=values[i*NDOF2+1];
    50984953        }
    50994954
    … …  
    51364991        xfree((void**)&doflist);
    51374992}
    5138 
    51394993/*}}}*/
    51404994/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticVert {{{1*/
    51414995void  Penta::InputUpdateFromSolutionDiagnosticVert(double* solution){
    51424996
    5143         int i;
    5144 
    5145         const int    numdofpervertex=1;
    5146         const int    numdof=numdofpervertex*NUMVERTICES;
    5147         int*         doflist=NULL;
    5148         double       values[numdof];
    5149         double       vx[NUMVERTICES];
    5150         double       vy[NUMVERTICES];
    5151         double       vz[NUMVERTICES];
    5152         double       vzmacayeal[NUMVERTICES];
    5153         double       vzpattyn[NUMVERTICES];
    5154         double       vzstokes[NUMVERTICES];
    5155         double       vel[NUMVERTICES];
    5156         int          dummy;
    5157         double       pressure[NUMVERTICES];
    5158         double       surface[NUMVERTICES];
    5159         double       rho_ice,g;
    5160         double       xyz_list[NUMVERTICES][3];
    5161 
    5162         double*      vx_ptr=NULL;
    5163         double*      vy_ptr=NULL;
    5164         double*      vzstokes_ptr=NULL;
    5165 
    5166         int         approximation;
     4997        const int numdof=NDOF1*NUMVERTICES;
     4998       
     4999        int      i,dummy;
     5000        int      approximation;
     5001        double   rho_ice,g;
     5002        double   values[numdof];
     5003        double   vx[NUMVERTICES];
     5004        double   vy[NUMVERTICES];
     5005        double   vz[NUMVERTICES];
     5006        double   vzmacayeal[NUMVERTICES];
     5007        double   vzpattyn[NUMVERTICES];
     5008        double   vzstokes[NUMVERTICES];
     5009        double   vel[NUMVERTICES];
     5010        double   pressure[NUMVERTICES];
     5011        double   surface[NUMVERTICES];
     5012        double   xyz_list[NUMVERTICES][3];
     5013        int*     doflist      = NULL;
     5014        double*  vx_ptr       = NULL;
     5015        double*  vy_ptr       = NULL;
     5016        double*  vzstokes_ptr = NULL;
     5017
    51675018
    51685019        /*Get the approximation and do nothing if the element in Stokes or None*/
    … …  
    51725023        }
    51735024
    5174         /*Get dof list: */
     5025        /*Get dof list and vertices coordinates: */
    51755026        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    5176 
    5177         /*Get node data: */
    51785027        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    51795028
    5180         /*Use the dof list to index into the solution vector: */
    5181         for(i=0;i<numdof;i++){
    5182                 values[i]=solution[doflist[i]];
    5183         }
    5184 
    5185         /*Ok, we have vz in values, fill in vz array: */
    5186         for(i=0;i<NUMVERTICES;i++){
    5187                 vz[i]=values[i*numdofpervertex+0];
    5188         }
     5029        /*Use the dof list to index into the solution vector vz: */
     5030        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     5031        for(i=0;i<NUMVERTICES;i++) vz[i]=values[i*NDOF1+0];
    51895032
    51905033        /*Get Vx and Vy*/
    … …  
    52345077
    52355078        /*Now Compute vel*/
    5236 
    52375079        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    52385080
    … …  
    52715113void  Penta::InputUpdateFromSolutionDiagnosticStokes(double* solution){
    52725114       
    5273         const int    numdof=NDOF4*NUMVERTICES;
    5274         int i;
    5275         int*         doflist=NULL;
    5276         double       values[numdof];
    5277         double       vx[NUMVERTICES];
    5278         double       vy[NUMVERTICES];
    5279         double       vz[NUMVERTICES];
    5280         double       vel[NUMVERTICES];
    5281         double       pressure[NUMVERTICES];
    5282         double       stokesreconditioning;
     5115        const int numdof=NDOF4*NUMVERTICES;
     5116
     5117        int     i;
     5118        double  values[numdof];
     5119        double  vx[NUMVERTICES];
     5120        double  vy[NUMVERTICES];
     5121        double  vz[NUMVERTICES];
     5122        double  vel[NUMVERTICES];
     5123        double  pressure[NUMVERTICES];
     5124        double  stokesreconditioning;
     5125        int*    doflist=NULL;
    52835126
    52845127        /*Get dof list: */
    … …  
    52865129
    52875130        /*Use the dof list to index into the solution vector: */
    5288         for(i=0;i<numdof;i++){
    5289                 values[i]=solution[doflist[i]];
    5290         }
     5131        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    52915132
    52925133        /*Ok, we have vx and vy in values, fill in all arrays: */
    … …  
    52985139        }
    52995140
    5300         /*Recondition pressure: */
     5141        /*Recondition pressure and compute vel: */
    53015142        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    5302         for(i=0;i<NUMVERTICES;i++){
    5303                 pressure[i]=pressure[i]*stokesreconditioning;
    5304         }
    5305 
    5306         /*Now Compute vel*/
     5143        for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*stokesreconditioning;
    53075144        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    53085145       
    … …  
    53245161        xfree((void**)&doflist);
    53255162}
    5326 
    53275163/*}}}*/
    53285164/*FUNCTION Penta::InputUpdateFromSolutionAdjointStokes {{{1*/
    53295165void  Penta::InputUpdateFromSolutionAdjointStokes(double* solution){
    53305166
    5331         int i;
    5332 
    5333         const int    numdofpervertex=4;
    5334         const int    numdof=numdofpervertex*NUMVERTICES;
    5335         int*         doflist=NULL;
    5336         double       values[numdof];
    5337         double       lambdax[NUMVERTICES];
    5338         double       lambday[NUMVERTICES];
    5339         double       lambdaz[NUMVERTICES];
    5340         double       lambdap[NUMVERTICES];
     5167        const int    numdof=NDOF4*NUMVERTICES;
     5168
     5169        int    i;
     5170        double values[numdof];
     5171        double lambdax[NUMVERTICES];
     5172        double lambday[NUMVERTICES];
     5173        double lambdaz[NUMVERTICES];
     5174        double lambdap[NUMVERTICES];
     5175        int*   doflist=NULL;
    53415176
    53425177        /*Get dof list: */
    … …  
    53445179
    53455180        /*Use the dof list to index into the solution vector: */
    5346         for(i=0;i<numdof;i++){
    5347                 values[i]=solution[doflist[i]];
    5348         }
     5181        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    53495182
    53505183        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    53515184        for(i=0;i<NUMVERTICES;i++){
    5352                 lambdax[i]=values[i*numdofpervertex+0];
    5353                 lambday[i]=values[i*numdofpervertex+1];
    5354                 lambdaz[i]=values[i*numdofpervertex+2];
    5355                 lambdap[i]=values[i*numdofpervertex+3];
     5185                lambdax[i]=values[i*NDOF4+0];
     5186                lambday[i]=values[i*NDOF4+1];
     5187                lambdaz[i]=values[i*NDOF4+2];
     5188                lambdap[i]=values[i*NDOF4+3];
    53565189        }
    53575190
    … …  
    53695202void  Penta::InputUpdateFromSolutionAdjointHoriz(double* solution){
    53705203
    5371         int i;
    5372 
    5373         const int    numdofpervertex=2;
    5374         const int    numdof=numdofpervertex*NUMVERTICES;
    5375         int*         doflist=NULL;
    5376         double       values[numdof];
    5377         double       lambdax[NUMVERTICES];
    5378         double       lambday[NUMVERTICES];
     5204        const int numdof=NDOF2*NUMVERTICES;
     5205
     5206        int    i;
     5207        double values[numdof];
     5208        double lambdax[NUMVERTICES];
     5209        double lambday[NUMVERTICES];
     5210        int*   doflist=NULL;
    53795211
    53805212        /*Get dof list: */
    … …  
    53825214
    53835215        /*Use the dof list to index into the solution vector: */
    5384         for(i=0;i<numdof;i++){
    5385                 values[i]=solution[doflist[i]];
    5386         }
     5216        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    53875217
    53885218        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    53895219        for(i=0;i<NUMVERTICES;i++){
    5390                 lambdax[i]=values[i*numdofpervertex+0];
    5391                 lambday[i]=values[i*numdofpervertex+1];
     5220                lambdax[i]=values[i*NDOF2+0];
     5221                lambday[i]=values[i*NDOF2+1];
    53925222        }
    53935223
    … …  
    54035233void  Penta::InputUpdateFromSolutionThermal(double* solution){
    54045234
    5405         int i;
    5406 
    54075235        const int    numdof=NDOF1*NUMVERTICES;
    5408         int*         doflist=NULL;
    5409         double       values[numdof];
    5410         double       B[numdof];
    5411         double       B_average;
    5412         bool         converged;
     5236
     5237        bool   converged;
     5238        int    i;
     5239        double values[numdof];
     5240        double B[numdof];
     5241        double B_average;
     5242        int*   doflist=NULL;
    54135243
    54145244        /*Get dof list: */
    … …  
    54165246
    54175247        /*Use the dof list to index into the solution vector: */
    5418         for(i=0;i<numdof;i++){
    5419                 values[i]=solution[doflist[i]];
    5420         }
     5248        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    54215249
    54225250        this->inputs->GetParameterValue(&converged,ConvergedEnum);
    … …  
    54415269void  Penta::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
    54425270
    5443         const int numdofpervertex = 1;
    5444         const int numdof          = numdofpervertex *NUMVERTICES;
    5445         int*         doflist=NULL;
    5446         double    values[numdof];
     5271        const int numdof = NDOF1*NUMVERTICES;
     5272
     5273        double values[numdof];
     5274        int*   doflist=NULL;
    54475275
    54485276        /*Get dof list: */
    … …  
    54505278
    54515279        /*Use the dof list to index into the solution vector: */
    5452         for(int i=0;i<numdof;i++){
    5453                 values[i]=solution[doflist[i]];
    5454         }
     5280        for(int i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    54555281
    54565282        /*Add input to the element: */
    … …  
    54645290void  Penta::InputUpdateFromSolutionOneDofCollapsed(double* solution,int enum_type){
    54655291
    5466         const int  numdofpervertex = 1;
    5467         const int  numdof          = numdofpervertex *NUMVERTICES;
    5468         const int  numdof2d        = numdof/2;
    5469         int*         doflist=NULL;
    5470         double     values[numdof];
    5471         Penta     *penta           = NULL;
     5292        const int  numdof   = NDOF1*NUMVERTICES;
     5293        const int  numdof2d = NDOF1*NUMVERTICES2D;
     5294
     5295        double  values[numdof];
     5296        int*    doflist = NULL;
     5297        Penta  *penta   = NULL;
    54725298
    54735299        /*If not on bed, return*/
    … …  
    54865312        penta=this;
    54875313        for(;;){
    5488 
    54895314                /*Add input to the element: */
    54905315                penta->inputs->AddInput(new PentaVertexInput(enum_type,values));
    … …  
    55385363        inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum);
    55395364        return onshelf;
    5540 
    55415365}
    55425366/*}}}*/
    … …  
    55475371        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    55485372        return onwater;
    5549 
    55505373}
    55515374/*}}}*/
    … …  
    55565379        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    55575380        return onsurface;
    5558 
    55595381}
    55605382/*}}}*/
    … …  
    55655387        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    55665388        return onbed;
    5567 
    55685389}
    55695390/*}}}*/
    … …  
    55745395        int    i,node0,node1;
    55755396        int    edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges
     5397        double length;
    55765398        double minlength=-1;
    5577         double length;
    55785399       
    55795400        for(i=0;i<9;i++){
    … …  
    55935414void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
    55945415
    5595         int i,j;
    5596 
     5416        int    i,j;
    55975417        double Kii[24][24];
    55985418        double Kib[24][3];
    … …  
    56305450
    56315451        /*Affect value to the reduced matrix */
    5632         for(i=0;i<24;i++){
    5633                 for(j=0;j<24;j++){
    5634                         *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];
    5635                 }
    5636         }
     5452        for(i=0;i<24;i++) for(j=0;j<24;j++) *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];
    56375453}
    56385454/*}}}*/
    … …  
    56405456void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
    56415457
    5642         int i,j;
    5643 
     5458        int    i,j;
    56445459        double Pi[24];
    56455460        double Pb[3];
    … …  
    56505465
    56515466        /*Create the four matrices used for reduction */
    5652         for(i=0;i<24;i++){
    5653                 Pi[i]=*(Pe_temp+i);
    5654         }
    5655         for(i=0;i<3;i++){
    5656                 Pb[i]=*(Pe_temp+i+24);
    5657         }
     5467        for(i=0;i<24;i++) Pi[i]=*(Pe_temp+i);
     5468        for(i=0;i<3;i++) Pb[i]=*(Pe_temp+i+24);
    56585469        for(j=0;j<3;j++){
    56595470                for(i=0;i<24;i++){
    … …  
    56745485
    56755486        /*Affect value to the reduced matrix */
    5676         for(i=0;i<24;i++){
    5677                 *(Pe_reduced+i)=Pi[i]-Pright[i];
    5678         }
     5487        for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i];
    56795488}
    56805489/*}}}*/
    … …  
    56885497Tria*  Penta::SpawnTria(int g0, int g1, int g2){
    56895498
    5690         int i;
    5691         int analysis_counter;
     5499        int   i,analysis_counter;
     5500        int   indices[3];
     5501        int   zero=0;
     5502        Tria*       tria            = NULL;
     5503        Inputs*     tria_inputs     = NULL;
     5504        Results*    tria_results    = NULL;
     5505        Parameters* tria_parameters = NULL;
    56925506       
    56935507        /*go into parameters and get the analysis_counter: */
    56945508        this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
    5695 
    5696         /*out of grids g0,g1 and g2 from Penta, build a tria element: */
    5697         Tria* tria=NULL;
    5698         int indices[3];
    5699         int zero=0;
    5700         Parameters* tria_parameters=NULL;
    5701         Inputs* tria_inputs=NULL;
    5702         Results* tria_results=NULL;
    57035509
    57045510        indices[0]=g0;
    … …  
    57345540void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
    57355541
    5736         int i;
    5737         double v13[3];
    5738         double v23[3];
     5542        int    i;
     5543        double v13[3],v23[3];
    57395544        double normal[3];
    57405545        double normal_norm;
Note: See TracChangeset for help on using the changeset viewer.