@dokken , here the results of the test I made. I think in this way it is not easy to identify the slaves from the masters with a geometric criteria. You can find here the code I am playing with (for simplicity, I have formulated the problem as 2 squares that are sliding rather than 2 concentric circles that rotates).
This is the mesh generation through gmsh:
import gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("2squares")
r = 0.1
cx, cy = 0.5, 0.6
# Create base surfaces
rect_left = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
disk = gmsh.model.occ.addDisk(cx, cy, 0.0, r, r)
rect_right = gmsh.model.occ.addRectangle(1, 0, 0, 1, 1)
gmsh.model.occ.synchronize()
# Fragment: object = left rect, tools = disk + right rect
outDimTags, outMap = gmsh.model.occ.fragment([(2, rect_left)], [(2, disk)])
gmsh.model.occ.synchronize()
for i, (dim, ctag) in enumerate(gmsh.model.getEntities(2), start=1):
pg = gmsh.model.addPhysicalGroup(2, [ctag], tag=i)
gmsh.model.setPhysicalName(2, pg, f"surf_{ctag}")
for i, (dim, ctag) in enumerate(gmsh.model.getEntities(1), start=1):
pg = gmsh.model.addPhysicalGroup(1, [ctag], tag=100 + i)
gmsh.model.setPhysicalName(1, pg, f"line_{ctag}")
# Mesh and write
# Keep left-square sizes unchanged; refine only in the right square.
def right_square_refine(dim, tag, x, y, z, lc):
if x > 1.0 and 0.0 <= y <= 1.0:
return min(lc, 0.02)
return lc
#make right square mesh finer to see the non conformality
field = gmsh.model.mesh.field.add("Constant")
gmsh.model.mesh.field.setNumber(field, "VIn", 0.19) # size inside selected entities
gmsh.model.mesh.field.setNumbers(field, "SurfacesList", [rect_right])
gmsh.model.mesh.field.setAsBackgroundMesh(field)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.refine()
gmsh.write("2squares.msh")
gmsh.finalize()
While this is the implementation in dolfinx:
import dolfinx
import dolfinx_mpc
from dolfinx.io import gmsh
import numpy as np
import pyvista
import ufl
from mpi4py import MPI
#cell tags are -> 1: disk, 2 right square, 3 left square
mesh,ct,ft,*_ = gmsh.read_from_msh("2squares.msh",MPI.COMM_WORLD,0)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
disk = 1
right_square = 2
left_square = 3
left_square_contact = 104
right_square_contact = 107
tdim = mesh.topology.dim
shift = 0.15
tags_to_shift = right_square
#shift upward domain with tag 2
mesh.topology.create_connectivity(tdim, 0)
c2v = mesh.topology.connectivity(tdim, 0)
mask = np.isin(ct.values, tags_to_shift)
cells = ct.indices[mask]
verts = np.unique(np.concatenate([c2v.links(c) for c in cells]))
sp = dolfinx.fem.functionspace(mesh, ("CG", 2))
nu0 = 1/(4*3.14*1e-7)
def shift_upward(x):
x[1,:]+=1
return x
#dirichlet bc on left edge
bc = dolfinx.fem.dirichletbc(0.0, dolfinx.fem.locate_dofs_geometrical(sp, lambda x: np.isclose(x[0], 0.0)), sp)
#mpc: antiperiodic upper and lower facets, contact between left and right square, antiperiodic between right and left edge for the unmatched parts
mpc = dolfinx_mpc.MultiPointConstraint(sp)
#antiperiodic condition on left square lower-upper edge
mpc.create_periodic_constraint_geometrical(sp, lambda x: np.bitwise_and(np.isclose(x[1], 0.0), x[0]<=1),shift_upward, [bc],scale=-1.0)
#antiperiodic condition right square lower-upper edge
mpc.create_periodic_constraint_geometrical(sp, lambda x: np.bitwise_and(np.isclose(x[1], 0.),x[0]>=1),shift_upward, [bc],scale=-1.0)
#combination of contact and antiperiodicity on the sliding face
mpc.create_periodic_constraint_geometrical(sp, lambda x: np.bitwise_and(np.isclose(x[0], 1.), x[1]<shift),shift_upward, [bc],scale=-1.0)
###
#VERSION WITH CONTACT CONDITION...
#mesh.geometry.x[verts, 1] += shift
#mpc.create_contact_inelastic_condition(ft, left_square_contact, right_square_contact, eps2=1e-11,allow_missing_masters=True)
#mesh.geometry.x[verts, 1] -= shift
#VERSION WITH GEOMETRICAL CONSTRAINTS(NOT WORKING)
mpc.create_periodic_constraint_geometrical(sp, lambda x: np.bitwise_and(np.isclose(x[0], 1.), x[1]>=shift,x[1]<=1.0), lambda x: np.vstack((x[0], x[1] - shift)), [bc],scale=1.0)
#finalize mpc
mpc.finalize()
#linear form: source inside domain disk, all materials are linear with nu0
a = nu0*ufl.inner(ufl.grad(ufl.TrialFunction(sp)), ufl.grad(ufl.TestFunction(sp)))*ufl.dx
L= 10000.0*ufl.TestFunction(sp)*dx(disk)
#solve linear problem
problem = dolfinx_mpc.LinearProblem(a,L, mpc=mpc,bcs=[bc])
sol_lin= problem.solve()
from dolfinx import plot
owned = mesh.topology.index_map(mesh.topology.dim).size_local
ghosts = mesh.topology.index_map(mesh.topology.dim).num_ghosts
topology, cell_types, geometry = plot.vtk_mesh(mpc.function_space,np.arange(owned + ghosts))
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["ulin"] = sol_lin.x.array
grid.set_active_scalars("ulin")
vmin = min(sol_lin.x.array.min(), sol_lin.x.array.min())
vmax = max(sol_lin.x.array.max(), sol_lin.x.array.max())
levels = np.linspace(vmin, vmax, 25)
plotter = pyvista.Plotter()
iso_lin = grid.contour(isosurfaces=levels, scalars="ulin")
plotter.add_mesh(grid, scalars="ulin", cmap="viridis", clim=[vmin, vmax], show_edges=False)
plotter.add_mesh(grid, style="wireframe", color="k", opacity=0.2)
plotter.add_mesh(iso_lin, color="white", line_width=2)
plotter.add_text("ulin", font_size=12)
plotter.view_xy()
plotter.show()
Using the part with #VERSION WITH CONTACT CONDITION makes the results I expect (with the hack of manual shift of mesh coords upwards and downwards to make the contact_inelastic work as I want) even if there is there is still the issues at the corners due to the multiple periodicities overapping. If I use #VERSION WITH GEOMETRICAL CONSTRAINTS, where the slaves (on left square) are mapped to the master (on right square) with a downward shift, I get an unexpected result, with a zero potential on both sides of the contact edge. I guess this happen because all dofs on the contact, i.e. the ones on left and right square, are identified as slaves with no master, and the only thing that is making a nonzero solution on the right square is the antiperiodicity on the vertical noncontact edges.
