Assemble is not working with fenics_adjoint and imported mesh with MESHIO

I can not get my subdomains (and boundaries) to work when I import fenics_adjoint. This post is close, but the subdomains are not included here. Dolfin_adjoint, overloaded

I have included a MWE, where I’m just trying to calculate the the volumes. The mesh is converted using msh2xdmf /github.com/floiseau/msh2xdmf, modifed to work with the newest version of meshio according to msh2xdmf issue.

Everything seems fine with dolfin_convert, but i prefer to use meshio.

from fenics import *
from fenics_adjoint import * # It is working without this
from msh2xdmf_modified import import_mesh
mesh, boundaries_mf, subdomains_mf, association_table = import_mesh(prefix='MM2_1_4' ,dim=2 ,subdomains=True)
# mesh = Mesh(mesh) # This is not working either!
dx = Measure('dx', domain = mesh, subdomain_data = subdomains_mf)
# Calculate Volumes of Different sections
voltot = assemble(Constant(1)*dx)      
volMat1 = assemble(Constant(1)*dx(10))   
volMat2 = assemble(Constant(1)*dx(11))

The MSH file

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
4
1 8 "Bottom"
1 9 "Top"
2 10 "Material1"
2 11 "Material2"
$EndPhysicalNames
$Entities
6 7 2 0
1 0 0 0 0 
2 0 0.5 0 0 
3 0 1 0 0 
4 1 1 0 0 
5 1 0.5 0 0 
6 1 0 0 0 
1 0 0 0 0 0.5 0 0 2 1 -2 
2 0 0.5 0 0 1 0 0 2 2 -3 
3 0 1 0 1 1 0 1 9 2 3 -4 
4 1 0.5 0 1 1 0 0 2 4 -5 
5 1 0 0 1 0.5 0 0 2 6 -5 
6 0 0 0 1 0 0 1 8 2 1 -6 
7 0 0.5 0 1 0.5 0 0 2 2 -5 
1 0 0 0 1 0.5 0 1 10 4 1 7 -5 -6 
2 0 0.5 0 1 1 0 1 11 4 2 3 4 -7 
$EndEntities
$Nodes
15 121 1 121
0 1 0 1
1
0 0 0
0 2 0 1
2
0 0.5 0
0 3 0 1
3
0 1 0
0 4 0 1
4
1 1 0
0 5 0 1
5
1 0.5 0
0 6 0 1
6
1 0 0
1 1 0 4
7
8
9
10
0 0.09999999999977846 0
0 0.1999999999994866 0
0 0.2999999999994725 0
0 0.3999999999997362 0
1 2 0 4
11
12
13
14
0 0.6 0
0 0.7 0
0 0.8 0
0 0.9 0
1 3 0 9
15
16
17
18
19
20
21
22
23
0.09999999999981414 1 0
0.1999999999995569 1 0
0.299999999999265 1 0
0.3999999999989731 1 0
0.4999999999986921 1 0
0.599999999998945 1 0
0.6999999999992088 1 0
0.7999999999994725 1 0
0.8999999999997362 1 0
1 4 0 4
24
25
26
27
1 0.9 0
1 0.8 0
1 0.7 0
1 0.6 0
1 5 0 4
28
29
30
31
1 0.09999999999977846 0
1 0.1999999999994866 0
1 0.2999999999994725 0
1 0.3999999999997362 0
1 6 0 9
32
33
34
35
36
37
38
39
40
0.09999999999981414 0 0
0.1999999999995569 0 0
0.299999999999265 0 0
0.3999999999989731 0 0
0.4999999999986921 0 0
0.599999999998945 0 0
0.6999999999992088 0 0
0.7999999999994725 0 0
0.8999999999997362 0 0
1 7 0 9
41
42
43
44
45
46
47
48
49
0.09999999999981414 0.5 0
0.1999999999995569 0.5 0
0.299999999999265 0.5 0
0.3999999999989731 0.5 0
0.4999999999986921 0.5 0
0.599999999998945 0.5 0
0.6999999999992088 0.5 0
0.7999999999994725 0.5 0
0.8999999999997362 0.5 0
2 1 0 36
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
0.09999999999981413 0.09999999999977846 0
0.1999999999995569 0.09999999999977847 0
0.299999999999265 0.09999999999977846 0
0.3999999999989731 0.09999999999977846 0
0.499999999998692 0.09999999999977843 0
0.5999999999989449 0.09999999999977846 0
0.6999999999992088 0.09999999999977846 0
0.7999999999994725 0.09999999999977846 0
0.8999999999997361 0.09999999999977846 0
0.09999999999981413 0.1999999999994866 0
0.1999999999995569 0.1999999999994866 0
0.299999999999265 0.1999999999994865 0
0.3999999999989731 0.1999999999994866 0
0.499999999998692 0.1999999999994865 0
0.5999999999989449 0.1999999999994866 0
0.6999999999992088 0.1999999999994865 0
0.7999999999994725 0.1999999999994866 0
0.8999999999997361 0.1999999999994866 0
0.09999999999981414 0.2999999999994724 0
0.1999999999995569 0.2999999999994726 0
0.299999999999265 0.2999999999994725 0
0.3999999999989731 0.2999999999994726 0
0.499999999998692 0.2999999999994726 0
0.599999999998945 0.2999999999994725 0
0.6999999999992088 0.2999999999994725 0
0.7999999999994725 0.2999999999994725 0
0.8999999999997362 0.2999999999994725 0
0.09999999999981415 0.3999999999997362 0
0.1999999999995569 0.3999999999997362 0
0.299999999999265 0.3999999999997362 0
0.3999999999989731 0.3999999999997362 0
0.4999999999986921 0.3999999999997362 0
0.599999999998945 0.3999999999997362 0
0.6999999999992088 0.3999999999997362 0
0.7999999999994725 0.3999999999997362 0
0.8999999999997362 0.3999999999997362 0
2 2 0 36
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
0.09999999999981414 0.6000000000000002 0
0.1999999999995569 0.5999999999999998 0
0.299999999999265 0.6000000000000001 0
0.3999999999989731 0.6000000000000001 0
0.4999999999986921 0.5999999999999998 0
0.599999999998945 0.6 0
0.6999999999992088 0.6 0
0.7999999999994725 0.6 0
0.8999999999997362 0.6 0
0.09999999999981414 0.7000000000000002 0
0.1999999999995569 0.7000000000000001 0
0.299999999999265 0.7000000000000001 0
0.3999999999989731 0.7000000000000001 0
0.499999999998692 0.7 0
0.599999999998945 0.7 0
0.6999999999992088 0.7 0
0.7999999999994725 0.7 0
0.8999999999997362 0.7 0
0.09999999999981414 0.8 0
0.1999999999995569 0.8 0
0.299999999999265 0.7999999999999998 0
0.3999999999989731 0.8000000000000002 0
0.499999999998692 0.8000000000000002 0
0.599999999998945 0.8 0
0.6999999999992088 0.8 0
0.7999999999994725 0.8 0
0.8999999999997362 0.8 0
0.09999999999981414 0.9000000000000001 0
0.1999999999995569 0.9 0
0.299999999999265 0.9000000000000001 0
0.3999999999989731 0.8999999999999998 0
0.4999999999986921 0.9000000000000001 0
0.599999999998945 0.9 0
0.6999999999992088 0.8999999999999998 0
0.7999999999994725 0.9 0
0.8999999999997362 0.8999999999999997 0
$EndNodes
$Elements
4 220 1 220
1 3 1 10
1 3 15 
2 15 16 
3 16 17 
4 17 18 
5 18 19 
6 19 20 
7 20 21 
8 21 22 
9 22 23 
10 23 4 
1 6 1 10
11 1 32 
12 32 33 
13 33 34 
14 34 35 
15 35 36 
16 36 37 
17 37 38 
18 38 39 
19 39 40 
20 40 6 
2 1 2 100
21 1 7 50 
22 50 32 1 
23 32 50 51 
24 51 33 32 
25 33 51 52 
26 52 34 33 
27 34 52 53 
28 53 35 34 
29 35 53 54 
30 54 36 35 
31 36 54 55 
32 55 37 36 
33 37 55 56 
34 56 38 37 
35 38 56 57 
36 57 39 38 
37 39 57 58 
38 58 40 39 
39 40 58 28 
40 28 6 40 
41 7 8 59 
42 59 50 7 
43 50 59 60 
44 60 51 50 
45 51 60 61 
46 61 52 51 
47 52 61 62 
48 62 53 52 
49 53 62 63 
50 63 54 53 
51 54 63 64 
52 64 55 54 
53 55 64 65 
54 65 56 55 
55 56 65 66 
56 66 57 56 
57 57 66 67 
58 67 58 57 
59 58 67 29 
60 29 28 58 
61 8 9 68 
62 68 59 8 
63 59 68 69 
64 69 60 59 
65 60 69 70 
66 70 61 60 
67 61 70 71 
68 71 62 61 
69 62 71 72 
70 72 63 62 
71 63 72 73 
72 73 64 63 
73 64 73 74 
74 74 65 64 
75 65 74 75 
76 75 66 65 
77 66 75 76 
78 76 67 66 
79 67 76 30 
80 30 29 67 
81 9 10 77 
82 77 68 9 
83 68 77 78 
84 78 69 68 
85 69 78 79 
86 79 70 69 
87 70 79 80 
88 80 71 70 
89 71 80 81 
90 81 72 71 
91 72 81 82 
92 82 73 72 
93 73 82 83 
94 83 74 73 
95 74 83 84 
96 84 75 74 
97 75 84 85 
98 85 76 75 
99 76 85 31 
100 31 30 76 
101 10 2 41 
102 41 77 10 
103 77 41 42 
104 42 78 77 
105 78 42 43 
106 43 79 78 
107 79 43 44 
108 44 80 79 
109 80 44 45 
110 45 81 80 
111 81 45 46 
112 46 82 81 
113 82 46 47 
114 47 83 82 
115 83 47 48 
116 48 84 83 
117 84 48 49 
118 49 85 84 
119 85 49 5 
120 5 31 85 
2 2 2 100
121 2 11 86 
122 86 41 2 
123 41 86 87 
124 87 42 41 
125 42 87 88 
126 88 43 42 
127 43 88 89 
128 89 44 43 
129 44 89 90 
130 90 45 44 
131 45 90 91 
132 91 46 45 
133 46 91 92 
134 92 47 46 
135 47 92 93 
136 93 48 47 
137 48 93 94 
138 94 49 48 
139 49 94 27 
140 27 5 49 
141 11 12 95 
142 95 86 11 
143 86 95 96 
144 96 87 86 
145 87 96 97 
146 97 88 87 
147 88 97 98 
148 98 89 88 
149 89 98 99 
150 99 90 89 
151 90 99 100 
152 100 91 90 
153 91 100 101 
154 101 92 91 
155 92 101 102 
156 102 93 92 
157 93 102 103 
158 103 94 93 
159 94 103 26 
160 26 27 94 
161 12 13 104 
162 104 95 12 
163 95 104 105 
164 105 96 95 
165 96 105 106 
166 106 97 96 
167 97 106 107 
168 107 98 97 
169 98 107 108 
170 108 99 98 
171 99 108 109 
172 109 100 99 
173 100 109 110 
174 110 101 100 
175 101 110 111 
176 111 102 101 
177 102 111 112 
178 112 103 102 
179 103 112 25 
180 25 26 103 
181 13 14 113 
182 113 104 13 
183 104 113 114 
184 114 105 104 
185 105 114 115 
186 115 106 105 
187 106 115 116 
188 116 107 106 
189 107 116 117 
190 117 108 107 
191 108 117 118 
192 118 109 108 
193 109 118 119 
194 119 110 109 
195 110 119 120 
196 120 111 110 
197 111 120 121 
198 121 112 111 
199 112 121 24 
200 24 25 112 
201 14 3 15 
202 15 113 14 
203 113 15 16 
204 16 114 113 
205 114 16 17 
206 17 115 114 
207 115 17 18 
208 18 116 115 
209 116 18 19 
210 19 117 116 
211 117 19 20 
212 20 118 117 
213 118 20 21 
214 21 119 118 
215 119 21 22 
216 22 120 119 
217 120 22 23 
218 23 121 120 
219 121 23 4 
220 4 24 121 
$EndElements

