Resolve #98.
MWE:
```python
import basix.ufl
from dolfinx.io import XDMFFile…
from dolfinx.io import VTXWriter
from dolfinx_mpc import LinearProblem
import dolfinx
import mpi4py
import ufl
import gmsh
import dolfinx_mpc
from dolfinx.io import gmshio
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):
# create the main group/node
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)
# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
# create the main group/node
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1)
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)
with XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(ft, mesh.geometry)
xdmf.write_meshtags(ct, mesh.geometry)
dom_idx_rot = 2 # rotating/circular domain
dom_idx_stat = 1 # stationary domain
# boundary where continuity is imposed, stationary and rotating sides respectively
bnd_idx_cont = (5, 6)
bnd_idx_l = (3, ) # left boundary
bnd_idx_r = (4, ) # right boundary
fem_degree = 2
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
el = basix.ufl.element(
"Lagrange", mesh.topology.cell_name(), fem_degree)
Omega = dolfinx.fem.functionspace(mesh, el)
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)
F = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx - \
ufl.inner(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]))
cl = dolfinx.fem.Constant(mesh, (1.0))
dbc_l = dolfinx.fem.dirichletbc(cl, dofs_l, Omega)
cr = dolfinx.fem.Constant(mesh, (0.0))
dofs_r = dolfinx.fem.locate_dofs_topological(
Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(cr, dofs_r, Omega)
mpc = dolfinx_mpc.MultiPointConstraint(Omega)
tol = float(1e-11)
bcs = [dbc_l, dbc_r]
mpc.create_contact_inelastic_condition(
ft, bnd_idx_cont[0], bnd_idx_cont[1], eps2=tol)
mpc.finalize()
# 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"}
problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
uh = problem.solve()
with VTXWriter(mesh.comm, "uh.bp", [uh], engine="BP4") as vtx:
vtx.write(0.0)
```