dolfin team mailing list archive

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

On 6 June 2011 12:33, Anders Logg <logg@xxxxxxxxx> wrote:
> On Mon, Jun 06, 2011 at 10:17:27AM -0000, Garth Wells wrote:
>> On 06/06/11 11:05, Martin Sandve Alnæs wrote:
>> > On 6 June 2011 11:54, Garth Wells <785874@xxxxxxxxxxxxxxxxxx> wrote:
>> >> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
>> >>> 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.
>> >>>
>> >>
>> >> We have the same issue of order blow-out for problems with lots of
>> >> coefficients. I'm therefore inclined not to include any heuristics, and
>> >> leave it up to the user.
>> >
>> > Do you mean we should actually crash and burn with this line?
>> >   f = assemble(sin(triangle.x[0]), mesh=mesh)
>> > With the current heuristic this will give degree 3,
>> > 1 from x and +2 from sin.
>> >
>>
>> We obviously need an approach for functions from non-polynomial spaces.
>> What I'm not inclined towards is arbitrary thresholds for integrating
>> polynomial products.
>
> What if we by default set that threshold to the maximal element degree
> + k, where k is say 1? It would be enough to retain expected
> convergence.
>
> --
> Anders

That's definitely not enough (e.g. L = x**8*v*dx) and would be basically
the same as throwing away the entire total polynomial degree algorithm.

It might be necessary/desirable to use different rules for bilinear
forms and linear forms/functionals?

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