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

I/O
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])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    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

2

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
3

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!

4

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
5

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!