Using dolfinx_mpc in generalised plain strain formulation

Gaurav_Verma · July 3, 2025, 10:41pm

Hello, I wish to implement create_slip_constraint from dolfinx_mpc on the generalized plane strain case. The left and bottom edges are allowed to slide tangentially but not normally.

I am stuck in the part dolfinx.fem.petsc.assemble_matrix_block implementation with mpc. A guidance will be quite helpful. Here is my code below.

import gmsh
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import fem, io, plot
import dolfinx.fem.petsc
from dolfinx.cpp.la.petsc import get_local_vectors
from scifem import create_real_functionspace
from dolfinx_mpc.utils import (create_point_to_point_constraint, determine_closest_block, rigid_motions_nullspace, facet_normal_approximation)
from dolfinx_mpc import LinearProblem, MultiPointConstraint

hsize = 0.2

Re = 11.0
Ri = 9.0

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
gdim = 2
model_rank = 0
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, [outer_circ], 3, name="outer")


gmsh.model.mesh.generate(gdim)

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

def eps(v, ezz):
    return ufl.sym(
        ufl.as_tensor(
            [
                [v[0].dx(0), v[0].dx(1), 0],
                [v[1].dx(0), v[1].dx(1), 0],
                [0, 0, ezz],
            ]
        )
    )


E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)


def sigma(v, ezz):
    return lmbda * ufl.tr(eps(v, ezz)) * ufl.Identity(3) + 2.0 * mu * eps(v, ezz)


def eps(v, ezz):
    return ufl.sym(
        ufl.as_tensor(
            [
                [v[0].dx(0), v[0].dx(1), 0],
                [v[1].dx(0), v[1].dx(1), 0],
                [0, 0, ezz],
            ]
        )
    )


E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)


def sigma(v, ezz):
    return lmbda * ufl.tr(eps(v, ezz)) * ufl.Identity(3) + 2.0 * mu * eps(v, ezz)


# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))
R = create_real_functionspace(domain)
W = MixedFunctionSpace(V, R)

# Define Variational Problem
du, dezz = TrialFunctions(W)
u_, ezz_ = TestFunctions(W)


#u = fem.Function(V, name = "Displacement")
a_form = inner(sigma(du, dezz), eps(u_, ezz_)) * dx


T = fem.Constant(domain, 1000000.0)


#Self-weight on the surface
n = FacetNormal(domain)

ds = Measure("ds", domain=domain, subdomain_data=facet_markers)
dx = Measure("dx", domain=domain, subdomain_data=cell_markers)
dS = Measure("dS", domain=domain, subdomain_data=facet_markers)

L_form = dot(T*n,u_) * ds(3)

a_blocked_compiled = fem.form(extract_blocks(a_form))
L_blocked_compiled = fem.form(extract_blocks(L_form))


# Get the facet markers for boundary edges (ps1 and ps2)

mpc = MultiPointConstraint(V)

for tag in [1, 2]:  # ps1 and ps2
    normal_vec = create_normal_approximation(V, facet_markers, tag)
    mpc.create_slip_constraint(V, (facet_markers, tag), normal_vec)

mpc.finalize()



a_blocked_compiled = fem.form(ufl.extract_blocks(a_form))
L_blocked_compiled = fem.form(ufl.extract_blocks(L_form))

after this part I am stuck and need your guidance.

A = dolfinx.fem.petsc.assemble_matrix_block(a_blocked_compiled, bcs=bcs)
A.assemble()
b = dolfinx.fem.petsc.assemble_vector_block(
    L_blocked_compiled, a_blocked_compiled, bcs=bcs
)

ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")

# solve system
xh = dolfinx.fem.petsc.create_vector_block(L_blocked_compiled)
ksp.solve(b, xh)
xh.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

maps = [(Wi.dofmap.index_map, Wi.dofmap.index_map_bs) for Wi in W.ufl_sub_spaces()]
u = fem.Function(V, name="Displacement")
x_local = get_local_vectors(xh, maps)
u.x.array[: len(x_local[0])] = x_local[0]
u.x.scatter_forward()
ezz = x_local[1][0]

dokken · July 4, 2025, 11:27am