Modified msh2xdmf

#!/usr/bin/env python

import argparse
import meshio
import os
import numpy as np
from configparser import ConfigParser
try:
    from dolfin import XDMFFile, Mesh, MeshValueCollection
    from dolfin.cpp.mesh import MeshFunctionSizet
except ImportError:
    print("Could not import dolfin. Continuing without Dolfin support.")


def msh2xdmf(mesh_name, dim=2, directory="."):
    """
    Function converting a MSH mesh into XDMF files.
    The XDMF files are:
        - "domain.xdmf": the domain;
        - "boundaries.xdmf": the boundaries physical groups from GMSH;
    """
    # Get the mesh name has prefix
    prefix = mesh_name.split('.')[0]
    # Read the input mesh
    msh = meshio.read("{}/{}".format(directory, mesh_name))
    # Generate the domain XDMF file
    export_domain(msh, dim, directory, prefix)
    # Generate the boundaries XDMF file
    export_boundaries(msh, dim, directory, prefix)
    # Export association table
    export_association_table(msh, prefix, directory)


def export_domain(msh, dim, directory, prefix):
    """
    Export the domain XDMF file as well as the subdomains values.
    """
    # Set cell type
    if dim == 2:
        cell_type = "triangle"
    elif dim == 3:
        cell_type = "tetra"
    # Generate the cell block for the domain cells
    data_array = [cell.data for cell in msh.cells if cell.type == cell_type]
    if len(data_array) == 0:
        print("WARNING: No domain physical group found.")
        return
    else:
        data = np.concatenate(data_array)
    cells = [
        meshio.CellBlock(
            cell_type=cell_type,
            data=data,
        )
    ]
    # Generate the domain cells data (for the subdomains)
    try:
        cell_data = {
            "subdomains": [
                np.concatenate(
                    [
                        msh.cell_data["gmsh:physical"][i]
                        for i, cellBlock in enumerate(msh.cells)
                        if cellBlock.type == cell_type
                    ]
                )
            ]
        }
    except KeyError:
        raise ValueError(
            """
            No physical group found for the domain.
            Define the domain physical group.
                - if dim=2, the domain is a surface
                - if dim=3, the domain is a volume
            """
        )

    # Generate a meshio Mesh for the domain
    domain = meshio.Mesh(
        points=msh.points[:, :dim],
        cells=cells,
        cell_data=cell_data,
    )
    # Export the XDMF mesh of the domain
    meshio.write(
        "{}/{}_{}".format(directory, prefix, "domain.xdmf"),
        domain,
        file_format="xdmf"
    )


