Changeset 2956


Ignore:
Timestamp:
02/04/10 08:04:08 (15 years ago)
Author:
seroussi
Message:

changed Params to Reference and Basic to nothing

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

Legend:

Unmodified
Added
Removed

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

    r2722 r2956  
    416416        int     ig1,ig2;
    417417        double  gauss_weight1,gauss_weight2;
    418         double  gauss_l1l2l3l4[4];
     418        double  gauss_coord[4];
    419419        int     order_area_gauss;
    420420        int     num_vert_gauss;
    … …  
    571571
    572572
    573                                 gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    574                                 gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    575                                 gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    576                                 gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     573                                gauss_coord[0]=*(first_gauss_area_coord+ig1);
     574                                gauss_coord[1]=*(second_gauss_area_coord+ig1);
     575                                gauss_coord[2]=*(third_gauss_area_coord+ig1);
     576                                gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    577577
    578578
    579579                                /*Get strain rate from velocity: */
    580                                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
    581                                 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
     580                                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_coord);
     581                                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_coord);
    582582
    583583                                /*Get viscosity: */
    … …  
    586586
    587587                                /*Get B and Bprime matrices: */
    588                                 GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
    589                                 GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     588                                GetB(&B[0][0], &xyz_list[0][0], gauss_coord);
     589                                GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
    590590
    591591                                /* Get Jacobian determinant: */
    592                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     592                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    593593
    594594                                /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    … …  
    960960        int     ig1,ig2;
    961961        double  gauss_weight1,gauss_weight2;
    962         double  gauss_l1l2l3l4[4];
     962        double  gauss_coord[4];
    963963        int     order_area_gauss;
    964964        int     num_vert_gauss;
    … …  
    10351035                        gauss_weight=gauss_weight1*gauss_weight2;
    10361036
    1037                         gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    1038                         gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    1039                         gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    1040                         gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     1037                        gauss_coord[0]=*(first_gauss_area_coord+ig1);
     1038                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
     1039                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
     1040                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    10411041
    10421042                        /*Get B and Bprime matrices: */
    1043                         GetB_vert(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
    1044                         GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     1043                        GetB_vert(&B[0][0], &xyz_list[0][0], gauss_coord);
     1044                        GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
    10451045
    10461046                        /* Get Jacobian determinant: */
    1047                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     1047                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    10481048                        DL_scalar=gauss_weight*Jdet;
    10491049
    … …  
    15271527        double* area_gauss_weights      =  NULL;
    15281528        double* vert_gauss_weights      =  NULL;
    1529         double  gauss_l1l2l3l4[4];
     1529        double  gauss_coord[4];
    15301530        int     order_area_gauss;
    15311531        int     num_vert_gauss;
    … …  
    15391539
    15401540        /*nodal functions: */
    1541         double l1l2l3l4l5l6[6];
     1541        double l1l6[6];
    15421542
    15431543        /*element vector at the gaussian points: */
    … …  
    16091609                                gauss_weight=gauss_weight1*gauss_weight2;
    16101610
    1611                                 gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    1612                                 gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    1613                                 gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    1614                                 gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     1611                                gauss_coord[0]=*(first_gauss_area_coord+ig1);
     1612                                gauss_coord[1]=*(second_gauss_area_coord+ig1);
     1613                                gauss_coord[2]=*(third_gauss_area_coord+ig1);
     1614                                gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    16151615
    16161616                                /*Compute thickness at gaussian point: */
    1617                                 GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
     1617                                GetParameterValue(&thickness, &h[0],gauss_coord);
    16181618
    16191619                                /*Compute slope at gaussian point: */
    1620                                 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3l4);
     1620                                GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_coord);
    16211621
    16221622                                /* Get Jacobian determinant: */
    1623                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     1623                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    16241624
    16251625                                /*Get nodal functions: */
    1626                                 GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
     1626                                GetNodalFunctions(l1l6, gauss_coord);
    16271627
    16281628                                /*Compute driving stress: */
    … …  
    16321632                                for (i=0;i<numgrids;i++){
    16331633                                        for (j=0;j<NDOF2;j++){
    1634                                                 pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3l4l5l6[i];
     1634                                                pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l6[i];
    16351635                                        }
    16361636                                }
    … …  
    19211921        double* area_gauss_weights      =  NULL;
    19221922        double* vert_gauss_weights      =  NULL;
    1923         double  gauss_l1l2l3l4[4];
     1923        double  gauss_coord[4];
    19241924        int     order_area_gauss;
    19251925        int     num_vert_gauss;
    … …  
    19351935        double  pe_g[numdof];
    19361936        double  pe_g_gaussian[numdof];
    1937         double l1l2l3l4l5l6[6];
     1937        double l1l6[6];
    19381938
    19391939        /*Spawning: */
    … …  
    20002000                        gauss_weight=gauss_weight1*gauss_weight2;
    20012001
    2002                         gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    2003                         gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    2004                         gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    2005                         gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     2002                        gauss_coord[0]=*(first_gauss_area_coord+ig1);
     2003                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
     2004                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
     2005                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    20062006
    20072007                        /*Get velocity derivative, with respect to x and y: */
    2008                         GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
    2009                         GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
     2008                        GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord);
     2009                        GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord);
    20102010                        dudx=du[0];
    20112011                        dvdy=dv[1];
    … …  
    20132013
    20142014                        /* Get Jacobian determinant: */
    2015                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     2015                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    20162016#ifdef _ISSM_DEBUG_
    20172017                        printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
    … …  
    20192019
    20202020                        /*Get nodal functions: */
    2021                         GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
     2021                        GetNodalFunctions(l1l6, gauss_coord);
    20222022
    20232023                        /*Build pe_g_gaussian vector: */
    20242024                        for (i=0;i<numgrids;i++){
    2025                                 pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l2l3l4l5l6[i];
     2025                                pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l6[i];
    20262026                        }
    20272027
    … …  
    25472547#undef __FUNCT__
    25482548#define __FUNCT__ "Penta::GetB"
    2549 void Penta::GetB(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     2549void Penta::GetB(double* B, double* xyz_list, double* gauss_coord){
    25502550
    25512551        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    2552          * For grid i, Bi can be expressed in the basic coordinate system
     2552         * For grid i, Bi can be expressed in the actual coordinate system
    25532553         * by:
    2554          *       Bi_basic=[ dh/dx          0      ]
     2554         *       Bi=[ dh/dx          0      ]
    25552555         *                [   0           dh/dy   ]
    25562556         *                [ 1/2*dh/dy  1/2*dh/dx  ]
    … …  
    25672567        const int NDOF2=2;
    25682568
    2569         double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
    2570 
    2571         /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    2572         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     2569        double dh1dh6[NDOF3][numgrids];
     2570
     2571        /*Get dh1dh6 in actual coordinate system: */
     2572        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
    25732573
    25742574#ifdef _ISSM_DEBUG_
    25752575        for (i=0;i<numgrids;i++){
    2576                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
     2576                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
    25772577        }
    25782578#endif
    … …  
    25802580        /*Build B: */
    25812581        for (i=0;i<numgrids;i++){
    2582                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i];
     2582                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
    25832583                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
    25842584
    25852585                *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
    2586                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
    2587 
    2588                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[1][i];
    2589                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[0][i];
    2590 
    2591                 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
     2586                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
     2587
     2588                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
     2589                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
     2590
     2591                *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
    25922592                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    25932593
    25942594                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    2595                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
     2595                *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
    25962596        }
    25972597
    … …  
    26042604
    26052605        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    2606          * For grid i, Bi' can be expressed in the basic coordinate system
     2606         * For grid i, Bi' can be expressed in the actual coordinate system
    26072607         * by:
    2608          *       Bi_artdiff_basic=[ dh/dx ]
     2608         *       Bi_artdiff=[ dh/dx ]
    26092609         *                       [ dh/dy ]
    26102610         * where h is the interpolation function for grid i.
    … …  
    26182618        int DOFPERGRID=1;
    26192619
    2620         /*Same thing in the basic coordinate system: */
    2621         double dh1dh6_basic[calculationdof][numgrids];
    2622 
    2623         /*Get dh1dh2dh3 in basic coordinates system : */
    2624         GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     2620        /*Same thing in the actual coordinate system: */
     2621        double dh1dh6[calculationdof][numgrids];
     2622
     2623        /*Get dh1dh2dh3 in actual coordinates system : */
     2624        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
    26252625
    26262626        /*Build B': */
    26272627        for (i=0;i<numgrids;i++){
    2628                 *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i];
    2629                 *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i];
     2628                *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     2629                *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    26302630        }
    26312631}
    … …  
    26372637
    26382638        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    2639          * For grid i, Bi' can be expressed in the basic coordinate system
     2639         * For grid i, Bi' can be expressed in the actual coordinate system
    26402640         * by:
    2641          *       Bi_advec_basic =[ h ]
     2641         *       Bi_advec =[ h ]
    26422642         *                       [ h ]
    26432643         *                       [ h ]
    … …  
    26522652        int DOFPERGRID=1;
    26532653
    2654         /*Same thing in the basic coordinate system: */
     2654        /*Same thing in the actual coordinate system: */
    26552655        double l1l6[6];
    26562656
    2657         /*Get dh1dh2dh3 in basic coordinates system : */
     2657        /*Get dh1dh2dh3 in actual coordinates system : */
    26582658        GetNodalFunctions(l1l6, gauss_coord);
    26592659
    … …  
    26722672
    26732673        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    2674          * For grid i, Bi' can be expressed in the basic coordinate system
     2674         * For grid i, Bi' can be expressed in the actual coordinate system
    26752675         * by:
    2676          *       Bi_conduct_basic=[ dh/dx ]
     2676         *       Bi_conduct=[ dh/dx ]
    26772677         *                       [ dh/dy ]
    26782678         *                       [ dh/dz ]
    … …  
    26872687        int DOFPERGRID=1;
    26882688
    2689         /*Same thing in the basic coordinate system: */
    2690         double dh1dh6_basic[calculationdof][numgrids];
    2691 
    2692         /*Get dh1dh2dh3 in basic coordinates system : */
    2693         GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     2689        /*Same thing in the actual coordinate system: */
     2690        double dh1dh6[calculationdof][numgrids];
     2691
     2692        /*Get dh1dh2dh3 in actual coordinates system : */
     2693        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
    26942694
    26952695        /*Build B': */
    26962696        for (i=0;i<numgrids;i++){
    2697                 *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i];
    2698                 *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i];
    2699                 *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i];
     2697                *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     2698                *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     2699                *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
    27002700        }
    27012701}
    … …  
    27042704#undef __FUNCT__
    27052705#define __FUNCT__ "Penta:GetB_vert"
    2706 void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     2706void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_coord){
    27072707
    27082708
    … …  
    27132713        const int NDOF3=3;
    27142714        const int numgrids=6;
    2715         double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
    2716 
    2717         /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    2718         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     2715        double dh1dh6[NDOF3][numgrids];
     2716
     2717        /*Get dh1dh6 in actual coordinate system: */
     2718        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
    27192719
    27202720#ifdef _ISSM_DEBUG_
    27212721        for (i=0;i<numgrids;i++){
    2722                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
     2722                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
    27232723        }
    27242724#endif
    … …  
    27262726        /*Build B: */
    27272727        for (i=0;i<numgrids;i++){
    2728                 B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
     2728                B[i]=dh1dh6[2][i]; 
    27292729        }
    27302730
    … …  
    27422742#undef __FUNCT__
    27432743#define __FUNCT__ "Penta::GetBPrime"
    2744 void Penta::GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     2744void Penta::GetBPrime(double* B, double* xyz_list, double* gauss_coord){
    27452745
    27462746        /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    2747          * For grid i, Bi can be expressed in the basic coordinate system
     2747         * For grid i, Bi can be expressed in the actual coordinate system
    27482748         * by:
    2749          *       Bi_basic=[ 2*dh/dx     dh/dy   ]
     2749         *       Bi=[ 2*dh/dx     dh/dy   ]
    27502750         *                [   dh/dx    2*dh/dy  ]
    27512751         *                [ dh/dy      dh/dx    ]
    … …  
    27622762        const int numgrids=6;
    27632763
    2764         double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
    2765 
    2766         /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    2767         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     2764        double dh1dh6[NDOF3][numgrids];
     2765
     2766        /*Get dh1dh6 in actual coordinate system: */
     2767        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
    27682768
    27692769#ifdef _ISSM_DEBUG_
    27702770        for (i=0;i<numgrids;i++){
    2771                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
     2771                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
    27722772        }
    27732773#endif
    … …  
    27752775        /*Build BPrime: */
    27762776        for (i=0;i<numgrids;i++){
    2777                 *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh2dh3dh4dh5dh6_basic[0][i];
    2778                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
    2779 
    2780                 *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i];
    2781                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh2dh3dh4dh5dh6_basic[1][i];
    2782 
    2783                 *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[1][i];
    2784                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[0][i];
    2785 
    2786                 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i];
     2777                *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];
     2778                *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
     2779
     2780                *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
     2781                *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
     2782
     2783                *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i];
     2784                *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i];
     2785
     2786                *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i];
    27872787                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    27882788
    27892789                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    2790                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i];
     2790                *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i];
    27912791        }
    27922792}
    … …  
    27992799
    28002800        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    2801          * For grid i, Bi' can be expressed in the basic coordinate system
     2801         * For grid i, Bi' can be expressed in the actual coordinate system
    28022802         * by:
    28032803         *       Biprime_advec=[ dh/dx ]
    … …  
    28142814        int DOFPERGRID=1;
    28152815
    2816         /*Same thing in the basic coordinate system: */
    2817         double dh1dh6_basic[calculationdof][numgrids];
    2818 
    2819         /*Get dh1dh2dh3 in basic coordinates system : */
    2820         GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     2816        /*Same thing in the actual coordinate system: */
     2817        double dh1dh6[calculationdof][numgrids];
     2818
     2819        /*Get dh1dh2dh3 in actual coordinates system : */
     2820        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
    28212821
    28222822        /*Build B': */
    28232823        for (i=0;i<numgrids;i++){
    2824                 *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6_basic[0][i];
    2825                 *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6_basic[1][i];
    2826                 *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6_basic[2][i];
     2824                *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     2825                *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     2826                *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
    28272827        }
    28282828}
    … …  
    28312831#undef __FUNCT__
    28322832#define __FUNCT__ "Penta:GetBPrime_vert"
    2833 void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     2833void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_coord){
    28342834
    28352835        // Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i
    … …  
    28372837        int i;
    28382838
    2839         GetNodalFunctions(B, gauss_l1l2l3l4);
     2839        GetNodalFunctions(B, gauss_coord);
    28402840
    28412841}
    … …  
    28472847
    28482848        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    2849          *      For grid i, Bi' can be expressed in the basic coordinate system
     2849         *      For grid i, Bi' can be expressed in the actual coordinate system
    28502850         *      by:
    2851          *                              Bi_basic'=[ dh/dx   0          0       0]
     2851         *                              Bi'=[ dh/dx   0          0       0]
    28522852         *                                       [   0      dh/dy      0       0]
    28532853         *                                       [   0      0         dh/dz    0]
    … …  
    28672867        int DOFPERGRID=4;
    28682868
    2869         double dh1dh7_basic[calculationdof][numgrids+1];
     2869        double dh1dh7[calculationdof][numgrids+1];
    28702870        double l1l6[numgrids];
    28712871
    2872         /*Get dh1dh7 in basic coordinate system: */
    2873         GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
     2872        /*Get dh1dh7 in actual coordinate system: */
     2873        GetNodalFunctionsDerivativesStokes(&dh1dh7[0][0],xyz_list, gauss_coord);
    28742874
    28752875        GetNodalFunctions(l1l6, gauss_coord);
    … …  
    28772877#ifdef _ISSM_DEBUG_
    28782878        for (i=0;i<6;i++){
    2879                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
     2879                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i]);
    28802880        }
    28812881
    … …  
    28842884        /*B_primeuild B_prime: */
    28852885        for (i=0;i<numgrids+1;i++){
    2886                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
     2886                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
    28872887                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    28882888                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    28892889                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    2890                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2890                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    28912891                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    28922892                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    28932893                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    2894                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
    2895                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7_basic[1][i];
    2896                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7_basic[0][i];
     2894                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
     2895                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];
     2896                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];
    28972897                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    2898                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7_basic[2][i];
     2898                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
    28992899                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    2900                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7_basic[0][i];
     2900                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
    29012901                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    2902                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7_basic[2][i];
    2903                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7_basic[1][i];
    2904                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7_basic[0][i];
    2905                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
    2906                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2902                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
     2903                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
     2904                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
     2905                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
     2906                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
    29072907                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
    29082908                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
    … …  
    29292929
    29302930        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
    2931          * For grid i, Bi can be expressed in the basic coordinate system
    2932          * by:                          Bi_basic=[ dh/dx          0             0       0  ]
     2931         * For grid i, Bi can be expressed in the actual coordinate system
     2932         * by:                          Bi=[ dh/dx          0             0       0  ]
    29332933         *                                      [   0           dh/dy           0       0  ]
    29342934         *                                      [   0             0           dh/dy     0  ]
    … …  
    29472947        int DOFPERGRID=4;
    29482948
    2949         double dh1dh7_basic[calculationdof][numgrids+1];
     2949        double dh1dh7[calculationdof][numgrids+1];
    29502950        double l1l6[numgrids];
    29512951
    29522952
    2953         /*Get dh1dh7 in basic coordinate system: */
    2954         GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
     2953        /*Get dh1dh7 in actual coordinate system: */
     2954        GetNodalFunctionsDerivativesStokes(&dh1dh7[0][0],xyz_list, gauss_coord);
    29552955
    29562956        GetNodalFunctions(l1l6, gauss_coord);
    … …  
    29582958#ifdef _ISSM_DEBUG_
    29592959        for (i=0;i<7;i++){
    2960                 printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
     2960                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i],dh1dh7[2][i]);
    29612961        }
    29622962
    … …  
    29652965        /*Build B: */
    29662966        for (i=0;i<numgrids+1;i++){
    2967                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7_basic[0][i]; //B[0][DOFPERGRID*i]=dh1dh6_basic[0][i];
     2967                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
    29682968                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    29692969                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    29702970                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    2971                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
     2971                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    29722972                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    29732973                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    29742974                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    2975                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
    2976                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7_basic[1][i];
    2977                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[0][i];
     2975                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
     2976                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
     2977                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
    29782978                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    2979                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7_basic[2][i];
     2979                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
    29802980                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    2981                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[0][i];
     2981                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
    29822982                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    2983                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7_basic[2][i];
    2984                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7_basic[1][i];
     2983                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
     2984                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
    29852985                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
    29862986                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
    29872987                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
    2988                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7_basic[0][i];
    2989                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7_basic[1][i];
    2990                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7_basic[2][i];
     2988                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
     2989                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
     2990                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
    29912991        }
    29922992
    … …  
    30413041#undef __FUNCT__
    30423042#define __FUNCT__ "Penta::GetJacobian"
    3043 void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4){
     3043void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_coord){
    30443044
    30453045        const int NDOF3=3;
    … …  
    30593059
    30603060        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    3061         A1=gauss_l1l2l3l4[0];
    3062         A2=gauss_l1l2l3l4[1];
    3063         A3=gauss_l1l2l3l4[2];
     3061        A1=gauss_coord[0];
     3062        A2=gauss_coord[1];
     3063        A3=gauss_coord[2];
    30643064
    30653065        xi=A2-A1;
    30663066        eta=sqrt3*A3;
    3067         zi=gauss_l1l2l3l4[3];
     3067        zi=gauss_coord[3];
    30683068
    30693069        x1=*(xyz_list+3*0+0);
    … …  
    31143114#undef __FUNCT__
    31153115#define __FUNCT__ "Penta::GetJacobianDeterminant"
    3116 void Penta::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4){
     3116void Penta::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_coord){
    31173117
    31183118        /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
    … …  
    31223122        double J[NDOF3][NDOF3];
    31233123
    3124         GetJacobian(&J[0][0],xyz_list,gauss_l1l2l3l4);
     3124        GetJacobian(&J[0][0],xyz_list,gauss_coord);
    31253125
    31263126        *Jdet= J[0][0]*J[1][1]*J[2][2]-J[0][0]*J[1][2]*J[2][1]-J[1][0]*J[0][1]*J[2][2]+J[1][0]*J[0][2]*J[2][1]+J[2][0]*J[0][1]*J[1][2]-J[2][0]*J[0][2]*J[1][1];
    … …  
    31373137        /*
    31383138         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    3139          * For grid i, Li can be expressed in the basic coordinate system
     3139         * For grid i, Li can be expressed in the actual coordinate system
    31403140         * by:
    3141          *       Li_basic=[ h    0    0   0]
     3141         *       Li=[ h    0    0   0]
    31423142         *               [ 0    h    0   0]
    31433143         *               [ 0    0    h   0]
    … …  
    31633163
    31643164
    3165         /*Get l1l2l3 in basic coordinate system: */
     3165        /*Get l1l2l3 in actual coordinate system: */
    31663166        l1l2l3[0]=gauss_coord_tria[0];
    31673167        l1l2l3[1]=gauss_coord_tria[1];
    … …  
    31773177        /*Build LStokes: */
    31783178        for (i=0;i<3;i++){
    3179                 *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     3179                *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
    31803180                *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    31813181                *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    … …  
    32443244        /*
    32453245         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    3246          * For grid i, Lpi can be expressed in the basic coordinate system
     3246         * For grid i, Lpi can be expressed in the actual coordinate system
    32473247         * by:
    3248          *       Lpi_basic=[ h    0    0   0]
     3248         *       Lpi=[ h    0    0   0]
    32493249         *               [ 0    h    0   0]
    32503250         *               [ h    0    0   0]
    … …  
    32683268
    32693269        double l1l2l3[numgrids2d];
    3270         double dh1dh6_basic[3][6];
    3271 
    3272 
    3273         /*Get l1l2l3 in basic coordinate system: */
     3270        double dh1dh6[3][6];
     3271
     3272
     3273        /*Get l1l2l3 in actual coordinate system: */
    32743274        l1l2l3[0]=gauss_coord_tria[0];
    32753275        l1l2l3[1]=gauss_coord_tria[1];
    32763276        l1l2l3[2]=gauss_coord_tria[2];
    32773277
    3278         GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
     3278        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list,gauss_coord);
    32793279
    32803280#ifdef _DELUG_
    … …  
    32863286        /*Build LprimeStokes: */
    32873287        for (i=0;i<3;i++){
    3288                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     3288                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    32893289                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    32903290                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    … …  
    33123312                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
    33133313                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
    3314                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6_basic[2][i];
     3314                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
    33153315                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
    33163316                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
    33173317                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
    3318                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6_basic[2][i];
     3318                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
    33193319                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
    33203320                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
    33213321                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
    3322                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6_basic[2][i];
     3322                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
    33233323                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
    3324                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6_basic[2][i];
     3324                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
    33253325                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
    3326                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6_basic[0][i];
     3326                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
    33273327                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
    33283328                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
    3329                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6_basic[2][i];
    3330                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6_basic[1][i];
     3329                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
     3330                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
    33313331                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
    33323332                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
    … …  
    33863386}
    33873387/*}}}*/
    3388 /*FUNCTION GetNodalFunctionsDerivativesBasic {{{1*/
     3388/*FUNCTION GetNodalFunctionsDerivatives {{{1*/
    33893389#undef __FUNCT__
    3390 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic"
    3391 void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
    3392 
    3393         /*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
     3390#define __FUNCT__ "Penta::GetNodalFunctionsDerivatives"
     3391void Penta::GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){
     3392
     3393        /*This routine returns the values of the nodal functions derivatives  (with respect to the actual coordinate system: */
    33943394
    33953395
    … …  
    33983398        const int numgrids=6;
    33993399
    3400         double dh1dh2dh3dh4dh5dh6_param[NDOF3][numgrids];
     3400        double dh1dh6_ref[NDOF3][numgrids];
    34013401        double Jinv[NDOF3][NDOF3];
    34023402
    34033403        /*Get derivative values with respect to parametric coordinate system: */
    3404         GetNodalFunctionsDerivativesParams(&dh1dh2dh3dh4dh5dh6_param[0][0], gauss_l1l2l3l4);
     3404        GetNodalFunctionsDerivativesReference(&dh1dh6_ref[0][0], gauss_coord);
    34053405
    34063406        /*Get Jacobian invert: */
    3407         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3l4);
    3408 
    3409         /*Build dh1dh2dh3_basic:
     3407        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
     3408
     3409        /*Build dh1dh3:
    34103410         *
    34113411         * [dhi/dx]= Jinv*[dhi/dr]
    … …  
    34153415
    34163416        for (i=0;i<numgrids;i++){
    3417                 *(dh1dh2dh3dh4dh5dh6_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[0][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[0][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
    3418                 *(dh1dh2dh3dh4dh5dh6_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[1][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[1][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
    3419                 *(dh1dh2dh3dh4dh5dh6_basic+numgrids*2+i)=Jinv[2][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[2][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[2][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
    3420         }
    3421 
    3422 }
    3423 /*}}}*/
    3424 /*FUNCTION GetNodalFunctionsDerivativesBasicStokes {{{1*/
    3425 #undef __FUNCT__
    3426 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
    3427 void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
     3417                *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
     3418                *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
     3419                *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
     3420        }
     3421
     3422}
     3423/*}}}*/
     3424/*FUNCTION GetNodalFunctionsDerivativesStokes {{{1*/
     3425#undef __FUNCT__
     3426#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesStokes"
     3427void Penta::GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord){
    34283428
    34293429        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    3430          * basic coordinate system: */
     3430         * actual coordinate system: */
    34313431
    34323432        int i;
    34333433
    34343434        const  int numgrids=7;
    3435         double dh1dh7_param[3][numgrids];
     3435        double dh1dh7_ref[3][numgrids];
    34363436        double Jinv[3][3];
    34373437
    34383438
    34393439        /*Get derivative values with respect to parametric coordinate system: */
    3440         GetNodalFunctionsDerivativesParamsStokes(&dh1dh7_param[0][0], gauss_coord);
     3440        GetNodalFunctionsDerivativesReferenceStokes(&dh1dh7_ref[0][0], gauss_coord);
    34413441
    34423442        /*Get Jacobian invert: */
    34433443        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_coord);
    34443444
    3445         /*Build dh1dh6_basic:
     3445        /*Build dh1dh6:
    34463446         *
    34473447         * [dhi/dx]= Jinv'*[dhi/dr]
    … …  
    34513451
    34523452        for (i=0;i<numgrids;i++){
    3453                 *(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[0][1]*dh1dh7_param[1][i]+Jinv[0][2]*dh1dh7_param[2][i];
    3454                 *(dh1dh7_basic+numgrids*1+i)=Jinv[1][0]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[1][2]*dh1dh7_param[2][i];
    3455                 *(dh1dh7_basic+numgrids*2+i)=Jinv[2][0]*dh1dh7_param[0][i]+Jinv[2][1]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
    3456         }
    3457 
    3458 }
    3459 /*}}}*/
    3460 /*FUNCTION GetNodalFunctionsDerivativesParams {{{1*/
     3453                *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
     3454                *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
     3455                *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
     3456        }
     3457
     3458}
     3459/*}}}*/
     3460/*FUNCTION GetNodalFunctionsDerivativesReference {{{1*/
    34613461#undef __FUNCT__
    3462 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams"
    3463 void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
     3462#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesReference"
     3463void Penta::GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord){
    34643464
    34653465        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    34703470        double sqrt3=sqrt(3.0);
    34713471
    3472         A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
    3473         A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
    3474         A3=gauss_l1l2l3l4[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
    3475         z=gauss_l1l2l3l4[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
     3472        A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
     3473        A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
     3474        A3=gauss_coord[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
     3475        z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
    34763476
    34773477
    34783478        /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
    3479         *(dl1dl2dl3dl4dl5dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
    3480         *(dl1dl2dl3dl4dl5dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
    3481         *(dl1dl2dl3dl4dl5dl6+numgrids*2+0)=-1.0/2.0*A1;
     3479        *(dl1dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
     3480        *(dl1dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
     3481        *(dl1dl6+numgrids*2+0)=-1.0/2.0*A1;
    34823482
    34833483        /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
    3484         *(dl1dl2dl3dl4dl5dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
    3485         *(dl1dl2dl3dl4dl5dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
    3486         *(dl1dl2dl3dl4dl5dl6+numgrids*2+1)=-1.0/2.0*A2;
     3484        *(dl1dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
     3485        *(dl1dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
     3486        *(dl1dl6+numgrids*2+1)=-1.0/2.0*A2;
    34873487
    34883488        /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
    3489         *(dl1dl2dl3dl4dl5dl6+numgrids*0+2)=0.0;
    3490         *(dl1dl2dl3dl4dl5dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
    3491         *(dl1dl2dl3dl4dl5dl6+numgrids*2+2)=-1.0/2.0*A3;
     3489        *(dl1dl6+numgrids*0+2)=0.0;
     3490        *(dl1dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
     3491        *(dl1dl6+numgrids*2+2)=-1.0/2.0*A3;
    34923492
    34933493        /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
    3494         *(dl1dl2dl3dl4dl5dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
    3495         *(dl1dl2dl3dl4dl5dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
    3496         *(dl1dl2dl3dl4dl5dl6+numgrids*2+3)=1.0/2.0*A1;
     3494        *(dl1dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
     3495        *(dl1dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
     3496        *(dl1dl6+numgrids*2+3)=1.0/2.0*A1;
    34973497
    34983498        /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
    3499         *(dl1dl2dl3dl4dl5dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
    3500         *(dl1dl2dl3dl4dl5dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
    3501         *(dl1dl2dl3dl4dl5dl6+numgrids*2+4)=1.0/2.0*A2;
     3499        *(dl1dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
     3500        *(dl1dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
     3501        *(dl1dl6+numgrids*2+4)=1.0/2.0*A2;
    35023502
    35033503        /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
    3504         *(dl1dl2dl3dl4dl5dl6+numgrids*0+5)=0.0;
    3505         *(dl1dl2dl3dl4dl5dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
    3506         *(dl1dl2dl3dl4dl5dl6+numgrids*2+5)=1.0/2.0*A3;
    3507 }
    3508 /*}}}*/
    3509 /*FUNCTION GetNodalFunctionsDerivativesParamsStokes {{{1*/
    3510 #undef __FUNCT__
    3511 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
    3512 void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
     3504        *(dl1dl6+numgrids*0+5)=0.0;
     3505        *(dl1dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
     3506        *(dl1dl6+numgrids*2+5)=1.0/2.0*A3;
     3507}
     3508/*}}}*/
     3509/*FUNCTION GetNodalFunctionsDerivativesReferenceStokes {{{1*/
     3510#undef __FUNCT__
     3511#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesReferenceStokes"
     3512void Penta::GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord){
    35133513
    35143514        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    35923592#undef __FUNCT__
    35933593#define __FUNCT__ "Penta::GetParameterDerivativeValue"
    3594 void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
    3595 
    3596         /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4:
     3594void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord){
     3595
     3596        /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
    35973597         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    35983598         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
    … …  
    36043604        const int NDOF3=3;
    36053605        const int numgrids=6;
    3606         double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
    3607 
    3608         /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    3609         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    3610 
    3611         *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
     3606        double dh1dh6[NDOF3][numgrids];
     3607
     3608        /*Get dh1dh6 in actual coordinate system: */
     3609        GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
     3610
     3611        *(p+0)=p_list[0]*dh1dh6[0][0]+p_list[1]*dh1dh6[0][1]+p_list[2]*dh1dh6[0][2]+p_list[3]*dh1dh6[0][3]+p_list[4]*dh1dh6[0][4]+p_list[5]*dh1dh6[0][5];
    36123612        ;
    3613         *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
    3614 
    3615         *(p+2)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[2][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[2][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[2][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[2][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[2][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[2][5];
     3613        *(p+1)=p_list[0]*dh1dh6[1][0]+p_list[1]*dh1dh6[1][1]+p_list[2]*dh1dh6[1][2]+p_list[3]*dh1dh6[1][3]+p_list[4]*dh1dh6[1][4]+p_list[5]*dh1dh6[1][5];
     3614
     3615        *(p+2)=p_list[0]*dh1dh6[2][0]+p_list[1]*dh1dh6[2][1]+p_list[2]*dh1dh6[2][2]+p_list[3]*dh1dh6[2][3]+p_list[4]*dh1dh6[2][4]+p_list[5]*dh1dh6[2][5];
    36163616
    36173617}
    … …  
    36203620#undef __FUNCT__
    36213621#define __FUNCT__ "Penta::GetParameterValue"
    3622 void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
     3622void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_coord){
    36233623
    36243624        const int numgrids=6;
    3625         double l1l2l3l4l5l6[numgrids];
    3626 
    3627         GetNodalFunctions(&l1l2l3l4l5l6[0], gauss_l1l2l3l4);
    3628 
    3629         *pvalue=l1l2l3l4l5l6[0]*v_list[0]+l1l2l3l4l5l6[1]*v_list[1]+l1l2l3l4l5l6[2]*v_list[2]+l1l2l3l4l5l6[3]*v_list[3]+l1l2l3l4l5l6[4]*v_list[4]+l1l2l3l4l5l6[5]*v_list[5];
     3625        double l1l6[numgrids];
     3626
     3627        GetNodalFunctions(&l1l6[0], gauss_coord);
     3628
     3629        *pvalue=l1l6[0]*v_list[0]+l1l6[1]*v_list[1]+l1l6[2]*v_list[2]+l1l6[3]*v_list[3]+l1l6[4]*v_list[4]+l1l6[5]*v_list[5];
    36303630}
    36313631/*}}}*/
    … …  
    36333633#undef __FUNCT__
    36343634#define __FUNCT__ "Penta::GetJacobianInvert"
    3635 void Penta::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4){
     3635void Penta::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_coord){
    36363636
    36373637        double Jdet;
    … …  
    36393639
    36403640        /*Call Jacobian routine to get the jacobian:*/
    3641         GetJacobian(Jinv, xyz_list, gauss_l1l2l3l4);
     3641        GetJacobian(Jinv, xyz_list, gauss_coord);
    36423642
    36433643        /*Invert Jacobian matrix: */
    … …  
    36533653#undef __FUNCT__
    36543654#define __FUNCT__ "Penta::GetNodalFunctions"
    3655 void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
     3655void Penta::GetNodalFunctions(double* l1l6, double* gauss_coord){
    36563656
    36573657        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    36583658
    3659         l1l2l3l4l5l6[0]=gauss_l1l2l3l4[0]*(1-gauss_l1l2l3l4[3])/2.0;
    3660 
    3661         l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
    3662 
    3663         l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
    3664 
    3665         l1l2l3l4l5l6[3]=gauss_l1l2l3l4[0]*(1+gauss_l1l2l3l4[3])/2.0;
    3666 
    3667         l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
    3668 
    3669         l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
     3659        l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0;
     3660
     3661        l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0;
     3662
     3663        l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0;
     3664
     3665        l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0;
     3666
     3667        l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0;
     3668
     3669        l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0;
    36703670
    36713671}
    … …  
    37303730#undef __FUNCT__
    37313731#define __FUNCT__ "Penta::GetStrainRate"
    3732 void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4){
     3732void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord){
    37333733
    37343734        int i;
    … …  
    37393739
    37403740        /*Get B matrix: */
    3741         GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
     3741        GetB(&B[0][0], xyz_list, gauss_coord);
    37423742
    37433743#ifdef _ISSM_DEBUG_

  • issm/trunk/src/c/objects/Penta.h

    r2722 r2956  
    9595                void*  SpawnTria(int g0, int g1, int g2);
    9696
    97                 void  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4);
    98                 void  GetB(double* pB, double* xyz_list, double* gauss_l1l2l3l4);
    99                 void  GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4);
    100                 void  GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
    101                 void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
    102                 void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4);
    103                 void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4);
    104                 void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4);
    105                 void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4);
    106                 void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
     97                void  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
     98                void  GetB(double* pB, double* xyz_list, double* gauss_coord);
     99                void  GetBPrime(double* B, double* xyz_list, double* gauss_coord);
     100                void  GetB_vert(double* B, double* xyz_list, double* gauss_coord);
     101                void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_coord);
     102                void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_coord);
     103                void  GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord);
     104                void  GetJacobian(double* J, double* xyz_list,double* gauss_coord);
     105                void  GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord);
     106                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_coord);
    107107                void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
    108108                void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
    109                 void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
    110                 void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
    111                 void  GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4);
     109                void  GetParameterValue(double* pvalue, double* v_list,double* gauss_coord);
     110                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord);
     111                void  GetNodalFunctions(double* l1l6, double* gauss_coord);
    112112                void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
    113113                void  ComputePressure(Vec p_g);
    … …  
    131131                void  GetLStokes(double* LStokes, double* gauss_coord_tria);
    132132                void  GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_coord_tria, double* gauss_coord);
    133                 void  GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord);
    134                 void  GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord);
     133                void  GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord);
     134                void  GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord);
    135135                void  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);
    136136                void  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);

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

    r2722 r2956  
    30543054
    30553055        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    3056          * For grid i, Bi can be expressed in the basic coordinate system
     3056         * For grid i, Bi can be expressed in the actual coordinate system
    30573057         * by:
    3058          *       Bi_basic=[ dh/dx    0    ]
     3058         *       Bi=[ dh/dx    0    ]
    30593059         *                [   0    dh/dy  ]
    30603060         *                [ 1/2*dh/dy  1/2*dh/dx  ]
    … …  
    30683068        const int numgrids=3;
    30693069
    3070         double dh1dh2dh3_basic[NDOF2][numgrids];
    3071 
    3072 
    3073         /*Get dh1dh2dh3 in basic coordinate system: */
    3074         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
     3070        double dh1dh3[NDOF2][numgrids];
     3071
     3072
     3073        /*Get dh1dh2dh3 in actual coordinate system: */
     3074        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
    30753075
    30763076        #ifdef _ISSM_DEBUG_
    30773077        for (i=0;i<3;i++){
    3078                 printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]);
     3078                printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh3[0][i],dh1dh3[1][i]);
    30793079        }
    30803080        #endif
    … …  
    30823082        /*Build B: */
    30833083        for (i=0;i<numgrids;i++){
    3084                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3_basic[0][i]; //B[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     3084                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
    30853085                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
    30863086                *(B+NDOF2*numgrids*1+NDOF2*i)=0;
    3087                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3_basic[1][i];
    3088                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3_basic[1][i];
    3089                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3_basic[0][i];
     3087                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
     3088                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];
     3089                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
    30903090        }
    30913091}
    … …  
    30983098
    30993099        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    3100          * For grid i, Bi can be expressed in the basic coordinate system
     3100         * For grid i, Bi can be expressed in the actual coordinate system
    31013101         * by:
    3102          *       Bi_basic=[ h ]
     3102         *       Bi=[ h ]
    31033103         *                [ h ]
    31043104         * where h is the interpolation function for grid i.
    … …  
    31143114
    31153115
    3116         /*Get dh1dh2dh3 in basic coordinate system: */
     3116        /*Get dh1dh2dh3 in actual coordinate system: */
    31173117        GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
    31183118
    … …  
    31453145
    31463146        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    3147          * For grid i, Bi' can be expressed in the basic coordinate system
     3147         * For grid i, Bi' can be expressed in the actual coordinate system
    31483148         * by:
    3149          *       Bi_prime__basic=[ 2*dh/dx dh/dy ]
     3149         *       Bi_prime=[ 2*dh/dx dh/dy ]
    31503150         *                       [ dh/dx  2*dh/dy]
    31513151         *                       [dh/dy dh/dx]
    … …  
    31593159        const int numgrids=3;
    31603160
    3161         /*Same thing in the basic coordinate system: */
    3162         double dh1dh2dh3_basic[NDOF2][numgrids];
    3163 
    3164 
    3165         /*Get dh1dh2dh3 in basic coordinates system : */
    3166         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
     3161        /*Same thing in the actual coordinate system: */
     3162        double dh1dh3[NDOF2][numgrids];
     3163
     3164
     3165        /*Get dh1dh2dh3 in actual coordinates system : */
     3166        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
    31673167
    31683168        /*Build B': */
    31693169        for (i=0;i<numgrids;i++){
    3170                 *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh2dh3_basic[0][i];
    3171                 *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3_basic[1][i];
    3172                 *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3_basic[0][i];
    3173                 *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh2dh3_basic[1][i];
    3174                 *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3_basic[1][i];
    3175                 *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3_basic[0][i];
     3170                *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i];
     3171                *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i];
     3172                *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i];
     3173                *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];
     3174                *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i];
     3175                *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i];
    31763176        }
    31773177}
    … …  
    31843184
    31853185        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    3186          * For grid i, Bi' can be expressed in the basic coordinate system
     3186         * For grid i, Bi' can be expressed in the actual coordinate system
    31873187         * by:
    3188          *       Bi_prime__basic=[ dh/dx ]
     3188         *       Bi_prime=[ dh/dx ]
    31893189         *                       [ dh/dy ]
    31903190         * where h is the interpolation function for grid i.
    … …  
    31983198        const int numgrids=3;
    31993199
    3200         /*Same thing in the basic coordinate system: */
    3201         double dh1dh2dh3_basic[NDOF2][numgrids];
    3202 
    3203         /*Get dh1dh2dh3 in basic coordinates system : */
    3204         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
     3200        /*Same thing in the actual coordinate system: */
     3201        double dh1dh3[NDOF2][numgrids];
     3202
     3203        /*Get dh1dh2dh3 in actual coordinates system : */
     3204        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
    32053205
    32063206        /*Build B': */
    32073207        for (i=0;i<numgrids;i++){
    3208                 *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh2dh3_basic[0][i];
    3209                 *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh2dh3_basic[1][i];
     3208                *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i];
     3209                *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i];
    32103210        }
    32113211}
    … …  
    33513351
    33523352        /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    3353          * For grid i, Li can be expressed in the basic coordinate system
     3353         * For grid i, Li can be expressed in the actual coordinate system
    33543354         * by:
    33553355         *       numdof=1:
    3356          *       Li_basic=h;
     3356         *       Li=h;
    33573357         *       numdof=2:
    3358          *       Li_basic=[ h    0    ]
     3358         *       Li=[ h    0    ]
    33593359         *                [   0   h  ]
    33603360         * where h is the interpolation function for grid i.
    … …  
    33703370
    33713371
    3372         /*Get l1l2l3 in basic coordinate system: */
     3372        /*Get l1l2l3 in actual coordinate system: */
    33733373        GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    33743374
    … …  
    33873387        else{
    33883388                for (i=0;i<numgrids;i++){
    3389                         *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
     3389                        *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
    33903390                        *(L+numdof*numgrids*0+numdof*i+1)=0;
    33913391                        *(L+numdof*numgrids*1+numdof*i)=0;
    … …  
    34233423}
    34243424/*}}}*/
    3425 /*FUNCTION GetNodalFunctionsDerivativesBasic {{{1*/
    3426 #undef __FUNCT__
    3427 #define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
    3428 void Tria::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3){
     3425/*FUNCTION GetNodalFunctionsDerivatives {{{1*/
     3426#undef __FUNCT__
     3427#define __FUNCT__ "Tria::GetNodalFunctionsDerivatives"
     3428void Tria::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3){
    34293429       
    34303430        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    3431          * basic coordinate system: */
     3431         * actual coordinate system: */
    34323432
    34333433        int i;
    … …  
    34353435        const int numgrids=3;
    34363436
    3437         double dh1dh2dh3_param[NDOF2][numgrids];
     3437        double dh1dh3_ref[NDOF2][numgrids];
    34383438        double Jinv[NDOF2][NDOF2];
    34393439
    34403440
    34413441        /*Get derivative values with respect to parametric coordinate system: */
    3442         GetNodalFunctionsDerivativesParams(&dh1dh2dh3_param[0][0], gauss_l1l2l3);
     3442        GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss_l1l2l3);
    34433443
    34443444        /*Get Jacobian invert: */
    34453445        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3);
    34463446
    3447         /*Build dh1dh2dh3_basic:
     3447        /*Build dh1dh3:
    34483448         *
    34493449         * [dhi/dx]= Jinv*[dhi/dr]
    … …  
    34523452
    34533453        for (i=0;i<numgrids;i++){
    3454                 *(dh1dh2dh3_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3_param[0][i]+Jinv[0][1]*dh1dh2dh3_param[1][i];
    3455                 *(dh1dh2dh3_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3_param[0][i]+Jinv[1][1]*dh1dh2dh3_param[1][i];
    3456         }
    3457 
    3458 }
    3459 /*}}}*/
    3460 /*FUNCTION GetNodalFunctionsDerivativesParams {{{1*/
    3461 #undef __FUNCT__
    3462 #define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
    3463 void Tria::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3){
     3454                *(dh1dh3+numgrids*0+i)=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
     3455                *(dh1dh3+numgrids*1+i)=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
     3456        }
     3457
     3458}
     3459/*}}}*/
     3460/*FUNCTION GetNodalFunctionsDerivativesReference {{{1*/
     3461#undef __FUNCT__
     3462#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesReference"
     3463void Tria::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
    34643464       
    34653465        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    34723472
    34733473        /*First nodal function: */
    3474         *(dl1dl2dl3+numgrids*0+0)=-1.0/2.0;
    3475         *(dl1dl2dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);
     3474        *(dl1dl3+numgrids*0+0)=-1.0/2.0;
     3475        *(dl1dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);
    34763476
    34773477        /*Second nodal function: */
    3478         *(dl1dl2dl3+numgrids*0+1)=1.0/2.0;
    3479         *(dl1dl2dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);
     3478        *(dl1dl3+numgrids*0+1)=1.0/2.0;
     3479        *(dl1dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);
    34803480
    34813481        /*Third nodal function: */
    3482         *(dl1dl2dl3+numgrids*0+2)=0;
    3483         *(dl1dl2dl3+numgrids*1+2)=1.0/sqrt3;
     3482        *(dl1dl3+numgrids*0+2)=0;
     3483        *(dl1dl3+numgrids*1+2)=1.0/sqrt3;
    34843484
    34853485}
    … …  
    35153515         */
    35163516       
    3517         double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system.
    3518 
    3519         /*Get dh1dh2dh3 in basic coordinate system: */
    3520         GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
    3521 
    3522         *(p+0)=plist[0]*dh1dh2dh3_basic[0][0]+plist[1]*dh1dh2dh3_basic[0][1]+plist[2]*dh1dh2dh3_basic[0][2];
    3523         *(p+1)=plist[0]*dh1dh2dh3_basic[1][0]+plist[1]*dh1dh2dh3_basic[1][1]+plist[2]*dh1dh2dh3_basic[1][2];
     3517        double dh1dh3[NDOF2][numgrids]; //nodal derivative functions in actual coordinate system.
     3518
     3519        /*Get dh1dh2dh3 in actual coordinate system: */
     3520        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
     3521
     3522        *(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
     3523        *(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
    35243524
    35253525}
    … …  
    36113611        double       xyz_list[numgrids][3];
    36123612        int          doflist1[numgrids];
    3613         double       dh1dh2dh3_basic[NDOF2][numgrids];
     3613        double       dh1dh3[NDOF2][numgrids];
    36143614
    36153615        /* grid data: */
    … …  
    37333733
    37343734                /*Get nodal functions derivatives*/
    3735                 GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     3735                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    37363736
    37373737                /*Get B derivative: dB/dx */
    … …  
    37453745
    37463746                        //Add regularization term
    3747                         grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]);
     3747                        grade_g_gaussian[i]-=numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);
    37483748
    37493749                        //min dampening
    … …  
    37873787        double       xyz_list[numgrids][3];
    37883788        int          doflist1[numgrids];
    3789         double       dh1dh2dh3_basic[NDOF2][numgrids];
     3789        double       dh1dh3[NDOF2][numgrids];
    37903790
    37913791        /* grid data: */
    … …  
    39533953
    39543954                /*Get nodal functions derivatives*/
    3955                 GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     3955                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    39563956
    39573957                /*Get k derivative: dk/dx */
    … …  
    39653965
    39663966                        //noise dampening d/dki(1/2*(dk/dx)^2)
    3967                         grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
     3967                        grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    39683968                       
    39693969                        //min dampening
    … …  
    40054005        double       xyz_list[numgrids][3];
    40064006        int          doflist1[numgrids];
    4007         double       dh1dh2dh3_basic[NDOF2][numgrids];
     4007        double       dh1dh3[NDOF2][numgrids];
    40084008
    40094009        /* grid data: */
    … …  
    41864186
    41874187                /*Get nodal functions derivatives*/
    4188                 GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3);
     4188                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    41894189
    41904190                /*Get k derivative: dk/dx */
    … …  
    42014201
    42024202                        //Add regularization term
    4203                         grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dk[0]+dh1dh2dh3_basic[1][i]*dk[1]);
     4203                        grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    42044204
    42054205                        //min dampening

  • issm/trunk/src/c/objects/Tria.h

    r2722 r2956  
    8989                void  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
    9090                void  GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);
    91                 void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3);
    92                 void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3);
     91                void  GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3);
     92                void  GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);
    9393                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
    9494                void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
Note: See TracChangeset for help on using the changeset viewer.