I am trying to rework the disconnected elasticiy demo of dolfinx_mpc (0.7) on a curved boundary. The mesh across the boundary is non-conformal on purpose. Is the mpc expected to work for this geometry? The result I am getting with my attempt is clearly not working. But then I might have done it wrong. Below is the code
import dolfinx, mpi4py, ufl, gmsh, petsc4py, pyvista, dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt
def create_geom(ro=50, fov=None, disp=False):
if fov is None:
fov = (400, 300)
fov_x, fov_y = fov
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gdim = 2
model_rank = 0
t_gap = 5 # thickness of airgap
occ = gmsh.model.occ
# first we create the outer domain
outer_box = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
occ.synchronize()
air_shell = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
dom1 = occ.cut([(gdim, outer_box)], [(gdim, air_shell)]) # one domain
# now the inner domain which is the wheel
dom2 = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
occ.synchronize()
# number all domains
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1) # create the main group/node
# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1) # create the main group/node
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()
gmsh.model.mesh.refine()
# gmsh.write('geom/mwe_mixed_domains.msh')
model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
if disp:
gmsh.fltk.run() # visible inspection, for debugging
gmsh.finalize()
return mesh, ct, ft
mesh, ct, ft = create_geom(disp=False)
dom_idx_rot = 2 # rotating/circular domain
dom_idx_stat = 1 # stationary domain
bnd_idx_cont = (1, 6) # boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_l = (3, ) # left boundary
bnd_idx_r = (4, ) # right boundary
fem_degree = 2
Omega = dolfinx.fem.functionspace(mesh, ("Lagrange", fem_degree))
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)
F = ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx - ufl.dot(dolfinx.fem.Constant(mesh, 0.0), v)*ufl.dx
a = ufl.lhs(F)
L = ufl.rhs(F)
# dirichlet boundary conditions at the left and right edges
dofs_l = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_l[0]))
dbc_l = dolfinx.fem.dirichletbc(1.0, dofs_l, Omega)
dofs_r = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(0., dofs_r, Omega)
def gather_dof_coordinates(V, dofs: np.ndarray):
"""
Distributes the dof coordinates of this subset of dofs to all processors
"""
x = V.tabulate_dof_coordinates()
local_dofs = dofs[dofs < V.dofmap.index_map.size_local * V.dofmap.index_map_bs]
coords = x[local_dofs]
num_nodes = len(coords)
glob_num_nodes = mpi4py.MPI.COMM_WORLD.allreduce(num_nodes, op=mpi4py.MPI.SUM)
recvbuf = None
if mpi4py.MPI.COMM_WORLD.rank == 0:
recvbuf = np.zeros(3 * glob_num_nodes, dtype=V.mesh.geometry.x.dtype)
sendbuf = coords.reshape(-1)
sendcounts = np.array(mpi4py.MPI.COMM_WORLD.gather(len(sendbuf), 0))
mpi4py.MPI.COMM_WORLD.Gatherv(sendbuf, (recvbuf, sendcounts), root=0)
glob_coords = mpi4py.MPI.COMM_WORLD.bcast(recvbuf, root=0).reshape((-1, 3))
return glob_coords
dofs_stat = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_cont[0]))
dofs_rot = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_cont[1]))
# Given the local coordinates of the dofs, distribute them on all processors
nodes = [gather_dof_coordinates(Omega, dofs_stat), gather_dof_coordinates(Omega, dofs_rot)]
pairs = []
# points are on the circumference across the interface
for x in nodes[0]:
theta1 = np.arctan2(x[1], x[0])
for y in nodes[1]:
theta2 = np.arctan2(y[1], y[0])
if np.isclose(theta1, theta2):
pairs.append([x, y])
break
from dolfinx_mpc.utils import create_point_to_point_constraint
mpc = dolfinx_mpc.MultiPointConstraint(Omega)
for i, pair in enumerate(pairs):
sl, ms, co, ow, off = create_point_to_point_constraint(
Omega, pair[0], pair[1])
mpc.add_constraint(Omega, sl, ms, co, ow, off)
mpc.finalize()
bcs = [dbc_l, dbc_r]
# petsc_options = {"ksp_type": "preonly", "pc_type": "lu"} # does not work without dolfinx_mpc
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
from dolfinx_mpc import LinearProblem
problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()
import matplotlib.tri as tri
import matplotlib.animation as animation
fig1, ax1 = plt.subplots(constrained_layout=True, figsize=(4.5, 3.5))
ax1.set_aspect('equal')
triang = tri.Triangulation(uh.function_space.tabulate_dof_coordinates()[:,0],
uh.function_space.tabulate_dof_coordinates()[:,1])
tcf = ax1.tricontourf(triang, uh.x.array, 100, cmap='hot')
plt.colorbar(tcf)