Per element quadrature

[tsokar]

Hi !
I have two questions concerning the definition of per-element quadratures...
In Element.gbasis it is stated that X can be shared by all elements (array((Ndim, Npoints))) or defined per element (array((Ndim, Nelems, Npoints))).
I tried to play with that a little bit: it works for the evaluation of the basis functions, but eventually fails when calculating the inverse of the jacobian. I must admit that the error is quite deep in the code, so that I cannot fully understand where the problem is :-/

Here is a MWE:

import numpy as np
import skfem as skf

mesh = skf.MeshQuad().refined(2)
e = skf.ElementQuad1()

basis = skf.Basis(mesh, e)
baseQuadratureX, baseQuadratureWeights = basis.quadrature
perElementQuadratureWeights = np.repeat(baseQuadratureWeights, mesh.nelements).reshape((mesh.nelements, baseQuadratureWeights.shape[0]))
perElementQuadratureX = np.repeat(baseQuadratureX[:,:,np.newaxis], mesh.nelements, axis=2).transpose((0,2,1))

basisPerElementQuadrature = skf.Basis(mesh, e, quadrature=(perElementQuadratureX, perElementQuadratureWeights))

My second question is the following: ultimately, I would like to consider a different amount of quadrature points in every element. Is it possible ? It would work using padding in the arrays defining the quadratures, but it will not be very efficient.

Thanks !
T.