Applying a constraint using a Lagrange multiplier dolfinx

Hi,

I’m trying to impose a global constraint using a Lagrange multiplier on dolfinx. I have created a minimal example to test the constraint, however, I’m receiving a different result compared with dolfin.

Here is a minimal working example (dolfin):

from dolfin import *

# Create mesh
mesh = UnitSquareMesh.create(50, 50, CellType.Type.quadrilateral)

# Build function space with Lagrange multiplier
V = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([V, R]))

bc1 = DirichletBC(W.sub(0), 0.0, "on_boundary")

bcs = [bc1]

# Define variational problem
(u, lm) = TrialFunctions(W)
(v, c) = TestFunctions(W)
g = Constant(0.0)
F = inner(grad(u), grad(v))*dx - g*v*dx
F += lm*v*dx + c*(u - Constant(0.5))*dx
a = lhs(F)
L = rhs(F)

# Compute solution
w = Function(W)
solve(a == L, w, bcs)
(u, lm) = w.split()

# Save solution in VTK format
file1 = File("lagrange_multipliers/u.pvd")
file1 << u
file2 = File("lagrange_multipliers/lm1.pvd")
file2 << lm

Here is a minimal working example (dolfinx):

import os
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile
import dolfinx
from dolfinx import (DirichletBC, Function, FunctionSpace, Constant)
from dolfinx.fem import locate_dofs_topological
from dolfinx.mesh import locate_entities_boundary
from ufl import (lhs, rhs, FiniteElement, MixedElement, TestFunctions, TrialFunctions, dx, 
grad, inner)
from dolfinx.la import solve


# Create mesh
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 50, 50, 
       dolfinx.cpp.mesh.CellType.quadrilateral)

def boundary(x):
    return np.logical_or(np.isclose(x[1], 1.0), np.logical_or(np.isclose(x[1], 0.0), 
    np.logical_or(np.isclose(x[0], 1.0), np.isclose(x[0], 0.0))))

# Build function space with Lagrange multiplier
V = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
R = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([V, R]))

facets = locate_entities_boundary(mesh, 1, boundary)

u_bc = Function(W.sub(0).collapse())
u_bc.vector.set(0.0)
u_bc.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
dofs = locate_dofs_topological((W.sub(0), FunctionSpace(mesh, V)), 1, facets)
bc1 = DirichletBC(u_bc, dofs, W.sub(0))
bcs = [bc1]

# Define variational problem
(u, lm) = TrialFunctions(W)
(v, c) = TestFunctions(W)
g = Constant(mesh, 0.0)
F = inner(grad(u), grad(v))*dx - g*v*dx
F += lm*v*dx + c*(u - Constant(mesh, 0.5))*dx
a = lhs(F)
L = rhs(F)

# Compute solution
w = Function(W)

A = dolfinx.fem.assemble_matrix(a, bcs)
A.assemble()
b = dolfinx.fem.assemble_vector(L)
dolfinx.fem.apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(b, bcs)
solve(A, w.vector, b)

with XDMFFile(MPI.COMM_WORLD, "u.xdmf", "w") as u_file:
    u_file.write_mesh(mesh)
    u_file.write_function(w.sub(0))

#with XDMFFile(MPI.COMM_WORLD, "lm.xdmf", "w") as lm_file:
#    lm_file.write_mesh(mesh)
#    lm_file.write_function(w.sub(1))

Another issue is that I’m receiving the following error when I try to save the Lagrange multiplier using XDMFFile

Loguru caught a signal: SIGSEGV
Segmentation fault (core dumped)

This is related to this issue around Real spaces not doing what they are meant to.

2 Likes

Hello,

as far as I’ve tracked the related issues within the git repository this issue hasn’t be solve yet or did I missed something?

Kind regards,
Maximilian

Good news: @jackhale has been working on Real spaces and there’s a now a PR open that adds support for them to DOLFINx (Add Real Space by jhale · Pull Request #2266 · FEniCS/dolfinx · GitHub)

4 Likes