Subject: Re: [Ffc] [Branch ~ffc-core/ffc/main] Rev 1498: Added support for CellVolume from UFL.

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 7 Jul 2010 18:26:20 +0100

Supporting CellVolume makes it possible to do:

CG = FiniteElement("Lagrange", triangle, 2)
DG = FiniteElement("DG", triangle, 0)
v = TestFunction(DG)
f = Coefficient(CG)
vol = triangle.v
L = 1.0/vol*v*f*dx

instead of using a Coefficient for 'vol' and then defining an
Expression in DOLFIN.

Currently, it will use the absolute value of the determinant of the
Jacobian in tabulate_tensor ( std::abs(detJ) ).

Kristian

On 7 July 2010 18:19,  <[email-masked]> wrote:
> ------------------------------------------------------------
> revno: 1498
> committer: Kristian B. Ølgaard <[email-masked]>
> branch nick: ffc
> timestamp: Wed 2010-07-07 18:17:07 +0100
> message:
>  Added support for CellVolume from UFL.
> modified:
>  ChangeLog
>  ffc/quadrature/optimisedquadraturetransformer.py
>  ffc/quadrature/quadraturetransformer.py
>  ffc/quadrature/quadraturetransformerbase.py
>
>
> --
> lp:ffc
> https://code.launchpad.net/~ffc-core/ffc/main
>
> Your team FFC Core Team is subscribed to branch lp:ffc.
> To unsubscribe from this branch go to https://code.launchpad.net/~ffc-core/ffc/main/+edit-subscription
>
> === modified file 'ChangeLog'
> --- ChangeLog   2010-07-01 17:43:36 +0000
> +++ ChangeLog   2010-07-07 17:17:07 +0000
> @@ -1,3 +1,4 @@
> + - Added support for new geometric quantity CellVolume in UFL.
>  0.9.3 [2010-07-01]
>  - Make global_dimension for Real return an int instead of double, bug # 592088
>  - Add support for facet normal in 1D.
>
> === modified file 'ffc/quadrature/optimisedquadraturetransformer.py'
> --- ffc/quadrature/optimisedquadraturetransformer.py    2010-07-01 17:41:40 +0000
> +++ ffc/quadrature/optimisedquadraturetransformer.py    2010-07-07 17:17:07 +0000
> @@ -187,7 +187,7 @@
>         return {():new_val}
>
>     # -------------------------------------------------------------------------
> -    # FacetNormal (geometry.py).
> +    # FacetNormal, CellVolume (geometry.py).
>     # -------------------------------------------------------------------------
>     def facet_normal(self, o,  *operands):
>         #print("Visiting FacetNormal:")
> @@ -216,6 +216,17 @@
>
>         return {(): create_symbol(normal_component, GEO)}
>
> +    def cell_volume(self, o,  *operands):
> +        # Safety check.
> +        ffc_assert(not operands, "Didn't expect any operands for FacetNormal: " + repr(operands))
> +
> +        # FIXME: KBO: This has to change for higher order elements
> +        detJ = format["det(J)"](self.restriction)
> +        volume = format["absolute value"](detJ)
> +        self.trans_set.add(detJ)
> +
> +        return {():create_symbol(volume, GEO)}
> +
>     def create_argument(self, ufl_argument, derivatives, component, local_comp,
>                   local_offset, ffc_element, transformation, multiindices):
>         "Create code for basis functions, and update relevant tables of used basis."
>
> === modified file 'ffc/quadrature/quadraturetransformer.py'
> --- ffc/quadrature/quadraturetransformer.py     2010-06-30 09:59:49 +0000
> +++ ffc/quadrature/quadraturetransformer.py     2010-07-07 17:17:07 +0000
> @@ -262,7 +262,7 @@
>         return {():f_abs(operands[0][()])}
>
>     # -------------------------------------------------------------------------
> -    # FacetNormal (geometry.py).
> +    # FacetNormal, CellVolume (geometry.py).
>     # -------------------------------------------------------------------------
>     def facet_normal(self, o,  *operands):
>         #print("Visiting FacetNormal:")
> @@ -291,6 +291,18 @@
>
>         return {():normal_component}
>
> +    def cell_volume(self, o,  *operands):
> +        # Safety check.
> +        ffc_assert(not operands, "Didn't expect any operands for CellVolume: " + repr(operands))
> +
> +        # FIXME: KBO: This has to change for higher order elements
> +        detJ = format["det(J)"](self.restriction)
> +        volume = format["absolute value"](detJ)
> +        self.trans_set.add(detJ)
> +
> +        return {():volume}
> +
> +
>     def create_argument(self, ufl_argument, derivatives, component, local_comp,
>                   local_offset, ffc_element, transformation, multiindices):
>         "Create code for basis functions, and update relevant tables of used basis."
>
> === modified file 'ffc/quadrature/quadraturetransformerbase.py'
> --- ffc/quadrature/quadraturetransformerbase.py 2010-07-01 17:41:40 +0000
> +++ ffc/quadrature/quadraturetransformerbase.py 2010-07-07 17:17:07 +0000
> @@ -250,12 +250,16 @@
>         error("This object should be implemented by the child class.")
>
>     # -------------------------------------------------------------------------
> -    # FacetNormal (geometry.py).
> +    # FacetNormal, CellVolume (geometry.py).
>     # -------------------------------------------------------------------------
>     def facet_normal(self, o,  *operands):
>         print "\n\nVisiting FacetNormal: ", repr(o)
>         error("This object should be implemented by the child class.")
>
> +    def cell_volume(self, o,  *operands):
> +        print "\n\nVisiting CellVolume: ", repr(o)
> +        error("This object should be implemented by the child class.")
> +
>     # -------------------------------------------------------------------------
>     # Things that can be handled by the base class.
>     # -------------------------------------------------------------------------
>
>
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 7 Jul 2010 21:14:33 +0200

