Projecting into a quadrature space

Legacy DOLFIN
1

I am struggling to project a function into a FunctionSpace of the “Quadrature” family in classic legacy dolfin. Has anyone had any experience with this?

import dolfin as df

mesh = df.UnitSquareMesh(10, 10)

# Step 1: Create a Quadrature element specifying the quadrature scheme
quad_element = df.FiniteElement("Quadrature",
                                 mesh.ufl_cell(),
                                 degree=4,
                                 quad_scheme="default")

# Step 2: Create a VectorElement from the Quadrature element
vector_element = df.VectorElement(quad_element)

# Step 3: Create the VectorFunctionSpace with the VectorElement
V = df.FunctionSpace(mesh, vector_element)

u = df.TrialFunction(V)
uf = df.Function(V)
v = df.TestFunction(V)

f = df.Expression("sin(x[0])*cos(x[1])", degree=2)

df.solve(df.inner(u,v)*df.dx == df.inner(f,v)*df.dx, uf)
#df.assemble(df.inner(u,v)*df.dx) #This one also fails
print("Quarature Test finished successfully!")

When I run the script I get the following error

Missing quad_scheme in quadrature element.
Traceback (most recent call last):
File “/home/gabriel/mywork/research_studies/CHESRA/chesra_gb/experiments/unload/quadrature_elem_test.py”, line 15, in
V = df.FunctionSpace(mesh, vector_element)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/functionspace.py”, line 31, in init
self._init_from_ufl(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/functionspace.py”, line 42, in _init_from_ufl
ufc_element, ufc_dofmap = ffc_jit(element, form_compiler_parameters=None,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 50, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 100, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 149, in compile_element
return compile_ufl_objects(elements, “element”, object_names,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File “/usr/lib/python3/dist-packages/ffc/analysis.py”, line 125, in analyze_ufl_objects
error(“Missing quad_scheme in quadrature element.”)
File “”, line 1, in
File “/usr/lib/python3/dist-packages/ufl_legacy/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
Exception: Missing quad_scheme in quadrature element.

2

Is there a reason not for using:

import dolfin as df

mesh = df.UnitSquareMesh(10, 10)

# Step 1: Create a Quadrature element specifying the quadrature scheme
quad_element = df.VectorElement(
    "Quadrature", mesh.ufl_cell(), degree=4, quad_scheme="default"
)
V = df.FunctionSpace(mesh, quad_element)
3

Thanks for the suggestion Jørgen, and for answering a dolfin-legacy question. I have updated the code as you suggested

import dolfin as df

mesh = df.UnitSquareMesh(10, 10)

# Step 1: Create a Quadrature element specifying the quadrature scheme

quad_element = df.VectorElement("Quadrature", mesh.ufl_cell(), degree=4, quad_scheme="default")

V = df.FunctionSpace(mesh, quad_element)

u = df.TrialFunction(V)

uf = df.Function(V)

v = df.TestFunction(V)

f = df.Expression(("sin(x[0])*cos(x[1])", "cos(x[0])*sin(x[1])"), degree=4)

df.solve(df.inner(u,v)*df.dx == df.inner(f,v)*df.dx, uf)

print("Quarature Test finished successfully!")

It seems to throw a different error this time “mismatch of quadrature points”.

Traceback (most recent call last):
File “/home/gabriel/mywork/research_studies/CHESRA/chesra_gb/experiments/unload/quadrature_elem_test.py”, line 26, in
df.solve(df.inner(u,v)*df.dx == df.inner(f,v)*df.dx, uf)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py”, line 268, in _solve_varproblem
problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py”, line 58, in init
L = Form(L, form_compiler_parameters=form_compiler_parameters)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py”, line 43, in init
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 50, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py”, line 100, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/lib/python3/dist-packages/ffc/jitcompiler.py”, line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 142, in compile_form
return compile_ufl_objects(forms, “form”, object_names,
File “/usr/lib/python3/dist-packages/ffc/compiler.py”, line 190, in compile_ufl_objects
ir = compute_ir(analysis, prefix, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 182, in compute_ir
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 182, in
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File “/usr/lib/python3/dist-packages/ffc/representation.py”, line 455, in _compute_integral_ir
ir = r.compute_integral_ir(itg_data,
File “/usr/lib/python3/dist-packages/ffc/uflacs/uflacsrepresentation.py”, line 132, in compute_integral_ir
uflacs_ir = build_uflacs_ir(itg_data.domain.ufl_cell(),
File “/usr/lib/python3/dist-packages/ffc/uflacs/build_uflacs_ir.py”, line 350, in build_uflacs_ir
build_optimized_tables(num_points, quadrature_rules,
File “/usr/lib/python3/dist-packages/ffc/uflacs/elementtables.py”, line 564, in build_optimized_tables
build_element_tables(num_points, quadrature_rules,
File “/usr/lib/python3/dist-packages/ffc/uflacs/elementtables.py”, line 398, in build_element_tables
name_ignored = add_table(subres)
File “/usr/lib/python3/dist-packages/ffc/uflacs/elementtables.py”, line 371, in add_table
tables[name] = get_ffc_table_values(
File “/usr/lib/python3/dist-packages/ffc/uflacs/elementtables.py”, line 194, in get_ffc_table_values
tbl = fiat_element.tabulate(deriv_order, entity_points)[derivative_counts]
File “/usr/lib/python3/dist-packages/FIAT/quadrature_element.py”, line 56, in tabulate
raise AssertionError(“Mismatch of quadrature points!”)
AssertionError: Mismatch of quadrature points!

4

The following works for me:

import dolfin as df

df.parameters["form_compiler"]["representation"] = "quadrature"
mesh = df.UnitSquareMesh(10, 10)

# Step 1: Create a Quadrature element specifying the quadrature scheme

quad_element = df.VectorElement(
    "Quadrature", mesh.ufl_cell(), degree=4, quad_scheme="default"
)

V = df.FunctionSpace(mesh, quad_element)

u = df.TrialFunction(V)

uf = df.Function(V)

v = df.TestFunction(V)

f = df.Expression(("sin(x[0])*cos(x[1])", "cos(x[0])*sin(x[1])"), degree=4)

dx = df.dx(metadata={"quadrature_degree": 4})
df.solve(
    df.inner(u, v) * dx == df.inner(f, v) * dx,
    uf,
)

print("Quarature Test finished successfully!")
5

Thanks again for the help Jørgen! I ran your version of the code successfully.
I did get a deprecation warning though:

  Ignoring precision in integral metadata compiled using quadrature representation. Not implemented.
/usr/lib/python3/dist-packages/ffc/jitcompiler.py:234: QuadratureRepresentationDeprecationWarning: 
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** representation instead.                               ***
*** ===================================================== ***
  issue_deprecation_warning()
Solving linear variational problem.
Quarature Test finished successfully!

Can the same code be run with the uflacs representation?

6

No, you cannot. Quadrature spaces were never finalized in “uflacs”.

I would strongly advice moving over to DOLFINx, as there has been alot of development going on since the last release of legacy dolfin in 2019.

7

After running the code, I had a question. @dokken

I saw that the code was running successfully, and I wanted to go through the visualisation process with ParaView.
After writing additional code as shown below,

fid = df.File("output/uf.pvd")
fid << uf

I got the error below, and confirmed that ParaView doesn’t run at all.

File "~/miniconda3/envs/dolfin-2019.1.0/lib/python3.11/site-packages/dolfin/io/__init__.py", line 14, in __lshift__
    self.write(u._cpp_object)
RuntimeError: interpolate_vertex_values: Function is not supported/implemented for QuadratureElement.

Is there any way to solve this error? Or is this an issue that is bound to exist in legacy FEniCS?

8

This is something you can’t avoid in legacy fenics without writing a point cloud writer, similar to what we have in dolfinx/scifem:
https://scientificcomputing.github.io/scifem/examples/xdmf_point_cloud.html

9

Thank you for reply! @dokken! :raised_hands:

If so, does this mean that we should use scifem.xdmf.XDMFFile rather than the df.XDMFFile built into legacy FEniCS?
(When I try to install with conda install -c conda-forge scifem in the conda envs, I get a dependency issue. Is scifem a tool developed for the dolfinx env?)

If it’s not a specific problem, it seems like it’s because I have no concept of point cloud writer.

10

It wont work with legacy fenics. You will have to write your own version of this XDMFFile handler for legacy dolfin in you want to visualize quadrature-spaces.

1 Like