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