Ufl.conditional doesn't return precise results when using ('CG', 1)

Hello everyone,

I am working on a non-local plasticity example. For the obtained solution of the plasticity I need to define a irreversability constraint, where the previous solution is supposed to be the lower threshhold (here an simple function called x0_func). This also called ad-hoc condition. That means at time n+1, we have:

alpha_n+1 > alpha_n

In order to realize that I am using ufl.conditional. When I define the function in the space of (‘DG’, 0), I get really good results. However, this is not the case when the functions are defined in (‘CG’, 1).

The code and the outputs are as follows:

``````import dolfinx
from dolfinx import mesh
from dolfinx.fem import (Function, FunctionSpace)
from petsc4py import PETSc
from mpi4py import MPI
import ufl as ufl
import numpy as np
from dolfinx.io import XDMFFile

L = 100.0
H = 100.0
dL = 30
dH = 30
domain = mesh.create_box(MPI.COMM_WORLD,[[0,0,0], [L, L, H]], [dL, dL, dH], mesh.CellType.tetrahedron)
x = ufl.SpatialCoordinate(domain)

# For plasticity problem
V_pl = FunctionSpace(domain, ('CG', 1))
#V_pl = FunctionSpace(domain, ('DG', 0)) #precise
pl = Function(V_pl)
with pl.vector.localForm() as pl_loc:
pl_loc.set(-1.0)

x0_func = Function(V_pl)
x0_ufl = x[0]
x0 = lambda x: eval(str(x0_ufl))
x0_func.interpolate(x0)
print("min of x0_func:", min(x0_func.vector.array), flush=True)
with XDMFFile(MPI.COMM_WORLD,  "test.xdmf", "w") as xdmf:#, encoding=XDMFFile.Encoding.HDF5) #FIXME
xdmf.write_mesh(domain)
xdmf.write_function(x0_func)

def project(v, target_func, bcs=[]):
# Ensure we have a mesh and attach to measure
V = target_func.function_space
dx = ufl.dx(V.mesh)

# Define variational problem for projection
w = ufl.TestFunction(V)
Pv = ufl.TrialFunction(V)
a = dolfinx.fem.form(ufl.inner(Pv, w) * dx)
L = dolfinx.fem.form(ufl.inner(v, w) * dx)

