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

Bug #862383 reported by Martin Sandve Alnæs
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
DOLFIN
Fix Released
High
Anders Logg

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)

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :

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.
*** -------------------------------------------------------------------------

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :
Revision history for this message
Kristian B. Ølgaard (k.b.oelgaard) wrote :

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?

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :

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.

Changed in dolfin:
milestone: none → 1.0-rc1
Changed in dolfin:
importance: Undecided → High
Anders Logg (logg)
Changed in dolfin:
assignee: nobody → Anders Logg (logg)
Revision history for this message
Anders Logg (logg) wrote :
Changed in dolfin:
status: New → In Progress
Anders Logg (logg)
Changed in dolfin:
status: In Progress → Fix Committed
Revision history for this message
Kristian B. Ølgaard (k.b.oelgaard) wrote : Re: [Bug 862383] Re: FacetArea evaluation fails if the form includes both a cell integral and an interior facet integral

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
>

Revision history for this message
Anders Logg (logg) wrote :

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

Anders Logg (logg)
Changed in dolfin:
status: Fix Committed → Fix Released
To post a comment you must log in.
This report contains Public information  
Everyone can see this information.

Other bug subscribers

Remote bug watches

Bug watches keep track of this bug in other bug trackers.