Source code for steps.API_1.utilities.meshio

####################################################################################
#
#    STEPS - STochastic Engine for Pathway Simulation
#    Copyright (C) 2007-2023 Okinawa Institute of Science and Technology, Japan.
#    Copyright (C) 2003-2006 University of Antwerp, Belgium.
#    
#    See the file AUTHORS for details.
#    This file is part of STEPS.
#    
#    STEPS is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License version 3,
#    as published by the Free Software Foundation.
#    
#    STEPS is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#    GNU General Public License for more details.
#    
#    You should have received a copy of the GNU General Public License
#    along with this program. If not, see <http://www.gnu.org/licenses/>.
#
#################################################################################   
###
from __future__ import  print_function

"""
.. Note::

    This module is preliminary, means some of the functions are still under development.
    Code modification / debugging is wellcomed.
    Please email steps.dev@gmail.com if you would like to share you changes with others.

Mesh Import and Export Utilities

The MeshIO utilities are designed for importing geometry data from outside mesh-generators,
as well as creating STEPS mesh object and exporting it using STEPS' own format.

Current importing support includes:

* Abaqus data format (http://www.simulia.com)

  Abaqus format is supported by CUBIT (http://cubit.sandia.gov/)
  and Netgen (http://sourceforge.net/projects/netgen-mesher/)

* Tengen's data formats (http://tetgen.berlios.de/)

  Currently, only the .node, .ele and .face files are supported,
  other Tetgen files is not used in STEPS.

Each import function (either importTetgen or importAbaqus) will read and import data
from given data file(s), returns a Tetmesh object for STEPS simulation, and three
ElementProxy objects where geometry data and id mapping are stored.

Once a Tetmesh object has been created the data can be saved to two separate files
by the saveTetmesh() function:

* An xml annotated file containing the nodes, triangles and tetrahedra.

* A text file containing further information needed by STEPS solvers found
  in Tetmesh object constructor.

These files can later be loaded by the loadTetmesh() method and reconstruct
the Tetmesh object. This is intened to drastically reduce mesh-loading time
for large meshes (over ~100,000 voxels). By storing all data required by STEPS
internallyin these two files this infomation does not have to be found each time
by the Tetmesh object constructor.

**IMPORTANT NOTICE**

User is recommanded to save the tetmesh objects using the saveTetmesh() method
and recreate the objects from the saved files, instead of creating the objects
via the importing functions each time, if the tetmesh objects are intented to
be used in multiple simulations. Since the importing functions require a massive
amount of time to create the Tetmesh object, comparing to the loadTetmesh() method.

"""

import time
import re
import steps.API_1.geom as stetmesh
import os.path as opath
from steps.API_1.utilities.steps_shadow import *

#############################################################################################

