← Back to team overview

dolfin team mailing list archive

Re: [Bug 785874] Re: Projection of x is not accurate

 

Ok, so NaN was because I manually specified too low degree.

So we're left with 'auto' which gives a big error:
0.00666666666667 0.993333333333 8.77744057992e-16

This is probably caused by ffc not considering
spatial coordinates in the expression. Fair enough.

So auto probably chooses 2 for the bilinear form
and 1 for the linear form, keeping the matrix nonsingular
but reducing the accuracy of the right hand side.

So the only bug here is that FFC could look at
spatial coordinates in the expression for better
autodetection, which would make project() of
expressions in PyDOLFIN more robust.
I'll close this bug and file a blueprint with this point.

Martin

On 30 May 2011 22:22, Garth Wells <785874@xxxxxxxxxxxxxxxxxx> wrote:
> The NaN is not caused by FFC autodetecting the degree.
>
> The numerical error in some cases is due to FFC not increasing the
> degree when the spatial coordinate is used as a coefficient function.
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/785874
>
> Title:
>  Projection of x is not accurate
>
> Status in DOLFIN:
>  New
>
> Bug description:
>  I've tested that projecting x works without the scaling bug that was
>  just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3
>  degrees. I print the max and min values of the vector of the
>  projection function, and the values are _close_ to 0 and 1 but not to
>  machine precision. The script is below.
>
>  There's up to 2.7% error in the 3D case. Is the projection form
>  integrated accurately enough? All but the DG0 function space should be
>  capable of representing x exactly. Not sure if this is a dolfin or ffc
>  bug.
>
>
>  from dolfin import *
>
>  def mcx(dim):
>      if dim == 1:
>          mesh = UnitInterval(20)
>          cell = interval
>          x = cell.x
>      if dim == 2:
>          mesh = UnitSquare(20, 20)
>          cell = triangle
>          x = cell.x[0]
>      if dim == 3:
>          mesh = UnitCube(20, 20, 20)
>          cell = tetrahedron
>          x = cell.x[0]
>      return mesh, cell, x
>
>  for dim in range(1, 4):
>      mesh, cell, x = mcx(dim)
>      minval, maxval = 1.0, 0.0
>      #print dim, "DG"
>      for degree in range(3):
>          V = FunctionSpace(mesh, "DG", degree)
>          u = project(x, V)
>          #print dim, degree, u.vector().min(), u.vector().max()
>          minval = min(minval, u.vector().min())
>          maxval = max(maxval, u.vector().max())
>      #print dim, "CG"
>      for degree in range(1, 4):
>          V = FunctionSpace(mesh, "CG", degree)
>          u = project(x, V)
>          #print dim, degree, u.vector().min(), u.vector().max()
>          minval = min(minval, u.vector().min())
>          maxval = max(maxval, u.vector().max())
>      print minval, maxval
>

-- 
You received this bug notification because you are a member of DOLFIN
Team, which is subscribed to DOLFIN.
https://bugs.launchpad.net/bugs/785874

Title:
  Projection of x is not accurate

Status in DOLFIN:
  New

Bug description:
  I've tested that projecting x works without the scaling bug that was
  just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3
  degrees. I print the max and min values of the vector of the
  projection function, and the values are _close_ to 0 and 1 but not to
  machine precision. The script is below.

  There's up to 2.7% error in the 3D case. Is the projection form
  integrated accurately enough? All but the DG0 function space should be
  capable of representing x exactly. Not sure if this is a dolfin or ffc
  bug.

  
  from dolfin import *

  def mcx(dim):
      if dim == 1:
          mesh = UnitInterval(20)
          cell = interval
          x = cell.x
      if dim == 2:
          mesh = UnitSquare(20, 20)
          cell = triangle
          x = cell.x[0]
      if dim == 3:
          mesh = UnitCube(20, 20, 20)
          cell = tetrahedron
          x = cell.x[0]
      return mesh, cell, x

  for dim in range(1, 4):
      mesh, cell, x = mcx(dim)
      minval, maxval = 1.0, 0.0
      #print dim, "DG"
      for degree in range(3):
          V = FunctionSpace(mesh, "DG", degree)
          u = project(x, V)
          #print dim, degree, u.vector().min(), u.vector().max()
          minval = min(minval, u.vector().min())
          maxval = max(maxval, u.vector().max())
      #print dim, "CG"
      for degree in range(1, 4):
          V = FunctionSpace(mesh, "CG", degree)
          u = project(x, V)
          #print dim, degree, u.vector().min(), u.vector().max()
          minval = min(minval, u.vector().min())
          maxval = max(maxval, u.vector().max())
      print minval, maxval


Follow ups

References