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


426
views
1
11 months ago by
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.
https://fenicsproject.org/qa/9633/quadrature-points?show=9633#q9633​
 
Community: FEniCS Project

1 Answer


2
11 months ago by
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()]
​
Thanks for the reply. This code works for finding the quadrature points. Will there be any way to find the associate quadrature weight? Like in this post:
https://fenicsproject.org/qa/12309/quadrature-weights (which does not work either).
written 11 months ago by Mohan Wu  
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.

Similar posts:
Search »