Enforcing vertex locations to implement periodic boundaries

Question

I would like to enforce periodic boundary conditions on a mildly irregular domain (rectangle with one wall replaced with an arbitrary curve). The periodic boundary conditions are enforced on the two parallel, linear walls. In order to do this, I infer from links (1) and (2) that I need to match up the vertical coordinates of mesh vertices along the periodic boundaries.

I have two questions. First: How do I do this? This question was asked previously in (2) but not answered.

Second: am I currently making my mesh in a silly way? I currently generate my mesh this way:

# Define lower left corner
boundary = [Point(0.0,0.0)]

# Lower right corner
boundary.append( Point(Lx,0.0) )

# Curved upper surface
for i in range(1,Nx+1):
    ThisX = float(Nx-i)/float(Nx) * Lx
    boundary.append(Point (  ThisX  , SomeFunction(ThisX) ))

# Close polygon
boundary.append( Point(0.0,0.0) )

# Build mesh
domain = Polygon ( boundary )
mesh = generate_mesh(domain,Nx)

(1) https://fenicsproject.org/qa/3868/clarification-periodic-boundary-conditions-matching-vertices
(2) https://fenicsproject.org/qa/6616/integration-periodic-boundary-conditions-unpredictable

Answer 1

I'll answer my own question, or at least the first of my questions. I didn't realize that the Polygon() / generate_mesh combination automatically adds vertices at the vertices of the polygon. So I can generate a mesh with the same boundary points by feeding these points to Polygon(). This seems kind of sketchy, but it works...

# Build a mesh
def AddVerticalBoundaryVertices(l,x,y,n,dr):
    for i in range(1,n):
        if dr == -1: # Values go from (near) y to 0 with n points
            val = y * float(n - i)/float(n)
        elif dr == 1: # values go from (near) 0 to y with n points
            val = y * float(i)/float(n)
        l.append( Point(x, val) )

boundary = [Point(0.0,0.0)]
boundary.append( Point(Lx,0.0) )
AddVerticalBoundaryVertices(boundary,Lx,hy[-1],Ny,1)
for i in range(0,Nx+1):
    boundary.append( Point( hx[i] , hy[i] ) )
AddVerticalBoundaryVertices(boundary,0.0,hy[-1],Ny,-1)
boundary.append( Point(0.0,0.0) )