diff --git a/seissolxdmf/seissolxdmf/seissolxdmf.py b/seissolxdmf/seissolxdmf/seissolxdmf.py index e83f378..29ac409 100644 --- a/seissolxdmf/seissolxdmf/seissolxdmf.py +++ b/seissolxdmf/seissolxdmf/seissolxdmf.py @@ -2,29 +2,33 @@ import os import xml.etree.ElementTree as ET + def find_line_number_endtag_xdmf(alines): for n, line in enumerate(alines): - if '' in line: - return n + 1 + if "" in line: + return n + 1 raise ValueError(" not found") + class seissolxdmf: def __init__(self, xdmfFilename): self.xdmfFilename = xdmfFilename - with open (xdmfFilename, "r") as fid: - lines=fid.readlines() + with open(xdmfFilename, "r") as fid: + lines = fid.readlines() # Remove potential extra content at the end of the file nlines = find_line_number_endtag_xdmf(lines) - if nlines!=len(lines): - print(f'Warning: extra content at the end of {xdmfFilename} detected') - file_txt = ' '.join([line for line in lines[0:nlines]]) + if nlines != len(lines): + print(f"Warning: extra content at the end of {xdmfFilename} detected") + file_txt = " ".join([line for line in lines[0:nlines]]) self.tree = ET.ElementTree(ET.fromstring(file_txt)) self.ndt = self.ReadNdt() self.nElements = self.ReadNElements() - def ReadHdf5DatasetChunk(self, absolute_path, hdf5var, firstElement, nchunk, idt=-1): - """ Read block of data in hdf5 format - idt!=-1 loads only one time step """ + def ReadHdf5DatasetChunk( + self, absolute_path, hdf5var, firstElement, nchunk, idt=-1 + ): + """Read block of data in hdf5 format + idt!=-1 loads only one time step""" import h5py lastElement = firstElement + nchunk @@ -53,11 +57,13 @@ def GetDtype(self, data_prec, isInt): else: return np.dtype("d") - def ReadSimpleBinaryFile(self, absolute_path, MemDimension, data_prec, isInt, idt=-1): + def ReadSimpleBinaryFile( + self, absolute_path, MemDimension, data_prec, isInt, idt=-1 + ): """Read block of data in binary format (posix) - idt!=-1 loads only one time step + idt!=-1 loads only one time step this function is a special case of ReadSimpleBinaryFileChunk - but kept for performance reasons """ + but kept for performance reasons""" oneDtMem = True if idt != -1 else False data_type = self.GetDtype(data_prec, isInt) @@ -72,10 +78,19 @@ def ReadSimpleBinaryFile(self, absolute_path, MemDimension, data_prec, isInt, id fid.close() return myData - def ReadSimpleBinaryFileChunk(self, absolute_path, MemDimension, data_prec, isInt, firstElement, nchunk, idt=-1): + def ReadSimpleBinaryFileChunk( + self, + absolute_path, + MemDimension, + data_prec, + isInt, + firstElement, + nchunk, + idt=-1, + ): """Read block of data in binary format (posix) same as ReadSimpleBinaryFile: but reads a subset of the second dimension - idt!=-1 loads only one time step """ + idt!=-1 loads only one time step""" oneDtMem = True if idt != -1 else False data_type = self.GetDtype(data_prec, isInt) @@ -93,7 +108,7 @@ def ReadSimpleBinaryFileChunk(self, absolute_path, MemDimension, data_prec, isIn return myData def GetDataLocationPrecisionNElementsMemDimension(self, attribute): - """ Common function called by ReadTopologyOrGeometry """ + """Common function called by ReadTopologyOrGeometry""" root = self.tree.getroot() for Property in root.findall(".//%s" % (attribute)): nElements = int(Property.get("NumberOfElements")) @@ -106,7 +121,7 @@ def GetDataLocationPrecisionNElementsMemDimension(self, attribute): return [dataLocation, data_prec, nElements, MemDimension] def GetDataLocationPrecisionMemDimension(self, dataName): - """ Common function called by ReadData """ + """Common function called by ReadData""" def get(prop): dataLocation = prop.text @@ -131,9 +146,14 @@ def get(prop): raise NameError("%s not found in dataset" % (dataName)) def ReadTopologyOrGeometry(self, attribute): - """ Common function to read either connect or geometry """ + """Common function to read either connect or geometry""" path = os.path.join(os.path.dirname(self.xdmfFilename), "") - dataLocation, data_prec, nElements, MemDimension = self.GetDataLocationPrecisionNElementsMemDimension(attribute) + ( + dataLocation, + data_prec, + nElements, + MemDimension, + ) = self.GetDataLocationPrecisionNElementsMemDimension(attribute) # 3 for surface, 4 for volume dim2 = MemDimension[1] splitArgs = dataLocation.split(":") @@ -143,21 +163,23 @@ def ReadTopologyOrGeometry(self, attribute): filename, hdf5var = splitArgs myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, 0, dim2) else: - myData = self.ReadSimpleBinaryFile(path + dataLocation, dim2, data_prec, isInt) + myData = self.ReadSimpleBinaryFile( + path + dataLocation, dim2, data_prec, isInt + ) # due to zero padding in SeisSol for memory alignement, the array read may be larger than the actual data myData = myData[0:nElements, :] return myData def ReadConnect(self): - """ Read the connectivity matrice defining the cells """ + """Read the connectivity matrice defining the cells""" return self.ReadTopologyOrGeometry("Topology") def ReadGeometry(self): - """ Read the connectivity matrice defining the cells """ + """Read the connectivity matrice defining the cells""" return self.ReadTopologyOrGeometry("Geometry") def ReadNdt(self): - """ read number of time steps in the file """ + """read number of time steps in the file""" root = self.tree.getroot() ndt = 0 for Property in root.findall(".//Grid"): @@ -189,25 +211,28 @@ def ReadAttributeValue(self, attribute, list_possible_location): raise NameError(f"{attribute} not found in {list_possible_location}") def ReadNElements(self): - """ read number of cell elements of the mesh """ + """read number of cell elements of the mesh""" list_possible_location = ["Domain/Grid/Grid/Topology", "Domain/Grid/Topology"] - value = self.ReadAttributeValue("NumberOfElements", list_possible_location) + value = self.ReadAttributeValue("NumberOfElements", list_possible_location) return int(value) def ReadNNodes(self): - """ read number of vertex of the mesh """ + """read number of vertex of the mesh""" list_possible_location = ["Domain/Grid/Grid/Geometry", "Domain/Grid/Geometry"] - value = self.ReadAttributeValue("NumberOfElements", list_possible_location) + value = self.ReadAttributeValue("NumberOfElements", list_possible_location) return int(value) def ReadNodesPerElement(self): - """ read number of nodes per elements of the mesh """ - list_possible_location = ["Domain/Grid/Grid/Topology/DataItem", "Domain/Grid/Topology/DataItem"] - value = self.ReadAttributeValue("Dimensions", list_possible_location) + """read number of nodes per elements of the mesh""" + list_possible_location = [ + "Domain/Grid/Grid/Topology/DataItem", + "Domain/Grid/Topology/DataItem", + ] + value = self.ReadAttributeValue("Dimensions", list_possible_location) return int(value.split()[1]) def ReadAvailableDataFields(self): - """ read all available data fields, e.g. SRs or P_n """ + """read all available data fields, e.g. SRs or P_n""" root = self.tree.getroot() availableDataFields = set() for Property in root.findall(".//Attribute"): @@ -215,7 +240,7 @@ def ReadAvailableDataFields(self): return availableDataFields def ReadTimeStep(self): - """ reading the time step (dt) in the xdmf file """ + """reading the time step (dt) in the xdmf file""" root = self.tree.getroot() i = 0 for Property in root.findall("Domain/Grid/Grid/Time"): @@ -228,52 +253,84 @@ def ReadTimeStep(self): raise NameError("time step could not be determined") def Read1dData(self, dataName, nElements, isInt=False): - """ Read 1 dimension array (used by ReadPartition) """ + """Read 1 dimension array (used by ReadPartition)""" path = os.path.join(os.path.dirname(self.xdmfFilename), "") - dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName) + ( + dataLocation, + data_prec, + MemDimension, + ) = self.GetDataLocationPrecisionMemDimension(dataName) splitArgs = dataLocation.split(":") isHdf5 = True if len(splitArgs) == 2 else False if isHdf5: filename, hdf5var = splitArgs - myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, 0, MemDimension) + myData = self.ReadHdf5DatasetChunk( + path + filename, hdf5var, 0, MemDimension + ) else: filename = dataLocation - myData = self.ReadSimpleBinaryFile(path + dataLocation, nElements, data_prec, isInt, 0) + # The Dimensions of the 1D arrays (e.g. partition) in the xdmf files + # are not not accounting for potential 0 padding (bug in xdmfwriter?). + # Therefore, we precise that there is only one "time-step" to read at location 0 + # by setting idt=0 + myData = self.ReadSimpleBinaryFile( + path + dataLocation, nElements, data_prec, isInt, 0 + ) return myData def ReadPartition(self): - """ Read partition array """ + """Read partition array""" partition = self.Read1dData("partition", self.nElements, isInt=True) return partition def ReadData(self, dataName, idt=-1): - """ Load a data array named 'dataName' (e.g. SRs) + """Load a data array named 'dataName' (e.g. SRs) if idt!=-1, only the time step idt is loaded - else all time steps are loaded """ + else all time steps are loaded""" - return self.ReadDataChunk(dataName, firstElement=0, nchunk=self.nElements, idt=idt) + return self.ReadDataChunk( + dataName, firstElement=0, nchunk=self.nElements, idt=idt + ) def ReadDataChunk(self, dataName, firstElement, nchunk, idt=-1): - """ Load a chunk of a data array named 'dataName' (e.g. SRs) + """Load a chunk of a data array named 'dataName' (e.g. SRs) That is instead of loading 0:nElements, load firstElement:firstElement+nchunk This function is used for generating in parallel Ground motion estimate maps if idt!=-1, only the time step idt is loaded - else all time steps are loaded """ + else all time steps are loaded""" path = os.path.join(os.path.dirname(self.xdmfFilename), "") - dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName) + ( + dataLocation, + data_prec, + MemDimension, + ) = self.GetDataLocationPrecisionMemDimension(dataName) splitArgs = dataLocation.split(":") isHdf5 = True if len(splitArgs) == 2 else False oneDtMem = True if idt != -1 else False if isHdf5: filename, hdf5var = splitArgs - myData = self.ReadHdf5DatasetChunk(path + filename, hdf5var, firstElement, nchunk, idt) + myData = self.ReadHdf5DatasetChunk( + path + filename, hdf5var, firstElement, nchunk, idt + ) else: - myData = self.ReadSimpleBinaryFileChunk(path + dataLocation, MemDimension, data_prec, isInt=False, firstElement=firstElement, nchunk=nchunk, idt=idt) + myData = self.ReadSimpleBinaryFileChunk( + path + dataLocation, + MemDimension, + data_prec, + isInt=False, + firstElement=firstElement, + nchunk=nchunk, + idt=idt, + ) return myData def LoadData(self, dataName, nElements, idt=0, oneDtMem=False, firstElement=-1): - """ Do the same as ReadDataChunk. here for backward compatibility """ - dataLocation, data_prec, MemDimension = self.GetDataLocationPrecisionMemDimension(dataName) + """Do the same as ReadDataChunk. here for backward compatibility""" + ( + dataLocation, + data_prec, + MemDimension, + ) = self.GetDataLocationPrecisionMemDimension(dataName) if not oneDtMem: idt = -1 if firstElement == -1: