1import numpy as np 2 3 4def buildmats(qref, qweight): 5 P, Q, dim = 6, 4, 2 6 interp = np.empty(P * Q, dtype="float64") 7 grad = np.empty(dim * P * Q, dtype="float64") 8 9 qref[0] = 0.2 10 qref[1] = 0.6 11 qref[2] = 1. / 3. 12 qref[3] = 0.2 13 qref[4] = 0.2 14 qref[5] = 0.2 15 qref[6] = 1. / 3. 16 qref[7] = 0.6 17 qweight[0] = 25. / 96. 18 qweight[1] = 25. / 96. 19 qweight[2] = -27. / 96. 20 qweight[3] = 25. / 96. 21 22 # Loop over quadrature points 23 for i in range(Q): 24 x1 = qref[0 * Q + i] 25 x2 = qref[1 * Q + i] 26 # Interp 27 interp[i * P + 0] = 2. * (x1 + x2 - 1.) * (x1 + x2 - 1. / 2.) 28 interp[i * P + 1] = -4. * x1 * (x1 + x2 - 1.) 29 interp[i * P + 2] = 2. * x1 * (x1 - 1. / 2.) 30 interp[i * P + 3] = -4. * x2 * (x1 + x2 - 1.) 31 interp[i * P + 4] = 4. * x1 * x2 32 interp[i * P + 5] = 2. * x2 * (x2 - 1. / 2.) 33 # Grad 34 grad[(i + 0) * P + 0] = 2. * \ 35 (1. * (x1 + x2 - 1. / 2.) + (x1 + x2 - 1.) * 1.) 36 grad[(i + Q) * P + 0] = 2. * \ 37 (1. * (x1 + x2 - 1. / 2.) + (x1 + x2 - 1.) * 1.) 38 grad[(i + 0) * P + 1] = -4. * (1. * (x1 + x2 - 1.) + x1 * 1.) 39 grad[(i + Q) * P + 1] = -4. * (x1 * 1.) 40 grad[(i + 0) * P + 2] = 2. * (1. * (x1 - 1. / 2.) + x1 * 1.) 41 grad[(i + Q) * P + 2] = 2. * 0. 42 grad[(i + 0) * P + 3] = -4. * (x2 * 1.) 43 grad[(i + Q) * P + 3] = -4. * (1. * (x1 + x2 - 1.) + x2 * 1.) 44 grad[(i + 0) * P + 4] = 4. * (1. * x2) 45 grad[(i + Q) * P + 4] = 4. * (x1 * 1.) 46 grad[(i + 0) * P + 5] = 2. * 0. 47 grad[(i + Q) * P + 5] = 2. * (1. * (x2 - 1. / 2.) + x2 * 1.) 48 49 return interp, grad 50