Coverage for pyxfem/core/gmsh_mesh_model_interface.py: 58%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on Tue Feb 20 10:14:29 2024 

5@author: nchevaug 

6This module is an interface to gmsh. gmsh must be initialized for the function to work. 

7""" 

8from __future__ import annotations 

9from typing import Literal 

10import gmsh 

11import numpy as np 

12from pyxfem.core.config import IndexType 

13if __name__ == "__main__": 

14 import pyxfem.core.gmsh_mesh_generator as meshbuilder 

15#import numpy.typing as npt 

16 

17gmshTypeTag2DimTypeNbv = {15:(0, 'point', 1), 1:(1, 'edge', 2), 

18 2:(2, 'tri', 3), 3:(2, 'quad', 4), 4:(3, 'tet', 4)} 

19""" a dictionnary to translate gmsh element type tag to a triplet (dim, elementtypename, nbvertex) 

20othergmsh elem type can be or in gmsh documentation or there : 

21https://docs.juliahub.com/GmshTools/9rYp5/0.4.2/element_types/ """ 

22 

23def set_model(func): 

24 ''' a decorator for function (func) that query gmsh 

25 

26 the first argument of these func is supposed to be the modelName 

27 the warper set gmsh current model to modelName, then call func 

28 then reset the model name to the previous one. 

29 ''' 

30 def wrapper_set_model(*args, **kwargs): 

31 gmsh_model_name = args[0] 

32 old_gmsh_model_name = gmsh.model.getCurrent() 

33 gmsh.model.setCurrent(gmsh_model_name) 

34 r = func(*args, **kwargs) 

35 gmsh.model.setCurrent(old_gmsh_model_name) 

36 return r 

37 return wrapper_set_model 

38 

39@set_model 

40def nb_vertices(gmsh_model_name): # pylint: disable=W0613 

41 ''' Return the number of vertices in gmsh_model_name ''' 

42 return len(gmsh.model.mesh.getNodes()[0]) 

43 

44@set_model 

45def nb_triangles(gmsh_model_name): # pylint: disable=W0613 

46 ''' Return the number of Linear Triangle in gmsh_model_name ''' 

47 return len(gmsh.model.mesh.getElementsByType(2)[0]) 

48 

49@set_model 

50def get_vertices_coord(gmsh_model_name:str, dim:Literal[2,3]=2): # pylint: disable=W0613 

51 ''' Return all the vertices coordinate (xy, 2d problems) in mesh 

52 gmsh_model_name as an np.array of shape (nvertex, 2) ''' 

53 if dim not in [2, 3]: 

54 raise ValueError("print dim should be Literal 1 or 2") 

55 tag, xyz, _ = gmsh.model.mesh.getNodes() 

56 xyz = np.array(xyz) 

57 if dim == 2 : 

58 xy_ = xyz.reshape((-1, 3))[:, :-1] 

59 return xy_[np.argsort(tag)] 

60 if dim == 3: 

61 xyz_ = xyz.reshape((-1, 3)) 

62 return xyz_[np.argsort(tag)] 

63 

64@set_model 

65def get_triangles(gmsh_model_name): # pylint: disable=W0613 

66 ''' Return all the linear triangles connectivities in mesh gmsh_model_name 

67 as an np.array of int of shape (ntris, 3), 

68 and array with the corresponding gmshtag of the triangles ''' 

69 gmshtags, tris2v = gmsh.model.mesh.getElementsByType(2) 

70 tris2v = np.array(tris2v).reshape((-1, 3)) -1 

71 return tris2v, np.array(gmshtags, dtype = IndexType)-1 

72 

73@set_model 

74def get_vertices_on_phys_tag(gmsh_model_name, dim, phystag): # pylint: disable=W0613 

75 ''' return the nodes id and coordinates of all the nodes of model gmsh_model_name 

76 classifyied on physical entity dim, phystag ''' 

77 nodes_tag, coord = gmsh.model.mesh.getNodesForPhysicalGroup(dim, phystag) 

78 nodes_tag = np.array(nodes_tag) 

79 coord = np.array(coord) 

80 return nodes_tag-1, coord.reshape((-1, 3))[:, :-1] 

81 

82@set_model 

83def get_toplevel_elements_on_phys_tag(gmsh_model_name, dim, phystag): 

84 ''' return all the element of highest possible dim on geometrical entity (dim ,phystag) ''' 

85 tags = gmsh.model.getEntitiesForPhysicalGroup(dim, phystag) 