# Assemble linear system
A = dolfinx.fem.petsc.assemble_matrix(a, bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(L)
dolfinx.fem.petsc.apply_lifting(b, [a], [bcs])
dolfinx.fem.petsc.set_bc(b, bcs)

# Solve linear system
solver = PETSc.KSP().create(A.getComm())
solver.setOperators(A)
solver.solve(b, target_func.vector)

def Pl_max(pl_, pln_):
return ufl.conditional(ufl.gt(pl_, pln_), pl_, pln_)

project(Pl_max(pl, x0_func), pl)
print("min of plasticity:", min(pl.vector.array), flush=True)
``````

Output for V_pl = FunctionSpace(domain, (‘CG’, 1))

``````min of x0_func: -5.244705219161143e-17
min of plasticity: -2.950809850125471e-05
``````

Output for V_pl = FunctionSpace(domain, (‘DG’, 0))

``````min of x0_func: 0.8333333333333334
min of plasticity: 0.8333333333334109
``````

Thank you

I’m not really sure what you are expecting here.

You have set the minimum bound og `pl` to -1, thus your condition states that `x0_func` is always greater than than `pl`, and thus the resulting plasticity will always be equal to `x0_func`.

Let’s consider why you get different results for `x0`.

Any degree of freedom in a DG 0 space is placed in the cell center, while the dofs in CG-1 is placed at the vertices. We can then look at what point `x0` is at its smallest:

``````import numpy as np
import dolfinx
from dolfinx import mesh
from dolfinx.fem import (Function, FunctionSpace)
from petsc4py import PETSc
from mpi4py import MPI
import ufl as ufl
from dolfinx.io import XDMFFile

def compute_x0(V_pl):
domain = V_pl.mesh
x = ufl.SpatialCoordinate(domain)

x0_func = Function(V_pl)
x0_ufl = x[0]
def x0(x): return eval(str(x0_ufl))

x0_func.interpolate(x0)
return x0_func

L = 100.0
H = 100.0
dL = 30
dH = 30
domain = mesh.create_box(MPI.COMM_WORLD, [[0, 0, 0], [L, L, H]], [
dL, dL, dH], mesh.CellType.tetrahedron)

# For plasticity problem
V_CG = FunctionSpace(domain, ('CG', 1))
x0 = compute_x0(V_CG)
coord_CG = V_CG.tabulate_dof_coordinates()
pos = np.argmin(x0.x.array)
print(coord_CG[pos], x0.x.array[pos])

with XDMFFile(domain.comm, "CG.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(x0)

V_DG = FunctionSpace(domain, ('DG', 0))  # precise
x0_dg = compute_x0(V_DG)
coord_DG = V_DG.tabulate_dof_coordinates()
pos = np.argmin(x0_dg.x.array)
print(coord_DG[pos], x0_dg.x.array[pos])

with XDMFFile(domain.comm, "DG.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(x0_dg)
``````

which gives you

``````[0.00000000e+00 1.85037171e-16 1.85037171e-16] 0.0
[0.83333333 2.5        1.66666667] 0.8333333333333334
``````

i.e. the CG 1 function is smallest at [0,0,0]
while DG 0 is smallest at [0.83,2.5,1.6667]
If you look at all the dof coordinates, you can check that the smallest x-value for each spaces is:

``````print(min(coord_CG[:,0]))
print(min(coord_DG[:,0]))
``````

giving

``````0.0
0.8333333333333334
``````

Thank you very much for this explanation. That would have been my next question

By setting `pl `to -1, I wanted to force the plasticity to be equal to `x0_func` so that I could test the accuracy of the ad-hoc condition by comparing `pl` and `x0_func`.

I found that the lowest values of `pl` and `x0_func` were much closer in DG-0 than in CG-1. I wondered why this was. I’m sorry if I didn’t make myself clear.

What happens in my full simulation is that I get a negative plasticity of up to `-4e-3`. This should not be possible because of the ad-hoc condition and the initial values (`pl_n` set to `0`). However, when I describe the functions in DG-0, I get a minimum plasticity value of 0, as expected.

Thanks a lot!

Without having an actual example of this behavior at hand, I cannot really give you any pointers.

You need to think about what `project` and `interpolate` does, and how they are affected by using a `ufl.conditional`. If you do a projection with a ufl-conditional, the condition is evalutated at a quadrature point when assembled into the RHS. This is different from what happens in an interpolation, where you would only evaluate the condition at a dof coordinate/interpolation point.

1 Like

Thank you! I just tested your suggestion on my full simulation and it works like a charm.

For future readers, the now working example:

``````import dolfinx
from dolfinx import mesh
from dolfinx.fem import (Expression, Function, FunctionSpace)
from petsc4py import PETSc
from mpi4py import MPI
import ufl as ufl
import numpy as np
from dolfinx.io import XDMFFile

L = 100.0
H = 100.0
dL = 30
dH = 30
domain = mesh.create_box(MPI.COMM_WORLD,[[0,0,0], [L, L, H]], [dL, dL, dH], mesh.CellType.tetrahedron)
x = ufl.SpatialCoordinate(domain)

# For plasticity problem
V_pl = FunctionSpace(domain, ('CG', 1))
#V_pl = FunctionSpace(domain, ('DG', 0))
pl = Function(V_pl)
with pl.vector.localForm() as pl_loc:
pl_loc.set(-1.0)

x0_func = Function(V_pl)
x0_ufl = x[0]
x0 = lambda x: eval(str(x0_ufl))
x0_func.interpolate(x0)
print("min of x0_func:", min(x0_func.vector.array), flush=True)

def Pl_max(pl_, pln_):
return ufl.conditional(ufl.gt(pl_, pln_), pl_, pln_)

expr = Expression(Pl_max(pl, x0_func), V_pl.element.interpolation_points)
pl.interpolate(expr)
print("min of plasticity:", min(pl.vector.array), flush=True)
``````

Output for V_pl = FunctionSpace(domain, (‘CG’, 1))

``````min of x0_func: -5.244705219161143e-17
min of plasticity: -5.244705219161143e-17
``````

Output for V_pl = FunctionSpace(domain, (‘DG’, 0))

``````min of x0_func: 0.8333333333333334
min of plasticity: 0.8333333333333334
``````

Have a nice day!