API: gmsh.model.mesh.getElementByCoordinates() crashes for higher order elements and meshes
Hello,
It seems there is a bug inside the definition of getElementByCoordinates() in the gmsh (Python)-API. This function only works for first order elements and somehow for second order tetrahedrons. In case of second order hexahedrons and second order prisms it crashes completely. I tested it with the code below that generates some minimal examples for the different second order element types. Maybe it is possible to fix this in future releases.
Regards
# -*- coding: utf-8 -*-
import gmsh
import numpy as np
model = gmsh.model
mesh = model.mesh
# Initialize a new gmsh model
gmsh.initialize()
entityTag = model.addDiscreteEntity(dim=3)
# Element type to test
# 'TET', 'HEX', 'PRISM'
test = 'PRISM'
if test == 'HEX':
# Minimum exaple for second order hex elements
# Add nodes to the entity
nodeCoords = np.array([[0,0,0], [2,0,0], [2,2,0], [0,2,0],
[0,0,2], [2,0,2], [2,2,2], [0,2,2],
[1,0,0], [0,1,0], [0,0,1], [2,1,0],
[2,0,1], [1,2,0], [2,2,1], [0,2,1],
[1,0,2], [0,1,2], [2,1,2], [1,2,2]], dtype=float)
nodeTags = np.arange(1,len(nodeCoords)+1)
mesh.addNodes(3, entityTag, nodeTags, nodeCoords.flatten())
# Add one 20-node second order hexahedron
mesh.addElementsByType(entityTag, 17, [1], nodeTags)
# Try getElementByCoordinates()
x, y, z = 1, 1, 1
returnStuff = mesh.getElementByCoordinates(x, y, z) # This will crash
elif test == 'TET':
# Minimum exaple for second order tet elements
# Add nodes to the entity
nodeCoords = np.array([[0,0,0], [2,0,0], [0,2,0], [0,0,2],
[1,0,0], [1,1,0], [0,1,0], [0,0,1],
[0,1,1], [1,0,1]], dtype=float)
nodeTags = np.arange(1,len(nodeCoords)+1)
mesh.addNodes(3, entityTag, nodeTags, nodeCoords.flatten())
# Add one 10-node second order tetrahedron
mesh.addElementsByType(entityTag, 11, [1], nodeTags)
# Try getElementByCoordinates()
x, y, z = 1, 1, 1
returnStuff = mesh.getElementByCoordinates(x, y, z) # This won't crash
elif test == 'PRISM':
# Minimum exaple for second order prism elements
# Add nodes to the entity
nodeCoords = np.array([[0,0,0], [2,0,0], [0,2,0], [0,0,2],
[2,0,2], [0,2,2], [1,0,0], [0,1,0],
[0,0,1], [1,1,0], [2,0,1], [0,2,1],
[1,0,2], [0,1,2], [1,1,2]], dtype=float)
nodeTags = np.arange(1,len(nodeCoords)+1)
mesh.addNodes(3, entityTag, nodeTags, nodeCoords.flatten())
# Add one 15-node second order prism
mesh.addElementsByType(entityTag, 18, [1], nodeTags)
# Try getElementByCoordinates()
x, y, z = 1, 1, 1
returnStuff = mesh.getElementByCoordinates(x, y, z) # This will crash
gmsh.fltk.run()
gmsh.finalize()
Edited by Jonathan Stollberg