MeshLine.load submeshes & boundaries

[gdmcbain]

MeshLine.load warns

UserWarning: Unable to load submeshes

for any input, and fails to read any subdomains and boundaries present.

The immediate cause is that 'line' is absent from the keys in

scikit-fem/skfem/mesh/mesh.py

Lines 286 to 291 in dadd4ba

	bnd_type = {
	'triangle' : 'line',
	'quad' : 'line',
	'tetra' : 'triangle',
	'hexahedron' : 'quad',
	}[self.meshio_type]

While the one-dimensional examples to date just use MeshLine.__init__, and this suffices for most purposes, it is occasionally useful to generate even one-dimensional meshes externally in Gmsh; e.g. for problems like the ‘composite wall composed of n slabs of thicknesses l1, …, ln and conductivities K1, …, Kn’ of Carslaw & Jaeger (1959, Conduction of Heat in Solids, §3.2 ‘Steady temperature’, Oxford University Press).

Minimal example: line.geo.gz.

Point(1) = {0, 0, 0};
Point(2) = {1, 0, 0};
Line(1) = {1, 2};

Generate mesh with

gmsh -1 line.geo

and attempt to load with

MeshLine.load('line.msh')

[gdmcbain]

Just adding 'line': 'vertex' seems to be enough to get MeshLine.load to not warn and to successfully read the subdomain defined by the addition line in line.geo

Physical Curve("domain") = {1};

[gdmcbain]

MeshLine.load is not loading the boundaries correctly; e.g. for

Point(1) = {0, 0, 0};
Point(2) = {1, 0, 0};
Line(1) = {1, 2};
Physical Point("left") = {1};
Physical Point("right") = {2};
Physical Curve("domain") = {1};

which produces line.msh.gz, containing

$PhysicalNames
3
0 1 "left"
0 2 "right"
1 3 "domain"
$EndPhysicalNames

but only 'left' and 'domain' are loaded, not 'right'.

I'm looking into Mesh._parse_submeshes now.

[gdmcbain]

I'm wondering whether the problem follows from MeshLine._build_mappings sorting .p whereas that method of other descendents, e.g. MeshTri, MeshTet, doesn't.

[gdmcbain]

Yes, this is the problem. In the above example, at

scikit-fem/skfem/mesh/mesh.py

Lines 314 to 316 in 9ff1fd5

	facets = self.external.cells[bnd_type]
	facets_tag = self.external.cell_data[bnd_type]['gmsh:physical']
	bndfacets = self.boundary_facets()

we have facets == array([[0], [1]]), but bndfacets == array([0, 5]), i.e. facets is numbered as in line.msh but bndfacets uses the reordered nodal indices.

Then dic == {{(0,): 1, (1,): 2} (input order) so in

scikit-fem/skfem/mesh/mesh.py

Lines 324 to 326 in 9ff1fd5

	index = np.array([[dic[tuple(np.sort(self.facets[:, i]))], i]
	for i in bndfacets
	if tuple(np.sort(self.facets[:, i])) in dic])

the if-clause does get 'left' but only by luck since node 0 survived the permutation but misses 'right' since self.external.points[1] at x = 1.0 has been renumbered as self.p[:, 5]; i.e. index only contains one entry

array([[1, 0]])

I don't know how to fix this. I don't think it can be done at this point of _parse_submeshes since the permutation of the points is not remembered.

[kinnala]

I can't recall now the exact reason for doing the sort in _build_mappings. There could have been some technical problem but maybe its only because mesh.p was easier to interpret this way.

[kinnala]

I think you could as a workaround (for 1D meshes) re-sort self.external.points in _parse_submeshes and by that way obtain the correspondence 0<->0, 1<->5 by setting the kwarg return_index or return_inverse to True in np.sort.

[gdmcbain]

I can't recall now the exact reason for doing the sort in _build_mappings. There could have been some technical problem but maybe its only because mesh.p was easier to interpret this way.

That sorting is relied upon in a couple of places:

• the argument t that is passed to MeshLine.__init__ is always silently ignored and
• MeshLine.t and MeshLine.facets are defined in terms of the sorted MeshLine.p.

scikit-fem/skfem/mesh/mesh_line.py

Lines 44 to 47 in 625ac42

	self.p = np.sort(self.p)
	self.t = np.array([np.arange(np.max(self.p.shape)-1), np.arange(np.max(self.p.shape)-1)+1])
	self.facets = np.array([np.arange(self.p.shape[1])])

[gdmcbain]

So if it's O. K. to drop the sorting of the p, there's a fix in #176.

[kinnala]

Nothing left here I presume.