# -*- 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 and result reader.
This module provides the :class:`RstReader` class and the :func:`ReadRst`
helper function to read meshes and solution fields (nodal, elemental and
elemental-nodal) from an Ansys ``.rst`` result file. The reading relies on the
``ansys-dpf-core`` and ``ansys-dpf-post`` packages and returns a Muscat
:class:`~Muscat.MeshContainers.Mesh.Mesh` instance.
"""
from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union
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, Warning
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
if TYPE_CHECKING: # pragma: no cover
from ansys.dpf.core import Model
#: Error message raised when the optional Ansys DPF dependencies are missing.
_ANSYS_IMPORT_ERROR = "To use this module ansys-dpf-core and ansys/pydpf-post must be installed"
# _APDL_CONTACT_TYPES = frozenset((169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179))
_APDL_CONTA_TYPES = frozenset((171, 172, 173, 174, 175, 176, 177, 178))
_APDL_TARGET_TYPES = frozenset((169, 170, 179))
[docs]
class Results:
"""Lightweight description of an available result of an rst file.
Attributes
----------
name : str
Name of the result.
operator_name : str
Name of the DPF operator producing the result.
native_location : str
Native location of the result (``"Nodal"``, ``"Elemental"`` or
``"ElementalNodal"``).
n_components : int
Number of components of the result.
"""
__slots__ = ("name", "operator_name", "native_location", "n_components")
def __init__(self, name: str, operator_name: str, native_location: str, n_components: int) -> None:
self.name = name
self.operator_name = operator_name
self.native_location = native_location
self.n_components = n_components
def __repr__(self) -> str:
return f"{self.name}({self.operator_name}, {self.native_location}, {self.n_components})"
[docs]
def ReadRst(fileName: str, time: Optional[MuscatFloat] = None, timeIndex: Optional[MuscatIndex] = None, fromElementalNodalToNodal: bool = True, dimensionality: Optional[Union[int, List[int]]] = None) -> Mesh:
"""Read the mesh and solution fields (nodal, elemental and elemental-nodal) from a rst file.
Parameters
----------
fileName : str
Name of the file to be read.
time : MuscatFloat, optional
Time at which the fields are read, by default None.
timeIndex : MuscatIndex, optional
Time index at which the fields are read, by default None.
fromElementalNodalToNodal : bool, optional
If True, elemental-nodal fields are averaged at the nodes, by default True.
Returns
-------
Mesh
Output Mesh object containing the results.
"""
reader = RstReader()
reader.SetFileName(fileName)
reader.ReadMetaData()
return reader.Read(time=time, timeIndex=timeIndex, fromElementalNodalToNodal=fromElementalNodalToNodal, dimensionality=dimensionality)
[docs]
class RstReader(ReaderBase):
"""Reader for Ansys ``.rst`` result files.
This class reads the mesh and the available solution fields stored in an
Ansys ``.rst`` file, supporting nodal, elemental and elemental-nodal fields
as well as temporal and modal simulations.
"""
def __init__(self) -> None:
"""Initialize the RstReader."""
super().__init__()
self.Reset()
self.canHandleTemporal = True
self.modal = False
@staticmethod
def _ImportDpf():
"""Import and return the ``ansys.dpf.core`` module.
Returns
-------
module
The imported ``ansys.dpf.core`` module.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-core`` and ``ansys-dpf-post`` are not installed.
"""
try:
import ansys.dpf.core as dpf
except ImportError:
raise ModuleNotFoundError(_ANSYS_IMPORT_ERROR)
return dpf
@staticmethod
def _ImportPost():
"""Import and return the ``ansys.dpf.post`` module.
Returns
-------
module
The imported ``ansys.dpf.post`` module.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-post`` is not installed.
"""
try:
from ansys.dpf import post
except ImportError:
raise ModuleNotFoundError(_ANSYS_IMPORT_ERROR)
return post
[docs]
def SetFileName(self, fileName: str) -> None:
"""Set the file to read.
Parameters
----------
fileName : str
Path to the simulation file.
"""
if self.fileName != fileName:
self.Reset()
self.fileName = fileName
[docs]
def Reset(self) -> None:
"""Reset the reader to its initial state."""
self.fileName = None
self.model = None
self.meshMetadata = None
self.nodal = {}
self.elementary = {}
self.elemental_nodal = {}
self.atIntegrationPoints = False
self.IdsToNodes = None
self.IdsToElements = None
self.timeToRead = -1
self.output = None
self.to_nodal = None
self.contactElementIds = np.empty(0, dtype=MuscatIndex)
self.targetElementIds = np.empty(0, dtype=MuscatIndex)
[docs]
def GetModel(self) -> "Model":
"""Recover the dpf Model associated with the rst file.
Returns
-------
Model
The dpf model associated with the file.
Raises
------
RuntimeError
If the fileName has not been set yet.
ModuleNotFoundError
If ``ansys-dpf-core`` and ``ansys-dpf-post`` are not installed.
"""
if self.model is None:
if self.fileName is None:
raise RuntimeError("Need to set the fileName first")
dpf = self._ImportDpf()
self.model = dpf.Model(self.fileName)
return self.model
[docs]
def MeshRead(self) -> Mesh:
"""Read the mesh stored in the rst file.
Returns
-------
Mesh
Mesh containing the reading result.
"""
res = Mesh()
rst = self.GetModel()
metadata = rst.metadata
ansys_mesh = metadata.meshed_region
oidToElementContainerAndIndex = {}
nbNodes = ansys_mesh.nodes.n_nodes
if nbNodes:
# ------------ 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 NODAL GROUPS/SELECTIONS ------------
for group_name in ansys_mesh.available_named_selections:
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:
element_types = ansys_mesh.elements.element_types_field.data
used_types, return_inverse = np.unique(element_types, return_inverse=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 used_type in used_types:
nameType = AnsysElementDescriptorToMuscatElementType[used_type]
numberOfNodesPerElement = ED.numberOfNodes[nameType]
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 = np.asarray(connectivity_buffer[mask], dtype=MuscatIndex, order="C")
local_ids = np.asarray(element_ids[index], dtype=MuscatIndex)
self._AddElementsFromAnsysConnectivity(
res, nameType, localConnectivity, local_ids, oidToElementContainerAndIndex
)
eolOriginalIds = res.GetElementsOriginalIDs()
self.IdsToElements = np.empty(np.max(eolOriginalIds)+1, dtype=MuscatIndex)
self.IdsToElements[eolOriginalIds] = np.arange(eolOriginalIds.shape[0])
# ------------ Reading the ELEMENTAL GROUPS/SELECTIONS ------------
for group_name in ansys_mesh.available_named_selections:
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
try:
apdlTypes = np.asarray(ansys_mesh.property_field("apdl_element_type").data)
contactMask = np.isin(apdlTypes, list(_APDL_CONTA_TYPES))
self.contactElementIds = np.asarray(element_ids, dtype=MuscatIndex)[contactMask]
targetMask = np.isin(apdlTypes, list(_APDL_TARGET_TYPES))
self.targetElementIds = np.asarray(element_ids, dtype=MuscatIndex)[targetMask]
except Exception:
self.contactElementIds = np.empty(0, dtype=MuscatIndex)
self.targetElementIds = np.empty(0, dtype=MuscatIndex)
return res
def _AddElementsFromAnsysConnectivity(
self,
res: Mesh,
elementType: ED.ElementType,
connectivity: np.ndarray,
originalIds: np.ndarray,
oidToElementContainerAndIndex: Dict,
) -> None:
"""Add Ansys/DPF connectivity to a Muscat mesh, handling dropped midside nodes.
Parameters
----------
res : Mesh
Mesh to which the elements are added.
elementType : ED.ElementType
Muscat element type of the elements to add.
connectivity : np.ndarray
Connectivity array of the elements. Negative entries denote dropped
midside nodes.
originalIds : np.ndarray
Original Ansys identifiers of the elements.
oidToElementContainerAndIndex : dict
Mapping (updated in place) from an element original id to the tuple
(element container, local index).
Raises
------
RuntimeError
If an element has a negative connectivity that cannot be completed.
"""
negativeConnectivity = np.any(connectivity < 0, axis=1)
if np.any(negativeConnectivity):
if elementType in self._GetMidNodeEdgesByElementType():
connectivity = self._CompleteMissingMidNodes(res, elementType, connectivity, originalIds)
self._AddElementsToContainer(res, elementType, connectivity, originalIds, oidToElementContainerAndIndex)
else:
badIndex = np.where(negativeConnectivity)[0][0]
badOriginalId = originalIds[badIndex]
badConnectivity = connectivity[badIndex]
raise RuntimeError(
f"RST element {badOriginalId} has negative connectivity for {elementType}: {badConnectivity}"
)
else:
self._AddElementsToContainer(
res, elementType, connectivity, originalIds, oidToElementContainerAndIndex
)
def _GetMidNodeEdgesByElementType(self) -> Dict[ED.ElementType, Dict[int, Tuple[int, int]]]:
"""Return, for each supported element type, the edge of each midside node.
Returns
-------
dict
Mapping from an element type to a dictionary that associates each
midside node position with the pair of corner node positions
defining the edge it lies on.
"""
return {
ED.Triangle_6: {
3: (0, 1),
4: (1, 2),
5: (2, 0),
},
ED.Quadrangle_8: {
4: (0, 1),
5: (1, 2),
6: (2, 3),
7: (3, 0),
},
ED.Tetrahedron_10: {
4: (0, 1),
5: (1, 2),
6: (2, 0),
7: (0, 3),
8: (1, 3),
9: (2, 3),
},
}
def _GetElementalNodalScopingIdsByDimensionality(self, fieldname: str, dimensionality: Optional[Union[int, List[int]]] = None) -> Dict[int, np.ndarray]:
"""Return the element ids per dimensionality (descending) for an elemental-nodal field.
Contact results (``contact_*``) are scoped on the contact elements only.
Structural results are grouped by dimensionality, highest first, excluding
contact and target elements. Reading dimensionalities separately (and
stopping at the first one carrying data) avoids the slow shell reads for
solid-carried fields while keeping shell-only results reachable.
"""
if fieldname.startswith("contact"):
if len(self.contactElementIds) == 0:
return {}
return {0: np.asarray(self.contactElementIds, dtype=MuscatIndex)}
dims = [3, 2, 1] if dimensionality is None else np.atleast_1d(dimensionality)[::-1].tolist()
excluded = np.concatenate((self.contactElementIds, self.targetElementIds))
res = {}
for dim in dims:
idArrays = [selection.elements.originalIds for selection in ElementFilter(dimensionality=int(dim))(self.output)]
if not idArrays:
continue
ids = np.concatenate(idArrays).astype(MuscatIndex)
if excluded.size:
ids = np.setdiff1d(ids, excluded.astype(MuscatIndex))
if ids.size:
res[int(dim)] = ids
return res
def _CompleteMissingMidNodes(
self,
res: Mesh,
elementType: ED.ElementType,
connectivity: np.ndarray,
originalIds: np.ndarray,
) -> np.ndarray:
"""Complete the missing midside nodes of a connectivity array.
Missing midside nodes (negative entries) are either retrieved from an
already known edge or created at the midpoint of the corresponding edge
and appended to the mesh nodes.
Parameters
----------
res : Mesh
Mesh providing the existing nodes and elements; new nodes are added
to it in place.
elementType : ED.ElementType
Muscat element type of the elements being completed.
connectivity : np.ndarray
Connectivity array possibly containing negative midside node entries.
originalIds : np.ndarray
Original Ansys identifiers of the elements.
Returns
-------
np.ndarray
Connectivity array with all midside nodes filled in.
Raises
------
RuntimeError
If a corner node connectivity is negative.
"""
resConnectivity = np.array(connectivity, dtype=MuscatIndex, order="C", copy=True)
edgeByMidNodeByElementType = self._GetMidNodeEdgesByElementType()
edgeByMidNode = edgeByMidNodeByElementType[elementType]
midNodePositions = np.array(list(edgeByMidNode), dtype=MuscatIndex)
numberOfCornerNodes = max(max(edge) for edge in edgeByMidNode.values())+1
invalidCorners = np.any(resConnectivity[:, :numberOfCornerNodes] < 0, axis=1)
if np.any(invalidCorners):
badIndex = np.where(invalidCorners)[0][0]
badOriginalId = originalIds[badIndex]
badConnectivity = connectivity[badIndex]
raise RuntimeError(
f"RST element {badOriginalId} has negative corner connectivity for {elementType}: {badConnectivity}"
)
missingElementIndices, missingMidPositionIndices = np.nonzero(resConnectivity[:, midNodePositions] < 0)
if len(missingElementIndices) == 0:
return resConnectivity
missingMidNodePositions = midNodePositions[missingMidPositionIndices]
missingEdges = self._ExtractEdges(resConnectivity, missingElementIndices, missingMidNodePositions, edgeByMidNode)
uniqueMissingEdges, uniqueMissingInverse = np.unique(missingEdges, axis=0, return_inverse=True)
knownEdges, knownMidNodes = self._CollectKnownMidNodes(
res, resConnectivity, edgeByMidNode, edgeByMidNodeByElementType
)
uniqueMissingMidNodes = np.empty(uniqueMissingEdges.shape[0], dtype=MuscatIndex)
newEdgeMask = np.ones(uniqueMissingEdges.shape[0], dtype=bool)
if len(knownEdges):
uniqueKnownEdges, uniqueKnownIndices = np.unique(knownEdges, axis=0, return_index=True)
uniqueKnownMidNodes = knownMidNodes[uniqueKnownIndices]
knownEdgesAsRecords = self._AsEdgeRecords(uniqueKnownEdges)
missingEdgesAsRecords = self._AsEdgeRecords(uniqueMissingEdges)
knownEdgeIndices = np.searchsorted(knownEdgesAsRecords, missingEdgesAsRecords)
knownEdgeMask = knownEdgeIndices < len(knownEdgesAsRecords)
knownEdgeMask[knownEdgeMask] = knownEdgesAsRecords[knownEdgeIndices[knownEdgeMask]] == missingEdgesAsRecords[knownEdgeMask]
uniqueMissingMidNodes[knownEdgeMask] = uniqueKnownMidNodes[knownEdgeIndices[knownEdgeMask]]
newEdgeMask[knownEdgeMask] = False
newEdges = uniqueMissingEdges[newEdgeMask]
if len(newEdges):
firstNewNodeIndex = res.nodes.shape[0]
newNodeIndices = np.arange(firstNewNodeIndex, firstNewNodeIndex+len(newEdges), dtype=MuscatIndex)
newNodes = np.mean(res.nodes[newEdges, :], axis=1)
res.nodes = np.vstack((res.nodes, newNodes))
nextOriginalNodeId = np.max(res.originalIDNodes)+1 if len(res.originalIDNodes) else 0
newOriginalIds = np.arange(nextOriginalNodeId, nextOriginalNodeId+len(newEdges), dtype=MuscatIndex)
res.originalIDNodes = np.hstack((res.originalIDNodes, newOriginalIds))
uniqueMissingMidNodes[newEdgeMask] = newNodeIndices
resConnectivity[missingElementIndices, missingMidNodePositions] = uniqueMissingMidNodes[uniqueMissingInverse]
return resConnectivity
def _CollectKnownMidNodes(
self,
res: Mesh,
connectivity: np.ndarray,
edgeByMidNode: Dict[int, Tuple[int, int]],
edgeByMidNodeByElementType: Dict[ED.ElementType, Dict[int, Tuple[int, int]]],
) -> Tuple[np.ndarray, np.ndarray]:
"""Collect the edges whose midside node is already known.
Parameters
----------
res : Mesh
Mesh providing the already created elements.
connectivity : np.ndarray
Connectivity array of the elements currently being processed.
edgeByMidNode : dict
Mapping from a midside node position to its edge for ``connectivity``.
edgeByMidNodeByElementType : dict
Mapping from an element type to its midside-node-to-edge mapping.
Returns
-------
tuple of np.ndarray
A pair ``(knownEdges, knownMidNodes)`` where ``knownEdges`` is an
array of sorted edges (corner node pairs) and ``knownMidNodes`` the
associated midside node indices.
"""
knownEdges = []
knownMidNodes = []
for existingElementType, existingEdgeByMidNode in edgeByMidNodeByElementType.items():
if existingElementType not in res.elements:
continue
elements = res.GetElementsOfType(existingElementType)
existingConnectivity = elements.connectivity[:elements.GetNumberOfElements(), :]
self._AppendKnownMidNodes(existingConnectivity, existingEdgeByMidNode, knownEdges, knownMidNodes)
self._AppendKnownMidNodes(connectivity, edgeByMidNode, knownEdges, knownMidNodes)
if len(knownEdges) == 0:
return np.empty((0, 2), dtype=MuscatIndex), np.empty((0,), dtype=MuscatIndex)
return np.vstack(knownEdges), np.hstack(knownMidNodes)
def _AppendKnownMidNodes(
self,
connectivity: np.ndarray,
edgeByMidNode: Dict[int, Tuple[int, int]],
knownEdges: List[np.ndarray],
knownMidNodes: List[np.ndarray],
) -> None:
"""Append the known edges and midside nodes of a connectivity array.
Parameters
----------
connectivity : np.ndarray
Connectivity array to inspect.
edgeByMidNode : dict
Mapping from a midside node position to its edge.
knownEdges : list
List of edge arrays, updated in place.
knownMidNodes : list
List of midside node arrays, updated in place.
"""
if len(connectivity) == 0:
return
for midNodePosition, edge in edgeByMidNode.items():
validMidNodes = connectivity[:, midNodePosition] >= 0
if not np.any(validMidNodes):
continue
edges = np.sort(connectivity[validMidNodes][:, edge], axis=1)
knownEdges.append(edges)
knownMidNodes.append(connectivity[validMidNodes, midNodePosition])
def _ExtractEdges(
self,
connectivity: np.ndarray,
elementIndices: np.ndarray,
midNodePositions: np.ndarray,
edgeByMidNode: Dict[int, Tuple[int, int]],
) -> np.ndarray:
"""Extract the sorted edges associated with a set of missing midside nodes.
Parameters
----------
connectivity : np.ndarray
Connectivity array of the elements.
elementIndices : np.ndarray
Indices of the elements owning the missing midside nodes.
midNodePositions : np.ndarray
Midside node positions of the missing nodes.
edgeByMidNode : dict
Mapping from a midside node position to its edge.
Returns
-------
np.ndarray
Array of shape ``(len(elementIndices), 2)`` containing the sorted
corner node pairs of the edges supporting the missing midside nodes.
"""
res = np.empty((len(elementIndices), 2), dtype=MuscatIndex)
for midNodePosition, edge in edgeByMidNode.items():
mask = midNodePositions == midNodePosition
if not np.any(mask):
continue
res[mask, :] = np.sort(connectivity[elementIndices[mask]][:, edge], axis=1)
return res
def _AsEdgeRecords(self, edges: np.ndarray) -> np.ndarray:
"""Convert an array of edges into a structured array of records.
Parameters
----------
edges : np.ndarray
Array of shape ``(n, 2)`` containing edges as corner node pairs.
Returns
-------
np.ndarray
A 1D structured array (fields ``n0`` and ``n1``) suitable for
sorting and searching operations on edges.
"""
return np.ascontiguousarray(edges).view([("n0", edges.dtype), ("n1", edges.dtype)]).ravel()
def _AddElementsToContainer(
self,
res: Mesh,
elementType: ED.ElementType,
connectivity: np.ndarray,
originalIds: np.ndarray,
oidToElementContainerAndIndex: Dict,
) -> None:
"""Add elements to the proper element container of a mesh.
Parameters
----------
res : Mesh
Mesh to which the elements are added.
elementType : ED.ElementType
Muscat element type of the elements to add.
connectivity : np.ndarray
Connectivity array of the elements.
originalIds : np.ndarray
Original Ansys identifiers of the elements.
oidToElementContainerAndIndex : dict
Mapping (updated in place) from an element original id to the tuple
(element container, local index).
"""
if len(originalIds) == 0:
return
elements = res.GetElementsOfType(elementType)
firstElementIndex = elements.GetNumberOfElements()
elements.Reserve(firstElementIndex+len(originalIds))
elements.AddNewElements(np.asarray(connectivity, dtype=MuscatIndex, order="C"), originalIds)
for localIndex, oid in enumerate(originalIds):
oidToElementContainerAndIndex[oid] = (elements, localIndex+firstElementIndex)
def _RescopeFields(self, dpf, fieldsContainer, wantedIds: np.ndarray) -> List:
"""Restrict each ElementalNodal field of a container to the wanted elements.
The rescope operates on in-memory fields and is fast, unlike passing a
mesh_scoping to the result operator, which extracts entity by entity
(~60 ms/element on some rst files).
"""
fields = []
for field in fieldsContainer:
if field.size == 0:
continue
presentIds = np.asarray(field.scoping.ids)
keepMask = np.isin(presentIds, wantedIds)
if not keepMask.any():
continue
if keepMask.all():
fields.append(field) # nothing to remove
continue
rescopeScoping = dpf.Scoping(location=dpf.locations.elemental)
rescopeScoping.ids = presentIds[keepMask] # array numpy, pas .tolist()
rescope = dpf.operators.scoping.rescope()
rescope.inputs.fields.connect(field)
rescope.inputs.mesh_scoping.connect(rescopeScoping)
fields.append(rescope.outputs.fields_as_field())
return fields
def _ReadElementalNodalDpfField(self, dpf, rst, fieldname: str, timeIndex: int, idsByDim: Dict[int, np.ndarray]):
"""Read an elemental-nodal DPF field, trying each dimensionality in turn.
Scoped reads are attempted per dimensionality, highest first, stopping at
the first non-empty result: solid-carried fields (stress, strain...) are
read on the fast solid-only scoping, while shell-only results are
recovered on lower dimensionalities. If a scoped read raises (e.g.
element_nodal_forces with degenerated contact elements), a single
unscoped read (which streams the rst) followed by an in-memory rescope
is used instead.
Returns
-------
Field or None
The DPF field, or None if the result could not be read or carries
no data on the requested elements.
"""
fields = []
for dim, ids in idsByDim.items():
try:
scoping = dpf.Scoping(location=dpf.locations.elemental)
scoping.ids = ids
fieldsContainer = getattr(rst.results, fieldname)(
mesh_scoping=scoping, time_scoping=int(timeIndex)
).outputs.fields_container()
fields = [f for f in fieldsContainer if f.size]
except Exception as e:
Warning(f"{fieldname}: scoped read failed on dimensionality {dim} ({e}); falling back to unscoped read")
try:
fieldsContainer = getattr(rst.results, fieldname)(
time_scoping=int(timeIndex)
).outputs.fields_container()
allIds = np.concatenate(list(idsByDim.values()))
fields = self._RescopeFields(dpf, fieldsContainer, allIds)
except Exception as e2:
Warning(f"{fieldname} could not be read: {e2}")
return None
break # the unscoped read already covers every dimensionality
if fields:
break
if not fields:
Warning(f"{fieldname}: no data on the requested elements")
return None
if len(fields) > 1:
Warning(f"{fieldname}: {len(fields)} fields returned, only the first one is used")
return fields[0]
[docs]
def Read(self, time: Optional[MuscatFloat] = None, timeIndex: Optional[MuscatIndex] = None, fromElementalNodalToNodal: bool = True, dimensionality: Optional[Union[int, List[int]]] = None) -> Mesh:
"""Read the mesh and the fields stored in the rst file.
Parameters
----------
time : MuscatFloat, optional
Time at which the fields are read, by default None.
timeIndex : MuscatIndex, optional
Time index at which the fields are read, by default None.
fromElementalNodalToNodal : bool, optional
If True, elemental-nodal fields are averaged at the nodes; otherwise
they are stored as discontinuous finite element fields, by default True.
dimensionality : int or list of int, optional
Dimensionality (or list of dimensionalities) of the elements on which
the elemental-nodal fields (e.g. stress, strain) are read. When
``None`` (default), the fields are read on every geometry present in
the mesh, whatever its dimensionality. This is important when several
geometries of different dimensions (e.g. a solid and a shell) coexist
in the same rst file and both carry results.
Returns
-------
Mesh
Output Mesh object containing the results.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-core`` and ``ansys-dpf-post`` are not installed.
"""
post = self._ImportPost()
# If the simulation is modal then it sets self.modal to True because it needs
# to accomodate to the 1 index of modal simulations
try:
self.modal = int(isinstance(post.load_simulation(self.fileName), post.ModalMechanicalSimulation))
except NotImplementedError:
self.modal = 0
mesh = self.MeshRead()
for nodefield, operators in self.nodal.items():
field = self.ReadField(fieldname=nodefield, time=time, timeIndex=timeIndex)
if field is not None:
self._StoreArrayFieldComponents(mesh.nodeFields, operators, field)
for elemField, operators in self.elementary.items():
field = self.ReadField(fieldname=elemField, time=time, timeIndex=timeIndex)
if field is not None:
self._StoreArrayFieldComponents(mesh.elemFields, operators, field)
if fromElementalNodalToNodal:
for elemField, operators in self.elemental_nodal.items():
field = self.ReadField(fieldname=elemField, time=time, timeIndex=timeIndex, dimensionality=dimensionality)
if field is not None:
self._StoreArrayFieldComponents(mesh.nodeFields, operators, field)
else:
for elemField, operators in self.elemental_nodal.items():
FEfield = self.ReadField(fieldname=elemField, time=time, timeIndex=timeIndex, fromElementalNodalToNodal=False, dimensionality=dimensionality)
if FEfield is not None:
self._StoreFEFieldComponents(mesh.elemFields, operators, FEfield)
return mesh
#: Component suffixes used when an array field exposes more components than
#: operator names. Indexed by the number of components.
_COMPONENT_SUFFIXES = {
1: ("",),
3: ("X", "Y", "Z"),
6: ("XX", "YY", "ZZ", "XY", "YZ", "XZ"),
}
def _StoreArrayFieldComponents(self, storage: Dict, operators: List[str], field: np.ndarray) -> None:
"""Store each component of an array field under its operator name.
When the number of operator names matches the number of components, each
component is stored under its operator name. When a single operator name
is provided for a multi-component field (e.g. a stress or strain tensor),
every component is stored under ``"{operatorName}{suffix}"`` where the
suffix follows the Ansys convention (``X``/``Y``/``Z`` for vectors,
``XX``/``YY``/``ZZ``/``XY``/``YZ``/``XZ`` for symmetric tensors, and the
component index otherwise).
Parameters
----------
storage : dict
Field container (e.g. ``mesh.nodeFields``) updated in place.
operators : list of str
Operator names used as keys. Either one per field component or a
single base name for the whole field.
field : np.ndarray
Field array whose columns are the components to store.
"""
numberOfComponents = field.shape[1]
if len(operators) >= numberOfComponents:
for componentIndex, op in enumerate(operators[:numberOfComponents]):
storage[op] = field[:, componentIndex]
return
baseName = operators[0]
suffixes = self._COMPONENT_SUFFIXES.get(
numberOfComponents, tuple(str(componentIndex) for componentIndex in range(numberOfComponents))
)
for componentIndex, suffix in enumerate(suffixes):
storage[f"{baseName}{suffix}"] = field[:, componentIndex]
def _StoreFEFieldComponents(self, storage: Dict, operators: List[str], fields: List[FEField]) -> None:
"""Store the cell representation of each finite element field component.
Parameters
----------
storage : dict
Field container (e.g. ``mesh.elemFields``) updated in place.
operators : list of str
Operator names used as keys, one per field component.
fields : list of FEField
Finite element fields, one per component.
"""
for op, field in zip(operators, fields):
storage[op] = field.GetCellRepresentation()
[docs]
def ReadField(self, fieldname: str, time: Optional[MuscatFloat] = None, timeIndex: Optional[MuscatIndex] = None, fromElementalNodalToNodal: bool = True, dimensionality: Optional[Union[int, List[int]]] = None) -> Union[np.ndarray, List[FEField], None]:
"""Read a single field from the rst file.
Parameters
----------
fieldname : str
Name of the field to be read.
time : MuscatFloat, optional
Time at which the field is read, by default None.
timeIndex : MuscatIndex, optional
Time index at which the field is read, by default None.
fromElementalNodalToNodal : bool, optional
For elemental-nodal fields, whether to average the field at the
nodes, by default True.
dimensionality : int or list of int, optional
For elemental-nodal fields, the dimensionality (or list of
dimensionalities) of the elements on which the field is read. When
``None`` (default), the field is read on every geometry present in
the mesh, whatever its dimensionality. This is required when several
geometries of different dimensions coexist in the same rst file and
more than one carries elemental-nodal results.
Returns
-------
np.ndarray or list of FEField
For a nodal or elemental field, a numpy array of the field values.
For an elemental-nodal field:
- if ``fromElementalNodalToNodal`` is True (default), a numpy array
of the averaged value at the nodes;
- if ``fromElementalNodalToNodal`` is False, a list of discontinuous
finite element fields (one for each component of the field).
``None`` is returned if the field could not be broadcast to the mesh.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-core`` and ``ansys-dpf-post`` are not installed.
KeyError
If the requested field cannot be found.
"""
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()
dpf = self._ImportDpf()
if fieldname in self.nodal:
scoping = rst.metadata.meshed_region.nodes.scoping
nbentities = self.output.GetNumberOfNodes()
mapping = self.IdsToNodes
elif fieldname in self.elementary:
scoping = rst.metadata.meshed_region.elements.scoping
nbentities = self.output.GetNumberOfElements()
mapping = self.IdsToElements
elif fieldname in self.elemental_nodal:
idsByDim = self._GetElementalNodalScopingIdsByDimensionality(fieldname, dimensionality)
if not idsByDim:
Warning(f"{fieldname}: empty element selection (dimensionality={dimensionality})")
return None
nbentities = self.output.GetNumberOfNodes()
mapping = self.IdsToNodes
else:
raise KeyError("unable to find field " + str(fieldname))
if fieldname in self.elemental_nodal:
field = self._ReadElementalNodalDpfField(dpf, rst, fieldname, timeIndex, idsByDim)
if field is None:
return None
else:
try:
field = getattr(rst.results, fieldname)(
mesh_scoping=scoping, time_scoping=int(timeIndex)
).outputs.fields_container()[0]
except Exception as e:
Warning(f"{fieldname} could not be read: {e}")
return None
if fieldname in self.elemental_nodal:
if not fromElementalNodalToNodal:
############
# generate a discontinuous space to hold the non-averaged elemental-nodal field
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, LagrangeSpaceGeo
if dimensionality is None:
dimensionality = [1, 2, 3]
elementFilterSolid = ElementFilter(dimensionality=dimensionality)
support_ids = np.asarray(field.scoping.ids)
data = field.data
if data.ndim == 1:
data = data[:, None]
muscat_offset = self.output.ComputeGlobalOffset()
# warning : depending of the nature of the field (stress, temperature...), the number of entries in field.data can correspond to the order of the mesh or to a linear version
count = 0
for selection in elementFilterSolid.IterOnTypesOnly(self.output):
count += ED.numberOfNodes[selection.elementType]*selection.elements.GetNumberOfElements()
if count == data.shape[0]:
space = LagrangeSpaceGeo
else:
space = LagrangeSpaceP1
numbering = ComputeDofNumbering(self.output, space=space, elementFilter=elementFilterSolid, discontinuous=True)
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
nodesPerElements_extraction = nodesPerElements[support_ids]
elementsTypes_extraction = elementsTypes[support_ids]
ansys_offset = np.empty(len(support_ids)+1, dtype=MuscatIndex)
ansys_offset[0] = 0
ansys_offset[1:] = np.add.accumulate(nodesPerElements_extraction)
res = []
for i in range(data.shape[1]):
res_i = FEField(name=fieldname, mesh=self.output, space=space, numbering=numbering)
res_i.Allocate()
cpt = 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)
ansys_mask = fillingMask + ansys_offset[index, None]
newData = data[ansys_mask, i]
targetIndices = (self.IdsToElements[support_ids[index]] - muscat_offset[selection.elementType])*numberOfNodesPerElement + cpt
muscat_mask = fillingMask + targetIndices[:, None]
res_i.data[muscat_mask.ravel()] = newData.ravel()
cpt += selection.elements.GetNumberOfElements()*numberOfNodesPerElement
res.append(res_i)
return res
#############################
else:
#############################
# construct the nodal averaged field
if field.size == 0:
Warning(f"{fieldname}: the DPF scoping selected no element (dimensionality={dimensionality}). "
"The field is probably carried by a geometry of another dimensionality; "
"try passing the appropriate 'dimensionality' argument.")
return None
to_nodal = dpf.operators.averaging.elemental_nodal_to_nodal()
to_nodal.inputs.field.connect(field)
field_nodal = to_nodal.outputs.field()
support_ids = np.asarray(field_nodal.scoping.ids)
data = field_nodal.data
#############################
else:
support_ids = np.asarray(field.scoping.ids)
data = field.data
if data.ndim == 1:
data = data[:, None]
res = np.zeros((nbentities, data.shape[1]))
try:
res[mapping[support_ids]] = data
except ValueError:
Warning(f"{fieldname} of shape {data.shape} could not be broadcast to indexing result of shape {support_ids.shape}")
return None
return res
[docs]
def SetTimeToRead(self, time: Optional[MuscatFloat] = None, timeIndex: Optional[MuscatIndex] = None) -> MuscatIndex:
"""Set the time at which the data is read.
Parameters
----------
time : MuscatFloat, optional
Time at which the data is read, by default None.
timeIndex : MuscatIndex, optional
Time index at which the data is read, by default None.
Returns
-------
MuscatIndex
Time index at which the data is read.
Raises
------
ValueError
If both ``time`` and ``timeIndex`` are specified.
"""
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][0]
if (time is not None) and (timeIndex is not None):
raise ValueError("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][0]
[docs]
def GetAvailableTimes(self) -> np.ndarray:
"""Return 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[Results]:
"""Return the available results together with their operator name.
Returns
-------
list of Results
List of :class:`Results` objects exposing the ``name``,
``operator_name``, ``native_location`` and ``n_components`` attributes.
"""
rst = self.GetModel()
return [Results(r.name, r.operator_name, r.native_location, r.n_components) for r in rst.metadata.result_info.available_results]
[docs]
def GetAvailablePostResults(self) -> List[str]:
"""Return the available post-processing results of the simulation.
Returns
-------
list
Names of the callable post-processing results exposed by the
``ansys-dpf-post`` simulation object.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-post`` is not installed.
"""
post = self._ImportPost()
simulation = post.load_simulation(self.fileName)
return [attr for attr in dir(simulation) if callable(getattr(simulation, attr)) and not attr.startswith("__")]
[docs]
def GetPostResults(self, post_names: Optional[List[str]] = None, locations: Optional[List[str]] = None) -> None:
"""Read post-processing results and store them in the output mesh.
Parameters
----------
post_names : list of str, optional
Names of the post-processing results to read, by default None.
locations : list of str, optional
Location (``"nodal"`` or ``"elemental"``) associated with each
requested result, by default None.
Raises
------
ModuleNotFoundError
If ``ansys-dpf-post`` is not installed.
"""
post_names = post_names if post_names is not None else []
locations = locations if locations is not None else []
post = self._ImportPost()
simulation = post.load_simulation(self.fileName)
for attribute, location in zip(post_names, locations):
result = getattr(simulation, attribute)(location=getattr(post.locations, location))
if location == "nodal":
storage = self.output.nodeFields
elif location == "elemental":
storage = self.output.elemFields
else:
continue
data = result.array
reordered = np.empty_like(data)
reordered[result.mesh_index.values-1] = data
self._StorePostResultComponents(storage, attribute, result, reordered)
def _StorePostResultComponents(self, storage: Dict, attribute: str, result, data: np.ndarray) -> None:
"""Store a post-processing result, splitting multi-component data per component.
A scalar result is stored under ``attribute``. A multi-component result
(e.g. a stress tensor) is split into one field per component, named
``"{attribute}_{componentLabel}"`` (for instance ``"stress_XX"``). The
component labels are read from the result when available and fall back to
the component index otherwise.
Parameters
----------
storage : dict
Field container (e.g. ``mesh.nodeFields``) updated in place.
attribute : str
Base name of the result.
result : DataFrame
The ``ansys-dpf-post`` result object, used to recover the component
labels.
data : np.ndarray
Result values, already reordered to match the mesh entities.
"""
if data.ndim == 1 or data.shape[1] == 1:
storage[attribute] = data.ravel()
return
componentLabels = self._GetComponentLabels(result, data.shape[1])
for componentIndex, label in enumerate(componentLabels):
storage[f"{attribute}_{label}"] = data[:, componentIndex]
@staticmethod
def _GetComponentLabels(result, numberOfComponents: int) -> List[str]:
"""Return the component labels of a post-processing result.
Parameters
----------
result : DataFrame
The ``ansys-dpf-post`` result object.
numberOfComponents : int
Number of components of the result.
Returns
-------
list of str
One label per component. When the labels cannot be recovered from
the result, the component indices (as strings) are returned.
"""
try:
componentIndex = result.columns.find_label("comp", "components")
labels = [str(label) for label in result.columns[componentIndex].values]
if len(labels) == numberOfComponents:
return labels
except Exception:
pass
return [str(componentIndex) for componentIndex in range(numberOfComponents)]
RegisterReaderClass(".rst", RstReader)
[docs]
def CheckIntegrity(GUI: bool = False) -> str:
"""Check the integrity of the RstReader module.
Parameters
----------
GUI : bool, optional
Whether to run the check in GUI mode, by default False.
Returns
-------
str
``"ok"`` if the check succeeds, or a ``"skip: ..."`` message when the
Ansys dependencies are not available.
"""
from Muscat.Helpers.CheckTools import MustFailFunctionWith
reader = RstReader()
storage = {}
reader._StoreArrayFieldComponents(storage, ["A_X", "A_Y", "A_Z"], np.ones((4, 3)))
assert list(storage.keys()) == ["A_X", "A_Y", "A_Z"]
storage = {}
reader._StoreArrayFieldComponents(storage, ["S"], np.arange(24, dtype=MuscatFloat).reshape(4, 6))
assert list(storage.keys()) == ["SXX", "SYY", "SZZ", "SXY", "SYZ", "SXZ"]
np.testing.assert_allclose(storage["SYY"], np.arange(24).reshape(4, 6)[:, 1])
mesh = Mesh()
mesh.nodes = np.array([
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[0.5, 0.0, 0.0],
[0.5, 0.5, 0.0],
[0.0, 0.5, 0.0],
[0.0, 0.0, 0.5],
[0.5, 0.0, 0.5],
[0.0, 0.5, 0.5],
], dtype=MuscatFloat)
mesh.originalIDNodes = np.arange(10, dtype=MuscatIndex)
oidToElementContainerAndIndex = {}
connectivity = np.array([
[0, 1, 2, 3, 4, 5, 6, 7, 8, -1],
], dtype=MuscatIndex)
originalIds = np.array([11], dtype=MuscatIndex)
reader._AddElementsFromAnsysConnectivity(mesh, ED.Tetrahedron_10, connectivity, originalIds, oidToElementContainerAndIndex)
mesh.PrepareForOutput()
assert mesh.GetElementsOfType(ED.Tetrahedron_10).GetNumberOfElements() == 1
assert mesh.nodes.shape[0] == 11
np.testing.assert_allclose(mesh.nodes[10], [0.0, 0.5, 0.5])
mesh = Mesh()
mesh.nodes = np.zeros((10, 3), dtype=MuscatFloat)
mesh.originalIDNodes = np.arange(10, dtype=MuscatIndex)
oidToElementContainerAndIndex = {}
connectivity = np.array([
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
[0, 1, 2, 3, -1, -1, -1, -1, -1, -1],
], dtype=MuscatIndex)
originalIds = np.array([12, 13], dtype=MuscatIndex)
reader._AddElementsFromAnsysConnectivity(mesh, ED.Tetrahedron_10, connectivity, originalIds, oidToElementContainerAndIndex)
mesh.PrepareForOutput()
assert mesh.GetElementsOfType(ED.Tetrahedron_10).GetNumberOfElements() == 2
assert mesh.nodes.shape[0] == 10
mesh = Mesh()
mesh.nodes = np.array([
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.5, 0.0, 0.0],
[0.5, 0.5, 0.0],
[0.0, 0.5, 0.0],
], dtype=MuscatFloat)
mesh.originalIDNodes = np.arange(6, dtype=MuscatIndex)
oidToElementContainerAndIndex = {}
connectivity = np.array([
[0, 1, 2, -1, 4, 5],
], dtype=MuscatIndex)
originalIds = np.array([15], dtype=MuscatIndex)
reader._AddElementsFromAnsysConnectivity(mesh, ED.Triangle_6, connectivity, originalIds, oidToElementContainerAndIndex)
mesh.PrepareForOutput()
assert mesh.GetElementsOfType(ED.Triangle_6).GetNumberOfElements() == 1
assert mesh.nodes.shape[0] == 7
np.testing.assert_allclose(mesh.nodes[6], [0.5, 0.0, 0.0])
invalidConnectivity = np.array([[-1, 1, 2, 3, 4, 5, 6, 7, 8, 9]], dtype=MuscatIndex)
MustFailFunctionWith(RuntimeError, reader._AddElementsFromAnsysConnectivity, Mesh(), ED.Tetrahedron_10, invalidConnectivity, np.array([14], dtype=MuscatIndex), {})
# ------------------------------------------------------------------
# Results dataclass-like container
# ------------------------------------------------------------------
result = Results("stress", "S", "ElementalNodal", 6)
assert result.name == "stress"
assert result.operator_name == "S"
assert result.native_location == "ElementalNodal"
assert result.n_components == 6
assert repr(result) == "stress(S, ElementalNodal, 6)"
# ------------------------------------------------------------------
# SetFileName and Reset
# ------------------------------------------------------------------
reader = RstReader()
reader.SetFileName("first.rst")
assert reader.fileName == "first.rst"
reader.model = "dummy" # simulate a loaded model
reader.SetFileName("first.rst") # same name: Reset must not be triggered
assert reader.model == "dummy"
reader.SetFileName("second.rst") # different name: Reset is triggered
assert reader.fileName == "second.rst"
assert reader.model is None
# GetModel must fail when the fileName has not been set.
emptyReader = RstReader()
MustFailFunctionWith(RuntimeError, emptyReader.GetModel)
# ------------------------------------------------------------------
# _GetMidNodeEdgesByElementType
# ------------------------------------------------------------------
edgesByType = reader._GetMidNodeEdgesByElementType()
assert edgesByType[ED.Triangle_6][3] == (0, 1)
assert edgesByType[ED.Quadrangle_8][7] == (3, 0)
assert edgesByType[ED.Tetrahedron_10][9] == (2, 3)
# ------------------------------------------------------------------
# _AsEdgeRecords
# ------------------------------------------------------------------
edges = np.array([[0, 1], [2, 3]], dtype=MuscatIndex)
records = reader._AsEdgeRecords(edges)
assert records.shape == (2,)
assert records[0]["n0"] == 0 and records[0]["n1"] == 1
assert records[1]["n0"] == 2 and records[1]["n1"] == 3
# ------------------------------------------------------------------
# _AddElementsToContainer (empty branch returns immediately)
# ------------------------------------------------------------------
emptyMesh = Mesh()
reader._AddElementsToContainer(
emptyMesh,
ED.Tetrahedron_4,
np.empty((0, 4), dtype=MuscatIndex),
np.empty(0, dtype=MuscatIndex),
{},
)
assert ED.Tetrahedron_4 not in emptyMesh.elements
# ------------------------------------------------------------------
# SetTimeToRead (every branch)
# ------------------------------------------------------------------
timeReader = RstReader()
timeReader.time = np.array([0.0, 1.0, 2.0], dtype=MuscatFloat)
# Both None with timeToRead == -1: the last time index is returned.
timeReader.timeToRead = -1
assert timeReader.SetTimeToRead() == 2
assert timeReader.timeToRead == 2.0
# Both None with a previously stored timeToRead.
timeReader.timeToRead = 1.0
assert timeReader.SetTimeToRead() == 1
# A time index is given directly.
assert timeReader.SetTimeToRead(timeIndex=0) == 0
assert timeReader.timeToRead == 0.0
# An explicit time is matched against the available times.
assert timeReader.SetTimeToRead(time=1.0) == 1
assert timeReader.timeToRead == 1.0
# Specifying both time and timeIndex is an error.
MustFailFunctionWith(ValueError, timeReader.SetTimeToRead, 1.0, 1)
# ------------------------------------------------------------------
# _StoreArrayFieldComponents (generic suffix branch)
# ------------------------------------------------------------------
storage = {}
reader._StoreArrayFieldComponents(storage, ["G"], np.arange(8, dtype=MuscatFloat).reshape(2, 4))
assert list(storage.keys()) == ["G0", "G1", "G2", "G3"]
# ------------------------------------------------------------------
# _GetElementalNodalScopingIdsByDimensionality (contact branch)
# ------------------------------------------------------------------
contactReader = RstReader()
contactReader.contactElementIds = np.empty(0, dtype=MuscatIndex)
assert contactReader._GetElementalNodalScopingIdsByDimensionality("contact_status") == {}
contactReader.contactElementIds = np.array([7, 8, 9], dtype=MuscatIndex)
contactIds = contactReader._GetElementalNodalScopingIdsByDimensionality("contact_status")
assert set(contactIds.keys()) == {0}
np.testing.assert_array_equal(contactIds[0], [7, 8, 9])
# ------------------------------------------------------------------
# _GetComponentLabels (fallback to indices when labels are unavailable)
# ------------------------------------------------------------------
labels = RstReader._GetComponentLabels(object(), 3)
assert labels == ["0", "1", "2"]
# ------------------------------------------------------------------
# _StorePostResultComponents (scalar and multi-component branches)
# ------------------------------------------------------------------
scalarStorage = {}
reader._StorePostResultComponents(scalarStorage, "temperature", object(), np.arange(3, dtype=MuscatFloat))
np.testing.assert_array_equal(scalarStorage["temperature"], [0.0, 1.0, 2.0])
tensorStorage = {}
reader._StorePostResultComponents(tensorStorage, "stress", object(), np.arange(6, dtype=MuscatFloat).reshape(2, 3))
assert list(tensorStorage.keys()) == ["stress_0", "stress_1", "stress_2"]
# ------------------------------------------------------------------
# _StoreFEFieldComponents (delegates to GetCellRepresentation)
# ------------------------------------------------------------------
class _FakeFEField:
def __init__(self, values):
self._values = values
def GetCellRepresentation(self):
return self._values
feStorage = {}
reader._StoreFEFieldComponents(feStorage, ["a", "b"], [_FakeFEField([1.0]), _FakeFEField([2.0])])
assert feStorage["a"] == [1.0]
assert feStorage["b"] == [2.0]
# ------------------------------------------------------------------
# Import helpers must fail with an explicit message when DPF is absent
# ------------------------------------------------------------------
import Muscat.TestData as MuscatTestData
from Muscat.Helpers.IO.Which import Which
try:
import ansys.dpf.core as dpf
except ImportError:
MustFailFunctionWith(ModuleNotFoundError, RstReader._ImportDpf)
MustFailFunctionWith(ModuleNotFoundError, RstReader._ImportPost)
return 'skip: ansys.dpf.core not available '
ansys_exec = "runwb2"
ansysExec = Which(ansys_exec)
if ansysExec == '' or ansysExec == None:
return "skip: ansys not available"
file = MuscatTestData.GetTestDataPath()+"file.rst"
ReadRst(file)
reader = RstReader()
reader.SetFileName(file)
times = reader.GetAvailableTimes()
reader.Read(time=times[-1], fromElementalNodalToNodal=False)
reader.GetAvailableResultsAndOperators()
reader.GetAvailablePostResults()
reader.GetPostResults(['stress'], ["nodal"])
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(True))# pragma: no cover