Expression conditional on subdomain

On 11/09/2012 05:41 PM, Nico Schlömer wrote:
> New question #213769 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/213769
>
> I'd like to define an expression that does (sorry for the ASCII art
>
> / 1.0/x[0] in subdomain 1
> f(x) = (
> \ 1.0 in subdomain 2
>
> (Subdomain 1 doesn't contain any points with x[0].)

So what does subdomain 1 contain then?

> I'm quite used to defining custom Expressions like
>
> class MyExpr(Expression):
> def eval_cell(self, values, x, cell):
> k = subdomains.array()[cell.index]
> if k==1:
> values[0] = 1.0 / x[0]
> else:
> values[0] = 0.0
>
> but x[0] isn't what I want here.

What do you want then? Looks good to me.

> (I'm not even sure what it does in eval_cell(). -- The barycenter?)

This depends on what Element the Expression has. For all I know it will
be the coordinates of the dofs in a local element. So if the element is
DG0 it will be the barry center.

> There are also UFL conditionals, but I couldn't find out how to make them depend on subdomains.array()[cell.index].
> Any hints?
>

>> (Subdomain 1 doesn't contain any points with x[0].)
> So what does subdomain 1 contain then?

That should have been x[0]==0. You picked it up already.

> What do you want then? Looks good to me.

I'm getting weird cases where x[0] is 0 where it shouldn't but this only occurs when I project() onto a subdomain. I'll try and construct a minimally failing example.

This issue is related to one or two other ones that came up recently.
The problem with defining subdomains in the way

    class MyExpr(Expression):
        def eval_cell(self, values, x, cell):
            k = subdomains.array()[cell.index]
            if k==1:
                values[0] = 1.0 / x[0]
            else:
                values[0] = 0.0

is that cell.index *must* be the global cell index. If MyExpr() is evaluated on the entire mesh, this is working fine. However, if MyExpr() is evaluated only on a subdomain, e.g.,

    integral = assemble(MyExpr()*dx, mesh = SubMesh(mesh, subdomains, 7))

the returned values will be bogus.

This makes me thing that accessing subdomains.array() isn't a very good way of handling things. Is there anything else that you'd recommend?

On 12/04/2012 03:21 PM, Nico Schlömer wrote:
> Question #213769 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/213769
>
> Status: Solved => Open
>
> Nico Schlömer is still having a problem:
> This issue is related to one or two other ones that came up recently.
> The problem with defining subdomains in the way
>
> class MyExpr(Expression):
> def eval_cell(self, values, x, cell):
> k = subdomains.array()[cell.index]
> if k==1:
> values[0] = 1.0 / x[0]
> else:
> values[0] = 0.0
>
> is that cell.index *must* be the global cell index. If MyExpr() is
> evaluated on the entire mesh, this is working fine. However, if MyExpr()
> is evaluated only on a subdomain, e.g.,
>
> integral = assemble(MyExpr()*dx, mesh = SubMesh(mesh, subdomains,
> 7))
>
> the returned values will be bogus.
>
> This makes me thing that accessing subdomains.array() isn't a very good
> way of handling things. Is there anything else that you'd recommend?

As it stands now your script does not make sense. You assemble an
Expression which is only non-zero for global cells of index 1. But then
you assemble over a submesh global index 7.

What are you trying to do?

Johan

Right, I better provide a complete example; get it at http://win.ua.ac.be/~nschloe/other/chi/.
You'll find a domain with 6 subdomains Omega_i and the characteristic function (expression) Chi for the union of the Omega_i. The code then integrates Chi() over each Omega_i. Wrong values come out.
The correct values are retrieved when projection Chi() into a fully qualified function space before integrating (e.g., FunctionSpace(mesh, 'DG', 0)).
Anyways, this is supposed to highlight how one can easily make a mistake when defining Expressions the way they are defined here.

I cannot get your example to run because of changes in the recent trunk.
uint->size_t. If you use subdomain integrals, dx(coilindex), instead of
SubMesh you will preserve the global numbering of the cells.

Johan

On 12/04/2012 11:06 PM, Nico Schlömer wrote:
> Question #213769 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/213769
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
> Right, I better provide a complete example; get it at http://win.ua.ac.be/~nschloe/other/chi/.
> You'll find a domain with 6 subdomains Omega_i and the characteristic function (expression) Chi for the union of the Omega_i. The code then integrates Chi() over each Omega_i. Wrong values come out.
> The correct values are retrieved when projection Chi() into a fully qualified function space before integrating (e.g., FunctionSpace(mesh, 'DG', 0)).
> Anyways, this is supposed to highlight how one can easily make a mistake when defining Expressions the way they are defined here.
>

I recommend splitting your form into multiple integrals and using different
coefficients on each subdomain:

  a = c1*u*v*dx(1) + c2*u*v*dx(2)
instead of
   a = c*u*v*dx

then you move the complexity out of the coefficient c.

Note that in latest trunk (of ufl) you can add measures like this
  dx12 = dx(1) + dx(2)
  a = c1*u*v*dx(1) + c2*u*v*dx(2) + alpha*u*v*dx12
so you don't have to repeat terms that are the same on two subdomains.

On 4 December 2012 23:06, Nico Schlömer <
<email address hidden>> wrote:

> Question #213769 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/213769
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
> Right, I better provide a complete example; get it at
> http://win.ua.ac.be/~nschloe/other/chi/.
> You'll find a domain with 6 subdomains Omega_i and the characteristic
> function (expression) Chi for the union of the Omega_i. The code then
> integrates Chi() over each Omega_i. Wrong values come out.
> The correct values are retrieved when projection Chi() into a fully
> qualified function space before integrating (e.g., FunctionSpace(mesh,
> 'DG', 0)).
> Anyways, this is supposed to highlight how one can easily make a mistake
> when defining Expressions the way they are defined here.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I suppose that'd work although I can't say I like the idea of moving complexity into the already complicated bilinear form.

I have one question remaining:
I would like to keep the number of subdomains variable, so ideally I would like to do something like

# define coefficients for the subdomains
coeffs = {1: Expression('x[0]'),
                  2: Expression('2*pi'),
                  3: Expression('x[0]*x[1]'),
                  ...}
...
a = <???>
for i in coeffs:
    a += inner(coeffs[i] * grad(u), grad(v)) * dx(i)

How do I sensibly initialize a? Is there such a thing as the "empty expression"?

Try:
  import ufl.constantvalue
  a = ufl.constantvalue.Zero()*dx

Johan

import ufl.constantvalue
a = ufl.constantvalue.Zero() * dx(0)

works. Thanks!

zero()*dx(0) would be the official UFL syntax.

On 01/14/2013 05:31 PM, Martin Sandve Alnæs wrote:
> Question #213769 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/213769
>
> Martin Sandve Alnæs posted a new comment:
> zero()*dx(0) would be the official UFL syntax.

Ahh, I missed zero...

Johan