Periodic Boundary Condition for Elasticity Equation

General
1

Dear FEniCS Users,
I do hope you are doing great, and Happy Canada Day.
I would like to implement the periodic boundary condition for a simple 2D elasticity equation, then calculating the averaged-volume stress and compare it with the analytical solution. A minimal of my code implementation has been given below for your consideration:

import numpy as np
from dolfin import *

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool((near(x[0], 0) or near(x[1], 0)) and on_boundary)

    def map(self, x, y):
        if near(x[0], 1) and near(x[1], 1):
            y[0] = x[0] - 1.0
            y[1] = x[1] - 1.0
        elif near(x[0], 1):
            y[0] = x[0] - 1.0
            y[1] = x[1]
        else:
            y[0] = x[0]
            y[1] = x[1] - 1.0

# Create the periodic boundary condition
periodic_bc = PeriodicBoundary()
W = 1.0
nFEx = 1
nFEy = nFEx
NFE = nFEx*nFEy
NFE *= 4 #if crossed, else *2 
E=1.0
nu = 0.0
mu = Constant(E/(2*(1+nu)))
lambda_ = Constant(E*nu/((1+nu)*(1-2*nu)))
mesh = RectangleMesh(Point(0.0, 0.0), Point(W, W), nFEx, nFEy, "crossed")
V = VectorFunctionSpace(mesh, 'P', 1, constrained_domain=periodic_bc)
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(0)
bc = []
boundaries = MeshFunction("size_t", mesh, 1)
boundaries.set_all(0)
        
ds = Measure("ds", domain = mesh, subdomain_data=boundaries)
dx = Measure("dx", domain = mesh)

def point(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0]) < tol) and (abs(x[1]) < tol)

# Define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
    return lambda_*div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx
# Define time stepping parameters
t = 1.0
num_steps = 1
dt = t/num_steps

for n in range(num_steps):
    # Update BCs
    bc.append(DirichletBC(V.sub(0), Constant((0.)), point, method="pointwise"))
    bc.append(DirichletBC(V.sub(1), Constant((0.)), point, method="pointwise"))
    # Solve variational problem
    solve(a == L, u, bc)
    sigmaTensor = project(sigma(u),S)
    epsTensor = project(epsilon(u),S)
    # Compute volume averages
    for i in range(4):
        sigmaAVG[n+1,i] = assemble(sigmaTensor.sub(i)*dx)/W
        epsAVG[n+1,i] = assemble(epsTensor.sub(i)*dx)/W
    sigmaGP[n+1] = sigmaTensor.vector().get_local().reshape(-1,4)
    epsGP[n+1] = epsTensor.vector().get_local().reshape(-1,4)
    # sig
    # Write result file
    xdmffile.write(u,n)
    # # Update time
    t += dt

However, the results are not correct in terms of using periodic boundary conditions. Any help regarding this would be highly appreciated.
Regards,
Benhour

2

Refer to Periodic BC Problem in 3D - #6 by conpierce8.

It looks like you need to modify your inside function so that it returns False for the points marked 2 and 4 in my illustrations, and that you need to add a fourth condition to map which handles all points that are not part of the periodic boundary.