Contact for mixed Elements

[JohannesMenning]

Hello everyone,

I would like to solve a coupled problem, e.g. thermo-mechanical with the displacement u and the temperature as the field variables. In the simulation, two bodies are in contact with each other, and therefore the necessary contact conditions have to be defined.

I am wondering if it is possible to define a mechanical contact, as is done here, and the necessary contact conditions for the temperature field, e.g. $T_1 = T_2$. For this I would then like to use the function "create_contact_inelastic_condition()", using a monolithic scheme.

Here is a "MME" in which I have tried to implement the idea. I have adopted the Code from the example here and just added the mixed Elements, and deleted everything besides the definition of the MPC.

I truly appreciate any advice you can share and look forward to hearing your tips!

from __future__ import annotations
import warnings
from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
from pathlib import Path

from mpi4py import MPI
from petsc4py import PETSc

import dolfinx.fem as fem
import numpy as np
import scipy.sparse.linalg
from dolfinx import default_real_type, default_scalar_type
from dolfinx.common import Timer, list_timings
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType
from ufl import Identity, TestFunction, TrialFunction, dx, grad, inner, sym, tr
import ufl

from create_and_export_mesh import gmsh_3D_stacked, mesh_3D_dolfin
from dolfinx_mpc import MultiPointConstraint, apply_lifting, assemble_matrix, assemble_vector
from dolfinx_mpc.utils import (
    compare_mpc_lhs,
    compare_mpc_rhs,
    create_normal_approximation,
    gather_PETScMatrix,
    gather_PETScVector,
    gather_transformation_matrix,
    log_info,
    rigid_motions_nullspace,
    rotation_matrix,
)

import dolfinx.fem.petsc
import dolfinx.nls.petsc
import dolfinx_mpc
import dolfinx.la as _la
from dolfinx import log
import basix



def demo_stacked_cubes(
    outfile_u: XDMFFile,
    outfile_phi: XDMFFile,
    theta: float,
    gmsh: bool = False,
    ct: CellType = CellType.tetrahedron,
    compare: bool = True,
    res: float = 0.1,
    noslip: bool = False,
):
    celltype = "hexahedron" if ct == CellType.hexahedron else "tetrahedron"
    type_ext = "no_slip" if noslip else "slip"
    mesh_ext = "_gmsh_" if gmsh else "_"
    log_info(f"Run theta:{theta:.2f}, Cell: {celltype}, GMSH {gmsh}, Noslip: {noslip}")

    
    if gmsh:
        mesh, mt = gmsh_3D_stacked(celltype, theta, res)
        tdim = mesh.topology.dim
        fdim = tdim - 1
        mesh.topology.create_connectivity(tdim, tdim)
        mesh.topology.create_connectivity(fdim, tdim)
    else:
        if default_real_type == np.float32:
            warnings.warn("Demo does not run for single float precision due to limited xdmf support")
            exit(0)
        mesh_3D_dolfin(theta, ct, celltype, res)
        MPI.COMM_WORLD.barrier()
        with XDMFFile(MPI.COMM_WORLD, f"meshes/mesh_{celltype}_{theta:.2f}.xdmf", "r") as xdmf:
            mesh = xdmf.read_mesh(name="mesh")
            tdim = mesh.topology.dim
            fdim = tdim - 1
            mesh.topology.create_connectivity(tdim, tdim)
            mesh.topology.create_connectivity(fdim, tdim)
            mt = xdmf.read_meshtags(mesh, "facet_tags")
    mesh.name = f"mesh_{celltype}_{theta:.2f}{type_ext}{mesh_ext}"

    
   
    Ele_u = basix.ufl.element("CG", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, )) # fuer Verschiebung
    Ele_phi = basix.ufl.element("CG", mesh.topology.cell_name(), 1) # fuer elektrisches Potential
    mixed_ele = basix.ufl.mixed_element([Ele_u, Ele_phi])

    V_ges = fem.functionspace(mesh, mixed_ele)

    w = fem.Function(V_ges)

    
    mpc = MultiPointConstraint(V_ges)

    tol = float(5e2 * np.finfo(default_scalar_type).resolution)
    with Timer("~Contact: Create contact constraint"):
        mpc.create_contact_inelastic_condition(mt, 4, 9, eps2=tol) # for scalar field
        nh_u = create_normal_approximation(V_ges, mt, 4) # for displacement, i.e. mechanical field
        mpc.create_contact_slip_condition(mt, 4, 9, nh_u, eps2=tol) # for displacement, i.e. mechanical field


    with Timer("~~Contact: Add data and finialize MPC"):
        mpc.finalize()

    
if __name__ == "__main__":
    parser = ArgumentParser(formatter_class=ArgumentDefaultsHelpFormatter)
    parser.add_argument("--res", default=0.1, type=np.float64, dest="res", help="Resolution of Mesh")
    parser.add_argument(
        "--theta",
        default=np.pi / 3,
        type=np.float64,
        dest="theta",
        help="Rotation angle around axis [1, 1, 0]",
    )
    hex = parser.add_mutually_exclusive_group(required=False)
    hex.add_argument("--hex", dest="hex", action="store_true", help="Use hexahedron mesh", default=False)
    slip = parser.add_mutually_exclusive_group(required=False)
    slip.add_argument(
        "--no-slip",
        dest="noslip",
        action="store_true",
        help="Use no-slip constraint",
        default=False,
    )
    gmsh = parser.add_mutually_exclusive_group(required=False)
    gmsh.add_argument(
        "--gmsh",
        dest="gmsh",
        action="store_true",
        help="Gmsh mesh instead of built-in grid",
        default=False,
    )
    comp = parser.add_mutually_exclusive_group(required=False)
    comp.add_argument(
        "--compare",
        dest="compare",
        action="store_true",
        help="Compare with global solution",
        default=False,
    )
    time = parser.add_mutually_exclusive_group(required=False)
    time.add_argument("--timing", dest="timing", action="store_true", help="List timings", default=False)

    args = parser.parse_args()
    outdir = Path("results")
    outdir.mkdir(exist_ok=True, parents=True)
    outfile_u = XDMFFile(MPI.COMM_WORLD, outdir / "demo_contact_3D_u.xdmf", "w")
    outfile_phi = XDMFFile(MPI.COMM_WORLD, outdir / "demo_contact_3D_phi.xdmf", "w")

    ct = CellType.hexahedron if args.hex else CellType.tetrahedron

    demo_stacked_cubes(
        outfile_u,
        outfile_phi,
        theta=args.theta,
        gmsh=args.gmsh,
        ct=ct,
        compare=args.compare,
        res=args.res,
        noslip=args.noslip,
    )

    outfile_u.close()
    outfile_phi.close()

    log_info("Simulation finished")
    if args.timing:
        list_timings(MPI.COMM_WORLD)

``

[jorgensd]

Hi @JohannesMenning,

I think a better approach to work on this would be to use ideas of block/nest systems, as shown in:
https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_stokes_nest.py
as it would more easily work with DOLFINx_MPC in its current state.

In the future, there might be a smoother interface to this using ufl.MixedFunctionSpace and ufl.extract_blocks (it should work right now, but I don't have the bandwidth to make a demo). See for instance: https://jsdokken.com/FEniCS-workshop/src/multiphysics/submeshes.html for the general idea of blocked systems with a ufl.MixedFunctionSpace.

[JohannesMenning]

Hello @jorgensd

Thank you very much, for your fast reply, I will look into the links and tutorials which you have recommended!