Proper Implemention of KKT conditions in energy minimization

[Joshua_Simpson-Gomez]

I am trying to simulate a linear elastic problem with a small amount of contact. To get the best results I was hoping to implement a Lagrange multiplier approach to solve for the displacements with constraints on the contact gap. When I read literature I typically find it is already in the matrix form and I am unsure of the best way to convert it into the variational form to solve in fenicsx.

Given a displacement u, a contact pressure K, and the normal to each node on the contact surface n the signorini contact conditions are as follows:

K\geq 0, \ u \cdot n - g \leq 0, \ K(u \cdot n - g) = 0

which also correspond to the KKT conditions where P is the lagrange multiplier and u \cdot n - g \leq 0 is the inequality constraint. The energy minimization is linear elastic and minimizes

E = \int \sigma (u) : \epsilon (u) \ d\Omega \

with u = u_{applied} on \Gamma_{u} with \sigma being the stress and \epsilon being the strain. I tried to follow the tutorial for constrained minimization of the poisson equation, but since the gap function is not a PDE it doesnt seem to form a linear problem when solving the adjoint equation, therefore I tried to derive the variational form as follows:

min \ E(u) \ s.t. \ G(u) = 0

with G(u) being the gap function seen as

G(u) = \int <u \cdot n - g>_{+} d\Gamma_{contact}

and <x>_{+} being x when x is positive and 0 when x is negative, automatically satisfying the inequality constraint. To satisfy the equality constraint that G(u) = 0, the previous and following equations are solved using a block solver for both u and k being a CG 2d vector and CG scalar values respectively.

\int \frac{\partial E}{\partial u}dx + \int K \frac{\partial G}{\partial u}d\Gamma_{contact}

The minimum working example code is included below written on dolfinx 0.10.0 installed on docker.

from mpi4py import MPI

import gmsh
import numpy as np

import dolfinx.fem.petsc
import ufl

from dolfinx import mesh

from BoltGroupComputer import BoltGroup

import pyvista

plotstress = True

# physical variables
N = 25
L = 1.0
D = 2.0*L

# modifying the applied load changes the max error (tested at finite difference step of 1e-6)
u_disp = 0.1

# material properties
E_val = 200.0
nu_val = 0.3
iniK = 50.0

# values for loading and contact
R = .5
d = 0.0

# create mesh
comm = MPI.COMM_WORLD
domain = mesh.create_rectangle(
    comm,
    [np.array([-L, -D]), np.array([L, 0])],
    [N, int(round(N*D/L))],
    cell_type=mesh.CellType.quadrilateral,
)

# locate all exterior faces (for later boolean operations)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim
boundary_facets_all = dolfinx.mesh.exterior_facet_indices(domain.topology)

# define the mesh marker for each facets on the mesh
num_facets = domain.topology.index_map(fdim).size_local
markers = np.zeros(num_facets, dtype=np.int32)
top_mark,bottom_mark,free_mark = 1,2,3

def topBnds(x):
    return np.isclose(x[1],0)

def bottomBnds(x):
    return np.isclose(x[1],-D)

# locate different regions of facets where the boundary conditions will be applied
top_facets, markers[top_facets] = mesh.locate_entities_boundary(domain, fdim, topBnds), top_mark
bottom_facets, markers[bottom_facets] = mesh.locate_entities_boundary(domain, fdim, bottomBnds), bottom_mark
free_facets, markers[free_facets] = np.setdiff1d(boundary_facets_all,np.union1d(top_facets,bottom_facets)), free_mark

facet_tags = mesh.meshtags(domain, fdim, np.arange(num_facets, dtype=np.int32), markers)
subdomain_tags = np.hstack(top_facets)

gamma, gamma_to_domain = dolfinx.mesh.create_submesh(domain, fdim, subdomain_tags)[0:2]

V = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
Q = dolfinx.fem.functionspace(gamma, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, Q)

dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) 

E = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(E_val))
nu = dolfinx.fem.Constant(domain, nu_val)

x = ufl.SpatialCoordinate(domain)
f = dolfinx.fem.Constant(domain, np.zeros(gdim, dtype=dolfinx.default_scalar_type))

v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
k = dolfinx.fem.Function(Q, name="lambda")
k.x.array[:] = iniK

mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(w):
    return ufl.sym(ufl.grad(w))

