Crashes when projecting to quadrature space of degree 2.

Issue #1084 new
Grant Bruer created an issue

I’m trying to project a function to a quadrature space, but it crashes if I try to use a quadrature degree higher than 1.

from fenics import *

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'P', 1)

u = Function(V)

# Create quadrature space.
quad_params = {
    'quadrature_degree': 2,        # 1 works here, but not 2 or higher.
    'quadrature_rule': 'default',
    'representation': 'quadrature' # Also tried uflacs here as warning suggested.
}
quad_el = TensorElement('Quadrature', mesh.ufl_cell(), shape=(),
    degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])
Q = FunctionSpace(mesh, quad_el)

qdx = dx(metadata=quad_params)

# Set up a system for projecting u onto Q.
v = TestFunction(Q)
q = TrialFunction(Q)

a = q * v * qdx
L = u * v * qdx

# Solve the system.
q_ = Function(Q)
solve(a == L, q_)

Output using ‘quadrature’ representation:

/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturerepresentation.py:58: QuadratureRepresentationDeprecationWarning: 
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** representation instead.                               ***
*** ===================================================== ***
  issue_deprecation_warning()
Traceback (most recent call last):
  File "mwe.py", line 29, in <module>
    solve(a == L, q_)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in <listcomp>
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 460, in _compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 100, in compute_integral_ir
    ir["optimise_parameters"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/optimisedquadraturetransformer.py", line 56, in __init__
    QuadratureTransformerBase.__init__(self, *args)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 148, in __init__
    self.entity_type)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadratureutils.py", line 49, in create_psi_tables
    name_map, unique_tables = unique_psi_tables(flat_tables, eliminate_zeros)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadratureutils.py", line 179, in unique_psi_tables
    if abs(vals[r][c]) < format_epsilon:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Output using ‘uflacs’ representation:

Traceback (most recent call last):
  File "mwe.py", line 29, in <module>
    solve(a == L, q_)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in <listcomp>
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 460, in _compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/uflacsrepresentation.py", line 139, in compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py", line 353, in build_uflacs_ir
    rtol=p["table_rtol"], atol=p["table_atol"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 553, in build_optimized_tables
    modified_terminals, rtol=rtol, atol=atol)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 385, in build_element_tables
    name_ignored = add_table(subres)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 362, in add_table
    entitytype, local_derivatives, flat_component)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 217, in get_ffc_table_values
    num_dofs, num_points = component_tables[0].shape
ValueError: too many values to unpack (expected 2)

Comments (2)

  1. Grant Bruer reporter

    Interestingly, it seems to work fine if I use vector spaces instead of scalar spaces.

    from fenics import *
    
    mesh = UnitSquareMesh(4, 4)
    V = VectorFunctionSpace(mesh, 'P', 1)
    
    u = Function(V)
    
    # Create quadrature space.
    quad_params = {
        'quadrature_degree': 2,
        'quadrature_rule': 'default',
        'representation': 'quadrature'
    }
    shape = V.ufl_element().value_shape()
    quad_el = TensorElement('Quadrature', mesh.ufl_cell(), shape=shape,
        degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])
    Q = FunctionSpace(mesh, quad_el)
    
    qdx = dx(metadata=quad_params)
    
    # Set up a system for projecting u onto Q.
    v = TestFunction(Q)
    q = TrialFunction(Q)
    
    a = inner(q, v) * qdx
    L = inner(u, v) * qdx
    
    # Solve the system.
    q_ = Function(Q)
    solve(a == L, q_) # Does not crash.
    
  2. Log in to comment