10ef72598Sjeremyltimport numpy as np 20ef72598Sjeremylt 30ef72598Sjeremylt 4*80a9ef05SNatalie Beamsdef buildmats(qref, qweight, mat_dtype="float64"): 50ef72598Sjeremylt P, Q, dim = 6, 4, 2 6*80a9ef05SNatalie Beams interp = np.empty(P * Q, dtype=mat_dtype) 7*80a9ef05SNatalie Beams grad = np.empty(dim * P * Q, dtype=mat_dtype) 80ef72598Sjeremylt 90ef72598Sjeremylt qref[0] = 0.2 100ef72598Sjeremylt qref[1] = 0.6 110ef72598Sjeremylt qref[2] = 1. / 3. 120ef72598Sjeremylt qref[3] = 0.2 130ef72598Sjeremylt qref[4] = 0.2 140ef72598Sjeremylt qref[5] = 0.2 150ef72598Sjeremylt qref[6] = 1. / 3. 160ef72598Sjeremylt qref[7] = 0.6 170ef72598Sjeremylt qweight[0] = 25. / 96. 180ef72598Sjeremylt qweight[1] = 25. / 96. 190ef72598Sjeremylt qweight[2] = -27. / 96. 200ef72598Sjeremylt qweight[3] = 25. / 96. 210ef72598Sjeremylt 220ef72598Sjeremylt # Loop over quadrature points 230ef72598Sjeremylt for i in range(Q): 240ef72598Sjeremylt x1 = qref[0 * Q + i] 250ef72598Sjeremylt x2 = qref[1 * Q + i] 260ef72598Sjeremylt # Interp 270ef72598Sjeremylt interp[i * P + 0] = 2. * (x1 + x2 - 1.) * (x1 + x2 - 1. / 2.) 280ef72598Sjeremylt interp[i * P + 1] = -4. * x1 * (x1 + x2 - 1.) 290ef72598Sjeremylt interp[i * P + 2] = 2. * x1 * (x1 - 1. / 2.) 300ef72598Sjeremylt interp[i * P + 3] = -4. * x2 * (x1 + x2 - 1.) 310ef72598Sjeremylt interp[i * P + 4] = 4. * x1 * x2 320ef72598Sjeremylt interp[i * P + 5] = 2. * x2 * (x2 - 1. / 2.) 330ef72598Sjeremylt # Grad 340ef72598Sjeremylt grad[(i + 0) * P + 0] = 2. * \ 350ef72598Sjeremylt (1. * (x1 + x2 - 1. / 2.) + (x1 + x2 - 1.) * 1.) 360ef72598Sjeremylt grad[(i + Q) * P + 0] = 2. * \ 370ef72598Sjeremylt (1. * (x1 + x2 - 1. / 2.) + (x1 + x2 - 1.) * 1.) 380ef72598Sjeremylt grad[(i + 0) * P + 1] = -4. * (1. * (x1 + x2 - 1.) + x1 * 1.) 390ef72598Sjeremylt grad[(i + Q) * P + 1] = -4. * (x1 * 1.) 400ef72598Sjeremylt grad[(i + 0) * P + 2] = 2. * (1. * (x1 - 1. / 2.) + x1 * 1.) 410ef72598Sjeremylt grad[(i + Q) * P + 2] = 2. * 0. 420ef72598Sjeremylt grad[(i + 0) * P + 3] = -4. * (x2 * 1.) 430ef72598Sjeremylt grad[(i + Q) * P + 3] = -4. * (1. * (x1 + x2 - 1.) + x2 * 1.) 440ef72598Sjeremylt grad[(i + 0) * P + 4] = 4. * (1. * x2) 450ef72598Sjeremylt grad[(i + Q) * P + 4] = 4. * (x1 * 1.) 460ef72598Sjeremylt grad[(i + 0) * P + 5] = 2. * 0. 470ef72598Sjeremylt grad[(i + Q) * P + 5] = 2. * (1. * (x2 - 1. / 2.) + x2 * 1.) 480ef72598Sjeremylt 490ef72598Sjeremylt return interp, grad 50