Poisson on a Moebius strip

Question

I am trying to run the Poisson demo on the Moebius strip as it was described here: http://fenicsproject.org/pub/workshops/fenics13/slides/Rognes.pdf.

The python script gives exactly the same results as it is in the pdf. However the cpp code produces a different result.
The ufl file:

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Constant(triangle)

a = inner(grad(u), grad(v))*dx
L = f*v*dx

and the cpp file:

#include <dolfin.h>
#include "Poisson.h"

using namespace dolfin;

class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  { return on_boundary; }
};

int main()
{
  Mesh mesh("moebius2.xml.gz");
  Poisson::FunctionSpace V(mesh);

  Constant u0(0.0);
  DirichletBoundary boundary;
  DirichletBC bc(V, u0, boundary);

  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);

  Constant f(1.0);
  L.f = f;

  Function u(V);
  solve(a == L, u, bc);

  plot(u);
  interactive();
  return 0;
}

This is the resulting image:

The python code:

from  dolfin  import *
mesh = Mesh("moebius2.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 1)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, "on_boundary")

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.0)
a = inner(grad(u), grad(v))*dx
L = f*v*dx

u = Function(V)
solve(a == L, u, bc)

plot(u, interactive=True)

The resulting image:

What is the reason for the different results? How should I modify the cpp code to get the same result as in the python version?

Answer 1

Hi, consider defining the finite element cell as

cell = Cell('triangle', 3)
element = FiniteElement("Lagrange", cell, 1)
# rest of your code

This way your FE space is setup from correct cells, i.e. triangles(topological dimension 2) in 3d(geometric dimension 3). You can see that this is the cell that is used in python code e.g. by print V.ufl_cell().