Fix function expression evaluation

- Fixed the function expression evaluation.
eval is unsafe, and ne.evaluate does not recognize
the heaviside function -> convert to where / np.where.

Author: dragos-ana

src/fourc_webviewer/fourc_webserver.py

@@ -23,16 +23,15 @@
 from fourc_webviewer.input_file_utils.io_utils import (
     create_file_object_for_browser,
     get_master_and_linked_material_indices,
+    get_variable_data_by_name_in_funct_item,
     read_fourc_yaml_file,
     write_fourc_yaml_file,
-    get_variable_data_by_name_in_funct_item,
 )
 from fourc_webviewer.python_utils import convert_string2number, find_value_recursively
 from fourc_webviewer.read_geometry_from_file import (
     FourCGeometry,
 )
 
-
 # always set pyvista to plot off screen with Trame
 pv.OFF_SCREEN = True
 
@@ -266,15 +265,15 @@ def update_pyvista_render_objects(self, init_rendering=False):
         )
 
         # get coords of node with prescribed result description
-        self._server_vars[
-            "pv_selected_result_description_node_coords"
-        ] = self._server_vars["pv_mesh"].points[
-            self.state.result_description_section[
-                self.state.selected_result_description_id
-            ]["PARAMETERS"]["NODE"]
-            - 1,
-            :,
-        ]
+        self._server_vars["pv_selected_result_description_node_coords"] = (
+            self._server_vars["pv_mesh"].points[
+                self.state.result_description_section[
+                    self.state.selected_result_description_id
+                ]["PARAMETERS"]["NODE"]
+                - 1,
+                :,
+            ]
+        )
 
         # update plotter / rendering
         pv_render.update_pv_plotter(
@@ -352,9 +351,9 @@ def init_general_sections_state_and_server_vars(self):
                         self.state.general_sections[main_section_name] = {}
 
                     # add subsection
-                    self.state.general_sections[main_section_name][
-                        section_name
-                    ] = section_data
+                    self.state.general_sections[main_section_name][section_name] = (
+                        section_data
+                    )
 
     def sync_general_sections_from_state(self):
         """Syncs the server-side general sections based on the current values
@@ -436,9 +435,9 @@ def init_materials_state_and_server_vars(self):
             ):
                 if mat_id in linked_material_indices_item:
                     # add linked material indices
-                    mat_item_val["RELATIONSHIPS"][
-                        "LINKED MATERIALS"
-                    ] = linked_material_indices_item
+                    mat_item_val["RELATIONSHIPS"]["LINKED MATERIALS"] = (
+                        linked_material_indices_item
+                    )
 
                     # add master material index
                     mat_item_val["RELATIONSHIPS"]["MASTER MATERIAL"] = material_indices[
@@ -498,9 +497,9 @@ def sync_materials_sections_from_state(self):
         # write to server-side content
         self._server_vars["fourc_yaml_content"]["MATERIALS"] = new_materials_section
         if new_cloning_material_map_section:
-            self._server_vars["fourc_yaml_content"][
-                "CLONING MATERIAL MAP"
-            ] = new_cloning_material_map_section
+            self._server_vars["fourc_yaml_content"]["CLONING MATERIAL MAP"] = (
+                new_cloning_material_map_section
+            )
 
     def init_design_conditions_state_and_server_vars(self):
         """Initialize the state and server variables for the design condition
@@ -692,9 +691,9 @@ def sync_result_description_section_from_state(self):
             new_result_description_section.append({field: params})
 
         # set result description section on the server
-        self._server_vars["fourc_yaml_content"][
-            "RESULT DESCRIPTION"
-        ] = new_result_description_section
+        self._server_vars["fourc_yaml_content"]["RESULT DESCRIPTION"] = (
+            new_result_description_section
+        )
 
     def init_funct_state_and_server_vars(self):
         """Initialize the state and server variables for the function

src/fourc_webviewer/gui_utils.py

@@ -540,7 +540,9 @@ def _functions_panel(server):
                         display_mode_bar="true",
                     )
                     server.controller.figure_update = figure.update
-                    if server.state.funct_section[server.state.selected_funct][
+                    if server.state.funct_section[
+                        server.state.selected_funct
+                    ][
                         server.state.selected_funct_item
                     ][
                         "VISUALIZATION"

src/fourc_webviewer/input_file_utils/fourc_yaml_file_visualization.py

@@ -8,17 +8,18 @@
 import numexpr as ne
 import numpy as np
 import plotly.express as px
+
 from fourc_webviewer.input_file_utils.io_utils import (
     get_variable_data_by_name_in_funct_item,
 )
 
-
 # functional expressions / constants known by 4C, that are replaced by the numpy counterpart during evaluation
 DEF_FUNCT = ["exp", "sqrt", "log", "sin", "cos", "tan", "heaviside", "pi"]
 
 
 def get_variable_names_in_funct_expression(funct_expression: str):
-    """Returns all variable names present in a functional expression, using regular expressions."""
+    """Returns all variable names present in a functional expression, using
+    regular expressions."""
     vars_found = re.findall(r"[A-Za-z_]+", funct_expression)
     return [
         v for v in vars_found if v not in DEF_FUNCT and v not in ["t", "x", "y", "z"]
@@ -159,48 +160,54 @@ def funct_using_eval(x, y, z, t):
         Returns:
             parsed object using ast.literal_eval
         """
-        # funct_string copy
+        # Create a safe environment
+        safe_dict = {
+            "x": x,
+            "y": y,
+            "z": z,
+            "t": t,
+            "sin": np.sin,
+            "cos": np.cos,
+            "exp": np.exp,
+            "log": np.log,
+            "sqrt": np.sqrt,
+            "where": np.where,
+            "pi": np.pi,
+            # "heaviside": np.heaviside, -> no heaviside, numexpr will
+            # not deal with this
+            # add other safe functions as needed
+        }
+
         funct_string_copy = funct_string
 
-        # replace variables by their functional expressions
+        # replace variables by their
         for k, v in variable_funct_strings.items():
             funct_string_copy = re.sub(
                 rf"(?<![A-Za-z]){k}(?![A-Za-z])", v, funct_string_copy
             )
 
-        # replace the defined functions in the funct_string with "<DEF_FUNCT>"
-        for i in range(len(DEF_FUNCT)):
-            funct_string_copy = funct_string_copy.replace(
-                DEF_FUNCT[i], f"np.{DEF_FUNCT[i]}"
-            )
-
-        # replace the used power sign
+        # replace heaviside functions with where / np.where
         funct_string_copy = funct_string_copy.replace("^", "**")
-
-        # replace variables
-        funct_string_copy = (
-            funct_string_copy.replace("x", str(x))
-            .replace("y", str(y))
-            .replace("z", str(z))
-            .replace("t", str(t))
+        funct_string_copy = re.sub(
+            r"heaviside\(([^),]+)\)", r"where(\1 >= 0, 1, 0)", funct_string_copy
         )
-
-        # for heaviside: np.heaviside takes two arguments -> second argument denotes the function value at the first argument -> we set it by default to 0
         funct_string_copy = re.sub(
-            r"heaviside\((.*?)\)", r"heaviside(\1,0)", funct_string_copy
-        )  # usage of raw strings, (.*?) is a non greedy capturing, and \1 replaces the captured value
+            r"heaviside\(([^),]+),\s*([^)]+)\)",
+            r"where(\1 > 0, 1, where(\1 == 0, \2, 0))",
+            funct_string_copy,
+        )
 
-        return eval(
-            funct_string_copy, {"np": np}, {}
-        )  # this parses string in as a function
+        # Numexpr evaluation (much safer)
+        return ne.evaluate(funct_string_copy, local_dict=safe_dict)
 
     return np.frompyfunc(funct_using_eval, 4, 1)
 
 
 def construct_funct_string_from_variable_data(
     variable_name: str, funct_section_item: dict
 ):
-    """Constructs a functional string from the given data for a function variable."""
+    """Constructs a functional string from the given data for a function
+    variable."""
 
     # retrieve variable data
     variable_data = get_variable_data_by_name_in_funct_item(
@@ -212,8 +219,9 @@ def construct_funct_string_from_variable_data(
     match variable_data["TYPE"]:
         case "linearinterpolation":
             # get times and values
-            times, values = np.array(variable_data["TIMES"]), np.array(
-                variable_data["VALUES"]
+            times, values = (
+                np.array(variable_data["TIMES"]),
+                np.array(variable_data["VALUES"]),
             )
 
             # consistency check: time should start with 0.0
@@ -226,7 +234,7 @@ def construct_funct_string_from_variable_data(
                 if time_instant_index != 0:
                     funct_string += "+"
 
-                funct_string += f"({values[time_instant_index]}+({values[time_instant_index+1]}-{values[time_instant_index]})/({times[time_instant_index+1]}-{time_instant})*(t-{time_instant}))*heaviside(t-{time_instant})*heaviside({times[time_instant_index+1]}-t)"
+                funct_string += f"({values[time_instant_index]}+({values[time_instant_index + 1]}-{values[time_instant_index]})/({times[time_instant_index + 1]}-{time_instant})*(t-{time_instant}))*heaviside(t-{time_instant})*heaviside({times[time_instant_index + 1]}-t)"
 
             funct_string += ")"
 
@@ -247,7 +255,7 @@ def construct_funct_string_from_variable_data(
                 if time_instant_index != 0:
                     funct_string += "+"
 
-                funct_string += f"({descriptions[time_instant_index]}*heaviside(t-{time_instant})*heaviside({times[time_instant_index+1]}-t))"
+                funct_string += f"({descriptions[time_instant_index]}*heaviside(t-{time_instant})*heaviside({times[time_instant_index + 1]}-t))"
 
             funct_string += ")"
 

src/fourc_webviewer/input_file_utils/io_utils.py

@@ -192,9 +192,7 @@ def get_main_and_clustered_section_names(sections_list):
             # append the main section "FUNCTIONS"
             main_sections.append("FUNCTIONS")
 
-            clustered_sections_to_be_added = (
-                []
-            )  # list of clustered sections to be added
+            clustered_sections_to_be_added = []  # list of clustered sections to be added
             # add current element to clustered sections and remove it from sections
             clustered_sections_to_be_added.append(sections.pop(0))
 
@@ -315,9 +313,7 @@ def get_master_and_linked_material_indices(materials_section):
     master_mat_indices = [mat_item["MAT"] for mat_item in materials_section]
 
     # check whether some of the master materials are actually related to others, and eliminate them from the master material index list
-    linked_mat_indices = (
-        []
-    )  # list of linked material indices for each of the "master" materials: [<list of linked materials for "MASTER" 1>, <list of linked materials for "MASTER" 2>,... ]
+    linked_mat_indices = []  # list of linked material indices for each of the "master" materials: [<list of linked materials for "MASTER" 1>, <list of linked materials for "MASTER" 2>,... ]
     for master_mat_index in master_mat_indices:
         linked_mat_indices.append(
             find_linked_materials(
@@ -339,7 +335,9 @@ def get_master_and_linked_material_indices(materials_section):
                 del master_mat_indices[next_master_list_ind]
                 del linked_mat_indices[next_master_list_ind]
             # next_master_ind stays the same
-            elif set(linked_mat_indices[curr_master_list_ind]).issubset(
+            elif set(
+                linked_mat_indices[curr_master_list_ind]
+            ).issubset(
                 set(linked_mat_indices[next_master_list_ind])
             ):  # this means that the current element is not truly a master -> has to be eliminated and we also break out of the for-loop
                 del master_mat_indices[curr_master_list_ind]
@@ -359,7 +357,8 @@ def get_master_and_linked_material_indices(materials_section):
 def get_variable_data_by_name_in_funct_item(
     funct_section_item: dict, variable_name: str
 ):
-    """Retrieves the entire dictionary for the variable called <variable_name> from the specified function section item, e.g. FUNCT1.
+    """Retrieves the entire dictionary for the variable called <variable_name>
+    from the specified function section item, e.g. FUNCT1.
 
     Args:
         funct_section_item (dict): specified function item as a

src/fourc_webviewer/read_geometry_from_file.py

@@ -4,24 +4,24 @@
 """
 
 import re
+from pathlib import Path
 
 import numpy as np
-from meshio.__about__ import __version__
-from meshio._common import warn
-from meshio._exceptions import ReadError
-from meshio._mesh import Mesh
 import pyvista as pv
 from fourcipp.fourc_input import FourCInput
 from lnmmeshio import read, write
-from lnmmeshio.meshio_to_discretization import mesh2Discretization
 from lnmmeshio.discretization import (
-    PointNodeset,
     LineNodeset,
+    PointNodeset,
     SurfaceNodeset,
     VolumeNodeset,
 )
 from lnmmeshio.fiber import Fiber
-from pathlib import Path
+from lnmmeshio.meshio_to_discretization import mesh2Discretization
+from meshio.__about__ import __version__
+from meshio._common import warn
+from meshio._exceptions import ReadError
+from meshio._mesh import Mesh
 
 exodus_to_meshio_type = {
     "SPHERE": "vertex",
@@ -106,9 +106,8 @@
 
 
 def switch_node_order(mesh_exo: Mesh) -> Mesh:
-    """Switch node orders for read-in Exodus mesh such that the node
-    orders are consistent upon conversion to lnmmeshio's Discretization
-    object.
+    """Switch node orders for read-in Exodus mesh such that the node orders are
+    consistent upon conversion to lnmmeshio's Discretization object.
 
     Args:
         mesh_exo (meshio.Mesh): read-in Exodus mesh
@@ -120,7 +119,6 @@ def switch_node_order(mesh_exo: Mesh) -> Mesh:
 
     # run through elements
     for cell_block in copy_mesh_exo.cells:
-
         cell_type = cell_block.type
 
         if cell_type in special_meshio_mesh_to_lnmmeshio_dis_to_vtu_node_order:
@@ -371,8 +369,10 @@ def read_exodus(filename, use_set_names=False):  # noqa: C901
 
 
 def check_for_geometry_files(fourc_yaml: FourCInput) -> list:
-    """Checks the content of a fourc yaml file for referenced geometry files, e.g., contained within
-    STRUCTURE GEOMETRY / FILE. Returns the identified file, whereby we previously verify whether all given files are the same file (required by 4C currently).
+    """Checks the content of a fourc yaml file for referenced geometry files,
+    e.g., contained within STRUCTURE GEOMETRY / FILE. Returns the identified
+    file, whereby we previously verify whether all given files are the same
+    file (required by 4C currently).
 
     Args:
         fourc_yaml(FourCInput): content of the fourc yaml file to be verified.
@@ -489,14 +489,16 @@ def get_element_ids_of_block(self, element_block_id):
         return np.where(self._dis.cell_data["GROUP_ID"] == element_block_id - 1)[0]
 
     def get_all_nodes_in_element_block_exo(self, element_block_id: int):
-        """Retrieve all unique node indices in a specified element block for Exodus geometry.
+        """Retrieve all unique node indices in a specified element block for
+        Exodus geometry.
+
         Args:
             element_block_id (int): id of the element block to retrieve the nodes from
         Returns:
             ndarray: array of node indices within the specified element block.
         """
 
-        # get cummulative element counts for each block
+        # get cumulative element counts for each block
         cum_el_counts = np.cumsum([len(cs) for cs in self._mesh_exo.cells])
 
         # get all element ids in the considered cell set
@@ -521,6 +523,8 @@ def get_all_nodes_in_element_block_exo(self, element_block_id: int):
         return np.unique(all_node_ids)
 
     def enhance_exo_dis_with_fourc_yaml_info(self):
+        """Enhance contained Discretization with further information from the
+        fourc yaml file -> read in nodesets, and material data."""
 
         # --> read in nodeset info (pointnodesets, linenodesets, surfacenodesets, volumenodesets) based on the design sections specified in the yaml file
         # get all design sections
@@ -675,7 +679,9 @@ def convert_dis_to_vtu(self, override=True):
         )
 
     def prepare_dis_for_vtu_output(self):
-        """Prepares discretization for vtu conversion by adding data contained within the yaml file (e.g. material id, design conditions) as nodal or element data."""
+        """Prepares discretization for vtu conversion by adding data contained
+        within the yaml file (e.g. material id, design conditions) as nodal or
+        element data."""
         self._dis.compute_ids(zero_based=False)
 
         # write node data

src/fourc_webviewer_default_files/tutorial_contact_3d.4C.yaml

@@ -64,7 +64,7 @@ FUNCT1:
     TYPE: "linearinterpolation"
     NUMPOINTS: 3
     TIMES: [0, 0.1, 2]
-    VALUES: [0,0.02,0.02]
+    VALUES: [0, 0.02, 0.02]
 DESIGN SURF MORTAR CONTACT CONDITIONS 3D:
   - E: 1
     ENTITY_TYPE: node_set_id