Changeset 5647 for issm/trunk/src/c/objects/Elements/PentaRef.cpp
- Timestamp:
- 09/01/10 10:45:00 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5203 r5647 981 981 } 982 982 /*}}}*/ 983 984 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/ 985 void 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*/ 1023 void 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*/ 1068 void 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*/ 1112 void 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*/ 1183 void 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*/ 1255 void 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*/ 1285 void 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*/ 1317 void 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*/ 1349 void 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*/ 1369 void 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*/ 1401 void 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*/ 1409 void 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*/ 1508 void 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*/ 1568 void 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*/ 1583 void 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*/ 1596 void 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*/ 1610 void 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*/ 1642 void 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*/ 1691 void 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*/ 1704 void 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*/ 1736 void 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*/ 1781 void 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*/ 1796 void 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.