[docs]class ElementProxy: """ Element Proxy Object The Element Proxy object contains data imported from outside files (Abaqus, Tetgen), as well as the mapping from the element's import id to its STEPS id, and vice versa. it also records block and/or group information contained in the imported file, which could be used for data tracing and compartment/patch creactions:: ElementProxy(type, unitlength) Construct an empty ElementProxy object. Parameters: * type (string): The type of the elements stored in this object. * unitlength (int): The length of vector units required for each element data. e.g. A node with 3 coordinates requires unitlength of 3, a tetrahedron with 4 nodes requires unitlength of 4. A triangle requires unitlength of 3. Example:: nodeproxy = ElementProxy('node', 3) """ def __init__(self, type, unitlength): self.type = type self.unitlength = unitlength self.idcounter = 0 self.data = [] self.importid = [] self.stepsid = {} self.blocks = {} self.groups = {} self.tempblockname = '' self.tempblockstart = -1
[docs] def insert(self, import_id, import_data): """ Insert an element to the Element Map object, a STEPS id will be assigned automatically. Parameters: * import_id (int): The id of the element given by the importing file. * import_data (list): a list of data belongs to the element, e.g. for a node, a 3 units list of its coordinates. Example: nodeproxy.insert(1, [0.1, 0.2, 0.4]) Insert a node with import id 1 and coordinates x = 0.1, y = 0.2, z = 0.4 to nodeproxy. type nodeproxy.getSTEPSID(1) will return the STEPS id of this element. """ self.data += import_data self.importid.append(import_id) self.stepsid[import_id] = self.idcounter self.idcounter += 1
[docs] def getType(self): """ Return the type of elements stored in this Element Map object. """ return self.type
[docs] def getSize(self): """ Return the number of elements stored in the Element Map object. """ return self.idcounter
[docs] def getDataFromSTEPSID(self, steps_id): """ Return the data of element with STEPS id steps_id. Parameters: * steps_id (int): The STEPS id of the element, this id is automatically allocated when inserting the element to the Element Map (see insert() method). """ return self.data[steps_id * self.unitlength : (steps_id + 1) * self.unitlength]
[docs] def getDataFromImportID(self, import_id): """ Return the data of element with import id import_id. Parameters: import_id (int): The import id of the element, this id is given by outside generator and contained in the imported file. """ return self.getDataFromSTEPSID(self.stepsid[import_id])
[docs] def getAllData(self): """ Return all data as an one dimentional list. This list can be used to construct the Tetmesh. Example:: nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() mesh = steps.geom.Tetmesh(nodedata, tetdata) """ return self.data
[docs] def getSTEPSID(self, import_id): """ Return the STEPS id of a element from its import id. Parameters: * import_id (int): The import id of the element. """ return self.stepsid[import_id]
[docs] def getImportID(self, steps_id): """ Return the import id of a element from its STEPS id. Parameters: * steps_id (int): The STEPS id of the element. """ return self.importid[steps_id]
[docs] def addGroup(self, group_name, group_ids): """ Add a group of elements in the Element Map object. An group is defined as a list of element ids. Unlike blocks (see blockBegin() method), the ids of a group need not be continual. For this reason, user should provide a list which contains all ids of elements belong to this group. Each group has a unique name for data access. No uniqueness or order is predefined for a group, this can be defined by the user for specific usage. Groups can be used to construct compartments and patches in STEPS, please refer to the user manual for examples about using groups and blocks. Parameters: * group_name (string): Name of the element group. * group_ids (list): A list of ids of elements belong to the group. """ self.groups[group_name] = group_ids
def addToGroup(self, group_name, id): if (group_name in self.groups): self.groups[group_name].append(id) else: self.groups[group_name] = [id]
[docs] def getGroups(self): """ Return all groups stored in the Element Map object. Example: The following group dictionary contains two groups:: "Group1" contains element ids 1,2,6,4,8 "Group2" contains element ids 3,1,2,1 groups = {"Group1":[1,2,6,4,8], "Group2":[3,1,2,1]} Note that element 1 appears twice in "Group2". For meshes imported from Abaqus files, the keys of the groups are fetched from the ELSET id. For meshes imported from Gmsh2.2 files, the keys of the groups are pairs of (tag_idx, tag_value), where tag_idx is the gmsh tag index position of the element, and tag_value is the tag value of the element. For example, if a tetrahedron is defined as below in Gmsh 1 4 2 2 3 16 17 13 18 It has two tags, tag[0] = 2 and tag[1] = 3, thus the element is contained in both groups[(0, 2)] and groups[(1, 3)]. Note that if a specfic name string is mapped to a value in the physical tag (tag[0]) in the file, the name string can also be used to access the same group, for example, if the following is defined in the file $PhysicalNames 1 3 2 "inner" $EndPhysicalNames groups[(0, 2)] can also be accessed using groups[(0, "inner")] """ return self.groups
[docs] def getGroup(self, group_key): """ Return the group stored in the Element Map object with key = group_key. If the key is not present in the Element Map object, an empty list is returned. Example: g1 = getGroup("Group1") g2 = getGroup("Group2") """ if group_key not in self.groups.keys(): return [] else: return self.groups[group_key]
[docs] def blockBegin(self, block_name): """ Notify the Element Map object that a new block with name block_name begins from the next element. A block is a special type of group whose element's STEPS ids is continual, i.e. can be represented by a begin id and an end id. An unique name should be given to each block, which can be used to access associated block information (i.e. a two unit list contains the block begin id and block end id). If another block was initialized before calling the blockEnd() method, this method will also end the previous block and the its data will be stored. A dictionary of blocks can be converted to a dictionary of group via the blocksToGroups() method. Notice: Although a block will be finished and stored when a new block is notified without calling the blockEnd() method, user is recommanded to call blockEnd() in the usual case. Parameters: * block_name (string): Name of the element block. """ if (self.tempblockstart != -1): self.blocks[self.tempblockname] = [self.tempblockstart, self.idcounter - 1] self.tempblockname = block_name self.tempblockstart = self.idcounter
[docs] def blockEnd(self): """ Notify the Element Map object that the current element block should be ended. Block information (i.e. a dictionary element blockname:[begin_id, end_id]) will be added to the ElementProxy object. Notice: If the blockEnd() method is called before any block begins, or a block has been ended and there is no new block started, no information will be added to the block dictionary. """ if (self.tempblockstart != -1): self.blocks[self.tempblockname] = [self.tempblockstart, self.idcounter - 1] self.tempblockstart = -1
[docs] def getBlocks(self): """ Return the block dictionary stored in the Element Map object. A block dictionary uses the unique block names as keys, and a list of the start and end ids as values. Example: The following block dictionary contains two blocks: "Block1" starts from element with id of 0 and end at element 100, "Block2" starts from element 101 and ends at element 110:: blocks = { "Block1":[0, 100], "Block2":[100, 110] } To access individual blocks, for example "Block1", type:: block1 = blocks["Block1"] User is recommanded to check Python's manual for the usage of dictionary. """ return self.blocks
[docs] def blocksToGroups(self): """ Convert the block dictionary to a group dictionary. Return: A dictionary of groups. """ it = iter(self.blocks) converted_groups = {} for key in it: blockrange = self.blocks[key] converted_groups[key] = list(range(blockrange[0], blockrange[1] + 1)) return converted_groups
#############################################################################################
[docs]def importTetGen(pathroot, scale, verbose = False): """ Read a TetGen-generated or refined mesh from a set of files. The TetGen file format for describing a mesh actually consists of multiple files with varying suffixes. Currently, this routine only reads meshes consisting of three files: * <input>.node: describing the tetrahedral mesh node points. * <input>.ele: describing tetrahedral elements, each of which consists of 4 pointers into the <input>.node file. (TetGen also supports 10-node elements; these 6 extra nodes are obviously not read by STEPS.) * <input>.face: describing triangular faces, each of which consists of 3 pointers into the <input>.node file. This file is optional. Other files are .vol (list of maximum volumes), .var (variant constraints) .neigh (list of neighbours), .smesh (simple PLC descriptions) and .edge (list of boundary edges) files. None of these seem relevant for our use cases, so we don't load them even when they are there. In particular, the .neigh file is computed by STEPS itself. Please refer to the TetGen manual (pages 31-40 in the last edition) for more information on these file formats tetgen.berlios.de/files/tetgen-manual.pdf (See the documentation for steps.geom.tetmesh for details about the mesh object.) PARAMETERS: * pathroot The root of the path name. E.g. mesh/torus would make this routine try to read files mesh/torus.node, mesh/torus.ele and optionally for mesh/torus.face * scale: LENGTH scale from the importing mesh to real geometry. e.g. a radius of 10 in the importing file to a radius of 1 micron in STEPS, scale is 1e-7. RETURNS: mesh, nodeproxy, tetproxy, triproxy mesh: The STEPS TetMesh object nodeproxy: Element Map for nodes tetproxy: Element Map for tetrahedrons triproxy: Element Map for triangles IMPORTANT NOTICE: User is recommanded to save the tetmesh objects using the saveTetmesh() method and recreate the objects from the saved files, instead of creating the objects via the importing functions each time, if the tetmesh objects are intented to be used in multiple simulations. Since the importing functions require a massive amount of time to create the Tetmesh object, comparing to the loadTetmesh() method. """ nodefname = pathroot + '.node' elefname = pathroot + '.ele' facefname = pathroot + '.face' nodeproxy = ElementProxy('node', 3) tetproxy = ElementProxy('tet', 4) triproxy = ElementProxy('tri', 3) # Is there a .node file? if not opath.isfile(nodefname): print(nodefname) return None if not opath.isfile(elefname): print(elefname) return None if not opath.isfile(facefname): facefname = '' # Try to read the node file. with open(nodefname, 'r') as nodefile: # First line is: <x> <y> <z> [att<# of points> <dimension (3)> <# of attributes> # <boundary marker (0 or 1)> line = nodefile.readline() while (line[0] == '#' or line[0] == '\n'): line = nodefile.readline() tokens = line.split() assert len(tokens) == 4 nnodes = int(tokens[0]) assert nnodes > 0 ndims = int(tokens[1]) assert ndims == 3 nattribs = int(tokens[2]) bmarkers = int(tokens[3]) while (nodeproxy.getSize() != nnodes): line = nodefile.readline() commentstart = line.find('#') if commentstart != -1: line = line[0:commentstart] # Remaing lines: <point #>ributes] # [boundary marker] tokens = line.split() if len(tokens) == 0: continue nodeid = int(tokens[0]) node = [0.0, 0.0, 0.0] node[0] = float(tokens[1]) * scale node[1] = float(tokens[2]) * scale node[2] = float(tokens[3]) * scale nodeproxy.insert(nodeid, node) # Try to read the tet file. with open(elefname, 'r') as elefile: line = elefile.readline() while (line[0] == '#' or line[0] == '\n'): line = elefile.readline() tokens = line.split() assert len(tokens) == 3 ntets = int(tokens[0]) assert ntets > 0 ndims = int(tokens[1]) assert ndims == 4 nattribs = int(tokens[2]) while (tetproxy.getSize() != ntets): line = elefile.readline() commentstart = line.find('#') if commentstart != -1: line = line[0:commentstart] # Remaing lines: <point #>ributes] # [boundary marker] tokens = line.split() if len(tokens) == 0: continue tetid = int(tokens[0]) tet = [0, 0, 0, 0] tet[0] = nodeproxy.getSTEPSID(int(tokens[1])) tet[1] = nodeproxy.getSTEPSID(int(tokens[2])) tet[2] = nodeproxy.getSTEPSID(int(tokens[3])) tet[3] = nodeproxy.getSTEPSID(int(tokens[4])) tetproxy.insert(tetid, tet) if nattribs == 1: tetproxy.addToGroup(tokens[5], tetproxy.getSTEPSID(tetid)) # Try to read the tri file. if (facefname != ''): with open(facefname, 'r') as facefile: line = facefile.readline() while (line[0] == '#' or line[0] == '\n'): line = facefile.readline() tokens = line.split() assert len(tokens) == 2 ntris = int(tokens[0]) assert ntris > 0 nattribs = int(tokens[1]) while (triproxy.getSize() != ntris): line = facefile.readline() commentstart = line.find('#') if commentstart != -1: line = line[0:commentstart] # Remaing lines: <point #>ributes] # [boundary marker] tokens = line.split() if len(tokens) == 0: continue triid = int(tokens[0]) tri = [0, 0, 0] tri[0] = nodeproxy.getSTEPSID(int(tokens[1])) tri[1] = nodeproxy.getSTEPSID(int(tokens[2])) tri[2] = nodeproxy.getSTEPSID(int(tokens[3])) triproxy.insert(triid, tri) if nattribs == 1: triproxy.addToGroup(tokens[4], triproxy.getSTEPSID(triid)) if (verbose): print("Read TetGen files succesfully") nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() tridata = triproxy.getAllData() if (verbose): print("creating Tetmesh object in STEPS...") mesh = stetmesh.Tetmesh(nodedata, tetdata, tridata) if (verbose): print("Tetmesh object created.") return mesh, nodeproxy, tetproxy, triproxy
############################################################################################# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ############################################################################################# def parseAbaqusLine(line): line = line.replace("\n", "") line = line.replace("\r", "") type = "undefined" result = {} if line == "": pass elif (line.find('**', 0, 2) == 0): type = "comment" # keyword elif (line.find('*', 0, 1) == 0): type = "keyword" line = line.split(",") for seg in line: if seg[0] == "*": result["keyword"] = seg[1:].upper() elif "=" in seg: seg = seg.replace(" ", "").split("=") result[seg[0].upper()] = seg[1] else: result[seg.upper()] = None # data else: type = "data" result["data"] = line.replace(" ", "").split(",") return type, result ############################################################################################# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #############################################################################################
[docs]def importAbaqus(filename, scale, ebs = None, shadow_mesh = None, verbose = False): """ Read a ABAQUS-formated mesh file, return the created steps.geom.Tetmesh object, the element mapping for nodes, tetraedrons and triangles. PARAMETERS: * filename: the Abaqus filename (or path) including any suffix. * scale: LENGTH scale from the importing mesh to real geometry. e.g. a radius of 10 in the importing file to a radius of 1 micron in STEPS, scale is 1e-7. * ebs: specify the names of selected element blocks which are included in the mesh. * shadow_mesh: name of the ShadowMesh file exported using the STEPS-CUBIT supporting toolkit, can also be the ShadowMesh object itself. RETURNS: mesh, nodeproxy, tetproxy, triproxy * mesh: The STEPS TetMesh object * nodeproxy: Element Map for nodes * tetproxy: Element Map for tetrahedrons * triproxy: Element Map for triangles IMPORTANT NOTICE: User is recommanded to save the tetmesh objects using the saveTetmesh() method and recreate the objects from the saved files, instead of creating the objects via the importing functions each time, if the tetmesh objects are intented to be used in multiple simulations. Since the importing functions require a massive amount of time to create the Tetmesh object, comparing to the loadTetmesh() method. """ if (verbose): print("Reading Abaqus file...") btime = time.time() with open(filename, 'r') as abaqusfile: nodeproxy = ElementProxy('node', 3) tetproxy = ElementProxy('tet', 4) triproxy = ElementProxy('tri', 3) line = abaqusfile.readline() # check we have the right kind of CUBIT output file here #assert(line.find('*HEADING')) if ebs == None: include = True else: include = False currmap = None # add tri data record for remove duplicated triangles recorded_tris = {} # read line one by one until the end while(line): type, result = parseAbaqusLine(line) # status check if type == "comment": pass elif type == "keyword" and (result["keyword"] == "NODE"): if (currmap != None): currmap.blockEnd() if "ELSET" in result.keys(): nodeproxy.blockBegin(result["ELSET"]) else: nodeproxy.blockBegin("AllNodes") currmap = nodeproxy elif type == "keyword" and (result["keyword"] == "ELEMENT"): if (currmap != None): currmap.blockEnd() if "ELSET" in result.keys(): blockname = result["ELSET"] else: blockname = "AllElements" if (ebs == None) or (blockname in ebs): include = True else: include = False if (result["TYPE"] == 'C3D4') and include: currmap = tetproxy tetproxy.blockBegin(blockname) elif (result["TYPE"] == 'STRI3') and include: currmap = triproxy triproxy.blockBegin(blockname) elif currmap != None: currmap.blockEnd() currmap = None elif (type == "keyword") and (currmap != None): currmap.blockEnd() currmap = None elif type == "data" and currmap != None: if (currmap == nodeproxy): nodeid = int(result["data"][0]) node = [0.0, 0.0, 0.0] # set to the right scale node[0] = float(result["data"][1])*scale node[1] = float(result["data"][2])*scale node[2] = float(result["data"][3])*scale currmap.insert(nodeid, node) elif (currmap == tetproxy) and include: nodeid = int(result["data"][0]) node = [0,0,0,0] node[0] = nodeproxy.getSTEPSID(int(result["data"][1])) node[1] = nodeproxy.getSTEPSID(int(result["data"][2])) node[2] = nodeproxy.getSTEPSID(int(result["data"][3])) node[3] = nodeproxy.getSTEPSID(int(result["data"][4])) currmap.insert(nodeid, node) elif (currmap == triproxy) and include: nodeid = int(result["data"][0]) node = [0,0,0] node[0] = nodeproxy.getSTEPSID(int(result["data"][1])) node[1] = nodeproxy.getSTEPSID(int(result["data"][2])) node[2] = nodeproxy.getSTEPSID(int(result["data"][3])) if (node[0], node[1], node[2]) in recorded_tris: if (verbose): print("Triangle: ", (node[0], node[1], node[2]), " with index ", recorded_tris[(node[0], node[1], node[2])], " has been imported, ignore duplicated triangle ", nodeid) else: currmap.insert(nodeid, node) recorded_tris[(node[0], node[1], node[2])] = nodeid line = abaqusfile.readline() if (currmap != None): currmap.blockEnd() nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() tridata = triproxy.getAllData() if (verbose): print("Number of nodes imported: ", nodeproxy.getSize()) if (verbose): print("Number of tetrahedrons imported: ", tetproxy.getSize()) if (verbose): print("Number of triangles imported: ", triproxy.getSize()) if (verbose): print("creating Tetmesh object in STEPS...") mesh = stetmesh.Tetmesh(nodedata, tetdata, tridata) if (verbose): print("Tetmesh object created.") if shadow_mesh != None: if isinstance(shadow_mesh, str): shadow_mesh = ShadowMesh.importFrom(shadow_mesh) if (verbose): print("Importing data from shadow mesh.") for c in shadow_mesh.comps.values(): if (verbose): print("Import compartment ", c.name) steps_indices = [tetproxy.getSTEPSID(i) for i in c.indices] comp = stetmesh.TmComp(c.name, mesh, steps_indices) for v in c.vsyss: comp.addVolsys(v) if len(shadow_mesh.patches.keys()) != 0: print("ImportAbaqus does not support patch importing, use ImportAbaqus2 instead.") for roi in shadow_mesh.rois.items(): roi_name = roi[0] roi_type = roi[1]["Type"] roi_import_indices = roi[1]["Indices"] if (verbose): print("Import ROI data ", roi[0]) if roi_type == ELEM_VERTEX: steps_indices = [nodeproxy.getSTEPSID(i) for i in roi_import_indices] elif roi_type == ELEM_TET: steps_indices = [tetproxy.getSTEPSID(i) for i in roi_import_indices] elif roi_type == ELEM_TRI: print("ImportAbaqus does not support triangle ROI importing, use ImportAbaqus2 instead.") else: steps_indices = roi_import_indices mesh.addROI(roi_name, roi_type, steps_indices) return mesh, nodeproxy, tetproxy,triproxy
[docs]def importAbaqus2(tetfilename, trifilename, scale, shadow_mesh = None, verbose = False): """ Read two ABAQUS-formated mesh files, one with tetrahedron data and the other with triangle data, return the created steps.geom.Tetmesh object, the element mapping for nodes, tetraedrons and triangles. PARAMETERS: * tetfilename: the Abaqus filename for tetrahedron data including any suffix. * trifilename: the Abaqus filename for triangle data including any suffix. * scale: LENGTH scale from the importing mesh to real geometry. e.g. a radius of 10 in the importing file to a radius of 1 micron in STEPS, scale is 1e-7. * shadow_mesh: name of the ShadowMesh file exported using the STEPS-CUBIT supporting toolkit, can also be the ShadowMesh object itself. RETURNS: mesh, nodeproxy, tetproxy, triproxy * mesh: The STEPS TetMesh object * nodeproxy: Element Map for nodes * tetproxy: Element Map for tetrahedrons * triproxy: Element Map for triangles IMPORTANT NOTICE: User is recommanded to save the tetmesh objects using the saveTetmesh() method and recreate the objects from the saved files, instead of creating the objects via the importing functions each time, if the tetmesh objects are intented to be used in multiple simulations. Since the importing functions require a massive amount of time to create the Tetmesh object, comparing to the loadTetmesh() method. """ if (verbose): print("Reading Abaqus file...") btime = time.time() with open(tetfilename, 'r') as tetfile: nodeproxy = ElementProxy('node', 3) tetproxy = ElementProxy('tet', 4) triproxy = ElementProxy('tri', 3) vert_idxs = {} line = tetfile.readline() # check we have the right kind of CUBIT output file here #assert(line.find('*HEADING')) currmap = None # read line one by one until the end while(line): type, result = parseAbaqusLine(line) # status check if type == "comment": pass elif type == "keyword" and (result["keyword"] == "NODE"): if (currmap != None): currmap.blockEnd() #if (verbose): print('Found *NODE section, start reading nodes.') if "ELSET" in result.keys(): nodeproxy.blockBegin(result["ELSET"]) else: nodeproxy.blockBegin("AllNodes") currmap = nodeproxy elif type == "keyword" and (result["keyword"] == "ELEMENT"): if (currmap != None): currmap.blockEnd() if "ELSET" in result.keys(): blockname = result["ELSET"] else: blockname = "AllElements" if result["TYPE"] == 'C3D4': currmap = tetproxy tetproxy.blockBegin(blockname) elif result["TYPE"] == 'STRI3': print("This importing function is not designed for file with both tetrahedrons and triangles, try importAbaqus instead.") return elif currmap != None: currmap.blockEnd() currmap = None elif (type == "keyword") and (currmap != None): currmap.blockEnd() currmap = None elif type == "data" and currmap != None: if (currmap == nodeproxy): nodeid = int(result["data"][0]) node = [0.0, 0.0, 0.0] # set to the right scale node[0] = float(result["data"][1])*scale node[1] = float(result["data"][2])*scale node[2] = float(result["data"][3])*scale currmap.insert(nodeid, node) vert_idxs[tuple(node)] = currmap.getSize() - 1 elif (currmap == tetproxy): nodeid = int(result["data"][0]) node = [0,0,0,0] # set to the right scale node[0] = nodeproxy.getSTEPSID(int(result["data"][1])) node[1] = nodeproxy.getSTEPSID(int(result["data"][2])) node[2] = nodeproxy.getSTEPSID(int(result["data"][3])) node[3] = nodeproxy.getSTEPSID(int(result["data"][4])) currmap.insert(nodeid, node) line = tetfile.readline() if (currmap != None): currmap.blockEnd() with open(trifilename, 'r') as trifile: line = trifile.readline() # check we have the right kind of CUBIT output file here #assert(line.find('*HEADING')) # add tri data record for remove duplicated triangles recorded_tris = {} currmap = None # read line one by one until the end while(line): # status check type, result = parseAbaqusLine(line) if type == "comment": pass elif type == "keyword" and (result["keyword"] == "NODE"): if (currmap != None): currmap.blockEnd() #if (verbose): print('Found *NODE section, start reading nodes.') if "ELSET" in result.keys(): nodeproxy.blockBegin(result["ELSET"]) else: nodeproxy.blockBegin("AllNodes") currmap = nodeproxy elif type == "keyword" and (result["keyword"] == "ELEMENT"): if (currmap != None): currmap.blockEnd() if "ELSET" in result.keys(): blockname = result["ELSET"] else: blockname = "AllElements" if result["TYPE"] == 'STRI3': currmap = triproxy triproxy.blockBegin(blockname) elif (type == "keyword") and (currmap != None): currmap.blockEnd() currmap = None elif type == "data": if (currmap == nodeproxy): nodeid = int(result["data"][0]) node = [0.0, 0.0, 0.0] # set to the right scale node[0] = float(result["data"][1])*scale node[1] = float(result["data"][2])*scale node[2] = float(result["data"][3])*scale if tuple(node) not in vert_idxs.keys(): print("coordinates ", node, " is not in the tet file, abort function.") raise elif (currmap == triproxy): nodeid = int(result["data"][0]) node = [0,0,0] # set to the right scale node[0] = nodeproxy.getSTEPSID(int(result["data"][1])) node[1] = nodeproxy.getSTEPSID(int(result["data"][2])) node[2] = nodeproxy.getSTEPSID(int(result["data"][3])) if (node[0], node[1], node[2]) in recorded_tris: if (verbose): print("Triangle: ", (node[0], node[1], node[2]), " with index ", recorded_tris[(node[0], node[1], node[2])], " has been imported, ignore duplicated triangle ", nodeid) else: currmap.insert(nodeid, node) recorded_tris[(node[0], node[1], node[2])] = nodeid line = trifile.readline() if (currmap != None): currmap.blockEnd() nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() tridata = triproxy.getAllData() if (verbose): print("Number of nodes imported: ", nodeproxy.getSize()) if (verbose): print("Number of tetrahedrons imported: ", tetproxy.getSize()) if (verbose): print("Number of triangles imported: ", triproxy.getSize()) if (verbose): print("creating Tetmesh object in STEPS...") mesh = stetmesh.Tetmesh(nodedata, tetdata, tridata) if (verbose): print("Tetmesh object created.") if shadow_mesh != None: if isinstance(shadow_mesh, str): shadow_mesh = ShadowMesh.importFrom(shadow_mesh) if (verbose): print("Importing data from shadow mesh.") temp_comps = {} for c in shadow_mesh.comps.values(): if (verbose): print("Import compartment ", c.name) steps_indices = [tetproxy.getSTEPSID(i) for i in c.indices] comp = stetmesh.TmComp(c.name, mesh, steps_indices) temp_comps[c] = comp for v in c.vsyss: comp.addVolsys(v) for p in shadow_mesh.patches.values(): if (verbose): print("Import patch ", p.name) steps_indices = [triproxy.getSTEPSID(i) for i in p.indices] icomp = temp_comps[p.icomp] ocomp = None if p.ocomp != None: ocomp = temp_comps[p.ocomp] patch = stetmesh.TmPatch(p.name, mesh, steps_indices, icomp, ocomp) for s in p.ssyss: patch.addSurfsys(s) for roi in shadow_mesh.rois.items(): roi_name = roi[0] roi_type = roi[1]["Type"] roi_import_indices = roi[1]["Indices"] if (verbose): print("Import ROI data ", roi[0]) if roi_type == ELEM_VERTEX: steps_indices = [nodeproxy.getSTEPSID(i) for i in roi_import_indices] elif roi_type == ELEM_TET: steps_indices = [tetproxy.getSTEPSID(i) for i in roi_import_indices] elif roi_type == ELEM_TRI: steps_indices = [triproxy.getSTEPSID(i) for i in roi_import_indices] else: steps_indices = roi_import_indices mesh.addROI(roi_name, roi_type, steps_indices) return mesh, nodeproxy, tetproxy,triproxy
############################################################################################# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #############################################################################################
[docs]def importGmsh(filename, scale, verbose = False): """ Read a Gmsh-formated mesh file, return the created steps.geom.Tetmesh object, the element mapping for nodes, tetraedrons and triangles. PARAMETERS: * filename: the Abaqus filename (or path) including any suffix. * scale: LENGTH scale from the importing mesh to real geometry. e.g. a radius of 10 in the importing file to a radius of 1 micron in STEPS, scale is 1e-7. RETURNS: mesh, nodeproxy, tetproxy, triproxy * mesh: The STEPS TetMesh object * nodeproxy: Element Map for nodes * tetproxy: Element Map for tetrahedrons * triproxy: Element Map for triangles IMPORTANT NOTICE: User is recommanded to save the tetmesh objects using the saveTetmesh() method and recreate the objects from the saved files, instead of creating the objects via the importing functions each time, if the tetmesh objects are intented to be used in multiple simulations. Since the importing functions require a massive amount of time to create the Tetmesh object, comparing to the loadTetmesh() method. """ if (verbose): print("Reading Gmsh file...") btime = time.time() nodeproxy = ElementProxy('node', 3) tetproxy = ElementProxy('tet', 4) triproxy = ElementProxy('tri', 3) physical_name_mapping = {} with open(filename, 'r') as meshfile: line = meshfile.readline() while line: if line == '$MeshFormat\n': line = meshfile.readline() linesec = line.split() if linesec[0] != "2.2" or linesec[1] != "0": raise Exception('This importer only supports Gmsh 2.2 ASCII format.') line = meshfile.readline() assert(line == '$EndMeshFormat\n') elif line == '$PhysicalNames\n': line = meshfile.readline() nphysical_names = int(line) for n in range(nphysical_names): line = meshfile.readline() linesec = line.split() physical_name_mapping[(int(linesec[0]), int(linesec[1]))] = linesec[2].replace('"', '') line = meshfile.readline() assert(line == '$EndPhysicalNames\n') elif line == '$Nodes\n': line = meshfile.readline() n_nodes = int(line) # read nodes for i in range(n_nodes): line = meshfile.readline() linesec = line.split() id = int(linesec[0]) node = [0.0,0.0,0.0] node[0] = float(linesec[1])*scale node[1] = float(linesec[2])*scale node[2] = float(linesec[3])*scale nodeproxy.insert(id, node) line = meshfile.readline() assert(line == '$EndNodes\n') elif line == '$Elements\n': line = meshfile.readline() n_elems = int(line) # read elem for i in range(n_elems): line = meshfile.readline() linesec = line.split() id = int(linesec[0]) type = int(linesec[1]) ntags = int(linesec[2]) # triangle if type == 2: node = [0,0,0] node[0] = nodeproxy.getSTEPSID(int(linesec[-3])) node[1] = nodeproxy.getSTEPSID(int(linesec[-2])) node[2] = nodeproxy.getSTEPSID(int(linesec[-1])) triproxy.insert(id, node) steps_id = triproxy.getSTEPSID(id) for tag in range(ntags): tag_id = int(linesec[3+tag]) triproxy.addToGroup((tag, tag_id), steps_id) # tet if type == 4: group_id = linesec[4] node = [0,0,0,0] node[0] = nodeproxy.getSTEPSID(int(linesec[-4])) node[1] = nodeproxy.getSTEPSID(int(linesec[-3])) node[2] = nodeproxy.getSTEPSID(int(linesec[-2])) node[3] = nodeproxy.getSTEPSID(int(linesec[-1])) tetproxy.insert(id, node) steps_id = tetproxy.getSTEPSID(id) for tag in range(ntags): tag_id = int(linesec[3+tag]) tetproxy.addToGroup((tag, tag_id), steps_id) line = meshfile.readline() assert(line == '$EndElements\n') line = meshfile.readline() for physical_tag in physical_name_mapping: dimension = physical_tag[0] tag = physical_tag[1] tag_name = physical_name_mapping[(dimension, tag)] if dimension == 0: nodeproxy.addGroup((0, tag_name), nodeproxy.getGroup((0, tag))) elif dimension == 2: triproxy.addGroup((0, tag_name), triproxy.getGroup((0, tag))) elif dimension == 3: tetproxy.addGroup((0, tag_name), tetproxy.getGroup((0, tag))) if (verbose): print("Read Msh file succesfully") nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() tridata = triproxy.getAllData() if (verbose): print("Number of nodes imported: ", nodeproxy.getSize()) if (verbose): print("Number of tetrahedrons imported: ", tetproxy.getSize()) if (verbose): print("Number of triangles imported: ", triproxy.getSize()) if (verbose): print("creating Tetmesh object in STEPS...") mesh = stetmesh.Tetmesh(nodedata, tetdata, tridata) if (verbose): print("Tetmesh object created.") return mesh, nodeproxy, tetproxy,triproxy
############################################################################################# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #############################################################################################
[docs]def importVTK(filename, scale, verbose = False): """ Read a VTK-formated mesh file, return the created steps.geom.Tetmesh object, the element mapping for nodes, tetraedrons and triangles. PARAMETERS: * filename: the VTK filename (or path) including any suffix. * scale: LENGTH scale from the importing mesh to real geometry. e.g. a radius of 10 in the importing file to a radius of 1 micron in STEPS, scale is 1e-7. RETURNS: mesh, nodeproxy, tetproxy, triproxy * mesh: The STEPS TetMesh object * nodeproxy: Element Map for nodes * tetproxy: Element Map for tetrahedrons * triproxy: Element Map for triangles IMPORTANT NOTICE: User is recommanded to save the tetmesh objects using the saveTetmesh() method and recreate the objects from the saved files, instead of creating the objects via the importing functions each time, if the tetmesh objects are intented to be used in multiple simulations. Since the importing functions require a massive amount of time to create the Tetmesh object, comparing to the loadTetmesh() method. """ if (verbose): print("Reading VTK file...") #btime = time.time() with open(filename, 'r') as vtkfile: nodeproxy = ElementProxy('node', 3) tetproxy = ElementProxy('tet', 4) triproxy = ElementProxy('tri', 3) elements = [] warningprinted = False line = vtkfile.readline() file_it = iter(vtkfile) for line in file_it: if line.strip() == "": continue tokens = line.split() if tokens[0] == "DATASET" and tokens[1] != "UNSTRUCTURED_GRID": print("Error: not an unstructured grid!") break elif tokens[0] == "POINTS": i = 0 n_points = int(tokens[1]) while i < n_points: line = next(file_it) tokens = line.split() for j in range(len(tokens) // 3): nodeproxy.insert(i, [float(tokens[j])*scale, float(tokens[j + 1])*scale, float(tokens[j + 2])*scale]) i += 1 elif tokens[0] == "CELLS": for i in range(int(tokens[1])): line = next(file_it) tokens = line.split() elements.append(tokens) elif tokens[0] == "CELL_TYPES": for i in range(int(tokens[1])): line = next(file_it) element = elements[i] if int(line) == 5 and int(element[0]) == 3: triproxy.insert(i, [int(element[1]), int(element[2]), int(element[3])]) elif int(line) == 10 and int(element[0]) == 4: tetproxy.insert(i, [int(element[1]), int(element[2]), int(element[3]), int(element[4])]) elif not warningprinted: print("Warning: at least one element is neither a triangle nor a tetrahedron") warningprinted = True break nodedata = nodeproxy.getAllData() tetdata = tetproxy.getAllData() tridata = triproxy.getAllData() if (verbose): print("creating Tetmesh object in STEPS...") mesh = stetmesh.Tetmesh(nodedata, tetdata, tridata) if (verbose): print("Tetmesh object created.") return mesh, nodeproxy, tetproxy, triproxy
############################################################################################# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #############################################################################################
[docs]def saveMesh(pathname, tetmesh): """ Save a STEPS Tetmesh in an XML file This file stores the basic information about the mesh which tends to be common information for any software that supports tetrahedral meshes. * NODES are stored by cartesian coordinates. * TRIANGLES are stored by the indices of their 3 nodes. * TETRAHEDRONS are stored by the indices of their 4 nodes. The XML file also stores infomation about any compartments or patches created in STEPS (class steps.geom.TmComp steps.geom.TmPatch respectively). * COMPARTMENT(S) are stored by: * their string identification. * a list of any volume systems added to the compartment at time of saving. * a list of tetrahedrons belonging to the compartment * PATCH(ES) are stored by: * their string identification. * a list of any surface systems added to the patch at time of saving. * the inner compartment id. * the outer compartment id (if it exists). * a list of trianlges belonging to this patch. PARAMETERS: * pathname: the root of the path to store the files. e.g. 'meshes/spine1' will save data in /meshes/spine1.xml * tetmesh: A valid STEPS Tetmesh object (of class steps.geom.Tetmesh). This mesh can be made in a variety of ways, e.g. to save a mesh loaded from tetgen:: >>> import meshio >>> ### Use cubit script to create steps.geom.Tetmesh object from tetgen output file ### >>> mymesh = meshio.readTetgen(tetgenfilename) >>> ### Save this mesh in XML (and ASCII) format for quick-loading in future ### >>> meshio.saveMesh('/meshes/spine1', mymesh[0]) """ # Perform a basic test on the object itself if not isinstance(tetmesh, stetmesh.Tetmesh): raise TypeError(f"Expected a steps.API_1.geom.Tetmesh object, got {tetmesh} instead.") # Following will raise IOError if pathname not a valid directory with open(pathname+'.xml', 'w') as xmlfile: xmlfile.write('<?xml version="1.0" encoding="ISO-8859-1"?>\n') xmlfile.write('<tetmesh>\n') ### Write the node information nverts = tetmesh.countVertices() xmlfile.write('\t<nodes size="'+str(nverts)+'">\n') for node in range(nverts): # Write indices and coordinates to xml file xmlfile.write('\t\t<node idx = "' + str(node) + '">\n') coords = str(tetmesh.getVertex(node))[1:-1] xmlfile.write('\t\t\t<coords>'+ coords + '</coords>\n') xmlfile.write('\t\t</node>\n') xmlfile.write('\t</nodes>\n') ### Write the triangle information ntris = tetmesh.countTris() xmlfile.write('\t<triangles size="' +str(ntris)+'">\n') for tri in range(ntris): # Write indices and nodes to xml file xmlfile.write('\t\t<tri idx="' + str(tri) + '">\n') nodes = str(tetmesh.getTri(tri))[1:-1] xmlfile.write('\t\t\t<nodes>' + nodes + '</nodes>\n') xmlfile.write('\t\t</tri>\n') xmlfile.write('\t</triangles>\n') ### Write the tetrahedron information ntets = tetmesh.countTets() xmlfile.write('\t<tetrahedrons size="' +str(ntets)+'">\n') for tet in range(ntets): # Write indices and nodes to xml file xmlfile.write ('\t\t<tet idx="' + str(tet) + '">\n') nodes = str(tetmesh.getTet(tet))[1:-1] xmlfile.write('\t\t\t<nodes>' + nodes + '</nodes>\n') xmlfile.write('\t\t</tet>\n') # Write the volume, barycenter, trineighbours, tetneighbours to xml file xmlfile.write('\t</tetrahedrons>\n') ### Write the comp and patch information. # TODO: Changes need to be made to steps code to make it easier to get the tet # and tri members. Currently the mesh returns base pointer (Comp or Patch) # which cannot be used to find the indices. comps = tetmesh.getAllComps() ncomps = comps.__len__() xmlfile.write('\t<compartments size="' + str(ncomps) + '">\n') if (ncomps > 0) : ids = [] vsys=[] tets = [] for c in range(ncomps): ids.append(comps[c].getID()) vsys.append(comps[c].getVolsys()) tets.append([]) assert(tets.__len__() == ncomps) # Only choice right now is to loop over all tets and compare comp to tet id for tet in range(ntets): comptemp = tetmesh.getTetComp(tet) if not comptemp: continue idtemp = comptemp.getID() for c in range(ncomps): if idtemp == ids[c]: tets[c].append(tet) break # Now we have the tet members of each comp, we can write this to xml for c in range(ncomps): xmlfile.write('\t\t<comp idx="' +str(c) + '">\n') xmlfile.write('\t\t\t<id>' + ids[c] + '</id>\n') xmlfile.write('\t\t\t<volsys>') for v in vsys[c]: if (v != vsys[c][0]) : xmlfile.write(',') xmlfile.write(v) xmlfile.write('</volsys>\n') xmlfile.write('\t\t\t<tets>') for t in tets[c]: if (t != tets[c][0]) : xmlfile.write(',') xmlfile.write(str(t)) xmlfile.write('</tets>\n') xmlfile.write('\t\t</comp>\n') xmlfile.write('\t</compartments>\n') patches = tetmesh.getAllPatches() npatches = patches.__len__() xmlfile.write('\t<patches size="' + str(npatches) + '">\n') if (npatches > 0) : ids = [] ssys = [] icomp=[] ocomp=[] tris = [] for p in range(npatches): ids.append(patches[p].getID()) ssys.append(patches[p].getSurfsys()) icomp.append(patches[p].getIComp().getID()) if (not patches[p].getOComp()): ocomp.append('null') else : ocomp.append(patches[p].getOComp().getID()) tris.append([]) assert(ids.__len__() == ssys.__len__() == icomp.__len__() == ocomp.__len__() == tris.__len__() == npatches) for tri in range(ntris): patchtemp = tetmesh.getTriPatch(tri) if not patchtemp: continue idtemp = patchtemp.getID() for p in range(npatches): if idtemp == ids[p]: tris[p].append(tri) break # Write all to xml for p in range(npatches): xmlfile.write('\t\t<patch idx = "'+str(p) + '">\n') xmlfile.write('\t\t\t<id>' + ids[p] + '</id>\n') xmlfile.write('\t\t\t<surfsys>') for s in ssys[p]: if (s != ssys[p][0]) : xmlfile.write(',') xmlfile.write(s) xmlfile.write('</surfsys>\n') xmlfile.write('\t\t\t<icomp>' + icomp[p] + '</icomp>\n') xmlfile.write('\t\t\t<ocomp>' + ocomp[p] + '</ocomp>\n') xmlfile.write('\t\t\t<tris>') for t in tris[p]: if (t != tris[p][0] ) : xmlfile.write(',') xmlfile.write(str(t)) xmlfile.write('</tris>\n') xmlfile.write('\t\t</patch>\n') xmlfile.write('\t</patches>\n') # add write ROI data roi_ids = tetmesh.getAllROINames() xmlfile.write('\t<ROI_records size="' + str(len(roi_ids)) + '">\n') for n in roi_ids: type = tetmesh.getROIType(n) xmlfile.write('\t\t<ROI id="'+ n + '" type="' + str(type) + '">\n') indices = tetmesh.getROIData(n) xmlfile.write('\t\t\t<indices>') for i in indices: if (i != indices[0] ) : xmlfile.write(',') xmlfile.write(str(i)) xmlfile.write('</indices>\n') xmlfile.write('\t\t</ROI>\n') xmlfile.write('\t</ROI_records>\n') xmlfile.write('</tetmesh>')
#############################################################################################
[docs]def loadMesh(pathname, scale=1, strict=False): """ Load a mesh in STEPS from the XML file. PARAMETERS: * pathname: the root of the path where the file(s) are stored. e.g. with 'meshes/spine1' this function will look for the file /meshes/spine1.xml * scale: optionally rescale the mesh on loading by given factor. * strict: apply strict(-er) checking to the input XML. RETURNS: A tuple (mesh, comps, patches) * mesh The STEPS Tetmesh object (steps.geom.Tetmesh) * comps A list of any compartment objects (steps.geom.TmComp) from XML file * patches A list of any patches objects (steps.geom.TmPatch) from XML file """ # More robust/relaxed attribute parsing: element_rx=re.compile(r'\s*<(\/?\w+)((?:\s+\w+\s*=\s*(?:\w+|"[^"]*"|\'[^\']*\'))*)\s*/?\s*>') attrib_rx=re.compile(r'\s*(\w+)\s*=\s*(?:(\w+)|"([^"]*)"|\'([^\']*)\')') def parse_xml_element(s): (tag,attrstr) = element_rx.match(s).groups() attrs = {} for a,u,d,q in attrib_rx.findall(attrstr): attrs[a] = u+d+q return (tag,attrs) # Try to open the XML file. An error will be raisen if it doesn't exist with open(pathname+'.xml', 'r') as xmlfile: # Perform a basic check to see if we have the expected kind of file which has not been altered. info = xmlfile.readline() if(xmlfile.readline().rstrip() != '<tetmesh>'): print('XML file is not a recognised STEPS mesh file') return # Collect basic node information and perform some checks on the data read from XML file nodeinfo = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(nodeinfo) assert(tag == 'nodes') assert('size' in attrs) nnodes = int(attrs['size']) nodes_out = [0.0]*(nnodes*3) for i in range(nnodes): idxtemp = xmlfile.readline().strip() assert(int(idxtemp[13:-2]) == i) coordtemp = xmlfile.readline().strip() assert(coordtemp[:8] == '<coords>' and coordtemp[-9:] == '</coords>') coordtemp = coordtemp[8:-9].split(', ') nodes_out[i*3], nodes_out[(i*3)+1], nodes_out[(i*3)+2] = float(coordtemp[0]), float(coordtemp[1]), float(coordtemp[2]) assert(xmlfile.readline().strip() == '</node>') assert(xmlfile.readline().strip() == '</nodes>') # Now read triangle information from xml file and text file if we have it triinfo = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(triinfo) assert(tag == 'triangles') assert('size' in attrs) ntris = int(attrs['size']) tris_out = [0]*(ntris*3) # numpy.zeros(ntris*3, dtype = 'int') for i in range(ntris): idxtemp = xmlfile.readline().strip() if strict: (tag,attrs) = parse_xml_element(idxtemp) assert(tag == 'tri') assert(int(attrs['idx']) == i) nodetemp = xmlfile.readline().strip() assert (nodetemp[:7] == '<nodes>' and nodetemp[-8:] == '</nodes>') nodetemp = nodetemp[7:-8].split(', ') tris_out[i*3], tris_out[(i*3)+1], tris_out[(i*3)+2] = int(nodetemp[0]), int(nodetemp[1]), int(nodetemp[2]) assert(xmlfile.readline().strip() == '</tri>') # Now read the text file, if it exists, and get the extra information assert(xmlfile.readline().strip() == '</triangles>') # Now read tet information from xml file and text file if we have it tetinfo = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(tetinfo) assert(tag == 'tetrahedrons') assert('size' in attrs) ntets = int(attrs['size']) tets_out = [0]*(ntets*4) # numpy.zeros(ntets*4, dtype = 'int') for i in range(ntets): idxtemp = xmlfile.readline().strip() if strict: (tag,attrs) = parse_xml_element(idxtemp) assert(tag == 'tet') assert(int(attrs['idx']) == i) nodetemp = xmlfile.readline().strip() assert (nodetemp[:7] == '<nodes>' and nodetemp[-8:] == '</nodes>') nodetemp = nodetemp[7:-8].split(', ') tets_out[i*4], tets_out[(i*4)+1], tets_out[(i*4)+2], tets_out[(i*4)+3] = int(nodetemp[0]), int(nodetemp[1]), int(nodetemp[2]), int(nodetemp[3]) assert(xmlfile.readline().strip() == '</tet>') assert(xmlfile.readline().strip() == '</tetrahedrons>') # Rescale coordinates if requested. if scale!=1: for i in range(len(nodes_out)): nodes_out[i] *= scale # We have all the information now. Time to make the Tetmesh object. mesh = stetmesh.Tetmesh(nodes_out, tets_out, tris_out) # Now fetch any comp and patch information from XML file compinfo = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(compinfo) assert(tag == 'compartments') assert('size' in attrs) ncomps = int(attrs['size']) comps_out = [] for i in range(ncomps): idxtemp = xmlfile.readline().strip() if strict: (tag,attrs) = parse_xml_element(idxtemp) assert(tag == 'comp') assert(int(attrs['idx']) == i) idtemp = xmlfile.readline().strip() assert(idtemp[:4] == '<id>' and idtemp[-5:] == '</id>') idtemp = idtemp[4:-5] volsystemp = xmlfile.readline().strip() assert(volsystemp[:8] == '<volsys>' and volsystemp[-9:] == '</volsys>') volsystemp = volsystemp[8:-9].split(',') if (volsystemp[0] == '') : volsystemp = [] tettemp = xmlfile.readline().strip() assert(tettemp[:6] == '<tets>' and tettemp[-7:] == '</tets>') tettemp = tettemp[6:-7].split(',') nctets = tettemp.__len__() ctets = [0]*nctets # numpy.zeros(nctets, dtype = 'int') for ct in range(nctets): ctets[ct] = int(tettemp[ct]) c_out = stetmesh.TmComp(idtemp, mesh, ctets) for v in volsystemp: c_out.addVolsys(v) comps_out.append(c_out) assert(xmlfile.readline().strip() == '</comp>') assert(xmlfile.readline().strip() == '</compartments>') # Retrieve patch info patchinfo = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(patchinfo) assert(tag == 'patches') assert('size' in attrs) npatches = int(attrs['size']) patches_out = [] for i in range(npatches): idxtemp = xmlfile.readline().strip() if strict: (tag,attrs) = parse_xml_element(idxtemp) assert(tag == 'patch') assert(int(attrs['idx']) == i) idtemp = xmlfile.readline().strip() assert(idtemp[:4] == '<id>' and idtemp[-5:] == '</id>') idtemp = idtemp[4:-5] surfsystemp = xmlfile.readline().strip() assert(surfsystemp[:9] == '<surfsys>' and surfsystemp[-10:] == '</surfsys>') surfsystemp = surfsystemp[9:-10].split(',') if (surfsystemp[0] == '') : surfsystemp = [] icomptemp = xmlfile.readline().strip() assert(icomptemp[:7] == '<icomp>' and icomptemp[-8:] == '</icomp>') icomptemp = icomptemp[7:-8] ocomptemp = xmlfile.readline().strip() assert(ocomptemp[:7] == '<ocomp>' and ocomptemp[-8:] == '</ocomp>') ocomptemp = ocomptemp[7:-8] tritemp = xmlfile.readline().strip() assert(tritemp[:6] == '<tris>' and tritemp[-7:] == '</tris>') tritemp = tritemp[6:-7].split(',') nptris = tritemp.__len__() ptris = [0]*nptris # numpy.zeros(nptris, dtype='int') for pt in range(nptris): ptris[pt] = int(tritemp[pt]) if (ocomptemp != 'null'): p_out = stetmesh.TmPatch(idtemp, mesh, ptris, mesh.getComp(icomptemp), mesh.getComp(ocomptemp)) else : p_out = stetmesh.TmPatch(idtemp, mesh, ptris, mesh.getComp(icomptemp)) for s in surfsystemp: p_out.addSurfsys(s) patches_out.append(p_out) assert(xmlfile.readline().strip() == '</patch>') assert(xmlfile.readline().strip() == '</patches>') # Retrieve ROI info info = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(info) if tag == 'ROI_records': ROIinfo = info assert('size' in attrs) nROI = int(attrs['size']) for r in range(nROI): ROItemp = xmlfile.readline().strip() (tag,attrs) = parse_xml_element(ROItemp) assert(tag == 'ROI') assert('id' in attrs) assert('type' in attrs) id = attrs['id'] type = int(attrs['type']) idxtemp = xmlfile.readline().strip() assert(idxtemp.find("<indices>") != -1) assert(idxtemp.find("</indices>") != -1) idxtemp = idxtemp.replace("<indices>", "").replace("</indices>", "").split(",") indices = [int(i) for i in idxtemp] assert(xmlfile.readline().strip() == '</ROI>') mesh.addROI(id, type, indices) assert(xmlfile.readline().strip() == '</ROI_records>') info = xmlfile.readline().strip() assert(info == '</tetmesh>') return (mesh,comps_out,patches_out)