project(): boundary oszillations, subdomain restrictions

Are you sure x[0] is not 0.0 in that domain?

If so this statement should be True which it is not:

  np.all(np.abs(submesh.coordinates()[:,0])>1e-10)

Johan

On 11/13/2012 07:41 PM, Nico Schlömer wrote:
> New question #214113 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/214113
>
> I discovered two ways in which Dolfin's project() behaves that I wouldn't have expected. Both are highlighted by the following setup.
>
> Suppose you have a domain Omega=[0,3]x[-1,1] that divides into a couple of subdomains. We would like to create the function
>
> ( 1/x[0] for the subdomains 1,2,3,
> f(x) = (
> ( 0 everywhere else.
>
> (The subdomains 1,2,3 don't contain points with x[0]==0.)
> The definition of f() can happen via an expression. When plotting the function, it has to be projected onto a meaningful space, e.g., FunctionSpace(mesh, 'CG', 1). When using CG1, oszillations around the boundary of the subdomains 1,2,3 occur.
> What is the reason for those?
> (Using DG0 avoids avoids such oszillations.)
>
> Now we would like to restrict f to subdomain 4 (which does contain points x[0]==0). We would expect it to be 0 throughout.
>
> Overever, when projecting onto
>
> submesh = SubMesh(mesh, subdomains, 4)
> Vsub = FunctionSpace(submesh, 'CG', 1)
> tmp = project(oneOverX, Vsub)
>
> the expression 1/x[0] is evaluated with x[0]==0.
> Why is that?
>
> The complete failing code is
>
> ====================== *snip* ======================
> from dolfin import *
>
> # -----------------------------------------------------------------------------
> def _main():
>
> # Read the mesh and its subdomains.
> mesh = Mesh('coils2d.xml')
> subdomains = MeshFunction('uint', mesh, 'coils2d_physical_region.xml')
> coil_indices = [1,2,3]
>
> class OneOverX(Expression):
> def eval_cell(self, values, x, cell):
> k = subdomains.array()[cell.index]
> try:
> index = coil_indices.index(k)
> assert abs(x[0]) > 1.0e-10, 'Division by 0.'
> values[0] = 1.0/ x[0]
> except ValueError:
> values[0] = 0.0
> oneOverX = OneOverX()
>
> # Plot it.
> V = FunctionSpace(mesh, 'CG', 1)
> File('oneOverX.pvd') << project(oneOverX, V)
>
> # Project onto subdomain.
> submesh = SubMesh(mesh, subdomains, 4)
> Vsub = FunctionSpace(submesh, 'CG', 1)
> tmp = project(oneOverX, Vsub)
> File('oneOverX-submesh.pvd') << tmp
>
> return
> # -----------------------------------------------------------------------------
> if __name__ == '__main__':
> _main()
> # -----------------------------------------------------------------------------
> ====================== *snap* ======================
>
> and the domain data can be retrieved from http://win.ua.ac.be/~nschloe/other/fenics/.
>

>> (which does contain points x[0]==0)
> Are you sure x[0] is not 0.0 in that domain?

The submesh onto which we project DOES contain x[0]==0, but that should not matter.

On 11/13/2012 08:35 PM, Nico Schlömer wrote:
> Question #214113 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214113
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
>>> (which does contain points x[0]==0)
>> Are you sure x[0] is not 0.0 in that domain?
>
> The submesh onto which we project DOES contain x[0]==0, but that should
> not matter.

Why?

Johan

Well, those subdomains 1,2,3 don't contain x[0]==0, so f() is meaningfully defined on the whole domain. Taking a peek at oneOverX.pvd (when plotted with DG0) tells you what things look like.
I would like to project f() onto a subdomain which happens to contain x[0]==0, but it should be f(x)==0 there anyways.

On 11/13/2012 08:45 PM, Nico Schlömer wrote:
> Question #214113 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214113
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
> Well, those subdomains 1,2,3 don't contain x[0]==0, so f() is meaningfully defined on the whole domain. Taking a peek at oneOverX.pvd (when plotted with DG0) tells you what things look like.
> I would like to project f() onto a subdomain which happens to contain x[0]==0, but it should be f(x)==0 there anyways.

I do not understand what you do within the eval method. k and index are
never used. Why use a try and except clause?

Johan

Johan: coil_indices.index(k) throws a value error if k is not in coil_indices. I guess the intent is that if k is not 1, 2, 3, then a value error will be thrown and caught and so values = 0.0

Nico: Could the problem be that you are assuming that the cell numbering in the mesh and the submesh are the same? Printing cell.index, k, index and x in the eval_cell, leads it to fail at cell 85 (relative to the submesh) where the corresponding subdomains marker (relative to the mesh) is 1, and for which x[0] is indeed 0.

> coil_indices.index(k) throws a value error [...]

Exactly that's what it does.

@Marie
That's interesting and likely the cause.
I suppose that the low-level access to subdomains.array() bites me here. Is there a more consistent way of finding out what subdomain I'm in?

On 11/13/2012 09:45 PM, Nico Schlömer wrote:
> Question #214113 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214113
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
>> coil_indices.index(k) throws a value error [...]
>
> Exactly that's what it does.

Ok but:

  if subdomains.array()[cell.index] in coil_indices:
    values[0] = 0.0
  else:
    values[0] = 1.0/ x[0]

is more explicit :)

> @Marie
> That's interesting and likely the cause. I suppose that the low-level
> access to subdomains.array() bites me here. Is there a more
> consistent way of finding out what subdomain I'm in?

I guess this works as intended for the full mesh, right?

For the submesh it really does not make sense to use the cell.index to
find the domain, as all cells have marker 4.

If you want to keep you logic you need to create a new "subdomains" for
the submesh, where cells with vertices x[0]=0 are marked with a new marker.

Johan

Johan

> if subdomains.array()[cell.index] in coil_indices:
> values[0] = 0.0
> else:
> values[0] = 1.0/ x[0]
>
> is more explicit :)

And how! This is way better, thanks.

> I guess this works as intended for the full mesh, right?

Right.

> For the submesh it really does not make sense to use the cell.index to
> find the domain, as all cells have marker 4.

Makes sense. It's clear now what the problem is, and I'll try to find an appropriate solution (maybe involving the definition of an Expression just on subdomain 4 or something like that).

The only open question right now would be the appearance of oszillations for the CG1-projection. I'm guessing the cause for this will be numerical since I'm trying to match a discontinuous function (the Expression) with continuous elements. You'd expect something to go wrong, but why does it oscillate?

On 11/14/2012 09:11 AM, Nico Schlömer wrote:
> Question #214113 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/214113
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
>> if subdomains.array()[cell.index] in coil_indices:
>> values[0] = 0.0
>> else:
>> values[0] = 1.0/ x[0]
>>
>> is more explicit :)
>
> And how! This is way better, thanks.
>
>> I guess this works as intended for the full mesh, right?
>
> Right.
>
>> For the submesh it really does not make sense to use the cell.index to
>> find the domain, as all cells have marker 4.
>
> Makes sense. It's clear now what the problem is, and I'll try to find an
> appropriate solution (maybe involving the definition of an Expression
> just on subdomain 4 or something like that).
>
> The only open question right now would be the appearance of oszillations
> for the CG1-projection. I'm guessing the cause for this will be
> numerical since I'm trying to match a discontinuous function (the
> Expression) with continuous elements. You'd expect something to go
> wrong, but why does it oscillate?

Not my main domain of expertise but have in mind that the expression
gets interpolated during assemble to a FiniteElement space, which means
it get evaluated at the dofs for that element. The element is by default
CG1 for Expressions.

You might want to try a Discontinuous element instead. DG1 if you want
nodal values, or DG0 if you are good with barry center values.

Johan

That does indeed work; see above.
I was just wondering *why* we have those funny effects there. (This is purely for my own enlightenment. -- I'll see if I can find something in liturature here.)

This question was expired because it remained in the 'Open' state without activity for the last 15 days.