# Error in defining Measure with subdomains

**URL:** <https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387>\
**Category:** variational formulation\
**Created:** [October 24, 2020, 11:07am UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387 "2020-10-24T11:07:23Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Iactor](https://avatars.discourse-cdn.com/v4/letter/i/6f9a4e/32.png) [@Iactor](https://fenicsproject.discourse.group/u/Iactor)\
**Post date:** [October 24, 2020, 11:07am UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/1 "2020-10-24T11:07:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 24, 2020, 11:16am UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/2 "2020-10-24T11:16:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Iactor](https://avatars.discourse-cdn.com/v4/letter/i/6f9a4e/32.png) [@Iactor](https://fenicsproject.discourse.group/u/Iactor)\
**Post date:** [October 24, 2020, 3:06pm UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/3 "2020-10-24T15:06:03Z")

</div>

Hi, after importing

```auto
import dolfin as dol
import meshio
import numpy as np

```

I convert a .msh file to XDMF using

```auto
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.

```auto
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,

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

```

does give the total area (5.75)  
but

```auto
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.

```auto
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};

```

---

<div class="post-metadata">

**Author:** ![Iactor](https://avatars.discourse-cdn.com/v4/letter/i/6f9a4e/32.png) [@Iactor](https://fenicsproject.discourse.group/u/Iactor)\
**Post date:** [October 25, 2020, 2:03pm UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/4 "2020-10-25T14:03:43Z")

</div>

Interestingly if instead of multiplying the measure by a dolfin constant

```auto
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)

```auto
dol.assemble(1 * dx(1))

```

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

```auto
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?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 25, 2020, 2:23pm UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/5 "2020-10-25T14:23:48Z")

</div>

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

```auto
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:

```auto
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

```

---

<div class="post-metadata">

**Author:** ![Iactor](https://avatars.discourse-cdn.com/v4/letter/i/6f9a4e/32.png) [@Iactor](https://fenicsproject.discourse.group/u/Iactor)\
**Post date:** [October 25, 2020, 5:03pm UTC](https://fenicsproject.discourse.group/t/error-in-defining-measure-with-subdomains/4387/6 "2020-10-25T17:03:40Z")

</div>

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