Dear All,
I am trying to reproduce a Onelab (GetDP) problem using Fenics. The problem to solve in the unit square is:
\begin{equation}
\begin{cases}\tag{1}
-\Delta u + u = f & \text{in } \Omega,\
\displaystyle{\frac{\partial u}{\partial \mathbf{n}} = 0} & \text{on }\partial\Omega,
\end{cases}
\end{equation}
with the corresponding weak form:
\begin{equation}\tag{2}
\int_{\Omega} \nabla u\cdot\nabla v \;{\rm d}\Omega + \int_{\Omega}u\, v \;{\rm d}\Omega - \int_{\Omega}fv\;{\rm d}\Omega = 0
\end{equation}
I thought that the corresponding Fenics code would be:
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(64, 64)
mesh.init()
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("(1+2pipi)cos(pix[0])cos(pix[1])")
a = (dot(grad(u), grad(v)) + cv + du) * dx
L = fvdx - uvdx
w = Function(W)
solve(a == L, w)
(u, c) = w.split()
plot(u, interactive=True)
output to Gmsh .pos file
coor = mesh.coordinates()
u_nodal_values = u.vector()
u_array = u_nodal_values.array()
fp = open('test_laplace.pos', 'w')
fp.write('View "utotm" {\n')
for facet in facets(mesh):
v0 = facet.entities(0)[0]
v1 = facet.entities(0)[1]
v2 = facet.entities(0)[2]
x0 = coor[v0][0]
y0 = coor[v0][1]
z0 = coor[v0][2]
x1 = coor[v1][0]
y1 = coor[v1][1]
z1 = coor[v1][2]
x2 = coor[v2][0]
y2 = coor[v2][1]
z2 = coor[v2][2]
fp.write("ST(% .13f,% .13f,% .13f,% .13f,% .13f,% .13f,% .13f,% .13f,% .13f){% .13f, % .13f, % .13f };\n" %(x0, y0, z0, x1, y1, z1, x2, y2, z2, u_array[v0], u_array[v1], u_array[v2] ) )
fp.write('};\n')
print 'Gmsh file test_laplace.pos has been written'
Unfortunately, the second term in (L = fvdx - uvdx) leads to the following error message (last line only):
ufl.log.UFLException: All terms in form must have same rank.
Could someone explain what is wrong in this formulation ?
Thanks.
