Changeset 6158


Ignore:
Timestamp:
10/05/10 15:13:08 (14 years ago)
Author:
seroussi
Message:

some cleaning in Tria

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6131 r6158  
    540540void  Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){
    541541
     542        bool      already=false;
    542543        int       i,j;
     544        int       partition[NUMVERTICES];
     545        int       offset[NUMVERTICES];
    543546        double    area;
    544547        double    mean;
    545         int       partition[NUMVERTICES];
    546         int       offset[NUMVERTICES];
    547548        double    values[3];
    548         bool      already=false;
    549549
    550550        /*First, get the area: */
    … …  
    900900        /*Intermediaries*/
    901901        int        i,ig;
    902         double     Jdet;
    903         double     viscosity_complement;
    904         double     vx,vy,lambda,mu,thickness;
    905         double     cm_noisedmp;
    906902        int        doflist[NUMVERTICES];
     903        double     vx,vy,lambda,mu,thickness,Jdet;
     904        double     cm_noisedmp,viscosity_complement;
    907905        double     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2];
    908906        double     xyz_list[NUMVERTICES][3];
    909         double     basis[3];
     907        double     basis[3],epsilon[3];
    910908        double     dbasis[NDOF2][NUMVERTICES];
     909        double     grad_g[NUMVERTICES];
    911910        double     grad[NUMVERTICES]={0.0};
    912         double     grad_g[NUMVERTICES];
    913         double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    914911        GaussTria *gauss = NULL;
    915912
    … …  
    954951                }
    955952                /*Add regularization term*/
    956                 for (i=0;i<NUMVERTICES;i++){
    957                         grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]);
    958                 }
    959 
     953                for (i=0;i<NUMVERTICES;i++) grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]);
    960954                for(i=0;i<NUMVERTICES;i++) grad[i]+=grad_g[i];
    961955        }
    962956
    963         /*Add grade_g to global vector gradient: */
    964957        VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES);
    965958
    … …  
    971964void  Tria::GradjDrag(Vec gradient){
    972965
    973         int i;
    974         int     ig;
    975 
    976         double       xyz_list[NUMVERTICES][3];
    977         int          doflist1[NUMVERTICES];
    978         double       dh1dh3[NDOF2][NUMVERTICES];
    979         GaussTria *gauss=NULL;
    980         double  dk[NDOF2];
    981         double  vx,vy;
    982         double  lambda,mu;
    983         double  bed,thickness,Neff;
    984         double  alpha_complement;
    985         int     drag_type;
    986         double  drag;
    987         Friction* friction=NULL;
    988         double  grade_g[NUMVERTICES]={0.0};
    989         double  grade_g_gaussian[NUMVERTICES];
    990         double Jdet;
    991         double l1l2l3[3];
    992         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    993         double  cm_noisedmp;
    994         int analysis_type;
     966        int        i,ig;
     967        int        drag_type,analysis_type;
     968        int        doflist1[NUMVERTICES];
     969        double     vx,vy,lambda,mu,alpha_complement,Jdet;
     970        double     bed,thickness,Neff,drag,cm_noisedmp;
     971        double     xyz_list[NUMVERTICES][3];
     972        double     dh1dh3[NDOF2][NUMVERTICES];
     973        double     dk[NDOF2];
     974        double     grade_g[NUMVERTICES]={0.0};
     975        double     grade_g_gaussian[NUMVERTICES];
     976        double     l1l2l3[3];
     977        double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     978        Friction*  friction=NULL;
     979        GaussTria  *gauss=NULL;
    995980
    996981        /*retrive parameters: */
    997982        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    998983
    999         /*retrieve some parameters: */
     984        /*retrieve some parameters ands return if iceshelf: */
    1000985        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    1001 
    1002         /*Get out if shelf*/
    1003986        if(IsOnShelf())return;
    1004 
    1005         /* Get node coordinates and dof list: */
    1006987        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1007988        GetDofList1(&doflist1[0]);
    … …  
    10241005                gauss->GaussPoint(ig);
    10251006
     1007                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1008                GetNodalFunctions(l1l2l3, gauss);
     1009                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     1010
    10261011                /*Build alpha_complement_list: */
    10271012                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
    10281013                else alpha_complement=0;
    10291014       
    1030                 /*Recover alpha_complement and k: */
    10311015                dragcoefficient_input->GetParameterValue(&drag, gauss);
    1032 
    1033                 /*recover lambda and mu: */
    10341016                adjointx_input->GetParameterValue(&lambda, gauss);
    10351017                adjointy_input->GetParameterValue(&mu, gauss);
    1036                        
    1037                 /*recover vx and vy: */
    10381018                vx_input->GetParameterValue(&vx,gauss);
    10391019                vy_input->GetParameterValue(&vy,gauss);
    1040 
    1041                 /* Get Jacobian determinant: */
    1042                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    1043                
    1044                 /* Get nodal functions value at gaussian point:*/
    1045                 GetNodalFunctions(l1l2l3, gauss);
    1046 
    1047                 /*Get nodal functions derivatives*/
    1048                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
    1049 
    1050                 /*Get k derivative: dk/dx */
    10511020                dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    10521021
    … …  
    10651034        }
    10661035
    1067         /*Add grade_g to global vector gradient: */
    10681036        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    10691037
    … …  
    10711039        delete gauss;
    10721040        delete friction;
    1073 
    10741041}
    10751042/*}}}*/
    … …  
    10821049        double gradient_g[NUMVERTICES];
    10831050
    1084         /*Retrieve dof list*/
    10851051        GetDofList1(&doflist1[0]);
    10861052
    … …  
    10891055        for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];
    10901056
    1091         /*Add grade_g to global vector gradient: */
    10921057        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)gradient_g,INSERT_VALUES);
    10931058}
    … …  
    10951060/*FUNCTION Tria::GradjVx{{{1*/
    10961061void  Tria::GradjVx(Vec gradient){
    1097 
    1098         /*Intermediaries*/
    1099         int        i,ig;
    1100         int        doflist1[NUMVERTICES];
    1101         double     l1l2l3[3];
    1102         double     thickness,Jdet;
    1103         double     Dlambda[2];
    1104         double     xyz_list[NUMVERTICES][3];
    1105         GaussTria *gauss                = NULL;
    1106         double     grade_g[NUMVERTICES] = {0.0};
    1107 
    1108         /* Get node coordinates and dof list: */
    1109         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1110         GetDofList1(&doflist1[0]);
    1111 
    1112         /*Retrieve all inputs we will be needing: */
    1113         Input* adjoint_input=inputs->GetInput(AdjointEnum);     ISSMASSERT(adjoint_input);
    1114         Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    1115 
    1116         /* Start  looping on the number of gaussian points: */
    1117         gauss=new GaussTria(2);
    1118         for (ig=gauss->begin();ig<gauss->end();ig++){
    1119 
    1120                 gauss->GaussPoint(ig);
    1121 
    1122                 adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    1123                 thickness_input->GetParameterValue(&thickness, gauss);
    1124 
    1125                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    1126                 GetNodalFunctions(l1l2l3, gauss);
    1127 
    1128                 for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
    1129         }
    1130 
    1131         VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    1132 
    1133         /*Clean up and return*/
    1134         delete gauss;
    1135 }
    1136 /*}}}*/
    1137 /*FUNCTION Tria::GradjVy{{{1*/
    1138 void  Tria::GradjVy(Vec gradient){
    11391062
    11401063        /*Intermediaries*/
    … …  
    11451068        double     Dlambda[2];
    11461069        double     xyz_list[NUMVERTICES][3];
     1070        double     grade_g[NUMVERTICES] = {0.0};
    11471071        GaussTria *gauss                = NULL;
    1148         double     grade_g[NUMVERTICES] = {0.0};
    11491072
    11501073        /* Get node coordinates and dof list: */
    … …  
    11621085                gauss->GaussPoint(ig);
    11631086
     1087                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1088                GetNodalFunctions(l1l2l3, gauss);
     1089               
    11641090                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    11651091                thickness_input->GetParameterValue(&thickness, gauss);
    11661092
     1093                for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
     1094        }
     1095
     1096        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     1097
     1098        /*Clean up and return*/
     1099        delete gauss;
     1100}
     1101/*}}}*/
     1102/*FUNCTION Tria::GradjVy{{{1*/
     1103void  Tria::GradjVy(Vec gradient){
     1104
     1105        /*Intermediaries*/
     1106        int        i,ig;
     1107        int        doflist1[NUMVERTICES];
     1108        double     thickness,Jdet;
     1109        double     l1l2l3[3];
     1110        double     Dlambda[2];
     1111        double     xyz_list[NUMVERTICES][3];
     1112        double     grade_g[NUMVERTICES] = {0.0};
     1113        GaussTria *gauss                = NULL;
     1114
     1115        /* Get node coordinates and dof list: */
     1116        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1117        GetDofList1(&doflist1[0]);
     1118
     1119        /*Retrieve all inputs we will be needing: */
     1120        Input* adjoint_input=inputs->GetInput(AdjointEnum);     ISSMASSERT(adjoint_input);
     1121        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     1122
     1123        /* Start  looping on the number of gaussian points: */
     1124        gauss=new GaussTria(2);
     1125        for (ig=gauss->begin();ig<gauss->end();ig++){
     1126
     1127                gauss->GaussPoint(ig);
     1128
     1129                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
     1130                thickness_input->GetParameterValue(&thickness, gauss);
     1131
    11671132                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    11681133                GetNodalFunctions(l1l2l3, gauss);
    … …  
    11811146
    11821147        /*Intermediary*/
     1148        int    control_type;
    11831149        Input* input=NULL;
    11841150        double cm_min,cm_max;
    1185         int    control_type;
    11861151
    11871152        /*retrieve some parameters: */
    … …  
    12491214                        ISSMERROR("control type %s not implemented yet",EnumToString(control_type));
    12501215        }
    1251 
    12521216}
    12531217/*}}}*/
    … …  
    12551219bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
    12561220
    1257         int i;
     1221        bool    converged=true;
     1222        int     i;
    12581223        Input** new_inputs=NULL;
    12591224        Input** old_inputs=NULL;
    1260         bool    converged=true;
    12611225
    12621226        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
    … …  
    12801244        xfree((void**)&old_inputs);
    12811245        return converged;
    1282 
    12831246}
    12841247/*}}}*/
    … …  
    13101273        else
    13111274         ISSMERROR("object %s not supported yet",EnumToString(object_enum));
    1312 
    13131275}
    13141276/*}}}*/
    … …  
    13541316
    13551317        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1356         if (enum_type==RheologyBbarEnum){
    1357                 input=this->matice->inputs->GetInput(enum_type);
    1358         }
    1359         else{
    1360                 input=this->inputs->GetInput(enum_type);
    1361         }
     1318        if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
     1319        else input=this->inputs->GetInput(enum_type);
     1320
    13621321        if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type));
    13631322
    … …  
    13651324         * object out of the input, with the additional step and time information: */
    13661325        this->results->AddObject((Object*)input->SpawnResult(step,time));
    1367 
    13681326}
    13691327/*}}}*/
    … …  
    13711329double Tria::MassFlux( double* segment,bool process_units){
    13721330
    1373         int i;
    1374 
    13751331        const int    numdofs=2;
    1376         double mass_flux=0;
    1377         double xyz_list[NUMVERTICES][3];
     1332
     1333        int        i;
     1334        double     mass_flux=0;
     1335        double     xyz_list[NUMVERTICES][3];
     1336        double     normal[2];
     1337        double     length,rho_ice;
     1338        double     x1,y1,x2,y2,h1,h2;
     1339        double     vx1,vx2,vy1,vy2;
    13781340        GaussTria* gauss_1=NULL;
    13791341        GaussTria* gauss_2=NULL;
    1380         double normal[2];
    1381         double length;
    1382         double x1,y1,x2,y2;
    1383         double h1,h2;
    1384         double vx1,vx2,vy1,vy2;
    1385         double rho_ice;
    13861342
    13871343        /*Get material parameters :*/
    … …  
    14031359        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    14041360
    1405         /*get normal of segment: */
    14061361        normal[0]=cos(atan2(x1-x2,y2-y1));
    14071362        normal[1]=sin(atan2(x1-x2,y2-y1));
    14081363
    1409         /*get length of segment: */
    14101364        length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
    14111365
    1412         /*get thickness and velocity at two segment extremities: */
    14131366        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     1367        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
     1368        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
     1369
    14141370        thickness_input->GetParameterValue(&h1, gauss_1);
    14151371        thickness_input->GetParameterValue(&h2, gauss_2);
    1416 
    1417         Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
    14181372        vx_input->GetParameterValue(&vx1,gauss_1);
    14191373        vx_input->GetParameterValue(&vx2,gauss_2);
    1420 
    1421         Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    14221374        vy_input->GetParameterValue(&vy1,gauss_1);
    14231375        vy_input->GetParameterValue(&vy2,gauss_2);
    … …  
    14311383        mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters);
    14321384
    1433         /*clean up*/
     1385        /*clean up and return:*/
    14341386        delete gauss_1;
    14351387        delete gauss_2;
    1436 
    14371388        return mass_flux;
    14381389}
    … …  
    14881439        /*Assign output pointers:*/
    14891440        *pmaxvel=maxvel;
    1490 
    14911441}
    14921442/*}}}*/
    … …  
    15021452        /*Assign output pointers:*/
    15031453        *pmaxvx=maxvx;
    1504 
    15051454}
    15061455/*}}}*/
    … …  
    15301479        /*Assign output pointers:*/
    15311480        *pmaxvz=maxvz;
    1532 
    15331481}
    15341482/*}}}*/
    … …  
    15441492        /*Assign output pointers:*/
    15451493        *pminvel=minvel;
    1546 
    15471494}
    15481495/*}}}*/
    … …  
    15581505        /*Assign output pointers:*/
    15591506        *pminvx=minvx;
    1560 
    15611507}
    15621508/*}}}*/
    … …  
    15721518        /*Assign output pointers:*/
    15731519        *pminvy=minvy;
    1574 
    15751520}
    15761521/*}}}*/
    … …  
    15861531        /*Assign output pointers:*/
    15871532        *pminvz=minvz;
    1588 
    15891533}
    15901534/*}}}*/
    … …  
    15931537
    15941538        /*intermediary: */
    1595         double C;
    1596         double maxabsvx;
    1597         double maxabsvy;
     1539        int    i;
     1540        double C,dt;
     1541        double dx,dy;
    15981542        double maxx,minx;
    15991543        double maxy,miny;
     1544        double maxabsvx,maxabsvy;
    16001545        double xyz_list[NUMVERTICES][3];
    1601         double dx,dy;
    1602         int    i;
    1603 
    1604         /*output: */
    1605         double dt;
    16061546
    16071547        /*get CFL coefficient:*/
    … …  
    16331573
    16341574        return dt;
    1635 
    16361575}
    16371576/*}}}*/
    … …  
    16871626double Tria::SurfaceAbsVelMisfit(bool process_units){
    16881627
    1689         int i;
    1690 
    1691         /* output: */
    1692         double Jelem=0;
    1693 
    1694         /* node data: */
    1695         const int    numdof=2*NUMVERTICES;
    1696         double       xyz_list[NUMVERTICES][3];
    1697 
    1698         /* grid data: */
    1699         double vx_list[NUMVERTICES];
    1700         double vy_list[NUMVERTICES];
    1701         double obs_vx_list[NUMVERTICES];
    1702         double obs_vy_list[NUMVERTICES];
    1703         double misfit_square_list[NUMVERTICES];
    1704         double misfit_list[NUMVERTICES];
    1705         double weights_list[NUMVERTICES];
    1706 
    1707         /* gaussian points: */
    1708         int     ig;
     1628        const int    numdof=NDOF2*NUMVERTICES;
     1629
     1630        int        i,ig;
     1631        double     Jelem=0;
     1632        double     velocity_mag,obs_velocity_mag;
     1633        double     meanvel, epsvel,misfit,Jdet;
     1634        double     xyz_list[NUMVERTICES][3];
     1635        double     vx_list[NUMVERTICES];
     1636        double     vy_list[NUMVERTICES];
     1637        double     obs_vx_list[NUMVERTICES];
     1638        double     obs_vy_list[NUMVERTICES];
     1639        double     misfit_square_list[NUMVERTICES];
     1640        double     misfit_list[NUMVERTICES];
     1641        double     weights_list[NUMVERTICES];
    17091642        GaussTria *gauss=NULL;
    1710 
    1711         /* parameters: */
    1712         double  velocity_mag,obs_velocity_mag;
    1713         double  misfit;
    1714 
    1715         /* Jacobian: */
    1716         double Jdet;
    1717         double  meanvel, epsvel;
    17181643
    17191644        /*If on water, return 0: */
    … …  
    17581683
    17591684        /*Apply weights to misfits*/
    1760         for (i=0;i<NUMVERTICES;i++){
    1761                 misfit_list[i]=weights_list[i]*misfit_list[i];
    1762         }
     1685        for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    17631686
    17641687        /* Start  looping on the number of gaussian points: */
    … …  
    17811704double Tria::SurfaceRelVelMisfit(bool process_units){
    17821705
    1783         int i;
    1784 
    1785         /* output: */
    1786         double Jelem=0;
    1787 
    1788         /* node data: */
    17891706        const int    numdof=2*NUMVERTICES;
    1790         double       xyz_list[NUMVERTICES][3];
    1791 
    1792         /* grid data: */
    1793         double vx_list[NUMVERTICES];
    1794         double vy_list[NUMVERTICES];
    1795         double obs_vx_list[NUMVERTICES];
    1796         double obs_vy_list[NUMVERTICES];
    1797         double misfit_square_list[NUMVERTICES];
    1798         double misfit_list[NUMVERTICES];
    1799         double weights_list[NUMVERTICES];
    1800 
    1801         /* gaussian points: */
    1802         int     ig;
     1707
     1708        int        i,ig;
     1709        double     Jelem=0;
     1710        double     scalex=1,scaley=1;
     1711        double     meanvel, epsvel,misfit,Jdet;
     1712        double     velocity_mag,obs_velocity_mag;
     1713        double     xyz_list[NUMVERTICES][3];
     1714        double     vx_list[NUMVERTICES];
     1715        double     vy_list[NUMVERTICES];
     1716        double     obs_vx_list[NUMVERTICES];
     1717        double     obs_vy_list[NUMVERTICES];
     1718        double     misfit_square_list[NUMVERTICES];
     1719        double     misfit_list[NUMVERTICES];
     1720        double     weights_list[NUMVERTICES];
    18031721        GaussTria *gauss=NULL;
    1804 
    1805         /* parameters: */
    1806         double  velocity_mag,obs_velocity_mag;
    1807         double  misfit;
    1808 
    1809         /* Jacobian: */
    1810         double Jdet;
    1811 
    1812         /*relative and logarithmic control method :*/
    1813         double scalex=1;
    1814         double scaley=1;
    1815         double  meanvel, epsvel;
    18161722
    18171723        /*If on water, return 0: */
    … …  
    18611767
    18621768        /*Apply weights to misfits*/
    1863         for (i=0;i<NUMVERTICES;i++){
    1864                 misfit_list[i]=weights_list[i]*misfit_list[i];
    1865         }
     1769        for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    18661770
    18671771        /* Start  looping on the number of gaussian points: */
    … …  
    18841788double Tria::SurfaceLogVelMisfit(bool process_units){
    18851789
    1886         int i;
    1887 
    1888         /* output: */
    1889         double Jelem=0;
    1890 
    1891         /* node data: */
    1892         const int    numdof=2*NUMVERTICES;
    1893         double       xyz_list[NUMVERTICES][3];
    1894 
    1895         /* grid data: */
    1896         double vx_list[NUMVERTICES];
    1897         double vy_list[NUMVERTICES];
    1898         double obs_vx_list[NUMVERTICES];
    1899         double obs_vy_list[NUMVERTICES];
    1900         double misfit_square_list[NUMVERTICES];
    1901         double misfit_list[NUMVERTICES];
    1902         double weights_list[NUMVERTICES];
    1903 
    1904         /* gaussian points: */
    1905         int     ig;
     1790        const int    numdof=NDOF2*NUMVERTICES;
     1791
     1792        int        i,ig;
     1793        double     Jelem=0;
     1794        double     scalex=1,scaley=1;
     1795        double     meanvel, epsvel,misfit,Jdet;
     1796        double     velocity_mag,obs_velocity_mag;
     1797        double     xyz_list[NUMVERTICES][3];
     1798        double     vx_list[NUMVERTICES];
     1799        double     vy_list[NUMVERTICES];
     1800        double     obs_vx_list[NUMVERTICES];
     1801        double     obs_vy_list[NUMVERTICES];
     1802        double     misfit_square_list[NUMVERTICES];
     1803        double     misfit_list[NUMVERTICES];
     1804        double     weights_list[NUMVERTICES];
    19061805        GaussTria *gauss=NULL;
    1907 
    1908         /* parameters: */
    1909         double  velocity_mag,obs_velocity_mag;
    1910         double  misfit;
    1911 
    1912         /* Jacobian: */
    1913         double Jdet;
    1914 
    1915         /*relative and logarithmic control method :*/
    1916         double scalex=1;
    1917         double scaley=1;
    1918         double  meanvel, epsvel;
    19191806
    19201807        /*If on water, return 0: */
    … …  
    19621849
    19631850        /*Apply weights to misfits*/
    1964         for (i=0;i<NUMVERTICES;i++){
    1965                 misfit_list[i]=weights_list[i]*misfit_list[i];
    1966         }
     1851        for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    19671852
    19681853        /* Start  looping on the number of gaussian points: */
    … …  
    19851870double Tria::SurfaceLogVxVyMisfit(bool process_units){
    19861871
    1987         int i;
    1988 
    1989         /* output: */
    1990         double Jelem=0;
    1991 
    1992         /* node data: */
    1993         const int    numdof=2*NUMVERTICES;
    1994         double       xyz_list[NUMVERTICES][3];
    1995 
    1996         /* grid data: */
    1997         double vx_list[NUMVERTICES];
    1998         double vy_list[NUMVERTICES];
    1999         double obs_vx_list[NUMVERTICES];
    2000         double obs_vy_list[NUMVERTICES];
    2001         double misfit_square_list[NUMVERTICES];
    2002         double misfit_list[NUMVERTICES];
    2003         double weights_list[NUMVERTICES];
    2004 
    2005         /* gaussian points: */
    2006         int     ig;
     1872        const int    numdof=NDOF2*NUMVERTICES;
     1873
     1874        int        i,ig;
     1875        int        fit=-1;
     1876        double     Jelem=0, S=0;
     1877        double     scalex=1,scaley=1;
     1878        double     meanvel, epsvel, misfit, Jdet;
     1879        double     velocity_mag,obs_velocity_mag;
     1880        double     xyz_list[NUMVERTICES][3];
     1881        double     vx_list[NUMVERTICES];
     1882        double     vy_list[NUMVERTICES];
     1883        double     obs_vx_list[NUMVERTICES];
     1884        double     obs_vy_list[NUMVERTICES];
     1885        double     misfit_square_list[NUMVERTICES];
     1886        double     misfit_list[NUMVERTICES];
     1887        double     weights_list[NUMVERTICES];
    20071888        GaussTria *gauss=NULL;
    2008 
    2009         /* parameters: */
    2010         double  velocity_mag,obs_velocity_mag;
    2011         double  misfit;
    2012 
    2013         /* Jacobian: */
    2014         double Jdet;
    2015 
    2016         /*relative and logarithmic control method :*/
    2017         double scalex=1;
    2018         double scaley=1;
    2019         double S=0;
    2020         int    fit=-1;
    2021         double  meanvel, epsvel;
    20221889
    20231890        /*If on water, return 0: */
    … …  
    20651932
    20661933        /*Apply weights to misfits*/
    2067         for (i=0;i<NUMVERTICES;i++){
    2068                 misfit_list[i]=weights_list[i]*misfit_list[i];
    2069         }
     1934        for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    20701935
    20711936        /* Start  looping on the number of gaussian points: */
    … …  
    20881953double Tria::SurfaceAverageVelMisfit(bool process_units){
    20891954
    2090         int i;
    2091 
    2092         /* output: */
    2093         double Jelem=0;
    2094 
    2095         /* node data: */
    20961955        const int    numdof=2*NUMVERTICES;
    2097         double       xyz_list[NUMVERTICES][3];
    2098 
    2099         /* grid data: */
    2100         double vx_list[NUMVERTICES];
    2101         double vy_list[NUMVERTICES];
    2102         double obs_vx_list[NUMVERTICES];
    2103         double obs_vy_list[NUMVERTICES];
    2104         double misfit_square_list[NUMVERTICES];
    2105         double misfit_list[NUMVERTICES];
    2106         double weights_list[NUMVERTICES];
    2107 
    2108         /* gaussian points: */
    2109         int     ig;
     1956
     1957        int        i,ig;
     1958        int        fit=-1;
     1959        double     Jelem=0,S=0;
     1960        double     scalex=1, scaley=1;
     1961        double     meanvel, epsvel,Jdet;
     1962        double     velocity_mag,obs_velocity_mag,misfit;
     1963        double     xyz_list[NUMVERTICES][3];
     1964        double     vx_list[NUMVERTICES];
     1965        double     vy_list[NUMVERTICES];
     1966        double     obs_vx_list[NUMVERTICES];
     1967        double     obs_vy_list[NUMVERTICES];
     1968        double     misfit_square_list[NUMVERTICES];
     1969        double     misfit_list[NUMVERTICES];
     1970        double     weights_list[NUMVERTICES];
    21101971        GaussTria *gauss=NULL;
    2111 
    2112         /* parameters: */
    2113         double  velocity_mag,obs_velocity_mag;
    2114         double  misfit;
    2115 
    2116         /* Jacobian: */
    2117         double Jdet;
    2118 
    2119         /*relative and logarithmic control method :*/
    2120         double scalex=1;
    2121         double scaley=1;
    2122         double S=0;
    2123         int    fit=-1;
    2124         double  meanvel, epsvel;
    21251972
    21261973        /*If on water, return 0: */
    … …  
    21672014
    21682015        /*Apply weights to misfits*/
    2169         for (i=0;i<NUMVERTICES;i++){
    2170                 misfit_list[i]=weights_list[i]*misfit_list[i];
    2171         }
     2016        for (i=0;i<NUMVERTICES;i++) misfit_list[i]=weights_list[i]*misfit_list[i];
    21722017
    21732018        /* Start  looping on the number of gaussian points: */
    … …  
    21832028
    21842029        /*clean-up and Return: */
    2185                 delete gauss;
     2030        delete gauss;
    21862031        return Jelem;
    21872032}
    … …  
    21902035void  Tria::PatchFill(int* prow, Patch* patch){
    21912036
    2192         int i;
    2193         int row;
     2037        int i,row;
    21942038        int vertices_ids[3];
    21952039
    … …  
    21972041        row=*prow;
    21982042               
    2199         /*will be needed later: */
    22002043        for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
    22012044
    … …  
    22202063
    22212064        int     i;
    2222        
    2223         /*output: */
    22242065        int     numrows     = 0;
    22252066        int     numnodes    = 0;
    … …  
    22382079        *pnumvertices=NUMVERTICES;
    22392080        *pnumnodes=numnodes;
    2240        
    22412081}
    22422082/*}}}*/
    … …  
    22502090                elementresult->ProcessUnits(this->parameters);
    22512091        }
    2252 
    22532092}
    22542093/*}}}*/
    … …  
    22562095double Tria::SurfaceArea(void){
    22572096
    2258         int i;
    2259 
    2260         /* output: */
     2097        int    i;
    22612098        double S;
    2262 
    2263         /* node data: */
     2099        double normal[3];
     2100        double v13[3],v23[3];
    22642101        double xyz_list[NUMVERTICES][3];
    2265         double v13[3];
    2266         double v23[3];
    2267         double normal[3];
    22682102
    22692103        /*If on water, return 0: */
    22702104        if(IsOnWater())return 0;
    22712105
    2272         /* Get node coordinates and dof list: */
    22732106        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    22742107
    … …  
    22952128        int    tria_node_ids[3];
    22962129        int    tria_vertex_ids[3];
     2130        int    tria_type;
    22972131        double nodeinputs[3];
    2298         int    tria_type;
    22992132
    23002133        /*Checks if debuging*/
    … …  
    25042337        /*If sheet: surface = bed + thickness*/
    25052338        else{
    2506 
     2339               
    25072340                /*The bed does not change, update surface only s = b + h*/
    25082341                this->inputs->DuplicateInput(BedEnum,SurfaceEnum);          //1: copy bed into surface
    25092342                this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum);     //2: surface = surface + 1 * thickness
    25102343        }
    2511 
    25122344}
    25132345/*}}}*/
    … …  
    46324464
    46334465        int i,j;
     4466        int count=0;
    46344467        int numberofdofs=0;
    4635         int count=0;
    4636 
    4637         /*output: */
    46384468        int* doflist=NULL;
    46394469
    4640         /*First, figure out size of doflist: */
    4641         for(i=0;i<3;i++){
    4642                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    4643         }
    4644 
    4645         /*Allocate: */
     4470        /*First, figure out size of doflist and create it: */
     4471        for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    46464472        doflist=(int*)xmalloc(numberofdofs*sizeof(int));
    46474473
    … …  
    46554481        /*Assign output pointers:*/
    46564482        *pdoflist=doflist;
    4657 
    46584483}
    46594484/*}}}*/
    … …  
    46624487
    46634488        int i;
    4664         for(i=0;i<3;i++){
    4665                 doflist[i]=nodes[i]->GetDofList1();
    4666         }
     4489        for(i=0;i<3;i++) doflist[i]=nodes[i]->GetDofList1();
    46674490
    46684491}
    … …  
    46964519void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
    46974520
    4698         /*Intermediaries*/
    46994521        double     value[NUMVERTICES];
    4700         GaussTria *gauss              = NULL;
    4701 
    4702         /*Recover input*/
    4703         Input* input=inputs->GetInput(enumtype);
     4522        GaussTria *gauss = NULL;
     4523        Input     *input = inputs->GetInput(enumtype);
    47044524
    47054525        /*Checks in debugging mode*/
    … …  
    47384558void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
    47394559
    4740         int i;
    4741         const int    numdofpervertex=2;
    4742         const int    numdof=numdofpervertex*NUMVERTICES;
     4560        const int    numdof=NDOF2*NUMVERTICES;
     4561
     4562        int          i;
    47434563        int*         doflist=NULL;
     4564        double       vx,vy;
    47444565        double       values[numdof];
    4745         double       vx;
    4746         double       vy;
    47474566        GaussTria*   gauss=NULL;
    47484567
    … …  
    47644583                vx_input->GetParameterValue(&vx,gauss);
    47654584                vy_input->GetParameterValue(&vy,gauss);
    4766                 values[i*numdofpervertex+0]=vx;
    4767                 values[i*numdofpervertex+1]=vy;
    4768         }
    4769 
    4770         /*Add value to global vector*/
     4585                values[i*NDOF2+0]=vx;
     4586                values[i*NDOF2+1]=vy;
     4587        }
     4588
    47714589        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    47724590
    … …  
    47744592        delete gauss;
    47754593        xfree((void**)&doflist);
    4776 
    47774594}
    47784595/*}}}*/
    … …  
    47804597void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
    47814598
    4782         int i;
    4783         const int    numdofpervertex=2;
    4784         const int    numdof=numdofpervertex*NUMVERTICES;
     4599        const int    numdof=NDOF2*NUMVERTICES;
     4600
     4601        int i,dummy;
    47854602        int*         doflist=NULL;
     4603        double       vx,vy;
    47864604        double       values[numdof];
    4787         double       vx;
    4788         double       vy;
    4789         int          dummy;
    47904605        GaussTria*   gauss=NULL;
    47914606
    … …  
    48074622                vx_input->GetParameterValue(&vx,gauss);
    48084623                vy_input->GetParameterValue(&vy,gauss);
    4809                 values[i*numdofpervertex+0]=vx;
    4810                 values[i*numdofpervertex+1]=vy;
    4811         }
    4812 
    4813         /*Add value to global vector*/
     4624                values[i*NDOF2+0]=vx;
     4625                values[i*NDOF2+1]=vy;
     4626        }
     4627
    48144628        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    48154629
    … …  
    48174631        delete gauss;
    48184632        xfree((void**)&doflist);
    4819 
    48204633}
    48214634/*}}}*/
    … …  
    48234636void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
    48244637        /*Compute the 2d Strain Rate (3 components):
    4825          *
    4826          * epsilon=[exx eyy exy]
    4827          */
     4638         * epsilon=[exx eyy exy] */
    48284639
    48294640        int i;
    4830 
    48314641        double epsilonvx[3];
    48324642        double epsilonvy[3];
    … …  
    48434653        /*Sum all contributions*/
    48444654        for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    4845 
    48464655}
    48474656/*}}}*/
    … …  
    48494658void  Tria::GradjDragStokes(Vec gradient){
    48504659
    4851         int i,ig;
    4852 
    4853         /* node data: */
    4854         double       xyz_list[NUMVERTICES][3];
    4855         int          doflist1[NUMVERTICES];
    4856         double       dh1dh3[NDOF2][NUMVERTICES];
    4857 
    4858         /* grid data: */
    4859         double drag;
    4860         double alpha_complement;
    4861         Friction* friction=NULL;
     4660        int        i,ig;
     4661        int        drag_type,analysis_type;
     4662        int        doflist1[NUMVERTICES];
     4663        double     bed,thickness,Neff;
     4664        double     lambda,mu,xi,Jdet,vx,vy,vz;
     4665        double     alpha_complement,drag,cm_noisedmp;
     4666        double     surface_normal[3],bed_normal[3];
     4667        double     xyz_list[NUMVERTICES][3];
     4668        double     dh1dh3[NDOF2][NUMVERTICES];
     4669        double     dk[NDOF2];
     4670        double     l1l2l3[3];
     4671        double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     4672        double     grade_g[NUMVERTICES]={0.0};
     4673        double     grade_g_gaussian[NUMVERTICES];
     4674        Friction*  friction=NULL;
    48624675        GaussTria* gauss=NULL;
    4863 
    4864         /* parameters: */
    4865         double  vx,vy,vz;
    4866         double  lambda,mu,xi;
    4867         double  bed,thickness,Neff;
    4868         double  surface_normal[3];
    4869         double  bed_normal[3];
    4870         double  dk[NDOF2];
    4871 
    4872         /*element vector at the gaussian points: */
    4873         double  grade_g[NUMVERTICES]={0.0};
    4874         double  grade_g_gaussian[NUMVERTICES];
    4875         double Jdet;
    4876         double l1l2l3[3];
    4877         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    4878         int  drag_type;
    4879         double  cm_noisedmp;
    4880         int analysis_type;
    48814676
    48824677        /*retrive parameters: */
    … …  
    49644759        }
    49654760
    4966         /*Add grade_g to global vector gradient: */
    49674761        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    49684762
    49694763        delete friction;
    49704764        delete gauss;
    4971 
    49724765}
    49734766/*}}}*/
    … …  
    49754768void  Tria::InputUpdateFromSolutionAdjointBalancedthickness(double* solution){
    49764769
    4977         int i;
    4978 
    4979         const int    numdofpervertex=1;
    4980         const int    numdof=numdofpervertex*NUMVERTICES;
    4981         int*         doflist=NULL;
    4982         double       values[numdof];
    4983         double       lambda[NUMVERTICES];
     4770        const int numdof=NDOF1*NUMVERTICES;
     4771
     4772        int       i;
     4773        int*      doflist=NULL;
     4774        double    values[numdof];
     4775        double    lambda[NUMVERTICES];
    49844776
    49854777        /*Get dof list: */
    … …  
    49874779
    49884780        /*Use the dof list to index into the solution vector: */
    4989         for(i=0;i<numdof;i++){
    4990                 values[i]=solution[doflist[i]];
    4991         }
     4781        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    49924782
    49934783        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4994         for(i=0;i<numdof;i++){
    4995                 lambda[i]=values[i];
    4996         }
     4784        for(i=0;i<numdof;i++) lambda[i]=values[i];
    49974785
    49984786        /*Add vx and vy as inputs to the tria element: */
    … …  
    50014789        /*Free ressources:*/
    50024790        xfree((void**)&doflist);
    5003 
    50044791}
    50054792/*}}}*/
    … …  
    50074794void  Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
    50084795
    5009         int i;
    5010 
    5011         const int    numdofpervertex=2;
    5012         const int    numdof=numdofpervertex*NUMVERTICES;
    5013         int*         doflist=NULL;
    5014         double       values[numdof];
    5015         double       lambdax[NUMVERTICES];
    5016         double       lambday[NUMVERTICES];
     4796        const int numdof=NDOF2*NUMVERTICES;
     4797
     4798        int       i;
     4799        int*      doflist=NULL;
     4800        double    values[numdof];
     4801        double    lambdax[NUMVERTICES];
     4802        double    lambday[NUMVERTICES];
    50174803
    50184804        /*Get dof list: */
    … …  
    50204806
    50214807        /*Use the dof list to index into the solution vector: */
    5022         for(i=0;i<numdof;i++){
    5023                 values[i]=solution[doflist[i]];
    5024         }
     4808        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    50254809
    50264810        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    50274811        for(i=0;i<NUMVERTICES;i++){
    5028                 lambdax[i]=values[i*numdofpervertex+0];
    5029                 lambday[i]=values[i*numdofpervertex+1];
     4812                lambdax[i]=values[i*NDOF2+0];
     4813                lambday[i]=values[i*NDOF2+1];
    50304814        }
    50314815
    … …  
    50364820        /*Free ressources:*/
    50374821        xfree((void**)&doflist);
    5038 
    50394822}
    50404823/*}}}*/
    … …  
    50424825void  Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
    50434826       
    5044         int i;
    5045 
    5046         const int    numdofpervertex=2;
    5047         const int    numdof=numdofpervertex*NUMVERTICES;
    5048         int*         doflist=NULL;
    5049         double       values[numdof];
    5050         double       vx[NUMVERTICES];
    5051         double       vy[NUMVERTICES];
    5052         double       vz[NUMVERTICES];
    5053         double       vel[NUMVERTICES];
    5054         double       pressure[NUMVERTICES];
    5055         double       thickness[NUMVERTICES];
    5056         double       rho_ice,g;
    5057         int          dummy;
     4827        const int numdof=NDOF2*NUMVERTICES;
     4828
     4829        int       i;
     4830        int       dummy;
     4831        int*      doflist=NULL;
     4832        double    rho_ice,g;
     4833        double    values[numdof];
     4834        double    vx[NUMVERTICES];
     4835        double    vy[NUMVERTICES];
     4836        double    vz[NUMVERTICES];
     4837        double    vel[NUMVERTICES];
     4838        double    pressure[NUMVERTICES];
     4839        double    thickness[NUMVERTICES];
    50584840       
    50594841        /*Get dof list: */
    … …  
    50614843
    50624844        /*Use the dof list to index into the solution vector: */
    5063         for(i=0;i<numdof;i++){
    5064                 values[i]=solution[doflist[i]];
    5065         }
     4845        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    50664846
    50674847        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    50684848        for(i=0;i<NUMVERTICES;i++){
    5069                 vx[i]=values[i*numdofpervertex+0];
    5070                 vy[i]=values[i*numdofpervertex+1];
    5071         }
    5072 
    5073         /*Get Vz*/
     4849                vx[i]=values[i*NDOF2+0];
     4850                vy[i]=values[i*NDOF2+1];
     4851        }
     4852
     4853        /*Get Vz and compute vel*/
    50744854        GetParameterListOnVertices(&vz[0],VzEnum,0);
    5075 
    5076         /*Now Compute vel*/
    50774855        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);
    50784856
    … …  
    51044882void  Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
    51054883       
    5106         int i;
    5107 
    5108         const int    numdofpervertex=2;
    5109         const int    numdof=numdofpervertex*NUMVERTICES;
    5110         int*         doflist=NULL;
    5111         double       values[numdof];
    5112         double       vx[NUMVERTICES];
    5113         double       vy[NUMVERTICES];
    5114         double       vz[NUMVERTICES];
    5115         double       vel[NUMVERTICES];
    5116         double       pressure[NUMVERTICES];
    5117         double       thickness[NUMVERTICES];
    5118         double       rho_ice,g;
    5119         Input*       vz_input=NULL;
    5120         double*      vz_ptr=NULL;
    5121         int          dummy;
     4884        const int numdof=NDOF2*NUMVERTICES;
     4885       
     4886        int       i;
     4887        int       dummy;
     4888        int*      doflist=NULL;
     4889        double    rho_ice,g;
     4890        double    values[numdof];
     4891        double    vx[NUMVERTICES];
     4892        double    vy[NUMVERTICES];
     4893        double    vz[NUMVERTICES];
     4894        double    vel[NUMVERTICES];
     4895        double    pressure[NUMVERTICES];
     4896        double    thickness[NUMVERTICES];
     4897        double*   vz_ptr=NULL;
     4898        Input*    vz_input=NULL;
    51224899       
    51234900        /*Get dof list: */
    … …  
    51254902
    51264903        /*Use the dof list to index into the solution vector: */
    5127         for(i=0;i<numdof;i++){
    5128                 values[i]=solution[doflist[i]];
    5129         }
     4904        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    51304905
    51314906        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    51324907        for(i=0;i<NUMVERTICES;i++){
    5133                 vx[i]=values[i*numdofpervertex+0];
    5134                 vy[i]=values[i*numdofpervertex+1];
     4908                vx[i]=values[i*NDOF2+0];
     4909                vy[i]=values[i*NDOF2+1];
    51354910        }
    51364911
    … …  
    51564931        g=matpar->GetG();
    51574932        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
    5158        
    5159         for(i=0;i<NUMVERTICES;i++){
    5160                 pressure[i]=rho_ice*g*thickness[i];
    5161         }
     4933        for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
    51624934
    51634935        /*Now, we have to move the previous Vx and Vy inputs  to old
    … …  
    51754947        /*Free ressources:*/
    51764948        xfree((void**)&doflist);
    5177 
    51784949}
    51794950/*}}}*/
    … …  
    51814952void  Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
    51824953
    5183         const int numdofpervertex = 1;
    5184         const int numdof          = numdofpervertex *NUMVERTICES;
    5185         int*         doflist=NULL;
     4954        const int numdof          = NDOF1*NUMVERTICES;
     4955
     4956        int*      doflist=NULL;
    51864957        double    values[numdof];
    51874958
    … …  
    51904961
    51914962        /*Use the dof list to index into the solution vector: */
    5192         for(int i=0;i<numdof;i++){
    5193                 values[i]=solution[doflist[i]];
    5194         }
     4963        for(int i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    51954964
    51964965        /*Add input to the element: */
    … …  
    51994968        /*Free ressources:*/
    52004969        xfree((void**)&doflist);
    5201 
    52024970}
    52034971/*}}}*/
    … …  
    52385006
    52395007        int i;
    5240         double v13[3];
    5241         double v23[3];
     5008        double v13[3],v23[3];
    52425009        double normal[3];
    52435010        double normal_norm;
Note: See TracChangeset for help on using the changeset viewer.