### How to find the quadrature points and weights of a mesh?

426

views

1

Hi,

Is there any way to find all the quadrature points and weights of a given function space? I tried the solution in this post, but it does not work for the current version of fenics.

Is there any way to find all the quadrature points and weights of a given function space? I tried the solution in this post, but it does not work for the current version of fenics.

`https://fenicsproject.org/qa/9633/quadrature-points?show=9633#q9633`

Community: FEniCS Project

### 1 Answer

2

Based on the post you mentioned and on https://fenicsproject.org/qa/9635/quadrature-element and https://fenicsproject.org/qa/3932/accessing-the-coordinates-of-a-degree-of-freedom, this should work in the current version

```
from fenics import *
mesh = UnitSquareMesh(10, 10)
gdim = mesh.geometry().dim()
quad_element = VectorElement(family = 'Quadrature',
cell = mesh.ufl_cell(),
degree = 2,
quad_scheme = 'default')
V = FunctionSpace(mesh, quad_element)
xq = V.tabulate_dof_coordinates().reshape((-1, gdim))
xq0 = xq[V.sub(0).dofmap().dofs()]
```

The code from the post works with

`from ffc.fiatinterface import create_quadrature`

written
11 months ago by
Adam Janecka

This will give the quadrature points and weights for a reference cell, how to get this relation for the global index?

written
11 months ago by
Mohan Wu

There surely must be a better solution but I was only able to come up with this ugly (hopefully correct) one:

```
from __future__ import print_function
from fenics import *
from ffc.fiatinterface import create_quadrature
mesh = UnitSquareMesh(2, 2)
gdim = mesh.geometry().dim()
shape = mesh.ufl_cell()
degree = 4
scheme = 'default'
quad_element = FiniteElement(family = 'Quadrature',
cell = shape,
degree = degree,
quad_scheme = scheme)
V = FunctionSpace(mesh, quad_element)
xq = V.tabulate_dof_coordinates().reshape((-1, gdim))
print('Quadrature points:\n', xq)
print('\n\t Point \t\t\t Weight')
for i in range(len(xq)):
u = Function(V)
u.vector()[i] = 1.0;
print(xq[i], ':\t', assemble(u*dx))
(points, weights) = create_quadrature(str(shape), degree, scheme)
print('\nReference element:')
print(' Points:\n', points)
print(' Weights:\n', weights)
```

written
11 months ago by
Adam Janecka

Please login to add an answer/comment or follow this question.

https://fenicsproject.org/qa/12309/quadrature-weights (which does not work either).