Hi all,
I have been trying to progress my code from a 2D Octagon with periodicity in x and y directions to a 3D truncated sphere with periodicity in x, y and z directions. The solution to my 2D code is indeed periodic in both directions but when I extend this to 3D, the solution is no longer periodic.
I have looked at this case:
https://fenicsproject.org/qa/262/possible-specify-more-than-one-periodic-boundary-condition
and this is how I implemented the 2D case which works.
I have looked at this case:
https://fenicsproject.org/qa/10700/modifying-periodicboundary3-class-to-be-periodic-in-x-y-z
but this simply brings me back to the previous post.
The way I define my boundary conditions for the 2D code is:
# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
# Making the negative parts the reference parts and the positive the slave parts
def inside(self, x, on_boundary):
return bool((near(x[0], -0.5*L) or near(x[1], -0.5*L)) and \
(not ((near(x[0], -0.5*L) and near(x[1], 0.5*L)) or \
(near(x[1], -0.5*L) and near(x[0], 0.5*L)))) and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
if (near(x[0], 0.5*L) and near(x[1], 0.5*L)):
y[0] = x[0] - 1.0*L
y[1] = x[1] - 1.0*L
elif near(x[0], 0.5*L):
y[0] = x[0] - 1.0*L
y[1] = x[1]
elif near(x[1], 0.5*L):
y[0] = x[0]
y[1] = x[1] - 1.0*L
else:
y[0] = 1000.*L
y[1] = 1000.*L
#-----------------
# Create geometry
#-----------------
L = 1.
# Unit Cell - (corner_1(x,y,z), corner_2(x,y,z))
square = Rectangle(Point(-.5*L, -.5*L), Point(.5*L, .5*L))
#Sphere - (Centre(x,y,z), radius)
circle = Circle(Point(0, 0), 0.6*L)
# Create anode geometry
quaffle = circle - (circle - square)
# Generate mesh
mesh = generate_mesh(quaffle, 32)
# Save mesh
File("2Dmesh.xml") << mesh
# Plot mesh and allow interaction
plot(mesh)
interactive()
#------------------------
# Boundary Conditions
#------------------------
bdry = FacetFunction("size_t", mesh)
bdry.set_all(30)
class Neumann(SubDomain):
def inside(self, x, on_boundary):
return bool(not (near(x[0], -0.5*L) or near(x[0], 0.5*L) or near(x[1], -0.5*L) or near(x[1], 0.5*L))) and on_boundary
Neumann().mark(bdry, 31);
ds = Measure("ds")[bdry]
bcs = []
and for the 3D code, this is done the following way:
# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
# Making the negative parts the reference parts and the positive the slave parts
def inside(self, x, on_boundary):
return bool((near(x[0], -0.5*L) or near(x[1], -0.5*L) or near(x[2], -0.5*L)) and \
(not ((near(x[0], -0.5*L) and near(x[1], 0.5*L)) or \
(near(x[0], -0.5*L) and near(x[2], 0.5*L)) or \
(near(x[1], -0.5*L) and near(x[0], 0.5*L)) or \
(near(x[1], -0.5*L) and near(x[2], 0.5*L)) or \
(near(x[2], -0.5*L) and near(x[0], 0.5*L)) or \
(near(x[2], -0.5*L) and near(x[1], 0.5*L)))) and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
if (near(x[0], 0.5*L) and near(x[1], 0.5*L) and near(x[2], 0.5*L)):
y[0] = x[0] - 1.0*L
y[1] = x[1] - 1.0*L
y[2] = x[2] - 1.0*L
elif (near(x[0], 0.5*L) and near(x[1], 0.5*L)):
y[0] = x[0] - 1.0*L
y[1] = x[1] - 1.0*L
y[2] = x[2]
elif (near(x[1], 0.5*L) and near(x[2], 0.5*L)):
y[0] = x[0]
y[1] = x[1] - 1.0*L
y[2] = x[2] - 1.0*L
elif (near(x[2], 0.5*L) and near(x[0], 0.5*L)):
y[0] = x[0] - 1.0*L
y[1] = x[1]
y[2] = x[2] - 1.0*L
elif near(x[0], 0.5*L):
y[0] = x[0] - 1.0*L
y[1] = x[1]
y[2] = x[2]
elif near(x[1], 0.5*L):
y[0] = x[0]
y[1] = x[1] - 1.0*L
y[2] = x[2]
elif near(x[2], 0.5*L):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - 1.0*L
else:
y[0] = 1000.*L
y[1] = 1000.*L
y[2] = 1000.*L
#-----------------
# Create geometry
#-----------------
L = 1.
# Unit Cell - (corner_1(x,y,z), corner_2(x,y,z))
box = Box(Point(-.5*L, -.5*L, -.5*L), Point(.5*L, .5*L, .5*L))
#Sphere - (Centre(x,y,z), radius)
sphere = Sphere(Point(0, 0, 0), 0.6*L)
# Create anode geometry
quaffle = sphere - (sphere - box)
# Generate mesh
mesh = generate_mesh(quaffle, 32)
# Save mesh
File("mesh.xml") << mesh
# Plot mesh and allow interaction
plot(mesh)
interactive()
#------------------------
# Boundary Conditions
#------------------------
bdry = FacetFunction("size_t", mesh)
bdry.set_all(30)
class electrolyte(SubDomain):
def inside(self, x, on_boundary):
return bool(not (near(x[0], -0.5*L) or near(x[0], 0.5*L) or near(x[1], -0.5*L) or \
near(x[1], 0.5*L) or near(x[2],-0.5*L) or near(x[2], 0.5*L)) and \
on_boundary)
electrolyte().mark(bdry, 31);
ds = Measure("ds")[bdry]
bcs = []
Is there something I am doing wrong to specify the periodic boundary condition or is it a badly posed problem and the periodicity is the first thing to go? Any help would be much appreciated.
For completeness, the code used to define the function spaces and constants for both codes is:
#------------------------
# Define function spaces
#------------------------
Vh = VectorFunctionSpace(mesh, "Lagrange", 1, constrained_domain=PeriodicBoundary())
u = TrialFunction(Vh)
v = TestFunction(Vh)
uh = Function(Vh)
#------------------------
# Constants
#------------------------
V_Si = 1.21E-5
V_C = 8.69E-6
c_maxSi = 3.65E5
c_maxC = 1.92E4
eta_Si = 0.1479
eta_C = 0.1933
nu_Si = 0.28
nu_C = 0.23
E_Si = 90.13
E_C = 27.6
lam_Si = (E_Si*nu_Si)/((1.+nu_Si)*(1.-2.*nu_Si))
lam_C = (E_C*nu_C)/((1.+nu_C)*(1.-2.*nu_C))
G_Si = E_Si/(2.*(1.+nu_Si))
G_C = E_C/(2.*(1.+nu_C))
lam = lam_Si/G_Si
G = G_Si/G_Si
and the solver is called by:
#------------------------
# Lam\'e Weak Form
#------------------------
f = Constant((-(lam+2.*G), -lam, -lam))
n = FacetNormal(mesh)
A1 = G * inner(grad(v),grad(u)+transpose(grad(u))) * dx
A2 = lam * div(u) * div(v) * dx
A = A1 + A2
b = inner(v,elem_mult(f,n)) * ds(31)
#------------------------
# Solve
#------------------------
solve(A == b, uh, bcs)
uh.rename("u","u")
File('output_periodic.pvd') << uh
for both of the codes.
Thank you!
