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

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!

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-03-03
Last reply:
2013-03-03
Anders Logg (logg) said : #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:
> https://answers.launchpad.net/dolfin/+question/223087
>
> 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 : #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 : #3

On Thu, Feb 28, 2013 at 08:45:59PM -0000, Xiaoxian Liu wrote:
> Question #223087 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223087
>
> Status: Answered => Open
>
> 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 : #4

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

Xiaoxian Liu (liuxiaox) said : #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 : #6

On Thu, Feb 28, 2013 at 09:16:02PM -0000, Xiaoxian Liu wrote:
> Question #223087 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223087
>
> Xiaoxian Liu gave more information on the question:
> 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 : #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 : #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:
https://launchpad.net/~fenics-packages/+archive/fenics-dev

The problem would come back if I upgrade the package following the instruction at:
http://fenicsproject.org/download/ubuntu_details.html#ubuntu-details

Is it a bug by any chance?

Garth Wells (garth-wells) said : #9

On 3 March 2013 01:41, Xiaoxian Liu
<email address hidden> wrote:
> Question #223087 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223087
>
> Xiaoxian Liu gave more information on the question:
> 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:
> https://launchpad.net/~fenics-packages/+archive/fenics-dev
>
> The problem would come back if I upgrade the package following the instruction at:
> http://fenicsproject.org/download/ubuntu_details.html#ubuntu-details
>
> 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 : #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 : #11

On 3 March 2013 14:25, Xiaoxian Liu
<email address hidden> wrote:
> Question #223087 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223087
>
> Status: Answered => Open
>
> 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.

Can you help with this problem?

Provide an answer of your own, or ask Xiaoxian Liu for more information if necessary.

To post a message you must log in.