Problem with projection (mesh domains ?)

Asked by Stepan Roucka on 2013-02-25

Dear all,
while trying to run a fenics model, which used to work about two years ago, I have found the following problem: I generate diffpack meshes using netgen and convert them to .xml using dolfin-convert. The meshes produced by the recent version of dolfin-convert contain a new section <domains>, which looks like this:
 ...
<tetrahedron index="32764" v0="972" v1="4739" v2="4754" v3="4788"/>
      <tetrahedron index="32765" v0="4747" v1="6939" v2="6941" v3="6942"/>
      <tetrahedron index="32766" v0="6161" v1="6201" v2="6202" v3="6203"/>
    </cells>
    <domains>
    <mesh_value_collection name="m" type="uint" dim="2" size="8936">
      <value cell_index="5" local_entity="3" value="3"/>
      <value cell_index="18" local_entity="3" value="3"/>
      <value cell_index="46" local_entity="3" value="3"/>
      <value cell_index="52" local_entity="3" value="3"/>
      <value cell_index="63" local_entity="3" value="3"/>
      <value cell_index="65" local_entity="0" value="3"/>
      <value cell_index="65" local_entity="1" value="3"/>
      <value cell_index="65" local_entity="2" value="3"/>
....
      <value cell_index="32766" local_entity="0" value="1"/>
    </mesh_value_collection>
    </domains>
  </mesh>
</dolfin>

The problem is that I am not able to define and project an Expression on such mesh. This can be demonstrated by the following minimal example, where pipe.xml and pipe2.xml are meshes without and with the <domains> section respectively. The projected field v2 is zero everywhere, while v1 is the expected parabolic velocity profile. The example meshes can be found at
http://www.roucka.eu/pub/pipe.xml
http://www.roucka.eu/pub/pipe2.xml

from dolfin import *

#pipe radius [m]
radius = 0.0094

# Load mesh and subdomains
mesh1 = Mesh("pipe.xml")
mesh2 = Mesh("pipe2.xml")

# Define function spaces
V1 = VectorFunctionSpace(mesh1, "CG", 2)
V2 = VectorFunctionSpace(mesh2, "CG", 2)

# Create velocity profile function
sqradius = radius**2
profile = "2.0*(1.0-(x[0]*x[0]+x[1]*x[1])/(" + str(sqradius) + "))"
v1 = project(Expression(("0.0", "0.0", profile)), V1)
v2 = project(Expression(("0.0", "0.0", profile)), V2)

# this shows that both meshes are loaded correctly?
plot(mesh1)
plot(mesh2)

# but the projections are not identical.
plot(v1)
plot(v2)
interactive()
exit()

I can solve the problem by deleting the <domains> from my meshes, but I believe that it is not the correct solution. Anyone knows what am I doing wrong?

Thanks in advance for any help.

Stepan

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
2013-02-25
Last query:
2013-02-25
Last reply:
2013-02-25
Stepan Roucka (rouckas) said : #1

I should add that I am running dolfin 1.1.0-1~ppa3~precise1 from the fenics ppa in ubuntu 12.04

Best Johan Hake (johan-hake) said : #2

Most probably there are cell markers, which are not equal to 0, stored
in the mesh. This in turn makes the dx integral return zeros. You can
either upgrade your dolfin installation, as this should have been fixed
in trunk, or you can try to clear the markers.

  mesh.domains().markers(3).clear()

Johan

On 02/25/2013 05:25 PM, Stepan Roucka wrote:
> The problem is that I am not able to define and project an Expression
> on such mesh. This can be demonstrated by the following minimal
> example, where pipe.xml and pipe2.xml are meshes without and with the
> <domains> section respectively. The projected field v2 is zero
> everywhere, while v1 is the expected parabolic velocity profile. The
> example meshes can be found at

Jan Blechta (blechta) said : #3

> I can solve the problem by deleting the <domains> from my meshes,
> but I believe that it is not the correct solution. Anyone
> knows what am I doing wrong?

This a good solution as long as you don't need these markers - it will save you some cost when loading mesh. The best solution would be to instruct your mesh generator not to write them to your original mesh files.

Stepan Roucka (rouckas) said : #4

Thanks Johan Hake, that solved my question.

Stepan Roucka (rouckas) said : #5

Thanks to both of you for the answers.

> This a good solution as long as you don't need these markers - it will save you some cost when
> loading mesh. The best solution would be to instruct your mesh generator not to write them
> to your original mesh files.

I am using netgen CSG and I don't know how disable (or change) the cell markers. I can only change the boundary markers AFAIK.