From 20907a036c0f27244e043fcb68211b87890191b7 Mon Sep 17 00:00:00 2001 From: Johnny Vogels <35307256+jmv2009@users.noreply.github.com> Date: Wed, 11 Dec 2024 20:35:09 +0100 Subject: [PATCH] Laplacian of trace kernel for internal dofs TNT quadrilateral According to https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py, we need to multiply by the primary bubble function, and take the laplacian. To this end tabulation which also generates the derivatives is used. --- python/demo/demo_tnt-elements.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/python/demo/demo_tnt-elements.py b/python/demo/demo_tnt-elements.py index 06eadf560b6..071c29fe016 100644 --- a/python/demo/demo_tnt-elements.py +++ b/python/demo/demo_tnt-elements.py @@ -218,10 +218,18 @@ def create_tnt_quad(degree): x[2].append(np.zeros([0, 2])) M[2].append(np.zeros([0, 1, 0, 1])) else: - pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree - 1) - poly = basix.tabulate_polynomials( - basix.PolynomialType.legendre, basix.CellType.quadrilateral, degree - 2, pts + pts, wts = basix.make_quadrature(basix.CellType.quadrilateral, 2 * degree) + u = pts[:, 0] + v = pts[:, 1] + pol_set = basix.polynomials.tabulate_polynomial_set( + basix.CellType.quadrilateral, basix.PolynomialType.legendre, degree-2, 2, pts ) + # this assumes the conventional [0 to 1][0 to 1] domain of the reference element, + # and takes the Laplacian of (1-u)*u*(1-v)*v*pol_set[0], + # cf https://github.com/mscroggs/symfem/blob/main/symfem/elements/tnt.py + poly = (pol_set[5]+pol_set[3])*(u-1)*u*(v-1)*v+ \ + 2*(pol_set[2]*(u-1)*u*(2*v-1)+pol_set[1]*(v-1)*v*(2*u-1)+ \ + pol_set[0]*((u-1)*u+(v-1)*v)) face_ndofs = poly.shape[0] x[2].append(pts) mat = np.zeros((face_ndofs, 1, len(pts), 1)) @@ -300,7 +308,7 @@ def poisson_error(V: fem.FunctionSpace): tnt_degrees.append(2) tnt_ndofs.append(V.dofmap.index_map.size_global) tnt_errors.append(poisson_error(V)) -print(f"TNT degree 2 error: {tnt_errors[-1]}") +print(f"TNT degree 1 error: {tnt_errors[-1]}") for degree in range(2, 9): V = fem.functionspace(msh, create_tnt_quad(degree)) tnt_degrees.append(degree + 1)