def export_boundaries(msh, dim, directory, prefix):
    """
    Export the boundaries XDMF file.
    """
    # Set the cell type
    if dim == 2:
        cell_type = "line"
    elif dim == 3:
        cell_type = "triangle"
    # Generate the cell block for the boundaries cells
    data_array = [cell.data for cell in msh.cells if cell.type == cell_type]
    if len(data_array) == 0:
        print("WARNING: No boundary physical group found.")
        return
    else:
        data = np.concatenate(data_array)
    boundaries_cells = [
        meshio.CellBlock(
            cell_type=cell_type,
            data=data,
        )
    ]
    # Generate the boundaries cells data
    cell_data = {
        "boundaries": [
            np.concatenate(
                [
                    msh.cell_data["gmsh:physical"][i]
                    for i, cellBlock in enumerate(msh.cells)
                    if cellBlock.type == cell_type
                ]
            )
        ]
    }
    # Generate the meshio Mesh for the boundaries physical groups
    boundaries = meshio.Mesh(
        points=msh.points[:, :dim],
        cells=boundaries_cells,
        cell_data=cell_data,
    )
    # Export the XDMF mesh of the lines boundaries
    meshio.write(
        "{}/{}_{}".format(directory, prefix, "boundaries.xdmf"),
        boundaries,
        file_format="xdmf"
    )


