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