1. FEniCS Project Avatar FEniCS Project
  2. Core
  3. FFC
  4. Issues

problem with Quadrature element

Issue #138 closed
Florian Bruckner created an issue

When using the quadrature element the following error occurs since we updated to the new FEniCS 2016.2:

from dolfin import *

mesh = UnitCubeMesh(1,1,1)

V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

element = TensorElement("Quadrature", "tetrahedron", degree=2)
g = Expression((("1.0","0.0","0.0"),("0.0","2.0","0.0"),("0.0","0.0","1.0")), element=element)

a = inner(grad(u),g*grad(v))*dx
L = Constant(1.)*v*dx

assembler = SystemAssembler(a,L)

          Calling FFC just-in-time (JIT) compiler, this may take some time.
Quadrature element must have specified quadrature scheme (None) equal to the integral (default).
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-19-b906d302dfc7> in <module>()
----> 1 execfile("test.py")

/home/florian/test.py in <module>()
     13 L = Constant(1.)*v*dx
     14 
---> 15 assembler = SystemAssembler(a,L)
     16 

/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.pyc in __init__(self, A_form, b_form, bcs, form_compiler_parameters)
    436         """
    437         # Create dolfin Form objects referencing all data needed by assembler
--> 438         A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
    439         b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
    440 

/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.pyc in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
     65         return Form(form,
     66                     form_compiler_parameters=form_compiler_parameters,
---> 67                     function_spaces=function_spaces)
     68     else:
     69         raise TypeError("Invalid form type %s" % (type(form),))

/usr/lib/python2.7/dist-packages/dolfin/fem/form.pyc in __init__(self, form, form_compiler_parameters, function_spaces)
     87         # Jit the module, this will take some time...
     88         jit_result = jit(form, form_compiler_parameters,
---> 89                          mpi_comm=mesh.mpi_comm())
     90         if jit_result is None:
     91             cpp.dolfin_error("form.py",

/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.pyc in mpi_jit(*args, **kwargs)
     66         # Just call JIT compiler when running in serial
     67         if MPI.size(mpi_comm) == 1:
---> 68             return local_jit(*args, **kwargs)
     69 
     70         # Compile first on process 0

/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.pyc in jit(ufl_object, form_compiler_parameters, mpi_comm)
    131         cpp.dolfin_error("jit.py",
    132                          "perform just-in-time compilation of form",
--> 133                          "ffc.jit failed with message:\n%s" % (tb_text,))
    134 
    135     if isinstance(ufl_object, ufl.Form):

/usr/lib/python2.7/dist-packages/dolfin/cpp/common.pyc in dolfin_error(location, task, reason)
   2971 
   2972     """
-> 2973     return _common.dolfin_error(location, task, reason)
   2974 
   2975 def deprecation(feature, version_deprecated, message):

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to perform just-in-time compilation of form.
*** Reason:  ffc.jit failed with message:
Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 128, in jit
    result = ffc.jit(ufl_object, parameters=p)
  File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 198, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 120, in jit_build
    generate=jit_generate)
  File "/usr/lib/python2.7/dist-packages/dijitso/jit.py", line 160, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 141, in compile_form
    prefix, parameters, jit)
  File "/usr/lib/python2.7/dist-packages/ffc/compiler.py", line 183, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    for form in forms)
  File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    for form in forms)
  File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 172, in _analyze_form
    _attach_integral_metadata(form_data, parameters)
  File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 370, in _attach_integral_metadata
    _validate_quadrature_schemes_of_elements(quad_schemes, form_data.unique_sub_elements)
  File "/usr/lib/python2.7/dist-packages/ffc/analysis.py", line 384, in _validate_quadrature_schemes_of_elements
    error("Quadrature element must have specified quadrature scheme ((null)) equal to the integral (�O�)." % (qs, scheme))
  File "<string>", line 1, in <lambda>
  File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 171, in error
    raise self._exception_type(self._format_raw(*message))
Exception: Quadrature element must have specified quadrature scheme (None) equal to the integral (default).
.
*** Where:   This error was encountered inside jit.py.
*** Process: 0
*** 
*** DOLFIN version: 2016.2.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

the quad_scheme option of the TensorElement seems to be ignored!?

Comments (8)

  2. Jan Blechta

    Certainly not a FFC issue. Your code is incorrect as you don't set quadrature scheme of quadrature element. The correct code including workaround for https://bitbucket.org/fenics-project/ufl/issues/93 is:

    from dolfin import *

mesh = UnitCubeMesh(1,1,1)

V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

element = TensorElement("Quadrature", "tetrahedron", degree=2, quad_scheme='default')

# Workaround to UFL #93
element._quad_scheme = 'default'
for e in element.sub_elements():
    e._quad_scheme = 'default'

g = Expression((("1.0","0.0","0.0"),("0.0","2.0","0.0"),("0.0","0.0","1.0")), element=element)

a_metadata = {'quadrature_degree': element.degree(), 'quadrature_scheme': element.quadrature_scheme()}
a = inner(grad(u),g*grad(v))*dx(metadata=a_metadata)
L = Constant(1.)*v*dx

assembler = SystemAssembler(a,L)
  4. Florian Bruckner reporter

    Hello Jan, thanks for your fast response.

    I tried setting the quad_scheme before, but without your UFL workaround, but it seems to be ignored (also passing 'asdf' does not lead to an error message).

    Your testcode now works for me for a piecewise linear function space, but unfortunately we are using quadratic elements: 

    from dolfin import *

mesh = UnitCubeMesh(1,1,1)

V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)

element = TensorElement("Quadrature", "tetrahedron", degree=2, quad_scheme='default')

# Workaround to UFL #93
element._quad_scheme = 'default'
for e in element.sub_elements():
    e._quad_scheme = 'default'

g = Expression((("1.0","0.0","0.0"),("0.0","2.0","0.0"),("0.0","0.0","1.0")), element=element)

a = inner(grad(u),g*grad(v))*dx
L = Constant(1.)*v*dx

assembler = SystemAssembler(a,L)

    This code still fails with an Exception: Points must be equal to coordinates of quadrature points. I'm not exactly sure what this means, but it worked under fenics-1.6.

    greetings Florian

  5. Jan Blechta

    Yes. The code works by accident because estimated quadrature degree in a is 2. The correct general code would do:

    a_metadata = {'quadrature_degree': element.degree(), 'quadrature_scheme': element.quadrature_scheme()}
a = inner(grad(u),g*grad(v))*dx(metadata=a_metadata)

    This should work for any combination of coefficient/argument elements in the form a. I've corrected my post above.

  8. Martin Sandve Alnæs

    @blechta @florian_bruckner you can also do

    a = inner(grad(u),g*grad(v))*dx(degree=d, scheme=s)

    Behind the scenes this currently just creates the metadata.

  9. Log in to comment
Assignee
–
Type
bug
Priority
major
Status
closed
Component
jit
Milestone
2017.1
Version
2016.2
Votes
0
Watchers
2
Jira: the preferred issue tracker for Bitbucket. Join the team!