Hi all,
I wanted to try imposing an integral constrain in a 3D Poisson problem. I tried to start from a simple case as demonstrated here for the legacy FEniCS.
Basically consider a 1D line in the interval [0,1] and the Laplace equation in this domain -\nabla^2u(x)=0 where we have the Dirichet Boundary Condition u(1)=1.5 and the integral constrain \int_{0}^{1}dxu(x)=C_0=1. The analytical answer is u(x)=x+0.5.
I defined the Lagrange multiplier \lambda and the Lagrangian L=\int_{0}^{1}dx\left(\frac{1}{2}\left|\nabla u(x)\right|^2-\lambda\left( \int_{0}^{1}dxu(x)-C_0\right)\right) and functional minimization with respect to u(x) and \lambda gives the following set of equations \nabla^2 u(x)=\lambda and the constrain \int_{0}^{1}dxu(x)=C_0.
So, seems that I need a Lagrange finite element for u(x) and a real one for \lambda. So, my questions for discussion are
- Does the above formulation seem correct?
- It seems, from what I read in the forum, that real finite elements are not incorporated in FEniCSx since the code I give below gives the error
Exception has occurred: ValueError Unknown element family: Real with cell type interval - Any ideas on if I can solve the problem by introducing Lagrange finite elements for the Lagrange multiplier? I tried it but the code outputs a wrong result.
Minimal Code for Real Element
from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import petsc4py.PETSc
domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
element1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
element2 = ufl.FiniteElement("Real", domain.ufl_cell(), 0)
mel = ufl.MixedElement([element1, element2])
V = dolfinx.fem.FunctionSpace(domain, mel)
V0, submap = V.sub(0).collapse()
# Test and Trial functions
v, r = ufl.TestFunctions(V)
uh = dolfinx.fem.Function(V)
u, m = ufl.TrialFunctions(V)
def right_boundary(x):
return np.isclose(x[0], 1)
right_facets = dolfinx.fem.locate_dofs_geometrical((V.sub(0),V0), right_boundary)
bc_r = dolfinx.fem.dirichletbc(ScalarType(1.5), right_facets[0], V.sub(0))
C0 = dolfinx.fem.Constant(domain, ScalarType(1))
f = dolfinx.fem.Constant(domain, ScalarType(0))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + m*v*ufl.dx + u * r * ufl.dx
L = f * v * ufl.dx + C0 * r * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc_r],u=uh)
uh = problem.solve()
u1, u2 = uh.split()
print(u1.x.array)
print(u2.x.array)
Minimal Code for Lagrange Element
from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import petsc4py.PETSc
domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
element1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
mel = ufl.MixedElement([element1, element1])
V = dolfinx.fem.FunctionSpace(domain, mel)
V0, submap = V.sub(0).collapse()
# Test and test functions
v, r = ufl.TestFunctions(V)
uh = dolfinx.fem.Function(V)
u, m = ufl.TrialFunctions(V)#
def right_boundary(x):
return np.isclose(x[0], 1)
right_facets = dolfinx.fem.locate_dofs_geometrical((V.sub(0),V0), right_boundary)
bc_r = dolfinx.fem.dirichletbc(ScalarType(1.5), right_facets[0], V.sub(0))
C0 = dolfinx.fem.Constant(domain, ScalarType(1))
f = dolfinx.fem.Constant(domain, ScalarType(0))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + m*v*ufl.dx + u * r * ufl.dx
L = f * v * ufl.dx + C0 * r * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc_r],u=uh)
uh = problem.solve()
u1, u2 = uh.split()
print(u1.x.array)
print(u2.x.array)