# creating DG-representation of a step function

Asked by Xiaoxian Liu on 2013-02-28

Hi all,

I was trying to define a step function over a rectangular region. It was supposed to be 0 on the left side and 1 on the right; these values coincide with my choice of mesh function indicating the two subdomains. However, when the result comes out, it's 1 on the left and 0 on the right... I couldn't solve it for a quite a while. Please take a look at the code below:

###### beginning of code
from dolfin import *

# a simple mesh consisting of 4 triangles
mesh = RectangleMesh(-1, 0, 1, 1, 2, 1)

class LeftRegion(SubDomain):
def inside(self, x, on_boundary):
return x[0]<= DOLFIN_EPS
class RightRegion(SubDomain):
def inside(self, x, on_boundary):
return x[0]>= -DOLFIN_EPS

# mark the left side as "0" and right side as "1"
marker = CellFunction('size_t', mesh, 0)
RightRegion().mark(marker, 1)

# DG-representation of marker
Vdg = FunctionSpace(mesh, 'DG', 0)
dx = Measure("dx")[marker]

u = TrialFunction(Vdg)
v = TestFunction(Vdg)
a = u * v * dx(0) + u*v*dx(1)
L = Constant(0)*v*dx(0) + Constant(1)*v*dx(1)

u = Function(Vdg)
solve(a==L, u)mesh, subdomains, domain markers

print "marker is", marker.array()
print "u is", u.vector().array()

########### end of code

marker is [0 0 1 1]
u is [ 1. 1. 0. 0.] (u is also supposed to be [0 0 1 1]!!!!)

Thank you!

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-03-03
2013-03-03
 Anders Logg (logg) said on 2013-02-28: #1

When I run the code I get

marker is [0 0 1 1]
u is [ 0. 0. 1. 1.]

--
Anders

On Thu, Feb 28, 2013 at 04:50:56PM -0000, Xiaoxian Liu wrote:
> New question #223087 on DOLFIN:
>
> Hi all,
>
> I was trying to define a step function over a rectangular region. It was supposed to be 0 on the left side and 1 on the right; these values coincide with my choice of mesh function indicating the two subdomains. However, when the result comes out, it's 1 on the left and 0 on the right... I couldn't solve it for a quite a while. Please take a look at the code below:
>
> ###### beginning of code
> from dolfin import *
>
> # a simple mesh consisting of 4 triangles
> mesh = RectangleMesh(-1, 0, 1, 1, 2, 1)
>
> class LeftRegion(SubDomain):
> def inside(self, x, on_boundary):
> return x[0]<= DOLFIN_EPS
> class RightRegion(SubDomain):
> def inside(self, x, on_boundary):
> return x[0]>= -DOLFIN_EPS
>
> # mark the left side as "0" and right side as "1"
> marker = CellFunction('size_t', mesh, 0)
> RightRegion().mark(marker, 1)
>
>
> # DG-representation of marker
> Vdg = FunctionSpace(mesh, 'DG', 0)
> dx = Measure("dx")[marker]
>
> u = TrialFunction(Vdg)
> v = TestFunction(Vdg)
> a = u * v * dx(0) + u*v*dx(1)
> L = Constant(0)*v*dx(0) + Constant(1)*v*dx(1)
>
> u = Function(Vdg)
> solve(a==L, u)mesh, subdomains, domain markers
>
> print "marker is", marker.array()
> print "u is", u.vector().array()
>
> ########### end of code
>
> Results of the code reads:
> marker is [0 0 1 1]
> u is [ 1. 1. 0. 0.] (u is also supposed to be [0 0 1 1]!!!!)
>
>
> Thank you!
>

 Xiaoxian Liu (liuxiaox) said on 2013-02-28: #2

