Allow integrals over entire domain mixed with integrals over subdomains

Registered by Martin Sandve Alnæs

(This is a crosscutting ufl-ffc-ufc-dolfin feature)

Today you can only integrate over the entire domain if you explicitly make sure to pass a meshfunction which has the same value everywhere. Related to ongoing work on multidomain support, I want to introduce the concept of integrals over the full domain.

This will involve changing the default meaning of f*dx, f*ds, etc. to mean integration over everything, and adding warnings to mixed dx and dx(1) notation such as "f*dx + f*dx(1)" if we want to play it safe. Most code that use such mixed notation (with and without domain ids) is probably a bug anyway. New notation will be something like dx(D) where D is a Domain (or a Mesh in dolfin), or dx("everywhere") or whatever. I can handle this smoothly on the UFL side but there are some issues that are more crosscutting.

In particular we need to handle this from the assembler point of view. The integration subdomains today are always disjoint, which means there is only one ufc integral object to call tabulate_tensor on for each cell in the mesh. With the upcoming concept of Region in UFL as a set of subdomain ids, we can preprocess integrals such that
  f*dx((1,2)) + g*dx((2,3))
becomes
  f*dx(1) + (f + g)*dx(2) + g*dx(3)
prior to form compilation (e.g. in ufl preprocess), and we therefore preserve the "disjointness" of the subdomains.

Now considering the integral
  f*dx("everywhere") + g*dx(1)
we cannot do the same transformation in UFL, because the set of all subdomain ids are not available before assembly time.

One solution is to add a special integral object such as "ufc::form::create_cell_integral_everywhere()" and let the assembler tabulate the local tensor for both integrals f*dx("everywhere") and g*dx(1) when we are on subdomain 1. This gives a loss of assembly efficiency.

To avoid this efficiency loss (which I would prefer), we can still preprocess the form on the ufl side if we add an integral object to ufc with different semantics. This integral is the integral over all subdomains where no other integral in the form contributes, e.g. "ufc::form::create_cell_integral_everywhere_else()". The assembler will then need to, for each cell:
  - Get subdomain id
  - Get integral for this subdomain id
  - NEW: If no integral, use default integral (if any)
  - Call tabulate tensor etc.
the cost of the additional "if (!integral) integral = default_integral;" should be neglible.

The integral over the entire domain then becomes a special case of a Region, and we need to do a similar tranformation in the preprocess stage (or maybe in FFC). This transformation requires a bit more care than the region->subdomain splitting shown above, because it cannot be applied twice. FFC currently calls preprocess twice, which is a bad idea for efficiency but not mathematically incorrect as of today. While we should avoid that in FFC, it is still a good property that preprocess _can_ be applied twice without changing the meaning of a form. Therefore I suggest UFL provides the algorithm and FFC applies it at an appropriate time.

The algorithm then looks simply like this (separately for each integral type dx,ds,etc.):
  find everywhere-integral
  for each integral in form:
      add everywhere-integral to integrand
for example,
  f*dx((1,2)) + g*dx((2,3)) + h*dx("everywhere")
becomes
  (h+f)*dx(1) + (h + f + g)*dx(2) + (h + g)*dx(3) + h*dx("everywhere else")
(Ok, so if we represent "everywhere" and "everywhere else" distinctly, we can apply this twice with no change the second time, but if you add e.g. q*dx(4) to the last form and then process again, it becomes muddy. So lets just don't do that.)

What do you think?

The first thing I want to do is to quickly and brutally deprecate the dx as dx(0) interpretation. This is already impossible if you use the dx[boundaries] feature in python, and code that depends on this behaviour is likely a bug.

Blueprint information

Status:
Complete
Approver:
None
Priority:
Essential
Drafter:
None
Direction:
Needs approval
Assignee:
Martin Sandve Alnæs
Definition:
Approved
Series goal:
None
Implementation:
Implemented
Milestone target:
None
Started by
Martin Sandve Alnæs
Completed by
Martin Sandve Alnæs

Related branches

Sprints

Whiteboard

This is the new behaviour (same for dx, ds, dS, just using dx here):

f*dx() is now the integral over the entire domain (may be removed)

f*dx() + g*dx(1) is allowed

f*dx() + g*dx((1,2)) + h*dx((2,3)) is allowed as well

f*dx is now ALSO the integral over the entire domain, changing previous behaviour, possibly breaking some old code

f*dx + g*dx(1) is currently allowed, and means the same as f*dx() + g*dx(1)

Old user code that specifies only *dx terms and no *dx(i) terms, and assembles without subdomains or with a meshfunction that is zero everywhere, will work as before.

Old user code that specifies no *dx terms and only *dx(i) terms, and assembles with subdomains will work as before.

Old user code that specifies *dx(i) terms but assemble with no subdomains may break.

Old user code that mixes *dx and *dx(i) terms and assembles with subdomains, or otherwise depends on the identity dx == dx(0) which is no longer valid, will break.

(?)

Work Items

Work items:
[martinal] UFL; Implement f*dx() + g*dx(1) syntax: DONE
[martinal] UFC; Introduce dummy ufc::form::has_*_integral() functions in UFC to check for both default and subdomain integrals easily: DONE
[martinal] UFC; Introduce dummy ufc::form::create_default_*_integral() functions: DONE
[martinal] DOLFIN; Use ufc::form::has_*_integral() functions: DONE
[martinal] FFC; Generate dummy ufc::form::has_*_integral() functions: DONE
[martinal] FFC; Generate actual ufc::form::has_*_integral() functions in FFC: DONE
[martinal] FFC; Generate dummy ufc::form::create_default_*_integral() functions: DONE
[martinal] UFL; Handle foo*dx as integral over everywhere: DONE
[martinal] UFL; Enable everywhere integrals and transform to otherwise integrals: DONE
[martinal] FFC; Generate actual ufc::form::create_default_*_integral() functions in FFC: DONE
[martinal] DOLFIN; Use ufc::form::create_default_*_integral() functions: DONE
[martinal] FFC; DOLFIN; Update error control code to use dx() syntax for generated forms: DONE
[martinal] DOLFIN; Update demos, tests, norms, etc. to use dx() syntax for everywhere integrals: DONE
[martinal] UFL; Canonicalize integrals from subdomains/regions/everywhere to subdomains/otherwise in preprocess: DONE
[martinal] DOLFIN; Make error control unit test pass: DONE

This blueprint contains Public information 
Everyone can see this information.