1099 double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
1100 if (edgeN != NULL) {
1106 double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
1107 *diff_edgeN30 = NULL;
1108 if (diff_edgeN != NULL) {
1109 diff_edgeN01 = diff_edgeN[0];
1110 diff_edgeN12 = diff_edgeN[1];
1111 diff_edgeN23 = diff_edgeN[2];
1112 diff_edgeN30 = diff_edgeN[3];
1116 for (; ee < 4; ee++)
1125 for (; ii < GDIM; ii++) {
1126 int node_shift = ii * 4;
1127 int node_diff_shift = 2 * node_shift;
1129 double shape0 =
N[node_shift + n0];
1130 double shape1 =
N[node_shift + n1];
1131 double shape2 =
N[node_shift + n2];
1132 double shape3 =
N[node_shift + n3];
1134 double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
1135 double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
1136 double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
1137 double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
1139 double extrude_zeta01 = shape0 + shape1;
1140 double extrude_ksi12 = shape1 + shape2;
1141 double extrude_zeta23 = shape2 + shape3;
1142 double extrude_ksi30 = shape0 + shape3;
1144 double bubble_ksi = extrude_ksi12 * extrude_ksi30;
1145 double bubble_zeta = extrude_zeta01 * extrude_zeta23;
1147 double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
1148 double diff_extrude_zeta01[2];
1149 double diff_extrude_ksi12[2];
1150 double diff_extrude_zeta23[2];
1151 double diff_extrude_ksi30[2];
1152 double diff_bubble_ksi[2];
1153 double diff_bubble_zeta[2];
1156 for (;
d < 2;
d++) {
1157 double diff_shape0 = diffN[node_diff_shift + 2 * n0 +
d];
1158 double diff_shape1 = diffN[node_diff_shift + 2 * n1 +
d];
1159 double diff_shape2 = diffN[node_diff_shift + 2 * n2 +
d];
1160 double diff_shape3 = diffN[node_diff_shift + 2 * n3 +
d];
1162 (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1164 (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1166 (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1168 (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1169 diff_extrude_zeta01[
d] = diff_shape0 + diff_shape1;
1170 diff_extrude_ksi12[
d] = diff_shape1 + diff_shape2;
1171 diff_extrude_zeta23[
d] = diff_shape2 + diff_shape3;
1172 diff_extrude_ksi30[
d] = diff_shape0 + diff_shape3;
1173 diff_bubble_ksi[
d] = diff_extrude_ksi12[
d] * extrude_ksi30 +
1174 extrude_ksi12 * diff_extrude_ksi30[
d];
1175 diff_bubble_zeta[
d] = diff_extrude_zeta01[
d] * extrude_zeta23 +
1176 extrude_zeta01 * diff_extrude_zeta23[
d];
1179 double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1180 double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1181 diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1182 ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1184 ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1186 ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1188 ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1192 if (edgeN != NULL) {
1194 shift = ii * (
P[0]);
1195 cblas_dcopy(
P[0], L01, 1, &edgeN01[shift], 1);
1196 cblas_dscal(
P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1198 shift = ii * (
P[1]);
1199 cblas_dcopy(
P[1], L12, 1, &edgeN12[shift], 1);
1200 cblas_dscal(
P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1202 shift = ii * (
P[2]);
1203 cblas_dcopy(
P[2], L23, 1, &edgeN23[shift], 1);
1204 cblas_dscal(
P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1206 shift = ii * (
P[3]);
1207 cblas_dcopy(
P[3], L30, 1, &edgeN30[shift], 1);
1208 cblas_dscal(
P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1210 if (diff_edgeN != NULL) {
1213 shift = ii * (
P[0]);
1214 bzero(&diff_edgeN01[2 * shift],
sizeof(
double) * 2 * (
P[0]));
1216 for (;
d != 2; ++
d) {
1217 cblas_daxpy(
P[0], bubble_ksi * extrude_zeta01,
1218 &diffL01[
d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift +
d],
1221 diff_bubble_ksi[
d] * extrude_zeta01 +
1222 bubble_ksi * diff_extrude_zeta01[
d],
1223 L01, 1, &diff_edgeN01[2 * shift +
d], 2);
1228 shift = ii * (
P[1]);
1229 bzero(&diff_edgeN12[2 * shift],
sizeof(
double) * 2 * (
P[1]));
1231 for (;
d != 2; ++
d) {
1232 cblas_daxpy(
P[1], bubble_zeta * extrude_ksi12,
1233 &diffL12[
d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift +
d],
1236 diff_bubble_zeta[
d] * extrude_ksi12 +
1237 bubble_zeta * diff_extrude_ksi12[
d],
1238 L12, 1, &diff_edgeN12[2 * shift +
d], 2);
1243 shift = ii * (
P[2]);
1244 bzero(&diff_edgeN23[2 * shift],
sizeof(
double) * 2 * (
P[2]));
1246 for (;
d != 2; ++
d) {
1247 cblas_daxpy(
P[2], bubble_ksi * extrude_zeta23,
1248 &diffL23[
d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift +
d],
1251 diff_bubble_ksi[
d] * extrude_zeta23 +
1252 bubble_ksi * diff_extrude_zeta23[
d],
1253 L23, 1, &diff_edgeN23[2 * shift +
d], 2);
1258 shift = ii * (
P[3]);
1259 bzero(&diff_edgeN30[2 * shift],
sizeof(
double) * 2 * (
P[3]));
1261 for (;
d != 2; ++
d) {
1262 cblas_daxpy(
P[3], bubble_zeta * extrude_ksi30,
1263 &diffL30[
d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift +
d],
1266 diff_bubble_zeta[
d] * extrude_ksi30 +
1267 bubble_zeta * diff_extrude_ksi30[
d],
1268 L30, 1, &diff_edgeN30[2 * shift +
d], 2);