Hi all,
I am trying to use Basix for building a custom anisotropic finite element defined, in a rectangle, as product of order 2 elements in the x variable and order 1 elements in the y variable. Specifically, this element will be defined by evaluating Lagrange polynomials on the following points:
o---o---o
| |
| |
o---o---o
We have the following set of polynomials spanned by each element: { 1, y, x, xy, x^2, x^2y }.
For that, I followed the steps in Creating a custom element — Basix 0.10.0.0 documentation.
First, we define the following 6x9 matrix to store the coefficients
import numpy as np
wcoeffs = np.zeros((6,9))
wcoeffs[0,0] = 1
wcoeffs[1,1] = 1
wcoeffs[2,3] = 1
wcoeffs[3,4] = 1
wcoeffs[4,6] = 1
wcoeffs[5,7] = 1
The Ciarlet functionals are, in this case, evaluations at the 4 corners and 2 edges of the quadrilateral.
#
# Nodes:
#
x = [ [], [], [], [] ]
# Corners
x[0].append(np.array([[0.0, 0.0]]))
x[0].append(np.array([[1.0, 0.0]]))
x[0].append(np.array([[0.0, 1.0]]))
x[0].append(np.array([[1.0, 1.0]]))
# Edges
x[1].append(np.array([[0.5, 0.0]]))
x[1].append(np.array([[0.5, 1.0]]))
x[1].append(np.zeros((0,2)))
x[1].append(np.zeros((0,2)))
# Empty array of points inside the 2d cell
x[2].append(np.zeros((0,2)))
As all the DOFs are point evaluations the matrices are all identity matrices for the entities that have a point:
M = [ [], [], [], [] ]
# 4 corners
for _ in range(4):
M[0].append(np.array([[[[1.0]]]]))
# 2 edges
M[1].append(np.array([[[[1.0]]]]))
M[1].append(np.array([[[[1.0]]]]))
M[1].append(np.zeros((0, 1, 0, 1)))
M[1].append(np.zeros((0, 1, 0, 1)))
# Inside cell
M[2].append(np.zeros((0, 1, 0, 1)))
Now we create the element
import dolfinx
import basix
from basix import CellType, MapType, SobolevSpace, PolysetType
import basix.ufl
P2xP1 = basix.ufl.custom_element(
cell_type=CellType.quadrilateral,
reference_value_shape=[],
wcoeffs=wcoeffs,
x=x,
M=M,
interpolation_nderivs=0,
map_type=MapType.identity,
sobolev_space=SobolevSpace.H1,
discontinuous=False,
embedded_subdegree=1,
embedded_superdegree=2,
polyset_type=PolysetType.standard,
)
And compute the matrix for a Poisson problem in a single cell (a 1x1 mesh)
from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
from ufl import dx, grad, inner
msh = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)),
n=(1,1),
cell_type=mesh.CellType.quadrilateral,
)
V = fem.functionspace(msh, P2xP1)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) * dx
# Matrix
import scipy
A=dolfinx.fem.assemble_matrix(dolfinx.fem.form(a))
As = scipy.sparse.csr_matrix((A.data, A.indices, A.indptr))
def my_print(A, float_formatter = "{:.1f}".format):
np.set_printoptions(formatter={'float_kind':float_formatter})
print(A)
A_arr = As.toarray()
print("Matrix shape:", A_arr.shape)
my_print(A_arr)
According to the number of degrees of freedom for this the element, I expected a 6x6 matrix. But, unexpectedly, I obtain a 8x6 matrix, where the two extra rows are composed by zeros. In particular, this the result of running the minimal script defined above:
Matrix shape: (8, 6)
[[0.9 0.1 0.3 0.1 -0.8 -0.5]
[0.1 0.9 0.1 0.3 -0.8 -0.5]
[0.3 0.1 0.9 0.1 -0.5 -0.8]
[0.1 0.3 0.1 0.9 -0.5 -0.8]
[-0.8 -0.8 -0.5 -0.5 2.3 0.4]
[-0.5 -0.5 -0.8 -0.8 0.4 2.3]
[0.0 0.0 0.0 0.0 0.0 0.0]
[0.0 0.0 0.0 0.0 0.0 0.0]]
Any help would be appreciated!
Thanks!