Computing the mesh volume using assemble(Constant(1)*dx, mesh) yields zero for certain meshes

You can also try:
vol = 0.
for cell in cells(b):
    volume += cell.volume()

Johan
On Oct 29, 2012 4:15 PM, "Maximilian Albert" <
<email address hidden>> wrote:

> New question #212697 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/212697
>
> I'm using dolfin version 1.0.0-2+bzr7035~ppa1~precise1 (provided via the
> development build PPA).
>
> In https://answers.launchpad.net/dolfin/+question/168647, Johan Hake
> suggests the following method to compute the volume of a mesh:
>
> volume = assemble(Constant(1), mesh=mesh)
>
> I have used this a lot with the stable release of FEniCS and it has worked
> fine for me so far. However, after a recent upgrade to the FEniCS
> development version, it has started to fail for certain meshes, with the
> volume being reported as zero. I have only seen this for meshes created by
> netgen so far (although I haven't investigated much, so this may be a red
> herring). Here is an example to reproduce this:
>
> ==>
> from dolfin import *
> brick = mesh.netgen.Brick(0,0,0,1,1,1)
> brick.generate_mesh("brick.xml", "verycoarse")
> b = Mesh("brick.xml")
> print assemble(Constant(1)*dx, mesh=b)
> <==
>
> Result: 0.0
>
> I don't know if the orientation of the cells is taken into account, and if
> this might have an impact on the result (e.g. some of the cells cancel each
> other because they have "negative volume"). However, "b.ordered()" yields
> "True" so it seems that this can't be the cause of the problem. Could
> someone please clarify whether the above is expected behaviour, and if so
> how I can ensure that I always get the correct volume for a mesh? Or is it
> a bug that was introduced in the development version? Many thanks!
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Many thanks for the quick reply, Johan! It looks like your new suggestion works. :) Do you have any idea why the other method fails? Is this expected, or is it a bug/regression?

On 10/30/2012 05:21 AM, Maximilian Albert wrote:
> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Maximilian Albert posted a new comment:
> Many thanks for the quick reply, Johan! It looks like your new
> suggestion works. :) Do you have any idea why the other method fails? Is
> this expected, or is it a bug/regression?

Well, it looks like some well hidden feature coming to surface here :)

The mesh generated from netgen is outputted as a diffpack mesh(!?). The
diffpack converter has gained some love recently, where all markers are
now stored within the mesh as MeshValueCollections. Have a look at the
end of the brick.xml file and you see the markers.

So now all cells are marked with marker 1. This means you need to define
your integral over this domain. dx is short for dx(0), which is a volume
integral over the zeroeth domain. So just change dx to dx(1) in your
form and you should be back in business.

Johan

Thanks again for the very helpful and detailed response, Johan! Much appreciated indeed. Just to make sure I understand correctly, are the "domains" you refer to a way of distinguishing different regions in the mesh (i.e. a logical grouping, e.g. because they belong to the same "material")? If so, does this mean that it will become easier to import meshes produced by Netgen which contain different materials/regions into Finmag? I seem to remember that this used to be difficult because the information wasn't retained.

Also, regarding your comment "dx is short for dx(0)": would it be possible to make "dx" apply to the entire mesh by default (i.e. to *all* domains, not just the 0-th domain)? This is what I would naturally expect as a naive user who doesn't know about the internal workings of FEniCS, nor even the existence of MeshValueCollections (until you pointed me to them).

Many thanks again!

On 10/30/2012 02:56 PM, Maximilian Albert wrote:
> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Maximilian Albert posted a new comment:
> Thanks again for the very helpful and detailed response, Johan! Much
> appreciated indeed. Just to make sure I understand correctly, are the
> "domains" you refer to a way of distinguishing different regions in the
> mesh (i.e. a logical grouping, e.g. because they belong to the same
> "material")?

Yes

> If so, does this mean that it will become easier to import
> meshes produced by Netgen which contain different materials/regions into
> Finmag? I seem to remember that this used to be difficult because the
> information wasn't retained.

Not sure, but it looks like it. But as I said it looks like the built in
python netgen mesh builder uses diffpack representation as an
intermediate. If the markers are converted correctly between netgen and
diffpack they should also be there in the stored mesh.

> Also, regarding your comment "dx is short for dx(0)": would it be
> possible to make "dx" apply to the entire mesh by default (i.e. to *all*
> domains, not just the 0-th domain)? This is what I would naturally
> expect as a naive user who doesn't know about the internal workings of
> FEniCS, nor even the existence of MeshValueCollections (until you
> pointed me to them).

It should be possible and we have discussed it previously. It keeps on
biting both frequent and new users of DOLFIN. However we would need to
change stuff in the ufc interface to get that going (if I remember
correctly), and we are pretty conservative to such changes.

