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)