Block-assembly is not implemented in DOLFINx_MPC, only nest assembly:

github.com/jorgensd/dolfinx_mpc

python/demos/demo_stokes_nest.py

v0.9.2

# Copyright (C) 2022 Nathan Sime
#
# This file is part of DOLFINX_MPC
#
# SPDX-License-Identifier:    MIT
#
# This demo illustrates how to apply a slip condition on an
# interface not aligned with the coordiante axis.
# The demos solves the Stokes problem using the nest functionality to
# avoid using mixed function spaces. The demo also illustrates how to use
#  block preconditioners with PETSc
from __future__ import annotations

from pathlib import Path

from mpi4py import MPI
from petsc4py import PETSc

import basix
import dolfinx.io

This file has been truncated. show original

which can use direct solvers, as illustrated in:

github.com/FEniCS/dolfinx

Deprecate block matrices/kind=MPI for PETSc

opened 12:51PM - 13 May 25 UTC

jorgensd

proposal

### Describe new/missing feature The main difference between PETSc NEST and PET…Sc block matrices **seems** to be the location of the ghosts (blocked matrix puts all ghosts at the end, while nest matrix puts them at the end of each block). @schnellerhase , @michalhabera and I tried to use a direct solver with a nest matrix, and it seems to work. Is there any performance study that indicates that our blocked matrices has a performance benefit compared to the nest matrices? If so, we should of course keep them. ```python from packaging.version import Version from mpi4py import MPI from petsc4py import PETSc from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors import dolfinx.fem.petsc import numpy as np from scifem import create_real_functionspace, assemble_scalar from scifem.petsc import apply_lifting_and_set_bc, zero_petsc_vector import ufl M = 40 mesh = dolfinx.mesh.create_unit_square( MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=np.float64 ) V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) def u_exact(x): return ufl.sin(2 * ufl.pi * x[0]) x = ufl.SpatialCoordinate(mesh) n = ufl.FacetNormal(mesh) g = ufl.dot(ufl.grad(u_exact(x)), n) f = -ufl.div(ufl.grad(u_exact(x))) h = assemble_scalar(u_exact(x) * ufl.dx) # We then create the Lagrange multiplier space R = create_real_functionspace(mesh) # Next, we can create a mixed-function space for our problem W = ufl.MixedFunctionSpace(V, R) u, lmbda = ufl.TrialFunctions(W) du, dl = ufl.TestFunctions(W) # We can now define the variational problem zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)) a00 = ufl.inner(ufl.grad(u), ufl.grad(du)) * ufl.dx a01 = ufl.inner(lmbda, du) * ufl.dx a10 = ufl.inner(u, dl) * ufl.dx L0 = ufl.inner(f, du) * ufl.dx + ufl.inner(g, du) * ufl.ds L1 = ufl.inner(zero, dl) * ufl.dx a = [[a00, a01], [a10, None]] L = [L0, L1] a_compiled = dolfinx.fem.form(a) L_compiled = dolfinx.fem.form(L) # Note that we have defined the variational form in a block form, and # that we have not included $h$ in the variational form. We will enforce this # once we have assembled the right hand side vector. # We can now assemble the matrix and vector A = dolfinx.fem.petsc.assemble_matrix(a_compiled, kind="nest", bcs=[]) A.assemble() # On the main branch of DOLFINx, the `assemble_vector` function for blocked spaces has been rewritten to reflect how # it works for standard assembly and `nest` assembly. This means that lifting is applied manually. # In this case, with no Dirichlet BC, we could skip those steps. # However, for clarity we include them here. bcs = [] b = dolfinx.fem.petsc.assemble_vector(L_compiled, kind="nest") apply_lifting_and_set_bc(b, a_compiled, bcs=bcs) # Next, we modify the second part of the block to contain `h` # We start by enforcing the multiplier constraint $h$ by modifying the right hand side vector. # On the main branch, this is greatly simplified uh = dolfinx.fem.Function(V, name="u") po = { "ksp_type":"preonly", "pc_type":"lu", "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True, } problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options=po, kind="nest" ) uh, rh = problem.solve() print(problem.solver.getConvergedReason()) diff = uh - u_exact(x) error = dolfinx.fem.form(ufl.inner(diff, diff) * ufl.dx) print(f"L2 error: {np.sqrt(assemble_scalar(error)):.2e}") ``` ### Suggested user interface ```python3 ```

