Skip to content

Isosurface example #1056

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 15 commits into from
Jul 26, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 113 additions & 0 deletions examples/12-fluids/03-fluids_isosurface.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
"""
.. _ref_fluids_isosurface:

Compute iso-surfaces on fluid models
------------------------------------------

This example demonstrates how to compute iso-surfaces on fluid models.
"""

###############################################################################
# Import the ``dpf-core`` module and its examples files.
# ~~~~~~~~~~~~~~~~~~

import ansys.dpf.core as dpf
from ansys.dpf.core import examples
from ansys.dpf.core.plotter import DpfPlotter

###############################################################################
# Specify the file path.
# ~~~~~~~~~~~~~~~~~~
# We work on a cas/dat.h5 file with only nodal variables.

path = examples.download_cfx_heating_coil()
ds = dpf.DataSources()
ds.set_result_file_path(path["cas"], "cas")
ds.add_file_path(path["dat"], "dat")
streams = dpf.operators.metadata.streams_provider(data_sources=ds)

###############################################################################
# Whole mesh scoping.
# ~~~~~~~~~~~~~~~~~~
# We evaluate the mesh with the mesh_provider operator to scope the mesh_cut operator
# with the whole mesh.

whole_mesh = dpf.operators.mesh.mesh_provider(streams_container=streams).eval()
print(whole_mesh)

pl = DpfPlotter()
pl.add_mesh(whole_mesh)
cpos_whole_mesh = [
(4.256160478475664, 4.73662111240005, 4.00410065817644),
(-0.0011924505233764648, 1.8596649169921875e-05, 1.125),
(-0.2738679385987956, -0.30771426079547065, 0.9112125360807675),
]
pl.show_figure(cpos=cpos_whole_mesh, show_axes=True)

###############################################################################
# Extract the physics variable
# ~~~~~~~~~~~~~~~~~
# Here we choose to work with the static pressure by default which is a scalar and
# nodal variable without multi-species/phases. With a multi-species case,
# select one using qualifier ellipsis pins and connecting a LabelSpace "species"/"phase".

P_S = dpf.operators.result.static_pressure(streams_container=streams, mesh=whole_mesh).eval()
print(P_S[0])

pl = DpfPlotter()
pl.add_field(P_S[0])
cpos_mesh_variable = [
(4.256160478475664, 4.73662111240005, 4.00410065817644),
(-0.0011924505233764648, 1.8596649169921875e-05, 1.125),
(-0.2738679385987956, -0.30771426079547065, 0.9112125360807675),
]
pl.show_figure(cpos=cpos_mesh_variable, show_axes=True)

###############################################################################
# Evaluate iso-surfaces
# ~~~~~~~~~~~~~~
# We can finally use the mesh_cut operator on this specific variable.
# We choose to cut the whole with 5 iso-surface equally spaced between min and max.

max_pressure = 361.8170 # Pa
min_pressure = -153.5356 # Pa
number_of_iso_surface = 5
step = (max_pressure - min_pressure) / number_of_iso_surface

pl = DpfPlotter()
c_pos_iso = [
(4.256160478475664, 4.73662111240005, 4.00410065817644),
(-0.0011924505233764648, 1.8596649169921875e-05, 1.125),
(-0.2738679385987956, -0.30771426079547065, 0.9112125360807675),
]
pl.add_mesh(
meshed_region=whole_mesh,
style="wireframe",
show_edges=True,
show_axes=True,
color="black",
opacity=0.3,
)

for i in range(number_of_iso_surface):
iso_surface = dpf.operators.mesh.mesh_cut(
field=P_S[0], iso_value=min_pressure, closed_surface=0, mesh=whole_mesh, slice_surfaces=True
).eval()
P_S_step = dpf.Field(location=dpf.locations.overall, nature=dpf.common.natures.scalar)
P_S_step.append([min_pressure], i)
P_S_step.name = "static pressure"
P_S_step.unit = "Pa"
pl.add_field(
field=P_S_step, meshed_region=iso_surface, style="surface", show_edges=False, show_axes=True
)
min_pressure += step

pl.show_figure(show_axes=True, cpos=c_pos_iso)