def export_association_table(msh, prefix='mesh', directory='.', verbose=True):
    """
    Display the association between the physical group label and the mesh
    value.
    """
    # Create association table
    association_table = {}

    # Display the correspondance
    formatter = "|{:^20}|{:^20}|"
    topbot = "+{:-^41}+".format("")
    separator = "+{:-^20}+{:-^20}+".format("", "")

    # Display
    if verbose:
        print('\n' + topbot)
        print(formatter.format("GMSH label", "MeshFunction value"))
        print(separator)

    for label, arrays in msh.cell_sets.items():
        # Get the index of the array in arrays
        for i, array in enumerate(arrays):
            if array.size != 0:
                index = i
        # Added check to make sure that the association table
        # doesn't try to import irrelevant information.
        if label != "gmsh:bounding_entities":
            value = msh.cell_data["gmsh:physical"][index][0]
            # Store the association table in a dictionnary
            association_table[label] = value
            # Display the association
            if verbose:
                print(formatter.format(label, value))
    if verbose:
        print(topbot)
    # Export the association table
    file_content = ConfigParser()
    file_content["ASSOCIATION TABLE"] = association_table
    file_name = "{}/{}_{}".format(directory, prefix, "association_table.ini")
    with open(file_name, 'w') as f:
        file_content.write(f)


