FunctionSpace(mesh, "Quadrature", 1) broken in v 2016.1

[Umberto Villa]

In the python interface, the creation of a "Quadrature" FunctionSpace using the backward compatible (Fenics 1.6 or earlier) constructor is broken.

The following code fragment

import dolfin as dl                                                                                                             
mesh = dl.UnitSquareMesh(1,1)                                                                                                   
Vh = dl.FunctionSpace(mesh, "Quadrature", 1)

gives this error:

Missing quad_scheme in quadrature element.                                                                                          
Traceback (most recent call last):                                                                                                  
  File "<stdin>", line 1, in <module>                                                                                               
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 206, in __init__                    
    self._init_convenience(*args, **kwargs)                                                                                         
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 258, in _init_convenience constrained_domain=constrained_domain)                                                                              
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 224, in _init_from_ufl              
    dolfin_element, dolfin_dofmap = _compile_dolfin_element(element, mesh, constrained_domain=constrained_domain)                   
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 89, in _compile_dolfin_element      
    ufc_element, ufc_dofmap = jit(element, mpi_comm=mesh.mpi_comm())                                                                
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 65, in mpi_jit                           
    return local_jit(*args, **kwargs)                                                                                               
  File "/home/fenics/build/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 124, in jit                              
    result = ffc.jit(ufl_object, parameters=p)                                                                                      
  File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 201, in jit                                        
    module = jit_build_with_instant(ufl_object, module_name, parameters)                                                            
  File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 98, in jit_build_with_instant                      
    code_h, code_c = jit_generate(ufl_object, module_name, parameters)                                                              
  File "/home/fenics/build/lib/python2.7/site-packages/ffc/jitcompiler.py", line 65, in jit_generate                                
    parameters=parameters, jit=True)                                                                                                
  File "/home/fenics/build/lib/python2.7/site-packages/ffc/compiler.py", line 209, in compile_element                               
    analysis = analyze_elements(elements, parameters)                                                                               
  File "/home/fenics/build/lib/python2.7/site-packages/ffc/analysis.py", line 99, in analyze_elements                               
    error("Missing quad_scheme in quadrature element.")                                                                             
  File "<string>", line 1, in <lambda>                                                                                              
  File "/home/fenics/build/lib/python2.7/site-packages/ufl/log.py", line 158, in error                                              
    raise self._exception_type(self._format_raw(*message))                                                                          
Exception: Missing quad_scheme in quadrature element. 

Thanks in advance.

Umberto

[Martin Sandve Alnæs]

I think this is the result of a bugfix where compiling a form had side effects by setting the quad_scheme of the function to match the form. Try setting the missing quad_scheme of the element.

[Jan Blechta]

I thought that it was intentional to require quad scheme/degree with Quadrature element.

[Martin Sandve Alnæs]

Yes, because it allowed fixing the abovementioned bug properly, and otherwise the element would not be fully defined on its own.

[Umberto Villa]

Hi Martin and Jan,

Thank you for looking into this issue.

If I understand correctly the 1.6-like construction of a FunctionSpace (i.e. mesh, "Family", degree) does not allow to specify the quadrature rule, but only the degree.

If the plan for the future is to maintain this interface, it may be a good idea to introduce more than one family for Quadrature finite element. Something like "GaussQuadrature", "GaussLobattoQuadrature", etc. Another option is to simply remove "Quadrature" from the list of finite element families and then force the user to adopt the new interface when working with Quadrature Finite Element Spaces.

Best,

Umberto

[Jan Blechta]

Yes. The name Quadrature is misleading. It means "point evaluation (without FE interpolation)". Effective quadrature scheme is naturally associated with an employed measure.

[Jan Blechta]

So the working code is

from dolfin import *
mesh = UnitSquareMesh(1, 1)
element = FiniteElement("Quadrature", mesh.ufl_cell(), 1, quad_scheme="default")
V = FunctionSpace(mesh, element)

Is it really necessary, @martinal?

[Martin Sandve Alnæs]

It doesn't make any sense to create a FunctionSpace for a 'quadrature element' if the scheme and degree is not determined. The previous behaviour had silent bugs that are now fixed. Feel free to make the interface prettier, but the scheme and degree must be part of the element/functionspace.

[Jan Blechta]

@martinal, thanks, it makes sense to me now.