Interpolate a function with periodic boundary condition

Hi all, I am new to fenics. Instead of solving a differential equation, I was trying to plot a gaussian function in a function space with 2D periodic boundary. But I was not getting the expected results when the center in the gaussian function changes eg.(0,0),(1,1),(0.5,0.5),(1,0),(0,1). In the code below, the gaussian function is applied at (1,0). Thanks for your help in advance.
My code:

from dolfin import *
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np



# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
    
    def __init__(self, b):
        SubDomain.__init__(self)
        self.b = b
        
    
    def inside(self, y, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool((near(y[0], 0) or near(y[1], 0)) and 
            (not ((near(y[0], self.b) and near(y[1], 0)) or 
                  (near(y[0], 0) and near(y[1], self.b)))) and on_boundary)

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

# Create mesh and finite element
mesh = UnitSquareMesh(32, 32)
b=1
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary(b))

#Gaussian function with center(1,0)
c_x=1
c_y=0
f = Expression("10*exp(-(pow(x[0] - c_x, 2) + pow(x[1] - c_y, 2)) / 0.02)", degree=2,c_x=c_x,c_y=c_y)

f = interpolate(f,V)
print(f)

plt.figure()
plot(f)

If you move the gaussian function to (0.8,0.9) it is easier to see whats going on (I have also plotted the mesh on top).
int

As you can observe here, the values are mapped correctly from the top right corner to the right side and bottom facets.

Note that that the interpolation is only occurring at the boundary dofs, and therefore the values rapidly decreases to the interpolated functions value for the dofs close to (but not on the boundary).

Code for this is supplied below:

from dolfin import *

import matplotlib.pyplot as plt
import numpy as np



# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    def __init__(self, b):
        SubDomain.__init__(self)
        self.b = b


    def inside(self, y, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool((near(y[0], 0) or near(y[1], 0)) and
            (not ((near(y[0], self.b) and near(y[1], 0)) or
                  (near(y[0], 0) and near(y[1], self.b)))) and on_boundary)

    def map(self, x, y):
        if (near(x[0], self.b) and near(x[1], self.b)):
            y[0] = x[0] - self.b
            y[1] = x[1] - self.b

        elif near(x[0], self.b):
            y[0] = x[0] - self.b
            y[1] = x[1]

        elif near(x[1], self.b):
            y[0] = x[0]
            y[1] = x[1] - self.b

        else:
            y[0] = -1000
            y[1] = -1000


# Create mesh and finite element
mesh = UnitSquareMesh(10, 10)
b=1
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary(b))

#Gaussian function with center(1,0)
c_x=0.8
c_y=0.9
f = Expression("10*exp(-(pow(x[0] - c_x, 2) + pow(x[1] - c_y, 2)) / 0.05)", degree=2,c_x=c_x,c_y=c_y)

f = interpolate(f,V)
plt.figure()
plot(f)
plot(mesh)
plt.savefig("int.png")

Thank you for your reply. But when I move the gaussian to (0.1,0.1), mapping is not happening. What could be done, if I want to interpolate a function in a periodic domain irrespective of the point of application of gaussian function?

Thank you for your help

That is because you are mapping the points (1,y)->(0,y), (x,1)->(x,0), (1,1)->(0,0). If the gaussian is centered at (0.1,0.1), the values that are mapped are close to zero (as (1,y), (x,1), (1,1) are far away from (0.1,0.1).

The direction of your periodic mapping will affect the result of an interpolation. If you want it to be independent of the direction, I would suggest using a projection (the fenics project function).

As you have rightly pointed out, I’m mapping (1,y)->(0,y), (x,1)->(x,0), but still I’m expecting periodicity independent of direction.

I’m seeing the following result after the interpolation of gaussian function at (0,0.5), which clearly shows no periodicity.
interpolate

Then when I do the projection of gaussian function at (0,0.5), I see the following result which shows periodicity but only at border as a flat structure.
project

But what I’m expecting is somewhat like the below picture, which shows some continuity interior to boundary
poisson

As the gaussian function at (0.5,0.5) shows the following structure.
center

Is there anyway where I can implement periodicity of a function independent of the direction of periodicity defined and getting the same structure of the function?

Thanks for the help.

I am not sure why you expect the third picture.
If we start by looking at the top picture, it is perfectly clear why it does not show periodicity, as you are mapping the values at (1,y) to (0, y) (evaluating your gaussian function), which is close to 0.

For the second figure, you observe periodicity, as for a projection, the direction of the mapping does not matter. However, as you are trying to fit a function that is not periodic to a periodic function, you will get some result.
If we just consider a 1D example, to simplify the problem:

from dolfin import *

import matplotlib.pyplot as plt
import numpy as np


class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        if np.isclose(x[0], 1):
            y[0] = x[0] - 1.0
        else:
            y[0] = x[0]

# Create mesh and finite element
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())

#Gaussian function with center(1,0)
c_x=0
exp = Expression("x[0]", degree=1)

# Interpolation
f = interpolate(exp,V)
f.rename("interpolation", "")
File("f_int.pvd") << f

# Projection
f = project(exp, V)
f.rename("projection", "")
File("f_project.pvd") << f

# Original function
V_ex = FunctionSpace(mesh, "CG", 1)
f = interpolate(exp, V_ex)
f.rename("exact", "")
File("f_ex.pvd") << f

which we can visualize with for instance “Plot over line” in Paraview, and obtain:


As you observe, when mapping from 1->0, the value at f(0) = f(1) and all other values are equal to the original function f(x)=x. For the projection, we minimize \int (u_h - x)\mathrm{d} x, we therefore observe overshoot and undershoots compared to the original function.

We can of course do this with another function, for instance 1-x, and obtain the following plot: