# -*- 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.
#
"""
rst based ansys mesh reader.
returns an unstructured mesh (class from BasicTools)
"""
from typing import List
import numpy as np
from Muscat.IO.ReaderBase import ReaderBase
from Muscat.IO.IOFactory import RegisterReaderClass
from Muscat.Types import MuscatFloat, MuscatIndex
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.Helpers.Logger import Debug
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.IO.AnsysTools import nbIntegrationsPoints, AnsysElementDescriptorToMuscatElementType
from Muscat.FE.Fields.FEField import FEField
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
[docs]
def ReadRst(fileName: str, time: MuscatFloat = None, timeIndex: MuscatIndex = None):
"""Function mesh and solution fields (nodal and elemental) from a rst file
Parameters
----------
fileName : str, optional
name of the file to be read, by default None
time : float, optional
time at which the field is read, by default None
timeIndex :int, optional
time index at which the field is read, by default None
Returns
-------
Mesh
output Mesh object containing the results
"""
reader = RstReader()
reader.SetFileName(fileName)
reader.ReadMetaData()
return reader.Read(time=time, timeIndex=timeIndex)
[docs]
class RstReader(ReaderBase):
"""Rst Reader class
"""
def __init__(self):
""" Initialize the RstReader
Raises:
ModuleNotFoundError: if the import ansys.dpf.core could not be load
"""
super().__init__()
self.Reset()
self.canHandleTemporal = True
try:
import ansys.dpf.core
except ImportError:
raise ModuleNotFoundError("To use this module ansys.pydpf.core must be installed")
[docs]
def SetFileName(self, fileName: str) -> None:
"""Function to set the file to read
If the simulation is modal then it sets self.modal to True because it needs
to accomodate to the 1 index of modal simulations
Parameters
----------
fileName : str
path to the simulation
"""
if self.fileName != fileName:
self.Reset()
from ansys.dpf import post
self.modal = int(isinstance(post.load_simulation(fileName), post.ModalMechanicalSimulation))
self.fileName = fileName
[docs]
def Reset(self):
"""Resets the reader
"""
self.fileName = None
self.model = None
self.meshMetadata = None
self.nodal = []
self.elemental_nodal = []
self.elementary = []
self.atIntegrationPoints = False
self.IdsToNodes = None
self.timeToRead = -1
self.output = None
[docs]
def GetModel(self):
"""Function to recover the dpf Model associated to the rst file
Raises:
RuntimeError: if the fileName is not already set
Returns:
Model: dpf model
"""
if self.model is None:
if self.fileName is None:
raise RuntimeError("Need to set the fileName first")
from ansys.dpf.core import Model
self.model = Model(self.fileName)
return self.model
[docs]
def MeshRead(self) -> Mesh:
"""Function that performs the reading of a mesh in a rst file
Returns
-------
Mesh
mesh containing reading result
"""
res = Mesh()
rst = self.GetModel()
metadata = rst.metadata
ansys_mesh = metadata.meshed_region
oidToElementContainerAndIndex = {}
nbNodes = ansys_mesh.nodes.n_nodes
if nbNodes:
# print(' ------------ Reading the NODES ------------')
res.originalIDNodes = np.asarray(ansys_mesh.nodes.scoping.get_ids(), dtype=MuscatIndex)
self.IdsToNodes = np.empty(np.max(res.originalIDNodes)+1, dtype=MuscatIndex)
self.IdsToNodes[res.originalIDNodes] = np.arange(res.originalIDNodes.shape[0])
res.nodes = np.asarray(ansys_mesh.nodes.coordinates_field.data, dtype=MuscatFloat)
# ------------ Reading the GROUPS/SELECTIONS ------------
group_names = ansys_mesh.available_named_selections
for group_name in group_names:
group_container = ansys_mesh.named_selection(group_name)
if group_container.location == 'Nodal':
nsetname = group_name
if len(nsetname[0]) and group_container.size:
tag = res.GetNodalTag(nsetname)
tag.SetIds(np.asarray(group_container.ids, dtype=MuscatIndex)-1)
nbElements = ansys_mesh.elements.n_elements
if nbElements:
j = 0
element_types = ansys_mesh.elements.element_types_field.data
used_types, return_inverse, return_counts = np.unique(element_types, return_inverse=True, return_counts=True)
connectivity_buffer = ansys_mesh.elements.connectivities_field.data
element_ids = ansys_mesh.elements.scoping.ids
nbnodes_per_type = np.array([ED.numberOfNodes[AnsysElementDescriptorToMuscatElementType[used_type]] for used_type in used_types])
offsets = np.empty(nbElements+1, dtype=MuscatIndex)
offsets[0] = 0
offsets[1:] = np.add.accumulate(nbnodes_per_type[return_inverse])
for i, used_type in enumerate(used_types):
nameType = AnsysElementDescriptorToMuscatElementType[used_type]
elements = res.GetElementsOfType(nameType)
cpt0 = elements.GetNumberOfElements()
elements.Reserve(cpt0+return_counts[i])
numberOfNodesPerElement = elements.GetNumberOfNodesPerElement()
index = np.where(element_types == used_type)[0]
mask = np.arange(numberOfNodesPerElement, dtype=MuscatIndex) * np.ones((len(index), numberOfNodesPerElement), dtype=MuscatIndex) + offsets[index, None]
localConnectivity = connectivity_buffer[mask]
local_ids = element_ids[index]
elements.AddNewElements(localConnectivity, local_ids)
for j, oid in enumerate(local_ids):
oidToElementContainerAndIndex[oid] = (elements, j+cpt0)
eolOriginalIds = res.GetElementsOriginalIDs()
self.IdsToElements = np.empty(np.max(eolOriginalIds)+1, dtype=MuscatIndex)
self.IdsToElements[eolOriginalIds] = np.arange(eolOriginalIds.shape[0])
# ------------ Reading the GROUPS/SELECTIONS ------------
group_names = ansys_mesh.available_named_selections
for group_name in group_names:
group_container = ansys_mesh.named_selection(group_name)
if group_container.location == 'Elemental':
if group_container.size != 0:
# res.AddElementToTagUsingOriginalId(group_container.ids, group_name)
for oid_elem in group_container.ids:
elementContainer, index = oidToElementContainerAndIndex[oid_elem]
elementContainer.tags.CreateTag(group_name, False).AddToTag(index)
res.PrepareForOutput()
self.output = res
return res
[docs]
def Read(self, time: MuscatFloat = None, timeIndex: MuscatIndex = None) -> Mesh:
"""Function that performs the reading of the mesh and the nodal/element fields in the rst
Parameters
----------
time : float, optional
time at which the field is read, by default None
timeIndex :int, optional
time index at which the field is read, by default None
Returns
-------
Mesh
output Mesh object containing the results
"""
mesh = self.MeshRead()
for nodefield in self.nodal:
mesh.nodeFields[nodefield] = self.ReadField(fieldname=nodefield, time=time, timeIndex=timeIndex)
for elemField in self.elementary:
mesh.elemFields[elemField] = self.ReadField(fieldname=elemField, time=time, timeIndex=timeIndex)
return mesh
[docs]
def ReadField(self, fieldname, time: MuscatFloat = None, timeIndex: MuscatIndex = None) -> np.ndarray:
"""Function that performs the reading of a field
Parameters
----------
fieldname : str, optional
name of the field to be read, by default None
time : float, optional
time at which the field is read, by default None
timeIndex : int, optional
time index at which the field is read, by default None
Returns
-------
np.ndarray or FEField
Depending on the location of the field in the rst this function can return
a 'Nodal' field (numpy array), "Elemental" field (numpy array), or a "ElementalNodal" field (FEField)
"""
self.ReadMetaData()
# need to shift the index in case of a modal study
timeIndex = self.SetTimeToRead(time, timeIndex) + self.modal
Debug("Reading timeIndex : " + str(timeIndex))
rst = self.GetModel()
results_ansys = rst.results
if fieldname in self.nodal:
scoping = rst.metadata.meshed_region.nodes.scoping
location = 'Nodal'
elif fieldname in self.elementary:
scoping = rst.metadata.meshed_region.elements.scoping
location = 'Elemental'
elif fieldname in self.elemental_nodal:
elementFilterSolid = ElementFilter(dimensionality=3)
solid_ids = []
for selection in elementFilterSolid(self.output):
solid_ids.extend(selection.elements.originalIds)
from ansys.dpf.core import mesh_scoping_factory
solid_ids = np.array(solid_ids, dtype=MuscatIndex)
scoping = mesh_scoping_factory.elemental_scoping(solid_ids)
location = 'ElementalNodal'
else:
raise(Exception("unable to find field " +str(fieldname) ))
avail_res = rst.metadata.result_info.available_results
# Extract the nodal data from the RST
name = None
for nm in avail_res:
if nm.native_location == location:
for component, sub_r in enumerate(nm.sub_results):
if fieldname == sub_r['operator name']:
name = nm.name
break
else:
if fieldname == nm.operator_name:
name = nm.name
component = None
if name is not None:
break
field_in_time = getattr(results_ansys, name)(mesh_scoping=scoping, time_scoping=timeIndex).outputs.fields_container()
data = np.array(field_in_time[0].data)
support_ids = np.asarray(field_in_time[0].scoping.get_ids())
if location == 'ElementalNodal':
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
# generate a discontinues space to hold the field "nodal non moyenné"
space = LagrangeSpaceP1
numbering = ComputeDofNumbering(self.output, space=space, elementFilter=elementFilterSolid, discontinuous=True)
res = FEField(name=name, mesh=self.output, space=space, numbering=numbering)
res.Allocate()
if component is not None:
data = data[:, component]
elementsOriginalIDs = self.output.GetElementsOriginalIDs()
nodesPerElements = np.zeros((np.max(elementsOriginalIDs)+1), dtype=MuscatIndex)
elementsTypes = np.zeros((np.max(elementsOriginalIDs)+1), dtype=MuscatIndex)-1
for selection in elementFilterSolid.IterOnTypesOnly(self.output):
nodesPerElements[selection.elements.originalIds] = space[selection.elementType].GetNumberOfShapeFunctions()
elementsTypes[selection.elements.originalIds] = selection.elementType.value
originalIdsToIndexGlobal = np.empty(np.max(elementsOriginalIDs)+1, dtype=MuscatIndex)
originalIdsToIndexGlobal[elementsOriginalIDs] = np.arange(len(elementsOriginalIDs), dtype=MuscatIndex)
nodesPerElements_extraction = nodesPerElements[support_ids]
elementsTypes_extraction = elementsTypes[support_ids]
offsets = np.empty(len(nodesPerElements_extraction)+1, dtype=MuscatIndex)
offsets[0] = 0
offsets[1:] = np.add.accumulate(nodesPerElements_extraction)
cpt = 0
off = 0
for selection in elementFilterSolid(self.output):
index = np.where(elementsTypes_extraction == selection.elementType.value)[0]
numberOfNodesPerElement = space[selection.elementType].GetNumberOfShapeFunctions()
fillingMask = np.arange(numberOfNodesPerElement, dtype=MuscatIndex) * np.ones((len(index), numberOfNodesPerElement), dtype=MuscatIndex)
mask = fillingMask + offsets[index, None]
localData = data[mask]
newData = np.reshape(localData, (selection.elements.GetNumberOfElements(), numberOfNodesPerElement))
targetIndices = (originalIdsToIndexGlobal[support_ids[index]] - cpt)*numberOfNodesPerElement + off
newMask = fillingMask + targetIndices[:, None]
res.data[newMask.ravel()] = newData.ravel('f')
off += selection.elements.GetNumberOfElements()*numberOfNodesPerElement
cpt += selection.elements.GetNumberOfElements()
return res
else:
if location == 'Nodal':
nbentities = self.output.GetNumberOfNodes()
mapping = self.IdsToNodes
elif location == 'Elemental':
nbentities = self.output.GetNumberOfElements()
mapping = self.IdsToElements
if component is None:
res = np.zeros((nbentities, data.shape[1]))
res[mapping[support_ids], :] = np.asarray(data)
else:
res = np.zeros((nbentities, 1))
res[mapping[support_ids], 0] = np.asarray(data[:, component])
return res
[docs]
def SetTimeToRead(self, time: MuscatFloat = None, timeIndex: MuscatIndex = None) -> MuscatIndex:
"""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 None) and (timeIndex is None):
if self.timeToRead == -1:
self.timeToRead = self.time[-1]
return len(self.time)-1
else:
return np.where(self.time == self.timeToRead)[0]
if (time is not None) and (timeIndex is not None):
raise (Exception("Cannot specify both time and timeIndex"))
if time is None:
self.timeToRead = self.time[timeIndex]
return timeIndex
elif time == -1:
self.timeToRead = self.time[-1]
return timeIndex
else:
self.timeToRead = time
return np.where(self.time == self.timeToRead)[0]
[docs]
def GetAvailableTimes(self) -> np.ndarray:
"""Returns the available times at which data can be read
Returns
-------
np.ndarray
available times at which data can be read
"""
self.ReadMetaData()
return self.time
[docs]
def GetAvailableResultsAndOperators(self) -> List:
"""Returns the available results and its operator name
Returns
-------
list
contains a namedTuple with the results. attributes are name, operator_name, native_location, n_components
"""
rst = self.GetModel()
class Results():
__slots__ = ("name", "operator_name", "native_location", "n_components")
def __init__(self, name, operator_name, native_location, n_components):
self.name = name
self.operator_name = operator_name
self.native_location = native_location
self.n_components = n_components
def __repr__(self):
return f"{self.name}({self.operator_name}, {self.native_location}, {self.n_components})"
return [Results(r.name, r.operator_name, r.native_location, r.n_components) for r in rst.metadata.result_info.available_results]
RegisterReaderClass(".rst", RstReader)
[docs]
def CheckIntegrity(GUI: bool = False):
try:
import ansys.dpf.core as dpf
except ImportError:
return 'skip: ansys.dpf.core not available '
return "ok"