# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
"""Xdmf file reader
"""
from __future__ import annotations
from typing import List, Union
import os
import sys
import numpy as np
import xml.sax
from Muscat.Helpers.TextFormatHelper import TFormat
from Muscat.Types import MuscatIndex, MuscatFloat
import Muscat.Containers.ElementsDescription as ED
from Muscat.IO.XdmfTools import FieldNotFound, HasHdf5Support, XdmfNumber
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
from Muscat.Containers.ElementsContainers import ElementsContainer, AllElements
from Muscat.Containers.Mesh import Mesh
from Muscat.IO.IOFactory import RegisterReaderClass
[docs]def ReadXdmf(fileName):
"""Function API for reading an xdmf file
Parameters
----------
fileName : str
name of the file to be read
Returns
-------
Mesh
output unstructured mesh object containing reading result
"""
obj = XdmfReader(filename=fileName)
return obj.Read()
[docs]class XdmfBase(object):
""" Base class for all the xdmf related objects """
[docs] def ReadAttribute(self, attrs, name, default=None):
"""" Helper to read attributes"""
if name in attrs:
return attrs.get(name)
else:
if default is None:
raise KeyError(f"Key :'{name}' not available (you can add a default value to bypass this error)")
else:
return default
[docs] def TreatCDATA(self) -> None:
"""Default implementation to read the heavy data"""
pass
def __AttributesToStr(self) -> str:
""" Helper function to make easier the introspection """
import inspect
res = ''
TFormat.II()
attributes = inspect.getmembers(self, lambda a: not inspect.isroutine(a))
attributes = [a for a in attributes if not (a[0].startswith('__') and a[0].endswith('__')) and a[0][0].isupper()]
# print(attributes)
for a in attributes:
# print(a)
res += TFormat.GetIndent() + str(a[0]) + ' : ' + str(a[1]) + '\n'
TFormat.DI()
return res
[docs]class Xdmf(XdmfBase):
""" Top class for the xdmf Document """
def __init__(self):
self.domains: List[XdmfDomain] = []
[docs] def GetDomain(self, nameOrNumber: Union[str,MuscatIndex]) -> XdmfDomain:
""" Get a Grid (XdmfGrid) using a name or a number """
if isinstance(nameOrNumber, str):
for g in self.domains:
if g.Name == nameOrNumber:
return g
raise Exception("Domain '"+nameOrNumber + "' not Found")
else:
return self.domains[nameOrNumber]
def __str__(self):
res = 'Xdmf\n'
TFormat.II()
for d in self.domains:
res += d.__str__()
TFormat.DI()
return res
[docs]class XdmfDomain(XdmfBase):
"""A Domain. Can contain many grids"""
def __init__(self):
self.grids: List[XdmfGrid] = []
self.Name: str = ''
self.information = []
[docs] def GetGrid(self, NameOrNumber: Union[str, MuscatIndex]) -> XdmfGrid:
""" Get a Grid (XdmfGrid) using a name or a number """
if isinstance(NameOrNumber, str):
for g in self.grids:
if g.Name == NameOrNumber:
return g
raise Exception(f" grid with name '{NameOrNumber}' not found")
else:
return self.grids[NameOrNumber]
[docs] def ReadAttributes(self, attrs):
self.Name = self.ReadAttribute(attrs, 'Name', '')
def __str__(self):
res = TFormat.GetIndent() + 'XdmfDomain\n'
TFormat.II()
for g in self.grids:
res += g.__str__()
TFormat.DI()
return res
[docs]class XdmfGrid(XdmfBase):
""" a Grid: contains a mesh (points and connectivity ) and fields by (element, nodes, grid)"""
def __init__(self):
self.information = []
self.topology: Union[None, XdmfTopology] = None
self.geometry: Union[None,XdmfGeometry] = None
self.attributes:List[XdmfAttribute] = []
self.Name: str = ''
self.GridType: str = 'Uniform'
self.time : Union[None,XdmfTime] = None
[docs] def GetSupport(self) -> Mesh:
"""Returns the support defined in the Xdmf file
Returns
-------
Mesh or ConstantRectilinearMesh
output mesh object containing reading result
"""
if self.geometry is None or self.topology is None:
raise RuntimeError("Cant recover support of an empty (uninitialized) XdmfGrid")
res = Mesh()
if self.geometry.Type == "ORIGIN_DXDYDZ":
swapper = [2, 1, 0]
res = CreateConstantRectilinearMesh(dimensions=self.topology.GetDimensions()[swapper], origin=self.geometry.GetOrigin(), spacing=self.geometry.GetSpacing() )
elif self.geometry.Type == "XYZ":
res.nodes = np.ascontiguousarray(self.geometry.GetNodes().reshape((-1, 3)), dtype=MuscatFloat)
res.elements = self.topology.GetConnectivity()
res.GenerateManufacturedOriginalIDs()
res.PrepareForOutput()
if np.linalg.norm(res.nodes[:, 2]) == 0:
from Muscat.Containers.MeshModificationTools import LowerNodesDimension
res = LowerNodesDimension(res)
elif self.geometry.Type == "XY":
res.nodes = self.geometry.GetNodes().reshape((-1, 2))
res.elements = self.topology.GetConnectivity()
res.GenerateManufacturedOriginalIDs()
res.PrepareForOutput()
else:
raise RuntimeError(f" geometry type unknown ({self.geometry.Type})")
nodeFields = {name: self.GetPointField(name)[0] for name in self.GetPointFieldsNames()}
elemFields = {name: self.GetCellField(name)[0] for name in self.GetCellFieldsNames()}
for name, data in nodeFields.items():
if data.shape[0] == res.GetNumberOfNodes() and data.dtype == np.int8 and np.min(data) >= 0 and np.max(data) <= 1:
res.nodesTags.CreateTag(name).SetIds(np.where(data == 1)[0])
else:
res.nodeFields[name] = data
for name, data in elemFields.items():
if data.shape[0] == res.GetNumberOfElements() and data.dtype == np.int8 and np.min(data) >= 0 and np.max(data) <= 1:
res.AddElementsToTag(np.where(data)[0], tagname=name)
else:
res.elemFields[name] = data
return res
[docs] def ReadAttributes(self, attrs):
self.Name = self.ReadAttribute(attrs, 'Name', default="")
self.GridType = self.ReadAttribute(attrs, 'CollectionType', default="Uniform")
[docs] def HasField(self, name):
for a in self.attributes:
if a.Name == name:
return True
return False
[docs] def GetTime(self):
"""Returns time at which data is read
Returns
-------
float
time at which data is read
"""
if self.time is None:
return None
else:
return self.time.Value[0]
[docs] def GetFieldsNames(self):
return [a.Name for a in self.attributes]
[docs] def GetPointFieldsNames(self) -> List[str]:
return [a.Name for a in self.attributes if a.Center == "Node"]
[docs] def GetCellFieldsNames(self) -> List[str]:
return [a.Name for a in self.attributes if a.Center == "Cell"]
[docs] def GetGridFieldsNames(self) -> List[str]:
return [a.Name for a in self.attributes if a.Center == "Grid"]
[docs] def GetPointFields(self) -> List[np.ndarray]:
return self.GetFieldsOfType("Node")
[docs] def GetPointField(self, name) -> List[np.ndarray]:
return self.GetFieldsOfType("Node", name)
[docs] def GetCellFields(self) -> List[np.ndarray]:
return self.GetFieldsOfType("Cell")
[docs] def GetCellField(self, name) -> List[np.ndarray]:
return self.GetFieldsOfType("Cell", name)
[docs] def GetGridFields(self) -> List[np.ndarray]:
return self.GetFieldsOfType("Grid")
[docs] def GetGridField(self, name) -> List[np.ndarray]:
return self.GetFieldsOfType("Grid", name)
[docs] def GetFieldsOfType(self, ftype, name=None):
res = []
for a in self.attributes:
if name is not None:
if a.Name != name:
continue
if a.Center == ftype:
data = a.dataitems[0].GetData()
if self.geometry.Type == "ORIGIN_DXDYDZ":
if a.Type == "Vector":
data = data.transpose(2, 1, 0, 3)
else:
data.shape = self.topology.GetDimensions()
data = data.transpose(2, 1, 0).ravel()
res.append(np.copy(data))
return res
[docs] def GetFieldData(self, name, noTranspose=False):
for att in self.attributes:
if att.Name == name:
data = att.dataitems[0].GetData()
if self.geometry.Type == "ORIGIN_DXDYDZ":
if att.Type == "Vector":
if noTranspose:
return data
else:
return data.transpose(2, 1, 0, 3)
else:
if noTranspose:
return data
else:
return data.transpose(2, 1, 0)
return data
raise FieldNotFound(name)
[docs] def GetFieldTermsAsColMatrix(self, fieldname, offset=0):
# first we check the padding max 8 zeros of padding
padding = 0
for i in range(8):
padding = i
ss = fieldname+str(offset).zfill(padding)
# print ss
if self.HasField(ss):
break
# print(padding)
else:
raise FieldNotFound(fieldname)
# first we check the number of terms
#print('padding ' + str(padding))
cpt = 0
while (self.HasField(fieldname + str(offset+cpt).zfill(padding))):
cpt += 1
# get the firt data to get the type and the size
d_0 = self.GetFieldData(fieldname+str(offset).zfill(padding))
# now we allocate the np matrix for the data
res = np.empty([d_0.size, cpt], dtype=d_0.dtype)
res[:, 0] = d_0.reshape(d_0.size)
for i in range(1, cpt):
res[:, i] = self.GetFieldData(fieldname+str(offset+i).zfill(padding)).reshape(d_0.size)
return res
[docs] def GetFieldTermsAsTensor(self, fieldname, sep='_', offseti=0, offsetj=0):
from itertools import product
# first we check the padding max 8 zeros of padding
paddingi = 0
paddingj = 0
for i, j in product(range(8), range(8)):
paddingi = i
paddingj = j
ss = fieldname + str(offseti).zfill(paddingi) + sep + str(offsetj).zfill(paddingj)
if self.HasField(ss):
break
else:
raise FieldNotFound(fieldname+"*" + str(offseti) + sep + "*" + str(offsetj))
# first we check the number of terms
cpti = 0
while (self.HasField(fieldname+str(offseti+cpti).zfill(paddingi) + sep + str(offsetj).zfill(paddingj))):
cpti += 1
cptj = 0
while (self.HasField(fieldname+str(offseti).zfill(paddingi) + sep + str(offsetj+cptj).zfill(paddingj))):
cptj += 1
# get the firt data to get the type and the size
d_0_0 = self.GetFieldData(fieldname+str(offseti).zfill(paddingi) + sep + str(offsetj).zfill(paddingj))
# now we allocate the np matrix for the data
res = np.empty([d_0_0.size, cpti, cptj], dtype=d_0_0.dtype)
for i in range(0, cpti):
for j in range(0, cptj):
#print(fieldname + str(offseti+i).zfill(paddingi) + sep + str(offsetj+j).zfill(paddingj))
res[:, i, j] = self.GetFieldData(fieldname + str(offseti+i).zfill(paddingi) + sep + str(offsetj+j).zfill(paddingj)).reshape(d_0_0.size)
return res
def __str__(self):
res = TFormat.GetIndent() + 'XdmfGrid\n'
TFormat.II()
for info in self.information:
res += info.__str__()
res += TFormat.GetIndent() + 'Name : ' + self.Name + '\n'
res += self.geometry.__str__()
res += self.topology.__str__()
for att in self.attributes:
res += att.__str__()
TFormat.DI()
return res
[docs]class XdmfTopology(XdmfBase):
""" XdmfTopology class: stores the connectivity of the Grid"""
def __init__(self):
self.dataitems = []
self.Dimensions = None
self.Type = None
def __str__(self):
res = TFormat.GetIndent() + 'XdmfTopology\n'
res += self._XdmfBase__AttributesToStr()
return res
[docs] def ReadAttributes(self, attrs):
self.Type = self.ReadAttribute(attrs, 'Type', default=-1)
if self.Type == -1:
self.Type = self.ReadAttribute(attrs, 'TopologyType')
# if self.Type != "Mixed":
try:
self.Dimensions = np.array(self.ReadAttribute(attrs, 'Dimensions').split(), dtype='int')[::-1]
except:
pass
# d = self.dataitems[0].ReadAttribute(attrs,'Dimensions').split()
# print(d)
# self.Dimensions = np.array(d, dtype='int')[::-1]
# print(self.Dimensions)
[docs] def GetConnectivity(self) -> AllElements:
from Muscat.IO.XdmfTools import XdmfNumberToEN
from Muscat.IO.XdmfTools import XdmfNameToEN
from Muscat.Containers.ElementsContainers import ElementsContainer, AllElements, StructuredElementsContainer
res = AllElements()
if self.Type == "Mixed":
rawdata = self.dataitems[0].GetData()
ts = len(rawdata)
cpt = 0
elcpt = 0
while cpt < ts:
# print(rawdata)
# print(XdmfNumberToEN)
elementName = XdmfNumberToEN[rawdata[cpt]]
if rawdata[cpt] == 0x2 or rawdata[cpt] == 0x1:
cpt += 1
cpt += 1
elements = res.GetElementsOfType(elementName)
connectivity = np.ascontiguousarray(rawdata[cpt:elements.GetNumberOfNodesPerElement()+cpt], dtype=MuscatIndex)
elements.AddNewElement(connectivity, elcpt)
elcpt += 1
cpt += elements.GetNumberOfNodesPerElement()
return res
elif self.Type == "3DCoRectMesh":
raise
elif self.Type == "3DSMesh":
if len(self.Dimensions) == 3:
data = StructuredElementsContainer(ED.Hexahedron_8)
else:
data = StructuredElementsContainer(ED.Quadrangle_4)
data.SetDimensions(self.Dimensions)
res.AddContainer(data)
return res
else:
en = XdmfNameToEN[self.Type]
data = ElementsContainer(en)
data.connectivity = np.ascontiguousarray(self.dataitems[0].GetData().ravel(), dtype=MuscatIndex)
size = (len(data.connectivity)//ED.numberOfNodes[en], ED.numberOfNodes[en])
data.connectivity.shape = size
data.cpt = size[0]
res.AddContainer(data)
return res
[docs] def GetDimensions(self) -> np.ndarray:
if self.Dimensions is None:
return self.dataitems[0].Dimensions
return self.Dimensions[::-1]
[docs]class XdmfGeometry(XdmfBase):
"""XdmfGeometry class: stores the point positions """
def __init__(self):
self.dataitems = []
self.Type = None
def __str__(self):
res = TFormat.GetIndent() + 'XdmfGeometry\n'
res += self._XdmfBase__AttributesToStr()
return res
[docs] def ReadAttributes(self, attrs):
if 'Type' in attrs:
self.Type = self.ReadAttribute(attrs, 'Type')
else:
self.Type = self.ReadAttribute(attrs, 'GeometryType')
[docs] def GetOrigin(self) -> np.ndarray:
if self.Type == "ORIGIN_DXDYDZ":
# self.Read()
return self.dataitems[0].GetData()[::-1]
else:
raise Exception # pragma: no cover
[docs] def GetSpacing(self) -> np.ndarray:
if self.Type == "ORIGIN_DXDYDZ":
# self.Read()
return self.dataitems[1].GetData()[::-1]
else:
raise Exception # pragma: no cover
[docs] def GetNodes(self) -> np.ndarray:
if self.Type == "XYZ" or self.Type == "XY":
return self.dataitems[0].GetData()
else:
raise Exception # pragma: no cover
[docs]class XdmfAttribute(XdmfBase):
""" XdmfAttribute class: to store the data over the grids """
def __init__(self):
self.dataitems = []
self.Name = ''
self.Type = ''
self.Center = ''
self.CDATA = ''
[docs] def ReadAttributes(self, attrs):
self.Name = self.ReadAttribute(attrs, 'Name')
try:
self.Type = self.ReadAttribute(attrs, 'Type')
except:
self.Type = self.ReadAttribute(attrs, 'AttributeType')
self.Center = self.ReadAttribute(attrs, 'Center')
def __str__(self):
res = TFormat.GetIndent() + 'XdmfAttribute\n'
TFormat.II()
res += TFormat.GetIndent() + 'Name : ' + self.Name + '\n'
res += TFormat.GetIndent() + 'Type : ' + self.Type + '\n'
res += TFormat.GetIndent() + 'Center : ' + self.Center + '\n'
for d in self.dataitems:
res += d.__str__()
TFormat.DI()
return res
[docs]class XdmfTime(XdmfBase):
""" XdmfTime class: to store the data over the grids """
def __init__(self):
self.Value = None
[docs] def ReadAttributes(self, attrs):
self.Value = np.array(self.ReadAttribute(attrs, 'Value').split(), dtype=MuscatFloat)
def __str__(self):
res = TFormat.GetIndent() + 'XdmfTime'
TFormat.II()
res += TFormat.GetIndent() + 'Value : ' + str(self.Value) + '\n'
TFormat.DI()
return res
[docs]class XdmfDataItem(XdmfBase):
""" XdmfDataItem class : class to manage the reading of the data Heavy and light """
def __init__(self):
self.Type = None
self.Dimensions = []
self.Precision = None
self.Format = None
self.Data = []
self.CDATA = ''
self.Seek = 0
self.Endian = "Native"
[docs] def ReadAttributes(self, attrs, path):
self.path = path
self.Dimensions = np.array(self.ReadAttribute(attrs, 'Dimensions').split(), dtype='int')
if 'DataType' in attrs:
self.Type = self.ReadAttribute(attrs, 'DataType')
elif 'NumberType' in attrs:
self.Type = self.ReadAttribute(attrs, 'NumberType')
else:
self.Type = 'float'
if ((self.Type.lower() == 'float' or self.Type.lower() == 'int') and 'Precision' in attrs):
self.Precision = int(self.ReadAttribute(attrs, 'Precision'))
self.Format = self.ReadAttribute(attrs, 'Format')
self.Seek = int(self.ReadAttribute(attrs, 'Seek', 0))
self.Endian = self.ReadAttribute(attrs, 'Endian', "Native")
def __str__(self):
res = TFormat.GetIndent() + 'XdmfDataItem \n'
TFormat.II()
res += TFormat.GetIndent() + 'Type : ' + self.Type + '\n'
res += TFormat.GetIndent() + 'Dimensions : ' + str(self.Dimensions) + '\n'
if self.Type.lower() == 'float':
res += TFormat.GetIndent() + 'Precision : ' + str(self.Precision) + '\n'
elif self.Type.lower() == 'int':
res += TFormat.GetIndent() + 'Precision : ' + str(self.Precision) + '\n'
res += TFormat.GetIndent() + 'Format : ' + str(self.Format) + '\n'
res += TFormat.GetIndent() + 'Data : \n ' + TFormat.GetIndent() + str(self.Data) + '\n'
res += TFormat.GetIndent() + 'CDATA : \n ' + TFormat.GetIndent() + str(self.CDATA) + '\n'
TFormat.DI()
return res
[docs] def TreatCDATA(self):
if len(self.CDATA) == 0:
return
if self.Format == 'XML':
if self.Type.lower() == 'float':
if self.Precision == 4:
numpytype = 'float32'
else:
numpytype = 'float_'
elif self.Type.lower() == 'int':
numpytype = 'int_'
elif self.Type.lower() == 'char':
numpytype = 'int8'
self.Data = np.array(self.CDATA.split(), dtype=numpytype)
self.Data = self.Data.reshape(self.Dimensions)
self.CDATA = ''
elif self.Format == 'HDF':
filename, dataSetPath = str(self.CDATA).strip().split(":",1)
# print(dataSetPath)
from h5py import File as __File
f = __File(os.path.join(self.path, filename), 'r')
# print(f[dataSetPath])
self.Data = np.array(f[dataSetPath])
self.Data = self.Data.reshape(self.Dimensions)
# print(self.Data)
self.CDATA = ''
elif self.Format == 'Binary':
if self.Type.lower() == 'float':
if self.Precision == 4:
numpytype = 'float32'
else:
numpytype = 'float_'
elif self.Type.lower() == 'int':
if self.Precision == 4:
numpytype = np.int32
elif self.Precision == 8:
numpytype = np.int64
else:
raise Exception(f"Dont know how to treat this type of Precision: {self.Precision}")
elif self.Type.lower() == 'char':
numpytype = 'int8'
binfilename = str(self.CDATA).strip()
binfile = open(os.path.join(self.path, binfilename), "rb")
binfile.seek(self.Seek)
self.Data = np.fromfile(binfile, dtype=numpytype, count=np.prod(self.Dimensions))
if self.Endian == "Native":
pass
elif self.Endian == "Big" and sys.byteorder == "little":
self.Data.byteswap(inplace=True)
elif self.Endian == "Little" and sys.byteorder == "big":
self.Data.byteswap(inplace=True)
self.Data.shape = self.Dimensions
binfile.close()
self.CDATA = ''
else:
raise Exception("Heavy data in format '" + self.Format + "' not suported yet") # pragma: no cover
[docs] def GetData(self):
self.TreatCDATA()
return self.Data
[docs]class XdmfReader(xml.sax.ContentHandler):
"""Xdmf Reader class
"""
def __init__(self, filename=''):
super().__init__()
self.xdmf = Xdmf()
self.pile = []
self.path = ''
self.filename = ''
self.SetFileName(filename)
self.readed = False
self.lazy = True
self.canHandleTemporal = True
self.time = np.array([])
self.timeToRead = -1
self.encoding = None
[docs] def Reset(self):
"""Resets the reader
"""
self.xdmf = Xdmf()
self.pile = []
self.readed = False
self.timeToRead = -1
self.time = np.array([])
[docs] def SetFileName(self, filename):
"""Sets the name of file to read
Parameters
----------
filename : str
file name to set
"""
# if same filename no need to read the file again
if self.filename == filename:
return
self.readed = False
self.filename = filename
self.path = os.path.dirname(filename)
[docs] def GetAvailableTimes(self):
"""Returns the available times at which data can be read
Returns
-------
np.ndarray
available times at which data can be read
"""
return self.time
[docs] def SetTimeToRead(self, time=None, timeIndex=None):
"""Sets the time at which the data is read
Parameters
----------
time : float, optional
time at which the data is read, by default None
timeIndex : int, optional
time index at which the data is read, by default None
Returns
-------
int
time index at which the data is read
"""
if time is not None:
self.timeToRead = time
else:
if timeIndex is not None:
self.timeToRead = self.time[timeIndex]
else:
self.timeToRead = -1
[docs] def Read(self, fileName=None):
"""Function that performs the reading of the data defined in an ut file
Parameters
----------
fileName : str, optional
name of the file to read, by default None
Returns
-------
Mesh or ConstantRectilinearMesh
output mesh object containing reading result
"""
if self.timeToRead == -1.:
timeIndex = len(self.time)-1
else:
if len(self.time) == 0:
timeIndex = -1
else:
timeIndex = np.argmin(abs(self.time - self.timeToRead))
if fileName is not None:
self.SetFileName(fileName)
self.readed = False
if len(self.filename) == 0:
raise Exception('Need a filename ')
# read only one time the filel
if self.readed:
return
self.Reset()
thefile = open(self.filename, "r")
# to deactivate the DTD validation
parser = xml.sax.make_parser()
parser.setContentHandler(self)
parser.setFeature(xml.sax.handler.feature_external_ges, False)
parser.parse(thefile)
thefile.close()
return self.xdmf.GetDomain(-1).GetGrid(timeIndex).GetSupport()
# this a a overloaded function (must start with lower case) !!!!!
[docs] def startElement(self, name, attrs):
if name == "Xdmf":
self.pile.append(self.xdmf)
return
father = self.pile[-1]
if name == "Domain":
res = XdmfDomain()
res.ReadAttributes(attrs)
father.domains.append(res)
elif name == 'Grid':
res = XdmfGrid()
res.ReadAttributes(attrs)
# for the moment we use a flat representation (no GridType collection)
if "GridType" in attrs and attrs.get("GridType").lower() == "collection":
res = father
else:
father.grids.append(res)
elif name == 'Information':
res = XdmfInformation()
res.ReadAttributes(attrs)
father.information.append(res)
elif name == 'Topology':
res = XdmfTopology()
res.ReadAttributes(attrs)
father.topology = res
elif name == 'Geometry':
res = XdmfGeometry()
res.ReadAttributes(attrs)
father.geometry = res
elif name == 'DataItem':
res = XdmfDataItem()
res.ReadAttributes(attrs, self.path)
father.dataitems.append(res)
elif name == 'Attribute':
res = XdmfAttribute()
res.ReadAttributes(attrs)
father.attributes.append(res)
elif name == "Time":
res = XdmfTime()
res.ReadAttributes(attrs)
father.time = res
else:
raise Exception("Unknown tag : '" + name + "' Sorry!") # pragma: no cover
self.pile.append(res)
# this a a overloaded function (must start with lower case) !!!!!
[docs] def characters(self, content):
#print("*** " + content + " +++ ")
father = self.pile[-1]
if isinstance(father, XdmfDataItem):
father.CDATA += content # Note: '+=', not '='
# this a a overloaded function (must start with lower case) !!!!!
[docs] def endElement(self, name):
if self.lazy:
self.pile.pop()
else:
self.pile.pop().TreatCDATA()
def __str__(self):
res = ''
for d in self.xdmf.domains:
res = d.__str__() + "\n"
return res
RegisterReaderClass(".xdmf", XdmfReader)
RegisterReaderClass(".pxdmf", XdmfReader)
[docs]def CheckIntegrity():
# All the try pasrt are to execute the code for error handling.
# The code in the try(s) must raise exceptions
# TODO Verified if the exceptions are the good ones
# Get TestData directory
from Muscat.TestData import GetTestDataPath
TestDataPath = GetTestDataPath()
# Create a Reader
res = XdmfReader()
try:
res.Read()
return 'Not OK' # pragma: no cover
except Exception as e:
pass
for filename in ["Unstructured.xmf", "UnstructuredBinary.xmf", "UnstructuredAscii.xmf"]:
if not HasHdf5Support() and filename == "Unstructured.xmf":
continue
res = XdmfReader(filename=TestDataPath + filename)
res.lazy = False
res.Read()
res = XdmfReader(filename=TestDataPath + filename)
# read only the xml part
res.Read()
rep = res.__str__()
# print(rep)
# Test xdmf part***********************
a = res.xdmf.__str__()
try:
res.xdmf.GetDomain("Imaginary Domain")
return 'Not OK' # pragma: no cover
except Exception as e:
pass
res.xdmf.GetDomain("Dom 1")
domain = res.xdmf.GetDomain(0)
# Test domain part***********************
try:
domain.GetGrid("imaginary grid")
return 'Not OK' # pragma: no cover
except Exception as e:
pass
domain.GetGrid(0)
grid = domain.GetGrid("Grid")
# Test Grid part**************************************************************
names = grid.GetFieldsNames()
if names[0] != 'RTData':
raise Exception()
grid.HasField(names[0])
grid.HasField('toto')
try:
data = res.xdmf.GetDomain("Dom 1").GetGrid("Grid").GetFieldData('ImaginaryField')
return 'Not OK' # pragma: no cover
except Exception as e:
pass
grid.GetFieldTermsAsColMatrix('term_')
try:
grid.GetFieldTermsAsColMatrix('ImaginaryField_')
return 'Not OK' # pragma: no cover
except Exception as e:
pass
grid.GetFieldTermsAsTensor('TensorField_', offsetj=1)
try:
grid.GetFieldTermsAsTensor('ImaginaryField_', offsetj=1)
return 'Not OK' # pragma: no cover
except Exception as e:
pass
grid.GetFieldData('IntField')
try:
grid.GetFieldData('UnknownField')
return 'Not OK' # pragma: no cover
except Exception as e:
pass
data = res.xdmf.domains[0].GetGrid("Grid").GetFieldData('RTData')
if data[49] != 260.0:
raise
geo = res.xdmf.domains[0].GetGrid("Grid").geometry.dataitems[0].GetData()
# print(geo)
if geo[0, 2] != -1:
raise
topo = res.xdmf.domains[0].GetGrid("Grid").topology.dataitems[0].GetData()
# print(topo)
if topo[0] != XdmfNumber[ED.Hexahedron_8]:
raise
######################### Structured #########################
res = XdmfReader(filename=TestDataPath + "Structured.xmf")
# read only the xml part
res.Read()
domain = res.xdmf.GetDomain(0)
# Test domain part***********************
grid = domain.GetGrid(0)
# print(grid.GetFieldsOfType("Node"))
# print("--")
# print(grid.GetFieldsOfType("Cell"))
grid.geometry.GetOrigin()
grid.geometry.GetSpacing()
##################################
Example1()
Example2()
return 'OK'
[docs]def Example1():
import Muscat.TestData as test
# Create a Reader
reader = XdmfReader(filename=test.GetTestDataPath() + "UnstructuredBinary.xmf")
# Do the reading (only the xml part, to read all the data set lazy to False)
#res.lazy = False
reader.Read()
# Get the domaine "Dom 1"
dom = reader.xdmf.GetDomain("Dom 1")
# Get the first Grid
grid = dom.GetGrid(0)
grid.topology.GetDimensions()
names = grid.GetPointFieldsNames()
allFields = grid.GetPointFields()
if len(allFields):
print(grid.GetPointField(names[0]))
names = grid.GetCellFieldsNames()
allFields = grid.GetCellFields()
if len(allFields):
print(grid.GetCellField(names[0]))
names = grid.GetGridFieldsNames()
allFields = grid.GetGridFields()
if len(allFields):
print(grid.GetGridField(names[0]))
# Get one field (or term)
dataField1 = grid.GetFieldData('RTData')
dataField1 = grid.GetFieldsOfType('Node')
# print(dataField1.shape)
# Get all the term as a matrix
dataField2 = grid.GetFieldTermsAsColMatrix('term_')
# print(dataField2.shape)
# Get all the term as a matrix
dataField3 = grid.GetFieldTermsAsTensor('TensorField_', offsetj=1)
# print(dataField3.shape)
[docs]def Example2():
import Muscat.TestData as test
reader = XdmfReader(filename=test.GetTestDataPath() + "TensorTestData.xmf")
# Do the reading (only the xml part, to read all the data set lazy to False)
#res.lazy = False
reader.Read()
# Get the domaine "Dom 1"
dom = reader.xdmf.GetDomain(0)
if __name__ == '__main__':
print(CheckIntegrity()) # pragma: no cover
import Muscat.Helpers.Logger as Logger
Logger.silent = "OK"