Defining a custom element by basix

1

I want define a element, the basis functions is (x-1/3, y-1/3), functional is (div u, 1), This DOF is associated with face 0 of the reference element. My code is as follows

mport basix
import symfem
import numpy as np
from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace

#  spans:
#
#  .. math::
#     \left\{(x-1/3, y-1/3))\right\}.
#     \left\{(1, 0),\; (x, 0),\; (y, 0),\; (0, 1),\; (0, x),\; (0, y)\right\}.

#  The degree of free is the integral over the cell of the divergence of function over the geometry dimension.
# polynomial coefficients
wcoeffs = np.zeros((1, 6))

pts, wts = basix.make_quadrature(CellType.triangle, 2)
poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.triangle, 1, pts)
x = pts[:, 0]
y = pts[:, 1]
for i in range(3):
    wcoeffs[0, i] = sum((x-1/3) * poly[i, :] * wts)
    wcoeffs[0, 3 + i] = sum((y-1/3) * poly[i, :] * wts)
# # Interpolation
# # -------------
pts1, wts1 = basix.make_quadrature(CellType.triangle, 1)
# # The points associated with the cell are calculated by mapping the quadrature points to the cell.
print(pts1, wts1)
x = [[], [], [], []]
for _ in range(3):
    x[0].append(np.zeros((0, 2)))
    x[1].append(np.zeros((0, 2)))
x[2].append(pts1)
dim = 2
M = [[], [], [], []]
for _ in range(3):
    M[0].append(np.zeros((0, 2, 0, 1)))
    M[1].append(np.zeros((0, 2, 0, 1)))
mat = np.zeros((1, 2, len(wts1), 2))
mat[0, 0, :, 0] = np.empty((len(wts1)))
mat[0, 1, :, 0] = np.empty((len(wts1)))
mat[0, 0, :, 1] = wts1/dim
mat[0, 1, :, 1] = wts1/dim
M[2].append(mat)
element = basix.create_custom_element(
    CellType.triangle, [2], wcoeffs, x, M, 1, MapType.identity, SobolevSpace.L2,
    True, 0, 1)

And the error is
Traceback (most recent call last):
File “/home/yxchen/Documents/python workspace/code/DG0.py”, line 51, in
element = basix.create_custom_element(
RuntimeError: M has the wrong shape (dimension 3 is wrong)

2

Your M has the wrong number of derivatives (dimension 3): you would need this to be 3 as you have the function value, x-derivative and y-derivative.

If you’re hoping to use this element with DOLFINx, though, elements whose M uses derivatives are not yet supported. You might instead like to try this alternative element which is uses an integral over the cell to give the function you want:

import basix

import numpy as np
from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace, PolysetType

wcoeffs = np.zeros((1, 6))

pts, wts = basix.make_quadrature(CellType.triangle, 2)
poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.triangle, 1, pts)
x = pts[:, 0]
y = pts[:, 1]
for i in range(3):
    wcoeffs[0, i] = sum((x-1/3) * poly[i, :] * wts)
    wcoeffs[0, 3 + i] = sum((y-1/3) * poly[i, :] * wts)

pts1, wts1 = basix.make_quadrature(CellType.triangle, 2)
x = [[], [], [], []]
for _ in range(3):
    x[0].append(np.zeros((0, 2)))
    x[1].append(np.zeros((0, 2)))
x[2].append(pts1)

M = [[], [], [], []]
for _ in range(3):
    M[0].append(np.zeros((0, 2, 0, 1)))
    M[1].append(np.zeros((0, 2, 0, 1)))
mat = np.zeros((1, 2, len(wts1), 1))   
mat[0, 0, :, 0] = wts1 * 18 * (pts1[:, 0] - 1/3)
mat[0, 1, :, 0] = wts1 * 18 * (pts1[:, 1] - 1/3)
M[2].append(mat)
element = basix.create_custom_element(
    CellType.triangle, [2], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.L2,
    True, 0, 1, PolysetType.standard)

I used this test code to check that the basis function is (x-1/3, y-1/3):

points = np.array([[i / 10, j / 10] for i in range(11) for j in range(11 - i)])
values = element.tabulate(0, points)[0, :, 0]
for (x, y), (v0, v1) in zip(points, values):
    assert np.isclose(v0, x - 1/3)
    assert np.isclose(v1, y - 1/3)
1 Like