Changeset 5629


Ignore:
Timestamp:
08/31/10 08:40:24 (15 years ago)
Author:
Mathieu Morlighem
Message:

Use GaussTria

Location:
issm/trunk/src/c/objects
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/PentaRef.h

    r5203 r5629  
    4646                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, double* gauss);
    4747                void GetParameterValue(double* pvalue,double* plist,double* gauss);
     48                void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    4849                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, double* gauss);
     50                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    4951
    5052};

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

    r5625 r5629  
    2929        this->parameters=NULL;
    3030        this->results=NULL;
     31
    3132}
    3233/*}}}*/
    … …  
    3435Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels)
    3536        :TriaRef(nummodels)
    36         ,TriaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id
    37                                                                   { //i is the element index
     37        ,TriaHook(nummodels,index+1,iomodel->numberofelements+1){
    3838        /*id: */
    3939        this->id=tria_id;
    … …  
    226226/*}}}*/
    227227/*FUNCTION Tria::Id {{{1*/
    228 int    Tria::Id(){ return id; }
     228int    Tria::Id(){
     229       
     230        return id;
     231
     232}
    229233/*}}}*/
    230234/*FUNCTION Tria::Marshall {{{1*/
    … …  
    356360void  Tria::InputUpdateFromSolution(double* solution){
    357361
     362        /*retrive parameters: */
    358363        int analysis_type;
    359 
    360         /*retrive parameters: */
    361364        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    362365
    363366        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    364         if (analysis_type==DiagnosticHorizAnalysisEnum){
    365                 InputUpdateFromSolutionDiagnosticHoriz( solution);
    366         }
    367         else if (analysis_type==DiagnosticHutterAnalysisEnum){
    368                 InputUpdateFromSolutionDiagnosticHoriz( solution);
    369         }
    370         else if (analysis_type==AdjointHorizAnalysisEnum){
    371                 InputUpdateFromSolutionAdjointHoriz( solution);
    372         }
    373         else if (analysis_type==BedSlopeXAnalysisEnum){
    374                 InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
    375         }
    376         else if (analysis_type==BedSlopeYAnalysisEnum){
    377                 InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
    378         }
    379         else if (analysis_type==SurfaceSlopeXAnalysisEnum){
    380                 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
    381         }
    382         else if (analysis_type==SurfaceSlopeYAnalysisEnum){
    383                 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
    384         }
    385         else if (analysis_type==PrognosticAnalysisEnum){
    386                 InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    387         }
    388         else if (analysis_type==BalancedthicknessAnalysisEnum){
    389                 InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    390         }
    391         else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){
    392                 InputUpdateFromSolutionAdjointBalancedthickness( solution);
    393         }
    394         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    395                 InputUpdateFromSolutionOneDof(solution,VelEnum);
    396         }
    397         else{
    398                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     367        switch(analysis_type){
     368                case DiagnosticHorizAnalysisEnum:
     369                        InputUpdateFromSolutionDiagnosticHoriz( solution);
     370                        break;
     371                case DiagnosticHutterAnalysisEnum:
     372                        InputUpdateFromSolutionDiagnosticHoriz( solution);
     373                        break;
     374                case AdjointHorizAnalysisEnum:
     375                        InputUpdateFromSolutionAdjointHoriz( solution);
     376                        break;
     377                case BedSlopeXAnalysisEnum:
     378                        InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
     379                        break;
     380                case BedSlopeYAnalysisEnum:
     381                        InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
     382                        break;
     383                case SurfaceSlopeXAnalysisEnum:
     384                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
     385                        break;
     386                case SurfaceSlopeYAnalysisEnum:
     387                        InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
     388                        break;
     389                case PrognosticAnalysisEnum:
     390                        InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     391                        break;
     392                case BalancedthicknessAnalysisEnum:
     393                        InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     394                        break;
     395                case AdjointBalancedthicknessAnalysisEnum:
     396                        InputUpdateFromSolutionAdjointBalancedthickness( solution);
     397                        break;
     398                case BalancedvelocitiesAnalysisEnum:
     399                        InputUpdateFromSolutionOneDof(solution,VelEnum);
     400                        break;
     401                default:
     402                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
    399403        }
    400404}
    … …  
    428432
    429433                default:
    430 
    431434                        ISSMERROR("type %i (%s) not implemented yet",type,EnumToString(type));
    432435        }
    … …  
    446449void  Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
    447450       
    448         const int numgrids=3;
    449         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     451        const int numvertices=3;
     452        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    450453        int i,j;
    451454
    … …  
    572575
    573576        int i;
    574         const int numgrids=3;
    575         int doflist[numgrids];
     577        const int numvertices=3;
     578        int doflist[numvertices];
    576579        double value;
    577580
    … …  
    615618void  Tria::ComputeStrainRate(Vec eps){
    616619
    617         int i;
    618         const int numgrids=3;
    619         int doflist[numgrids];
    620         double value;
     620        int       i;
     621        const int numvertices = 3;
     622        int       doflist[numvertices];
     623        double    value;
    621624
    622625        /*plug local pressure values into global pressure vector: */
    … …  
    628631/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
    629632void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
    630 
    631         int analysis_counter;
    632633       
    633634        /*go into parameters and get the analysis_counter: */
     635        int analysis_counter;
    634636        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
    635637
    … …  
    645647/*FUNCTION Tria::Configure {{{1*/
    646648void  Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
    647 
    648         int analysis_counter;
    649649       
    650650        /*go into parameters and get the analysis_counter: */
     651        int analysis_counter;
    651652        parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
    652653
    … …  
    674675double Tria::RegularizeInversion(){
    675676
    676         int i;
    677 
    678         /* output: */
    679         double Jelem=0;
    680 
    681         /* node data: */
    682         const int    numgrids=3;
    683         const int    numdof=2*numgrids;
     677        /*constants: */
     678        const int    numvertices=3;
    684679        const int    NDOF2=2;
    685         double       xyz_list[numgrids][3];
    686 
    687         /* grid data: */
    688         double B[numgrids];
    689 
    690         /* gaussian points: */
    691         int     num_gauss,ig;
    692         double* first_gauss_area_coord  =  NULL;
    693         double* second_gauss_area_coord =  NULL;
    694         double* third_gauss_area_coord  =  NULL;
    695         double* gauss_weights           =  NULL;
    696         double  gauss_weight;
    697         double  gauss_l1l2l3[3];
    698         double  B_gauss;
    699         double  dhdt_gauss;
    700 
    701         /* parameters: */
    702         double  dk[NDOF2];
    703         double  dB[NDOF2];
    704         int    control_type;
    705         double  cm_noisedmp;
    706 
    707         /* Jacobian: */
    708         double Jdet;
    709 
    710         /*inputs: */
    711         bool onwater;
    712         bool shelf;
    713 
    714         /*retrieve inputs :*/
     680        const int    numdof=NDOF2*numvertices;
     681
     682        /* Intermediaries */
     683        int        i,j,ig;
     684        int        control_type;
     685        bool       onwater,shelf;
     686        double     Jelem = 0;
     687        double     cm_noisedmp;
     688        double     Jdet;
     689        double     xyz_list[numvertices][3];
     690        double     dk[NDOF2],dB[NDOF2];
     691        GaussTria *gauss = NULL;
     692
     693        /*retrieve parameters and inputs*/
     694        this->parameters->FindParam(&control_type,ControlTypeEnum);
     695        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    715696        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    716697        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     698        Input* drag_input=inputs->GetInput(DragCoefficientEnum);
     699        Input* dhdt_input=inputs->GetInput(DhDtEnum);
     700        Input* Bbar_input=matice->inputs->GetInput(RheologyBbarEnum);
    717701
    718702        /*If on water, return 0: */
    719         if(onwater)return 0;
    720 
    721         /*retrieve some parameters: */
    722         this->parameters->FindParam(&control_type,ControlTypeEnum);
    723         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    724 
    725         /* Get node coordinates and dof list: */
    726         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    727 
    728         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    729         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    730 
    731         /* Start  looping on the number of gaussian points: */
    732         for (ig=0; ig<num_gauss; ig++){
    733                 /*Pick up the gaussian point: */
    734                 gauss_weight=*(gauss_weights+ig);
    735                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    736                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    737                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    738 
    739                 /* Get Jacobian determinant: */
    740                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     703        if(onwater) return 0;
     704
     705        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
     706
     707        /* Start looping on the number of gaussian points: */
     708        gauss=new GaussTria(2);
     709        for (ig=gauss->begin();ig<gauss->end();ig++){
     710
     711                gauss->GaussPoint(ig);
     712
     713                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    741714
    742715                /*Add Tikhonov regularization term to misfit*/
    743716                if (control_type==DragCoefficientEnum){
    744                         if (shelf){
    745                                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    746                                 Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;
    747 
    748                         }
     717                        if (!shelf) return 0; //only shelf?... TO BE CHECKED
     718                        ISSMASSERT(drag_input);
     719                        drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     720                        Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
    749721                }
    750722                else if (control_type==RheologyBbarEnum){
    751                         matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBbarEnum);
    752                         Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     723                        ISSMASSERT(Bbar_input);
     724                        Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
     725                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
    753726                }
    754727                else if (control_type==DhDtEnum){
    755                         this->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],DhDtEnum);
    756                         Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     728                        ISSMASSERT(dhdt_input);
     729                        dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
     730                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
    757731                }
    758732                else{
    759733                        ISSMERROR("%s%i","unsupported control type: ",control_type);
    760734                }
    761 
    762         }
    763 
    764         xfree((void**)&first_gauss_area_coord);
    765         xfree((void**)&second_gauss_area_coord);
    766         xfree((void**)&third_gauss_area_coord);
    767         xfree((void**)&gauss_weights);
    768 
    769         /*Return: */
     735        }
     736
     737        /*clean-up and return: */
     738        delete gauss;
    770739        return Jelem;
    771740}
    … …  
    774743void  Tria::CreateKMatrix(Mat Kgg){
    775744
     745        /*retrive parameters: */
    776746        int analysis_type;
    777 
    778         /*retrive parameters: */
    779747        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    780748
    … …  
    884852void  Tria::GetBedList(double* bedlist){
    885853       
    886         const int numgrids=3;
    887         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     854        const int numvertices=3;
     855        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    888856       
    889857        inputs->GetParameterValues(bedlist,&gaussgrids[0][0],3,BedEnum);
    … …  
    953921void Tria::GetThicknessList(double* thickness_list){
    954922
    955         const int numgrids=3;
    956         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     923        const int numvertices=3;
     924        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    957925        inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],3,ThicknessEnum);
    958926
    … …  
    1008976
    1009977        /* node data: */
    1010         const int    numgrids=3;
     978        const int    numvertices=3;
    1011979        const int    NDOF1=1;
    1012980        const int    NDOF2=2;
    1013         const int    numdof=NDOF2*numgrids;
    1014         double       xyz_list[numgrids][3];
    1015         int          doflist1[numgrids];
    1016         double       dh1dh3[NDOF2][numgrids];
     981        const int    numdof=NDOF2*numvertices;
     982        double       xyz_list[numvertices][3];
     983        int          doflist1[numvertices];
     984        double       dh1dh3[NDOF2][numvertices];
    1017985
    1018986        /* grid data: */
    1019         double B[numgrids];
     987        double B[numvertices];
    1020988
    1021989
    … …  
    1030998
    1031999        /*element vector at the gaussian points: */
    1032         double  grade_g[numgrids]={0.0};
    1033         double  grade_g_gaussian[numgrids];
     1000        double  grade_g[numvertices]={0.0};
     1001        double  grade_g_gaussian[numvertices];
    10341002
    10351003        /* Jacobian: */
    … …  
    10771045
    10781046        /* Get node coordinates and dof list: */
    1079         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1047        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    10801048        GetDofList1(&doflist1[0]);
    10811049
    … …  
    11281096
    11291097                /*Build gradje_g_gaussian vector (actually -dJ/dB): */
    1130                 for (i=0;i<numgrids;i++){
     1098                for (i=0;i<numvertices;i++){
    11311099                        //standard gradient dJ/dki
    11321100                        grade_g_gaussian[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss_weight*l1l2l3[i];
    … …  
    11481116
    11491117                /*Add grade_g_gaussian to grade_g: */
    1150                 for( i=0; i<numgrids;i++) grade_g[i]+=grade_g_gaussian[i];
     1118                for( i=0; i<numvertices;i++) grade_g[i]+=grade_g_gaussian[i];
    11511119        }
    11521120
    11531121        /*Add grade_g to global vector gradient: */
    1154         VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
     1122        VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
    11551123
    11561124        xfree((void**)&first_gauss_area_coord);
    … …  
    11661134
    11671135        /* node data: */
    1168         const int    numgrids=3;
     1136        const int    numvertices=3;
    11691137        const int    NDOF2=2;
    1170         const int    numdof=NDOF2*numgrids;
    1171         double       xyz_list[numgrids][3];
    1172         int          doflist1[numgrids];
    1173         double       dh1dh3[NDOF2][numgrids];
     1138        const int    numdof=NDOF2*numvertices;
     1139        double       xyz_list[numvertices][3];
     1140        int          doflist1[numvertices];
     1141        double       dh1dh3[NDOF2][numvertices];
    11741142
    11751143        /* gaussian points: */
    … …  
    11931161
    11941162        /*element vector at the gaussian points: */
    1195         double  grade_g[numgrids]={0.0};
    1196         double  grade_g_gaussian[numgrids];
     1163        double  grade_g[numvertices]={0.0};
     1164        double  grade_g_gaussian[numvertices];
    11971165
    11981166        /* Jacobian: */
    … …  
    12011169        /*nodal functions: */
    12021170        double l1l2l3[3];
    1203         double    alpha2complement_list[numgrids];                          //TO BE DELETED
    1204         double    gauss[numgrids][numgrids] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     1171        double    alpha2complement_list[numvertices];                          //TO BE DELETED
     1172        double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    12051173
    12061174        /* strain rate: */
    … …  
    12411209
    12421210        /* Get node coordinates and dof list: */
    1243         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1211        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    12441212        GetDofList1(&doflist1[0]);
    12451213
    … …  
    12491217
    12501218        /*COMPUT alpha2_list (TO BE DELETED)*/
    1251         for(i=0;i<numgrids;i++){
     1219        for(i=0;i<numvertices;i++){
    12521220                friction->GetAlphaComplement(&alpha2complement_list[i],&gauss[i][0],VxEnum,VyEnum);
    12531221        }
    … …  
    13001268
    13011269                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    1302                 for (i=0;i<numgrids;i++){
     1270                for (i=0;i<numvertices;i++){
    13031271
    13041272                        //standard term dJ/dki
    … …  
    13201288               
    13211289                /*Add gradje_g_gaussian vector to gradje_g: */
    1322                 for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
     1290                for( i=0; i<numvertices; i++)grade_g[i]+=grade_g_gaussian[i];
    13231291        }
    13241292
    13251293        /*Add grade_g to global vector gradient: */
    1326         VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
     1294        VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
    13271295
    13281296        xfree((void**)&first_gauss_area_coord);
    … …  
    13401308
    13411309        /* node data: */
    1342         const int    numgrids=3;
     1310        const int    numvertices=3;
    13431311        const int    NDOF1=1;
    1344         int          doflist1[numgrids];
     1312        int          doflist1[numvertices];
    13451313
    13461314        /*Gauss*/
    1347         double  gauss[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     1315        double  gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    13481316
    13491317        /* grid data: */
    … …  
    13511319
    13521320        /*element vector at the gaussian points: */
    1353         double  gradient_g[numgrids];
     1321        double  gradient_g[numvertices];
    13541322
    13551323        /*Inputs*/
    … …  
    13631331
    13641332        /* Start  looping on the vertices: */
    1365         for(i=0; i<numgrids;i++){
     1333        for(i=0; i<numvertices;i++){
    13661334                adjoint_input->GetParameterValue(&lambda,&gauss[i][0]);
    13671335                gradient_g[i]=-lambda;
    … …  
    13691337
    13701338        /*Add grade_g to global vector gradient: */
    1371         VecSetValues(gradient,numgrids,doflist1,(const double*)gradient_g,INSERT_VALUES);
     1339        VecSetValues(gradient,numvertices,doflist1,(const double*)gradient_g,INSERT_VALUES);
    13721340}
    13731341/*}}}*/
    … …  
    15681536        int i;
    15691537
    1570         const int    numgrids=3;
     1538        const int    numvertices=3;
    15711539        const int    numdofs=2;
    15721540        double mass_flux=0;
    1573         double xyz_list[numgrids][3];
     1541        double xyz_list[numvertices][3];
    15741542        double gauss_1[3];
    15751543        double gauss_2[3];
    … …  
    15911559
    15921560        /*Get xyz list: */
    1593         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1561        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    15941562
    15951563        /*get area coordinates of 0 and 1 locations: */
    … …  
    16301598        int i;
    16311599        int dim;
    1632         const int numgrids=3;
    1633         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1634         double  vx_values[numgrids];
     1600        const int numvertices=3;
     1601        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1602        double  vx_values[numvertices];
    16351603        double  maxabsvx;
    16361604
    … …  
    16391607
    16401608        /*retrive velocity values at nodes */
    1641         inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
     1609        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum);
    16421610
    16431611        /*process units if requested: */
    1644         if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters);
     1612        if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters);
    16451613
    16461614        /*now, compute maximum:*/
    16471615        maxabsvx=fabs(vx_values[0]);
    1648         for(i=1;i<numgrids;i++){
     1616        for(i=1;i<numvertices;i++){
    16491617                if (fabs(vx_values[i])>maxabsvx)maxabsvx=fabs(vx_values[i]);
    16501618        }
    … …  
    16591627        int i;
    16601628        int dim;
    1661         const int numgrids=3;
    1662         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1663         double  vy_values[numgrids];
     1629        const int numvertices=3;
     1630        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1631        double  vy_values[numvertices];
    16641632        double  maxabsvy;
    16651633
    … …  
    16681636
    16691637        /*retrive velocity values at nodes */
    1670         inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
     1638        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum);
    16711639
    16721640        /*process units if requested: */
    1673         if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters);
     1641        if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters);
    16741642
    16751643        /*now, compute maximum:*/
    16761644        maxabsvy=fabs(vy_values[0]);
    1677         for(i=1;i<numgrids;i++){
     1645        for(i=1;i<numvertices;i++){
    16781646                if (fabs(vy_values[i])>maxabsvy)maxabsvy=fabs(vy_values[i]);
    16791647        }
    … …  
    16881656        int i;
    16891657        int dim;
    1690         const int numgrids=3;
    1691         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1692         double  vz_values[numgrids];
     1658        const int numvertices=3;
     1659        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1660        double  vz_values[numvertices];
    16931661        double  maxabsvz;
    16941662
    … …  
    16971665
    16981666        /*retrive velocity values at nodes */
    1699         inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
     1667        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum);
    17001668
    17011669        /*process units if requested: */
    1702         if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters);
     1670        if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters);
    17031671
    17041672        /*now, compute maximum:*/
    17051673        maxabsvz=fabs(vz_values[0]);
    1706         for(i=1;i<numgrids;i++){
     1674        for(i=1;i<numvertices;i++){
    17071675                if (fabs(vz_values[i])>maxabsvz)maxabsvz=fabs(vz_values[i]);
    17081676        }
    … …  
    17171685        int i;
    17181686        int dim;
    1719         const int numgrids=3;
    1720         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1721         double  vel_values[numgrids];
     1687        const int numvertices=3;
     1688        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1689        double  vel_values[numvertices];
    17221690        double  maxvel;
    17231691
    … …  
    17261694
    17271695        /*retrive velocity values at nodes */
    1728         inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numgrids,VelEnum);
     1696        inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numvertices,VelEnum);
    17291697
    17301698        /*process units if requested: */
    1731         if(process_units)UnitConversion(&vel_values[0],numgrids,IuToExtEnum,VelEnum,this->parameters);
     1699        if(process_units)UnitConversion(&vel_values[0],numvertices,IuToExtEnum,VelEnum,this->parameters);
    17321700
    17331701        /*now, compute maximum:*/
    17341702        maxvel=vel_values[0];
    1735         for(i=1;i<numgrids;i++){
     1703        for(i=1;i<numvertices;i++){
    17361704                if (vel_values[i]>maxvel)maxvel=vel_values[i];
    17371705        }
    … …  
    17471715        int i;
    17481716        int dim;
    1749         const int numgrids=3;
    1750         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1751         double  vx_values[numgrids];
     1717        const int numvertices=3;
     1718        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1719        double  vx_values[numvertices];
    17521720        double  maxvx;
    17531721
    … …  
    17561724
    17571725        /*retrive velocity values at nodes */
    1758         inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
     1726        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum);
    17591727
    17601728        /*process units if requested: */
    1761         if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters);
     1729        if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters);
    17621730
    17631731        /*now, compute maximum:*/
    17641732        maxvx=vx_values[0];
    1765         for(i=1;i<numgrids;i++){
     1733        for(i=1;i<numvertices;i++){
    17661734                if (vx_values[i]>maxvx)maxvx=vx_values[i];
    17671735        }
    … …  
    17771745        int i;
    17781746        int dim;
    1779         const int numgrids=3;
    1780         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1781         double  vy_values[numgrids];
     1747        const int numvertices=3;
     1748        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1749        double  vy_values[numvertices];
    17821750        double  maxvy;
    17831751
    … …  
    17861754
    17871755        /*retrive velocity values at nodes */
    1788         inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
     1756        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum);
    17891757
    17901758        /*process units if requested: */
    1791         if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters);
     1759        if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters);
    17921760
    17931761        /*now, compute maximum:*/
    17941762        maxvy=vy_values[0];
    1795         for(i=1;i<numgrids;i++){
     1763        for(i=1;i<numvertices;i++){
    17961764                if (vy_values[i]>maxvy)maxvy=vy_values[i];
    17971765        }
    … …  
    18071775        int i;
    18081776        int dim;
    1809         const int numgrids=3;
    1810         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1811         double  vz_values[numgrids];
     1777        const int numvertices=3;
     1778        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1779        double  vz_values[numvertices];
    18121780        double  maxvz;
    18131781
    … …  
    18161784
    18171785        /*retrive velocity values at nodes */
    1818         inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
     1786        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum);
    18191787
    18201788        /*process units if requested: */
    1821         if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters);
     1789        if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters);
    18221790
    18231791        /*now, compute maximum:*/
    18241792        maxvz=vz_values[0];
    1825         for(i=1;i<numgrids;i++){
     1793        for(i=1;i<numvertices;i++){
    18261794                if (vz_values[i]>maxvz)maxvz=vz_values[i];
    18271795        }
    … …  
    18371805        int i;
    18381806        int dim;
    1839         const int numgrids=3;
    1840         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1841         double  vel_values[numgrids];
     1807        const int numvertices=3;
     1808        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1809        double  vel_values[numvertices];
    18421810        double  minvel;
    18431811
    … …  
    18461814
    18471815        /*retrive velocity values at nodes */
    1848         inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numgrids,VelEnum);
     1816        inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numvertices,VelEnum);
    18491817
    18501818        /*process units if requested: */
    1851         if(process_units)UnitConversion(&vel_values[0],numgrids,IuToExtEnum,VelEnum,this->parameters);
     1819        if(process_units)UnitConversion(&vel_values[0],numvertices,IuToExtEnum,VelEnum,this->parameters);
    18521820
    18531821        /*now, compute minimum:*/
    18541822        minvel=vel_values[0];
    1855         for(i=1;i<numgrids;i++){
     1823        for(i=1;i<numvertices;i++){
    18561824                if (vel_values[i]<minvel)minvel=vel_values[i];
    18571825        }
    … …  
    18671835        int i;
    18681836        int dim;
    1869         const int numgrids=3;
    1870         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1871         double  vx_values[numgrids];
     1837        const int numvertices=3;
     1838        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1839        double  vx_values[numvertices];
    18721840        double  minvx;
    18731841
    … …  
    18761844
    18771845        /*retrive velocity values at nodes */
    1878         inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum);
     1846        inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum);
    18791847
    18801848        /*process units if requested: */
    1881         if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters);
     1849        if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters);
    18821850
    18831851        /*now, compute minimum:*/
    18841852        minvx=vx_values[0];
    1885         for(i=1;i<numgrids;i++){
     1853        for(i=1;i<numvertices;i++){
    18861854                if (vx_values[i]<minvx)minvx=vx_values[i];
    18871855        }
    … …  
    18971865        int i;
    18981866        int dim;
    1899         const int numgrids=3;
    1900         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1901         double  vy_values[numgrids];
     1867        const int numvertices=3;
     1868        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1869        double  vy_values[numvertices];
    19021870        double  minvy;
    19031871
    … …  
    19061874
    19071875        /*retrive velocity values at nodes */
    1908         inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum);
     1876        inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum);
    19091877
    19101878        /*process units if requested: */
    1911         if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters);
     1879        if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters);
    19121880
    19131881        /*now, compute minimum:*/
    19141882        minvy=vy_values[0];
    1915         for(i=1;i<numgrids;i++){
     1883        for(i=1;i<numvertices;i++){
    19161884                if (vy_values[i]<minvy)minvy=vy_values[i];
    19171885        }
    … …  
    19271895        int i;
    19281896        int dim;
    1929         const int numgrids=3;
    1930         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    1931         double  vz_values[numgrids];
     1897        const int numvertices=3;
     1898        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
     1899        double  vz_values[numvertices];
    19321900        double  minvz;
    19331901
    … …  
    19361904
    19371905        /*retrive velocity values at nodes */
    1938         inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum);
     1906        inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum);
    19391907
    19401908        /*process units if requested: */
    1941         if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters);
     1909        if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters);
    19421910
    19431911        /*now, compute minimum:*/
    19441912        minvz=vz_values[0];
    1945         for(i=1;i<numgrids;i++){
     1913        for(i=1;i<numvertices;i++){
    19461914                if (vz_values[i]<minvz)minvz=vz_values[i];
    19471915        }
    … …  
    19611929
    19621930        /* node data: */
    1963         const int    numgrids=3;
    1964         const int    numdof=1*numgrids;
     1931        const int    numvertices=3;
     1932        const int    numdof=1*numvertices;
    19651933        const int    NDOF=1;
    1966         double       xyz_list[numgrids][3];
     1934        double       xyz_list[numvertices][3];
    19671935
    19681936        /* gaussian points: */
    … …  
    19741942        double  gauss_weight;
    19751943        double  gauss_l1l2l3[3];
    1976         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     1944        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    19771945
    19781946        /* parameters: */
    … …  
    19951963
    19961964        /* Get node coordinates and dof list: */
    1997         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1965        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    19981966
    19991967        /*Retrieve all inputs we will be needing: */
    … …  
    20432011
    20442012        /* node data: */
    2045         const int    numgrids=3;
    2046         const int    numdof=2*numgrids;
     2013        const int    numvertices=3;
     2014        const int    numdof=2*numvertices;
    20472015        const int    NDOF2=2;
    2048         double       xyz_list[numgrids][3];
     2016        double       xyz_list[numvertices][3];
    20492017
    20502018        /* grid data: */
    2051         double vx_list[numgrids];
    2052         double vy_list[numgrids];
    2053         double obs_vx_list[numgrids];
    2054         double obs_vy_list[numgrids];
    2055         double misfit_square_list[numgrids];
    2056         double misfit_list[numgrids];
    2057         double weights_list[numgrids];
     2019        double vx_list[numvertices];
     2020        double vy_list[numvertices];
     2021        double obs_vx_list[numvertices];
     2022        double obs_vy_list[numvertices];
     2023        double misfit_square_list[numvertices];
     2024        double misfit_list[numvertices];
     2025        double weights_list[numvertices];
    20582026
    20592027        /* gaussian points: */
    … …  
    20652033        double  gauss_weight;
    20662034        double  gauss_l1l2l3[3];
    2067         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2035        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    20682036
    20692037        /* parameters: */
    … …  
    20862054
    20872055        /* Get node coordinates and dof list: */
    2088         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2056        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    20892057
    20902058        /* Recover input data: */
    … …  
    21162084         *
    21172085         */
    2118         for (i=0;i<numgrids;i++){
     2086        for (i=0;i<numvertices;i++){
    21192087                misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2));
    21202088        }
    21212089        /*Process units: */
    2122         if(process_units)UnitConversion(&misfit_list[0],numgrids,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
     2090        if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters);
    21232091
    21242092        /*Apply weights to misfits*/
    2125         for (i=0;i<numgrids;i++){
     2093        for (i=0;i<numvertices;i++){
    21262094                misfit_list[i]=weights_list[i]*misfit_list[i];
    21272095        }
    … …  
    21662134
    21672135        /* node data: */
    2168         const int    numgrids=3;
    2169         const int    numdof=2*numgrids;
     2136        const int    numvertices=3;
     2137        const int    numdof=2*numvertices;
    21702138        const int    NDOF2=2;
    2171         double       xyz_list[numgrids][3];
     2139        double       xyz_list[numvertices][3];
    21722140
    21732141        /* grid data: */
    2174         double vx_list[numgrids];
    2175         double vy_list[numgrids];
    2176         double obs_vx_list[numgrids];
    2177         double obs_vy_list[numgrids];
    2178         double misfit_square_list[numgrids];
    2179         double misfit_list[numgrids];
    2180         double weights_list[numgrids];
     2142        double vx_list[numvertices];
     2143        double vy_list[numvertices];
     2144        double obs_vx_list[numvertices];
     2145        double obs_vy_list[numvertices];
     2146        double misfit_square_list[numvertices];
     2147        double misfit_list[numvertices];
     2148        double weights_list[numvertices];
    21812149
    21822150        /* gaussian points: */
    … …  
    21882156        double  gauss_weight;
    21892157        double  gauss_l1l2l3[3];
    2190         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2158        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    21912159
    21922160        /* parameters: */
    … …  
    22132181
    22142182        /* Get node coordinates and dof list: */
    2215         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2183        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    22162184
    22172185        /* Recover input data: */
    … …  
    22432211         *              obs                        obs                     
    22442212         */
    2245         for (i=0;i<numgrids;i++){
     2213        for (i=0;i<numvertices;i++){
    22462214                scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2);
    22472215                scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2);
    … …  
    22522220
    22532221        /*Process units: */
    2254         if(process_units)UnitConversion(&misfit_list[0],numgrids,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
     2222        if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters);
    22552223
    22562224        /*Apply weights to misfits*/
    2257         for (i=0;i<numgrids;i++){
     2225        for (i=0;i<numvertices;i++){
    22582226                misfit_list[i]=weights_list[i]*misfit_list[i];
    22592227        }
    … …  
    22982266
    22992267        /* node data: */
    2300         const int    numgrids=3;
    2301         const int    numdof=2*numgrids;
     2268        const int    numvertices=3;
     2269        const int    numdof=2*numvertices;
    23022270        const int    NDOF2=2;
    2303         double       xyz_list[numgrids][3];
     2271        double       xyz_list[numvertices][3];
    23042272
    23052273        /* grid data: */
    2306         double vx_list[numgrids];
    2307         double vy_list[numgrids];
    2308         double obs_vx_list[numgrids];
    2309         double obs_vy_list[numgrids];
    2310         double misfit_square_list[numgrids];
    2311         double misfit_list[numgrids];
    2312         double weights_list[numgrids];
     2274        double vx_list[numvertices];
     2275        double vy_list[numvertices];
     2276        double obs_vx_list[numvertices];
     2277        double obs_vy_list[numvertices];
     2278        double misfit_square_list[numvertices];
     2279        double misfit_list[numvertices];
     2280        double weights_list[numvertices];
    23132281
    23142282        /* gaussian points: */
    … …  
    23202288        double  gauss_weight;
    23212289        double  gauss_l1l2l3[3];
    2322         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2290        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    23232291
    23242292        /* parameters: */
    … …  
    23452313
    23462314        /* Get node coordinates and dof list: */
    2347         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2315        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    23482316
    23492317        /* Recover input data: */
    … …  
    23752343         *                            obs
    23762344         */
    2377         for (i=0;i<numgrids;i++){
     2345        for (i=0;i<numvertices;i++){
    23782346                velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil.
    23792347                obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil.
    … …  
    23822350
    23832351        /*Process units: */
    2384         if(process_units)UnitConversion(&misfit_list[0],numgrids,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
     2352        if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters);
    23852353
    23862354        /*Apply weights to misfits*/
    2387         for (i=0;i<numgrids;i++){
     2355        for (i=0;i<numvertices;i++){
    23882356                misfit_list[i]=weights_list[i]*misfit_list[i];
    23892357        }
    … …  
    24282396
    24292397        /* node data: */
    2430         const int    numgrids=3;
    2431         const int    numdof=2*numgrids;
     2398        const int    numvertices=3;
     2399        const int    numdof=2*numvertices;
    24322400        const int    NDOF2=2;
    2433         double       xyz_list[numgrids][3];
     2401        double       xyz_list[numvertices][3];
    24342402
    24352403        /* grid data: */
    2436         double vx_list[numgrids];
    2437         double vy_list[numgrids];
    2438         double obs_vx_list[numgrids];
    2439         double obs_vy_list[numgrids];
    2440         double misfit_square_list[numgrids];
    2441         double misfit_list[numgrids];
    2442         double weights_list[numgrids];
     2404        double vx_list[numvertices];
     2405        double vy_list[numvertices];
     2406        double obs_vx_list[numvertices];
     2407        double obs_vy_list[numvertices];
     2408        double misfit_square_list[numvertices];
     2409        double misfit_list[numvertices];
     2410        double weights_list[numvertices];
    24432411
    24442412        /* gaussian points: */
    … …  
    24502418        double  gauss_weight;
    24512419        double  gauss_l1l2l3[3];
    2452         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2420        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    24532421
    24542422        /* parameters: */
    … …  
    24772445
    24782446        /* Get node coordinates and dof list: */
    2479         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2447        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    24802448
    24812449        /* Recover input data: */
    … …  
    25072475         *                              obs                       obs
    25082476         */
    2509         for (i=0;i<numgrids;i++){
     2477        for (i=0;i<numvertices;i++){
    25102478                misfit_list[i]=0.5*pow(meanvel,(double)2)*(
    25112479                                        pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) +
    … …  
    25142482
    25152483        /*Process units: */
    2516         if(process_units)UnitConversion(&misfit_list[0],numgrids,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
     2484        if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters);
    25172485
    25182486        /*Apply weights to misfits*/
    2519         for (i=0;i<numgrids;i++){
     2487        for (i=0;i<numvertices;i++){
    25202488                misfit_list[i]=weights_list[i]*misfit_list[i];
    25212489        }
    … …  
    25602528
    25612529        /* node data: */
    2562         const int    numgrids=3;
    2563         const int    numdof=2*numgrids;
     2530        const int    numvertices=3;
     2531        const int    numdof=2*numvertices;
    25642532        const int    NDOF2=2;
    2565         double       xyz_list[numgrids][3];
     2533        double       xyz_list[numvertices][3];
    25662534
    25672535        /* grid data: */
    2568         double vx_list[numgrids];
    2569         double vy_list[numgrids];
    2570         double obs_vx_list[numgrids];
    2571         double obs_vy_list[numgrids];
    2572         double misfit_square_list[numgrids];
    2573         double misfit_list[numgrids];
    2574         double weights_list[numgrids];
     2536        double vx_list[numvertices];
     2537        double vy_list[numvertices];
     2538        double obs_vx_list[numvertices];
     2539        double obs_vy_list[numvertices];
     2540        double misfit_square_list[numvertices];
     2541        double misfit_list[numvertices];
     2542        double weights_list[numvertices];
    25752543
    25762544        /* gaussian points: */
    … …  
    25822550        double  gauss_weight;
    25832551        double  gauss_l1l2l3[3];
    2584         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     2552        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    25852553
    25862554        /* parameters: */
    … …  
    26092577
    26102578        /* Get node coordinates and dof list: */
    2611         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2579        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    26122580
    26132581        /* Recover input data: */
    … …  
    26392607         *      S                obs            obs
    26402608         */
    2641         for (i=0;i<numgrids;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
     2609        for (i=0;i<numvertices;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2);
    26422610
    26432611        /*Process units: */
    2644         if(process_units)UnitConversion(&misfit_square_list[0],numgrids,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
     2612        if(process_units)UnitConversion(&misfit_square_list[0],numvertices,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters);
    26452613
    26462614        /*Take the square root, and scale by surface: */
    2647         for (i=0;i<numgrids;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
     2615        for (i=0;i<numvertices;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S;
    26482616
    26492617        /*Apply weights to misfits*/
    2650         for (i=0;i<numgrids;i++){
     2618        for (i=0;i<numvertices;i++){
    26512619                misfit_list[i]=weights_list[i]*misfit_list[i];
    26522620        }
    … …  
    27602728
    27612729        /* node data: */
    2762         int numgrids=3;
    2763         double xyz_list[numgrids][3];
     2730        int numvertices=3;
     2731        double xyz_list[numvertices][3];
    27642732        double v13[3];
    27652733        double v23[3];
    … …  
    27762744
    27772745        /* Get node coordinates and dof list: */
    2778         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2746        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    27792747
    27802748        for (i=0;i<3;i++){
    … …  
    30262994
    30272995        /* node data: */
    3028         const int    numgrids=3;
     2996        const int    numvertices=3;
    30292997        const int    NDOF1=1;
    3030         const int    numdof=NDOF1*numgrids;
    3031         double       xyz_list[numgrids][3];
     2998        const int    numdof=NDOF1*numvertices;
     2999        double       xyz_list[numvertices][3];
    30323000        int*         doflist=NULL;
    30333001
    … …  
    30423010
    30433011        /* matrices: */
    3044         double L[numgrids];
    3045         double B[2][numgrids];
    3046         double Bprime[2][numgrids];
     3012        double L[numvertices];
     3013        double B[2][numvertices];
     3014        double Bprime[2][numvertices];
    30473015        double DL[2][2]={0.0};
    30483016        double DLprime[2][2]={0.0};
    … …  
    30733041
    30743042        /* Get node coordinates and dof list: */
    3075         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3043        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    30763044        GetDofList(&doflist);
    30773045
    … …  
    31923160
    31933161        /* node data: */
    3194         const int    numgrids=3;
     3162        const int    numvertices=3;
    31953163        const int    NDOF1=1;
    3196         const int    numdof=NDOF1*numgrids;
    3197         double       xyz_list[numgrids][3];
     3164        const int    numdof=NDOF1*numvertices;
     3165        double       xyz_list[numvertices][3];
    31983166        int*         doflist=NULL;
    31993167
    … …  
    32083176
    32093177        /* matrices: */
    3210         double B[2][numgrids];
    3211         double Bprime[2][numgrids];
     3178        double B[2][numvertices];
     3179        double Bprime[2][numvertices];
    32123180        double DL[2][2]={0.0};
    32133181        double DLprime[2][2]={0.0};
    … …  
    32273195
    32283196        /* Get node coordinates and dof list: */
    3229         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3197        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    32303198        GetDofList(&doflist);
    32313199        this->parameters->FindParam(&dim,DimEnum);
    … …  
    32933261
    32943262        /* node data: */
    3295         const int    numgrids=3;
     3263        const int    numvertices=3;
    32963264        const int    NDOF1=1;
    3297         const int    numdof=NDOF1*numgrids;
    3298         double       xyz_list[numgrids][3];
     3265        const int    numdof=NDOF1*numvertices;
     3266        double       xyz_list[numvertices][3];
    32993267        int*         doflist=NULL;
    3300         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     3268        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    33013269
    33023270        /* gaussian points: */
    … …  
    33103278
    33113279        /* matrices: */
    3312         double L[numgrids];
    3313         double B[2][numgrids];
    3314         double Bprime[2][numgrids];
     3280        double L[numvertices];
     3281        double B[2][numvertices];
     3282        double Bprime[2][numvertices];
    33153283        double DL[2][2]={0.0};
    33163284        double DLprime[2][2]={0.0};
    … …  
    33473315
    33483316        /* Get node coordinates and dof list: */
    3349         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3317        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    33503318        GetDofList(&doflist);
    33513319
    … …  
    33633331        /*Modify z so that it reflects the surface*/
    33643332        surface_input->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3);
    3365         for(i=0;i<numgrids;i++) xyz_list[i][2]=surface_list[i];
     3333        for(i=0;i<numvertices;i++) xyz_list[i][2]=surface_list[i];
    33663334
    33673335        /*Get normal vector to the surface*/
    … …  
    34733441
    34743442        /* node data: */
    3475         const int    numgrids=3;
    3476         const int    numdof=2*numgrids;
    3477         double       xyz_list[numgrids][3];
     3443        const int    numvertices=3;
     3444        const int    numdof=2*numvertices;
     3445        double       xyz_list[numvertices][3];
    34783446        int*         doflist=NULL;
    34793447
    … …  
    35333501
    35343502        /* Get node coordinates and dof list: */
    3535         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3503        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    35363504        GetDofList(&doflist,MacAyealApproximationEnum);
    35373505
    … …  
    40414009        int    i;
    40424010        int    connectivity;
    4043         const int numgrids=3;
     4011        const int numvertices=3;
    40444012        const int NDOF2=2;
    4045         const int numdofs=numgrids*NDOF2;
     4013        const int numdofs=numvertices*NDOF2;
    40464014        int*         doflist=NULL;
    40474015
    … …  
    40784046
    40794047        /* node data: */
    4080         const int    numgrids=3;
     4048        const int    numvertices=3;
    40814049        const int    NDOF1=1;
    4082         const int    numdof=NDOF1*numgrids;
    4083         double       xyz_list[numgrids][3];
     4050        const int    numdof=NDOF1*numvertices;
     4051        double       xyz_list[numvertices][3];
    40844052        int*         doflist=NULL;
    40854053
    … …  
    41144082
    41154083        /* Get node coordinates and dof list: */
    4116         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4084        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    41174085        GetDofList(&doflist);
    41184086
    … …  
    41954163        int i,j;
    41964164
    4197         const int  numgrids=3;
     4165        const int  numvertices=3;
    41984166        const int  NDOF1=1;
    4199         const int  numdof=numgrids*NDOF1;
     4167        const int  numdof=numvertices*NDOF1;
    42004168        int*       doflist=NULL;
    42014169
    42024170        /*Grid data: */
    4203         double     xyz_list[numgrids][3];
     4171        double     xyz_list[numvertices][3];
    42044172
    42054173        /*Material constants */
    … …  
    42284196
    42294197        /* Get node coordinates and dof list: */
    4230         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4198        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    42314199        GetDofList(&doflist);
    42324200
    … …  
    42544222                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
    42554223
    4256                 for(i=0;i<numgrids;i++){
    4257                         for(j=0;j<numgrids;j++){
     4224                for(i=0;i<numvertices;i++){
     4225                        for(j=0;j<numvertices;j++){
    42584226                                K_terms[i][j]+=Ke_gaussian[i][j];
    42594227                        }
    … …  
    42784246
    42794247        /* node data: */
    4280         const int    numgrids=3;
     4248        const int    numvertices=3;
    42814249        const int    NDOF1=1;
    4282         const int    numdof=NDOF1*numgrids;
    4283         double       xyz_list[numgrids][3];
     4250        const int    numdof=NDOF1*numvertices;
     4251        double       xyz_list[numvertices][3];
    42844252        int*         doflist=NULL;
    42854253        int          numberofdofspernode=1;
    … …  
    42954263
    42964264        /* matrices: */
    4297         double L[numgrids];
    4298         double B[2][numgrids];
    4299         double Bprime[2][numgrids];
     4265        double L[numvertices];
     4266        double B[2][numvertices];
     4267        double Bprime[2][numvertices];
    43004268        double DL[2][2]={0.0};
    43014269        double DLprime[2][2]={0.0};
    … …  
    43304298
    43314299        /* Get node coordinates and dof list: */
    4332         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4300        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    43334301        GetDofList(&doflist);
    43344302
    … …  
    44594427
    44604428        /* node data: */
    4461         const int    numgrids=3;
     4429        const int    numvertices=3;
    44624430        const int    NDOF1=1;
    4463         const int    numdof=NDOF1*numgrids;
    4464         double       xyz_list[numgrids][3];
     4431        const int    numdof=NDOF1*numvertices;
     4432        double       xyz_list[numvertices][3];
    44654433        int*         doflist=NULL;
    44664434        int          numberofdofspernode=1;
    … …  
    44764444
    44774445        /* matrices: */
    4478         double L[numgrids];
    4479         double B[2][numgrids];
    4480         double Bprime[2][numgrids];
     4446        double L[numvertices];
     4447        double B[2][numvertices];
     4448        double Bprime[2][numvertices];
    44814449        double DL[2][2]={0.0};
    44824450        double DLprime[2][2]={0.0};
    … …  
    45024470
    45034471        /* Get node coordinates and dof list: */
    4504         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4472        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    45054473        GetDofList(&doflist);
    45064474
    … …  
    46324600       
    46334601        /* node data: */
    4634         const int    numgrids=3;
     4602        const int    numvertices=3;
    46354603        const int    NDOF1=1;
    4636         const int    numdof=NDOF1*numgrids;
    4637         double       xyz_list[numgrids][3];
     4604        const int    numdof=NDOF1*numvertices;
     4605        double       xyz_list[numvertices][3];
    46384606        int*         doflist=NULL;
    46394607
    … …  
    46564624        double  K_terms[numdof][numdof]={0.0};
    46574625        double  Ke_gaussian[numdof][numdof]={0.0};
    4658         double  l1l2l3[numgrids];
     4626        double  l1l2l3[numvertices];
    46594627        double     tl1l2l3D[3];
    46604628        double  D_scalar;
    … …  
    46674635
    46684636        /* Get node coordinates and dof list: */
    4669         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4637        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    46704638        GetDofList(&doflist);
    46714639
    … …  
    47264694
    47274695        /* node data: */
    4728         const int    numgrids=3;
     4696        const int    numvertices=3;
    47294697        const int    NDOF1=1;
    4730         const int    numdof=NDOF1*numgrids;
    4731         double       xyz_list[numgrids][3];
     4698        const int    numdof=NDOF1*numvertices;
     4699        double       xyz_list[numvertices][3];
    47324700        int*         doflist=NULL;
    47334701        int          numberofdofspernode=1;
    … …  
    47434711
    47444712        /* matrix */
    4745         double pe_g[numgrids]={0.0};
    4746         double L[numgrids];
     4713        double pe_g[numvertices]={0.0};
     4714        double L[numvertices];
    47474715        double Jdettria;
    47484716
    … …  
    47584726
    47594727        /* Get node coordinates and dof list: */
    4760         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4728        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    47614729        GetDofList(&doflist);
    47624730
    … …  
    48104778
    48114779        /* node data: */
    4812         const int    numgrids=3;
     4780        const int    numvertices=3;
    48134781        const int    NDOF1=1;
    4814         const int    numdof=NDOF1*numgrids;
    4815         double       xyz_list[numgrids][3];
     4782        const int    numdof=NDOF1*numvertices;
     4783        double       xyz_list[numvertices][3];
    48164784        int*         doflist=NULL;
    48174785        int          numberofdofspernode=1;
    … …  
    48274795
    48284796        /* matrix */
    4829         double pe_g[numgrids]={0.0};
    4830         double L[numgrids];
     4797        double pe_g[numvertices]={0.0};
     4798        double L[numvertices];
    48314799        double Jdettria;
    48324800
    … …  
    48424810
    48434811        /* Get node coordinates and dof list: */
    4844         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4812        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    48454813        GetDofList(&doflist);
    48464814
    … …  
    48944862
    48954863        /* node data: */
    4896         const int    numgrids=3;
     4864        const int    numvertices=3;
    48974865        const int    NDOF1=1;
    4898         const int    numdof=NDOF1*numgrids;
    4899         double       xyz_list[numgrids][3];
     4866        const int    numdof=NDOF1*numvertices;
     4867        double       xyz_list[numvertices][3];
    49004868        int*         doflist=NULL;
    49014869        int          numberofdofspernode=1;
    … …  
    49114879
    49124880        /* matrix */
    4913         double pe_g[numgrids]={0.0};
    4914         double L[numgrids];
     4881        double pe_g[numvertices]={0.0};
     4882        double L[numvertices];
    49154883        double Jdettria;
    49164884
    … …  
    49244892
    49254893        /* Get node coordinates and dof list: */
    4926         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4894        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    49274895        GetDofList(&doflist);
    49284896
    … …  
    49744942
    49754943        /* node data: */
    4976         const int    numgrids=3;
     4944        const int    numvertices=3;
    49774945        const int    NDOF1=1;
    4978         const int    numdof=NDOF1*numgrids;
    4979         double       xyz_list[numgrids][3];
     4946        const int    numdof=NDOF1*numvertices;
     4947        double       xyz_list[numvertices][3];
    49804948        int*         doflist=NULL;
    49814949
    … …  
    50004968
    50014969        /* matrices: */
    5002         double L[numgrids];
     4970        double L[numvertices];
    50034971
    50044972        /*input parameters for structural analysis (diagnostic): */
    … …  
    50154983
    50164984        /* Get node coordinates and dof list: */
    5017         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4985        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    50184986        GetDofList(&doflist);
    50194987
    … …  
    50565024
    50575025                /*Build gaussian vector: */
    5058                 for(i=0;i<numgrids;i++){
     5026                for(i=0;i<numvertices;i++){
    50595027                        pe_g_gaussian[i]=-Jdet*gauss_weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
    50605028                }
    … …  
    50825050
    50835051        /* node data: */
    5084         const int    numgrids=3;
    5085         const int    numdof=2*numgrids;
     5052        const int    numvertices=3;
     5053        const int    numdof=2*numvertices;
    50865054        const int    NDOF2=2;
    5087         double       xyz_list[numgrids][3];
     5055        double       xyz_list[numvertices][3];
    50885056        int*         doflist=NULL;
    50895057       
    … …  
    51335101
    51345102        /* Get node coordinates and dof list: */
    5135         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5103        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    51365104        GetDofList(&doflist,MacAyealApproximationEnum);
    51375105
    … …  
    51695137                /*Build pe_g_gaussian vector: */
    51705138                if(drag_type==1){
    5171                         for (i=0;i<numgrids;i++){
     5139                        for (i=0;i<numvertices;i++){
    51725140                                for (j=0;j<NDOF2;j++){
    51735141                                        pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i];
    … …  
    51765144                }
    51775145                else {
    5178                         for (i=0;i<numgrids;i++){
     5146                        for (i=0;i<numvertices;i++){
    51795147                                for (j=0;j<NDOF2;j++){
    51805148                                        pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
    … …  
    52055173
    52065174        /* node data: */
    5207         const int    numgrids=3;
    5208         const int    numdof=1*numgrids;
    5209         double       xyz_list[numgrids][3];
     5175        const int    numvertices=3;
     5176        const int    numdof=1*numvertices;
     5177        double       xyz_list[numvertices][3];
    52105178        int*         doflist=NULL;
    52115179
    … …  
    52385206
    52395207        /* Get node coordinates and dof list: */
    5240         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5208        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    52415209        GetDofList(&doflist);
    52425210
    … …  
    52685236                weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]);
    52695237
    5270                 for (i=0;i<numgrids;i++){
     5238                for (i=0;i<numvertices;i++){
    52715239                        pe_g_gaussian[i]=(thicknessobs-thickness)*weight*Jdet*gauss_weight*l1l2l3[i];
    52725240                }
    … …  
    52955263
    52965264        /* node data: */
    5297         const int    numgrids=3;
    5298         const int    numdof=2*numgrids;
     5265        const int    numvertices=3;
     5266        const int    numdof=2*numvertices;
    52995267        const int    NDOF2=2;
    5300         double       xyz_list[numgrids][3];
     5268        double       xyz_list[numvertices][3];
    53015269        int*         doflist=NULL;
    5302         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     5270        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    53035271
    53045272        /* grid data: */
    5305         double vx_list[numgrids];
    5306         double vy_list[numgrids];
    5307         double obs_vx_list[numgrids];
    5308         double obs_vy_list[numgrids];
    5309         double dux_list[numgrids];
    5310         double duy_list[numgrids];
    5311         double weights_list[numgrids];
     5273        double vx_list[numvertices];
     5274        double vy_list[numvertices];
     5275        double obs_vx_list[numvertices];
     5276        double obs_vy_list[numvertices];
     5277        double dux_list[numvertices];
     5278        double duy_list[numvertices];
     5279        double weights_list[numvertices];
    53125280
    53135281        /* gaussian points: */
    … …  
    53475315
    53485316        /* Get node coordinates and dof list: */
    5349         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5317        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    53505318        GetDofList(&doflist);
    53515319
    … …  
    53885356                 *        du     obs
    53895357                 */
    5390                 for (i=0;i<numgrids;i++){
     5358                for (i=0;i<numvertices;i++){
    53915359                        dux_list[i]=obs_vx_list[i]-vx_list[i];
    53925360                        duy_list[i]=obs_vy_list[i]-vy_list[i];
    … …  
    54065374                 *               obs
    54075375                 */
    5408                 for (i=0;i<numgrids;i++){
     5376                for (i=0;i<numvertices;i++){
    54095377                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
    54105378                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    … …  
    54285396                 *           
    54295397                 */
    5430                 for (i=0;i<numgrids;i++){
     5398                for (i=0;i<numvertices;i++){
    54315399                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
    54325400                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
    … …  
    54475415                 *        du      S  2 sqrt(...)           obs
    54485416                 */
    5449                 for (i=0;i<numgrids;i++){
     5417                for (i=0;i<numvertices;i++){
    54505418                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    54515419                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    … …  
    54645432                 *        du                         |u| + eps  |u|                           u + eps
    54655433                 */
    5466                 for (i=0;i<numgrids;i++){
     5434                for (i=0;i<numvertices;i++){
    54675435                        dux_list[i] = - pow(meanvel,(double)2)*(
    54685436                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
    … …  
    54775445
    54785446        /*Apply weights to DU*/
    5479         for (i=0;i<numgrids;i++){
     5447        for (i=0;i<numvertices;i++){
    54805448                dux_list[i]=weights_list[i]*dux_list[i];
    54815449                duy_list[i]=weights_list[i]*duy_list[i];
    … …  
    55065474
    55075475                /*compute Du*/
    5508                 for (i=0;i<numgrids;i++){
     5476                for (i=0;i<numvertices;i++){
    55095477                        pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i];
    55105478                        pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i];
    … …  
    55345502
    55355503        /* node data: */
    5536         const int    numgrids=3;
     5504        const int    numvertices=3;
    55375505        const int    NDOF4=4;
    5538         const int    numdof=NDOF4*numgrids;
    5539         double       xyz_list[numgrids][3];
     5506        const int    numdof=NDOF4*numvertices;
     5507        double       xyz_list[numvertices][3];
    55405508        int*         doflist=NULL;
    5541         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     5509        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    55425510
    55435511        /* grid data: */
    5544         double vx_list[numgrids];
    5545         double vy_list[numgrids];
    5546         double obs_vx_list[numgrids];
    5547         double obs_vy_list[numgrids];
    5548         double dux_list[numgrids];
    5549         double duy_list[numgrids];
    5550         double weights_list[numgrids];
     5512        double vx_list[numvertices];
     5513        double vy_list[numvertices];
     5514        double obs_vx_list[numvertices];
     5515        double obs_vy_list[numvertices];
     5516        double dux_list[numvertices];
     5517        double duy_list[numvertices];
     5518        double weights_list[numvertices];
    55515519
    55525520        /* gaussian points: */
    … …  
    55865554
    55875555        /* Get node coordinates and dof list: */
    5588         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5556        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    55895557        GetDofList(&doflist);
    55905558
    … …  
    56275595                 *        du     obs
    56285596                 */
    5629                 for (i=0;i<numgrids;i++){
     5597                for (i=0;i<numvertices;i++){
    56305598                        dux_list[i]=obs_vx_list[i]-vx_list[i];
    56315599                        duy_list[i]=obs_vy_list[i]-vy_list[i];
    … …  
    56455613                 *               obs
    56465614                 */
    5647                 for (i=0;i<numgrids;i++){
     5615                for (i=0;i<numvertices;i++){
    56485616                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
    56495617                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    … …  
    56675635                 *           
    56685636                 */
    5669                 for (i=0;i<numgrids;i++){
     5637                for (i=0;i<numvertices;i++){
    56705638                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
    56715639                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
    … …  
    56865654                 *        du      S  2 sqrt(...)           obs
    56875655                 */
    5688                 for (i=0;i<numgrids;i++){
     5656                for (i=0;i<numvertices;i++){
    56895657                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    56905658                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    … …  
    57035671                 *        du                         |u| + eps  |u|                           u + eps
    57045672                 */
    5705                 for (i=0;i<numgrids;i++){
     5673                for (i=0;i<numvertices;i++){
    57065674                        dux_list[i] = - pow(meanvel,(double)2)*(
    57075675                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
    … …  
    57165684
    57175685        /*Apply weights to DU*/
    5718         for (i=0;i<numgrids;i++){
     5686        for (i=0;i<numvertices;i++){
    57195687                dux_list[i]=weights_list[i]*dux_list[i];
    57205688                duy_list[i]=weights_list[i]*duy_list[i];
    … …  
    57455713
    57465714                /*compute Du*/
    5747                 for (i=0;i<numgrids;i++){
     5715                for (i=0;i<numvertices;i++){
    57485716                        pe_g[i*NDOF4+0]+=dux*Jdet*gauss_weight*l1l2l3[i];
    57495717                        pe_g[i*NDOF4+1]+=duy*Jdet*gauss_weight*l1l2l3[i];
    … …  
    57675735        /*Collapsed formulation: */
    57685736        int       i;
    5769         const int numgrids=3;
     5737        const int numvertices=3;
    57705738        const int NDOF2=2;
    5771         const int numdofs=NDOF2*numgrids;
     5739        const int numdofs=NDOF2*numvertices;
    57725740        int*         doflist=NULL;
    57735741        double    constant_part,ub,vb;
    … …  
    57845752        Input* slopey_input=NULL;
    57855753        Input* thickness_input=NULL;
    5786         double gauss[numgrids][numgrids] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     5754        double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    57875755
    57885756        /*recover some inputs: */
    … …  
    58405808
    58415809        /* node data: */
    5842         const int    numgrids=3;
     5810        const int    numvertices=3;
    58435811        const int    NDOF1=1;
    5844         const int    numdof=NDOF1*numgrids;
    5845         double       xyz_list[numgrids][3];
     5812        const int    numdof=NDOF1*numvertices;
     5813        double       xyz_list[numvertices][3];
    58465814        int*         doflist=NULL;
    58475815        int          numberofdofspernode=1;
    … …  
    58575825
    58585826        /* matrix */
    5859         double pe_g[numgrids]={0.0};
    5860         double L[numgrids];
     5827        double pe_g[numvertices]={0.0};
     5828        double L[numvertices];
    58615829        double Jdettria;
    58625830
    … …  
    58765844
    58775845        /* Get node coordinates and dof list: */
    5878         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5846        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    58795847        GetDofList(&doflist);
    58805848
    … …  
    59305898
    59315899        /* node data: */
    5932         const int    numgrids=3;
     5900        const int    numvertices=3;
    59335901        const int    NDOF1=1;
    5934         const int    numdof=NDOF1*numgrids;
    5935         double       xyz_list[numgrids][3];
     5902        const int    numdof=NDOF1*numvertices;
     5903        double       xyz_list[numvertices][3];
    59365904        int*         doflist=NULL;
    59375905        int          numberofdofspernode=1;
    … …  
    59475915
    59485916        /* matrix */
    5949         double pe_g[numgrids]={0.0};
    5950         double L[numgrids];
     5917        double pe_g[numvertices]={0.0};
     5918        double L[numvertices];
    59515919        double Jdettria;
    59525920
    … …  
    59645932
    59655933        /* Get node coordinates and dof list: */
    5966         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     5934        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    59675935        GetDofList(&doflist);
    59685936
    … …  
    60165984
    60175985        /* node data: */
    6018         const int    numgrids=3;
     5986        const int    numvertices=3;
    60195987        const int    NDOF1=1;
    6020         const int    numdof=NDOF1*numgrids;
    6021         double       xyz_list[numgrids][3];
     5988        const int    numdof=NDOF1*numvertices;
     5989        double       xyz_list[numvertices][3];
    60225990        int*         doflist=NULL;
    60235991       
    … …  
    60506018
    60516019        /* Get node coordinates and dof list: */
    6052         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6020        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    60536021        GetDofList(&doflist);
    60546022
    … …  
    61096077        int i,found;
    61106078       
    6111         const int  numgrids=3;
     6079        const int  numvertices=3;
    61126080        const int  NDOF1=1;
    6113         const int  numdof=numgrids*NDOF1;
     6081        const int  numdof=numvertices*NDOF1;
    61146082        int*         doflist=NULL;
    6115         double       xyz_list[numgrids][3];
     6083        double       xyz_list[numvertices][3];
    61166084
    61176085        double mixed_layer_capacity;
    … …  
    61396107        double  Jdet;
    61406108        double  P_terms[numdof]={0.0};
    6141         double  l1l2l3[numgrids];
     6109        double  l1l2l3[numvertices];
    61426110
    61436111        double  t_pmp;
    … …  
    61486116
    61496117        /* Get node coordinates and dof list: */
    6150         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6118        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    61516119        GetDofList(&doflist);
    61526120
    … …  
    62146182        int i,found;
    62156183       
    6216         const int  numgrids=3;
     6184        const int  numvertices=3;
    62176185        const int  NDOF1=1;
    6218         const int  numdof=numgrids*NDOF1;
     6186        const int  numdof=numvertices*NDOF1;
    62196187        int*       doflist=NULL;
    6220         double     xyz_list[numgrids][3];
     6188        double     xyz_list[numvertices][3];
    62216189
    62226190        double rho_ice;
    … …  
    62326200        double alpha2,vx,vy;
    62336201        double geothermalflux_value;
    6234         double    alpha2_list[numgrids];                                       //TO BE DELETED
    6235         double    gauss[numgrids][numgrids] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    6236         double vx_list[numgrids]; //TO BE DELETED
    6237         double vy_list[numgrids]; //TO BE DELETED
    6238         double basalfriction_list[numgrids]; //TO BE DELETED
     6202        double    alpha2_list[numvertices];                                       //TO BE DELETED
     6203        double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     6204        double vx_list[numvertices]; //TO BE DELETED
     6205        double vy_list[numvertices]; //TO BE DELETED
     6206        double basalfriction_list[numvertices]; //TO BE DELETED
    62396207
    62406208        /* gaussian points: */
    … …  
    62506218        double  Jdet;
    62516219        double  P_terms[numdof]={0.0};
    6252         double  l1l2l3[numgrids];
     6220        double  l1l2l3[numvertices];
    62536221        double  scalar;
    62546222
    … …  
    62656233       
    62666234        /* Get node coordinates and dof list: */
    6267         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6235        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    62686236        GetDofList(&doflist);
    62696237
    … …  
    62876255
    62886256        /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/
    6289         for(i=0;i<numgrids;i++){
     6257        for(i=0;i<numvertices;i++){
    62906258                friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
    62916259        }
    62926260        vx_input->GetParameterValues(&vx_list[0],&gauss[0][0],3); //TO BE DELETED
    62936261        vy_input->GetParameterValues(&vy_list[0],&gauss[0][0],3); //TO BE DELETED
    6294         for(i=0;i<numgrids;i++){
     6262        for(i=0;i<numvertices;i++){
    62956263                basalfriction_list[i]=alpha2_list[i]*(pow(vx_list[i],(double)2.0)+pow(vy_list[i],(double)2.0)); //TO BE DELETED
    62966264        }
    … …  
    63476315
    63486316        double area=0;
    6349         const int    numgrids=3;
    6350         double xyz_list[numgrids][3];
     6317        const int    numvertices=3;
     6318        double xyz_list[numvertices][3];
    63516319        double x1,y1,x2,y2,x3,y3;
    63526320
    63536321        /*Get xyz list: */
    6354         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6322        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    63556323        x1=xyz_list[0][0]; y1=xyz_list[0][1];
    63566324        x2=xyz_list[1][0]; y2=xyz_list[1][1];
    … …  
    63656333        /*Intermediaries*/
    63666334        double    area = 0;
    6367         const int numgrids = 3;
    6368         double    xyz_list[numgrids][3];
     6335        const int numvertices = 3;
     6336        double    xyz_list[numvertices][3];
    63696337        double    x1,y1,x2,y2,x3,y3;
    63706338
    … …  
    63736341
    63746342        /*Get xyz list: */
    6375         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6343        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    63766344        x1=xyz_list[0][0]; y1=xyz_list[0][1];
    63776345        x2=xyz_list[1][0]; y2=xyz_list[1][1];
    … …  
    67356703
    67366704        /* node data: */
    6737         const int    numgrids=3;
     6705        const int    numvertices=3;
    67386706        const int    NDOF2=2;
    6739         double       xyz_list[numgrids][3];
    6740         int          doflist1[numgrids];
    6741         double       dh1dh3[NDOF2][numgrids];
     6707        double       xyz_list[numvertices][3];
     6708        int          doflist1[numvertices];
     6709        double       dh1dh3[NDOF2][numvertices];
    67426710
    67436711        /* grid data: */
    … …  
    67546722        double  gauss_weight;
    67556723        double  gauss_l1l2l3[3];
    6756         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     6724        double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    67576725
    67586726        /* parameters: */
    … …  
    67656733
    67666734        /*element vector at the gaussian points: */
    6767         double  grade_g[numgrids]={0.0};
    6768         double  grade_g_gaussian[numgrids];
     6735        double  grade_g[numvertices]={0.0};
     6736        double  grade_g_gaussian[numvertices];
    67696737
    67706738        /* Jacobian: */
    … …  
    68086776
    68096777        /* Get node coordinates and dof list: */
    6810         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     6778        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    68116779        GetDofList1(&doflist1[0]);
    68126780
    … …  
    68616829
    68626830                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    6863                 for (i=0;i<numgrids;i++){
     6831                for (i=0;i<numvertices;i++){
    68646832                        //standard gradient dJ/dki
    68656833                        grade_g_gaussian[i]=(
    … …  
    68846852
    68856853                /*Add gradje_g_gaussian vector to gradje_g: */
    6886                 for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
     6854                for( i=0; i<numvertices; i++)grade_g[i]+=grade_g_gaussian[i];
    68876855        }
    68886856
    68896857        /*Add grade_g to global vector gradient: */
    6890         VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
     6858        VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
    68916859
    68926860        /*Add grade_g to the inputs of this element: */

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

    r5625 r5629  
    359359}
    360360/*}}}*/
    361 /*FUNCTION TriaRef::GetJacobianInvert {{{1*/
     361/*FUNCTION TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss) {{{1*/
    362362void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
     363
     364        /*Jacobian*/
     365        double J[2][2];
     366
     367        /*Call Jacobian routine to get the jacobian:*/
     368        GetJacobian(&J[0][0], xyz_list, gauss);
     369
     370        /*Invert Jacobian matrix: */
     371        Matrix2x2Invert(Jinv,&J[0][0]);
     372
     373}
     374/*}}}*/
     375/*FUNCTION TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss) {{{1*/
     376void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss){
    363377
    364378        /*Jacobian*/
    … …  
    393407}
    394408/*}}}*/
    395 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives {{{1*/
     409/*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss) {{{1*/
    396410void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
    397411
    … …  
    422436}
    423437/*}}}*/
    424 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference {{{1*/
     438/*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss) {{{1*/
     439void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){
     440
     441        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     442         * actual coordinate system): */
     443        int       i;
     444        const int NDOF2    = 2;
     445        const int numgrids = 3;
     446        double    dh1dh3_ref[NDOF2][numgrids];
     447        double    Jinv[NDOF2][NDOF2];
     448
     449        /*Get derivative values with respect to parametric coordinate system: */
     450        GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss);
     451
     452        /*Get Jacobian invert: */
     453        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
     454
     455        /*Build dh1dh3:
     456         *
     457         * [dhi/dx]= Jinv*[dhi/dr]
     458         * [dhi/dy]       [dhi/ds]
     459         */
     460        for (i=0;i<numgrids;i++){
     461                dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
     462                dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
     463        }
     464
     465}
     466/*}}}*/
     467/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss) {{{1*/
    425468void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss){
    426469
    … …  
    445488}
    446489/*}}}*/
    447 /*FUNCTION TriaRef::GetParameterDerivativeValue {{{1*/
     490/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss) {{{1*/
     491void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){
     492
     493        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     494         * natural coordinate system) at the gaussian point. */
     495
     496        const int NDOF2=2;
     497        const int numgrids=3;
     498
     499        /*First nodal function: */
     500        *(dl1dl3+numgrids*0+0)=-0.5;
     501        *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
     502
     503        /*Second nodal function: */
     504        *(dl1dl3+numgrids*0+1)=0.5;
     505        *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
     506
     507        /*Third nodal function: */
     508        *(dl1dl3+numgrids*0+2)=0;
     509        *(dl1dl3+numgrids*1+2)=1.0/SQRT3;
     510
     511}
     512/*}}}*/
     513/*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss) {{{1*/
    448514void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss){
    449515
    … …  
    468534}
    469535/*}}}*/
    470 /*FUNCTION TriaRef::GetParameterValue{{{1*/
     536/*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss) {{{1*/
     537void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){
     538
     539        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
     540         * point specified by gauss_l1l2l3:
     541         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
     542         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
     543         *
     544         * p is a vector of size 2x1 already allocated.
     545         */
     546
     547        /*Nodal Derivatives*/
     548        double dh1dh3[2][3]; //nodal derivative functions in actual coordinate system.
     549
     550        /*Get dh1dh2dh3 in actual coordinate system: */
     551        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss);
     552
     553        /*Assign values*/
     554        *(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
     555        *(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
     556
     557}
     558/*}}}*/
     559/*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, double* gauss){{{1*/
    471560void TriaRef::GetParameterValue(double* p, double* plist, double* gauss){
    472561
    … …  
    484573}
    485574/*}}}*/
     575/*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){{{1*/
     576void TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){
     577
     578        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
     579         * point specifie by gauss: */
     580
     581        /*nodal functions annd output: */
     582        double l1l2l3[3];
     583
     584        /*Get nodal functions*/
     585        GetNodalFunctions(l1l2l3, gauss);
     586
     587        /*Get parameter*/
     588        *p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];
     589}
     590/*}}}*/

  • issm/trunk/src/c/objects/Elements/TriaRef.h

    r5625 r5629  
    3737                void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
    3838                void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
     39                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss);
    3940                void GetNodalFunctions(double* l1l2l3, double* gauss);
    4041                void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
    4142                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
     43                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss);
    4244                void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);
     45                void GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss);
    4346                void GetParameterValue(double* pp, double* plist, double* gauss);
     47                void GetParameterValue(double* pp, double* plist, GaussTria* gauss);
    4448                void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, double* gauss);
     49                void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, GaussTria* gauss);
    4550
    4651};

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r5103 r5629  
    169169void BoolInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    170170/*}}}*/
     171/*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     172void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
     173/*}}}*/
    171174/*FUNCTION BoolInput::GetParameterValues{{{1*/
    172175void BoolInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    173176/*}}}*/
    174 /*FUNCTION BoolInput::GetParameterDerivativeValue{{{1*/
     177/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    175178void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
     179/*}}}*/
     180/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
     181void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
    176182/*}}}*/
    177183/*FUNCTION BoolInput::ChangeEnum{{{1*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r5578 r5629  
    1111#include "./Input.h"
    1212#include "../../include/include.h"
     13class GaussTria;
    1314/*}}}*/
    1415
    … …  
    4748                void GetParameterValue(double* pvalue);
    4849                void GetParameterValue(double* pvalue,double* gauss);
     50                void GetParameterValue(double* pvalue,GaussTria* gauss);
    4951                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5052                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
     53                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5154                void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");};
    5255                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};

  • issm/trunk/src/c/objects/Inputs/DoubleInput.cpp

    r5103 r5629  
    182182void DoubleInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    183183/*}}}*/
     184/*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     185void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
     186/*}}}*/
    184187/*FUNCTION DoubleInput::GetParameterValues{{{1*/
    185188void DoubleInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    186189/*}}}*/
    187 /*FUNCTION DoubleInput::GetParameterDerivativeValue{{{1*/
     190/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    188191void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
     192/*}}}*/
     193/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
     194void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
    189195/*}}}*/
    190196/*FUNCTION DoubleInput::ChangeEnum{{{1*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r5578 r5629  
    1111#include "./Input.h"
    1212#include "../../include/include.h"
     13class GaussTria;
    1314/*}}}*/
    1415
    … …  
    4647                void GetParameterValue(double* pvalue);
    4748                void GetParameterValue(double* pvalue,double* gauss);
     49                void GetParameterValue(double* pvalue,GaussTria* gauss);
    4850                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    4951                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
     52                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5053                void GetParameterAverage(double* pvalue);
    5154                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};

  • issm/trunk/src/c/objects/Inputs/Input.h

    r5578 r5629  
    1212class Node;
    1313class ElementResult;
     14class GaussTria;
    1415/*}}}*/
    1516
    … …  
    2526                virtual void GetParameterValue(double* pvalue)=0;
    2627                virtual void GetParameterValue(double* pvalue,double* gauss)=0;
     28                virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0;
    2729                virtual void GetParameterValues(double* values,double* gauss_pointers, int numgauss)=0;
    2830                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0;
     31                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0;
    2932                virtual void GetParameterAverage(double* pvalue)=0;
    3033                virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss)=0;

  • issm/trunk/src/c/objects/Inputs/IntInput.cpp

    r5103 r5629  
    170170void IntInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    171171/*}}}*/
     172/*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     173void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
     174/*}}}*/
    172175/*FUNCTION IntInput::GetParameterValues{{{1*/
    173176void IntInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    174177/*}}}*/
    175 /*FUNCTION IntInput::GetParameterDerivativeValue{{{1*/
     178/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    176179void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
     180/*}}}*/
     181/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
     182void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
    177183/*}}}*/
    178184/*FUNCTION IntInput::ChangeEnum{{{1*/

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r5578 r5629  
    1111#include "./Input.h"
    1212#include "../../include/include.h"
     13class GaussTria;
    1314/*}}}*/
    1415
    … …  
    4748                void GetParameterValue(double* pvalue);
    4849                void GetParameterValue(double* pvalue,double* gauss);
     50                void GetParameterValue(double* pvalue,GaussTria* gauss);
    4951                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5052                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
     53                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5154                void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");};
    5255                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r5529 r5629  
    176176
    177177/*Object functions*/
    178 /*FUNCTION PentaVertexInput::GetParameterValue{{{1*/
     178/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    179179void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
     180
     181        /*Call PentaRef function*/
     182        PentaRef::GetParameterValue(pvalue,&values[0],gauss);
     183
     184}
     185/*}}}*/
     186/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     187void PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){
    180188
    181189        /*Call PentaRef function*/
    … …  
    201209}
    202210/*}}}*/
    203 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue{{{1*/
     211/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
    204212void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
     213
     214        /*Call PentaRef function*/
     215        PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
     216}
     217/*}}}*/
     218/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
     219void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
    205220
    206221        /*Call PentaRef function*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r5578 r5629  
    1111#include "./Input.h"
    1212#include "../Elements/PentaRef.h"
     13class GaussTria;
    1314/*}}}*/
    1415
    … …  
    4748                void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");};
    4849                void GetParameterValue(double* pvalue,double* gauss);
     50                void GetParameterValue(double* pvalue,GaussTria* gauss);
    4951                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5052                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
     53                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5154                void GetParameterAverage(double* pvalue);
    5255                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r5578 r5629  
    165165
    166166/*Object functions*/
    167 /*FUNCTION TriaVertexInput::GetParameterValue{{{1*/
     167/*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    168168void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){
     169
     170        /*Call TriaRef function*/
     171        TriaRef::GetParameterValue(pvalue,&values[0],gauss);
     172
     173}
     174/*}}}*/
     175/*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     176void TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){
    169177
    170178        /*Call TriaRef function*/
    … …  
    191199}
    192200/*}}}*/
    193 /*FUNCTION TriaVertexInput::GetParameterDerivativeValue{{{1*/
     201/*FUNCTION TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
    194202void TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
     203
     204        /*Call TriaRef function*/
     205        TriaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
     206}
     207/*}}}*/
     208/*FUNCTION TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
     209void TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
    195210
    196211        /*Call TriaRef function*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r5578 r5629  
    1111#include "./Input.h"
    1212#include "../Elements/TriaRef.h"
     13class GaussTria;
    1314/*}}}*/
    1415
    … …  
    4748                void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");}
    4849                void GetParameterValue(double* pvalue,double* gauss);
     50                void GetParameterValue(double* pvalue,GaussTria* gauss);
    4951                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5052                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
     53                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5154                void GetParameterAverage(double* pvalue);
    5255                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss);
Note: See TracChangeset for help on using the changeset viewer.