1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
|
# Copyright (c) 2020 Chris Richardson
# FEniCS Project
# SPDX-License-Identifier: MIT
import basix
import numpy as np
def test_gll_warped_pyramid():
n = 8
# Check that all the surface points of the pyramid match up with the same points on
# quad and triangle
tri_pts = basix.create_lattice(basix.CellType.triangle, n, basix.LatticeType.gll_warped, True)
quad_pts = basix.create_lattice(basix.CellType.quadrilateral, n, basix.LatticeType.gll_warped, True)
pyr_pts = basix.create_lattice(basix.CellType.pyramid, n, basix.LatticeType.gll_warped, True)
# Remove any near-zero values to make sorting robust
pyr_pts[np.where(abs(pyr_pts) < 1e-12)] = 0.0
tri_pts[np.where(abs(tri_pts) < 1e-12)] = 0.0
quad_pts[np.where(abs(quad_pts) < 1e-12)] = 0.0
idx = np.where(np.isclose(pyr_pts[:, 0], 0.0))
pyr_x0 = pyr_pts[idx][:, 1:]
assert np.allclose(np.sort(tri_pts), np.sort(pyr_x0))
idx = np.where(np.isclose(pyr_pts[:, 0] + pyr_pts[:, 2], 1.0))
pyr_xz = pyr_pts[idx][:, 1:]
assert np.allclose(np.sort(tri_pts), np.sort(pyr_xz))
idx = np.where(np.isclose(pyr_pts[:, 1], 0.0))
pyr_y0 = pyr_pts[idx][:, 0::2]
assert np.allclose(np.sort(tri_pts), np.sort(pyr_y0))
idx = np.where(np.isclose(pyr_pts[:, 1] + pyr_pts[:, 2], 1.0))
pyr_yz = pyr_pts[idx][:, 0::2]
assert np.allclose(np.sort(tri_pts), np.sort(pyr_yz))
idx = np.where(np.isclose(pyr_pts[:, 2], 0.0))
pyr_z0 = pyr_pts[idx][:, :2]
assert np.allclose(np.sort(quad_pts), np.sort(pyr_z0))
def test_gll_warped_tetrahedron():
n = 8
# Check that all the surface points of the tet match up with the same points on
# triangle
tri_pts = basix.create_lattice(basix.CellType.triangle, n, basix.LatticeType.gll_warped, True)
tet_pts = basix.create_lattice(basix.CellType.tetrahedron, n, basix.LatticeType.gll_warped, True)
tet_pts[np.where(abs(tet_pts) < 1e-12)] = 0.0
tri_pts[np.where(abs(tri_pts) < 1e-12)] = 0.0
idx = np.where(np.isclose(tet_pts[:, 0], 0.0))
tet_x0 = tet_pts[idx][:, 1:]
assert np.allclose(np.sort(tri_pts), np.sort(tet_x0))
idx = np.where(np.isclose(tet_pts[:, 1], 0.0))
tet_y0 = tet_pts[idx][:, 0::2]
assert np.allclose(np.sort(tri_pts), np.sort(tet_y0))
idx = np.where(np.isclose(tet_pts[:, 2], 0.0))
tet_z0 = tet_pts[idx][:, :2]
assert np.allclose(np.sort(tri_pts), np.sort(tet_z0))
# Project x+y+z=1 onto x=0
idx = np.where(np.isclose(tet_pts[:, 0] + tet_pts[:, 1] + tet_pts[:, 2], 1.0))
tet_xyz = tet_pts[idx][:, 1:]
assert np.allclose(np.sort(tri_pts), np.sort(tet_xyz))
|