Error in defining Measure with subdomains

variational formulation
1

Hi
I defined a volume measure

measure

using subdomain data from an external mesh.
If I check the subdomains by plotting

dol.plot(measure.subdomain_data())

I obtain the plot I would expect with the subdomains marked correctly (one with id 1 and another one with id 2)
However if I try

print(dol.assemble(dol.Constant(1) * measure(1)))

I get 0, while

print(dol.assemble(dol.Constant(1) * measure))

gives the total volume.
Does anyone have an idea of what can be wrong? (i am using dolfin 2019.1.0)
Thanks in advance

2

Please create a minimal working code example, that reproduces your error.

3

Hi, after importing

import dolfin as dol
import meshio
import numpy as np

I convert a .msh file to XDMF using

msh = meshio.read('device.msh')
volume_cell_type = 'triangle' 
volume_cells = []

for cell in msh.cells:
    if cell.type == volume_cell_type:
        if len(volume_cells) == 0:
            volume_cells = cell.data
        else:
            volume_cells = np.vstack([volume_cells, cell.data])
        
for key in msh.cell_data_dict['gmsh:physical'].keys():
    if key == volume_cell_type:
        volume_data = msh.cell_data_dict['gmsh:physical'][key]

volume_mesh = meshio.Mesh(points = msh.points[:,:2],
                          cells = [(volume_cell_type, volume_cells)],
                          cell_data = {'marker' : [volume_data]})

meshio.write('device_2D_mesh.xdmf', volume_mesh)

I then open the file, read the mesh and the marker, and define the measure.

mesh = dol.Mesh()
with dol.XDMFFile('device_2D_mesh.xdmf') as infile:
    infile.read(mesh)
mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = 2) 
with dol.XDMFFile('device_2D_mesh.xdmf') as infile:
    infile.read(mvc, 'marker')
volume_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
del(mvc)
dx = dol.Measure('dx', domain = mesh, subdomain_data = volume_marker)

I plotted volume_marker using dolfin plot function and it is correct,

dol.assemble(dol.Constant(1) * dx )

does give the total area (5.75)
but

dol.assemble(dol.Constant(1) * dx(1))
dol.assemble(dol.Constant(1) * dx(2))

Both return 0 instead of 4.5 and 1.25.
I am using
dolfin 2019.1.0
meshio 4.2.0
gmsh 4.5.5 for producing the msh file
For completeness I also attach the geometry file I used to produce the .msh one.

SetFactory("OpenCASCADE");

lc = .1;

Point(1) = {0, 0, 0, lc};
Point(2) = {4, 0, 0, lc};
Point(3) = {3, 2, 0, lc/2};
Point(4) = {1, 1, 0, lc/2};
Point(5) = {0.5, 0.5, 0, lc};

Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,5};
Line(5) = {5,1};

Line Loop(1) = {1,2,3,4,5};
Plane Surface(1) = {1};

Point(6) = {1, 2, 0, lc};
Point(7) = {0, 2, 0, lc};
Point(8) = {0, 1, 0, lc};

Line(6) = {4,6};
Line(7) = {6,7};
Line(8) = {7,8};
Line(9) = {8,5};

Line Loop(2) = {4,6,7,8,9};
Plane Surface(2) = {2};

Point(9) = {2, 1, 0, lc / 5};
Point(10) = {3, 0.5, 0, lc / 5};
Line(10) = {9,10};
Line{10} In Surface {1};

v() = BooleanFragments{Surface{1}; Delete;}{ Surface{2}; Delete; };

Physical Surface('blg', 1) = {1};
Physical Surface('slg', 2) = {2};

Physical Line('source', 1) = {7,8};
Physical Line('drain', 2) = {2};
Physical Line('wall', 3) = {10};
4

Interestingly if instead of multiplying the measure by a dolfin constant

dol.assemble(dol.Constant(1)  * dx(1))

I multiply it by a python int or float (I wasn’t aware of this possibility but I saw it reading the Measure code)

dol.assemble(1  * dx(1))

I get the correct result.
Moreover, while the second options returns the message

Calling FFC just-in-time (JIT) compiler, this may take some time.

and takes about 0.018 s on my pc the first one doesn’t display messages and takes much less (0.005 s).
Any idea of what is going wrong?

5

I cannot reproduce this with your mesh using the latest docker image: quay.io/fenicsproject/stable:latest

import dolfin as dol
mesh = dol.Mesh()
mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = 2)
with dol.XDMFFile('device_2D_mesh.xdmf') as infile:
    infile.read(mesh)
    infile.read(mvc, 'marker')
volume_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)

del(mvc)
dx = dol.Measure('dx', domain = mesh, subdomain_data = volume_marker)

print(dol.assemble(1*dx(1)), dol.assemble(dol.Constant(1)*dx(1)))
print(dol.assemble(1*dx(2)), dol.assemble(dol.Constant(1)*dx(2)))

With output:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
4.499999999999997 4.499999999999997
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
1.2500000000000009 1.2500000000000009
6

Indeed the issue disappears with 2019.2.0.dev0 I just reinstalled fro apt-get