def import_mesh(
        prefix="mesh",
        subdomains=False,
        dim=2,
        directory=".",
):
    """Function importing a dolfin mesh.

    Arguments:
        prefix (str, optional): mesh files prefix (eg. my_mesh.msh,
            my_mesh_domain.xdmf, my_mesh_bondaries.xdmf). Defaults to "mesh".
        subdomains (bool, optional): True if there are subdomains. Defaults to
            False.
        dim (int, optional): dimension of the domain. Defaults to 2.
        directory (str, optional): directory of the mesh files. Defaults to ".".

    Output:
        - dolfin Mesh object containing the domain;
        - dolfin MeshFunction object containing the physical lines (dim=2) or
            surfaces (dim=3) defined in the msh file and the sub-domains;
        - association table
    """
    # Set the file name
    domain = "{}_domain.xdmf".format(prefix)
    boundaries = "{}_boundaries.xdmf".format(prefix)

    # create 2 xdmf files if not converted before
    if not os.path.exists("{}/{}".format(directory, domain)) or \
       not os.path.exists("{}/{}".format(directory, boundaries)):
        msh2xdmf("{}.msh".format(prefix), dim=dim, directory=directory)

    # Import the converted domain
    mesh = Mesh()
    with XDMFFile("{}/{}".format(directory, domain)) as infile:
        infile.read(mesh)
    # Import the boundaries
    boundaries_mvc = MeshValueCollection("size_t", mesh, dim=dim)
    with XDMFFile("{}/{}".format(directory, boundaries)) as infile:
        infile.read(boundaries_mvc, 'boundaries')
    boundaries_mf = MeshFunctionSizet(mesh, boundaries_mvc)
    # Import the subdomains
    if subdomains:
        subdomains_mvc = MeshValueCollection("size_t", mesh, dim=dim)
        with XDMFFile("{}/{}".format(directory, domain)) as infile:
            infile.read(subdomains_mvc, 'subdomains')
        subdomains_mf = MeshFunctionSizet(mesh, subdomains_mvc)
    # Import the association table
    association_table_name = "{}/{}_{}".format(
        directory, prefix, "association_table.ini")
    file_content = ConfigParser()
    file_content.read(association_table_name)
    association_table = dict(file_content["ASSOCIATION TABLE"])
    # Convert the value from string to int
    for key, value in association_table.items():
        association_table[key] = int(value)
    # Return the Mesh and the MeshFunction objects
    if not subdomains:
        return mesh, boundaries_mf, association_table
    else:
        return mesh, boundaries_mf, subdomains_mf, association_table


if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "msh_file",
        help="input .msh file",
        type=str,
    )
    parser.add_argument(
        "-d",
        "--dimension",
        help="dimension of the domain",
        type=int,
        default=2,
    )
    args = parser.parse_args()
    # Get current directory
    current_directory = os.getcwd()
    # Conert the mesh
    msh2xdmf(args.msh_file, args.dimension, directory=current_directory)```