Apply MPC to enforced “non-homogeneous” periodic BCs

[Tcadart]

Hello everyone,

I'm working on a model where I try to impose some specific boundary conditions. Those are close to periodic boundary conditions but instead of setting an equality between the solutions over the periodic interfaces, I am looking for applying a specific jump:

$$
\mathbf{u}{|\Gamma^+} - \mathbf{u}{\Gamma^-} = \mathbf{d}.
$$

Such BCs are of interest in computational mechanics for multi-scale homogenization, see for instance equation (1) in [P. Henyš et al., (2019), "Comparison of current methods for implementing periodic boundary conditions in multi-scale homogenisation".](https://doi.org/10.1016/j.euromechsol.2019.103825) This is actually what I am trying to do, on a unit cell modeled via timoshenko beams (see the piece of code below).

This seems to be trackable via dolfinx_mpc, but I didn’t succeed in doing so for the moment. In some sense, I want to enforce MPC of the following form (u) = [K] (û) + (d).
The function dolfinx_mpc.apply_lifting seems to do the trick. To my understanding, the imposed jump might be added via the input argument x0.

However, when attempting to use dolfinx_mpc.apply_lifting(b, [l_fem_form], [bcs], constraint=mpc), I encounter the following error: AttributeError: 'MultiPointConstraint' object has no attribute '_cpp_object'.

Here's a snippet of my setup:

# Setup mpc and periodicity
mpc = MultiPointConstraint(V)
mpc.add_constraint(V, dof_master, dof_slave, coeff, owner, offset)
# Problematic call
apply_lifting(b, [l_fem_form], [bcs], constraint=mpc)

Does anyone have an idea of what might be going wrong, or can suggest an alternative approach to resolve this issue?

Here is a functioning code without the “non-zero” jump (i.e. standard periodic boundary condition):

from mpi4py import MPI
from petsc4py import PETSc
import math
import numpy as np
import ufl
from ufl import sqrt, cross, dot, diag, as_vector, split, grad, dx
from dolfinx import fem, io, mesh, default_scalar_type
from dolfinx.fem import Constant, Expression
import dolfinx.fem.petsc
from dolfinx.common import Timer
from dolfinx_mpc import MultiPointConstraint, LinearProblem, apply_lifting
import gmsh
from Beam_parameters import *

#Geometry definition
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 0.1
gdim = 1
gmsh.initialize()
if mesh_comm.rank == model_rank:
    geom = gmsh.model.geo
    P0 = geom.addPoint(0, 0, 0, meshSize=h)
    P1 = geom.addPoint(1, 0, 0, meshSize=h)
    P2 = geom.addPoint(0, 1, 0, meshSize=h)
    P3 = geom.addPoint(1, 1, 0, meshSize=h)
    P4 = geom.addPoint(0.5, 0.5, 0, meshSize=h)
    L0 = geom.addLine(P0, P4)
    L1 = geom.addLine(P1, P4)
    L2 = geom.addLine(P2, P4)
    L3 = geom.addLine(P3, P4)
    geom.synchronize()
    gmsh.model.addPhysicalGroup(gdim, [L0,L1,L2,L3], 1)
    gmsh.model.setPhysicalName(gdim, 1, "Beam")
    gmsh.model.addPhysicalGroup(gdim - 1, [P0], 0)
    gmsh.model.addPhysicalGroup(gdim - 1, [P1], 1)
    gmsh.model.addPhysicalGroup(gdim - 1, [P2], 2)
    gmsh.model.addPhysicalGroup(gdim - 1, [P3], 3)
    gmsh.model.mesh.generate(gdim)
    gmsh.write("test.msh")
    gmsh.finalize()

domain, markers, facets = io.gmshio.read_from_msh("test.msh", MPI.COMM_WORLD)

# Material data
E = Constant(domain, 70e3)
nu = Constant(domain, 0.3)
G = E / 2 / (1 + nu)
rho = Constant(domain, 2.7e-3)

#Geometry data
radius = Constant(domain, 0.075)  # Radius
S = math.pi * radius ** 2  # Surface
I = (math.pi * (radius * 2) ** 4) / 64  # Bending inertia
ES = E * S
EI = E * I
J = (math.pi * (radius * 2) ** 4) / 32  # Torsional inertia
GJ = G * J
kappa = Constant(domain, 0.9)  # Shear area
GS = G * kappa * math.pi * radius ** 2

# Define Function Space on mesh P1/P1-Lagrange interpolation
element = ufl.VectorElement("CG", domain.ufl_cell(), 1, dim=domain.geometry.dim)
element_mixed = ufl.MixedElement([element, element])
V = fem.functionspace(domain, element_mixed)

# Define Test Function
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

def tangent(mesh): # Compute tangent for all beams in the mesh
    t = ufl.Jacobian(mesh)
    return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(ufl.inner(t, t))

t = tangent(domain)
ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))

def tgrad(u,t):
    return dot(grad(u), t)

def generalized_strains(u,t,a1,a2):
    (w, theta) = split(u)
    return as_vector([dot(tgrad(w,t), t),
                      dot(tgrad(w,t), a1)-dot(theta, a2),
                      dot(tgrad(w,t), a2)+dot(theta, a1),
                      dot(tgrad(theta,t), t),
                      dot(tgrad(theta,t), a1),
                      dot(tgrad(theta,t), a2)])
def generalized_stress(u,ES, GS, GJ, EI,t,a1,a2):
    return dot(diag(as_vector([ES, GS, GS, GJ, EI, EI])), generalized_strains(u,t,a1,a2))

Sig = generalized_stress(du,ES, GS, GJ, EI,t,a1,a2)
Eps =  generalized_strains(u_,t,a1,a2)

#Define Dirichlet Boundary Condition
pointList = np.array([facets.find(0),facets.find(1)], dtype=np.int32)
fixed_dofs = fem.locate_dofs_topological(V, domain.topology.dim-1, pointList)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, fixed_dofs)]

# Define z-axis force
force_magnitude = 0.0001
force_direction = np.array([0, -1, 0])
force_vector = force_magnitude * force_direction
force = fem.Constant(domain, force_vector)

# Define application domain
ds = ufl.Measure("ds", domain=domain, subdomain_data=facets)
dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 1})

# Define Variational form
k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
l_form = dot(force, w_)*ds(2)
l_fem_form = dolfinx.fem.form(l_form)

# Define periodic nodes
dof_master = fem.locate_dofs_topological(V, domain.topology.dim-1, facets.find(2))
dof_slave = fem.locate_dofs_topological(V, domain.topology.dim-1, facets.find(3))
coeff = np.ones(dof_slave.size)
owner = np.zeros(dof_slave.size,dtype=np.int32)
offset = np.arange(0,dof_slave.size+1)

b = PETSc.Vec().createSeq(dof_master.size)
b.set(0)
b.setValue(dof_slave[0], 1.0, addv=PETSc.InsertMode.INSERT_VALUES)
b.assemblyBegin()
b.assemblyEnd()
l_fem_form = dolfinx.fem.form(l_form)

with dolfinx.common.Timer("Multi constraint"):
    mpc = MultiPointConstraint(V)
    mpc.add_constraint(V,dof_master,dof_slave,coeff,owner,offset)
    # apply_lifting(b, [l_fem_form], [bcs], constraint=mpc)
mpc.finalize()

problem = LinearProblem(k_form, l_form, mpc, bcs=bcs, petsc_options={"ksp_type": "preonly","pc_type": "lu"})
u = problem.solve()

saveDataLinearBeam(domain, "reformated_file.pvd", u,ES, GS, GS, GJ, EI, EI,t,a1,a2)