Cannot align integratoin points and weights in a mesh cell

General
1

With dolfin and ffc, I found the post to generate physical quadrature points in the following post.

quadrature points

Since there is no weights, I use create_quadrature trying to get the quadrature weights. These two block of codes are run as follows:

# fristly get the quad points
import ffc
from dolfin import *

mesh = UnitSquareMesh(1, 1)

el = FiniteElement("Quadrature", mesh.ufl_cell(),  degree=1, quad_scheme="default")
V = FunctionSpace(mesh, el)
print(V.tabulate_dof_coordinates())

points, weights = ffc.fiatinterface.create_quadrature(mesh.cell_name(), degree=1, scheme="default")
print(points, weights)

Get the result:

Calling FFC just-in-time (JIT) compiler, this may take some time.
[[0.66666667 0.33333333]
[0.33333333 0.66666667]]
[[0.33333333 0.33333333]] [0.5]

With the same mesh cell, they generate different number of quadrature points, let alone match them precisely. How can I get the quadrature points and relevant weights of a physical mesh cell ?

2

This quadrature is based on the reference element.

The quadrature extracted with:

Gives you the quadrature points in physical space, i.e.

# fristly get the quad points
import ffc
from dolfin import *

mesh = UnitSquareMesh(1, 1)

el = FiniteElement("Quadrature", mesh.ufl_cell(),  degree=2, quad_scheme="default")
V = FunctionSpace(mesh, el)
x = V.tabulate_dof_coordinates()


for cell in cells(mesh):
    print(V.dofmap().cell_dofs(cell.index()), x[V.dofmap().cell_dofs(cell.index())])

gives you

[1 0 2] [[0.33333333 0.16666667]
 [0.83333333 0.66666667]
 [0.83333333 0.16666667]]
[4 3 5] [[0.16666667 0.33333333]
 [0.66666667 0.83333333]
 [0.16666667 0.83333333]]

which are the three quadrature points in each cell.

3

Thanks. I see that the mesh function gives two triangular elements for the square. With the following code

# fristly get the quad points
import ffc
from dolfin import *

degree=4

mesh = UnitSquareMesh(1, 1)

el = FiniteElement("Quadrature", mesh.ufl_cell(),  degree=degree, quad_scheme="default")
V = FunctionSpace(mesh, el)
x = V.tabulate_dof_coordinates()

print("\n physical cel")

for cell in cells(mesh):
    print("cell",cell.index(),"cell nodes",V.dofmap().cell_dofs(cell.index()),'\n', x[V.dofmap().cell_dofs(cell.index())],'\n')

print("\nreference cel, cell name is  " ,mesh.cell_name())
points, weights = ffc.fiatinterface.create_quadrature(mesh.cell_name(), degree=degree, scheme="default")
print("weights of ref cel: \n", weights)

I get the results:

 physical cel
cell 0 cell nodes [1 0 2 3 4 5] 
 [[0.90842379 0.09157621]
 [0.90842379 0.81684757]
 [0.18315243 0.09157621]
 [0.55405151 0.44594849]
 [0.55405151 0.10810302]
 [0.89189698 0.44594849]] 

cell 1 cell nodes [ 7  6  8  9 10 11] 
 [[0.09157621 0.90842379]
 [0.81684757 0.90842379]
 [0.09157621 0.18315243]
 [0.44594849 0.55405151]
 [0.10810302 0.55405151]
 [0.44594849 0.89189698]] 


reference cel, cell name is   triangle
weights of ref cel: 
 [0.05497587 0.05497587 0.05497587 0.11169079 0.11169079 0.11169079]

Does it mean that I can multiply the weights form the ref cel accordingly with the quad points above? Do the quad points from the physical cell automatically align with the weights from the ref cell?

4

The sum of weights should equal to the volume of the physical cell. If I just use the quad weights from the ref cell, it surely won’t get correct integration because the volume of general physical cell != volume of the ref cell.

5

Now you are mixing quadrature weights and the absolute value of the determinant of the Jacobian of the mapping from the physical element to the reference element.

That determinant can be found as shown here: Spatial coordinate local values - #4 by dokken

6

Thanks, got it. So the order of weights in ref cell will exactly match with the order of quad points in physical cell, right?

7

Yes, they should match.

1 Like