def sigma(w, mu, lmbda):
    ew = epsilon(w)
    gdim = ew.ufl_shape[0]
    return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)

F = ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= ufl.inner(f, v) * dx # body force

# define the obstacle norm and gap functions
centroid = ufl.as_vector([0.0,R-d])
x = ufl.SpatialCoordinate(domain)
dist = ufl.sqrt(ufl.dot(centroid-x,centroid-x))
gap = dist-R
norm = (centroid-x)/dist
deformedGap = ufl.dot(norm,u)-gap

# derived equations
G = ufl.conditional(ufl.ge(deformedGap,0),deformedGap,0)
dGdu = ufl.derivative(G,u,v)
F += k*dGdu*ds(top_mark)
F += ufl.inner(G,w)*ds(top_mark)

residual = ufl.extract_blocks(F)

du, dpsi = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, k, dpsi)
J = ufl.extract_blocks(jac)


bc_bottom = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((u_disp))),
                    dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(bottom_mark)), V.sub(1))
bc_fixed = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((0))),
                    dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(bottom_mark)), V.sub(0))
bcs = [bc_bottom, bc_fixed]

V_DG = dolfinx.fem.functionspace(domain, ("DG", 1, (domain.geometry.dim,)))
stress_space, stress_to_disp = V_DG.sub(0).collapse()
von_mises = dolfinx.fem.Function(stress_space, name="von_Mises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, stress_space.element.interpolation_points)

entity_maps = [gamma_to_domain]
solver = dolfinx.fem.petsc.NonlinearProblem(
    residual,
    u=[u, k],
    J=J,
    bcs=bcs,
    entity_maps=entity_maps,
    petsc_options={
        "snes_monitor": None,
        "ksp_type": "preonly",
        "pc_type": "hypre",
        "ksp_error_if_not_converged": True,
        # "snes_error_if_not_converged": True,
        "snes_max_it": 1000,
    },
    petsc_options_prefix="signorini_",
)

solver.solve()

von_mises.interpolate(stress_expr)
u_dg.interpolate(u)

print("max deformation: " + str(np.max(np.abs(u.x.array))))
print("avg pressure value: " + str(np.mean(k.x.array)) + " and std: " + str(np.std(k.x.array)))

if (plotstress):
    grid = dolfinx.plot.vtk_mesh(u_dg.function_space)
    pyvista_grid = pyvista.UnstructuredGrid(*grid)
    values = u_dg.x.array.reshape(-1, domain.geometry.dim)
    values_padded = np.zeros((values.shape[0], 3))
    values_padded[:, : domain.geometry.dim] = values
    pyvista_grid.point_data["u"] = values_padded
    warped = pyvista_grid.warp_by_vector("u")
    stresses = np.zeros_like(u_dg.x.array)
    stresses[stress_to_disp] = von_mises.x.array
    stresses = stresses.reshape(-1, domain.geometry.dim)[:, 0]
    warped.point_data["von_Mises"] = stresses

    plotter = pyvista.Plotter()
    plotter.add_mesh(pyvista_grid, style="wireframe", color="b")
    plotter.add_mesh(warped, scalars="von_Mises", lighting=True, show_edges=True)
    plotter.show_bounds()
    plotter.view_xy()
    plotter.show()

After running the simulation, the displacement is incorrect, but has a correct trend. One of the main concerns is during solving the value of K is not modified from the initial guess, and only u changes, making it behave like a worse performing penalty method but I’m not sure why since 2 equations are solved for. I think it has to do with how the inequality constraint is modified, but Im not sure how else to do it. In the problem, the top surface is the contact surface and has a circular indenter fixed above the domain. the bottom of the domain has a displacement applied of 0.1 upwards to make it contact the circular indenter (hertzian contact). The correct solution is included below.

Before trying the lagrange multiplier method I implemented a penalty approach and the latent variable approach from the fenicsx examples, and while they both work well at larger displacements, it seems like the smaller contact distances in the problem I am simulating cause convergence issues and have less accurate results than desired. This should be solved by the lagrange multiplier method since that directly solves for the contact pressure, but I just dont know how to formulate it, any help or tips would be greatly appreciated!

Thank you!

[dokken]

@Joshua_Simpson-Gomez

latent variable approach from the fenicsx examples, and while they both work well at larger displacements, it seems like the smaller contact distances in the problem I am simulating cause convergence issues and have less accurate results than desired.

It surprises me that LVPP gives non optimal results for small deformation. Do you have the code for that?

Anyhow, for your current approach the only thing I spot at a quick glance is that you solve with preonly and hypre.
Have you tried preonly and lu (direct solver)?

I am also uncertain if you can use the conditional to enforce the variational inequality for the Lagrange multiplier, ref: 10.1007/s00211-009-0273-z

[Joshua_Simpson-Gomez]

Thank you for your quick reply! For the LVPP, on a second look into the code and you are right, It can give good results, I had the value of alpha too low but now it converges below 1e-10. The other reason I wanted to use a lagrange multiplier approach is that I am trying to implement topology optimization with contact, and it seems like the LVPP method would require a history dependent sensitivity analysis whereas the lagrange method would not.

I have tried lu and it gives me this error

I also tried solving the first and second equations sequentially and for that I had to use hypre to solve for G(u) = 0, but could use lu for solving the second equation above involving the derivatives of E and G with respect to u, however it resulted in very similar results to the above code.

thank you for the source, I think you are right that applying the conditional is not correct, looking at their formulation when solving they have:

in which the second equation is greater than or equal to 0. To try and impose this in fenicsx I added a slack variable so that there are 3 equations to solve:

G = \int u \cdot n - g + m^{2} \ d\Gamma_{contact} = 0

\int \frac{\partial E}{\partial u} d\Omega + \int K\frac{\partial G}{\partial u} d\Gamma_{contact}

\int Km \ d\Gamma_{contact} = 0

and while this runs in fenicsx using the hypre solver (gives same error as before using lu), the values of K still dont change as well as the values of m from the initial values. I think it might have to do with the integral enforcing the complementary slackness; is using a slack variable and solving for u, K, and m a valid approach?

I included the code with the slack variable below in case it is useful, thank you!

from mpi4py import MPI

import gmsh
import numpy as np

import dolfinx.fem.petsc
import ufl

from dolfinx import mesh

from BoltGroupComputer import BoltGroup

import pyvista

plotstress = True

# physical variables
N = 25
L = 1.0
D = 2.0*L

# modifying the applied load changes the max error (tested at finite difference step of 1e-6)
u_disp = 0.1

# material properties
E_val = 200.0
nu_val = 0.3
iniK = 50.0

# values for loading and contact
R = .5
d = 0.0

# create mesh
comm = MPI.COMM_WORLD
domain = mesh.create_rectangle(
    comm,
    [np.array([-L, -D]), np.array([L, 0])],
    [N, int(round(N*D/L))],
    cell_type=mesh.CellType.quadrilateral,
)

# locate all exterior faces (for later boolean operations)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim
boundary_facets_all = dolfinx.mesh.exterior_facet_indices(domain.topology)

# define the mesh marker for each facets on the mesh
num_facets = domain.topology.index_map(fdim).size_local
markers = np.zeros(num_facets, dtype=np.int32)
top_mark,bottom_mark,free_mark = 1,2,3

def topBnds(x):
    return np.isclose(x[1],0)

def bottomBnds(x):
    return np.isclose(x[1],-D)

# locate different regions of facets where the boundary conditions will be applied
top_facets, markers[top_facets] = mesh.locate_entities_boundary(domain, fdim, topBnds), top_mark
bottom_facets, markers[bottom_facets] = mesh.locate_entities_boundary(domain, fdim, bottomBnds), bottom_mark
free_facets, markers[free_facets] = np.setdiff1d(boundary_facets_all,np.union1d(top_facets,bottom_facets)), free_mark

facet_tags = mesh.meshtags(domain, fdim, np.arange(num_facets, dtype=np.int32), markers)
subdomain_tags = np.hstack(top_facets)

gamma, gamma_to_domain = dolfinx.mesh.create_submesh(domain, fdim, subdomain_tags)[0:2]

V = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
Q = dolfinx.fem.functionspace(gamma, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, Q, Q)

dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) 

E = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(E_val))
nu = dolfinx.fem.Constant(domain, nu_val)

x = ufl.SpatialCoordinate(domain)
f = dolfinx.fem.Constant(domain, np.zeros(gdim, dtype=dolfinx.default_scalar_type))

v, w, z = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
k = dolfinx.fem.Function(Q, name="lambda")
m = dolfinx.fem.Function(Q, name="slack")
k.x.array[:] = iniK

mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(w):
    return ufl.sym(ufl.grad(w))

def sigma(w, mu, lmbda):
    ew = epsilon(w)
    gdim = ew.ufl_shape[0]
    return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)

F = ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= ufl.inner(f, v) * dx # body force

# define the obstacle norm and gap functions
centroid = ufl.as_vector([0.0,R-d])
x = ufl.SpatialCoordinate(domain)
dist = ufl.sqrt(ufl.dot(centroid-x,centroid-x))
gap = dist-R
norm = (centroid-x)/dist
deformedGap = ufl.dot(norm,u)-gap

# derived equations
G = deformedGap+m**2
dGdu = ufl.derivative(G,u,v)
F += k*dGdu*ds(top_mark)
F += ufl.inner(G,w)*ds(top_mark)
F += ufl.inner(ufl.inner(k,m),z)*ds(top_mark)

