Subject: Re: [Ffc] [Branch ~ffc-core/ffc/dev] Rev 1586: Return values for QuadratrueElement.tabulate().

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 29 Jan 2010 16:13:13 +0100 (CET)

I think it is almost there, but we need to rethink the way we manipulate the degrees of finite elements and forms, I think we have discussed this before but I don't recall if we ever reached a conclusion.
Try running the QuadratureElement.ufl demo.

Kristian

On 29 January 2010 16:04,  <[email-masked]> wrote:
> ------------------------------------------------------------
> revno: 1586
> committer: Kristian B. Ølgaard <[email-masked]>
> branch nick: dev
> timestamp: Fri 2010-01-29 16:02:09 +0100
> message:
>  Return values for QuadratrueElement.tabulate().
> modified:
>  ffc/quadratureelement.py
>
>
> --
> lp:~ffc-core/ffc/dev
> https://code.launchpad.net/~ffc-core/ffc/dev
>
> You are subscribed to branch lp:~ffc-core/ffc/dev.
> To unsubscribe from this branch go to https://code.launchpad.net/~ffc-core/ffc/dev/+edit-subscription.
>
> === modified file 'ffc/quadratureelement.py'
> --- ffc/quadratureelement.py    2010-01-29 14:51:42 +0000
> +++ ffc/quadratureelement.py    2010-01-29 15:02:09 +0000
> @@ -13,14 +13,7 @@
>  from FIAT.functional import PointEvaluation
>
>  # FFC modules.
> -from log import error
> -#from finiteelement import ufl_domain2fiat_domain
> -#from dofrepresentation import DofRepresentation
> -#from ffcquadraturerules import make_quadrature
> -#from finiteelement import FiniteElement
> -#from finiteelement import AFFINE
> -#from finiteelement import CONTRAVARIANT_PIOLA
> -#from finiteelement import COVARIANT_PIOLA
> +from log import error, info_red
>
>  # Default quadrature element degree
>  default_quadrature_degree = 1
> @@ -62,7 +55,7 @@
>
>         # Save the geometric dimension.
>         # FIXME: KBO: Do we need to change this in order to integrate on facets?
> -        self.geometric_dimension = ufl_element.cell().geometric_dimension
> +        self._geometric_dimension = ufl_element.cell().geometric_dimension()
>
>     def space_dimension(self):
>         "The element space dimension is simply the number of quadrature points"
> @@ -101,8 +94,8 @@
>         # RuntimeError again.
>         if order:
>             # error("Derivatives are not defined on a QuadratureElement")
> -            print "\n*** WARNING: Derivatives are not defined on a QuadratureElement,"
> -            print   "             returning values of basisfunction.\n"
> +            info_red("\n*** WARNING: Derivatives are not defined on a QuadratureElement,")
> +            info_red("             returning values of basisfunction.\n")
>
>         # Check that incoming points are equal to the quadrature points.
>         if len(points) != len(self._points) or abs(numpy.array(points) - self._points).max() > 1e-12:
> @@ -110,11 +103,10 @@
>             print "\nquad points:\n", self._points
>             error("Points must be equal to coordinates of quadrature points")
>
> -#        # Return the identity matrix of size __num_quad_points in a
> -#        # suitable format for monomialintegration.
> -#        values = numpy.identity(self._num_quad_points, float)
> -#        table = [{(0,)*self.cell().topological_dimension(): values}]
> -#        return table
> +        # Return the identity matrix of size len(self._points) in a
> +        # suitable format for tensor and quadrature representations.
> +        values = numpy.eye(len(self._points))
> +        return {(0,)*self._geometric_dimension: values}
>
>  def _create_entity_dofs(cell, num_dofs):
>     "This function is ripped from FIAT/discontinuous_lagrange.py"
>
>
>

---

From: Anders Logg <[email-masked]>
Date: Fri, 29 Jan 2010 20:21:23 +0100

On Fri, Jan 29, 2010 at 04:13:13PM +0100, Kristian Oelgaard wrote:
> I think it is almost there, but we need to rethink the way we
> manipulate the degrees of finite elements and forms, I think we have
> discussed this before but I don't recall if we ever reached a
> conclusion.  Try running the QuadratureElement.ufl demo.

I don't see exactly what the problem is. Is it the automatic degree
selection?

I've added some more debug output and the following is now printed for
the QuadratureElement demo:

  Found quadrature element(s) with the following degree(s): [3]
  Selecting quadrature degree based on quadrature element: 3
  quadrature_degree: auto --> 3

That looks correct to me. Shouldn't this result in the same points
being chosen for both the QuadratureElement and the quadrature loop?

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 29 Jan 2010 21:08:44 +0100 (CET)

> being chosen for both the QuadratureElement and the quadrature loop?

I see now that it is only a problem for tensor representation. For the linear form the first term v*u*dx makes it a form of degree 4, and tensor representation will use that degree for the entire form if I remember correctly?

Kristian

> --
> Anders
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAktjNTMACgkQTuwUCDsYZdHwxACcCTzBr0hzp1nXowRvQKnIsrmg
> /zAAn36IpBEEvbxC2lVEh1y0TP+F7eDZ
> =/PfD
> -----END PGP SIGNATURE-----
>
>

---

From: Anders Logg <[email-masked]>
Date: Fri, 29 Jan 2010 21:12:07 +0100

On Fri, Jan 29, 2010 at 09:08:44PM +0100, Kristian Oelgaard wrote:
> I see now that it is only a problem for tensor representation. For
> the linear form the first term v*u*dx makes it a form of degree 4,
> and tensor representation will use that degree for the entire form
> if I remember correctly?

Yes, and that's not much of a problem since it happens at compile
time.

But for the FFC demo the representation is chosen as quadrature since
it involves quadrature elements, so if it fails it shouldn't be
related to tensor representation?

--
Anders

---

From: Kristian Oelgaard <[email-masked]>
Date: Fri, 29 Jan 2010 21:22:51 +0100 (CET)

> Yes, and that's not much of a problem since it happens at compile

It is a problem if you try to tabulate a QuadratureElement with the wrong number of points :)

> related to tensor representation?

Yes, for the bilinear form, for the linear form tensor is chosen, the degree is set to 4 and the problems begin.

 File "/home/oelgaard/software/fenics_branches/dev/ffc/tensor/monomialintegration.py", line 101, in _init_table
    table[(ufl_element, None)] = fiat_element.tabulate(order, points)
  File "/home/oelgaard/software/fenics_branches/dev/ffc/mixedelement.py", line 77, in tabulate
    table = element.tabulate(order, points)
  File "/home/oelgaard/software/fenics_branches/dev/ffc/quadratureelement.py", line 104, in tabulate
    error("Points must be equal to coordinates of quadrature points")
  File "<string>", line 1, in <lambda>
  File "/home/oelgaard/software/fenics/ufl/ufl/log.py", line 124, in error
    raise UFLException(self._format_raw(*message))
ufl.log.UFLException: Points must be equal to coordinates of quadrature points

Kristian

> iEYEARECAAYFAktjQRYACgkQTuwUCDsYZdHxBQCfUtAbBw0wXjBsiyIzu8wbzc8p
> HzQAn1I8j9IG7LUPVtK7IgShyIHEtvie
> =4nxQ

---

From: Anders Logg <[email-masked]>
Date: Fri, 29 Jan 2010 21:31:18 +0100

> ufl.log.UFLException: Points must be equal to coordinates of quadrature points

ok, I missed that! The automatic selection of representation should
check for quadrature elements and in that case use quadrature.

I'll fix it.

--
Anders