Interpolation into Quadrature Elements

[bhaveshshrimali]

Hi @michalhabera

Many thanks for creating this wrapper. I wanted to check the feasibility of using dolfinx to simulate nonlinear viscoelasticity (a bit similar to plasticity but with some differences). I started by taking a look at the plasticity example and https://github.com/FEniCS/dolfinx/blob/c9f3c459a75f02abcf255a9dd05185350eb1edde/python/test/unit/fem/test_expression.py#L170-L253

but wanted to check if it is possible to interpolate into a quadrature function. I do see an open issue. Has there been any work on the same? In any case, the following should give an idea of what I want to do

from dolfinx.mesh import create_unit_cube
import dolfinx.fem as fem
from mpi4py import MPI
from ufl import TensorElement

metadata = {"quadrature_degree":5}
mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5)
Wq = TensorElement("Quadrature", mesh.ufl_cell(), degree=metadata["quadrature_degree"], symmetry=True, quad_scheme="default")
Vq = fem.FunctionSpace(mesh, Wq)
Cv = fem.Function(Vq, name="Cv")

# now something like?
Cv.interpolate(lambda x: np.eye(x))  # the initial state is identity

Please let me know if you see a way to implement it. Thanks in advance!

[bhaveshshrimali]

Alright, this does seem to do the trick for me at least:

from dolfinx.mesh import create_unit_cube
import dolfinx.fem as fem
from mpi4py import MPI
from ufl import TensorElement, SpatialCoordinate, grad, dx
import numpy as np
from petsc4py.PETSc import ScalarType as st
from dolfiny.interpolation import interpolate as dlinterp
metadata = {"quadrature_degree":5}
mesh = create_unit_cube(MPI.COMM_WORLD, 5, 5, 5)
Wq = TensorElement("Quadrature", mesh.ufl_cell(), degree=metadata["quadrature_degree"], symmetry=True, quad_scheme="default")
Vq = fem.FunctionSpace(mesh, Wq)
Vl = fem.TensorFunctionSpace(mesh, ("CG", 1), symmetry=True)
Cv_l = fem.Function(Vl)
xs = SpatialCoordinate(mesh)
expr = grad(xs)
Cv = fem.Function(Vq)
dlinterp(expr, Cv)

# verify
for i in range(3):
    for j in range(3):
        Cv_ij = fem.assemble_scalar(fem.form(Cv[i, j]*dx))
        print(f"Cv_{i+1}_{j+1} = {Cv_ij}")

# Cv_1_1 = 1.0000000000000389
# Cv_1_2 = 0.0
# Cv_1_3 = 0.0
# Cv_2_1 = 0.0
# Cv_2_2 = 1.0000000000000389
# Cv_2_3 = 0.0
# Cv_3_1 = 0.0
# Cv_3_2 = 0.0
# Cv_3_3 = 1.0000000000000389

Let me know if you see a better way to do this. Thanks!

[michalhabera]

Hi @bhaveshshrimali !

Our monolithic plasticity demo indeed uses interpolation into Quadrature function space, see https://github.com/michalhabera/dolfiny/blob/master/demo/plasticity/solid_plasticity_monolithic.py#L255 Such step is necessary to avoid discretisation error.

In dolfinx, interpolation into spaces which are not based on basix element is an open issue, as you correctly pointed out. Feel free to use dolfiny interpolation in the meantime. I will close this issue, but let us know if there still is a missing functionality.

[bhaveshshrimali]

Hi @michalhabera ,

The above example workaround fails with dolfinx version 0.4.2.0 (basix version 0.4.3.0). Is there a way to assign a function belonging to a tensor function space a value of identity ? Cv = I, in context of the above example. I also tried the following but failed:

def f(x):
    one = np.ones_like(x)
    zero = np.zeros_like(x)
    return np.array(
        [[one, zero, zero], [zero, one, zero], [zero, zero, one]]
    )
dlinterp(f, Cv)

[michalhabera]

Hi @bhaveshshrimali, we have our own "versioning" system for FEniCSx libraries, as you can see at the top of dockerfile https://github.com/michalhabera/dolfiny/blob/master/Dockerfile

Arguably, we should use standard versions as you mentions. Will think about it.

[bhaveshshrimali]

THanks @michalhabera , I hadn't seen the dockerfile. Just pulled the master branch of dolfiny and checked if everything was working as expected but apparently not.

[marchirschvogel]

Hi, I'm also following this for quite a while, cf. the minimal example https://fenicsproject.discourse.group/t/quadrature-function-space/5117/6 I posted. I'm hesitating of porting my code to dolfiny since I'm not sure if everything written using dolfinx would work in y straight away? Is there any rough estimate on if and when the interpolation to quadrature elements will be added to x? Thanks! Best, Marc