diff --git a/doc/source/conf.py b/doc/source/conf.py index 6ba5d1f60..312555382 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -3,7 +3,7 @@ from datetime import datetime import os -from ansys_sphinx_theme import ansys_favicon, get_version_match, pyansys_logo_black +from ansys_sphinx_theme import ansys_favicon, get_version_match import numpy as np import pyvista from pyvista.plotting.utilities.sphinx_gallery import DynamicScraper @@ -28,7 +28,6 @@ release = version = __version__ # Select desired logo, theme, and declare the html title -html_logo = pyansys_logo_black html_favicon = ansys_favicon html_theme = "ansys_sphinx_theme" html_short_title = html_title = "PyDPF Composites" @@ -41,6 +40,7 @@ cname = os.environ.get("DOCUMENTATION_CNAME", "composites.dpf.docs.pyansys.com") html_theme_options = { + "logo": "pyansys", "github_url": "https://github.com/ansys/pydpf-composites", "show_prev_next": False, "show_breadcrumbs": True, diff --git a/examples/013_thermal_example.py b/examples/013_thermal_example.py new file mode 100644 index 000000000..53b98b683 --- /dev/null +++ b/examples/013_thermal_example.py @@ -0,0 +1,128 @@ +# Copyright (C) 2023 - 2024 ANSYS, Inc. and/or its affiliates. +# SPDX-License-Identifier: MIT +# +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +""" +.. _thermal_example: + +Thermal analysis +---------------- + +PyDPF Composites can also be used to post-process thermal analyses. +In this case, the simulation is a two-step analysis where the results of +a thermal analysis are an input of the structural analysis. So, the RST +contains temperature and structural results. +The example mimics a PCB which was modeled with Ansys Composites PrePost (ACP). +where the solid model feature of ACP is used to generate the volume mesh. + +In detail, the example shows how to extract the temperatures for a specific ply, +and a specific material. + +.. note:: + + When using a Workbench project, + use the :func:`.get_composite_files_from_workbench_result_folder` + method to obtain the input files. + +""" + +# %% +# Set up analysis +# ~~~~~~~~~~~~~~~ +# Setting up the analysis consists of loading the required modules, connecting to the +# DPF server, and retrieving the example files. +# +import ansys.dpf.core as dpf +import numpy as np + +from ansys.dpf.composites.composite_model import CompositeModel +from ansys.dpf.composites.constants import TEMPERATURE_COMPONENT +from ansys.dpf.composites.example_helper import get_continuous_fiber_example_files +from ansys.dpf.composites.layup_info import get_all_analysis_ply_names +from ansys.dpf.composites.ply_wise_data import SpotReductionStrategy, get_ply_wise_data +from ansys.dpf.composites.select_indices import get_selected_indices_by_dpf_material_ids +from ansys.dpf.composites.server_helpers import connect_to_or_start_server + +server = connect_to_or_start_server() +composite_files = get_continuous_fiber_example_files(server, "thermal_solid") + +# %% +# Initialize the model +# ~~~~~~~~~~~~~~~~~~~~ +# The composite model is initialized with the composite files and the server. +# It provides access to the mesh, results, lay-up and materials +composite_model = CompositeModel(composite_files, server) + +# %% +# Get Results - Temperatures +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# The temperatures are stored under structural_temperature +temp_op = composite_model.core_model.results.structural_temperature() +temperatures_fc = temp_op.outputs.fields_container() + +# %% +# Ply-wise results +# ~~~~~~~~~~~~~~~~ +# Ply-wise results can be easily extracted using the function +# :func:`.get_ply_wise_data` and by passing the ply name. + +all_ply_names = get_all_analysis_ply_names(composite_model.get_mesh()) +print(all_ply_names) + +nodal_values = get_ply_wise_data( + field=temperatures_fc, + ply_name="P1L1__ModelingPly.8", + mesh=composite_model.get_mesh(), + component=TEMPERATURE_COMPONENT, + spot_reduction_strategy=SpotReductionStrategy.MAX, + requested_location=dpf.locations.nodal, +) + +composite_model.get_mesh().plot(nodal_values) + +# %% +# Material-wise results +# ~~~~~~~~~~~~~~~~~~~~~ +# It is also possible to filter the results by material. +# In this example the element-wise maximum temperature +# is extracted for the material `Honeycomb Aluminum Alloy`. +print(composite_model.material_names) +material_id = composite_model.material_names["Honeycomb Aluminum Alloy"] + +# get the last result field +temperatures_field = temperatures_fc[-1] + +material_result_field = dpf.field.Field(location=dpf.locations.elemental, nature=dpf.natures.scalar) +# performance optimization: use a local field instead of a field which is pushed to the server +with material_result_field.as_local_field() as local_result_field: + element_ids = temperatures_field.scoping.ids + + for element_id in element_ids: + element_info = composite_model.get_element_info(element_id) + assert element_info is not None + if material_id in element_info.dpf_material_ids: + temp_data = temperatures_field.get_entity_data_by_id(element_id) + selected_indices = get_selected_indices_by_dpf_material_ids(element_info, [material_id]) + + value = np.max(temp_data[selected_indices]) + local_result_field.append([value], element_id) + +composite_model.get_mesh().plot(material_result_field) diff --git a/pyproject.toml b/pyproject.toml index 6ca55c195..1231a4489 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ build-backend = "poetry.core.masonry.api" # Check https://python-poetry.org/docs/pyproject/ for all available sections name = "ansys-dpf-composites" # Switch to released version of dpf core releasing pydpf-composites! -version = "0.6.0" +version = "0.6.1" description = "Post-processing of composite structures based on Ansys DPF" license = "MIT" authors = ["ANSYS, Inc. "] diff --git a/src/ansys/dpf/composites/_indexer.py b/src/ansys/dpf/composites/_indexer.py index 193f9e111..6cedb4bbb 100644 --- a/src/ansys/dpf/composites/_indexer.py +++ b/src/ansys/dpf/composites/_indexer.py @@ -22,7 +22,7 @@ """Indexer helper classes.""" from dataclasses import dataclass -from typing import Optional, Protocol, cast +from typing import Optional, Protocol, Union, cast from ansys.dpf.core import Field, PropertyField, Scoping import numpy as np @@ -52,18 +52,56 @@ def setup_index_by_id(scoping: Scoping) -> IndexToId: return IndexToId(mapping=indices, max_id=len(indices) - 1) -class PropertyFieldIndexerSingleValue(Protocol): +class PropertyFieldIndexerProtocol(Protocol): """Protocol for single value property field indexer.""" def by_id(self, entity_id: int) -> Optional[np.int64]: - """Get index by id.""" + """ + Get index by id. + + Note: An exception is thrown if the entry has multiple values. + """ + + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.int64]]: + """Get indices by id.""" + + +def _has_data_pointer(field: Union[PropertyField, Field]) -> bool: + if ( + field._data_pointer is not None # pylint: disable=protected-access + and field._data_pointer.any() # pylint: disable=protected-access + ): + return True + return False + + +def get_property_field_indexer( + field: PropertyField, no_bounds_check: bool +) -> PropertyFieldIndexerProtocol: + """Get indexer for a property field. + + Parameters + ---------- + field: property field + no_bounds_check: whether to get the indexer w/o bounds check. More performant but less safe. + """ + if no_bounds_check: + if _has_data_pointer(field): + return PropertyFieldIndexerWithDataPointerNoBoundsCheck(field) + return PropertyFieldIndexerNoDataPointerNoBoundsCheck(field) + if _has_data_pointer(field): + return PropertyFieldIndexerWithDataPointer(field) + return PropertyFieldIndexerNoDataPointer(field) -class PropertyFieldIndexerArrayValue(Protocol): - """Protocol for array valued property field indexer.""" +class FieldIndexexProtocol(Protocol): + """Protocol for single value field indexer.""" - def by_id(self, entity_id: int) -> Optional[NDArray[np.int64]]: - """Get index by id.""" + def by_id(self, entity_id: int) -> Optional[np.double]: + """Get value by id.""" + + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.double]]: + """Get values by id.""" # General comment for all Indexer: @@ -81,10 +119,15 @@ class PropertyFieldIndexerNoDataPointer: def __init__(self, field: PropertyField): """Create indexer and get data.""" - index_by_id = setup_index_by_id(field.scoping) - self._indices = index_by_id.mapping - self._max_id = index_by_id.max_id - self._data: NDArray[np.int64] = np.array(field.data, dtype=np.int64) + if field.scoping.size > 0: + index_by_id = setup_index_by_id(field.scoping) + self._indices = index_by_id.mapping + self._max_id = index_by_id.max_id + self._data: NDArray[np.int64] = np.array(field.data, dtype=np.int64) + else: + self._indices = np.array([], dtype=np.int64) + self._data = np.array([], dtype=np.int64) + self._max_id = 0 def by_id(self, entity_id: int) -> Optional[np.int64]: """Get index by id. @@ -101,35 +144,17 @@ def by_id(self, entity_id: int) -> Optional[np.int64]: return None return cast(np.int64, self._data[idx]) - -class FieldIndexerNoDataPointer: - """Indexer for a dpf field with no data pointer.""" - - def __init__(self, field: Field): - """Create indexer and get data.""" - if field.scoping.size > 0: - index_by_id = setup_index_by_id(field.scoping) - self._indices = index_by_id.mapping - self._max_id = index_by_id.max_id - self._data: NDArray[np.double] = np.array(field.data, dtype=np.double) - else: - self._indices = np.array([], dtype=np.int64) - self._max_id = 0 - self._data = np.array([], dtype=np.double) - - def by_id(self, entity_id: int) -> Optional[np.double]: - """Get index by ID. + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.int64]]: + """Get indices by id. Parameters ---------- entity_id """ - if entity_id > self._max_id: - return None - idx = self._indices[entity_id] - if idx < 0: + value = self.by_id(entity_id) + if value is None: return None - return cast(np.double, self._data[idx]) + return np.array([value], dtype=np.int64) class PropertyFieldIndexerNoDataPointerNoBoundsCheck: @@ -159,6 +184,18 @@ def by_id(self, entity_id: int) -> Optional[np.int64]: else: return cast(np.int64, self._data[self._indices[entity_id]]) + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.int64]]: + """Get indices by id. + + Parameters + ---------- + entity_id + """ + value = self.by_id(entity_id) + if value is None: + return None + return np.array([value], dtype=np.int64) + class PropertyFieldIndexerWithDataPointer: """Indexer for a property field with data pointer.""" @@ -183,9 +220,18 @@ def __init__(self, field: PropertyField): self._n_components = 0 self._data_pointer = np.array([], dtype=np.int64) - def by_id(self, entity_id: int) -> Optional[NDArray[np.int64]]: + def by_id(self, entity_id: int) -> Optional[np.int64]: """Get index by ID. + Parameters + ---------- + entity_id + """ + raise NotImplementedError("PropertyFieldIndexerWithDataPointer does not support by_id.") + + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.int64]]: + """Get indices by ID. + Parameters ---------- entity_id @@ -206,44 +252,50 @@ def by_id(self, entity_id: int) -> Optional[NDArray[np.int64]]: ) -class FieldIndexerWithDataPointer: - """Indexer for a dpf field with data pointer.""" +class PropertyFieldIndexerWithDataPointerNoBoundsCheck: + """Indexer for a property field with data pointer and no bounds checks.""" - def __init__(self, field: Field): + def __init__(self, field: PropertyField): """Create indexer and get data.""" if field.scoping.size > 0: index_by_id = setup_index_by_id(field.scoping) self._indices = index_by_id.mapping - self._max_id = index_by_id.max_id - self._data: NDArray[np.double] = np.array(field.data, dtype=np.double) + self._data: NDArray[np.int64] = np.array(field.data, dtype=np.int64) self._n_components = field.component_count self._data_pointer: NDArray[np.int64] = np.append( field._data_pointer, len(self._data) * self._n_components ) else: - self._max_id = 0 self._indices = np.array([], dtype=np.int64) - self._data = np.array([], dtype=np.double) + self._data = np.array([], dtype=np.int64) self._n_components = 0 self._data_pointer = np.array([], dtype=np.int64) - def by_id(self, entity_id: int) -> Optional[NDArray[np.double]]: + def by_id(self, entity_id: int) -> Optional[np.int64]: """Get index by ID. Parameters ---------- entity_id """ - if entity_id > self._max_id: - return None + raise NotImplementedError( + "PropertyFieldIndexerWithDataPointerNoBoundsCheck does not support by_id." + ) + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.int64]]: + """Get index by ID. + + Parameters + ---------- + entity_id + """ idx = self._indices[entity_id] if idx < 0: return None return cast( - NDArray[np.double], + NDArray[np.int64], self._data[ self._data_pointer[idx] // self._n_components : self._data_pointer[idx + 1] @@ -252,39 +304,128 @@ def by_id(self, entity_id: int) -> Optional[NDArray[np.double]]: ) -class PropertyFieldIndexerWithDataPointerNoBoundsCheck: - """Indexer for a property field with data pointer and no bounds checks.""" +# DPF does not set the data pointers if a field has just +# one value per entity. Therefore, it is unknown if +# data pointer are set for some field of layered elements, for +# instance angles, thickness etc. +def get_field_indexer(field: Field) -> FieldIndexexProtocol: + """Get field indexer based on data pointer. - def __init__(self, field: PropertyField): + Parameters + ---------- + field + """ + if _has_data_pointer(field): + return FieldIndexerWithDataPointer(field) + return FieldIndexerNoDataPointer(field) + + +class FieldIndexerNoDataPointer: + """Indexer for a dpf field with no data pointer.""" + + def __init__(self, field: Field): """Create indexer and get data.""" if field.scoping.size > 0: index_by_id = setup_index_by_id(field.scoping) self._indices = index_by_id.mapping + self._max_id = index_by_id.max_id + self._data: NDArray[np.double] = np.array(field.data, dtype=np.double) + else: + self._indices = np.array([], dtype=np.int64) + self._max_id = 0 + self._data = np.array([], dtype=np.double) - self._data: NDArray[np.int64] = np.array(field.data, dtype=np.int64) + def by_id(self, entity_id: int) -> Optional[np.double]: + """Get value by ID. + + Parameters + ---------- + entity_id + """ + if entity_id > self._max_id: + return None + idx = self._indices[entity_id] + if idx < 0: + return None + return cast(np.double, self._data[idx]) + + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.double]]: + """Get values by id. + + Parameters + ---------- + entity_id + """ + value = self.by_id(entity_id) + if value is None: + return None + return np.array([value], dtype=np.double) + + +class FieldIndexerWithDataPointer: + """Indexer for a dpf field with data pointer.""" + + def __init__(self, field: Field): + """Create indexer and get data.""" + if field.scoping.size > 0: + index_by_id = setup_index_by_id(field.scoping) + self._indices = index_by_id.mapping + self._max_id = index_by_id.max_id + + self._data: NDArray[np.double] = np.array(field.data, dtype=np.double) self._n_components = field.component_count self._data_pointer: NDArray[np.int64] = np.append( field._data_pointer, len(self._data) * self._n_components ) else: + self._max_id = 0 self._indices = np.array([], dtype=np.int64) - self._data = np.array([], dtype=np.int64) + self._data = np.array([], dtype=np.double) self._n_components = 0 self._data_pointer = np.array([], dtype=np.int64) - def by_id(self, entity_id: int) -> Optional[NDArray[np.int64]]: - """Get index by ID. + def by_id(self, entity_id: int) -> Optional[np.double]: + """Get value by ID. Parameters ---------- entity_id """ + values = self.by_id_as_array(entity_id) + if values is None or len(values) == 0: + return None + if len(values) == 1: + return cast(np.double, values[0]) + + # There is an issue with the DPF server 2024r1_pre0 and before. + # Values of the laminate offset field does not have length 1. + # In this case the format of values is [offset, 0, 0., ...] + offset = values[0] + if all([v == 0 for v in values[1:]]): + return cast(np.double, offset) + + raise RuntimeError( + f"Cannot extract value for entity {entity_id}. " + "Use the latest version of the DPF server to get the correct value. " + f"Values: {values}" + ) + + def by_id_as_array(self, entity_id: int) -> Optional[NDArray[np.double]]: + """Get values by ID. + + Parameters + ---------- + entity_id + """ + if entity_id > self._max_id: + return None + idx = self._indices[entity_id] if idx < 0: return None return cast( - NDArray[np.int64], + NDArray[np.double], self._data[ self._data_pointer[idx] // self._n_components : self._data_pointer[idx + 1] diff --git a/src/ansys/dpf/composites/constants.py b/src/ansys/dpf/composites/constants.py index fbe20c507..4047bb867 100644 --- a/src/ansys/dpf/composites/constants.py +++ b/src/ansys/dpf/composites/constants.py @@ -30,10 +30,13 @@ "REF_SURFACE_NAME", "FAILURE_LABEL", "TIME_LABEL", + "TEMPERATURE_COMPONENT", ) FAILURE_LABEL = "failure_label" TIME_LABEL = "time" +TEMPERATURE_COMPONENT = 0 +REF_SURFACE_NAME = "Reference Surface" class Spot(IntEnum): @@ -69,6 +72,3 @@ class FailureOutput(IntEnum): MAX_GLOBAL_LAYER_IN_STACK = 5 MAX_LOCAL_LAYER_IN_ELEMENT = 6 MAX_SOLID_ELEMENT_ID = 7 - - -REF_SURFACE_NAME = "Reference Surface" diff --git a/src/ansys/dpf/composites/example_helper/__init__.py b/src/ansys/dpf/composites/example_helper/__init__.py index 9dc4c152a..8d85bb3a5 100644 --- a/src/ansys/dpf/composites/example_helper/__init__.py +++ b/src/ansys/dpf/composites/example_helper/__init__.py @@ -40,6 +40,10 @@ EXAMPLE_REPO = "https://github.com/ansys/example-data/raw/master/pydpf-composites/" +# Example URL to run the examples locally +# EXAMPLE_REPO = "file:////D:/Development/pyansys-example-data/pydpf-composites/" + + @dataclass class _ContinuousFiberCompositeFiles: definition: str @@ -62,7 +66,7 @@ class _ShortFiberCompositesExampleFilenames: @dataclass class _ContinuousFiberExampleLocation: - """Location of the a given continuous fiber example in the example_data repo. + """Location of a given continuous fiber example in the example_data repo. Parameters ---------- @@ -151,6 +155,16 @@ class _ShortFiberExampleLocation: }, ), ), + "thermal_solid": _ContinuousFiberExampleLocation( + directory="thermal_solid", + files=_ContinuousFiberCompositesExampleFilenames( + rst=["file.rst"], + engineering_data="MatML.xml", + composite={ + "shell": _ContinuousFiberCompositeFiles(definition="ACPSolidModel_SM.h5"), + }, + ), + ), } _short_fiber_examples: dict[str, _ShortFiberExampleLocation] = { diff --git a/src/ansys/dpf/composites/layup_info/_layup_info.py b/src/ansys/dpf/composites/layup_info/_layup_info.py index 44bf6b5b3..cbe1e215c 100644 --- a/src/ansys/dpf/composites/layup_info/_layup_info.py +++ b/src/ansys/dpf/composites/layup_info/_layup_info.py @@ -34,16 +34,7 @@ import numpy as np from numpy.typing import NDArray -from .._indexer import ( - FieldIndexerNoDataPointer, - FieldIndexerWithDataPointer, - PropertyFieldIndexerArrayValue, - PropertyFieldIndexerNoDataPointer, - PropertyFieldIndexerNoDataPointerNoBoundsCheck, - PropertyFieldIndexerSingleValue, - PropertyFieldIndexerWithDataPointer, - PropertyFieldIndexerWithDataPointerNoBoundsCheck, -) +from .._indexer import get_field_indexer, get_property_field_indexer from ..server_helpers import version_equal_or_later, version_older_than from ._enums import LayupProperty @@ -211,7 +202,7 @@ def __init__(self, mesh: MeshedRegion, name: str): """Initialize AnalysisPlyProvider.""" self.name = name self.property_field = _get_analysis_ply(mesh, name) - self._layer_indices = PropertyFieldIndexerNoDataPointer(self.property_field) + self._layer_indices = get_property_field_indexer(self.property_field, False) def get_layer_index_by_element_id(self, element_id: int) -> Optional[np.int64]: """Get the layer index for the analysis ply in a given element. @@ -404,31 +395,22 @@ def __init__( # focused on the most important properties. We can add different providers # for other properties (such as thickness and angles) - def get_indexer_with_data_pointer(field: PropertyField) -> PropertyFieldIndexerArrayValue: - if no_bounds_checks: - return PropertyFieldIndexerWithDataPointerNoBoundsCheck(field) - else: - return PropertyFieldIndexerWithDataPointer(field) - - def get_indexer_no_data_pointer(field: PropertyField) -> PropertyFieldIndexerSingleValue: - if no_bounds_checks: - return PropertyFieldIndexerNoDataPointerNoBoundsCheck(field) - else: - return PropertyFieldIndexerNoDataPointer(field) - # Has to be always with bounds checks because it does not contain # data for all the elements - self.layer_indices = PropertyFieldIndexerWithDataPointer(layer_indices) - self.layer_materials = get_indexer_with_data_pointer(material_ids) - self.apdl_element_types = get_indexer_no_data_pointer(element_types_apdl) - self.dpf_element_types = get_indexer_no_data_pointer(element_types_dpf) - self.keyopt_8_values = get_indexer_no_data_pointer(keyopt_8_values) - self.keyopt_3_values = get_indexer_no_data_pointer(keyopt_3_values) + self.layer_indices = get_property_field_indexer(layer_indices, no_bounds_checks) + self.layer_materials = get_property_field_indexer(material_ids, no_bounds_checks) + + self.apdl_element_types = get_property_field_indexer(element_types_apdl, no_bounds_checks) + self.dpf_element_types = get_property_field_indexer(element_types_dpf, no_bounds_checks) + self.keyopt_8_values = get_property_field_indexer(keyopt_8_values, no_bounds_checks) + self.keyopt_3_values = get_property_field_indexer(keyopt_3_values, no_bounds_checks) self.mesh = mesh self.corner_nodes_by_element_type = _get_corner_nodes_by_element_type_array() - self.apdl_material_indexer = get_indexer_no_data_pointer(self.mesh.elements.materials_field) + self.apdl_material_indexer = get_property_field_indexer( + self.mesh.elements.materials_field, no_bounds_checks + ) self.solver_material_to_dpf_id = {} if solver_material_ids is not None: @@ -473,9 +455,10 @@ def get_element_info(self, element_id: int) -> Optional[ElementInfo]: if element_type is None: raise IndexError(f"No DPF element type for element with id {element_id}.") - layer_data = self.layer_indices.by_id(element_id) + layer_data = self.layer_indices.by_id_as_array(element_id) if layer_data is not None: - dpf_material_ids = self.layer_materials.by_id(element_id) + # can be of type int for single layer elements or array for multilayer materials + dpf_material_ids = self.layer_materials.by_id_as_array(element_id) assert dpf_material_ids is not None assert layer_data[0] + 1 == len(layer_data), "Invalid size of layer data" n_layers = layer_data[0] @@ -605,24 +588,24 @@ def __init__(self, layup_provider: Operator, mesh: MeshedRegion): layup_outputs_container = layup_provider.outputs.section_data_container() composite_label = layup_outputs_container.labels[0] angle_field = layup_outputs_container.get_field({composite_label: LayupProperty.ANGLE}) - self._angle_indexer = FieldIndexerWithDataPointer(angle_field) + self._angle_indexer = get_field_indexer(angle_field) thickness_field = layup_outputs_container.get_field( {composite_label: LayupProperty.THICKNESS} ) - self._thickness_indexer = FieldIndexerWithDataPointer(thickness_field) + self._thickness_indexer = get_field_indexer(thickness_field) shear_angle_field = layup_outputs_container.get_field( {composite_label: LayupProperty.SHEAR_ANGLE} ) - self._shear_angle_indexer = FieldIndexerWithDataPointer(shear_angle_field) + self._shear_angle_indexer = get_field_indexer(shear_angle_field) offset_field = layup_outputs_container.get_field( {composite_label: LayupProperty.LAMINATE_OFFSET} ) - self._offset_indexer = FieldIndexerNoDataPointer(offset_field) + self._offset_indexer = get_field_indexer(offset_field) self._index_to_name_map = get_analysis_ply_index_to_name_map(mesh) - self._analysis_ply_indexer = PropertyFieldIndexerWithDataPointer( - mesh.property_field("layer_to_analysis_ply") + self._analysis_ply_indexer = get_property_field_indexer( + mesh.property_field("layer_to_analysis_ply"), False ) def get_layer_angles(self, element_id: int) -> Optional[NDArray[np.double]]: @@ -633,7 +616,7 @@ def get_layer_angles(self, element_id: int) -> Optional[NDArray[np.double]]: element_id: Element Id/Label """ - return self._angle_indexer.by_id(element_id) + return self._angle_indexer.by_id_as_array(element_id) def get_layer_thicknesses(self, element_id: int) -> Optional[NDArray[np.double]]: """Get thicknesses for all layers. Returns None if element is not layered. @@ -644,7 +627,7 @@ def get_layer_thicknesses(self, element_id: int) -> Optional[NDArray[np.double]] Element Id/Label """ - return self._thickness_indexer.by_id(element_id) + return self._thickness_indexer.by_id_as_array(element_id) def get_layer_shear_angles(self, element_id: int) -> Optional[NDArray[np.double]]: """Get shear angle for all layers. Returns None if element is not layered. @@ -654,7 +637,7 @@ def get_layer_shear_angles(self, element_id: int) -> Optional[NDArray[np.double] element_id: Element Id/Label """ - return self._shear_angle_indexer.by_id(element_id) + return self._shear_angle_indexer.by_id_as_array(element_id) def get_element_laminate_offset(self, element_id: int) -> Optional[np.double]: """Get laminate offset of element. Returns None if element is not layered. @@ -676,7 +659,7 @@ def get_analysis_plies(self, element_id: int) -> Optional[Sequence[str]]: Element Id/Label """ - indexes = self._analysis_ply_indexer.by_id(element_id) - if indexes is None: + indices = self._analysis_ply_indexer.by_id_as_array(element_id) + if indices is None: return None - return [self._index_to_name_map[index] for index in indexes] + return [self._index_to_name_map[index] for index in indices] diff --git a/src/ansys/dpf/composites/select_indices.py b/src/ansys/dpf/composites/select_indices.py index 2e4c8798e..28ea98148 100644 --- a/src/ansys/dpf/composites/select_indices.py +++ b/src/ansys/dpf/composites/select_indices.py @@ -142,10 +142,17 @@ def get_selected_indices( ) # Todo: Use numpy. Probably use ravel_multi_index method. current_index = 0 + + num_nodes_per_spot = ( + element_info.number_of_nodes_per_spot_plane + if element_info.is_layered + else element_info.n_corner_nodes + ) + for layer_index in layer_indices: - layer_start_index = layer_index * element_info.n_corner_nodes * element_info.n_spots + layer_start_index = layer_index * num_nodes_per_spot * element_info.n_spots for spot_index in spot_indices: - spot_start_index = layer_start_index + spot_index * element_info.n_corner_nodes + spot_start_index = layer_start_index + spot_index * num_nodes_per_spot for corner_index in node_indices: all_indices[current_index] = spot_start_index + corner_index current_index = current_index + 1 diff --git a/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer.rst b/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer.rst new file mode 100644 index 000000000..9173fe602 Binary files /dev/null and b/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer.rst differ diff --git a/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer_material.xml b/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer_material.xml new file mode 100644 index 000000000..815ab7248 --- /dev/null +++ b/tests/data/all_element_types/model_with_all_element_types_all_output_1_layer_material.xml @@ -0,0 +1,107 @@ + + + + + + + 1 + + - + ACP + Isotropic + + + - + Isotropic + + 70000 + + + 0.30000000000000004 + + + 26923.076923076922 + + + + + + + 2 + + - + ACP + Isotropic + + + - + Isotropic + + 210000 + + + 0.30000000000000004 + + + 80769.230769230766 + + + + + + + Ply Type + + + + Elasticity + + + + Young's Modulus + + + [Length] + + + [Mass] + + + [Time] + + + + + Poisson's Ratio + + + + Shear Modulus + + + [Length] + + + [Mass] + + + [Time] + + + + + + + + + + 1 + 2d0cf88a-4ac7-4d14-b949-21b9f3a3f932 + + + 2 + af1cc58e-261a-40eb-b06c-95c53f017823 + + + + diff --git a/tests/element_info_output_all_element_types_test.py b/tests/element_info_output_all_element_types_test.py index dc9cdf336..921e36e0b 100644 --- a/tests/element_info_output_all_element_types_test.py +++ b/tests/element_info_output_all_element_types_test.py @@ -25,16 +25,19 @@ import pathlib import ansys.dpf.core as dpf -from ansys.dpf.core import Field, MeshedRegion, PropertyField +from ansys.dpf.core import Field, MeshedRegion, PropertyField, unit_systems +import numpy as np import pytest +from ansys.dpf.composites.composite_model import CompositeModel from ansys.dpf.composites.constants import Spot +from ansys.dpf.composites.data_sources import ContinuousFiberCompositesFiles from ansys.dpf.composites.layup_info import ElementInfo, get_element_info_provider from ansys.dpf.composites.select_indices import ( get_selected_indices, get_selected_indices_by_dpf_material_ids, ) -from ansys.dpf.composites.server_helpers import upload_file_to_unique_tmp_folder +from ansys.dpf.composites.server_helpers import upload_file_to_unique_tmp_folder, version_older_than @dataclass(frozen=True) @@ -47,6 +50,7 @@ class ExpectedOutput: # otherwise it is 2 or 3 for top/bot or top/bot/mid output n_spots: int is_layered: bool + number_of_nodes_per_spot: int @dataclass(frozen=True) @@ -111,25 +115,25 @@ def get_expected_output( n_spots_homogeneous_solid = 0 return { - 1: ExpectedOutput(3, 4, 181, n_spots_shell, True), - 2: ExpectedOutput(3, 3, 181, n_spots_shell, True), - 3: ExpectedOutput(3, 4, 281, n_spots_shell, True), - 4: ExpectedOutput(3, 3, 281, n_spots_shell, True), - 10: ExpectedOutput(1, 8, 185, n_spots_homogeneous_solid, False), - 11: ExpectedOutput(1, 6, 185, n_spots_homogeneous_solid, False), - 12: ExpectedOutput(1, 5, 185, n_spots_homogeneous_solid, False), - 13: ExpectedOutput(1, 4, 185, n_spots_homogeneous_solid, False), - 20: ExpectedOutput(1, 8, 186, n_spots_homogeneous_solid, False), - 21: ExpectedOutput(1, 6, 186, n_spots_homogeneous_solid, False), - 22: ExpectedOutput(1, 5, 186, n_spots_homogeneous_solid, False), - 23: ExpectedOutput(1, 4, 186, n_spots_homogeneous_solid, False), - 24: ExpectedOutput(1, 4, 187, n_spots_homogeneous_solid, False), - 30: ExpectedOutput(3, 8, 185, n_spots_layered_solid, True), - 31: ExpectedOutput(3, 6, 185, n_spots_layered_solid, True), - 40: ExpectedOutput(3, 8, 186, n_spots_layered_solid, True), - 41: ExpectedOutput(3, 6, 186, n_spots_layered_solid, True), - 50: ExpectedOutput(3, 8, 190, n_spots_layered_solid, True), - 51: ExpectedOutput(3, 6, 190, n_spots_layered_solid, True), + 1: ExpectedOutput(3, 4, 181, n_spots_shell, True, 4), + 2: ExpectedOutput(3, 3, 181, n_spots_shell, True, 3), + 3: ExpectedOutput(3, 4, 281, n_spots_shell, True, 4), + 4: ExpectedOutput(3, 3, 281, n_spots_shell, True, 3), + 10: ExpectedOutput(1, 8, 185, n_spots_homogeneous_solid, False, -1), + 11: ExpectedOutput(1, 6, 185, n_spots_homogeneous_solid, False, -1), + 12: ExpectedOutput(1, 5, 185, n_spots_homogeneous_solid, False, -1), + 13: ExpectedOutput(1, 4, 185, n_spots_homogeneous_solid, False, -1), + 20: ExpectedOutput(1, 8, 186, n_spots_homogeneous_solid, False, -1), + 21: ExpectedOutput(1, 6, 186, n_spots_homogeneous_solid, False, -1), + 22: ExpectedOutput(1, 5, 186, n_spots_homogeneous_solid, False, -1), + 23: ExpectedOutput(1, 4, 186, n_spots_homogeneous_solid, False, -1), + 24: ExpectedOutput(1, 4, 187, n_spots_homogeneous_solid, False, -1), + 30: ExpectedOutput(3, 8, 185, n_spots_layered_solid, True, 4), + 31: ExpectedOutput(3, 6, 185, n_spots_layered_solid, True, 3), + 40: ExpectedOutput(3, 8, 186, n_spots_layered_solid, True, 4), + 41: ExpectedOutput(3, 6, 186, n_spots_layered_solid, True, 3), + 50: ExpectedOutput(3, 8, 190, n_spots_layered_solid, True, 4), + 51: ExpectedOutput(3, 6, 190, n_spots_layered_solid, True, 3), } def check_output(rst_file, expected_output): @@ -184,6 +188,8 @@ def check_output(rst_file, expected_output): corner_nodes_per_layer = element_info.n_corner_nodes else: corner_nodes_per_layer = element_info.n_corner_nodes / 2 + + assert corner_nodes_per_layer == element_info.number_of_nodes_per_spot_plane assert ( num_elementary_data == corner_nodes_per_layer * element_info.n_spots * element_info.n_layers @@ -205,32 +211,32 @@ def check_output(rst_file, expected_output): ) -def test_document_error_cases_indices(dpf_server): +def get_element_info_provider_for_rst(rst_file, server): TEST_DATA_ROOT_DIR = pathlib.Path(__file__).parent / "data" + rst_path = TEST_DATA_ROOT_DIR / "all_element_types" / rst_file - def get_layup_info_for_rst(rst_file): - rst_path = TEST_DATA_ROOT_DIR / "all_element_types" / rst_file + if not server.local_server: + rst_path = upload_file_to_unique_tmp_folder(rst_path, server=server) - if not dpf_server.local_server: - rst_path = upload_file_to_unique_tmp_folder(rst_path, server=dpf_server) + rst_data_source = dpf.DataSources(rst_path) - rst_data_source = dpf.DataSources(rst_path) + mesh_provider = dpf.Operator("MeshProvider") + mesh_provider.inputs.data_sources(rst_data_source) + mesh: MeshedRegion = mesh_provider.outputs.mesh() - mesh_provider = dpf.Operator("MeshProvider") - mesh_provider.inputs.data_sources(rst_data_source) - mesh: MeshedRegion = mesh_provider.outputs.mesh() + with pytest.raises(RuntimeError) as exc_info: + layup_info = get_element_info_provider(mesh, stream_provider_or_data_source=rst_data_source) + assert str(exc_info.value).startswith("Missing property field in mesh") + material_property_field, layer_indices_property_field = get_layup_property_fields() + mesh.set_property_field("element_layered_material_ids", material_property_field) + mesh.set_property_field("element_layer_indices", layer_indices_property_field) + return get_element_info_provider(mesh, stream_provider_or_data_source=rst_data_source) - with pytest.raises(RuntimeError) as exc_info: - layup_info = get_element_info_provider( - mesh, stream_provider_or_data_source=rst_data_source - ) - assert str(exc_info.value).startswith("Missing property field in mesh") - material_property_field, layer_indices_property_field = get_layup_property_fields() - mesh.set_property_field("element_layered_material_ids", material_property_field) - mesh.set_property_field("element_layer_indices", layer_indices_property_field) - return get_element_info_provider(mesh, stream_provider_or_data_source=rst_data_source) - layup_info = get_layup_info_for_rst("model_with_all_element_types_minimal_output.rst") +def test_document_error_cases_indices(dpf_server): + layup_info = get_element_info_provider_for_rst( + "model_with_all_element_types_minimal_output.rst", dpf_server + ) for element_id in get_element_ids().layered: with pytest.raises(RuntimeError) as exc_info: @@ -248,7 +254,9 @@ def get_layup_info_for_rst(rst_file): "Computation of indices is not supported for non-layered elements." ) - layup_info = get_layup_info_for_rst("model_with_all_element_types_all_output.rst") + layup_info = get_element_info_provider_for_rst( + "model_with_all_element_types_all_output.rst", dpf_server + ) for element_id in get_element_ids().non_layered: with pytest.raises(RuntimeError) as exc_info: @@ -279,7 +287,9 @@ def get_layup_info_for_rst(rst_file): selected_indices = get_selected_indices_by_dpf_material_ids(element_info, [5]) assert len(selected_indices) == 0 - layup_info = get_layup_info_for_rst("model_with_all_element_types_all_except_mid_output.rst") + layup_info = get_element_info_provider_for_rst( + "model_with_all_element_types_all_except_mid_output.rst", dpf_server + ) for element_id in get_element_ids().layered: with pytest.raises(RuntimeError) as exc_info: @@ -288,3 +298,112 @@ def get_layup_info_for_rst(rst_file): assert str(exc_info.value).startswith( "Spot index 2 is greater or equal to the number of spots" ) + + +def test_select_indices_all_element_types(dpf_server): + """ + Test get_selected_indices for all types of layered elements. + + The test verifies the indices for the first layer, 2nd layer, and + 2nd layer in combination with spot TOP. + + Note: Non-layered elements are not supported by get_selected_indices + """ + ref_indices_layer_0 = { + 1: np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]), # 4 node shell181, 3 layers, 3 spots + 2: np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]), # 3 node shell181, 3 layers, 2 spots + 3: np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]), # 8 node shell281, 3 layers, 3 spots + 4: np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]), # 6 node shell281, 3 layers, 2 spots + 30: np.array([0, 1, 2, 3, 4, 5, 6, 7]), # 8 node solid185 + 31: np.array([0, 1, 2, 3, 4, 5]), # 6 node solid185 + 40: np.array([0, 1, 2, 3, 4, 5, 6, 7]), # 20 node solid186 + 41: np.array([0, 1, 2, 3, 4, 5]), # 15 node solid186 + 50: np.array([0, 1, 2, 3, 4, 5, 6, 7]), # 8 node solid190 + 51: np.array([0, 1, 2, 3, 4, 5]), # 6 node solid190 + } + + element_info_provider = get_element_info_provider_for_rst( + "model_with_all_element_types_all_output.rst", dpf_server + ) + element_ids = get_element_ids() + for elem_id in element_ids.all: + element_info = element_info_provider.get_element_info(elem_id) + # get_selected_indices is only supported for layered elements + if element_info.is_layered: + # All indices of the first layer + indices = get_selected_indices(element_info, layers=[0]) + assert ( + indices == ref_indices_layer_0[elem_id] + ).all(), f"{element_info}, {indices} != {ref_indices_layer_0[elem_id]}" + + # All indices of the second layer via it's material ID + material_id = element_info.dpf_material_ids[1] # this is equivalent to the second layer + indices = get_selected_indices_by_dpf_material_ids(element_info, list([material_id])) + + # Offset indices for the second layer + ref_2nd_layer = ref_indices_layer_0[elem_id] + max(ref_indices_layer_0[elem_id]) + 1 + assert ( + indices == ref_2nd_layer + ).all(), f"{element_info}, i{indices} != {ref_2nd_layer}" + + # Indices of the second layer and the top spot + indices = get_selected_indices(element_info, layers=[1], spots=[Spot.TOP]) + num_indices = element_info.number_of_nodes_per_spot_plane + if element_info.is_shell: + # The order of the spots is bot, top, middle for shells + # So the indices in the middle are used to retrieve the data at the top + ref_2nd_layer_top = ref_2nd_layer[num_indices:-num_indices] + else: + # Layered solids have only bottom and top. + ref_2nd_layer_top = ref_2nd_layer[-num_indices:] + assert ( + indices == ref_2nd_layer_top + ).all(), f"{element_info}, {indices} != {ref_2nd_layer_top}" + + +def test_get_element_info_all_element_types(dpf_server): + """ + Test get_element_info for all element types. + + The layered elements have one layer only. In this case, the dpf fields do not have data + pointers. In addition, the analysis_ply_layer_indices property field is not available + because the model is loaded from an RST file without a lay-up definition file. + """ + + if version_older_than(dpf_server, "8.0"): + pytest.xfail( + "Not supported because section data from RST is not implemented before version 8.0." + ) + + model_path = pathlib.Path(__file__).parent / "data" / "all_element_types" + + rst_file = model_path / "model_with_all_element_types_all_output_1_layer.rst" + mat_xml_file = model_path / "model_with_all_element_types_all_output_1_layer_material.xml" + + composite_files = ContinuousFiberCompositesFiles( + rst=[rst_file], + composite={}, + engineering_data=mat_xml_file, + ) + + composite_model = CompositeModel( + composite_files, server=dpf_server, default_unit_system=unit_systems.solver_mks + ) + + expected_indices = { + 1: np.array([4, 5, 6, 7], dtype=np.int64), + 2: np.array([3, 4, 5], dtype=np.int64), + 3: np.array([4, 5, 6, 7], dtype=np.int64), + 4: np.array([3, 4, 5], dtype=np.int64), + 30: np.array([4, 5, 6, 7], dtype=np.int64), + 31: np.array([3, 4, 5], dtype=np.int64), + 40: np.array([4, 5, 6, 7], dtype=np.int64), + 41: np.array([3, 4, 5], dtype=np.int64), + 50: np.array([4, 5, 6, 7], dtype=np.int64), + 51: np.array([3, 4, 5], dtype=np.int64), + } + + for element_id, ref_indices in expected_indices.items(): + element_info = composite_model.get_element_info(element_id) + indices = get_selected_indices(element_info, layers=[0], spots=[Spot.TOP]) + assert (indices == ref_indices).all(), f"{element_info}, {indices} != {ref_indices}" diff --git a/tests/performance_test.py b/tests/performance_test.py index bd53b06a6..b8d62cc30 100644 --- a/tests/performance_test.py +++ b/tests/performance_test.py @@ -27,7 +27,7 @@ import numpy as np import pytest -from ansys.dpf.composites._indexer import FieldIndexerWithDataPointer +from ansys.dpf.composites._indexer import get_field_indexer from ansys.dpf.composites.data_sources import CompositeDefinitionFiles, get_composites_data_sources from ansys.dpf.composites.layup_info import ( LayupPropertiesProvider, @@ -335,7 +335,7 @@ def test_performance_flat(dpf_server): all_data = np.full((setup_result.field.elementary_data_count, 11), -1.0) start_index = 0 - indexer_data = FieldIndexerWithDataPointer(setup_result.field) + indexer_data = get_field_indexer(setup_result.field) timer.add("indexer") with setup_result.mesh.elements.connectivities_field.as_local_field() as local_connectivity: @@ -368,7 +368,7 @@ def test_performance_flat(dpf_server): flat_layer_indices = unraveled_indices[0] flat_spot_indices = unraveled_indices[1] flat_node_indices = unraveled_indices[2] - element_data = indexer_data.by_id(element_id) + element_data = indexer_data.by_id_as_array(element_id) nodes = np.array(local_connectivity.get_entity_data_by_id(element_id)) num_elementary_data = ( diff --git a/tests/test_metadata.py b/tests/test_metadata.py index 37ff358d7..c989b27a7 100644 --- a/tests/test_metadata.py +++ b/tests/test_metadata.py @@ -24,4 +24,4 @@ def test_pkg_version(): - assert __version__ == "0.6.0" + assert __version__ == "0.6.1"