> DG = FiniteElement("DG", triangle, 0)

Would it be better to call it vol or volume instead of v? Or does it
have to be a one-letter word?

--
Anders

> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 07 Jul 2010 20:22:33 +0100

> have to be a one-letter word?

. . . or call it 'volume'.

Garth

>> On 7 July 2010 18:19,<[email-masked]>  wrote:
>>> committer: Kristian B. Ølgaard<[email-masked]>
>>> branch nick: ffc

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 7 Jul 2010 22:34:39 +0100

> . . . or call it 'volume'.

It can be whatever we want, I just followed what was already there.
Should we then rename 'd', 'n' and 'x' to 'geometric_dimension',
'facet_normal', and 'spatial_coordinate' while we're at it?

Kristian

>>>> To unsubscribe from this branch go to
>>>> https://code.launchpad.net/~ffc-core/ffc/main/+edit-subscription
>>>>
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>     #
>>>> -------------------------------------------------------------------------
>>>>
>>>>
>>>>
>>>

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Jul 2010 08:58:04 +0200

> 'facet_normal', and 'spatial_coordinate' while we're at it?

I think d, n, x are fine, but v does not necessarily look like a
volume to me (it looks like a test function).

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: 08 Jul 2010 08:22:36 +0100

>volume to me (it looks like a test function).

I agree - d, n and x are all commonly used, but v for volume isn't.

Garth

