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

Added GaussPenta for Pattyn

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.