Johan

> Many thanks again!
>

On Wed, Oct 31, 2012 at 06:31:09AM -0000, Johan Hake wrote:
> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Johan Hake proposed the following answer:
> On 10/30/2012 02:56 PM, Maximilian Albert wrote:
> > Question #212697 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/212697
> >
> > Maximilian Albert posted a new comment:
> > Thanks again for the very helpful and detailed response, Johan! Much
> > appreciated indeed. Just to make sure I understand correctly, are the
> > "domains" you refer to a way of distinguishing different regions in the
> > mesh (i.e. a logical grouping, e.g. because they belong to the same
> > "material")?
>
> Yes
>
> > If so, does this mean that it will become easier to import
> > meshes produced by Netgen which contain different materials/regions into
> > Finmag? I seem to remember that this used to be difficult because the
> > information wasn't retained.
>
> Not sure, but it looks like it. But as I said it looks like the built in
> python netgen mesh builder uses diffpack representation as an
> intermediate. If the markers are converted correctly between netgen and
> diffpack they should also be there in the stored mesh.

The Netgen mesh generator in DOLFIN will soon be removed and replaced
by the CGAL-based mesh generator in the CSG branch.

> > Also, regarding your comment "dx is short for dx(0)": would it be
> > possible to make "dx" apply to the entire mesh by default (i.e. to *all*
> > domains, not just the 0-th domain)? This is what I would naturally
> > expect as a naive user who doesn't know about the internal workings of
> > FEniCS, nor even the existence of MeshValueCollections (until you
> > pointed me to them).
>
> It should be possible and we have discussed it previously. It keeps on
> biting both frequent and new users of DOLFIN. However we would need to
> change stuff in the ufc interface to get that going (if I remember
> correctly), and we are pretty conservative to such changes.

I think we concluded that it wasn't worth the trouble. The current
logic is:

0. dx = dx(0)

1. If a user doesn't mark the mesh, dx = dx(0) will integrate over the
entire mesh.

2. If a user marks the mesh, those markers will be used. In particular
dx = dx(0) will integrate only over domain 0.

4. We don't have a shortcut for "entire mesh".

We could discuss if it's worth to add such a shortcut. It would
require introducing a label in UFL for "all domains", something like
dx("all") and letting dx = dx("all"). It would also require UFC to
have an additional function, something like
create_cell_integral_default or create_cell_integral_all.

--
Anders

We should definitely handle this case when designing the improved
multidomain support next year.

Martin
Den 31. okt. 2012 08:35 skrev "Anders Logg" <
<email address hidden>> følgende:

> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Anders Logg proposed the following answer:
> On Wed, Oct 31, 2012 at 06:31:09AM -0000, Johan Hake wrote:
> > Question #212697 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/212697
> >
> > Johan Hake proposed the following answer:
> > On 10/30/2012 02:56 PM, Maximilian Albert wrote:
> > > Question #212697 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/212697
> > >
> > > Maximilian Albert posted a new comment:
> > > Thanks again for the very helpful and detailed response, Johan! Much
> > > appreciated indeed. Just to make sure I understand correctly, are the
> > > "domains" you refer to a way of distinguishing different regions in the
> > > mesh (i.e. a logical grouping, e.g. because they belong to the same
> > > "material")?
> >
> > Yes
> >
> > > If so, does this mean that it will become easier to import
> > > meshes produced by Netgen which contain different materials/regions
> into
> > > Finmag? I seem to remember that this used to be difficult because the
> > > information wasn't retained.
> >
> > Not sure, but it looks like it. But as I said it looks like the built in
> > python netgen mesh builder uses diffpack representation as an
> > intermediate. If the markers are converted correctly between netgen and
> > diffpack they should also be there in the stored mesh.
>
> The Netgen mesh generator in DOLFIN will soon be removed and replaced
> by the CGAL-based mesh generator in the CSG branch.
>
> > > Also, regarding your comment "dx is short for dx(0)": would it be
> > > possible to make "dx" apply to the entire mesh by default (i.e. to
> *all*
> > > domains, not just the 0-th domain)? This is what I would naturally
> > > expect as a naive user who doesn't know about the internal workings of
> > > FEniCS, nor even the existence of MeshValueCollections (until you
> > > pointed me to them).
> >
> > It should be possible and we have discussed it previously. It keeps on
> > biting both frequent and new users of DOLFIN. However we would need to
> > change stuff in the ufc interface to get that going (if I remember
> > correctly), and we are pretty conservative to such changes.
>
> I think we concluded that it wasn't worth the trouble. The current
> logic is:
>
> 0. dx = dx(0)
>
> 1. If a user doesn't mark the mesh, dx = dx(0) will integrate over the
> entire mesh.
>
> 2. If a user marks the mesh, those markers will be used. In particular
> dx = dx(0) will integrate only over domain 0.
>
> 4. We don't have a shortcut for "entire mesh".
>
> We could discuss if it's worth to add such a shortcut. It would
> require introducing a label in UFL for "all domains", something like
> dx("all") and letting dx = dx("all"). It would also require UFC to
> have an additional function, something like
> create_cell_integral_default or create_cell_integral_all.
>
> --
> Anders
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Thanks Johan Hake, that solved my question.

