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(