dolfin team mailing list archive

[Martin Sandve Alnæs <785874@xxxxxxxxxxxxxxxxxx>]

On 31 May 2011 00:24, Anders Logg <logg@xxxxxxxxx> wrote:
> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>> There's two, don't remember what they do:
>>   def estimate_max_polynomial_degree(e, default_degree=1):
>>   def estimate_total_polynomial_degree(e, default_degree=1):
>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>>
>> ** Changed in: dolfin
>>        Status: New => Invalid
>
> And those include spatial coordinates?

Turns out they didn't. Just checked the code.
But it was easy to add. I'm commiting changes
to estimate_total_polynomial_degree now which
incorporate the spatial degree. Maybe this should
be used for assembling rhs and functionals, while
looking at elements are enough for the bilinear form?

PyDOLFIN could do something like

d = estimate_total_polynomial_degree(expr)
d = max(d, 1)
d = min(d, 8)

to limit the degree to some reasonable range in cases such as
expr = sin(x**5)*cos(y**5)
which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
Look at the code and tests in the last commit for more details, it's
quite short.

Martin

-- 
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:
  Invalid

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