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)