Hi,

Am using dolfinx_mpc for the first time and am stuck with the implementation part. I wanted to couple the nodes along facet 3( all the dofs on this facet should have the same temperature). Can anyone guide me through how it should be implemented in a 2D case? I have attached my code where I haven’t implemented the multi-point constraint correctly:

$PhysicalNames

3

1 2 “convection”

1 3 “coupling”

2 1 “domain”

$EndPhysicalNames

import dolfinx

from dolfinx.fem import (Constant,dirichletbc, Function, FunctionSpace, LinearProblem,assemble_scalar, locate_dofs_geometrical, locate_dofs_topological)

from dolfinx.fem import (assemble_matrix, assemble_vector, form, apply_lifting, locate_dofs_geometrical, set_bc)

import ufl

from ufl import (SpatialCoordinate, TestFunction, TrialFunction,dx, ds, grad, inner)

from mpi4py import MPI

import dolfinx.io

import numpy as np

from petsc4py import PETSc

from petsc4py.PETSc import ScalarType

import dolfinx_mpc

import dolfinx_mpc.utils

from dolfinx.io import XDMFFile

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, “Mesh.xdmf”, “r”) as xdmf:

mesh = xdmf.read_mesh(name=“Grid”)

mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)

m_subdomains = xdmf.read_meshtags(mesh, name=“Grid”)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, “mt.xdmf”, “r”) as xdmf:

m_boundaries = xdmf.read_meshtags(mesh, name=“Grid”)

# function space

V = FunctionSpace(mesh, (“CG”, 1))

num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs

print(f"Number of dofs global: {num_dofs_global}")

dx = ufl.Measure(“dx”, domain=mesh, subdomain_data=m_subdomains)

ds = ufl.Measure(“ds”, domain=mesh, subdomain_data=m_boundaries)

# Solve Problem without MPC for reference

# trial and test functions

T = TrialFunction(V)

v = TestFunction(V)

# parameters

T_amb = Constant(mesh, 50.0) # ambient temperature

kappa = Constant(mesh, 52.0) # thermal conductivity [W/(m*K)]

h = Constant(mesh, 5.0) # heat Transfer Co-efficient [W/(m^2 deg C)]

# Robin boundary conditions

r1 = h

Integral_Ra = []

Integral_RL = []

Integral_Ra.append(r1 * T * v * ds(2))

Integral_RL.append(r1 * T_amb * v * ds(2))

# JZ: reformulate for steady state

a = kappa * inner(grad(T),grad(v)) * dx + sum(Integral_Ra)

L = sum(Integral_RL)

problem = LinearProblem(a, L, petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})

uh = problem.solve()

# Create multipoint constraint

def l2b(li):

return np.array(li, dtype=np.float64).tobytes()

n = dolfinx_mpc.utils.create_normal_approximation(V, m_boundaries, 3)

mpc = dolfinx_mpc.MultiPointConstraint(V)

mpc.finalize()

problem = dolfinx_mpc.LinearProblem(a, L, mpc ,petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})

u_h = problem.solve(