How to choose a specific quadrature rule?

General
1

I have a question about the codes as follows which is in Fenics tutorial of hyperelasticity

metadata={"quadrature_degree":4}
ds=ufl.Measure('ds',domain=domain,subdomain_data=facet_tag,metadata=metadata)
dx=ufl.Measure('dx',domain=domain,metadata=metadata)

I want to know how metadata is implemented in ufl.Measure?And what is the quadrature rule used here?(a 4th degree rule?)
Is there a way to change to different quadrature rules in Dolfinx?

In addition, I try to use the following codes to try to compare the degree of the qudrature rules but the result is the same for degree from 1 to 4, so it seems that I should also change the functionspace where the function f lives ?

from mpi4py import MPI
from dolfinx import fem,mesh
import ufl

domain=mesh.create_interval(MPI.COMM_WORLD,1,[0,1])
V=fem.FunctionSpace(domain,("CG",1))

#I want to intergate f=x^2 from 0 to 1
f=fem.Function(V)
u=lambda x: x[0]**2
f.interpolate(u)

metadata={"qudrature_degree":4} 
a=f*ufl.dx(domain=domain,metadata=metadata)
a_=fem.form(a)
b=fem.assemble_scalar(a_)

print(b)
2

For the first set of questions, I would recommend Total number of gauss point - #2 by dokken

Yes, to see an effect of the quadrature choice, you cannot just integrate a CG-1 function, as it is exact with degree=1. If you instead use

x=ufl.SpatialCoordinate(domain)
f = x[0]**2

you should see an effect with degree=1 and all the others.

3

Thanks a lot.Now I change the code to

x=ufl.SpatialCoordinate(domain)
f=x[0]**2
for k in rang(3):
    meta={"qudrature_degree":k}
    dx=ufl.Measure('dx',domain=domain,metadata=meta)
    a=f*ufl.dx
    b=fem.assemble_scalar(fem.form(a))
    print(b)

but I get the same results, is there something wrong?

4

Im not at a computer, but it might be super-convergence due to the mesh being structured. Either try with a higher order polynomial or an unstructured mesh.

5

Thanks a lot.I will have a try

6

For example:

from mpi4py import MPI
from dolfinx import fem, mesh
import ufl

domain = mesh.create_interval(MPI.COMM_WORLD, 1, [0, 1])

# I want to intergate f=x^2 from 0 to 1
x = ufl.SpatialCoordinate(domain)
f = x[0]**2
for i in range(5):
    metadata = {"quadrature_degree": i}
    a = f*ufl.dx(domain=domain, metadata=metadata)
    a_ = fem.form(a)
    b = fem.assemble_scalar(a_)
    print(i, f"{b:.5e}")

gives:

0 2.50000e-01
1 2.50000e-01
2 3.33333e-01
3 3.33333e-01
4 3.33333e-01

Note that you have a typo in

7

Great thanks!
Sorry for questions about this typo :smiling_face_with_tear:
when I use “qudrature_degree” it sends no error ,so what is the reason here?It seems that something is implicitly implemented here.

8

Why should it? DOLFINx only parses the metadata for known keys, if they are not found it uses default values. Potentially your dict could have thousands of entries, that we would not want to parse to check the validity of each key, or making educated guesses on what a key should be interpreted as.

9

Thanks.I am more clear about it.