Thanks for the reply, Anders. But this is impossible... I still got the wrong result... What could possibly go wrong in this case?
(I could define the DG function cell by cell, but I couldn't avoid this issue when I project a function to a different function space...)

 Anders Logg (logg) said on 2013-02-28: #3

On Thu, Feb 28, 2013 at 08:45:59PM -0000, Xiaoxian Liu wrote:
> Question #223087 on DOLFIN changed:
>
>
> Xiaoxian Liu is still having a problem:
> Thanks for the reply, Anders. But this is impossible... I still got
> the wrong result... What could possibly go wrong in this case?

Try plotting the function. You can't be sure that the dofs are
numbered from left to right.

--
Anders

> (I could define the DG function cell by cell, but I couldn't avoid this issue when I project a function to a different function space...)
>

 Xiaoxian Liu (liuxiaox) said on 2013-02-28: #4

The graph shows "1" on the left square and "0" on the right square

 Xiaoxian Liu (liuxiaox) said on 2013-02-28: #5

I recently upgraded the package to 1.1.0.
Could it be something went wrong when I did the upgrade? I basically followed the instruction on fenicsproject.org
I shall try this code on a different computer to see what the result is.

 Anders Logg (logg) said on 2013-03-01: #6

On Thu, Feb 28, 2013 at 09:16:02PM -0000, Xiaoxian Liu wrote:
> Question #223087 on DOLFIN changed:
>
> I recently upgraded the package to 1.1.0.
> Could it be something went wrong when I did the upgrade? I basically followed the instruction on fenicsproject.org
> I shall try this code on a different computer to see what the result is.

I'm using the development version (1.1+).

--
Anders

 Xiaoxian Liu (liuxiaox) said on 2013-03-02: #7

I've tested the code on another computer and I received the same wrong result... Any possible solutions? I will not be able to do anything if the marker of the domain is not reliable... Thanks!

 Xiaoxian Liu (liuxiaox) said on 2013-03-03: #8

After I degraded my fenics to 1.0.0(the default one coming with Ubuntu and I'm using Ubuntu 12.10), the problem is gone and I received the result I expected, i.e. [0,0,1,1]. I followed the instruction at:

The problem would come back if I upgrade the package following the instruction at:

Is it a bug by any chance?

 Garth Wells (garth-wells) said on 2013-03-03: #9

On 3 March 2013 01:41, Xiaoxian Liu
> Question #223087 on DOLFIN changed:
>
> After I degraded my fenics to 1.0.0(the default one coming with Ubuntu and I'm using Ubuntu 12.10), the problem is gone and I received the result I expected, i.e. [0,0,1,1]. I followed the instruction at:
>
> The problem would come back if I upgrade the package following the instruction at:
>
> Is it a bug by any chance?
>

No.

You're making an assumption on the dof ordering, which is not
guaranteed and has changed.

Garth

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Xiaoxian Liu (liuxiaox) said on 2013-03-03: #10

Thank you for the reply, Garth.
I certainly shouldn't make any assumption on the dof ordering as you said, but when I plot the solution, it should give me the right answer, i.e. 0 on the left and 1 on the right. But FEniCs 1.1.0 gave me the wrong solution which is 1 on the left and 0 on the right, whereas FEniCS 1.0.0 gave me the correct plot with 0 on the left and 1 on the right.

 Garth Wells (garth-wells) said on 2013-03-03: #11

On 3 March 2013 14:25, Xiaoxian Liu
> Question #223087 on DOLFIN changed:
>
>
> Xiaoxian Liu is still having a problem:
> Thank you for the reply, Garth.
> I certainly shouldn't make any assumption on the dof ordering as you said, but when I plot the solution, it should give me the right answer, i.e. 0 on the left and 1 on the right. But FEniCs 1.1.0 gave me the wrong solution which is 1 on the left and 0 on the right, whereas FEniCS 1.0.0 gave me the correct plot with 0 on the left and 1 on the right.
>

Your question was on the output of:

print "u is", u.vector().array()

which does depend on the dof ordering. In your output, the array is
flipped, which is exactly what I would expect the new ordering to do
for P0 elements.

Garth

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.