# Non periodic corner - eigenvalue problem with Slepc

**URL:** <https://fenicsproject.discourse.group/t/non-periodic-corner-eigenvalue-problem-with-slepc/19580>\
**Category:** dolfinx\
**Tags:** boundary-conditions, slepc\
**Created:** [March 6, 2026, 12:34pm UTC](https://fenicsproject.discourse.group/t/non-periodic-corner-eigenvalue-problem-with-slepc/19580 "2026-03-06T12:34:40Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Felipe\_Rocha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/felipe_rocha/32/7071_2.png) [@Felipe\_Rocha](https://fenicsproject.discourse.group/u/Felipe_Rocha)\
**Post date:** [March 6, 2026, 12:34pm UTC](https://fenicsproject.discourse.group/t/non-periodic-corner-eigenvalue-problem-with-slepc/19580/1 "2026-03-06T12:34:40Z")

</div>

```auto
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, plot
import ufl
import pyvista as pv
from petsc4py import PETSc
from slepc4py import SLEPc
from dolfinx_mpc import MultiPointConstraint, assemble_matrix

def slave_boundary_x(x):
    return np.isclose(x[0], 1.0)

def slave_boundary_y(x):
    return np.logical_and(np.isclose(x[1], 1.0), np.logical_not(np.isclose(x[0], 1.0)))
    # return np.isclose(x[1], 1.0) # first try

# Second try
# def slave_corner(x):
# return np.logical_and(np.isclose(x[0],1.0),
# np.isclose(x[1],1.0))

# mapping from slave -> master
def periodic_relation_x(x):
    out = np.zeros_like(x)
    out[0] = x[0] - 1.0
    out[1] = x[1]
    return out

# mapping from slave -> master
def periodic_relation_y(x):
    out = np.zeros_like(x)
    out[0] = x[0]
    out[1] = x[1] - 1.0
    return out

# Second try
# def periodic_relation_corner(x):
# out = np.zeros_like(x)
# out[0] = x[0] - 1.0
# out[1] = x[1] - 1.0
# return out

# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary_x, periodic_relation_x, [])
mpc.create_periodic_constraint_geometrical(V, slave_boundary_y, periodic_relation_y, [])
# mpc.create_periodic_constraint_geometrical(V, slave_corner, periodic_relation_corner, []) # second try
mpc.finalize()

W = mpc.function_space

# Variational forms
u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
m = ufl.inner(u, v) * ufl.dx

# Assemble matrices
A = assemble_matrix(fem.form(a), mpc)
A.assemble()
B = assemble_matrix(fem.form(m), mpc)
B.assemble()

# Slepc settings
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setDimensions(10) # number of eigenvalues
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)

# Solving and recovering eigenvectors
eps.solve()
nconv = eps.getConverged()
eigvec_list = []
for i in range(nconv):
    eigval = eps.getEigenvalue(i)
    xr = A.createVecRight()
    eps.getEigenvector(i, xr)
    print("Eigenvalue:", eigval)
    eigvh = fem.Function(W, name = "eigvec{0}".format(i))
    eigvh.x.array[:] = xr[:]
    mpc.backsubstitution(eigvh)
    eigvec_list.append(eigvh)

# Plot pyvista
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

plotter = pv.Plotter()

grid.point_data["eigv"] = eigvec_list[8].x.array # also in 12, ...  
plotter.add_mesh(grid, scalars="eigv")

plotter.view_xy()
plotter.enable_parallel_projection()
plotter.add_scalar_bar("eigv")
plotter.show()

```

Dear all,

Here you have MWE of an eigenvalue problem solved by Slepc. The periodic boundary conditions are imposed by dolfinx\_mpc. I’ve noticed some weird behaviour on the upper-right corner for some eigenmodes. I’m using dolfinx 0.10, serial mode. Do you any idea how to fix that? I have left in comments some wrong fixes that I tried.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 7, 2026, 1:50pm UTC](https://fenicsproject.discourse.group/t/non-periodic-corner-eigenvalue-problem-with-slepc/19580/2 "2026-03-07T13:50:00Z")

</div>

What about enforce the constraint in one go, i.e.

```python
def slave_boundary(x):
    return np.isclose(x[0], 1.0) | np.isclose(x[1], 1.0)

# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], 1.0)
    out[1] = x[1] - np.isclose(x[1], 1.0)
    return out

# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])
mpc.finalize()

```

resulting in the complete code:

```python
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, plot
import ufl
import pyvista as pv
from slepc4py import SLEPc
from dolfinx_mpc import MultiPointConstraint, assemble_matrix

def slave_boundary(x):
    return np.isclose(x[0], 1.0) | np.isclose(x[1], 1.0)

# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], 1.0)
    out[1] = x[1] - np.isclose(x[1], 1.0)
    return out

# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])
mpc.finalize()

W = mpc.function_space

# Variational forms
u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
m = ufl.inner(u, v) * ufl.dx

# Assemble matrices
A = assemble_matrix(fem.form(a), mpc)
A.assemble()
B = assemble_matrix(fem.form(m), mpc)
B.assemble()

# Slepc settings
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setDimensions(10) # number of eigenvalues
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)

# Solving and recovering eigenvectors
eps.solve()
nconv = eps.getConverged()
eigvec_list = []
for i in range(nconv):
    eigval = eps.getEigenvalue(i)
    eigvh = fem.Function(W, name = "eigvec{0}".format(i))
    eps.getEigenvector(i, eigvh.x.petsc_vec)
    eigvh.x.scatter_forward()
    print("Eigenvalue:", eigval)
    mpc.backsubstitution(eigvh)
    eigvec_list.append(eigvh)

# Plot pyvista
topology, cell_types, geometry = plot.vtk_mesh(W)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

plotter = pv.Plotter()

grid.point_data["eigv"] = eigvec_list[8].x.array # also in 12, ...  
plotter.add_mesh(grid, scalars="eigv")

plotter.view_xy()
plotter.enable_parallel_projection()
plotter.add_scalar_bar("eigv")
plotter.show()

```

which to me yields:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b4effc861ea74e36cd4e0cbbd5a684a2d5a78318.jpeg)

---

<div class="post-metadata">

**Author:** ![Felipe\_Rocha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/felipe_rocha/32/7071_2.png) [@Felipe\_Rocha](https://fenicsproject.discourse.group/u/Felipe_Rocha)\
**Post date:** [March 10, 2026, 3:06pm UTC](https://fenicsproject.discourse.group/t/non-periodic-corner-eigenvalue-problem-with-slepc/19580/3 "2026-03-10T15:06:43Z")

</div>

Thanks! It has solved the issue.