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