How to generate FEniCS friendly periodic meshes?

1

I need help trying to solve a periodic poisson problem on nonsquare mesh. I used my own mesher to create a symmetric mesh. The problem runs on fine on a square mesh and gives a periodic solution however when I use my own mesh the solution is not periodic.

Thanks!

MWE:

import matplotlib.pyplot as plt
from dolfin import *

# Source term
class Source(UserExpression):
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = x[0]*sin(5.0*DOLFIN_PI*x[1]) \
                    + 1.0*exp(-(dx*dx + dy*dy)/0.02)
class PeriodicBoundary(SubDomain):
    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
        return bool((near(x[0], 0) or near(x[1], 0,)) and 
                (not ((near(x[0], 0) and near(x[1], 1)) or 
                        (near(x[0], 1) and 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.
            y[1] = x[1] - 1.
        elif near(x[0], 1):
            y[0] = x[0] - 1.
            y[1] = x[1]
        else:   # near(x[1], 1)
            y[0] = x[0]
            y[1] = x[1] - 1.


# Create mesh and finite element
mesh = UnitSquareMesh(64,64)
#mesh = Mesh("wiggle_mesh.xml.gz")
W = FiniteElement("Lagrange",mesh.ufl_cell(),1)
R = FiniteElement("R",mesh.ufl_cell(),0)
TH = W * R
V = FunctionSpace(mesh, TH, constrained_domain=PeriodicBoundary())

# No Dirichlet Boundary Condtions
bcs= []

# Define variational problem
(u,c) = TrialFunction(V)
(v,d) = TestFunction(V)
f = Source(degree=1)
# Bilinear with integral constraints on the solution
a = dot(grad(u), grad(v))*dx+c*v*dx+u*d*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)
(ss,c) = u.split()
# Save solution to file
file = File("periodic.pvd")
file << ss

# Plot solution
plot(ss)
plt.show()

holes
holes.png1085Ă—806 12.9 KB
square
square.png1085Ă—806 17.5 KB

2

Consider using gmsh. An example geometry to get you started is

// periodic_.geo
size = 0.1;

SetFactory("Built-in");

Point(1) = {0, 0, 0, size};
Point(2) = {1, 0, 0, size/2};
Point(3) = {1, 1, 0, size/4};
Point(4) = {0, 1, 0, size/6};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

cx = 0.8;
cy = 0.6;
Point(5) = {cx, cy, 0, size};
Point(6) = {cx-0.1, cy, 0, size/10};
Point(7) = {cx, cy-0.1, 0, size/10};
Point(8) = {cx+0.1, cy, 0, size/10};
Point(9) = {cx, cy+0.1, 0, size/10};

Circle(5) = {6, 5, 7};
Circle(6) = {7, 5, 8};
Circle(7) = {8, 5, 9};
Circle(8) = {9, 5, 6};

// Build mapping left-right   4|   ^
//                             v   |2
Periodic Line {2} = {-4};

Line Loop(1) = {1, 2, 3, 4};
Line Loop(2) = {5, 6, 7, 8};

Plane Surface(1) = {1, 2};

Physical Line(1) = {2};  // master-slave entities
Physical Line(2) = {4};
Physical Surface(1) = {1};

Run gmsh -2 periodic_.geo -format msh2 to generate the mesh. Below I assume you will use dolfin-convert periodic_.msh periodic_.xml to get XML files for mesh and entity values. Mesh can be checked for periodicity e.g. as follows

from dolfin import Mesh, MeshFunction, SubsetIterator
from scipy.spatial.distance import cdist
import numpy as np

mesh = Mesh('periodic_.xml')
surfs = MeshFunction('size_t', mesh, 'periodic__facet_region.xml')

mesh.init(1)
mesh.init(1, 0)
# Get all master vertices for marked surfaces
master = set(sum(map(list, (f.entities(0) for f in SubsetIterator(surfs, 1))), []))
# Get all slave vertices
slave = set(sum(map(list, (f.entities(0) for f in SubsetIterator(surfs, 2))), []))

x = mesh.coordinates()
# For every master vertex there is one master == shift(slave)
master = x[list(master)]
slave = x[list(slave)]
# So when we compute distance between points there will be exactly one
# zero in the distance matrix
shift = np.array([1, 0])

C = cdist(master, shift + slave)
assert all(np.count_nonzero(row) == len(row)-1 for row in C)
2 Likes
3

Great! Thanks.
I have used gmsh a little bit in the past and the example will help.

A couple of quick questions on the test script:
Where does “mesh.init(1)” come into play?
From what I have read the code is intializing all edges, and then determining what vertices these edges are connected to but I don’t see it being used.

Also without having a facet_region file I can check for periodicity by creating my own MeshFunction and applying it to the different subdomains on my mesh and in theory the assertion should fail for my mesh, right?

Once again thanks MiroK, I appreciate the help.

4

Edge-vertex connectivity is needed in f.entities(0) (get entities of 0 topological dimension of f).

Correct.

2 Likes