On the main branch of DOLFINx_MPC the nest interface has gotten a serious face-lift:

github.com/jorgensd/dolfinx_mpc

python/demos/demo_stokes_nonlinear_nest.py

main

# Copyright (C) 2022-2025 Nathan Sime and Jørgen S. Dokken
#
# This file is part of DOLFINX_MPC
#
# SPDX-License-Identifier:    MIT
#
# This demo illustrates how to apply a slip condition on an
# interface not aligned with the coordiante axis.
# The demos solves the Stokes problem using the nest functionality to
# avoid using mixed function spaces and as a non-linear problem.

from __future__ import annotations

from pathlib import Path

from mpi4py import MPI
from petsc4py import PETSc

import basix
import dolfinx.io

This file has been truncated. show original

which makes it way easier to use than it was previously.

Gaurav_Verma · July 5, 2025, 1:42am

I tried implement your suggestions for this geometry but I am getting some errors. I think there is some syntax mistake. Your intervention will be helpful here. Here is the geometry and my code.

// Gmsh project created on Tue May  6 11:35:50 2025
//+
Point(1) = {0, 0.509795579104159, 0, 1.0};
//+
Point(2) = {0.195090322016128, 0.470989701299072, 0, 1.0};
//+
Point(3) = {0, 1.01959115820832, 0, 1.0};
//+
Point(4) = {0.390180644032256, 0.941979402598143, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 4};
//+
Line(3) = {4, 3};
//+
Line(4) = {3, 1};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Surface(1) = {1};
//+
Physical Curve("ps1", 5) = {2};
//+
Physical Curve("ps2", 6) = {4};
//+
Physical Curve("bottom", 7) = {1};
//+
Physical Curve("top", 8) = {3};
//+
Physical Surface("my_surf", 9) = {1};
//+
Show "*";

my code

import gmsh
import os
import math
import numpy as np
import matplotlib.pyplot as plt
import sys
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx import (mesh, fem, io, common, default_scalar_type, geometry, plot)
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx.plot import vtk_mesh
from dolfinx.io import XDMFFile, VTKFile
from ufl import (FacetNormal, as_matrix, as_tensor, extract_blocks, Identity, Measure, TestFunctions, tr, TrialFunctions, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedFunctionSpace)
from dolfinx_mpc.utils import (create_point_to_point_constraint, determine_closest_block, rigid_motions_nullspace, facet_normal_approximation, create_normal_approximation, gather_transformation_matrix, gather_PETScMatrix, gather_PETScVector)
from dolfinx_mpc import LinearProblem, MultiPointConstraint
from dolfinx.cpp.la.petsc import get_local_vectors
from packaging.version import Version
from dolfinx.cpp.la.petsc import scatter_local_vectors, get_local_vectors
from scifem import create_real_functionspace, assemble_scalar
from scifem.petsc import apply_lifting_and_set_bc, zero_petsc_vector

# Material properties
E = fem.Constant(domain, 200e9)       # Young's modulus
nu = fem.Constant(domain, 0.3)        # Poisson's ratio

lmbda = E * nu / ((1 + nu)*(1 - 2*nu))
mu = E / (2*(1 + nu))

# Strain with plane strain + extra ezz
def eps(u, ezz):
    return sym(as_tensor([
        [u[0].dx(0), u[0].dx(1), 0],
        [u[1].dx(0), u[1].dx(1), 0],
        [0, 0, ezz]
    ]))

def sigma(u, ezz):
    return lmbda * tr(eps(u, ezz)) * Identity(3) + 2 * mu * eps(u, ezz)


# Function spaces: V for displacement, R for scalar ezz
degree = 2
V = fem.functionspace(domain, ("P", degree, (2,)))
R = create_real_functionspace(domain)
W = MixedFunctionSpace(V, R)


# Trial and test functions
(du, dezz) = TrialFunctions(W)
(v, vezz) = TestFunctions(W)

zero_vec = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0)))