###############################################################################
# Important note
# ------------------------------
# Iso-surfaces computation through the `mesh_cut` operator are only supported for Nodal Fields.
# For Elemental variables, you must perform an averaging operation on the Nodes before
# running the `mesh_cut` operator. This can be done by chaining the `elemental_to_nodal` operator
# output with the `mesh_cut` operator input.
23 changes: 22 additions & 1 deletion src/ansys/dpf/core/field.py
Original file line number Diff line number Diff line change
Expand Up @@ -575,7 +575,6 @@ def dimensionality(self, value):
@property
def name(self):
"""Name of the field."""
# return self._api.csfield_get_name(self)
from ansys.dpf.gate import integral_types

size = integral_types.MutableInt32()
Expand All @@ -585,6 +584,28 @@ def name(self):
)
return str(name)

@name.setter
def name(self, value):
"""Change the name of the field

Parameters
----------
value : str
Name of the field.

Examples
--------
Units for a displacement field.

>>> from ansys.dpf import core as dpf
>>> my_field = dpf.Field(10, dpf.natures.vector,dpf.locations.nodal)
>>> my_field.name = "my-field"
>>> my_field.name
'my-field'

"""
self._field_definition._api.csfield_definition_set_name(self._field_definition, name=value)

def _set_field_definition(self, field_definition):
"""Set the field definition.

Expand Down
51 changes: 51 additions & 0 deletions src/ansys/dpf/core/fields_factory.py
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,57 @@ def create_vector_field(num_entities, num_comp, location=locations.nodal, server
return _create_field(server, natures.vector, num_entities, location, ncomp_n=num_comp)


def create_overall_field(
value, nature, num_entities, num_comp, location=locations.overall, server=None
):
"""Create a specific `:class:`ansys.dpf.core.Field` with entities that have an
overall location.

Regarding the nature of the entity contained in the field, we set the same value
for all elements.

Parameters
----------
value : float
Value of the entity
nature : str
Nature of the field entity data. For example:

- :class:`ansys.dpf.core.natures.matrix`
- :class:`ansys.dpf.core.natures.scalar`
num_entities : int
Number of entities to reserve.
num_comp : int
Number of vector components.
location : str, optional
Location of the field. Options are in :class:`locations <ansys.dpf.core.common.locations>`.
The default is ``dpf.locations.nodal``.

server : ansys.dpf.core.server, optional
Server with the channel connected to the remote or local instance.
The default is ``None``, in which case an attempt is made to use the
global server.

Returns
-------
field : Field
DPF field in the requested format.

Examples
--------
Create a field containing 10 scalar entities of 1 component each with an
overall location (default). Same value (1.0) is set for all element of the field.

>>> from ansys.dpf.core import fields_factory
>>> field = fields_factory.create_overall_field(1.0, natures.scalar, 10, 1)

"""
overall_field = _create_field(server, nature, num_entities, location, ncomp_n=num_comp)
for i in range(num_entities):
overall_field.append(value, i)
return overall_field


def _create_field(server, nature, nentities, location=locations.nodal, ncomp_n=0, ncomp_m=0):
"""Create a specific :class:`ansys.dpf.core.Field`.

Expand Down
12 changes: 9 additions & 3 deletions src/ansys/dpf/core/plotter.py
Original file line number Diff line number Diff line change
Expand Up @@ -269,20 +269,26 @@ def add_field(
show_min = False
elif location == locations.faces:
mesh_location = meshed_region.faces
if len(mesh_location) == 0:
raise ValueError("No faces found to plot on")
if show_max or show_min:
warnings.warn("`show_max` and `show_min` is only supported for Nodal results.")
show_max = False
show_min = False
elif location == locations.overall:
mesh_location = meshed_region.elements
else:
raise ValueError("Only elemental, nodal or faces location are supported for plotting.")
component_count = field.component_count
if component_count > 1:
overall_data = np.full((len(mesh_location), component_count), np.nan)
else:
overall_data = np.full(len(mesh_location), np.nan)
ind, mask = mesh_location.map_scoping(field.scoping)
overall_data[ind] = field.data[mask]

if location != locations.overall:
ind, mask = mesh_location.map_scoping(field.scoping)
overall_data[ind] = field.data[mask]
else:
overall_data[:] = field.data[0]
# Filter kwargs for add_mesh
kwargs_in = _sort_supported_kwargs(bound_method=self._plotter.add_mesh, **kwargs)
# Have to remove any active scalar field from the pre-existing grid object,
Expand Down