Project function in dolfinx

Regarding the FEniCSx demo: Elasto-plastic analysis of a 2D von Mises material - Elasto-plastic analysis of a 2D von Mises material — Computational Mechanics Numerical Tours with FEniCSx, I have two questions:

  1. In the legacy demo (Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation), it is said the plastic strain (in quadrature function space) must first be projected onto the previously defined DG FunctionSpace: p_avg.assign(project(p, P0)); file_results.write(p_avg, t), how to do this in FEniCSx?
  2. Following item1, I defined a project function, but when trying to save the vector function in VTK file, the kernel died. please see the following code.
import numpy as np
import matplotlib.pyplot as plt

import gmsh
from mpi4py import MPI
import ufl
import basix
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from petsc4py import PETSc

hsize = 0.2

Re = 1.3
Ri = 1.0

gmsh.initialize()
gdim = 2
model_rank = 0
if MPI.COMM_WORLD.rank == 0:
    gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
    gmsh.model.add("Model")

    geom = gmsh.model.geo
    center = geom.add_point(0, 0, 0)
    p1 = geom.add_point(Ri, 0, 0)
    p2 = geom.add_point(Re, 0, 0)
    p3 = geom.add_point(0, Re, 0)
    p4 = geom.add_point(0, Ri, 0)

    x_radius = geom.add_line(p1, p2)
    outer_circ = geom.add_circle_arc(p2, center, p3)
    y_radius = geom.add_line(p3, p4)
    inner_circ = geom.add_circle_arc(p4, center, p1)

    boundary = geom.add_curve_loop([x_radius, outer_circ, y_radius, inner_circ])
    surf = geom.add_plane_surface([boundary])

    geom.synchronize()

    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", hsize)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hsize)

    gmsh.model.addPhysicalGroup(gdim, [surf], 1)
    gmsh.model.addPhysicalGroup(gdim - 1, [x_radius], 1, name="bottom")
    gmsh.model.addPhysicalGroup(gdim - 1, [y_radius], 2, name="left")
    gmsh.model.addPhysicalGroup(gdim - 1, [inner_circ], 3, name="inner")

    gmsh.model.mesh.generate(gdim)

domain, _, facets = io.gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim
)
gmsh.finalize()

E = fem.Constant(domain, 70e3)  # in MPa
nu = fem.Constant(domain, 0.3)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2.0 / (1 + nu)
sig0 = fem.Constant(domain, 250.0)  # yield strength in MPa
Et = E / 100.0  # tangent modulus
H = E * Et / (E - Et)  # hardening modulus


deg_u = 2
shape = (gdim,)
V = fem.functionspace(domain, ("P", deg_u, shape))

Vx, _ = V.sub(0).collapse()
Vy, _ = V.sub(1).collapse()
bottom_dofsy = fem.locate_dofs_topological((V.sub(1), Vy), gdim - 1, facets.find(1))
top_dofsx = fem.locate_dofs_topological((V.sub(0), Vx), gdim - 1, facets.find(2))


# used for post-processing
def bottom_inside(x):
    return np.logical_and(np.isclose(x[0], Ri), np.isclose(x[1], 0))


bottom_inside_dof = fem.locate_dofs_geometrical((V.sub(0), Vx), bottom_inside)[0]

u0x = fem.Function(Vx)
u0y = fem.Function(Vy)
bcs = [
    fem.dirichletbc(u0x, top_dofsx, V.sub(0)),
    fem.dirichletbc(u0y, bottom_dofsy, V.sub(1)),
]

n = ufl.FacetNormal(domain)
q_lim = float(2 / np.sqrt(3) * np.log(Re / Ri) * sig0)

loading = fem.Constant(domain, 0.0)

deg_quad = 2  # quadrature degree for internal state variable representation
W0e = basix.ufl.quadrature_element(
    domain.basix_cell(), value_shape=(), scheme="default", degree=deg_quad
)
We = basix.ufl.quadrature_element(
    domain.basix_cell(), value_shape=(4,), scheme="default", degree=deg_quad
)
W = fem.functionspace(domain, We)
W0 = fem.functionspace(domain, W0e)

sig = fem.Function(W)
sig_old = fem.Function(W)
n_elas = fem.Function(W)
beta = fem.Function(W0)
p = fem.Function(W0, name="Cumulative_plastic_strain")
dp = fem.Function(W0)
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="Iteration_correction")
Du = fem.Function(V, name="Current_increment")
v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

P0 = fem.functionspace(domain, ("DG", 0))
p_avg = fem.Function(P0, name="Plastic_strain")

def eps(v):
    e = ufl.sym(ufl.grad(v))
    return ufl.as_tensor([[e[0, 0], e[0, 1], 0], [e[0, 1], e[1, 1], 0], [0, 0, 0]])


def elastic_behavior(eps_el):
    return lmbda * ufl.tr(eps_el) * ufl.Identity(3) + 2 * mu * eps_el


def as_3D_tensor(X):
    return ufl.as_tensor([[X[0], X[3], 0], [X[3], X[1], 0], [0, 0, X[2]]])


def to_vect(X):
    return ufl.as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1]])

ppos = lambda x: ufl.max_value(x, 0)


def constitutive_update(Δε, old_sig, old_p):
    sig_n = as_3D_tensor(old_sig)
    sig_elas = sig_n + elastic_behavior(Δε)
    s = ufl.dev(sig_elas)
    sig_eq = ufl.sqrt(3 / 2.0 * ufl.inner(s, s))
    f_elas = sig_eq - sig0 - H * old_p
    dp = ppos(f_elas) / (3 * mu + H)
    n_elas = s / sig_eq * ppos(f_elas) / f_elas
    beta = 3 * mu * dp / sig_eq
    new_sig = sig_elas - beta * s
    return to_vect(new_sig), to_vect(n_elas), beta, dp

def sigma_tang(eps):
    N_elas = as_3D_tensor(n_elas)
    return (
        elastic_behavior(eps)
        - 3 * mu * (3 * mu / (3 * mu + H) - beta) * ufl.inner(N_elas, eps) * N_elas
        - 2 * mu * beta * ufl.dev(eps)
    )

ds = ufl.Measure("ds", domain=domain, subdomain_data=facets)
dx = ufl.Measure(
    "dx",
    domain=domain,
    metadata={"quadrature_degree": deg_quad, "quadrature_scheme": "default"},
)

Residual = ufl.inner(eps(u_), as_3D_tensor(sig)) * dx - ufl.inner(
    -loading * n, u_
) * ds(3)
tangent_form = ufl.inner(eps(v), sigma_tang(eps(u_))) * dx

basix_celltype = getattr(basix.CellType, domain.topology.cell_type.name)
quadrature_points, weights = basix.make_quadrature(basix_celltype, deg_quad)

map_c = domain.topology.index_map(domain.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)


def interpolate_quadrature(ufl_expr, function):
    expr_expr = fem.Expression(ufl_expr, quadrature_points)
    expr_eval = expr_expr.eval(domain, cells)
    function.x.array[:] = expr_eval.flatten()[:]

class CustomLinearProblem(fem.petsc.LinearProblem):
    def assemble_rhs(self, u=None):
        """Assemble right-hand side and lift Dirichlet bcs.

        Parameters
        ----------
        u : dolfinx.fem.Function, optional
            For non-zero Dirichlet bcs u_D, use this function to assemble rhs with the value u_D - u_{bc}
            where u_{bc} is the value of the given u at the corresponding. Typically used for custom Newton methods
            with non-zero Dirichlet bcs.
        """

        # Assemble rhs
        with self._b.localForm() as b_loc:
            b_loc.set(0)
        fem.petsc.assemble_vector(self._b, self._L)

        # Apply boundary conditions to the rhs
        x0 = [] if u is None else [u.vector]
        fem.petsc.apply_lifting(self._b, [self._a], bcs=[self.bcs], x0=x0, scale=1.0)
        self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        x0 = None if u is None else u.vector
        fem.petsc.set_bc(self._b, self.bcs, x0, scale=1.0)

    def assemble_lhs(self):
        self._A.zeroEntries()
        fem.petsc.assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
        self._A.assemble()

    def solve_system(self):
        # Solve linear system and update ghost values in the solution
        self._solver.solve(self._b, self._x)
        self.u.x.scatter_forward()