# Variational forms
a00 = inner(sigma(du, 0), eps(v, 0)) * dx
a01 = inner(sigma(zero_vec, dezz), eps(v, 0)) * dx
a10 = inner(sigma(du, 0), eps(zero_vec, vezz)) * dx
a11 = inner(sigma(zero_vec, dezz), eps(zero_vec, vezz)) * dx


T = fem.Constant(domain, PETSc.ScalarType((0.0, 1e6)))
n = FacetNormal(domain)
L0 = inner(T, v) * ds
L1 = fem.Constant(domain, PETSc.ScalarType(0.0)) * vezz * dx


# Assemble nested matrix and vector
a = [[a00, a01], [a10, a11]]
L = [L0, L1]
a_compiled = fem.form(a)
L_compiled = fem.form(L)


# Slip constraints on selected facets (example facets 5 and 6)

mpc = MultiPointConstraint(V)

for tag in [5, 6]:  # ps1 and ps2
    normal_vec = create_normal_approximation(V, facet_markers, tag)
    mpc.create_slip_constraint(V, (facet_markers, tag), normal_vec)

mpc.finalize()

# Assemble matrix WITHOUT mpc argument
A = assemble_matrix(a_compiled, bcs=[])
A.assemble()

# Assemble RHS vector
b = assemble_vector(L_compiled)
apply_lifting(b, [a_compiled], bcs=[])
set_bc(b, [])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)


# Create and setup solver
ksp = PETSc.KSP().create(A.getComm())
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")

# Solve system
X = A.createVecRight()
ksp.solve(b, X)

# Split solution vector
x_u, x_ezz = X.getNestSubVecs()

# Distribute solution to function space
u_sol = fem.Function(mpc.function_space)
mpc.distribute_solution(x_u, u_sol.vector)

ezz_sol = x_ezz.array

dokken · July 5, 2025, 8:15am

Depending on your DOLFINx version, you need to:

Stable/0.9.0

Load your geometry as usual (you have omitted this part of the code)
Create the two function spaces V, R
Create MPC on V
Use the assembly procedure from dolfinx_mpc/python/demos/demo_stokes_nest.py at v0.9.2 · jorgensd/dolfinx_mpc · GitHub
Create your ksp object as in: dolfinx_mpc/python/demos/demo_stokes_nest.py at v0.9.2 · jorgensd/dolfinx_mpc · GitHub

Main branch
Step 1-3 are the same.
4-5. Use dolfinx_mpc.LinearProblem or dolfinx_mpc.NonLinearProblem as described in: dolfinx_mpc/python/demos/demo_stokes_nonlinear_nest.py at main · jorgensd/dolfinx_mpc · GitHub

Ergonsa · July 11, 2025, 8:00pm

Hello @dokken, I have some doubts regarding this topic. I want to combine the “lagrange multipliers” from the real spaces (scifem package) with the periodic boundary conditions (dolfinx_mpc package). I have seen that in some recent posts (Non-linear blocked newton solver with real spaces #93) you managed to make it work for a non-linear problem but using block-assembly. Can this be done with nest assembly since

Block-assembly is not implemented in DOLFINx_MPC, only nest assembly

?
Thank you in advance.

dokken · July 12, 2025, 12:52pm

Yes, it can be done. One has to make a custom implementation of the assemble_residual required by the Newton-solver, but it is possible.

Ergonsa · July 12, 2025, 5:39pm

Ok, thank you very much for the information. So, to make it clear, I should combine the assemble_residual_mpc from dolfinx_mpc/problem.py with the _assemble_residual of scifem/solvers.py, right?

Topic		Replies	Views
Mixing dolfinx_mpc and scifem packages for a nonlinear problem dolfinx boundary-conditions , dolfinx	2	45	July 8, 2025
Explicitly assembling matrices with dolfinx_mpc General	31	1052	August 11, 2023
Assembly of matrices for use in external solver General boundary-conditions , assembly , dolfinx	4	101	February 18, 2025
Issues Solving Pure Neumann Problems in dolfinx dolfinx	9	1546	May 1, 2024
Nested Newton solver for mpc dolfinx	2	95	December 5, 2024

Using dolfinx_mpc in generalised plain strain formulation

Related topics