Coverage for pyxfem/postpro/gmsh_post_pro_interface.py: 49%
215 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3# pylint: disable=C0103
4"""
5Interface with gmsh-postpro
6Created on Thu Nov 23 14:05:34 2023
8@author: nchevaug
9"""
10from __future__ import annotations
11import gmsh
12import numpy as np
13import numpy.typing as npt
15FNDArray = npt.NDArray[np.float64]
17nbData2DataType = {1:'scal', 2:'vec2d', 3:'vec3d', 4:'tens2d', 9:'tens3d' }
18entityType2nbVert = {'point':1,'edge':2, 'tri':3, 'quad':4, 'tet':4}
19entityType2minDim = {'point':1,'edge':1, 'tri':2, 'quad':2, 'tet':3}
20entityType2gmsh = {'point':'P','edge':'L', 'tri':'T'} #, 'quad':2, 'tet':3}
21dataType2gmsh = {'scal':'S', 'vec2d':'V', 'vec3d':'V', 'tens2d': 'T','tens3d':'T' }
22gmshVectorDisplay = {'Line':1, 'Arrow':2, 'Pyramid':3, '3d Arrow':4, 'Displacement':5, 'Comet':6}
23gmshIntervalType = {'Iso-Value':1, 'Continuous map':2, 'Filled iso-values':3, 'Numeric values':4}
24gmshPointType = {'Color dot':4, '3D sphere':1, 'Scaled dot':2, 'Scaled sphere':3}
25gmshLineType = {'Color segment':3, '3D cylinder':1, 'Tapered cylinder':2}
28#views = []
31def gmshExportsViews(basefilename, extension, views = None) :
32 ''' export view; either as mpg or .pos '''
33 if views is None :
34 views = []
35 if extension in ['mpg']:
36 gmsh.write(basefilename+'.'+extension)
37 return
38 if extension == 'pos':
39 if len(views)> 0:
40 gmsh.view.write(views[0], basefilename+'.pos', append = False)
41 for view in views[1:]:
42 gmsh.view.write(view, basefilename+'.pos', append = True)
43 gmsh.write(basefilename+'.opt')
44 return
45 for view in views:
46 gmsh.view.setVisibilityPerWindow(view, 0)
47 gmsh.option.setNumber('General.SmallAxes' , 0)
48 digits = 1 if len(views) == 0 else int(np.log10(len(views))+1)
49 for i, view in enumerate(views) :
50 gmsh.view.setVisibilityPerWindow(view, 1)
51 gmsh.write(basefilename+str(i).zfill(digits)+'.'+extension)
52 gmsh.view.setVisibilityPerWindow(view, 0)
53 for view in views:
54 gmsh.view.setVisibilityPerWindow(view, 1)
55 gmsh.option.setNumber('General.SmallAxes' , 1)
57def gmshListPlotFormatEntities(xyz_:FNDArray, entity_type:str, nEntities = None) :
58 ''' set xyz and nEntities to correct format '''
59 in_shapexyz_ = xyz_.shape
60 fname = gmshListPlotFormatEntities.__name__
61 valueerror_str="Can't interprete the xyz in "+fname
62 nbv = entityType2nbVert.get(entity_type)
63 if nbv is None :
64 raise ValueError("entityType "+entity_type+"is unknown in "+fname)
65 if nEntities is None :
66 if len(xyz_.shape) == 1:
67 raise ValueError(valueerror_str)
68 nEntities = xyz_.shape[0]
69 if len(xyz_.shape) == 1 :
70 if xyz_.shape[0]%nEntities :
71 raise ValueError(valueerror_str)
72 xyz_ = xyz_.reshape((nEntities, -1))
73 if len(xyz_.shape) == 2 :
74 min_dim = entityType2minDim[entity_type]
75 if xyz_.shape[1] == nbv :
76 if min_dim > 1 :
77 raise ValueError(valueerror_str)
78 xyz_ = xyz_.reshape((nEntities, nbv, 1))
79 elif xyz_.shape[1] == 2*nbv :
80 if min_dim > 2 :
81 raise ValueError(valueerror_str)
82 xyz_ = xyz_.reshape((nEntities, nbv, 2))
83 elif xyz_.shape[1] == 3*nbv:
84 xyz_= xyz_.reshape((nEntities, nbv, 3))
85 else : raise ValueError(valueerror_str)
86 if xyz_.shape[:2] != (nEntities, nbv) :
87 raise ValueError(valueerror_str)
88 if len(xyz_.shape) != 3 :
89 raise ValueError("xyz should be of ndim atleast 3")
90 ncoord = xyz_.shape[2]
91 if ncoord < 3 :
92 #filling up yz coordinates with zeros if needed
93 z = np.zeros(xyz_.shape[:-1]+(3-ncoord,))
94 xyz = np.concatenate((xyz_,z),2)
95 elif ncoord == 3 :
96 xyz = xyz_.copy()
97 else :
98 raise ValueError(valueerror_str)
99 xyz = np.moveaxis(xyz, 1,2)
100 # return shapes to there input values
101 xyz_ = xyz_.reshape(in_shapexyz_)
102 # return result
103 return nEntities, xyz
105def gmshListPlotFormatDataElem2ElemNode(constant_elem_data:FNDArray, nEntities, entityType) \
106 -> FNDArray:
107 ''' Transform data per element to data per element per node, as needed by gmsh list format '''
108 nbv = entityType2nbVert[entityType]
109 data = np.repeat(constant_elem_data.reshape((nEntities, -1)), nbv, axis = 0)
110 return data.reshape((nEntities, nbv, -1))
112def gmshListPlotFormatData(data_:FNDArray, nEntities:int, entityType:str, P0 =False) \
113 -> tuple[str, int, FNDArray] :
114 fname = gmshListPlotFormatData.__name__
115 errorstr = "can't interprete the data_ in "+fname
116 inShapedata_= data_.shape
117 indata_ = data_
118 if P0:
119 data_ = gmshListPlotFormatDataElem2ElemNode(data_, nEntities, entityType)
120 data_type = 'unknown'
121 nbv = entityType2nbVert[entityType]
122 if data_.ndim > 4:
123 raise ValueError(errorstr)
124 if data_.ndim == 1 :
125 if data_.size%nEntities == 0 :
126 data_ =data_.reshape((nEntities, -1))
127 else :
128 raise ValueError(errorstr)
129 if data_.shape[0] != nEntities :
130 raise ValueError(errorstr)
131 if data_.ndim == 2 :
132 if data_.shape[1]%nbv != 0:
133 raise ValueError(errorstr)
134 data_ = data_.reshape((nEntities, nbv, -1))
135 if data_.ndim == 3 :
136 nbdata = data_.shape[2]
137 data_type = nbData2DataType.get(nbdata, None)
138 if data_type is None:
139 raise ValueError(errorstr)
140 if data_type == 'tens2d' :
141 data_ = data_.reshape((nEntities, nbv, 2, 2))
142 if data_type == 'tens3d' :
143 data_ = data_.reshape((nEntities, nbv, 3, 3))
144 elif data_.ndim == 4:
145 if data_.shape[2:] == (2,2) :
146 data_type = 'tens2d'
147 elif data_.shape[2:] == (3,3) :
148 data_type = 'tens3d'
149 else :
150 raise ValueError(errorstr)
151 # filling up data with zeros where needed
152 if data_type == 'vec2d':
153 data = np.zeros((nEntities, nbv, 3))
154 data[:,:,0:-1] = data_
155 elif data_type == 'tens2d':
156 data = np.zeros((nEntities, nbv, 3, 3))
157 data[:,:,0:-1, 0:-1] = data_
158 else :
159 data = data_.copy()
160 # return shapes to there input values
161 indata_ = indata_.reshape(inShapedata_)
162 # return result
163 return data_type, nEntities, data
165def gmshListPlotFormatEntityData(entityVertCoo, data, entityType, nEntities = None, P0 = False):
166 nEntities, entityVertCoo = gmshListPlotFormatEntities(entityVertCoo, entityType, nEntities)
167 dataType, __, data = gmshListPlotFormatData(data, nEntities, entityType, P0)
168 data = np.column_stack((entityVertCoo.reshape((nEntities,-1)), data.reshape((nEntities,-1)))).flatten()
169 gmshDataTypeEntityType = dataType2gmsh[dataType] + entityType2gmsh[entityType]
170 return gmshDataTypeEntityType, nEntities, data
172def setViewOption(gv, IntervalsType = "Filled iso-values", NbIso = 20, VectorType = "Arrow", DisplacementFactor = 1., Range = None) :
173 gmsh.view.option.setNumber(gv, "IntervalsType", gmshIntervalType.get(IntervalsType, 3))
174 gmsh.view.option.setNumber(gv, "NbIso", NbIso)
175 gmsh.view.option.setNumber(gv, "VectorType", gmshVectorDisplay.get(VectorType, 2))
176 gmsh.view.option.setNumber(gv, "DisplacementFactor", DisplacementFactor)
177 if Range is not None :
178 gmsh.view.option.setNumber(gv, "RangeType", 2)
179 gmsh.view.option.setNumber(gv, "CustomMin", Range[0])
180 gmsh.view.option.setNumber(gv, "CustomMax", Range[1])
181 return gv
183def listPlotFieldPoint(xyz, data, viewname = 'T', nPoints = None, PointType = "3D sphere", PointSize = 5,
184 IntervalsType = "Filled iso-values", NbIso = 20, VectorType = "Arrow", DisplacementFactor = 0.,
185 Range = None, gv = None):
186 if gv is None :
187 gv = gmsh.view.add(viewname)
188 #views.append(gv)
189 gmsh.view.addListData(gv, *gmshListPlotFormatEntityData(xyz, data, 'point', nPoints) )
190 gmsh.view.option.setNumber(gv,"PointType", gmshPointType.get(PointType, 2))
191 gmsh.view.option.setNumber(gv,"PointSize", PointSize)
192 setViewOption(gv, IntervalsType = IntervalsType, NbIso = NbIso,
193 VectorType = VectorType, DisplacementFactor = DisplacementFactor, Range = Range)
194 return gv
196def listPlotFieldLine(lineVerticesCoord, data, viewname = 'T', P0 = False, nLines = None,
197 LineType = '3D cylinder', LineWidth = 2, IntervalsType = "Filled iso-values",
198 NbIso = 20, VectorType = "Arrow", DisplacementFactor = 0., Range = None, gv = None):
199 if gv is None :
200 gv = gmsh.view.add(viewname)
201 #views.append(gv)
202 gmsh.view.addListData(gv, *gmshListPlotFormatEntityData(lineVerticesCoord, data, 'edge',
203 nLines, P0))
204 gmsh.view.option.setNumber(gv, "LineType", gmshLineType.get(LineType, 2))
205 gmsh.view.option.setNumber(gv, "LineWidth", LineWidth)
206 gmsh.view.option.setNumber(gv, "GlyphLocation", 1 if P0 else 2)
207 setViewOption(gv, IntervalsType = IntervalsType, NbIso = NbIso, VectorType = VectorType,
208 DisplacementFactor = DisplacementFactor, Range = Range)
209 return gv
211def listPlotFieldTri(trisVerticesCoord, trisData, viewname ='T', P0 = False, ntri = None,
212 IntervalsType = "Filled iso-values", NbIso = 20, VectorType = "Arrow",
213 DisplacementFactor = 0., Range = None, gv =None):
214 if gv is None :
215 gv = gmsh.view.add(viewname)
216 #views.append(gv)
217 gmsh_data_entity_type, ntri, data = gmshListPlotFormatEntityData(trisVerticesCoord,
218 trisData, 'tri', ntri, P0)
219 gmsh.view.addListData(gv, gmsh_data_entity_type, ntri, data)
220 gmsh.view.option.setNumber(gv,"ShowElement", 0 if ntri > 10000 else 1)
221 gmsh.view.option.setNumber(gv, "GlyphLocation", 1 if P0 else 2)
222 setViewOption(gv, IntervalsType = IntervalsType, NbIso = NbIso,
223 VectorType = VectorType, DisplacementFactor = DisplacementFactor, Range = Range)
224 return gv
226if __name__ == '__main__':
227 testTri = True
228 testLine = True
229 testPoint = True
230 popUpGmsh = True
231 gmsh.initialize()
232 vertices_xy = np.array([[0.,0.], [1.,0.], [1.,1.], [0.,1.]])
233 vertices_xyz = np.array([[0.,0.,0.], [1.,0.,0.], [1.,1.,1.], [0.,1.,0.]])
234 if testPoint :
235 vertex_data = np.array([0.,1.,2.,3.])
236 listPlotFieldPoint(vertices_xy, vertex_data, viewname = 'Point xy Scalar',
237 nPoints = None, PointType = '3D sphere', PointSize = 5,
238 IntervalsType = "Filled iso-values", VectorType = "Arrow", gv = None)
239 listPlotFieldPoint(vertices_xyz, vertex_data, viewname = 'Point xyz, Scalar',
240 nPoints = None, PointType = 'Color dot', PointSize = 5,
241 IntervalsType = "Filled iso-values", VectorType = "Arrow", gv = None)
242 vertex_data = np.array([1.,0.,0.,1., -1.,0., 0., -1.])
243 listPlotFieldPoint(vertices_xy, vertex_data, viewname = 'Point xy vector2d',
244 nPoints = None, PointType = 'Scaled dot', PointSize = 5,
245 IntervalsType = "Filled iso-values", VectorType = "Arrow", gv = None)
246 listPlotFieldPoint(vertices_xyz, vertex_data, viewname = 'Point xyz, vector2d',
247 nPoints = None, PointSize = 5, IntervalsType = "Filled iso-values",
248 VectorType = "Arrow", gv = None)
249 vertex_data = np.array([1.,0.,0., 0.,1.,1., -1.,0., -1., 0., -1., 0.])
250 listPlotFieldPoint(vertices_xy, vertex_data, viewname = 'Point xy vector3d',
251 nPoints = None, PointSize = 5, IntervalsType = "Filled iso-values",
252 VectorType = "Arrow", gv = None)
253 listPlotFieldPoint(vertices_xyz, vertex_data, viewname = 'Point xyz, vector3d',
254 nPoints = None, PointSize = 5, IntervalsType = "Filled iso-values",
255 VectorType = "Comet", gv = None)
256 if testLine :
257 lines = np.array([[0,1],[1,2],[2,3], [3,0]])
258 xyLine = vertices_xy[lines]
259 xyzLine = vertices_xy[lines]
260 scalDataP1 = np.array([[0,1],[1,2], [2,3], [3,4]])
261 scalDataP0 = np.array([0, 1, 2, 3])
262 vectDataP0 = np.array([[0,-1],[1,0], [0,1], [-1,0]])
263 listPlotFieldLine(xyzLine, scalDataP1, viewname = 'Line xyz scalar P1')
264 listPlotFieldLine(xyzLine, scalDataP0, P0 =True,
265 viewname = 'Line xyz scalar P0', LineType="Tapered cylinder")
266 listPlotFieldLine(xyzLine, vectDataP0, P0 =True,
267 viewname = 'Line xyz vector P0', LineType="Tapered cylinder")
268 if testTri :
269 triCoord = np.array([[[0.,0.],[1.,0.], [0., 1.]],[[1.,0.], [1.,1.], [0.,1.]]])
270 #test tri scalar field
271 triScalData = np.array([[0.,1.,2], [1.,3.,2]])
272 listPlotFieldTri(triCoord, triScalData, viewname ='scalar Tri test', NbIso = 10,
273 Range = None, gv =None)
274 #test tri vector 2d field
275 triVec2dData = np.array([[[0.,0.], [1.,0.], [0.,1]], [[1.,0.], [1.,1.], [0.,1.]]])
276 listPlotFieldTri(triCoord, triVec2dData, viewname ='vecteur2D Tri test', NbIso = 10,
277 Range = None, gv =None)
278 #test tri vector 3d field
279 triVec3dData = np.array([[[0.,0.,1.], [1.,0.,1.], [0.,1., 1.]],
280 [[1.,0.,1.], [1.,1.,1.], [0.,1.,1.]]])
281 listPlotFieldTri(triCoord, triVec3dData, viewname ='vecteur3D Tri test', NbIso = 10,
282 Range = None, gv =None)
283 #test tri tensor 2d Flatten field
284 triTens2dData = np.array([[[1.,0., 0., 0.], [0.,1.,0.,0.], [0.,0.,1.,0.]],
285 [[0.,1.,0.,0.], [0.,0.,0.,1.], [0.,0.,1.,0.]]])
286 listPlotFieldTri(triCoord, triTens2dData, viewname ='tensor2D flat Tri test', NbIso = 10,
287 Range = None, gv =None)
288 #test tri tensor 2d Flatten field
289 triTens2dData = np.array([[[[1.,0.], [0., 0.]], [[0.,1.],[0.,0.]], [[0.,0.],[1.,0.]]],
290 [[[0.,1.],[0.,0.]], [[0.,0.],[0.,1.]], [[0.,0.],[1.,0.]]]])
291 listPlotFieldTri(triCoord, triTens2dData, viewname ='tensor2D Tri test', NbIso = 10,
292 Range = None, gv =None)
293 gmsh.option.setNumber('PostProcessing.AnimationCycle', 1)
294 if popUpGmsh:
295 gmsh.fltk.run() #start gmsh gui
296 gmsh.finalize() #close gmsh