tangent_problem = CustomLinearProblem(
    tangent_form,
    -Residual,
    u=du,
    bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)

def project(source_function, target_space):
    """
    Project a function from its current function space onto a target function space.

    Parameters:
    - source_function: fem.Function
      The function to be projected.
    - target_space: fem.FunctionSpace
      The function space onto which the source function should be projected.

    Returns:
    - fem.Function
      The projected function in the target function space.
    """
    # Create the projected function in the target space
    projected_function = fem.Function(target_space)

    # Define the test function and trial function in the target space
    v = ufl.TestFunction(target_space)
    u = ufl.TrialFunction(target_space)

    # Define the variational problem (a is the bilinear form, L is the linear form)
    a = ufl.inner(u, v) * ufl.dx
    L = ufl.inner(source_function, v) * ufl.dx

    # Solve the variational problem to project the source function onto the target space
    problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    projected_function = problem.solve()

    return projected_function
# alternative
def project1(v, V):
    """
    Project a function or expression `v` onto a function space `V`.
    
    Parameters:
    v : dolfinx.fem.Function or ufl.Expression
        The function or expression to project.
    V : dolfinx.fem.FunctionSpace
        The target function space to project onto.
    
    Returns:
    u_proj : dolfinx.fem.Function
        The projected function in the function space `V`.
    """
    # Define trial and test functions in the target function space
    u = ufl.TrialFunction(V)
    w = ufl.TestFunction(V)
    
    # Define the bilinear form (left-hand side)
    a = ufl.inner(u, w) * dx
    
    # Define the linear form (right-hand side) using the function to project
    L = ufl.inner(v, w) * dx
    
    # Convert the forms to dolfinx.fem.Form objects
    a_form = dolfinx.fem.form(a)
    L_form = dolfinx.fem.form(L)
    
    # Assemble the system matrix A and the right-hand side vector b
    A = dolfinx.fem.petsc.assemble_matrix(a_form)
    A.assemble()
    b = dolfinx.fem.petsc.assemble_vector(L_form)
    
    # Create a function in the target space to hold the projected result
    u_proj = dolfinx.fem.Function(V)
    
    # Solve the linear system to compute the projection
    solver = PETSc.KSP().create(V.mesh.comm)
    solver.setOperators(A)
    solver.solve(b, u_proj.vector)
    u_proj.x.scatter_forward()
    
    return u_proj

SigOut = fem.functionspace(domain, ("DG", 0, (4, )))
sigmaout=fem.Function(SigOut, name="Stresses")

vtk = io.VTKFile(domain.comm, "results1/elastoplastic.pvd", "w")
Nitermax, tol = 200, 1e-6  # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr + 1)[1:] ** 0.5
results = np.zeros((Nincr + 1, 3))

# we set all functions to zero before entering the loop in case we would like to reexecute this code cell
sig.vector.set(0.0)
sig_old.vector.set(0.0)
p.vector.set(0.0)
u.vector.set(0.0)
n_elas.vector.set(0.0)
beta.vector.set(0.0)

Δε = eps(Du)
sig_, n_elas_, beta_, dp_ = constitutive_update(Δε, sig_old, p)


