Generate second degree Gaussian quadrature points for tets.
(mesh)
| 37 | for c in charges], axis=0); |
| 38 | |
| 39 | def get_quadrature_pts(mesh): |
| 40 | """ Generate second degree Gaussian quadrature points for tets. |
| 41 | """ |
| 42 | vertices = mesh.vertices; |
| 43 | tets = mesh.voxels; |
| 44 | |
| 45 | b = 0.1381966011250105; |
| 46 | a = 1.0 - 3 * b; |
| 47 | |
| 48 | q_pts = np.ndarray((mesh.num_voxels * 4, 3)); |
| 49 | q_pts[0:mesh.num_voxels*4:4] = \ |
| 50 | a*vertices[tets[:,0]] + \ |
| 51 | b*vertices[tets[:,1]] + \ |
| 52 | b*vertices[tets[:,2]] + \ |
| 53 | b*vertices[tets[:,3]]; |
| 54 | q_pts[1:mesh.num_voxels*4:4] = \ |
| 55 | b*vertices[tets[:,0]] + \ |
| 56 | a*vertices[tets[:,1]] + \ |
| 57 | b*vertices[tets[:,2]] + \ |
| 58 | b*vertices[tets[:,3]]; |
| 59 | q_pts[2:mesh.num_voxels*4:4] = \ |
| 60 | b*vertices[tets[:,0]] + \ |
| 61 | b*vertices[tets[:,1]] + \ |
| 62 | a*vertices[tets[:,2]] + \ |
| 63 | b*vertices[tets[:,3]]; |
| 64 | q_pts[3:mesh.num_voxels*4:4] = \ |
| 65 | b*vertices[tets[:,0]] + \ |
| 66 | b*vertices[tets[:,1]] + \ |
| 67 | b*vertices[tets[:,2]] + \ |
| 68 | a*vertices[tets[:,3]]; |
| 69 | return q_pts; |
| 70 | |
| 71 | def interpolate_at_quadrature_pts(mesh, sol_values, sol_grads): |
| 72 | tets = mesh.voxels; |