Hi,
This is a follow up from the discussion here. It is a semi-simplification of a more complex 3D problem.
Consider the following Dirichlet BC problem with a global integral constraint.
- \partial _x^2u(x) = 0,~~x \in [0,1],~~u(0) = 0.5,~~u(1) = 1.5, ~~\int\limits_0^1 {dx} u(x)=2
This is the same as introducing a Lagrange multiplier \lambda (real constant since the constraint is global) and minimizing the “action”
J[u,\lambda ] = \int\limits_0^1 {dx} {\left( {{\partial _x}u(x)} \right)^2} - \lambda \left( {\int\limits_0^1 {dx} u(x) - 2} \right)
Minimization to \lambda gives the integral constraint and towards u gives - \partial _x^2u(x) = \lambda. Direct analytical solution gives u(x)=-6x^2+7x+0.5. Now since FEniCSx does not support real elements for \lambda in the aforementioned topic it was suggested by @conpierce8 to use @dokken’s dolfinx_mpc library to constrain all the dofs in that subspace to a single dof, and it works perfectly. I will post slightly modified minimal code since the problem is different from the original topic.
from mpi4py import MPI
import dolfinx
import dolfinx_mpc
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
comm = MPI.COMM_WORLD
### Create mesh and function spaces
domain = dolfinx.mesh.create_unit_interval(comm, 1000)
element1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, ufl.MixedElement([element1, element1]))
V0c, submap = V.sub(0).collapse()
V1 = V.sub(1)
### Create test and test functions
v, r = ufl.TestFunctions(V)
uh = dolfinx.fem.Function(V)
u, m = ufl.TrialFunctions(V)#
points = domain.geometry.x
# Export the mesh points and the boundary points
length=len(points)
### Mark boundaries
def right_boundary(x): return np.isclose(x[0], 1.0)
def left_boundary(x): return np.isclose(x[0], 0.0)
def domain_nodes(x): return ~np.logical_or(left_boundary(x), right_boundary(x))
node_left = dolfinx.mesh.locate_entities_boundary(domain, 0, left_boundary)
node_right = dolfinx.mesh.locate_entities_boundary(domain, 0, right_boundary)
node_middle = dolfinx.mesh.locate_entities_boundary(domain, 0, domain_nodes)
tags = {"left": 2, "right": 3, "middle": 1}
marker_left = np.full_like(node_left, tags["left"])
marker_right = np.full_like(node_right, tags["right"])
marker_middle = np.full_like(node_middle, tags["middle"])
nodes = np.concatenate((node_left, node_right, node_middle))
markers = np.concatenate((marker_left, marker_right, marker_middle))
idx = np.argsort(nodes)
ft = dolfinx.mesh.meshtags(domain, 0, nodes[idx], markers[idx])
### Create Dirichlet bcs
right_facets = dolfinx.fem.locate_dofs_geometrical((V.sub(0), V0c), right_boundary)
left_facets = dolfinx.fem.locate_dofs_geometrical((V.sub(0), V0c), left_boundary)
bc_r = dolfinx.fem.dirichletbc(ScalarType(1.5), right_facets[0], V.sub(0))
bc_l = dolfinx.fem.dirichletbc(ScalarType(0.5), left_facets[0], V.sub(0))
### MultiPointConstraint
mpc = dolfinx_mpc.MultiPointConstraint(V)
# Find the master dof (x = 0.0) and owning process
nodes = dolfinx.mesh.locate_entities_boundary(domain, 0, lambda x: np.isclose(x[0], 0.0))
master_dof = dolfinx.fem.locate_dofs_topological(V1, 0, nodes, remote=False)
master_owner = np.full_like(master_dof, comm.rank)
master_dof = np.concatenate(comm.allgather(master_dof), dtype=np.int64)
master_owner = np.concatenate(comm.allgather(master_owner), dtype=np.int32)
# Check for correct number of master dofs (should be 1)
if master_dof.size != 1:
raise RuntimeError(f"Expected to find 1 master dof; found {master_dof.size}")
else:
print(f"[{comm.rank}]: master_dof = {master_dof}")
print(f"[{comm.rank}]: master_owner = {master_owner}")
# Find the slave nodes
nodes = dolfinx.mesh.locate_entities(domain, 0, lambda x: ~np.isclose(x[0], 0.0))
slave_dof = dolfinx.fem.locate_dofs_topological(V1, 0, nodes, remote=False)
#print(f"[{comm.rank}]: slave_dof = {slave_dof}")
# Create arrays for use by mpc.add_constraint
masters = np.broadcast_to(master_dof, slave_dof.shape)
owners = np.broadcast_to(master_owner, slave_dof.shape)
coeffs = np.ones_like(masters, dtype=ScalarType)
offsets = np.arange(masters.size + 1, dtype=np.int32)
# Create constraint and finalize mpc
mpc.add_constraint(V, slave_dof, masters, coeffs, owners, offsets)
mpc.finalize()
### Create variational form and solve
C0 = dolfinx.fem.Constant(domain, ScalarType(2))
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
n = ufl.FacetNormal(domain)
ds = ufl.Measure("ds", domain, subdomain_data=ft)
a = (
ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
- ufl.inner(m, v)*ufl.dx
+ ufl.inner(u, r)*ufl.dx
#- ufl.inner(ufl.grad(u), n * v)*ds(tags["left"])
)
L = ufl.inner(f, v)*ufl.dx + ufl.inner(C0, r)*ufl.dx
problem = dolfinx_mpc.LinearProblem(
a, L, mpc, bcs=[bc_l,bc_r],
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
}
)
uh = problem.solve()
uh1=uh.x.array
new_array1 = [uh1[i] for i in range(len(uh1)) if i not in slave_dof]
new_array = [new_array1[i] for i in range(len(new_array1)) if i not in master_dof]
def u(x):
return -6 * (x**2) + 7 * x + 0.5
nn=len(points)
for i in range (nn):
print(points[i,0],new_array[i],u(points[i,0]))
I want to try to go this problem a step further by introducing a local constraint in a subdomain, \int\limits_0^{0.5} {dx} u(x)=2. If I am not terribly mistaken this means that \lambda takes a constant value in [0,0.5] and has zero contribution otherwise. The action now is J[u,\lambda ] = \int\limits_0^1 {dx} {\left( {{\partial _x}u(x)} \right)^2} - \lambda \left( {\int\limits_0^{0.5} {dx} u(x) - 2} \right) = \int\limits_0^1 {dx} \left[ {{{\left( {{\partial _x}u(x)} \right)}^2} - \lambda \left( {\theta (0.5 - x)u(x) - 2} \right)} \right] and the resulting PDE is - \partial _x^2u(x) = \lambda \theta (0.5 - x) where the step function \theta complicates the problem.
My questions and discussion topics are
- The above Lagrangian formalism does not seem elegant. Is there another way to formulate the problem?
- How could you suggest to tackle the - \partial _x^2u(x) = \lambda \theta (0.5 - x) and \int\limits_0^{0.5} {dx} u(x)=2 problem in FEniCSx?
- Could I maybe create a custom integral measure in FEniCSx and use code similar to the above? How can I do this while working on the dolfinx_mpc library? Or is there a way to split the calculation in two regions?
Any input is obviously highly appreciated.
PS. To solve the problem analytically you split the region in 0.5 and impose continuity in u and \partial _xu in this point and find a piecewise solution. Analytical solution: u(x) = \begin{cases} -\frac{312}{10}{x^2} + \frac{122}{5}x + {\frac{1}{2}}, ~~~ x \in [0,0.5] \\ -\frac{34}{5}x + {\frac{83}{10}},~~~~~~~~~~~~~~ x \in [0.5,1] \end{cases}