How to apply periodic BC with a constant shift in dolfinx-mpc for a linear elasticity problem?

General
1

dolfinx-mpc: v0.7.2
I found tutorial that explains how we apply u(x=1) = u(x=0) and u(y=1) = u(y=0) constraints here. How can we achieve a constant shift with periodic BC, i.e., u_x(x=1) = u_x(x=0) + C, with u(y=1) = u(y=0) and C being a constant scalar?
Here is an illustration to express what I want to achieve: solve displacement field (blue) for this problem.
image

Here is what I was trying, but this is not correct.

from __future__ import annotations

from pathlib import Path
from typing import List, Tuple
import numpy as np

from mpi4py import MPI
import ufl

from petsc4py import PETSc

import dolfinx
import dolfinx.fem as fem
#from dolfinx.mesh import create_unit_square


import dolfinx_mpc.utils
from dolfinx_mpc import MultiPointConstraint, assemble_matrix, LinearProblem

import pyvista


#def assemble_and_solve(boundary_condition: List[str] = ["dirichlet", "periodic"], Dapp: float 1.0):
def assemble_and_solve(boundary_condition, Dapp):
    # Create mesh and finite element
    comm = MPI.COMM_WORLD
    N = 20
    Lx = 1
    Ly = 1
    mesh = dolfinx.mesh.create_unit_square(comm, 10, 10),
    mesh = mesh[0]
    V = fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
    fdim = mesh.topology.dim - 1


    locate_bc_fixed = lambda x: np.isclose(x[0], 0.0)
    fixed_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, locate_bc_fixed)
    u_fixed = np.array((0,) * mesh.geometry.dim, dtype=dolfinx.default_scalar_type)
    fixed_dofs_xy = fem.locate_dofs_topological(V, fdim, fixed_facets)
    bc_fixed = fem.dirichletbc(u_fixed, fixed_dofs_xy, V)

    '''u_bc0 = fem.Constant(mesh, PETSc.ScalarType(Dapp))
    facet_R = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0],Lx))
    dofsX_R = fem.locate_dofs_topological(V.sub(1), fdim, facet_R)
    facet_B = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[1],0.0))
    dofsY_B = fem.locate_dofs_topological(V.sub(1), fdim, facet_B)
    dof_RB = np.intersect1d(dofsX_R, dofsY_B)
    bc1 = fem.dirichletbc(u_bc0, dof_RB, V.sub(1))'''
    

    bcs = [bc_fixed]  
    pbc_directions = []
    pbc_slave_tags = []
    pbc_is_slave = []
    pbc_is_master = []
    pbc_meshtags = []
    pbc_slave_to_master_maps = []

    def generate_pbc_slave_to_master_map(i):
        def pbc_slave_to_master_map(x):
            out_x = x.copy()
            out_x[i] = x[i] - 1
            if i == 0:
                out_x[i] = out_x[i] + Dapp
            return out_x

        return pbc_slave_to_master_map

    def generate_pbc_is_slave(i):
        return lambda x: np.isclose(x[i], 1)

    def generate_pbc_is_master(i):
        return lambda x: np.isclose(x[i], 0)

    # Parse boundary conditions
    for i, bc_type in enumerate(boundary_condition):
        if bc_type == "dirichlet":
            u_bc = fem.Function(V)
            u_bc.x.array[:] = 0

            def dirichletboundary(x):
                return np.logical_or(np.isclose(x[i], 0), np.isclose(x[i], 1))

            facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, dirichletboundary)
            topological_dofs = fem.locate_dofs_topological(V, fdim, facets)
            bcs.append(fem.dirichletbc(u_bc, topological_dofs))

        elif bc_type == "periodic":
            pbc_directions.append(i)
            pbc_slave_tags.append(i + 2)
            pbc_is_slave.append(generate_pbc_is_slave(i))
            pbc_is_master.append(generate_pbc_is_master(i))
            pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

            facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, pbc_is_slave[-1])
            arg_sort = np.argsort(facets)
            pbc_meshtags.append(
                dolfinx.mesh.meshtags(
                    mesh,
                    fdim,
                    facets[arg_sort],
                    np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32),
                )
            )

    # Create MultiPointConstraint object
    mpc = MultiPointConstraint(V)
    N_pbc = len(pbc_directions)
    for i in range(N_pbc):
        if N_pbc > 1:

            def pbc_slave_to_master_map(x):
                out_x = pbc_slave_to_master_maps[i](x)
                idx = pbc_is_slave[(i + 1) % N_pbc](x)
                out_x[pbc_directions[i]][idx] = np.nan
                return out_x
        else:
            pbc_slave_to_master_map = pbc_slave_to_master_maps[i]

        mpc.create_periodic_constraint_topological(V, pbc_meshtags[i], pbc_slave_tags[i], pbc_slave_to_master_map, bcs)
    if len(pbc_directions) > 1:
        # Map intersection(slaves_x, slaves_y) to intersection(masters_x, masters_y),
        # i.e. map the slave dof at (1, 1) to the master dof at (0, 0)
        def pbc_slave_to_master_map(x):
            out_x = x.copy()
            out_x[0] = x[0] - 1 + Dapp
            out_x[1] = x[1] - 1
            idx = np.logical_and(pbc_is_slave[0](x), pbc_is_slave[1](x))
            out_x[0][~idx] = np.nan
            out_x[1][~idx] = np.nan
            return out_x

        mpc.create_periodic_constraint_topological(V, pbc_meshtags[1], pbc_slave_tags[1], pbc_slave_to_master_map, bcs)
    mpc.finalize()
    # Define variational problem
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    # Elasticity parameters
    E = 1.0e4
    nu = 0.3
    mu = fem.Constant(mesh, dolfinx.default_scalar_type(E / (2.0 * (1.0 + nu))))
    lmbda = fem.Constant(mesh, dolfinx.default_scalar_type(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))))
    T = fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
    f = fem.Constant(mesh, dolfinx.default_scalar_type((0, 0)))
    def sigma(v):
        return 2.0 * mu * ufl.sym(ufl.grad(v)) + lmbda * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(len(v))

    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ufl.ds

    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
    problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options)
    uh = problem.solve()

assemble_and_solve(["periodic", "periodic"], 0.5)
2

This is currently not supported in DOLFINx_MPC. It should be possible to implement within the MPC framework with some minor modifications in the source code (The key is to associate a single constant per constraint, and apply this to the RHS of the equation through lifting).

I’ve unfortunately not had time to implement this in the last year, even if I did make a prototype of it a few years ago.

1 Like
3

I see, will try out your suggestions about the required modifications. Thank you dokken!

4

What is your targeted application ? From the sketch I am assuming elasticity with periodic boundary conditions ?
If this is the case you can simply replace the total displacement \boldsymbol{u} with the periodic fluctuation \boldsymbol{v} as the main unknown such that:
\boldsymbol{u} = \mathbf{E}\mathbf{x} + \boldsymbol{v}
where \mathbf{E} is the imposed total strain which would be \mathbf{E} = \begin{bmatrix} c & 0 \\ 0 & 0 \end{bmatrix} in your example.
In this setting, \boldsymbol{v} is subject to periodic BC without any constant shift. For more details, you can check this old demo: Periodic homogenization of linear elastic materials — Numerical tours of continuum mechanics using FEniCS master documentation

2 Likes
5

Thank you so much, bleyerj!! Yes, the targeted application is elasticity with periodic boundary conditions, and this approach works great with dolfinx-mpc. The demo was very helpful!