Get quadrature points and weights

[IvanBioli]

Hi,

is there a way to get the xy coordinates of quadrature points, and the corresponding weights, employed in integrals?

To provide a MWE: to integrate the function x_0^2 on a square using a quadrature rule of degree 2 I can write:

from nutils import mesh
from nutils.expression_v2 import Namespace

nelems = 10
topo, x = mesh.unitsquare(nelems, etype='square')
ns = Namespace()
ns.x = x
ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))

I = topo.integral('x_0^2 dV' @ ns, degree=2)
I.eval()

Is there any way I can get the quadrature points and weights?

Best,
Ivan

[gertjanvanzwieten]

The command topo.integral('x_0^2 dV' @ ns, degree=2) is short for topo.sample('gauss', 2).integral('x_0^2 dV' @ ns); from there it's a small step to topo.sample('gauss', 2).eval('x_0' @ ns) to get a vector of point coordinates. The weights are a bit more involved because there is no ready made evaluable for it, so we will have to make one. I modified your MWE below to include such a construction.

from nutils import mesh, function
from nutils.expression_v2 import Namespace

nelems = 10
topo, x = mesh.unitsquare(nelems, etype='square')
ns = Namespace()
ns.x = x
ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))

gauss = topo.sample('gauss', 2)
I = gauss.integral('x_0^2 dV' @ ns)

class Weights(function.Array):
    def __init__(self, sample):
        self.sample = sample
        super().__init__(shape=(), dtype=float, spaces=sample.spaces, arguments={})
    def lower(self, lowerargs):
        _, index = lowerargs.transform_chains[self.sample.space]
        weights = self.sample.get_evaluable_weights(index)
        if lowerargs.points_shape != weights.shape:
            raise ValueError('the Weights evaluable can only be evaluated on the sample for which it was created')
        return weights
        
v = gauss.eval('x_0^2 dV' @ ns)
w = gauss.eval(Weights(gauss))

assert v @ w == I.eval()

[gertjanvanzwieten]

That is strange, for me the shape is (400,) which is what I would expect. Are you running the code precisely as above or did you make modifications?

[MichaelLoibl]

Yes, sorry. I am not used to Python syntax. I meant (400,) when I said 400x1. But still I don't understand where I see the x- and where the y-coordinates?

[joostvanzwieten]

With gauss.eval(ns.x) you can evaluate the geometry in the quadrature points.

[gertjanvanzwieten]

With gauss.eval(ns.x) you can evaluate the geometry in the quadrature points.

Or, equivalently, gauss.eval('x_i' @ ns). The reason you do not get the y-coordinate is that you evaluate x_0, so only the first component.

[MichaelLoibl]

Thank you a lot. That was exactly what I was look for :)