Add a layer of hexahedrons to a 2D quad mesh
I have the following example to add a layer of hexahedrons to a structured mesh of quadrangles. Everything looks fine, but my physical group does not include the mesh:
import gmsh
import numpy as np
import sys
gmsh.initialize(sys.argv)
# Choose kernel
option = gmsh.option
model = gmsh.model
occ = model.occ
mesh = model.mesh
def split_list_by(l, N):
return [l[i : i + N] for i in range(0, len(l), N)]
def _get_nodes(dim, tag):
node_tags, coord, _ = mesh.getNodes(dim, tag, True, False)
return dict(zip(node_tags, split_list_by(coord, 3)))
def _get_elements(dim, tag):
element_types, element_tags, node_tags = mesh.getElements(dim, tag)
return {
t: zip(et, split_list_by(nt, element_type_map[t]["nodes"]))
for t, et, nt in zip(element_types, element_tags, node_tags)
}
def get_entitiy_data(dim=-1):
return {
e: {
"boundary": model.getBoundary([e]),
"nodes": _get_nodes(e[0], e[1]),
"elements": _get_elements(e[0], e[1]),
"type": model.getType(e[0], e[1]),
}
for e in model.getEntities(dim)
}
element_type_map = {
15: {"higher": 1, "nodes": 1, "dim": 0, "name": "point"},
1: {"higher": 3, "nodes": 2, "dim": 1, "name": "line"},
3: {"higher": 5, "nodes": 4, "dim": 2, "name": "quadrangle"},
5: {"higher": None, "nodes": 8, "dim": 3, "name": "hexahedron"},
}
def extrude_hexahedron(dx, dy, dz):
m = get_entitiy_data()
faces = model.getEntities(2)
occ.extrude(faces, dx, dy, dz)
occ.synchronize()
offset_node = 1e6
offset_element = 1e6
# cycle all entities sorted by dimension
for entity, v in sorted(m.items()):
up, down = model.getAdjacencies(*entity)
# add new higher dimensional entity
entity_dim = entity[0] + 1
entity_tag = up[0]
# cycle all elements in this entity
for element_type, tags in v["elements"].items():
if element_type == 3:
for element_tag, node_tags in tags:
element_coords = [v["nodes"][t] for t in node_tags]
shifted_element_coords = [
[n[0] + dx, n[1] + dy, n[2] + dz] for n in element_coords
]
new_coords = element_coords + shifted_element_coords
new_node_tags = np.concatenate(
(node_tags, [n + offset_node for n in node_tags])
)
mesh.addNodes(
entity_dim,
entity_tag,
new_node_tags,
[ce for c in new_coords for ce in c],
)
mesh.addElements(
entity_dim,
entity_tag,
[element_type_map[element_type]["higher"]], # higher order type
[[element_tag + offset_element]],
[new_node_tags],
)
occ.addRectangle(0, 0, 0, 1, 1)
occ.synchronize()
mesh.setTransfiniteAutomatic()
mesh.generate(2)
extrude_hexahedron(0, 0, 0.1)
model.addPhysicalGroup(2, [2], name="test")
gmsh.fltk.run()
gmsh.finalize()
What am I doing wrong?