>> >>>> === modified file 'ChangeLog' --- ChangeLog   2010-07-01 17:43:36 
>> >>>> +0000 +++ ChangeLog   2010-07-07 17:17:07 +0000 @@ -1,3 +1,4 @@ + - 
>> >>>> Added support for new geometric quantity CellVolume in UFL.  0.9.3 
>> >>>> [2010-07-01]  - Make global_dimension for Real return an int 
>> >>>> instead of double, bug # 592088  - Add support for facet normal in 
>> >>>> 1D.
>> >>>>
>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume 
>> >>>> (geometry.py).     # 
>> >>>>     def facet_normal(self, o,  *operands):         #print("Visiting 
>> >>>> FacetNormal:") @@ -216,6 +216,17 @@
>> >>>>
>> >>>> +    def cell_volume(self, o,  *operands): +        # Safety 
>> >>>> check. +        ffc_assert(not operands, "Didn't expect any 
>> >>>> operands for FacetNormal: " + repr(operands)) + +        # FIXME: 
>> >>>> KBO: This has to change for higher order elements +        detJ = 
>> >>>> format["det(J)"](self.restriction) +        volume = 
>> >>>> format["absolute value"](detJ) +        self.trans_set.add(detJ) + 
>> >>>> +        return {():create_symbol(volume, GEO)} +     def 
>> >>>> create_argument(self, ufl_argument, derivatives, component, 
>> >>>> local_comp,                   local_offset, ffc_element, 
>> >>>> transformation, multiindices):         "Create code for basis 
>> >>>> functions, and update relevant tables of used basis."
>> >>>>
>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume 
>> >>>> (geometry.py).     # 
>> >>>>     def facet_normal(self, o,  *operands):         #print("Visiting 
>> >>>> FacetNormal:") @@ -291,6 +291,18 @@
>> >>>>
>> >>>> +    def cell_volume(self, o,  *operands): +        # Safety 
>> >>>> check. +        ffc_assert(not operands, "Didn't expect any 
>> >>>> operands for CellVolume: " + repr(operands)) + +        # FIXME: 
>> >>>> KBO: This has to change for higher order elements +        detJ = 
>> >>>> format["det(J)"](self.restriction) +        volume = 
>> >>>> format["absolute value"](detJ) +        self.trans_set.add(detJ) + 
>> >>>> +        return {():volume} + +     def create_argument(self, 
>> >>>> ufl_argument, derivatives, component, local_comp,                   
>> >>>> --- ffc/quadrature/quadraturetransformerbase.py 2010-07-01 17:41:40 
>> >>>> +0000 +++ ffc/quadrature/quadraturetransformerbase.py 2010-07-07 
>> >>>> 17:17:07 +0000 @@ -250,12 +250,16 @@         error("This object 
>> >>>> should be implemented by the child class.")
>> >>>>
>> >>>>     # 
>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume 
>> >>>> (geometry.py).     # 
>> >>>>     def facet_normal(self, o,  *operands):         print 
>> >>>> "\n\nVisiting FacetNormal: ", repr(o)         error("This object 
>> >>>> should be implemented by the child class.")
>> >>>>
>> >>>> +    def cell_volume(self, o,  *operands): +        print 
>> >>>> "\n\nVisiting CellVolume: ", repr(o) +        error("This object 
>> >>>> should be implemented by the child class.") +     # 
>> >>>>     # Things that can be handled by the base class.     #

---

From: Mehdi Nikbakht <[email-masked]>
Date: Thu, 08 Jul 2010 09:33:24 +0200

> I agree - d, n and x are all commonly used, but v for volume isn't.

I have another request related to geometric computation inside UFL. Is
it possible to add a class called SurfaceNormal to geometry.py?

Mehdi

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 08:50:39 +0100

On 08/07/10 08:33, Mehdi Nikbakht wrote:
>>>> On 7 July 2010 20:22, Garth N. Wells<[email-masked]>  wrote:
>>>>>
>>>>>

I don't see how this is possible - it  can't be computed from the cell 
geometry. Shouldn't it be a Constant/Coefficient computed by the user?

Garth

---

From: Mehdi Nikbakht <[email-masked]>
Date: Thu, 08 Jul 2010 10:02:15 +0200

> geometry. Shouldn't it be a Constant/Coefficient computed by the user?

Yes, it can. However, we need to define an extra finite element space
for it. Now we have,

  elem = VectorElement("DG", "triangle", 0)
  n = Coefficient(elem)                       # Surface normal

I would expect to have just one line for this,

  n = SurfaceNormal("triangle")

SurfaceNormal is just a coefficient defined on DG space with zero order.
The same can be done for surface tangent.

Mehdi

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 09:08:57 +0100

> The same can be done for surface tangent.

Making it constant is a big restriction. It should be up to the user.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 8 Jul 2010 09:42:13 +0100

> The same can be done for surface tangent.

This is something completely different, the FacetNormal and CellVolume
are geometric quantities which can be computed from cell information
alone (vertex coordinates) it does not rely on coefficients from any
finite element space. To compute the surface normal and surface
tangent, you would need extra information about the surface which is
not available in tabulate_tensor.

Kristian

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 8 Jul 2010 11:49:14 +0100

> I agree - d, n and x are all commonly used, but v for volume isn't.

I changed 'v' to 'volume' and fixed a few bugs along the way.

On a related note, should we implement CellDiagonal too?
We use cell.diagonal() for MeshSize in SpecialFunctions.h of DOLFIN
and since we often use MeshSize for the 'h' coefficient in DG forms
e.g., PoissonDG we could remove the need for this function evaluation.
The code to compute the CellDiagonal could just be copied from
IntervalCell, TriangleCell and TetrahedronCell of DOLFIN and dumped in
codesnippets.py.

Kristian

>>> >>>> >>>> Added support for new geometric quantity CellVolume in UFL.  0.9.3 >>>>
>>> >>>> [2010-07-01]  - Make global_dimension for Real return an int >>>> instead of
>>> >>>> double, bug # 592088  - Add support for facet normal in >>>> 1D.
>>> >>>>
>>> >>>>     # >>>>
>>> >>>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume >>>>
>>> >>>> (geometry.py).     # >>>>
>>> >>>> +    def cell_volume(self, o,  *operands): +        # Safety >>>>
>>> >>>> check. +        ffc_assert(not operands, "Didn't expect any >>>> operands
>>> >>>> for FacetNormal: " + repr(operands)) + +        # FIXME: >>>> KBO: This has
>>> >>>> to change for higher order elements +        detJ = >>>>
>>> >>>> format["det(J)"](self.restriction) +        volume = >>>> format["absolute
>>> >>>> value"](detJ) +        self.trans_set.add(detJ) + >>>> +        return
>>> >>>> {():create_symbol(volume, GEO)} +     def >>>> create_argument(self,
>>> >>>> ufl_argument, derivatives, component, >>>> local_comp,
>>> >>>> local_offset, ffc_element, >>>> transformation, multiindices):
>>> >>>> "Create code for basis >>>> functions, and update relevant tables of used
>>> >>>> basis."
>>> >>>>
>>> >>>>     # >>>>
>>> >>>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume >>>>
>>> >>>> (geometry.py).     # >>>>
>>> >>>> +    def cell_volume(self, o,  *operands): +        # Safety >>>>
>>> >>>> check. +        ffc_assert(not operands, "Didn't expect any >>>> operands
>>> >>>> for CellVolume: " + repr(operands)) + +        # FIXME: >>>> KBO: This has
>>> >>>> to change for higher order elements +        detJ = >>>>
>>> >>>> format["det(J)"](self.restriction) +        volume = >>>> format["absolute
>>> >>>> value"](detJ) +        self.trans_set.add(detJ) + >>>> +        return
>>> >>>> {():volume} + +     def create_argument(self, >>>> ufl_argument,
>>> >>>> derivatives, component, local_comp,                   >>>> local_offset,
>>> >>>> ffc_element, transformation, multiindices):         >>>> "Create code for
>>> >>>> basis functions, and update relevant tables of >>>> used basis."
>>> >>>>
>>> >>>> === modified file 'ffc/quadrature/quadraturetransformerbase.py' >>>>
>>> >>>> --- ffc/quadrature/quadraturetransformerbase.py 2010-07-01 17:41:40 >>>>
>>> >>>> +0000 +++ ffc/quadrature/quadraturetransformerbase.py 2010-07-07 >>>>
>>> >>>> 17:17:07 +0000 @@ -250,12 +250,16 @@         error("This object >>>> should
>>> >>>> be implemented by the child class.")
>>> >>>>
>>> >>>>     # >>>>
>>> >>>> >>>> -    # FacetNormal (geometry.py). +    # FacetNormal, CellVolume >>>>
>>> >>>> (geometry.py).     # >>>>
>>> >>>> >>>>     def facet_normal(self, o,  *operands):         print >>>>
>>> >>>> "\n\nVisiting FacetNormal: ", repr(o)         error("This object >>>> should
>>> >>>> be implemented by the child class.")
>>> >>>>
>>> >>>> +    def cell_volume(self, o,  *operands): +        print >>>>
>>> >>>> "\n\nVisiting CellVolume: ", repr(o) +        error("This object >>>> should
>>> >>>> be implemented by the child class.") +     # >>>>
>>> >>>> >>>>     # Things that can be handled by the base class.     # >>>>

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 12:00:38 +0100

On 08/07/10 11:49, Kristian Oelgaard wrote:
> On 8 July 2010 08:22, Garth N. Wells<[email-masked]>  wrote:

What if we call it 'Circumradius'?

I recall that the agreed not have have CellSize because of the ambiguity 
in its definition, but circumradius is unambiguous, so I vote to add it.

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Jul 2010 13:59:20 +0200

On Thu, Jul 08, 2010 at 12:00:38PM +0100, Garth N. Wells wrote:
> I recall that the agreed not have have CellSize because of the
> ambiguity in its definition, but circumradius is unambiguous, so I
> vote to add it.

Isn't 'diameter' better? It is shorter and well-defined
(2*circumradius).

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 13:02:43 +0100

> Isn't 'diameter' better? It is shorter and well-defined

'diameter' isn't well defined. We can use it, but we would need to 
define it as being 2*circumradius.

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Jul 2010 14:13:08 +0200

> define it as being 2*circumradius.

Why isn't it well-defined? The diameter of a set T is

  d = sup_{x,y \in T} ||x - y||

and that should be the same as 2*circumradius unless I'm mistaken.

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 14:36:09 +0100

> and that should be the same as 2*circumradius unless I'm mistaken.

What's the set T? At a glance, the above is true if T is not the 
polygon, but the set of all points inside the circle inside which the 
polygonal can be inscribed.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Thu, 8 Jul 2010 15:28:42 +0100

On 8 July 2010 14:36, Garth N. Wells <[email-masked]> wrote:
> What's the set T? At a glance, the above is true if T is not the polygon,
> but the set of all points inside the circle inside which the polygonal can
> be inscribed.

I added cell.circumradius, if we want cell.diameter too we can easily add it.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Jul 2010 20:42:09 +0200

On Thu, Jul 08, 2010 at 02:36:09PM +0100, Garth N. Wells wrote:
> polygon, but the set of all points inside the circle inside which
> the polygonal can be inscribed.

Yes, don't know what I was thinking. Let's try again... For a
polygon, the diameter should be the same as the longest edge which
should be very easy to compute so why not use that?

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 08 Jul 2010 20:03:35 +0100

> should be very easy to compute so why not use that?

Longest edge and circumradius (which is also called the radius of a 
polygon, just to make things confusing) are not the same.

Garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 8 Jul 2010 21:25:47 +0200

> polygon, just to make things confusing) are not the same.

I see that but I was suggesting having 'diameter' in addition to
'circumradius' since 'diameter' is shorter and simpler and easier to
compute. It also seems to me like the difference between the two is
small, in particular for thin triangles where the two should be the
same (diameter = 2*circumradius).

--
Anders