Changeset 5647


Ignore:
Timestamp:
09/01/10 10:45:00 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added GaussPenta for Pattyn

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

Legend:

Unmodified
Added
Removed

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

    r5640 r5647  
    27132713
    27142714        /* 3d gaussian points: */
    2715         int     num_gauss,ig;
    2716         double* first_gauss_area_coord  =  NULL;
    2717         double* second_gauss_area_coord =  NULL;
    2718         double* third_gauss_area_coord  =  NULL;
    2719         double* fourth_gauss_vert_coord  =  NULL;
    2720         double* area_gauss_weights           =  NULL;
    2721         double* vert_gauss_weights           =  NULL;
    2722         int     ig1,ig2;
    2723         double  gauss_weight1,gauss_weight2;
    2724         double  gauss_coord[4];
    2725         int     order_area_gauss;
    2726         int     num_vert_gauss;
    2727         int     num_area_gauss;
    2728         double  gauss_weight;
    2729 
    2730         /* 2d gaussian point: */
    2731         int     num_gauss2d;
    2732         double* first_gauss_area_coord2d  =  NULL;
    2733         double* second_gauss_area_coord2d =  NULL;
    2734         double* third_gauss_area_coord2d  =  NULL;
    2735         double* gauss_weights2d=NULL;
    2736         double  gauss_l1l2l3[3];
     2715        int     ig;
     2716        GaussPenta *gauss=NULL;
    27372717
    27382718        /* material data: */
    … …  
    28022782        GetDofList(&doflist,PattynApproximationEnum);
    28032783
    2804         /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    2805           get tria gaussian points as well as segment gaussian points. For tria gaussian
    2806           points, order of integration is 2, because we need to integrate the product tB*D*B'
    2807           which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    2808           points, same deal, which yields 3 gaussian points.*/
    2809 
    2810         order_area_gauss=5;
    2811         num_vert_gauss=5;
    2812 
    2813         gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    2814 
    28152784        /*Retrieve all inputs we will be needing: */
    28162785        vx_input=inputs->GetInput(VxEnum);
    … …  
    28202789
    28212790        /* Start  looping on the number of gaussian points: */
    2822         for (ig1=0; ig1<num_area_gauss; ig1++){
    2823                 for (ig2=0; ig2<num_vert_gauss; ig2++){
    2824 
    2825                         /*Pick up the gaussian point: */
    2826                         gauss_weight1=*(area_gauss_weights+ig1);
    2827                         gauss_weight2=*(vert_gauss_weights+ig2);
    2828                         gauss_weight=gauss_weight1*gauss_weight2;
    2829 
    2830 
    2831                         gauss_coord[0]=*(first_gauss_area_coord+ig1);
    2832                         gauss_coord[1]=*(second_gauss_area_coord+ig1);
    2833                         gauss_coord[2]=*(third_gauss_area_coord+ig1);
    2834                         gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    2835 
    2836 
    2837                         /*Get strain rate from velocity: */
    2838                         this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
    2839                         this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
    2840 
    2841                         /*Get viscosity: */
    2842                         matice->GetViscosity3d(&viscosity, &epsilon[0]);
    2843                         matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    2844 
    2845                         /*Get B and Bprime matrices: */
    2846                         GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord);
    2847                         GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
    2848 
    2849                         /* Get Jacobian determinant: */
    2850                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    2851 
    2852                         /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
    2853                           onto this scalar matrix, so that we win some computational time: */
    2854 
    2855                         newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    2856                         D_scalar=2*newviscosity*gauss_weight*Jdet;
    2857                         for (i=0;i<5;i++){
    2858                                 D[i][i]=D_scalar;
    2859                         }
    2860 
    2861                         /*  Do the triple product tB*D*Bprime: */
    2862                         TripleMultiply( &B[0][0],5,numdof,1,
    2863                                                 &D[0][0],5,5,0,
    2864                                                 &Bprime[0][0],5,numdof,0,
    2865                                                 &Ke_gg_gaussian[0][0],0);
    2866 
    2867                         /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    2868                         for( i=0; i<numdof; i++){
    2869                                 for(j=0;j<numdof;j++){
    2870                                         Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2871                                 }
    2872                         }
    2873                 } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    2874         } //for (ig1=0; ig1<num_area_gauss; ig1++)
     2791        gauss=new GaussPenta(5,5);
     2792        for (ig=gauss->begin();ig<gauss->end();ig++){
     2793
     2794                gauss->GaussPoint(ig);
     2795
     2796                /*Get strain rate from velocity: */
     2797                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     2798                this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     2799
     2800                /*Get viscosity: */
     2801                matice->GetViscosity3d(&viscosity, &epsilon[0]);
     2802                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
     2803
     2804                /*Get B and Bprime matrices: */
     2805                GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
     2806                GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
     2807
     2808                /* Get Jacobian determinant: */
     2809                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     2810
     2811                /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
     2812                  onto this scalar matrix, so that we win some computational time: */
     2813
     2814                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     2815                D_scalar=2*newviscosity*gauss->weight*Jdet;
     2816                for (i=0;i<5;i++) D[i][i]=D_scalar;
     2817
     2818                /*  Do the triple product tB*D*Bprime: */
     2819                TripleMultiply( &B[0][0],5,numdof,1,
     2820                                        &D[0][0],5,5,0,
     2821                                        &Bprime[0][0],5,numdof,0,
     2822                                        &Ke_gg_gaussian[0][0],0);
     2823
     2824                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
     2825                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2826        }
    28752827
    28762828
    … …  
    28802832        //Deal with 2d friction at the bedrock interface
    28812833        if((onbed && !shelf)){
    2882 
    28832834                /*Build a tria element using the 3 grids of the base of the penta. Then use
    28842835                 * the tria functionality to build a friction stiffness matrix on these 3
    … …  
    28902841        }
    28912842
    2892         xfree((void**)&first_gauss_area_coord);
    2893         xfree((void**)&second_gauss_area_coord);
    2894         xfree((void**)&third_gauss_area_coord);
    2895         xfree((void**)&fourth_gauss_vert_coord);
    2896         xfree((void**)&area_gauss_weights);
    2897         xfree((void**)&vert_gauss_weights);
    2898         xfree((void**)&first_gauss_area_coord2d);
    2899         xfree((void**)&second_gauss_area_coord2d);
    2900         xfree((void**)&third_gauss_area_coord2d);
    2901         xfree((void**)&gauss_weights2d);
     2843        /*Clean up and return*/
     2844        delete gauss;
    29022845        xfree((void**)&doflist);
    29032846}
    … …  
    52085151/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
    52095152void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){
     5153        /*Compute the 3d Blatter/PattynStrain Rate (5 components):
     5154         *
     5155         * epsilon=[exx eyy exy exz eyz]
     5156         *
     5157         * with exz=1/2 du/dz
     5158         *      eyz=1/2 dv/dz
     5159         *
     5160         * the contribution of vz is neglected
     5161         */
     5162
     5163        int i;
     5164
     5165        double epsilonvx[5];
     5166        double epsilonvy[5];
     5167
     5168        /*Check that both inputs have been found*/
     5169        if (!vx_input || !vy_input){
     5170                ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
     5171        }
     5172
     5173        /*Get strain rate assuming that epsilon has been allocated*/
     5174        vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
     5175        vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     5176
     5177        /*Sum all contributions*/
     5178        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     5179
     5180}
     5181/*}}}*/
     5182/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
     5183void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
    52105184        /*Compute the 3d Blatter/PattynStrain Rate (5 components):
    52115185         *

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

    r5635 r5647  
    168168                void      GetSolutionFromInputsThermal(Vec solutiong);
    169169                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
     170                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    170171                void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    171172                Penta*  GetUpperElement(void);

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

    r5203 r5647  
    981981}
    982982/*}}}*/
     983
     984/*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/
     985void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){
     986        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     987         * For grid i, Bi can be expressed in the actual coordinate system
     988         * by:
     989         *       Bi=[ dh/dx          0      ]
     990         *          [   0           dh/dy   ]
     991         *          [ 1/2*dh/dy  1/2*dh/dx  ]
     992         * where h is the interpolation function for grid i.
     993         *
     994         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
     995         */
     996
     997        int i;
     998        const int numgrids=6;
     999        const int NDOF3=3;
     1000        const int NDOF2=2;
     1001
     1002        double dh1dh6[NDOF3][numgrids];
     1003
     1004        /*Get dh1dh6 in actual coordinate system: */
     1005        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     1006
     1007        /*Build B: */
     1008        for (i=0;i<numgrids;i++){
     1009                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
     1010                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
     1011
     1012                *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
     1013                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
     1014
     1015                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
     1016                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
     1017
     1018        }
     1019
     1020}
     1021/*}}}*/
     1022/*FUNCTION PentaRef::GetBPattyn {{{1*/
     1023void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){
     1024        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     1025         * For grid i, Bi can be expressed in the actual coordinate system
     1026         * by:
     1027         *       Bi=[ dh/dx          0      ]
     1028         *          [   0           dh/dy   ]
     1029         *          [ 1/2*dh/dy  1/2*dh/dx  ]
     1030         *          [ 1/2*dh/dz      0      ]
     1031         *          [  0         1/2*dh/dz  ]
     1032         * where h is the interpolation function for grid i.
     1033         *
     1034         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
     1035         */
     1036
     1037        int i;
     1038        const int numgrids=6;
     1039        const int NDOF3=3;
     1040        const int NDOF2=2;
     1041
     1042        double dh1dh6[NDOF3][numgrids];
     1043
     1044        /*Get dh1dh6 in actual coordinate system: */
     1045        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     1046
     1047        /*Build B: */
     1048        for (i=0;i<numgrids;i++){
     1049                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
     1050                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
     1051
     1052                *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
     1053                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
     1054
     1055                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
     1056                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
     1057
     1058                *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
     1059                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
     1060
     1061                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
     1062                *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
     1063        }
     1064
     1065}
     1066/*}}}*/
     1067/*FUNCTION PentaRef::GetBprimePattyn {{{1*/
     1068void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){
     1069        /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     1070         * For grid i, Bi can be expressed in the actual coordinate system
     1071         * by:
     1072         *       Bi=[ 2*dh/dx     dh/dy   ]
     1073         *                [   dh/dx    2*dh/dy  ]
     1074         *                [ dh/dy      dh/dx    ]
     1075         *                [ dh/dz         0     ]
     1076         *                [  0         dh/dz    ]
     1077         * where h is the interpolation function for grid i.
     1078         *
     1079         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
     1080         */
     1081
     1082        int i;
     1083        const int NDOF3=3;
     1084        const int NDOF2=2;
     1085        const int numgrids=6;
     1086
     1087        double dh1dh6[NDOF3][numgrids];
     1088
     1089        /*Get dh1dh6 in actual coordinate system: */
     1090        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
     1091
     1092        /*Build BPrime: */
     1093        for (i=0;i<numgrids;i++){
     1094                *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];
     1095                *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
     1096
     1097                *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
     1098                *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
     1099
     1100                *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i];
     1101                *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i];
     1102
     1103                *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i];
     1104                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
     1105
     1106                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
     1107                *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i];
     1108        }
     1109}
     1110/*}}}*/
     1111/*FUNCTION PentaRef::GetBStokes {{{1*/
     1112void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){
     1113
     1114        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
     1115         * For grid i, Bi can be expressed in the actual coordinate system
     1116         * by:          Bi=[ dh/dx          0              0       0  ]
     1117         *                                      [   0           dh/dy           0       0  ]
     1118         *                                      [   0             0           dh/dy     0  ]
     1119         *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
     1120         *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
     1121         *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
     1122         *                                      [   0             0             0       h  ]
     1123         *                                      [ dh/dx         dh/dy         dh/dz     0  ]
     1124         *      where h is the interpolation function for grid i.
     1125         *      Same thing for Bb except the last column that does not exist.
     1126         */
     1127
     1128        int i;
     1129        const int calculationdof=3;
     1130        const int numgrids=6;
     1131        int DOFPERGRID=4;
     1132
     1133        double dh1dh7[calculationdof][numgrids+1];
     1134        double l1l6[numgrids];
     1135
     1136
     1137        /*Get dh1dh7 in actual coordinate system: */
     1138        GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
     1139        GetNodalFunctionsP1(l1l6, gauss);
     1140
     1141        /*Build B: */
     1142        for (i=0;i<numgrids+1;i++){
     1143                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
     1144                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
     1145                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
     1146                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
     1147                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
     1148                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
     1149                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
     1150                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
     1151                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
     1152                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
     1153                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
     1154                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
     1155                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
     1156                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
     1157                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
     1158                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
     1159                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
     1160                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
     1161                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
     1162                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
     1163                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
     1164                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
     1165                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
     1166                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
     1167        }
     1168
     1169        for (i=0;i<numgrids;i++){ //last column not for the bubble function
     1170                *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
     1171                *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
     1172                *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
     1173                *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
     1174                *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
     1175                *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
     1176                *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
     1177                *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
     1178        }
     1179
     1180}
     1181/*}}}*/
     1182/*FUNCTION PentaRef::GetBprimeStokes {{{1*/
     1183void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){
     1184        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     1185         *      For grid i, Bi' can be expressed in the actual coordinate system
     1186         *      by:
     1187         *                              Bi'=[  dh/dx   0          0       0]
     1188         *                                       [   0      dh/dy      0       0]
     1189         *                                       [   0      0         dh/dz    0]
     1190         *                                       [  dh/dy   dh/dx      0       0]
     1191         *                                       [  dh/dz   0        dh/dx     0]
     1192         *                                       [   0      dh/dz    dh/dy     0]
     1193         *                                       [  dh/dx   dh/dy    dh/dz     0]
     1194         *                                       [   0      0          0       h]
     1195         *      where h is the interpolation function for grid i.
     1196         *
     1197         *      Same thing for the bubble fonction except that there is no fourth column
     1198         */
     1199
     1200        int i;
     1201        const int calculationdof=3;
     1202        const int numgrids=6;
     1203        int DOFPERGRID=4;
     1204
     1205        double dh1dh7[calculationdof][numgrids+1];
     1206        double l1l6[numgrids];
     1207
     1208        /*Get dh1dh7 in actual coordinate system: */
     1209        GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
     1210
     1211        GetNodalFunctionsP1(l1l6, gauss);
     1212
     1213        /*B_primeuild B_prime: */
     1214        for (i=0;i<numgrids+1;i++){
     1215                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
     1216                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
     1217                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
     1218                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
     1219                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
     1220                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
     1221                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
     1222                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
     1223                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
     1224                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];
     1225                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];
     1226                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
     1227                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
     1228                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
     1229                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
     1230                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
     1231                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
     1232                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
     1233                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
     1234                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
     1235                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
     1236                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
     1237                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
     1238                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;
     1239        }
     1240
     1241        for (i=0;i<numgrids;i++){ //last column not for the bubble function
     1242                *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
     1243                *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
     1244                *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
     1245                *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
     1246                *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
     1247                *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
     1248                *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;
     1249                *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];
     1250        }
     1251
     1252}
     1253/*}}}*/
     1254/*FUNCTION PentaRef::GetBArtdiff {{{1*/
     1255void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){
     1256        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     1257         * For grid i, Bi' can be expressed in the actual coordinate system
     1258         * by:
     1259         *       Bi_artdiff=[ dh/dx ]
     1260         *                 [ dh/dy ]
     1261         * where h is the interpolation function for grid i.
     1262         *
     1263         * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
     1264         */
     1265
     1266        int i;
     1267        const int calculationdof=3;
     1268        const int numgrids=6;
     1269        int DOFPERGRID=1;
     1270
     1271        /*Same thing in the actual coordinate system: */
     1272        double dh1dh6[calculationdof][numgrids];
     1273
     1274        /*Get dh1dh6 in actual coordinates system : */
     1275        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     1276
     1277        /*Build B': */
     1278        for (i=0;i<numgrids;i++){
     1279                *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     1280                *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     1281        }
     1282}
     1283/*}}}*/
     1284/*FUNCTION PentaRef::GetBAdvec{{{1*/
     1285void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){
     1286        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     1287         * For grid i, Bi' can be expressed in the actual coordinate system
     1288         * by:
     1289         *       Bi_advec =[ h ]
     1290         *                 [ h ]
     1291         *                 [ h ]
     1292         * where h is the interpolation function for grid i.
     1293         *
     1294         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
     1295         */
     1296
     1297        int i;
     1298        int calculationdof=3;
     1299        int numgrids=6;
     1300        int DOFPERGRID=1;
     1301
     1302        /*Same thing in the actual coordinate system: */
     1303        double l1l6[6];
     1304
     1305        /*Get dh1dh2dh3 in actual coordinates system : */
     1306        GetNodalFunctionsP1(l1l6, gauss);
     1307
     1308        /*Build B': */
     1309        for (i=0;i<numgrids;i++){
     1310                *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i];
     1311                *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i];
     1312                *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i];
     1313        }
     1314}
     1315/*}}}*/
     1316/*FUNCTION PentaRef::GetBConduct{{{1*/
     1317void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){
     1318        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     1319         * For grid i, Bi' can be expressed in the actual coordinate system
     1320         * by:
     1321         *       Bi_conduct=[ dh/dx ]
     1322         *                  [ dh/dy ]
     1323         *                  [ dh/dz ]
     1324         * where h is the interpolation function for grid i.
     1325         *
     1326         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
     1327         */
     1328
     1329        int i;
     1330        const int calculationdof=3;
     1331        const int numgrids=6;
     1332        int DOFPERGRID=1;
     1333
     1334        /*Same thing in the actual coordinate system: */
     1335        double dh1dh6[calculationdof][numgrids];
     1336
     1337        /*Get dh1dh2dh3 in actual coordinates system : */
     1338        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     1339
     1340        /*Build B': */
     1341        for (i=0;i<numgrids;i++){
     1342                *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     1343                *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     1344                *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
     1345        }
     1346}
     1347/*}}}*/
     1348/*FUNCTION PentaRef::GetBVert{{{1*/
     1349void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){
     1350        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
     1351                where hi is the interpolation function for grid i.*/
     1352
     1353        int i;
     1354        const int NDOF3=3;
     1355        const int numgrids=6;
     1356        double dh1dh6[NDOF3][numgrids];
     1357
     1358        /*Get dh1dh6 in actual coordinate system: */
     1359        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     1360
     1361        /*Build B: */
     1362        for (i=0;i<numgrids;i++){
     1363                B[i]=dh1dh6[2][i]; 
     1364        }
     1365
     1366}
     1367/*}}}*/
     1368/*FUNCTION PentaRef::GetBprimeAdvec{{{1*/
     1369void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){
     1370        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
     1371         * For grid i, Bi' can be expressed in the actual coordinate system
     1372         * by:
     1373         *       Biprime_advec=[ dh/dx ]
     1374         *                     [ dh/dy ]
     1375         *                     [ dh/dz ]
     1376         * where h is the interpolation function for grid i.
     1377         *
     1378         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
     1379         */
     1380
     1381        int i;
     1382        const int calculationdof=3;
     1383        const int numgrids=6;
     1384        int DOFPERGRID=1;
     1385
     1386        /*Same thing in the actual coordinate system: */
     1387        double dh1dh6[calculationdof][numgrids];
     1388
     1389        /*Get dh1dh2dh3 in actual coordinates system : */
     1390        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     1391
     1392        /*Build B': */
     1393        for (i=0;i<numgrids;i++){
     1394                *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
     1395                *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
     1396                *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
     1397        }
     1398}
     1399/*}}}*/
     1400/*FUNCTION PentaRef::GetBprimeVert{{{1*/
     1401void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){
     1402        /* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/
     1403
     1404        GetNodalFunctionsP1(B, gauss);
     1405
     1406}
     1407/*}}}*/
     1408/*FUNCTION PentaRef::GetLprimeStokes {{{1*/
     1409void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss){
     1410
     1411        /*
     1412         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     1413         * For grid i, Lpi can be expressed in the actual coordinate system
     1414         * by:
     1415         *       Lpi=[ h    0    0   0]
     1416         *                     [ 0    h    0   0]
     1417         *                     [ h    0    0   0]
     1418         *                     [ 0    h    0   0]
     1419         *                     [ 0    0    h   0]
     1420         *                     [ 0    0    h   0]
     1421         *                     [ 0    0  dh/dz 0]
     1422         *                     [ 0    0  dh/dz 0]
     1423         *                     [ 0    0  dh/dz 0]
     1424         *                     [dh/dz 0  dh/dx 0]
     1425         *                     [ 0 dh/dz dh/dy 0]
     1426         *           [ 0    0    0   h]
     1427         *           [ 0    0    0   h]
     1428         *           [ 0    0    0   h]
     1429         * where h is the interpolation function for grid i.
     1430         */
     1431        int i;
     1432        const int numgrids2d=3;
     1433        int num_dof=4;
     1434
     1435        double l1l2l3[numgrids2d];
     1436        double dh1dh6[3][6];
     1437
     1438        /*Get l1l2l3 in actual coordinate system: */
     1439        l1l2l3[0]=gauss_tria[0];
     1440        l1l2l3[1]=gauss_tria[1];
     1441        l1l2l3[2]=gauss_tria[2];
     1442
     1443        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
     1444
     1445        /*Build LprimeStokes: */
     1446        for (i=0;i<3;i++){
     1447                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
     1448                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
     1449                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
     1450                *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
     1451                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;
     1452                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
     1453                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
     1454                *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
     1455                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];
     1456                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
     1457                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;
     1458                *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
     1459                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;
     1460                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];
     1461                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;
     1462                *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
     1463                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;
     1464                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
     1465                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];
     1466                *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
     1467                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;
     1468                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;
     1469                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];
     1470                *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
     1471                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
     1472                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
     1473                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
     1474                *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
     1475                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
     1476                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
     1477                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
     1478                *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
     1479                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
     1480                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
     1481                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
     1482                *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
     1483                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
     1484                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
     1485                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
     1486                *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
     1487                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
     1488                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
     1489                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
     1490                *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
     1491                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
     1492                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
     1493                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
     1494                *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];
     1495                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;
     1496                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;
     1497                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
     1498                *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];
     1499                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;
     1500                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
     1501                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
     1502                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
     1503
     1504        }
     1505}
     1506/*}}}*/
     1507/*FUNCTION PentaRef::GetJacobian {{{1*/
     1508void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){
     1509
     1510        const int NDOF3=3;
     1511        int i,j;
     1512
     1513        /*The Jacobian is constant over the element, discard the gaussian points.
     1514         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1515
     1516        double A1,A2,A3; //area coordinates
     1517        double xi,eta,zi; //parametric coordinates
     1518
     1519        double x1,x2,x3,x4,x5,x6;
     1520        double y1,y2,y3,y4,y5,y6;
     1521        double z1,z2,z3,z4,z5,z6;
     1522
     1523        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
     1524        A1=gauss->coord1;
     1525        A2=gauss->coord2;
     1526        A3=gauss->coord3;
     1527
     1528        xi=A2-A1;
     1529        eta=SQRT3*A3;
     1530        zi=gauss->coord4;
     1531
     1532        x1=*(xyz_list+3*0+0);
     1533        x2=*(xyz_list+3*1+0);
     1534        x3=*(xyz_list+3*2+0);
     1535        x4=*(xyz_list+3*3+0);
     1536        x5=*(xyz_list+3*4+0);
     1537        x6=*(xyz_list+3*5+0);
     1538
     1539        y1=*(xyz_list+3*0+1);
     1540        y2=*(xyz_list+3*1+1);
     1541        y3=*(xyz_list+3*2+1);
     1542        y4=*(xyz_list+3*3+1);
     1543        y5=*(xyz_list+3*4+1);
     1544        y6=*(xyz_list+3*5+1);
     1545
     1546        z1=*(xyz_list+3*0+2);
     1547        z2=*(xyz_list+3*1+2);
     1548        z3=*(xyz_list+3*2+2);
     1549        z4=*(xyz_list+3*3+2);
     1550        z5=*(xyz_list+3*4+2);
     1551        z6=*(xyz_list+3*5+2);
     1552
     1553        *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
     1554        *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
     1555        *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
     1556
     1557        *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
     1558        *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
     1559        *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
     1560
     1561        *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
     1562        *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
     1563        *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
     1564
     1565}
     1566/*}}}*/
     1567/*FUNCTION PentaRef::GetJacobianDeterminant {{{1*/
     1568void PentaRef::GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss){
     1569        /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
     1570         * the determinant of it: */
     1571        double J[3][3];
     1572
     1573        /*Get Jacobian*/
     1574        GetJacobian(&J[0][0],xyz_list,gauss);
     1575
     1576        /*Get Determinant*/
     1577        Matrix3x3Determinant(Jdet,&J[0][0]);
     1578        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     1579
     1580}
     1581/*}}}*/
     1582/*FUNCTION PentaRef::GetJacobianInvert {{{1*/
     1583void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){
     1584
     1585        /*Jacobian*/
     1586        double J[3][3];
     1587
     1588        /*Call Jacobian routine to get the jacobian:*/
     1589        GetJacobian(&J[0][0], xyz_list, gauss);
     1590
     1591        /*Invert Jacobian matrix: */
     1592        Matrix3x3Invert(Jinv,&J[0][0]);
     1593}
     1594/*}}}*/
     1595/*FUNCTION PentaRef::GetNodalFunctionsMINI{{{1*/
     1596void PentaRef::GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss){
     1597        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     1598
     1599        l1l7[0]=gauss->coord1*(1.0-gauss->coord4)/2.0;
     1600        l1l7[1]=gauss->coord2*(1.0-gauss->coord4)/2.0;
     1601        l1l7[2]=gauss->coord3*(1.0-gauss->coord4)/2.0;
     1602        l1l7[3]=gauss->coord1*(1.0+gauss->coord4)/2.0;
     1603        l1l7[4]=gauss->coord2*(1.0+gauss->coord4)/2.0;
     1604        l1l7[5]=gauss->coord3*(1.0+gauss->coord4)/2.0;
     1605        l1l7[6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(1.0+gauss->coord4)*(1.0-gauss->coord4);
     1606
     1607}
     1608/*}}}*/
     1609/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/
     1610void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss){
     1611
     1612        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     1613         * actual coordinate system): */
     1614
     1615        int       i;
     1616        const int numgrids = 7;
     1617        double    dh1dh7_ref[3][numgrids];
     1618        double    Jinv[3][3];
     1619
     1620        /*Get derivative values with respect to parametric coordinate system: */
     1621        GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss);
     1622
     1623        /*Get Jacobian invert: */
     1624        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
     1625
     1626        /*Build dh1dh6:
     1627         *
     1628         * [dhi/dx]= Jinv'*[dhi/dr]
     1629         * [dhi/dy]        [dhi/ds]
     1630         * [dhi/dz]        [dhi/dzeta]
     1631         */
     1632
     1633        for (i=0;i<numgrids;i++){
     1634                *(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];
     1635                *(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];
     1636                *(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];
     1637        }
     1638
     1639}
     1640/*}}}*/
     1641/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{1*/
     1642void PentaRef::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss){
     1643
     1644        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     1645         * natural coordinate system) at the gaussian point. */
     1646
     1647        int    numgrids=7; //six plus bubble grids
     1648
     1649        double r=gauss->coord2-gauss->coord1;
     1650        double s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
     1651        double zeta=gauss->coord4;
     1652
     1653        /*First nodal function: */
     1654        *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
     1655        *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
     1656        *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
     1657
     1658        /*Second nodal function: */
     1659        *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
     1660        *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
     1661        *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
     1662
     1663        /*Third nodal function: */
     1664        *(dl1dl7+numgrids*0+2)=0;
     1665        *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
     1666        *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
     1667
     1668        /*Fourth nodal function: */
     1669        *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
     1670        *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
     1671        *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
     1672
     1673        /*Fith nodal function: */
     1674        *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
     1675        *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
     1676        *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
     1677
     1678        /*Sixth nodal function: */
     1679        *(dl1dl7+numgrids*0+5)=0;
     1680        *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
     1681        *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
     1682
     1683        /*Seventh nodal function: */
     1684        *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
     1685        *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
     1686        *(dl1dl7+numgrids*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
     1687
     1688}
     1689/*}}}*/
     1690/*FUNCTION PentaRef::GetNodalFunctionsP1 {{{1*/
     1691void PentaRef::GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss){
     1692        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     1693
     1694        l1l6[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     1695        l1l6[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     1696        l1l6[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     1697        l1l6[3]=gauss->coord1*(1+gauss->coord4)/2.0;
     1698        l1l6[4]=gauss->coord2*(1+gauss->coord4)/2.0;
     1699        l1l6[5]=gauss->coord3*(1+gauss->coord4)/2.0;
     1700
     1701}
     1702/*}}}*/
     1703/*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{1*/
     1704void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss){
     1705
     1706        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     1707         * actual coordinate system): */
     1708        int       i;
     1709        const int NDOF3    = 3;
     1710        const int numgrids = 6;
     1711        double    dh1dh6_ref[NDOF3][numgrids];
     1712        double    Jinv[NDOF3][NDOF3];
     1713
     1714        /*Get derivative values with respect to parametric coordinate system: */
     1715        GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss);
     1716
     1717        /*Get Jacobian invert: */
     1718        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
     1719
     1720        /*Build dh1dh3:
     1721         *
     1722         * [dhi/dx]= Jinv*[dhi/dr]
     1723         * [dhi/dy]       [dhi/ds]
     1724         * [dhi/dz]       [dhi/dn]
     1725         */
     1726
     1727        for (i=0;i<numgrids;i++){
     1728                *(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];
     1729                *(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];
     1730                *(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];
     1731        }
     1732
     1733}
     1734/*}}}*/
     1735/*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{1*/
     1736void PentaRef::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss){
     1737
     1738        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     1739         * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
     1740
     1741        const int numgrids=6;
     1742        double A1,A2,A3,z;
     1743
     1744        A1=gauss->coord1; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
     1745        A2=gauss->coord2; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
     1746        A3=gauss->coord3; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value  In term of xi and eta: A3=y/SQRT3;
     1747        z =gauss->coord4; ISSMASSERT(z>=-1 &&  z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
     1748
     1749        /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
     1750        *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
     1751        *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
     1752        *(dl1dl6+numgrids*2+0)=-0.5*A1;
     1753
     1754        /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
     1755        *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
     1756        *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
     1757        *(dl1dl6+numgrids*2+1)=-0.5*A2;
     1758
     1759        /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
     1760        *(dl1dl6+numgrids*0+2)=0.0;
     1761        *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
     1762        *(dl1dl6+numgrids*2+2)=-0.5*A3;
     1763
     1764        /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
     1765        *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
     1766        *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
     1767        *(dl1dl6+numgrids*2+3)=0.5*A1;
     1768
     1769        /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
     1770        *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
     1771        *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
     1772        *(dl1dl6+numgrids*2+4)=0.5*A2;
     1773
     1774        /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
     1775        *(dl1dl6+numgrids*0+5)=0.0;
     1776        *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
     1777        *(dl1dl6+numgrids*2+5)=0.5*A3;
     1778}
     1779/*}}}*/
     1780/*FUNCTION PentaRef::GetParameterValue{{{1*/
     1781void PentaRef::GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss){
     1782        /*P1 interpolation on Gauss point*/
     1783
     1784        /*intermediary*/
     1785        double l1l6[6];
     1786
     1787        /*nodal functions: */
     1788        GetNodalFunctionsP1(&l1l6[0],gauss);
     1789
     1790        /*Assign output pointers:*/
     1791        *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];
     1792
     1793}
     1794/*}}}*/
     1795/*FUNCTION PentaRef::GetParameterDerivativeValue{{{1*/
     1796void PentaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussPenta* gauss){
     1797        /*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:
     1798         *   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;
     1799         *   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;
     1800         *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
     1801         *
     1802         *   p is a vector of size 3x1 already allocated.
     1803         */
     1804        double dh1dh6[3][6];
     1805
     1806        /*Get nodal funnctions derivatives in actual coordinate system: */
     1807        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
     1808
     1809        /*Assign output*/
     1810        p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];
     1811        p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];
     1812        p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];
     1813
     1814}
     1815/*}}}*/

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

    r5629 r5647  
    5050                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    5151
     52                void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss);
     53                void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss);
     54                void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, GaussPenta* gauss);
     55                void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, GaussPenta* gauss);
     56                void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss);
     57                void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss);
     58                void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss);
     59                void GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
     60                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussPenta* gauss);
     61                void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss);
     62                void GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss);
     63                void GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss);
     64                void GetBStokes(double* B, double* xyz_list, GaussPenta* gauss);
     65                void GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss);
     66                void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss);
     67                void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss);
     68                void GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss);
     69                void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss);
     70                void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
     71                void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss);
     72                //void GetLStokes(double* LStokes, double* gauss_tria);
     73                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss);
     74                void GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss);
     75                //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
     76                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss);
     77                //void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
     78
    5279};
    5380#endif

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

    r5629 r5647  
    172172void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    173173/*}}}*/
     174/*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
     175void BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
     176/*}}}*/
    174177/*FUNCTION BoolInput::GetParameterValues{{{1*/
    175178void BoolInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    … …  
    180183/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
    181184void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
     185/*}}}*/
     186/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/
     187void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    182188/*}}}*/
    183189/*FUNCTION BoolInput::ChangeEnum{{{1*/

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

    r5631 r5647  
    4949                void GetParameterValue(double* pvalue,double* gauss);
    5050                void GetParameterValue(double* pvalue,GaussTria* gauss);
     51                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    5152                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5253                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5354                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
     55                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);
    5456                void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");};
    5557                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};
    … …  
    6264                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    6365                void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    64                 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    65                 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    66                 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    67                 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    68                 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     66                void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     67                void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     68                void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     69                void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     70                void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
    6971                void ChangeEnum(int newenumtype);
    7072                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);

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

    r5629 r5647  
    185185void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    186186/*}}}*/
     187/*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
     188void DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
     189/*}}}*/
    187190/*FUNCTION DoubleInput::GetParameterValues{{{1*/
    188191void DoubleInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    … …  
    193196/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
    194197void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
     198/*}}}*/
     199/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/
     200void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    195201/*}}}*/
    196202/*FUNCTION DoubleInput::ChangeEnum{{{1*/

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

    r5631 r5647  
    4848                void GetParameterValue(double* pvalue,double* gauss);
    4949                void GetParameterValue(double* pvalue,GaussTria* gauss);
     50                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    5051                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5152                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5253                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
     54                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);
    5355                void GetParameterAverage(double* pvalue);
    5456                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};
    … …  
    6163                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    6264                void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    63                 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    64                 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    65                 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    66                 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    67                 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     65                void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     66                void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     67                void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     68                void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     69                void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
    6870                void ChangeEnum(int newenumtype);
    6971                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);

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

    r5631 r5647  
    2727                virtual void GetParameterValue(double* pvalue,double* gauss)=0;
    2828                virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0;
     29                virtual void GetParameterValue(double* pvalue,GaussPenta* gauss)=0;
    2930                virtual void GetParameterValues(double* values,double* gauss_pointers, int numgauss)=0;
    3031                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0;
    3132                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0;
     33                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss)=0;
    3234                virtual void GetParameterAverage(double* pvalue)=0;
    3335                virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss)=0;
    … …  
    4042                virtual void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0;
    4143                virtual void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0;
    42                 virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0;
    43                 virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0;
    44                 virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss)=0;
    45                 virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss)=0;
    46                 virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss)=0;
     44                virtual void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0;
     45                virtual void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0;
     46                virtual void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss)=0;
     47                virtual void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss)=0;
     48                virtual void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss)=0;
    4749                virtual void ChangeEnum(int newenumtype)=0;
    4850

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

    r5629 r5647  
    173173void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    174174/*}}}*/
     175/*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
     176void IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
     177/*}}}*/
    175178/*FUNCTION IntInput::GetParameterValues{{{1*/
    176179void IntInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){ISSMERROR(" not supported yet!");}
    … …  
    181184/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/
    182185void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR(" not supported yet!");}
     186/*}}}*/
     187/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){{{1*/
     188void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    183189/*}}}*/
    184190/*FUNCTION IntInput::ChangeEnum{{{1*/

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

    r5631 r5647  
    4949                void GetParameterValue(double* pvalue,double* gauss);
    5050                void GetParameterValue(double* pvalue,GaussTria* gauss);
     51                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    5152                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5253                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5354                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
     55                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);
    5456                void GetParameterAverage(double* pvalue){ISSMERROR("not implemented yet");};
    5557                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};
    … …  
    6264                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    6365                void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    64                 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    65                 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    66                 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    67                 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    68                 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     66                void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     67                void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     68                void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     69                void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     70                void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
    6971                void ChangeEnum(int newenumtype);
    7072                void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);

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

    r5629 r5647  
    184184}
    185185/*}}}*/
    186 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
    187 void PentaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){
     186/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
     187void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){
    188188
    189189        /*Call PentaRef function*/
    … …  
    216216}
    217217/*}}}*/
    218 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
    219 void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
     218/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/
     219void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){
    220220
    221221        /*Call PentaRef function*/
    … …  
    386386/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
    387387void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){
     388
     389        int i;
     390        const int numgrids=6;
     391        const int NDOF2=2;
     392        double B[5][NDOF2*numgrids];
     393        double velocity[numgrids][NDOF2];
     394
     395        /*Get B matrix: */
     396        GetBPattyn(&B[0][0], xyz_list, gauss);
     397
     398        /*Here, we are computing the strain rate of (0,vy)*/
     399        for(i=0;i<numgrids;i++){
     400                velocity[i][0]=0.0;
     401                velocity[i][1]=this->values[i];
     402        }
     403
     404        /*Multiply B by velocity, to get strain rate: */
     405        MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
     406                                &velocity[0][0],NDOF2*numgrids,1,0,
     407                                epsilonvy,0);
     408
     409}
     410/*}}}*/
     411/*FUNCTION PentaVertexInput::GetVxStrainRate3d{{{1*/
     412void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
     413        int i,j;
     414
     415        const int numgrids=6;
     416        const int DOFVELOCITY=3;
     417        double B[8][27];
     418        double B_reduced[6][DOFVELOCITY*numgrids];
     419        double velocity[numgrids][DOFVELOCITY];
     420
     421        /*Get B matrix: */
     422        GetBStokes(&B[0][0], xyz_list, gauss);
     423        /*Create a reduced matrix of B to get rid of pressure */
     424        for (i=0;i<6;i++){
     425                for (j=0;j<3;j++){
     426                        B_reduced[i][j]=B[i][j];
     427                }
     428                for (j=4;j<7;j++){
     429                        B_reduced[i][j-1]=B[i][j];
     430                }
     431                for (j=8;j<11;j++){
     432                        B_reduced[i][j-2]=B[i][j];
     433                }
     434                for (j=12;j<15;j++){
     435                        B_reduced[i][j-3]=B[i][j];
     436                }
     437                for (j=16;j<19;j++){
     438                        B_reduced[i][j-4]=B[i][j];
     439                }
     440                for (j=20;j<23;j++){
     441                        B_reduced[i][j-5]=B[i][j];
     442                }
     443        }
     444
     445        /*Here, we are computing the strain rate of (vx,0,0)*/
     446        for(i=0;i<numgrids;i++){
     447                velocity[i][0]=this->values[i];
     448                velocity[i][1]=0.0;
     449                velocity[i][2]=0.0;
     450        }
     451        /*Multiply B by velocity, to get strain rate: */
     452        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);
     453
     454}
     455/*}}}*/
     456/*FUNCTION PentaVertexInput::GetVyStrainRate3d{{{1*/
     457void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
     458        int i,j;
     459
     460        const int numgrids=6;
     461        const int DOFVELOCITY=3;
     462        double B[8][27];
     463        double B_reduced[6][DOFVELOCITY*numgrids];
     464        double velocity[numgrids][DOFVELOCITY];
     465
     466        /*Get B matrix: */
     467        GetBStokes(&B[0][0], xyz_list, gauss);
     468        /*Create a reduced matrix of B to get rid of pressure */
     469        for (i=0;i<6;i++){
     470                for (j=0;j<3;j++){
     471                        B_reduced[i][j]=B[i][j];
     472                }
     473                for (j=4;j<7;j++){
     474                        B_reduced[i][j-1]=B[i][j];
     475                }
     476                for (j=8;j<11;j++){
     477                        B_reduced[i][j-2]=B[i][j];
     478                }
     479                for (j=12;j<15;j++){
     480                        B_reduced[i][j-3]=B[i][j];
     481                }
     482                for (j=16;j<19;j++){
     483                        B_reduced[i][j-4]=B[i][j];
     484                }
     485                for (j=20;j<23;j++){
     486                        B_reduced[i][j-5]=B[i][j];
     487                }
     488        }
     489
     490        /*Here, we are computing the strain rate of (0,vy,0)*/
     491        for(i=0;i<numgrids;i++){
     492                velocity[i][0]=0.0;
     493                velocity[i][1]=this->values[i];
     494                velocity[i][2]=0.0;
     495        }
     496        /*Multiply B by velocity, to get strain rate: */
     497        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);
     498
     499}
     500/*}}}*/
     501/*FUNCTION PentaVertexInput::GetVzStrainRate3d{{{1*/
     502void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){
     503        int i,j;
     504
     505        const int numgrids=6;
     506        const int DOFVELOCITY=3;
     507        double B[8][27];
     508        double B_reduced[6][DOFVELOCITY*numgrids];
     509        double velocity[numgrids][DOFVELOCITY];
     510
     511        /*Get B matrix: */
     512        GetBStokes(&B[0][0], xyz_list, gauss);
     513        /*Create a reduced matrix of B to get rid of pressure */
     514        for (i=0;i<6;i++){
     515                for (j=0;j<3;j++){
     516                        B_reduced[i][j]=B[i][j];
     517                }
     518                for (j=4;j<7;j++){
     519                        B_reduced[i][j-1]=B[i][j];
     520                }
     521                for (j=8;j<11;j++){
     522                        B_reduced[i][j-2]=B[i][j];
     523                }
     524                for (j=12;j<15;j++){
     525                        B_reduced[i][j-3]=B[i][j];
     526                }
     527                for (j=16;j<19;j++){
     528                        B_reduced[i][j-4]=B[i][j];
     529                }
     530                for (j=20;j<23;j++){
     531                        B_reduced[i][j-5]=B[i][j];
     532                }
     533        }
     534
     535        /*Here, we are computing the strain rate of (0,0,vz)*/
     536        for(i=0;i<numgrids;i++){
     537                velocity[i][0]=0.0;
     538                velocity[i][1]=0.0;
     539                velocity[i][2]=this->values[i];
     540        }
     541
     542        /*Multiply B by velocity, to get strain rate: */
     543        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);
     544
     545}
     546/*}}}*/
     547/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn{{{1*/
     548void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
     549
     550        int i;
     551        const int numgrids=6;
     552        const int NDOF2=2;
     553        double B[5][NDOF2*numgrids];
     554        double velocity[numgrids][NDOF2];
     555
     556        /*Get B matrix: */
     557        GetBPattyn(&B[0][0], xyz_list, gauss);
     558
     559        /*Here, we are computing the strain rate of (vx,0)*/
     560        for(i=0;i<numgrids;i++){
     561                velocity[i][0]=this->values[i];
     562                velocity[i][1]=0.0;
     563        }
     564
     565        /*Multiply B by velocity, to get strain rate: */
     566        MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
     567                                &velocity[0][0],NDOF2*numgrids,1,0,
     568                                epsilonvx,0);
     569
     570}
     571/*}}}*/
     572/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn{{{1*/
     573void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
    388574
    389575        int i;

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

    r5631 r5647  
    4848                void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");};
    4949                void GetParameterValue(double* pvalue,double* gauss);
    50                 void GetParameterValue(double* pvalue,GaussTria* gauss);
     50                void GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR("not implemented yet");};
     51                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    5152                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5253                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    53                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
     54                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     55                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);
    5456                void GetParameterAverage(double* pvalue);
    5557                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss){ISSMERROR("not implemented yet");};
    … …  
    6264                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    6365                void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    64                 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    65                 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    66                 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    67                 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    68                 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     66                void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss);
     67                void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss);
     68                void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss);
     69                void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss);
     70                void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss);
    6971                void ChangeEnum(int newenumtype);
    7072

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

    r5631 r5647  
    4949                void GetParameterValue(double* pvalue,double* gauss);
    5050                void GetParameterValue(double* pvalue,GaussTria* gauss);
     51                void GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR("not implemented yet");};
    5152                void GetParameterValues(double* values,double* gauss_pointers, int numgauss);
    5253                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5354                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
     55                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
    5456                void GetParameterAverage(double* pvalue);
    5557                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, double* gauss);
    … …  
    6264                void GetVxStrainRate2d(double* epsilonvx,double* xyz_list, GaussTria* gauss);
    6365                void GetVyStrainRate2d(double* epsilonvy,double* xyz_list, GaussTria* gauss);
    64                 void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    65                 void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    66                 void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    67                 void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    68                 void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
     66                void GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     67                void GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     68                void GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     69                void GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
     70                void GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};
    6971                void ChangeEnum(int newenumtype);
    7072
Note: See TracChangeset for help on using the changeset viewer.