residual = ufl.extract_blocks(F)

du, dk, dm = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, k, dk) + ufl.derivative(F, m, dm)
J = ufl.extract_blocks(jac)


bc_bottom = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((u_disp))),
                    dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(bottom_mark)), V.sub(1))
bc_fixed = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((0))),
                    dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(bottom_mark)), V.sub(0))
bcs = [bc_bottom, bc_fixed]

V_DG = dolfinx.fem.functionspace(domain, ("DG", 1, (domain.geometry.dim,)))
stress_space, stress_to_disp = V_DG.sub(0).collapse()
von_mises = dolfinx.fem.Function(stress_space, name="von_Mises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, stress_space.element.interpolation_points)

entity_maps = [gamma_to_domain]
solver = dolfinx.fem.petsc.NonlinearProblem(
    residual,
    u=[u, k, m],
    J=J,
    bcs=bcs,
    entity_maps=entity_maps,
    petsc_options={
        "snes_monitor": None,
        "ksp_type": "preonly",
        "pc_type": "hypre",
        "ksp_error_if_not_converged": True,
        # "snes_error_if_not_converged": True,
        "snes_max_it": 1000,
    },
    petsc_options_prefix="signorini_",
)

solver.solve()

von_mises.interpolate(stress_expr)
u_dg.interpolate(u)

print("max deformation: " + str(np.max(np.abs(u.x.array))))
print("avg pressure value: " + str(np.mean(k.x.array)) + " and std: " + str(np.std(k.x.array)))
print("avg slack value: " + str(np.mean(m.x.array)) + " and std: " + str(np.std(m.x.array)))

if (plotstress):
    grid = dolfinx.plot.vtk_mesh(u_dg.function_space)
    pyvista_grid = pyvista.UnstructuredGrid(*grid)
    values = u_dg.x.array.reshape(-1, domain.geometry.dim)
    values_padded = np.zeros((values.shape[0], 3))
    values_padded[:, : domain.geometry.dim] = values
    pyvista_grid.point_data["u"] = values_padded
    warped = pyvista_grid.warp_by_vector("u")
    stresses = np.zeros_like(u_dg.x.array)
    stresses[stress_to_disp] = von_mises.x.array
    stresses = stresses.reshape(-1, domain.geometry.dim)[:, 0]
    warped.point_data["von_Mises"] = stresses

    plotter = pyvista.Plotter()
    plotter.add_mesh(pyvista_grid, style="wireframe", color="b")
    plotter.add_mesh(warped, scalars="von_Mises", lighting=True, show_edges=True)
    plotter.show_bounds()
    plotter.view_xy()
    plotter.show()

[Joshua_Simpson-Gomez]

I also tried implementing the approach from the paper

@dokken

ref: 10.1007/s00211-009-0273-z

directly, but it doesnt seem to compile, can equation 12 from the source (picture below) be directly implemented in fenicsx?

I know that having a >= condition in equation 2 would need special consideration, but I thought that while it wouldnt converge properly, the form in the second equation should compile since its the same as the paper, but it gives me the following error when I try it:

Are there any examples of how to solve a mixed FEA approach where one of the equations is >= 0 like the second equation?

I included the code for the approach from the paper in case the error is something unrelated to the variational formulation. I also tried solving it iteratively and while the first equation in line 1 of equation 12 compiles and solves, the second nonlinearproblem has the same error as the block solver which makes me think it is something with the formulation and not the code.

from mpi4py import MPI

import gmsh
import numpy as np

import dolfinx.fem.petsc
import ufl

from dolfinx import mesh

from BoltGroupComputer import BoltGroup

import pyvista

plotstress = True

# physical variables
N = 25
L = 1.0
D = 2.0*L

# modifying the applied load changes the max error (tested at finite difference step of 1e-6)
u_disp = 0.1

# material properties
E_val = 200.0
nu_val = 0.3
iniK = 50.0

# values for loading and contact
R = .5
d = 0.0

# create mesh
comm = MPI.COMM_WORLD
domain = mesh.create_rectangle(
    comm,
    [np.array([-L, -D]), np.array([L, 0])],
    [N, int(round(N*D/L))],
    cell_type=mesh.CellType.quadrilateral,
)

# locate all exterior faces (for later boolean operations)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
gdim = domain.geometry.dim
boundary_facets_all = dolfinx.mesh.exterior_facet_indices(domain.topology)

# define the mesh marker for each facets on the mesh
num_facets = domain.topology.index_map(fdim).size_local
markers = np.zeros(num_facets, dtype=np.int32)
top_mark,bottom_mark,free_mark = 1,2,3

def topBnds(x):
    return np.isclose(x[1],0)

def bottomBnds(x):
    return np.isclose(x[1],-D)

# locate different regions of facets where the boundary conditions will be applied
top_facets, markers[top_facets] = mesh.locate_entities_boundary(domain, fdim, topBnds), top_mark
bottom_facets, markers[bottom_facets] = mesh.locate_entities_boundary(domain, fdim, bottomBnds), bottom_mark
free_facets, markers[free_facets] = np.setdiff1d(boundary_facets_all,np.union1d(top_facets,bottom_facets)), free_mark

facet_tags = mesh.meshtags(domain, fdim, np.arange(num_facets, dtype=np.int32), markers)
subdomain_tags = np.hstack(top_facets)

gamma, gamma_to_domain = dolfinx.mesh.create_submesh(domain, fdim, subdomain_tags)[0:2]

V = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
Q = dolfinx.fem.functionspace(gamma, ("CG", 1))
W = ufl.MixedFunctionSpace(V, Q)

dx = ufl.Measure("dx", domain=domain)
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) 

E = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(E_val))
nu = dolfinx.fem.Constant(domain, nu_val)

x = ufl.SpatialCoordinate(domain)
f = dolfinx.fem.Constant(domain, np.zeros(gdim, dtype=dolfinx.default_scalar_type))

v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
k = dolfinx.fem.Function(Q, name="lambda")
k.x.array[:] = iniK

mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(w):
    return ufl.sym(ufl.grad(w))

def sigma(w, mu, lmbda):
    ew = epsilon(w)
    gdim = ew.ufl_shape[0]
    return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)

def sigmaN(n,w,mu,lmbda):
    return ufl.dot(n,ufl.dot(sigma(w,mu,lmbda),n))

F = ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= ufl.inner(f, v) * dx # body force

# define the obstacle norm and gap functions
centroid = ufl.as_vector([0.0,R-d])
x = ufl.SpatialCoordinate(domain)
dist = ufl.sqrt(ufl.dot(centroid-x,centroid-x))
gap = dist-R
norm = (centroid-x)/dist

# # DIRECT EQUATIONS FROM EQ12
alpha = dolfinx.fem.Constant(domain, 0.1)
F -= k*ufl.dot(v,norm)*ds(top_mark)
F += alpha*(k-sigmaN(norm,u,mu,lmbda))*sigmaN(norm,v,mu,lmbda)*ds(top_mark)
# second equation
F += (w-k)*ufl.dot(u,norm)*ds(top_mark)
F += alpha*(w-k)*(k-sigmaN(norm,u,mu,lmbda))*ds(top_mark)

residual = ufl.extract_blocks(F)

du, dk = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, k, dk)
J = ufl.extract_blocks(jac)


bc_bottom = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((u_disp))),
                    dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, facet_tags.find(bottom_mark)), V.sub(1))
bc_fixed = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(domain, dolfinx.default_scalar_type((0))),
                    dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, facet_tags.find(bottom_mark)), V.sub(0))
bcs = [bc_bottom, bc_fixed]

