1import numpy as np 2 3 4def buildmats(qref, qweight, mat_dtype="float64"): 5 P, Q, dim = 6, 4, 2 6 interp = np.empty(P * Q, dtype=mat_dtype) 7 grad = np.empty(dim * P * Q, dtype=mat_dtype) 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 51 52def buildmatshdiv(qref, qweight, mat_dtype="float64"): 53 P, Q, dim = 4, 4, 2 54 interp = np.empty(dim * P * Q, dtype=mat_dtype) 55 div = np.empty(P * Q, dtype=mat_dtype) 56 57 qref[0] = -1. / np.sqrt(3.) 58 qref[1] = qref[0] 59 qref[2] = qref[0] 60 qref[3] = -qref[0] 61 qref[4] = -qref[0] 62 qref[5] = -qref[0] 63 qref[6] = qref[0] 64 qref[7] = qref[0] 65 qweight[0] = 1. 66 qweight[1] = 1. 67 qweight[2] = 1. 68 qweight[3] = 1. 69 70 # Loop over quadrature points 71 for i in range(Q): 72 x1 = qref[0 * Q + i] 73 x2 = qref[1 * Q + i] 74 # Interp 75 interp[(i + 0) * P + 0] = 0. 76 interp[(i + Q) * P + 0] = 1. - x2 77 interp[(i + 0) * P + 1] = x1 - 1. 78 interp[(i + Q) * P + 1] = 0. 79 interp[(i + 0) * P + 2] = -x1 80 interp[(i + Q) * P + 2] = 0. 81 interp[(i + 0) * P + 3] = 0. 82 interp[(i + Q) * P + 3] = x2 83 # Div 84 div[i * P + 0] = -1. 85 div[i * P + 1] = 1. 86 div[i * P + 2] = -1. 87 div[i * P + 3] = 1. 88 89 return interp, div 90 91 92def buildmatshcurl(qref, qweight, mat_dtype="float64"): 93 P, Q, dim = 3, 4, 2 94 interp = np.empty(dim * P * Q, dtype=mat_dtype) 95 curl = np.empty(P * Q, dtype=mat_dtype) 96 97 qref[0] = 0.2 98 qref[1] = 0.6 99 qref[2] = 1. / 3. 100 qref[3] = 0.2 101 qref[4] = 0.2 102 qref[5] = 0.2 103 qref[6] = 1. / 3. 104 qref[7] = 0.6 105 qweight[0] = 25. / 96. 106 qweight[1] = 25. / 96. 107 qweight[2] = -27. / 96. 108 qweight[3] = 25. / 96. 109 110 # Loop over quadrature points 111 for i in range(Q): 112 x1 = qref[0 * Q + i] 113 x2 = qref[1 * Q + i] 114 # Interp 115 interp[(i + 0) * P + 0] = -x2 116 interp[(i + Q) * P + 0] = x1 117 interp[(i + 0) * P + 1] = x2 118 interp[(i + Q) * P + 1] = 1. - x1 119 interp[(i + 0) * P + 2] = 1. - x2 120 interp[(i + Q) * P + 2] = x1 121 # Curl 122 curl[i * P + 0] = 2. 123 curl[i * P + 1] = -2. 124 curl[i * P + 2] = -2. 125 126 return interp, curl 127