Unexpected stress distribution on domain with two materials


I am attempting to solve the following problem:

A 2D beam is clamped on its left end subjected a tensile load in positive x-direction on its right end. An inner square in the middle of the beam consists of a stiffer material.

The displacement solution is in line with expectation and other finite element software. The stress distribution however shows some unexpected behaviour. Below is a plot of the normal stress in x-dimension. The lower material boundary behaves as expected, with the inner material showing significantly higher stress values due to its higher stiffness. However, the upper boundary does not reflect this trend.

I am not sure what the reason for this is, but I would suspect that it has something to do with the interpolation of the stresses? Below is my code to reproduce the plot above. I appreciate any help.

import numpy as np
import ufl
import matplotlib.pyplot as plt

from dolfinx import fem, mesh, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI

# structural parameters
Lx = 2
Ly = 1
# number of elements per unit length in each dimension
N = 100
# lower left and upper right corner of inner square
ll = [0.75, 0.25]
ur = [1.25, 0.75]
# force on right end
f_right = (0.5, 0.0)
tol = 0
element_degree = 1
quadrature_degree = 4

# material parameters
matmod = 'Hooke'
E1 = default_scalar_type(1.)
E2 = default_scalar_type(10.)
nu = default_scalar_type(0.3)

# mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [Lx, Ly]], [Lx*N, Ly*N], mesh.CellType.quadrilateral)
Q = fem.functionspace(domain, ("DG", 0))
V = fem.functionspace(domain, ("CG", element_degree, (domain.geometry.dim, )))

# left end of beam
def left(x):
    return np.isclose(x[0], 0)

# right end of beam
def right(x):
    return np.isclose(x[0], Lx)

# outer material
def domain1(x):
    return np.logical_or(np.logical_or(x[0] <= ll[0] + tol, x[0] >= ur[0] - tol), np.logical_or(x[1] <= ll[1] + tol, x[1] >= ur[1] - tol))

# inner material
def domain2(x):
    return np.logical_and(np.logical_and(ll[0] - tol <= x[0], x[0] <= ur[0] + tol), np.logical_and(ll[1] - tol <= x[1], x[1] <= ur[1] + tol))

# mark cells for bcs
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, left)
right_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag_bc = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])

# define mu and lambda differently for each material
cells1 = mesh.locate_entities(domain, domain.topology.dim, domain1)
cells2 = mesh.locate_entities(domain, domain.topology.dim, domain2)
mu = fem.Function(Q)
lam = fem.Function(Q)
mu.x.array[cells1] = np.full_like(cells1, E1 / (2 * (1 + nu)), dtype=default_scalar_type)
mu.x.array[cells2] = np.full_like(cells2, E2 / (2 * (1 + nu)), dtype=default_scalar_type)
lam.x.array[cells1] = np.full_like(cells1, E1 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
lam.x.array[cells2] = np.full_like(cells2, E2 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)

# dirchlet bc
u_zero = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
dirichlet_bc = [fem.dirichletbc(u_zero, left, V)]

# test and solution functions
v = ufl.TestFunction(V)
u = fem.Function(V)
# traction
t = fem.Constant(domain, default_scalar_type(f_right))

# spatial dimension
d = len(u)
# identity tensor
I = ufl.variable(ufl.Identity(d))
# lagrange strain
eps = ufl.sym(ufl.grad(u))

# strain energy
sig = 2 * mu * eps + lam * ufl.tr(eps) * I
psi = 0.5 * ufl.inner(eps, sig)

# set gauß quadrature degree
metadata = {'quadrature_degree': quadrature_degree}
# measure inner domain
dx = ufl.Measure('dx', domain=domain, metadata=metadata)
# measure right end
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_bc, metadata=metadata)

# potential energy (integrate strain energy on inner and outer material, traction energy on right boundary)
Pi = psi*dx - ufl.dot(t, u)*ds(2)
# variation of Pi
dPidu = ufl.derivative(Pi, u, v)

# define problem F(u) = 0
problem = NonlinearProblem(dPidu, u, dirichlet_bc)
# set solver and options
solver = NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

# solve

# get stress
stress = 2 * mu * eps + lam * ufl.tr(eps) * I

# interpolate normal stress in x
V = fem.functionspace(domain, ("CG", element_degree))
sxx_expr = fem.Expression(stress[0, 0], V.element.interpolation_points())
sxx = fem.Function(V)

# plot normal stress in x
x = np.array(V.tabulate_dof_coordinates())
plt.scatter(x[:, 0], x[:, 1], c=sxx.vector.array, marker='s')

You probably want to interpolate your stress over a DG space rather than CG
Taking your code and replacing

V = fem.functionspace(domain, ("DG", 0))

yields the plot that, I think, you are trying to get:

Hope this helps,


Thank you for your help! Could you also tell me how to interpolate the displacement over a DG space? It doesn’t seem as intuitive because I cannot index the DOFs of the displacement solution separately as far as I know.

I’m not sure I understand where you have a problem. You can interpolate in the same way the displacement over a vector DG space, or interpolate one of its component on the scalar DG space by indexing (u[0] and u[1]).

This works, thank you. Sorry for the basic questions, I am quite new to Fenics.

If you have the time, maybe you can help me more. To give some context, I am using Fenics to set up an FEM reference solution to judge the quality of a meshfree method I am investigating. In order to calculate the error in each collocation point, it would be convenient to have the FEM solution on the same coordinates as the collocation points. From my understanding, using your approach, I obtain values for stress and displacement in the middle of each element. It would be better to interpolate between nodes to assure the possibility of the FEM evaluation points and colocation points sharing the same coordinates. However, if I use the DG space with degree 1 to interpolate between nodes, I get the same problem with the stress field as before.

I see what you are trying to achieve. But your stress is discontinuous at the interface between the two materials, right, so if you have an interpolation point at this interface, what should be the value of the stress at that point? Maybe this is the core of the issue.

The FEM calculates stresses in the integration points of each element right? So I thought in each element the stress is interpolated to the nodes using the shape functions, and then each node takes the average of all neighbouring elements interpolating to that node. I do see that discontinuous stress can pose an issue here, but I don’t understand how it leads to the asymmetric distribution in the plot in my first post.

Stress is indeed calculated at the integration points when computing the integral over the domain. However, as far as I understand you are trying to get values at the nodes (that’s what you get by using a projection on ("CG",1)). The problem is that stress is not well defined at the nodes, because the value at the node is dependent on the cell.

That is not how standard interpolation works. Standard interpolation of a finite element function only happens for a coordinate within a specific cell. You evaluate the stress at the interpolation points of a given element to get the stress from the point of view of this cell.

You could project the stress into a continuous space, but you would suffer from Gibbs phenomenon, as described in: Projection and interpolation — FEniCS Tutorial @ Sorbonne

Thank for the info, I wasn’t aware of that. In that case, I think it makes the most sense for me to just project into a discontinuous space as suggested by Maxime and evaluate the meshfree method accordingly.
One more question remains for me: I don’t understand why the stress on the left Dirichlet boundary shows asymmetric behaviour. In the upper corner, two elements are affected by high stress concentration, while in the lower corner it is only one. This also occurs on a homogeneous domain (reproducible by setting E2=E1 in my code).

Please provide a minimal reproducible example, including the plotting script that exhibits this behavior. In general, i would use matplotlib for plotting stresses, but pyvista or paraview, as it has support for CG/DG elements on unstructured grids

Thanks for the hint about the plotting. The problem was indeed simply due to an unsuitable marker size in the scatter plot.