Defining FunctionSpace over a Quadrature Element in dolfinx

dolfinx
1

Hi,
I am trying to implement some plasticity models in dolfinx and i am following the approach in https://comet-fenics.readthedocs.io/en/latest/demo/2D_plasticity/vonMises_plasticity.py.html . However I run into issues when I try to define a FunctionSpace over the quadrature element. I tried,

Ue = ufl.VectorElement("Lagrange",mesh.ufl_cell(),1)
## Plastic strain Equivalent
Pe = ufl.FiniteElement("Quadrature",mesh.ufl_cell(),degree=2,quad_scheme='default')
## To store stress in the element
Se = ufl.VectorElement("Quadrature",mesh.ufl_cell(),degree=2,dim=6,quad_scheme='default')

U = FunctionSpace(mesh,Ue)
P = FunctionSpace(mesh,Pe)
S = FunctionSpace(mesh,Se)

I got the following error,

Traceback (most recent call last):
  File "J2_plasticity.py", line 111, in <module>
    P = FunctionSpace(mesh,Pe)
  File "/home/potaman/.local/lib/python3.7/site-packages/dolfinx/function.py", line 256, in __init__
    self.ufl_element(), form_compiler_parameters=None, mpi_comm=mesh.mpi_comm())
  File "/home/potaman/.local/lib/python3.7/site-packages/dolfinx/jit.py", line 48, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/potaman/.local/lib/python3.7/site-packages/dolfinx/jit.py", line 115, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_elements([ufl_object], parameters=p, cache_dir=cache_dir, **cffi_options)
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/codegeneration/jit.py", line 130, in compile_elements
    cffi_extra_compile_args, cffi_verbose, cffi_debug)
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/codegeneration/jit.py", line 287, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix="JIT", parameters=parameters)
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/compiler.py", line 111, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, parameters, visualise)
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/ir/representation.py", line 119, in compute_ir
    for e in analysis.unique_elements
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/ir/representation.py", line 119, in <listcomp>
    for e in analysis.unique_elements
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/ir/representation.py", line 186, in _compute_element_ir
    ir["evaluate_basis"] = _evaluate_basis(ufl_element, fiat_element, epsilon)
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/ir/representation.py", line 624, in _evaluate_basis
    "max_degree": max([e.degree() for e in elements])
  File "/home/potaman/.local/lib/python3.7/site-packages/ffcx/ir/representation.py", line 624, in <listcomp>
    "max_degree": max([e.degree() for e in elements])
  File "/usr/lib/python3/dist-packages/FIAT/tensor_product.py", line 394, in degree
    return self.element.degree()
AttributeError: 'QuadratureElement' object has no attribute 'degree'

Is there anything I am doing wrong? or is it not implemented yet in dolfinx?

Thanks,
Subramanya.

Update : This is getting curious, I changed the code to

e = ufl.VectorElement("Lagrange",mesh.ufl_cell(),1)
## Plastic strain Equivalent
Pe = ufl.FiniteElement("Q",mesh.ufl_cell(),degree=2,quad_scheme='default')
## To store stress in the element
Se = ufl.VectorElement("Q",mesh.ufl_cell(),degree=2,dim=6,quad_scheme='default')

And this works.

2

I cannot reproduce this with the latest docker image of dolfinx/real or dolfinx/complex.
The code i ran for reference:

from dolfinx import UnitSquareMesh, FunctionSpace
import ufl
from mpi4py import MPI
mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
Ue = ufl.VectorElement("Lagrange",mesh.ufl_cell(),1)
## Plastic strain Equivalent
Pe = ufl.FiniteElement("Quadrature",mesh.ufl_cell(),degree=2,quad_scheme='default')
## To store stress in the element
Se = ufl.VectorElement("Quadrature",mesh.ufl_cell(),degree=2,dim=6,quad_scheme='default')

U = FunctionSpace(mesh,Ue)
P = FunctionSpace(mesh,Pe)
S = FunctionSpace(mesh,Se)
3

I tried some more experiments. The thing that is causing the issue is my mesh definition.

mesh = dolfinx.generation.BoxMesh(dolfinx.MPI.comm_world,[p0,p1],[int(l/el_size), int(b/el_size),int(h/el_size)]
                         ,cell_type=CellType.hexahedron)

With this, the code only works with “Q” elements, but with tetrahedra, it only works with “Quadrature”. Is there a difference between how the two are defined?

Thanks!
Subramanya.

4

See http://femtable.org/

5

Oh. So this overrides the ufl documentation for Quadrature elements.

image
image829×183 28.1 KB

Should I use the “dQ” elements to approximate a quadrature element for holding the stresses for the quadrature?

Thanks!
Subramanya.

6

I see that Quadrature Elements are not supported for quadrilateral/hexahedral elements even in dolfin. Is there a plan to add them?

7

It is being worked on by @mscroggs in dolfinx. I guess he can elaborate on what is currently supported. See: https://github.com/FEniCS/dolfinx/pull/946/

1 Like
8

To elaborate, Seems like Q elements work according to the following lines:

parametrize_cell_types_tp
@pytest.mark.parametrize("family", ["Q"])
@pytest.mark.parametrize("degree", [2, 3, 4])
def test_P_tp(family, degree, cell_type, datadir):
    mesh = get_mesh(cell_type, datadir)
    V = FunctionSpace(mesh, (family, degree))
    run_scalar_test(mesh, V, degree)
1 Like