for i, t in enumerate(load_steps):
    loading.value = t * q_lim

    # compute the residual norm at the beginning of the load step
    tangent_problem.assemble_rhs()
    nRes0 = tangent_problem._b.norm()
    nRes = nRes0
    Du.x.array[:] = 0

    niter = 0
    while nRes / nRes0 > tol and niter < Nitermax:
        # solve for the displacement correction
        tangent_problem.assemble_lhs()
        tangent_problem.solve_system()

        # update the displacement increment with the current correction
        Du.vector.axpy(1, du.vector)  # Du = Du + 1*du
        Du.x.scatter_forward()

        # interpolate the new stresses and internal state variables
        interpolate_quadrature(sig_, sig)
        interpolate_quadrature(n_elas_, n_elas)
        interpolate_quadrature(beta_, beta)

        # compute the new residual
        tangent_problem.assemble_rhs()
        nRes = tangent_problem._b.norm()

        niter += 1

    # Update the displacement with the converged increment
    u.vector.axpy(1, Du.vector)  # u = u + 1*Du
    u.x.scatter_forward()

    # Update the previous plastic strain
    interpolate_quadrature(dp_, dp)
    p.vector.axpy(1, dp.vector)

    # Update the previous stress
    sig_old.x.array[:] = sig.x.array[:]
    
    p_avg.interpolate(project(p, P0))
    vtk.write_function(p_avg,t)

    sigmaout.interpolate(project(sig,SigOut))
    vtk.write_function(sigmaout,t)
    
vtk.close()

You seem to have implemented the projection already, so I dont really understand this question. See for instance Where to find 'project'-function in dolfinx? - #2 by nate

Try to run your code as a Python script (rather than a notebook) as it might give you a error traceback you could share.

Please also consider simplifying your problem.
Start by removing one of the project functions, remove the loading loop and deduce what line of code is causing your crash.

Hi dokken,

Thank you for your prompt response. I came up with two project functions (project and project1) as I am not sure the correct/proper way to project a function from quadrature function space to FEM function space. But it looks like neither of them work when running the script

p_avg.interpolate(project(p, P0))
vtk.write_function(p_avg,t)
sigmaout.interpolate(project(sig,SigOut))
vtk.write_function(sigmaout,t)

for project, it show the following error:

corrupted size vs. prev_size
[C23WYGT:197508] *** Process received signal ***
[C23WYGT:197508] Signal: Aborted (6)
[C23WYGT:197508] Signal code:  (-6)
[C23WYGT:197508] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7ff67db7c520]
[C23WYGT:197508] [ 1] /lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7ff67dbd09fc]
[C23WYGT:197508] [ 2] /lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7ff67db7c476]
[C23WYGT:197508] [ 3] /lib/x86_64-linux-gnu/libc.so.6(abort+0xd3)[0x7ff67db627f3]
[C23WYGT:197508] [ 4] /lib/x86_64-linux-gnu/libc.so.6(+0x89676)[0x7ff67dbc3676]
[C23WYGT:197508] [ 5] /lib/x86_64-linux-gnu/libc.so.6(+0xa0cfc)[0x7ff67dbdacfc]
[C23WYGT:197508] [ 6] /lib/x86_64-linux-gnu/libc.so.6(+0xa17e2)[0x7ff67dbdb7e2]
[C23WYGT:197508] [ 7] /lib/x86_64-linux-gnu/libc.so.6(+0xa2d2b)[0x7ff67dbdcd2b]
[C23WYGT:197508] [ 8] /lib/x86_64-linux-gnu/libc.so.6(free+0x73)[0x7ff67dbdf453]
[C23WYGT:197508] [ 9] /lib/x86_64-linux-gnu/libpetsc_real.so.3.15(+0xea080)[0x7ff60ed6c080]
[C23WYGT:197508] [10] /usr/lib/petsc/lib/python3/dist-packages/petsc4py/lib/PETSc.cpython-310-x86_64-linux-gnu.so(+0xf152b)[0x7ff669b5d52b]
[C23WYGT:197508] [11] python3(+0x15adae)[0x560974d72dae]
[C23WYGT:197508] [12] /usr/lib/petsc/lib/python3/dist-packages/petsc4py/lib/PETSc.cpython-310-x86_64-linux-gnu.so(+0x29685b)[0x7ff669d0285b]
[C23WYGT:197508] [13] python3(_PyEval_EvalFrameDefault+0x17d7)[0x560974d5ca37]
[C23WYGT:197508] [14] python3(_PyFunction_Vectorcall+0x7c)[0x560974d736ac]
[C23WYGT:197508] [15] python3(_PyObject_FastCallDictTstate+0x16d)[0x560974d6876d]
[C23WYGT:197508] [16] python3(+0x1657a4)[0x560974d7d7a4]
[C23WYGT:197508] [17] python3(_PyObject_MakeTpCall+0x1fc)[0x560974d694cc]
[C23WYGT:197508] [18] python3(_PyEval_EvalFrameDefault+0x7611)[0x560974d62871]
[C23WYGT:197508] [19] python3(_PyFunction_Vectorcall+0x7c)[0x560974d736ac]
[C23WYGT:197508] [20] python3(_PyEval_EvalFrameDefault+0x6d5)[0x560974d5b935]
[C23WYGT:197508] [21] python3(+0x140096)[0x560974d58096]
[C23WYGT:197508] [22] python3(PyEval_EvalCode+0x86)[0x560974e4df66]
[C23WYGT:197508] [23] python3(+0x260e98)[0x560974e78e98]
[C23WYGT:197508] [24] python3(+0x25a79b)[0x560974e7279b]
[C23WYGT:197508] [25] python3(+0x260be5)[0x560974e78be5]
[C23WYGT:197508] [26] python3(_PyRun_SimpleFileObject+0x1a8)[0x560974e780c8]
[C23WYGT:197508] [27] python3(_PyRun_AnyFileObject+0x43)[0x560974e77d13]
[C23WYGT:197508] [28] python3(Py_RunMain+0x2be)[0x560974e6a70e]
[C23WYGT:197508] [29] python3(Py_BytesMain+0x2d)[0x560974e40dfd]
[C23WYGT:197508] *** End of error message ***
Aborted

for project1

Traceback (most recent call last):
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py", line 492, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py", line 455, in _interpolate
    self._cpp_object.interpolate(u, cells, nmm_interpolation_data)  # type: ignore
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. interpolate(self, f: ndarray[dtype=float64, writable=False, shape=(*), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    2. interpolate(self, f: ndarray[dtype=float64, writable=False, shape=(*, *), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    3. interpolate(self, u: dolfinx.cpp.fem.Function_float64, cells: ndarray[dtype=int32, writable=False, shape=(*), order='C'], cell_map: ndarray[dtype=int32, writable=False, shape=(*), order='C'], nmm_interpolation_data: tuple[ndarray[dtype=int32, writable=False, shape=(*), order='C'], ndarray[dtype=int32, writable=False, shape=(*), order='C'], ndarray[dtype=float64, writable=False, shape=(*), order='C'], ndarray[dtype=int32, writable=False, shape=(*), order='C']]) -> None
    4. interpolate(self, expr: dolfinx.cpp.fem.Expression_float64, cells: ndarray[dtype=int32, writable=False, order='C'], expr_mesh: dolfinx.cpp.mesh.Mesh_float64, cell_map: ndarray[dtype=int32, writable=False, order='C']) -> None

Invoked with types: dolfinx.cpp.fem.Function_float64, NoneType, ndarray, dolfinx.fem.function.PointOwnershipData

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/mnt/c/Work/FEniCSx/plasticity1.py", line 557, in <module>
    p_avg.interpolate(project(p, P0))
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py", line 495, in interpolate
    assert callable(u)
AssertionError

Please note that there is so much unrelated code in your script, that I do not have the bandwidth to reduce it (there are over 350 lines of code), and you are asking about just the project function. Please make an example only projecting a function p into P0, rather than setting up a custom mesh, with a related PDE.