Sorry to revive this old thread, but it is possible that I am being bitten by the "dx = dx(0)" issue again. I'm not entirely sure that this is the cause of the problem, though, so in order to rule it out I'd like to make sure that I'm actually integrating over the whole mesh even if it contains multiple regions. What is the correct idiom to do that?

P.S.: I hope it is appropriate to ask this in this thread. Or should I open a new question for this so that it is easier to find for others?

With the latest trunk you can do f*(dx(0)+dx(2)), but we still do not have
a natural notation for the entire domain. Maybe try something like
f*dx[mf](0), where mf is a CellFunction with the value 0 everywhere.

Martin
Den 2. des. 2012 17:45 skrev "Maximilian Albert" <
<email address hidden>> følgende:

> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Status: Solved => Open
>
> Maximilian Albert is still having a problem:
> Sorry to revive this old thread, but it is possible that I am being
> bitten by the "dx = dx(0)" issue again. I'm not entirely sure that this
> is the cause of the problem, though, so in order to rule it out I'd like
> to make sure that I'm actually integrating over the whole mesh even if
> it contains multiple regions. What is the correct idiom to do that?
>
> P.S.: I hope it is appropriate to ask this in this thread. Or should I
> open a new question for this so that it is easier to find for others?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Hi again,

apologies for being a nag, but we were just hit again by this problem and by the fact that there isn't a very good workaround. Part of the blame goes to Netgen because it adds the marker '1' to mesh cells by default.

But imho it's a big problem that there simply isn't anything about the syntax of 'df.dx' to imply that it only operates on part of the mesh. A few colleages of mine also stumbled over this and it took them quite a while to find the error. What's more, there are situations where we need to operate on the entire mesh even if it contains multiple regions. So it would be very much appreciated if at least an efficient shortcut for this case could be added. Ideally, it would be possible to use such a shortcut as the default argument for a function, e.g. like this:

def f(mesh, dx=df.dx('all')):
    assemble(...*dx, mesh=mesh)

I don't know anything about the internals of UFL, but if there is anything I can do to help nonetheless please let me know. If someone points me to the relevant bits of code I can also try to submit a patch (although I can't promise that I'll be able to get up to speed quickly enough to be able to do that).

Many thanks for everyone's help and great work!

I think that Martin may be looking at this. It is a problem that

    dx == dx(0)

and it is a frequent source of errors. It should be fixed.

Resolving this problem/feature is currently active work in progress, see:

  https://blueprints.launchpad.net/dolfin/+spec/integrals-over-everywhere

I'm very much working on this. Will post an explanation to the list when it is working.

Just to put some pressure on myself: I _hope_ to finish this tomorrow.

On Mon, Jan 28, 2013 at 03:21:05PM -0000, Martin Sandve Alnæs wrote:
> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Martin Sandve Alnæs posted a new comment:
> Just to put some pressure on myself: I _hope_ to finish this tomorrow.

Related update on the "redesign" of UFC I started in Cambridge which I
will need to merge at some point: currently compiling but changes spread
throughout UFC/FFC/DOLFIN and I'm looking forward to a painful
merge...

--
Anders

By "redesign", do you mean the cell related stuff or have you done more? I've mostly touched integral creation related parts. Hopefully the merge can be fairly smooth...

On Tue, Jan 29, 2013 at 08:41:29AM -0000, Martin Sandve Alnæs wrote:
> Question #212697 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/212697
>
> Martin Sandve Alnæs posted a new comment:
> By "redesign", do you mean the cell related stuff or have you done more?
> I've mostly touched integral creation related parts. Hopefully the merge
> can be fairly smooth...

I hope so too. What I've done is to add utility functions to
ufc_geometry.h and use them from the code generation. At the same
time, I started to replace ufc::cell in tabulate_tensor with a
flattened array of vertex coordinates, but the flattening and the
change of arguments spread throughout UFC, FFC and DOLFIN.

--
Anders

Awesome to see so much activity on this. And thanks again for everyone's hard work! I appreciate how difficult it is to put a software as huge and powerful as FEniCS together and how many issues need to be tackled at once. Great job, everyone!