V_DG = dolfinx.fem.functionspace(domain, ("DG", 1, (domain.geometry.dim,)))
stress_space, stress_to_disp = V_DG.sub(0).collapse()
von_mises = dolfinx.fem.Function(stress_space, name="von_Mises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, stress_space.element.interpolation_points)

entity_maps = [gamma_to_domain]
solver = dolfinx.fem.petsc.NonlinearProblem(
    residual,
    u=[u, k],
    J=J,
    bcs=bcs,
    entity_maps=entity_maps,
    petsc_options={
        "snes_monitor": None,
        "ksp_type": "preonly",
        "pc_type": "lu",
        "ksp_error_if_not_converged": True,
        # "snes_error_if_not_converged": True,
        "snes_max_it": 1000,
    },
    petsc_options_prefix="signorini_",
)

solver.solve()

von_mises.interpolate(stress_expr)
u_dg.interpolate(u)

print("max deformation: " + str(np.max(np.abs(u.x.array))))
print("avg pressure value: " + str(np.mean(k.x.array)) + " and std: " + str(np.std(k.x.array)))

if (plotstress):
    grid = dolfinx.plot.vtk_mesh(u_dg.function_space)
    pyvista_grid = pyvista.UnstructuredGrid(*grid)
    values = u_dg.x.array.reshape(-1, domain.geometry.dim)
    values_padded = np.zeros((values.shape[0], 3))
    values_padded[:, : domain.geometry.dim] = values
    pyvista_grid.point_data["u"] = values_padded
    warped = pyvista_grid.warp_by_vector("u")
    stresses = np.zeros_like(u_dg.x.array)
    stresses[stress_to_disp] = von_mises.x.array
    stresses = stresses.reshape(-1, domain.geometry.dim)[:, 0]
    warped.point_data["von_Mises"] = stresses

    plotter = pyvista.Plotter()
    plotter.add_mesh(pyvista_grid, style="wireframe", color="b")
    plotter.add_mesh(warped, scalars="von_Mises", lighting=True, show_edges=True)
    plotter.show_bounds()
    plotter.view_xy()
    plotter.show()

Thank you!

[Stein]

The second line is not linear in the test function. Your (w-k)*… term (with w a test function and k some function) leads to a condition that is impossible to satisfy (both as an equality and as an inequality) if k is not linearly dependent on w. In terms of FEniCs language, FEniCS cannot even assemble a residual vector if not all terms are linear in w.

I didn’t look at the paper, but I’d expect they have an explicit expression for k.

[Joshua_Simpson-Gomez]

Thank you for the reply! I was also uncertain about that, in the formulation from the paper it seems like \mu (w in the code) is the test function while \lambda (k in the code) is the second value in the formulation which is solved for, where the spaces are


and it is also mentioned that the mixed formulation of the contact problem is finding

such that

and given they are only solving for \lambda and u, and not \mu I thought \mu had to be the test function, but I am not sure, is \mu something else or a value to be solved for as well?

I am especially confused because they say the variational form

can then be used to compile the following matrix form,

so given that fenics cant assemble the residual vector if not all terms are linear in \lambda as you mentioned and it seems like the above equation would lead to that result, would the best approach be to try and compile each individual form into the matrix form and solve that instead of solving the variational form directly using fenicsx?

[Stein]

Hmmm, maybe I jumped the gun a bit. I now see that the space for \mu is not a linear space, so then my point about inability to satisfy maybe doesn’t quite apply.

But then this cannot be assembled into vectors and matrices, as they are not linear forms (my point about FEniCS not being able to do that wasn’t meant as a comment about FEniCS; this simply cannot be done). I suppose you could manually assemble K C F L B D U and M from (35), but note that \bar{L} is a vector of unknowns that are only permitted to be positive. So it doesn’t have a vector assembled form. The path forward to a solution scheme is thus obscured to me. One would have to dig into that paper a bit deeper to understand how they actually solve that system then.

B.t.w, are you aware of this nice demo? Coupling PDEs of multiple dimensions — FEniCS Workshop

[Joshua_Simpson-Gomez]

Thank you for the response, I will look into assembling the vectors manually to solve. Since the end goal is to do topology optimization, I was hoping to implement the lagrange multiplier approach since its more accurate than a penalty approach but does not require a history dependent sensitivity analysis needed for latent variable approach. It is a good method if I didnt have to do the sensitivity analysis.