...
Code Block | ||
---|---|---|
| ||
proc MeshIo::ReadPreUnstructuredMeshInit { filename } { set fail 0variable meshio_num_nodes_per_cell variable meshio_gid_element #linear and serendipit variableset meshio_num_nodes_per_cell(vertex) 1 MeshIo::Import_gid_meshio_py ;#to load in python the file gid_meshio.py to define its python functions before be called set meshio_num_nodes_per_cell(line) 2 set meshio_num_nodes_per_cell(triangle) 3 set meshio_num_nodes_per_cell(quad) 4 set m [GiD_Python_Call gid_meshio.my_meshio_read_mesh $filename]meshio_num_nodes_per_cell(quad8) 8 set meshio_num_nodes_coordinates [lindex [MeshIo::PythonArrayToTclList [lindex $m 0]] 0]per_cell(tetra) 4 set meshio_num_nodes_per_cell(hexahedron) 8 #e.g.[lindex $m 1] == {<meshio CellBlock, type: triangle, num cells: 156, tags: []>}set meshio_num_nodes_per_cell(hexahedron20) 20 set meshio_num_nodes_per_cell(hexahedron24) 24 set meshio_elementnum_typesnodes_and_connectivities [MeshIo::PythonArrayToTclList [lindex $m 2]] set layer [GiD_Layers get to_use]per_cell(wedge) 6 set meshio_num_nodes_per_cell(pyramid) 5 #quadratic set offsetmeshio_num_nodes [GiD_Info mesh MaxNumNodes]_per_cell(line3) 3 set offset_elements [GiD_Info mesh MaxNumElements]meshio_num_nodes_per_cell(triangle6) 6 set lastmeshio_num_elementnodes_id $offset_elements per_cell(quad9) 9 #better use GiD_MeshPre_Create with same syintax as GiD_MeshPost (and some day could be implemented to be faster in C++)set meshio_num_nodes_per_cell(tetra10) 10 set meshio_num_nodes_per_cell(hexahedron27) 27 set meshio_num_nodes_per_cell(wedge15) 15 set meshio_num_nodes [llength $nodes_coordinates]_per_cell(wedge18) 18 set node_ids [objarray new_from_to intarray [expr $offset_nodes+1] [expr $offset_nodes+$num_nodes]] meshio_num_nodes_per_cell(pyramid13) 13 set meshio_num_nodes_per_cell(pyramid14) 14 set vertices [objarray new doublearray [expr $num_nodes*3]] #degree 3 set i 0meshio_num_nodes_per_cell(line4) 4 foreach node $nodes_coordinates { set meshio_num_nodes_per_cell(triangle10) 10 set meshio_num_nodes_per_cell(quad16) 16 foreach value $node { set meshio_num_nodes_per_cell(tetra20) 20 objarray set $vertices $i $valuemeshio_num_nodes_per_cell(wedge40) 40 set meshio_num_nodes_per_cell(hexahedron64) 64 #degree 4 incr i }set meshio_num_nodes_per_cell(line5) 5 set }meshio_num_nodes_per_cell(triangle15) 15 foreachset {meshio_elementnum_type meshio_connectivities} $meshio_element_types_and_connectivities { nodes_per_cell(quad25) 25 set element_type "" meshio_num_nodes_per_cell(tetra35) 35 set elementmeshio_num_nodes 0_per_cell(wedge75) 75 if { [info existsset meshio_gidnum_element($meshionodes_elementper_typecell(hexahedron125)] }125 { #degree 5 set element_type $meshio_gid_element($meshio_element_type) meshio_num_nodes_per_cell(line6) 6 set meshio_num_nodes_per_cell(triangle21) 21 set elementmeshio_num_nodes $meshio_num_nodes_per_cell($meshio_element_typequad36) 36 } else {set meshio_num_nodes_per_cell(tetra56) 56 set meshio_num_nodes_per_cell(wedge126) 126 W "element $meshio_element_type not supported"set meshio_num_nodes_per_cell(hexahedron216) 216 #degree 6 set meshio_num_nodes_per_cell(line7) 7 continue set meshio_num_nodes_per_cell(triangle28) 28 } set meshio_num_nodes_per_cell(quad49) 49 set elements [lindex $meshio_connectivities 0] meshio_num_nodes_per_cell(tetra84) 84 set meshio_num_nodes_per_cell(wedge196) 196 set meshio_num_elements [llength $elements]nodes_per_cell(hexahedron343) 343 #degree 7 set element_ids [objarray new_from_to intarray [expr $last_element_id+1] [expr $last_element_id+$num_elements]] meshio_num_nodes_per_cell(line8) 8 set meshio_num_nodes_per_cell(triangle36) 36 set element_vertex_indices [objarray new intarray [expr $num_elements*$element_num_nodes]] meshio_num_nodes_per_cell(quad64) 64 set meshio_num_nodes_per_cell(tetra120) 120 set meshio_num_nodes_per_cell(wedge288) 288 set i 0 set meshio_num_nodes_per_cell(hexahedron512) 512 foreach#degree element8 $elements { set meshio_num_nodes_per_cell(line9) 9 foreach node_id $element {set meshio_num_nodes_per_cell(triangle45) 45 set meshio_num_nodes_per_cell(quad81) 81 set meshio_num_nodes_per_cell(tetra165) 165 objarray set $elementmeshio_vertex_indices $i [expr $node_id+$offset_nodes+1] num_nodes_per_cell(wedge405) 405 set meshio_num_nodes_per_cell(hexahedron729) 729 #degree 9 incr i set meshio_num_nodes_per_cell(line10) 10 set meshio_num_nodes_per_cell(triangle55) 55 } set meshio_num_nodes_per_cell(quad100) 100 } set meshio_num_nodes_per_cell(tetra220) 220 set zero_based_array 0meshio_num_nodes_per_cell(wedge550) 550 set meshio_num_nodes_per_cell(hexahedron1000) 1000 GiD_MeshPre_Create $element_type $element set meshio_num_nodes $node_ids $vertices $element_ids $element_vertex_indices $zero_based_array $layer_per_cell(hexahedron1331) 1331 #degree 10 incrset last_element_id $num_elementsmeshio_num_nodes_per_cell(line11) 11 set meshio_num_nodes_per_cell(triangle66) 66 } set meshio_num_nodes_per_cell(quad121) 121 return $fail } |
For the import feature GiD is invoking from its Tcl code functions of Python files and modules like meshio and numpy.
For the export feature it is implemented in two alternative ways:
For Tcl-like programmers, with most code in Tcl
For Pyhton-like programmers: with most code in Python
Note that only the second approach is bi-directional and will invoke from Python Tcl commands (of the GiD interpreter) importing in Python the tohil module
Ask to GiD its current mesh data with Tcl, and process this data with Tcl code to reach the format expected by the Python meshio function and then call a Python function that create a meshio.Mesh and write it to file with the desired format.
From Tcl call some Python code that ask GiD mesh information calling Tcl commands and process this data with Python to create a meshio.Mesh and write it to file with the desired format.
The implementation of 1. is something like this
meshio.tcl
(big code)
Code Block | ||
---|---|---|
| ||
proc MeshIo::Init { } { variable meshio_num_nodes_per_cell variable meshio_gid_element variable meshio_element_name variable meshio_element_num_nodes #linear and serendipit set meshio_num_nodes_per_cell(vertex) 1set meshio_num_nodes_per_cell(tetra286) 286 #linear and serendipit set meshio_gid_element(vertex) point set meshio_gid_element(line) line set meshio_gid_element(triangle) triangle set meshio_gid_element(quad) quadrilateral set meshio_gid_element(tetra) tetrahedra set meshio_gid_element(hexahedron) hexahedra set meshio_gid_element(wedge) prism set meshio_gid_element(pyramid) pyramid #quadratic set meshio_gid_element(line3) line set meshio_num_nodes_per_cell(linegid_element(triangle6) 2triangle set meshio_num_nodes_per_cell(trianglegid_element(quad8) 3quadrilateral set meshio_num_nodes_per_cell(quadgid_element(quad9) 4quadrilateral set meshio_num_nodes_per_cell(quad8gid_element(tetra10) 8tetrahedra set meshio_num_nodes_per_cell(tetragid_element(hexahedron20) 4hexahedra set meshio_num_nodes_per_cell(hexahedrongid_element(hexahedron27) 8hexahedra set meshio_num_nodes_per_cell(hexahedron20gid_element(wedge15) 20prism set meshio_num_nodes_per_cell(hexahedron24gid_element(pyramid13) 24pyramid set meshio_num_nodes_per_cell(wedge) 6 set meshio_num_nodes_per_cell(pyramid) 5 #quadratic ... } proc MeshIo::ReadPreUnstructuredMesh { filename } { set fail 0 setvariable meshio_num_nodes_per_cell(line3) 3gid_element setvariable meshio_num_nodes_per_cell(triangle6) 6 set meshio_num_nodes_per_cell(quad9) 9 set meshio_num_nodes_per_cell(tetra10) 10 set meshio_num_nodes_per_cell(hexahedron27) 27MeshIo::Import_gid_meshio_py ;#to load in python the file gid_meshio.py to define its python functions before be called set m [GiD_Python_Call gid_meshio.my_nummeshio_nodes_per_cell(wedge15) 15read_mesh $filename] set meshio_num_nodes_per_cell(wedge18) 18 set meshio_num_nodes_per_cell(pyramid13) 13 set meshio_num_nodes_per_cell(pyramid14) 14 #degree 3 set meshio_num_nodes_per_cell(line4) 4coordinates [lindex [MeshIo::PythonArrayToTclList [lindex $m 0]] 0] #e.g.[lindex $m 1] == {<meshio CellBlock, type: triangle, num cells: 156, tags: []>} set meshio_numelement_nodestypes_per_cell(triangle10) 10and_connectivities [MeshIo::PythonArrayToTclList [lindex $m 2]] set meshio_num_nodes_per_cell(quad16) 16layer [GiD_Layers get to_use] set meshio_numoffset_nodes_per_cell(tetra20) 20 [GiD_Info mesh MaxNumNodes] set meshio_num_nodes_per_cell(wedge40) 40offset_elements [GiD_Info mesh MaxNumElements] set meshiolast_numelement_nodes_per_cell(hexahedron64) 64id $offset_elements #degree 4 #better setuse meshioGiD_num_nodes_per_cell(line5) 5 set meshio_num_nodes_per_cell(triangle15) 15 set meshio_num_nodes_per_cell(quad25) 25MeshPre_Create with same syintax as GiD_MeshPost (and some day could be implemented to be faster in C++) set meshio_num_nodes_per_cell(tetra35) 35 [llength $nodes_coordinates] set meshio_num_nodes_per_cell(wedge75) 75 set meshio_num_nodes_per_cell(hexahedron125) 125node_ids [objarray new_from_to intarray [expr $offset_nodes+1] [expr $offset_nodes+$num_nodes]] set #degreevertices 5[objarray new doublearray set meshio_num_nodes_per_cell(line6) 6[expr $num_nodes*3]] set meshio_num_nodes_per_cell(triangle21) 21 i 0 foreach setnode meshio_num_nodes_per_cell(quad36) 36 set meshio_num_nodes_per_cell(tetra56) 56$nodes_coordinates { foreach value $node { set meshio_num_nodes_per_cell(wedge126) 126 objarray set meshio_num_nodes_per_cell(hexahedron216) 216 $vertices $i $value #degree 6 set meshio_num_nodes_per_cell(line7) 7 incr i set meshio_num_nodes_per_cell(triangle28) 28 } set meshio_num_nodes_per_cell(quad49) 49 } setforeach {meshio_element_type meshio_connectivities} $meshio_numelement_nodestypes_per_cell(tetra84) 84and_connectivities { set meshio_num_nodes_per_cell(wedge196) 196element_type "" set meshioelement_num_nodes_per_cell(hexahedron343) 343 0 #degree 7 if { [info setexists meshio_numgid_nodeselement($meshio_perelement_cell(line8) 8type)] } { set meshio_num_nodes_per_cell(triangle36) 36 set meshio_num_nodes_per_cell(quad64) 64element_type $meshio_gid_element($meshio_element_type) set meshio_num_nodes_per_cell(tetra120) 120 set meshioelement_num_nodes $meshio_pernum_nodes_per_cell(wedge288) 288$meshio_element_type) } else set meshio_num_nodes_per_cell(hexahedron512) 512{ #degree 8 W set"element meshio$meshio_num_nodes_per_cell(line9) 9element_type not supported" set meshio_num_nodes_per_cell(triangle45) 45 set meshio_num_nodes_per_cell(quad81) 81 continue set meshio_num_nodes_per_cell(tetra165) 165 } set meshio_num_nodes_per_cell(wedge405) 405 set elements set[lindex meshio_num_nodes_per_cell(hexahedron729) 729$meshio_connectivities 0] #degree 9 set meshio_num_nodes_per_cell(line10) 10 set meshio_num_nodes_per_cell(triangle55) 55elements [llength $elements] set meshio_num_nodes_per_cell(quad100) 100element_ids [objarray new_from_to intarray [expr $last_element_id+1] [expr $last_element_id+$num_elements]] set meshioelement_num_nodes_per_cell(tetra220) 220 set meshio_num_nodes_per_cell(wedge550) 550vertex_indices [objarray new intarray [expr $num_elements*$element_num_nodes]] set meshio_num_nodes_per_cell(hexahedron1000) 1000 set meshio_num_nodes_per_cell(hexahedron1331) 1331i 0 foreach element $elements { #degree 10 setforeach meshio_num_nodes_per_cell(line11) 11 node_id $element { set meshio_num_nodes_per_cell(triangle66) 66 set meshio_num_nodes_per_cell(quad121) 121 objarray set meshio$element_num_nodes_per_cell(tetra286) 286vertex_indices $i [expr $node_id+$offset_nodes+1] #linear and serendipit set meshio_gid_element(vertex) point incr i set meshio_gid_element(line) line set meshio_gid_element(triangle) triangle } set meshio_gid_element(quad) quadrilateral } set meshio_gid_element(tetra) tetrahedra set meshiozero_gid_element(hexahedron) hexahedrabased_array 0 set meshioGiD_gid_element(wedge) prism set meshio_gid_element(pyramid) pyramidMeshPre_Create $element_type $element_num_nodes $node_ids $vertices $element_ids $element_vertex_indices $zero_based_array $layer #quadratic incr set meshiolast_gid_element(line3) lineelement_id $num_elements set meshio_gid_element(triangle6) triangle} return set meshio_gid_element(quad8) quadrilateral set meshio_gid_element(quad9) quadrilateral set meshio_gid_element(tetra10) tetrahedra set meshio_gid_element(hexahedron20) hexahedra set meshio_gid_element(hexahedron27) hexahedra set meshio_gid_element(wedge15) prism set meshio_gid_element(pyramid13) pyramid $fail } |
For the import feature GiD is invoking from its Tcl code functions of Python files and modules like meshio and numpy.
For the export feature it is implemented in two alternative ways:
For Tcl-like programmers, with most code in Tcl
For Pyhton-like programmers: with most code in Python
Note that only the second approach is bi-directional and will invoke from Python Tcl commands (of the GiD interpreter) importing in Python the tohil module
Ask to GiD its current mesh data with Tcl, and process this data with Tcl code to reach the format expected by the Python meshio function and then call a Python function that create a meshio.Mesh and write it to file with the desired format.
From Tcl call some Python code that ask GiD mesh information calling Tcl commands and process this data with Python to create a meshio.Mesh and write it to file with the desired format.
The implementation of 1. is something like this
meshio.tcl
(big code)
Code Block | ||
---|---|---|
| ||
proc MeshIo::Init { } {
variable meshio_element_name
...
set meshio_element_name(point,1) vertex
set meshio_element_name(line,2) line
set meshio_element_name(triangle,3) triangle
set meshio_element_name(quadrilateral,4) quad
set meshio_element_name(tetrahedra,4) tetra
set meshio_element_name(hexahedra,8) hexahedron
set meshio_element_name(prism,6) wedge
set meshio_element_name(pyramid,5) pyramid
#quadratic
set meshio_element_name(line,3) line3
set meshio_element_name(triangle,6) triangle6
set meshio_element_name(quadrilateral,8) quad8
set meshio_element_name(quadrilateral,9) quad9
set meshio_element_name(tetrahedra,10) tetra10
set meshio_element_name(hexahedra,20) hexahedron20
set meshio_element_name(hexahedra,27) hexahedron27
set meshio_element_name(prism,15) wedge15
set meshio_element_name(pyramid,13) pyramid13
}
proc MeshIo::TclObjarrayToPythonArrayPoints { node_xyzs } {
set points ""
lassign $node_xyzs xs ys zs
set num_nodes [objarray length $xs]
for {set i_node 0} {$i_node<$num_nodes} {incr i_node} {
set x [objarray get $xs $i_node]
set y [objarray get $ys $i_node]
set z [objarray get $zs $i_node]
append points "\[$x,$y,$z\],"
}
return "\[$points\]"
}
proc MeshIo::TclObjarrayToPythonArrayConnectivities { element_num_nodes connectivities } {
set cells ""
set num_elements [expr [objarray length $connectivities]/$element_num_nodes]
set i 0
for {set i_element 0} {$i_element<$num_elements} {incr i_element} {
set node_ids [list]
for {set i_node 0} {$i_node<$element_num_nodes} {incr i_node} {
lappend node_ids [objarray get $connectivities $i]
incr i
}
append cells "\[[join $node_ids ,]\],"
}
return "\[$cells\]"
}
proc MeshIo::WritePreUnstructuredMesh { filename } {
variable meshio_element_name
MeshIo::Import_gid_meshio_py ;#to load in python the file gid_meshio.py to define its python functions before be called
# coordinates
lassign [GiD_Info mesh nodes -array] node_ids node_xyzs
set max_id_nodes [objarray get $node_ids end]
set num_nodes [objarray length [lindex $node_xyzs 0]]
set new_node_ids ""
set nodes_renumbered_for_meshio 0
if { $max_id_nodes != $num_nodes } {
set nodes_renumbered_for_meshio 1
set new_node_ids [objarray new_from_to intarray 0 [expr $num_nodes-1]]
}
set points [MeshIo::TclObjarrayToPythonArrayPoints $node_xyzs]
set cells "\["
foreach element_type {linear triangle quadrilateral tetrahedra pyramid prism hexahedra } {
set elements_data [lindex [GiD_Info mesh elements $element_type -array2] 0]
if { [llength $elements_data] } {
lassign $elements_data element_type_ret element_ids connectivities materials
set num_elements_block [objarray length $element_ids]
if { $num_elements_block } {
if { $nodes_renumbered_for_meshio } {
objarray renumber $connectivities $node_ids $new_node_ids
} else {
objarray incr $connectivities -1 ;#meshio is zero based
}
set some_element_id [objarray get $element_ids 0]
set element_num_nodes [llength [GiD_Mesh get element $some_element_id connectivities]]
set meshio_etype $meshio_element_name($element_type,$element_num_nodes)
append cells "(\"$meshio_etype\",[MeshIo::TclObjarrayToPythonArrayConnectivities $element_num_nodes $connectivities]),"
}
}
}
append cells "\]"
set result [GiD_Python_Call gid_meshio.my_meshio_write_mesh $points $cells $filename]
return 0
} |
...
Code Block | ||
---|---|---|
| ||
import numpy as np import tohil import meshio #to create functions and variables for all tcl available ones tcl=tohil.import_tcl() gid_to_meshio_type = { "sphere,1":"vertex", "point,1":"vertex", "line,2":"line", "triangle,3":"triangle", "quadrilateral,4":"quad", "tetrahedra,4":"tetra", "hexahedra,8":"hexahedron", "prism,6":"wedge", "pyramid,5":"pyramid", #quadratic "line,3":"line3", "triangle,6":"triangle6", "quadrilateral,8":"quad8", "quadrilateral,9":"quad9", "tetrahedra,10":"tetra10", "hexahedra,20":"hexahedron20", "hexahedra,27":"hexahedron27", "prism,15":"wedge15", "pyramid,13":"pyramid13", } def gid_points_to_meshio_points(node_xyzs): xs,ys,zs=node_xyzs num_nodes=len(xs) points=np.empty((num_nodes,3)) for i_node in range(num_nodes): points[i_node]=(float(xs[i_node]),float(ys[i_node]),float(zs[i_node])) #points[:, 0] = xs[:] #points[:, 1] = ys[:] #points[:, 2] = zs[:] return points def gid_elements_to_meshio_cells(element_num_nodes,connectivities): cells=[] #operator // is for integer division num_elements=len(connectivities)//element_num_nodes i=0 for i_element in range(num_elements): node_ids=[] for i_node in range(element_num_nodes): node_ids.append(connectivities[i]) i+=1 cells.append(node_ids) return cells def tohil_obj_array_int_to_numpy(items): num_items=len(items) numpy_array=np.empty(num_items,np.int64) for i in range(num_items) : numpy_array[i]=int(items[i]) return numpy_array def numpy_renumber(connectivities,old_node_ids,new_node_ids): fail=0 length_connectivities=len(connectivities) length_old_ids=len(old_node_ids) length_new_ids=len(new_node_ids) if(length_old_ids==length_new_ids): max_old_id=np.max(old_node_ids) new_number=np.empty(max_old_id,np.int32) for i in range(length_old_ids) : new_number[old_node_ids[i]]=new_node_ids[i] for i in range(length_connectivities) : connectivities[i]=new_number[connectivities[i]] else: fail=1 return connectivities def numpy_incr(connectivities,increment): connectivities=connectivities+increment return connectivities #similar to my_meshio_write_mesh but asking GiD data from python and processing this data here def my_meshio_write_mesh2(filename): info_nodes=tuple(tcl.GiD_Info('mesh','nodes','-array')) node_ids,node_xyzs=info_nodes #tcl.W(node_ids) #tcl.W(node_xyzs) max_id_nodes=int(node_ids[-1]) num_nodes=len(node_ids) nodes_renumbered_for_meshio=False if (max_id_nodes != num_nodes): nodes_renumbered_for_meshio=True new_node_ids=np.arange(num_nodes) points=gid_points_to_meshio_points(node_xyzs) cells = [] for element_type in ['linear','triangle','quadrilateral','tetrahedra','pyramid','prism','hexahedra']: info_elements=tuple(tcl.GiD_Info('mesh','elements',element_type,'-array2')) if (len(info_elements)): #tcl.W(info_elements) elements_data=info_elements[0] element_type_ret,element_ids_original,connectivities_original,materials=elements_data element_ids=tohil_obj_array_int_to_numpy(tuple(element_ids_original)) connectivities=tohil_obj_array_int_to_numpy(tuple(connectivities_original)) num_elements_block=len(element_ids) if (num_elements_block): if (nodes_renumbered_for_meshio): connectivities=numpy_renumber(connectivities,node_ids,new_node_ids) else: connectivities=numpy_incr(connectivities,-1) some_element_id=element_ids[0] else: element_num_nodes=len(tcl.GiD_Mesh('get','element',some_element_id,'connectivities')) connectivitieskey=numpy_incr(connectivities,-1element_type+','+str(element_num_nodes) some_element_id=element_ids[0meshio_etype=gid_to_meshio_type[key] element_num_nodes=len(tcl.GiD_Mesh('get','element',some_element_id,'connectivities'cells.append((meshio_etype,gid_elements_to_meshio_cells(element_num_nodes,connectivities))) mesh=meshio.Mesh(points,cells) result=mesh.write(filename) key=element_type+','+str(element_num_nodes) meshio_etype=gid_to_meshio_type[key] cells.append((meshio_etype,gid_elements_to_meshio_cells(element_num_nodes,connectivities))) mesh=meshio.Mesh(points,cells) result=mesh.write(filename) return resultreturn result |
Note that in this case use python commands to call Tcl like these:
import tohil
tcl=tohil.import_tcl()
info_nodes=tuple(tcl.GiD_Info('mesh','nodes','-array'))
info_elements=tuple(tcl.GiD_Info('mesh','elements',element_type,'-array2'))
#tcl.W(node_ids)
#to show information in a GiD message window for debug
element_num_nodes=len(tcl.GiD_Mesh('get','element',some_element_id,'connectivities'))