Question about quadratures rules

Question

Hello,

Is the Gauss quadrature the only quadrature rule in FEniCS?

I would like to use Newton-Cotes rule. These method gives diagonal mass matrix and avoids any other lumping method.

Many Thanks for answers,

Erwan

Comments

That's interesting. I can see this for 1D and linear hat functions as FEM basis. Do you have further insight you can share or references?

---

I found references on the following paper :

http://www.code-aster.org/V2/doc/v10/en/man_r/r3/r3.06.07.pdf

---

Thank you. However, as there are no concrete computations in there, I have digged a little further and shared my findings in an 'answer' below.

Answer 1

In development version of FFC vertex scheme is now available in addition to Gauss rules. Docstring in ffc/ffc/quadrature_scheme.py says

    # The vertex scheme, i.e., averaging the function value in the vertices
    # and multiplying with the simplex volume, is only of order 1 and
    # inferior to other generic schemes in terms of error reduction.
    # Equation systems generated with the vertex scheme have some
    # properties that other schemes lack, e.g., the mass matrix is
    # a simple diagonal matrix. This may be prescribed in certain cases.

This command

parameters['form_compiler']['quadrature_scheme'] = 'vertex'

should utilize vertex scheme globally in your script. Alternatively you apply vertex scheme on single form like this

m = u*v*dx(None, form_compiler_parameters={'quadrature_degree': 1,
                                           'quadrature_scheme': 'vertex'}))

or pass it to assembler

M = assemble(u*v*dx, 
             form_compiler_parameters={'quadrature_degree': 1,
                                       'quadrature_scheme': 'vertex'})

or to whatever accepts form_compiler_parameters dictionary.

You can also consider implementing new rule you need into ffc/ffc/quadrature_scheme.py. Probably it is only file that needs a change to implement a new rule now.

Comments

Yeah, I know. I just suggested "implement new rule into ffc/ffc/quadrature_scheme.py, the rest of the FEniCS will probably handle it". Now I was just checking the second part:)

---

This is what I understood

Thanks

Erwan

---

Jan,

The argument form_compiler_parameters is not being accepted.

When I try, for example,

dx(None, form_compiler_parameters={'quadrature_degree': 5})

I get

TypeError: __call__() got an unexpected keyword argument 'form_compiler_parameters'

Cheers!

---

Measures (and other UFL stuff) have been substantially redesigned. I believe now it is

dx(None, metadata={'quadrature_degree': 5})

(not tested).

---

Thanks Jan. It works.

Answer 2

This is more a comment than an answer, but I thought it may help you and hopefully others, when thinking about enforcing diagonal mass matrices.

As I understand now, in 2D, the Newton-Cotes formulae give diagonal mass matrices

• for $P_1$ and $P_2^+ $ elements
• at the expense of a consistency error of order $2$ and $4$, respectively,

see, e.g., this lecture by Kuzmin for the formulae (lecture 5) and the definition of the elements (lecture 7).

Since the $P_2$ elements are never zero at the center of gravity, see the symbolic computations in this python script, the higher order Newton-Cotes formulae do NOT give a diagonal mass matrix.

One may combine the $P_1$ and $P_2$ exact formuale, but this leads inevitably to an consistency error of order $2$ which is worse than the consistency error that comes with the $P_2$ approximation.