Add scalar constraint to inelastic condition

Resolve #98. MWE:
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)
jorgensd / dolfinx_mpc

Add scalar constraint to inelastic condition #99

Quality Gate passed
