FacetArea evaluation fails if the form includes both a cell integral and an interior facet integral

Traceback (most recent call last):
  File "facetareabug.py", line 17, in <module>
    A2 = assemble(a2)
  File "/opt/fenics/dorsal-unstable/lib/python2.7/site-packages/dolfin/fem/assembling.py", line 180, in assemble
    add_values)
  File "/opt/fenics/dorsal-unstable/lib/python2.7/site-packages/dolfin/cpp.py", line 19792, in assemble
    return _cpp.assemble(*args)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to evaluate FacetArea expression.
*** Reason: Facet area is only defined on mesh boundaries.
*** Where: This error was encountered inside SpecialFunctions.cpp.
*** -------------------------------------------------------------------------

Is this is also a performance issue? Consider the form:

M = c*d*e*f*g*h*dx + c*ds + c('+')*dS

where c,...,h are Expressions, then the eval method of d,...,h will be called for each facet in the mesh although it's only necessary to evaluate them for each cell?

Yes, that would be negative for performance. But I guess that is actually the same for all kinds of subdomains? For example f*dx(1)+g*dx(2) would cause evaluation of both f and g for both domains 1 and 2, because there is no propagation of dependencies through UFC. I would settle for a hack making this case work correctly (a warning would be fine), and the performance issue should be addressed in the future through some design involving UFC.

Hack added in 1.0.x-logg.

See also https://blueprints.launchpad.net/ufc/+spec/missing-coefficient-domain-dependency

Would it be an idea to also add 'FacetArea' as a GeometricQuantity to
Cell in UFL?
Then h_F in the bug report can be defined as:

h_F = cell.facet_area

i.e., the area of the facet currently being integrated.

One could perhaps also consider defining SurfaceArea, the sum of all
FacetAreas of a Cell.

Note that:

M = h_F*dx

should be illegal, but this can be checked in UFL at compile time.

Kristian

On 17 November 2011 01:24, Anders Logg <email address hidden> wrote:
> ** Changed in: dolfin
>       Status: In Progress => Fix Committed
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/862383
>
> Title:
>  FacetArea evaluation fails if the form includes both a cell integral
>  and an interior facet integral
>
> Status in DOLFIN:
>  Fix Committed
>
> Bug description:
>  Reproduce with this script:
>
>  from dolfin import *
>
>  mesh = UnitSquare(3,3)
>  V = VectorFunctionSpace(mesh, "DG", 1)
>  v = TestFunction(V)
>  u = TrialFunction(V)
>
>  h_F = FacetArea(mesh)('+')
>
>  # Ok:
>  a1 = h_F*dot(jump(u),jump(v))*dS
>  A1 = assemble(a1)
>
>  # Fails:
>  a2 = a1 + inner(grad(u), grad(v))*dx
>  A2 = assemble(a2)
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/862383/+subscriptions
>

Yes, that would be a better solution. Many other SpecialFunctions have
gone that way.

--
Anders

On Thu, Nov 17, 2011 at 08:01:53AM -0000, Kristian B. Ølgaard wrote:
> Would it be an idea to also add 'FacetArea' as a GeometricQuantity to
> Cell in UFL?
> Then h_F in the bug report can be defined as:
>
> h_F = cell.facet_area
>
> i.e., the area of the facet currently being integrated.
>
> One could perhaps also consider defining SurfaceArea, the sum of all
> FacetAreas of a Cell.
>
> Note that:
>
> M = h_F*dx
>
> should be illegal, but this can be checked in UFL at compile time.
>
> Kristian
>
>
> On 17 November 2011 01:24, Anders Logg <email address hidden> wrote:
> > ** Changed in: dolfin
> >       Status: In Progress => Fix Committed
> >
> >
> > Title:
> >  FacetArea evaluation fails if the form includes both a cell integral
> >  and an interior facet integral
> >
> > Status in DOLFIN:
> >  Fix Committed
> >
> > Bug description:
> >  Reproduce with this script:
> >
> >  from dolfin import *
> >
> >  mesh = UnitSquare(3,3)
> >  V = VectorFunctionSpace(mesh, "DG", 1)
> >  v = TestFunction(V)
> >  u = TrialFunction(V)
> >
> >  h_F = FacetArea(mesh)('+')
> >
> >  # Ok:
> >  a1 = h_F*dot(jump(u),jump(v))*dS
> >  A1 = assemble(a1)
> >
> >  # Fails:
> >  a2 = a1 + inner(grad(u), grad(v))*dx
> >  A2 = assemble(a2)
> >
> > To manage notifications about this bug go to:
> > https://bugs.launchpad.net/dolfin/+bug/862383/+subscriptions
> >
>