Using ufl.cell_avg for B-bar correction of strain

General
1

I am trying to implement a B-bar formulation similar to the one described in B-bar, where the ufl expression describing the strain would need to contain the mean value of the voulmetric strain across the cell.

\bar \varepsilon_\mathrm{vol} = \frac{1}{V_e}\int_{V_e}\frac{1}{3} \mathrm {tr}(\varepsilon)\mathrm d V\\ \bar\varepsilon = \varepsilon+(\bar \varepsilon_\mathrm{vol}-\frac{1}{3}\mathrm{tr}(\varepsilon))\cdot I

I have tried the following code using ufl.cell_avg to calculate the mean volumetric strain


import numpy as np
import dolfinx as dfx
import ufl

from mpi4py import MPI

def b_bar_strain(u):
    eps = ufl.sym(ufl.grad(u))
    vol = (1/3) * ufl.cell_avg(ufl.tr(eps))
    return eps +(vol - (1/3) * ufl.tr(eps)) * ufl.Identity(2)

mesh = dfx.mesh.create_rectangle(MPI.COMM_WORLD,np.array([[0,0],[1., 1.]]), [1,1], cell_type=dfx.mesh.CellType.quadrilateral)

displacement_space = dfx.fem.VectorFunctionSpace(mesh,("CG",1))
stress_space = dfx.fem.TensorFunctionSpace(mesh,("CG",1),shape=(2,2))

stress = dfx.fem.Function(stress_space)
u_= ufl.TestFunction(displacement_space)

stress.vector.array[:] = np.random.random(stress.vector.array.size)

avg = dfx.fem.form(ufl.inner(stress,b_bar_strain(u_))*ufl.dx)
f = dfx.fem.petsc.assemble_vector(avg)

This does not work and I assume it is because of ufl.cell_avg, because a reduced MWE throws the same error


import numpy as np
import dolfinx as dfx
import ufl

from mpi4py import MPI

mesh = dfx.mesh.create_rectangle(MPI.COMM_WORLD,np.array([[0,0],[1., 1.]]), [1,1], cell_type=dfx.mesh.CellType.quadrilateral)

V = dfx.fem.FunctionSpace(mesh,("CG",1))
v = dfx.fem.Function(V)
v_= ufl.TestFunction(V)

v.vector.array[:] = np.random.random(v.vector.array.size)

avg = dfx.fem.form(ufl.inner(v,ufl.cell_avg(v_))*ufl.dx)
f = dfx.fem.petsc.assemble_vector(avg)

I get the following error message

Traceback (most recent call last):
  File "/home/srosenbu/git-packages/fenicsX-constitutive/test/test_b_bar.py", line 16, in <module>
    avg = dfx.fem.form(ufl.inner(v,ufl.cell_avg(v_))*ufl.dx)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 166, in form
    return _create_form(form)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 161, in _create_form
    return _form(form)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 135, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 168, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 232, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, parameters=parameters)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, parameters, visualise)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/ir/representation.py", line 198, in compute_ir
    irs = [
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/ir/representation.py", line 199, in <listcomp>
    _compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/ir/representation.py", line 480, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ffcx/ir/integral.py", line 65, in compute_integral_ir
    expression = balance_modifiers(expression)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ufl/algorithms/balancing.py", line 76, in balance_modifiers
    return map_expr_dag(mf, expr)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 99, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ufl/algorithms/balancing.py", line 63, in _modifier
    return balance_modified_terminal(expr)
  File "/home/srosenbu/mambaforge/envs/dolfinx/lib/python3.10/site-packages/ufl/algorithms/balancing.py", line 31, in balance_modified_terminal
    assert expr._ufl_is_terminal_modifier_
AssertionError

I am using dolfinx 0.5.0 which was installed from conda-forge.
I have already looked into the ufl-documentation, but I am not sure what exactly a terminal modifier is and why this would matter.

2

cell_avg is currently not supported. See:

Which tried to fix it but failed.

3

Thank you for looking into it. I will just watch if the pull request gets any updates in the future

4

Just out of curiosity, were you able to implement your own version of this or any kind of work around? It seems that cell_avg is still not supported

5

Hi @driley,

In short: I did not find a workaround for cell_avg, but I tried to solve the \bar B correction the following way:
In my simulation, I evaluate strains at the quadrature points and write them to a numpy array. Then I took this large vector of strains and wrote a function that replaces the volumetric part of the strain with an average of the volumetric strain of all quadrature points in a cell. With this new strain I calculated the stresses and wrote them to another quadrature space which I used to determine internal forces (alternatively you could build a stiffness matrix, but I only do explicit dynamics).

It would look somewhat like this:

# u represents the displacements, cells an array with the indices of all cells,
# and quadrature_points comes from the basix quadrature rule
# stress is a quadrature function
# measure is the measure corresponding to the used quadrature rule
strain_expr = dolfinx.fem.Expression(ufl.sym(ufl.grad(u)), quadrature_points)
strain_array = strain_expr.eval(cells)
strain_array = b_bar_correction(strain_array)
stress.vector.array[:] = constitutive_model(strain_array)
f_int_form = ufl.inner(ufl.sym(grad(test_function)), stress) * measure

I left a lot of things out. I don’t know how useful my approach is, since I just take the mean of all quadrature points in a cell which is not exactly the \bar B approach, and also I don’t use the correction later in the weak form. But I think I read somewhere in the LS-DYNA Theory manual, that the correction is not actually needed in the weak form, but I didn’t prove it for myself.
But I didn’t actually need the correction, therefore all my trys stop there.
If you are interested in some of the code I would share it. (It is already on github, but it is really messy, parts of it are in C++ with python bindings and it needs a lot of context, which is why I don’t want to just link it here)