86 etype2ids = {} 

87 for tag in tags: 

88 elem_types, elem_ids, _ = gmsh.model.mesh.getElements(dim, tag) 

89 elem_types = np.array(elem_types) 

90 elem_ids = np.array(elem_ids) 

91 for typetag, elem_id in zip(elem_types, elem_ids): 

92 edim, name, __ = gmshTypeTag2DimTypeNbv[typetag] 

93 if edim == dim: 

94 etype2ids[name] = np.append( 

95 etype2ids.get(name, np.zeros(0, dtype=np.uint64)), 

96 elem_id -1) 

97 return etype2ids 

98 

99@set_model 

100def get_toplevel_elements_vertices_on_phys_tag(gmsh_model_name, dim, phystag): # pylint: disable=W0613 

101 ''' get all the element of gmsh_model_name of dimension dim classifyed on dim phystag 

102 the result is given as dictionnary, indexed by elementype name 

103 ('point', 'edge', 'tri', 'quad', 'tet' ... etc) 

104 for each entry of the dictionnary, an np array of shape (n, nv) 

105 containing for each of the n elements there nv node. 

106 ''' 

107 tags = gmsh.model.getEntitiesForPhysicalGroup(dim, phystag) 

108 etype2nodes = {} 

109 for tag in tags: 

110 elem_types, ___, elem_node_tags = gmsh.model.mesh.getElements(dim, tag) 

111 elem_types = np.array(elem_types) 

112 elem_node_tags = np.array(elem_node_tags) 

113 for typetag, elem_nodes in zip(elem_types, elem_node_tags): 

114 edim, name, nv = gmshTypeTag2DimTypeNbv[typetag] 

115 if edim == dim: 

116 etype2nodes[name] = np.vstack(( 

117 etype2nodes.get(name, np.zeros((0, nv), dtype=np.uint64)), 

118 elem_nodes.reshape((-1, nv)) -1)) 

119 return etype2nodes 

120 

121@set_model 

122def get_group_name2dim_phys_tag(gmsh_model_name, groupe_names=None): # pylint: disable=W0613 

123 """ return a dictionary mapping group Name to dim and phystag """ 

124 dic = {gmsh.model.getPhysicalName(*pg): pg for pg in gmsh.model.getPhysicalGroups()} 

125 if groupe_names is None: 

126 return dic 

127 return {name: dic[name] for name in groupe_names} 

128 

129@set_model 

130def renumber(gmsh_model_name, method = 'RCMK'): 

131 """ renumber the msj nod to improve memory bandwith 

132 method must be one off: 

133 - 'RCMK' Reverse Cuthill McKee 

134 - 'Hilbert' : space filling hilbert curves 

135 - 'Metis" : well metis ... 

136 """ 

137 old, new = gmsh.model.mesh.computeRenumbering(method) 

138 gmsh.model.mesh.renumberNodes(old, new) 

139 

140if __name__ == '__main__': 

141 def main(): 

142 ''' a simple test checking how set the set_model decorator works''' 

143 gmsh.initialize() 

144 gmsh_model_name1, phys_dict1 = meshbuilder.square(relh=2.) 

145 gmsh_model_name2, phys_dict2 = meshbuilder.square(relh=1.) 

146 for gmsh_model_name, phys_dict in zip([gmsh_model_name1, gmsh_model_name2], 

147 [phys_dict1, phys_dict2]): 

148 print("### model: ", gmsh_model_name) 

149 print(f'nv: {nb_vertices(gmsh_model_name):d},\ 

150 nf: {nb_triangles(gmsh_model_name):d} ') 

151 print('vertices: ') 

152 print(get_vertices_coord(gmsh_model_name)) 

153 print('Triangles to vertices: ') 

154 print(get_triangles(gmsh_model_name)[0]) 

155 print('vertices On bottom physical edge') 

156 bottom_vertex, bottom_xy = get_vertices_on_phys_tag(gmsh_model_name, 

157 *phys_dict['bottom']) 

158 print(bottom_vertex) 

159 print(bottom_xy) 

160 print('edges On bottom physical edge') 

161 e2v = get_toplevel_elements_vertices_on_phys_tag(gmsh_model_name, 

162 *phys_dict['bottom'])['edge'] 

163 print(e2v) 

164 print("Dictionnary: GroupeName To dim phystag ") 

165 print(get_group_name2dim_phys_tag(gmsh_model_name)) 

166 gmsh.finalize() 

167 

168 main()