Dear All,
I’m new to fenics and I’m struggling with apparently simple problem. I need to solve 2D Poisson equation with periodic boundary conditions on both X and Y and with the rhs loaded from file.
I assembled together several examples and snippets from the forum. Here is my code:
from dolfin import *
import numpy as np
from scipy.interpolate import interp2d
# Load data stored as NxN matrix
data = np.loadtxt('q_distrib.dat')
f = np.array(data)
N = f.shape
print('N=',N)
# Need the format x,y,value for interpolation
# I know that this code is not optimal, nevr mind
x = np.zeros(N[0]*N[0])
y = np.zeros(N[0]*N[0])
val = np.zeros(N[0]*N[0])
dx = 1.0/float(N[0])
n = 0
for i in range(N[0]):
for j in range(N[1]):
x[n] = 0.5*dx + dx*i
y[n] = 0.5*dx + dx*j
val[n] = data[i,j]
n+=1
interpolant = interp2d(x, y, val, kind='linear', copy=False, bounds_error=True)
# Define source function (taken from forums)
class Source(Expression):
def __init__(self, f, **kwargs):
self._f = f
UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
values[:] = self._f(*x)
# Sub domain for Periodic boundary conditions over x and y (taken from forums)
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 define function space
mesh = UnitSquareMesh(32,32)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())
f = Source(interpolant, element=V.ufl_element())
f = interpolate(f, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, [])
# Save solution in VTK format
file = File("poisson.pvd")
file << u
# Plot solution
import matplotlib.pyplot as plt
plot(u)
plt.show()
When I run this code it shows no error but it then runs forever (waited for an hour) and produces no output.
I can’t figure out what is wrong since no errors are reported… Any suggestions?