dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Is there some UFL function that could be used instead to select the
degree (something that includes spatial coordinates)?

--
Anders


On Mon, May 30, 2011 at 08:44:15PM -0000, Martin Sandve Alnæs wrote:
> 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.
> >
> >
> > 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