Added routine to convert tet4 mesh to tet10 mesh by adding midside nodes

Author: calleman21

exodus_helper/__init__.py

@@ -13,6 +13,7 @@
 __version_info__ = parse(__version__).release
 
 from .core import CONNECTIVITY_SIDES
+from .core import SIDES_CONNECTIVITY
 from .core import Exodus
 from .core import add_variable
 from .core import get_data_exodus
@@ -27,6 +28,7 @@
 
 from .topology import RectangularPrism
 
+from .reconfigure_mesh import convert_tet4_tet10
 from .reconfigure_mesh import create_sets_canonical
 from .reconfigure_mesh import scale_mesh
 

exodus_helper/core.py

@@ -36,7 +36,12 @@
         (1, 2, 4): 1,
         (2, 3, 4): 2,
         (1, 3, 4): 3,
-        (1, 2, 3): 4}}
+        (1, 2, 3): 4},
+    'TETRA10': {
+        (1, 5, 2, 9, 4, 8): 1,
+        (2, 6, 3, 10, 4, 9): 2,
+        (1, 7, 3, 10, 4, 8): 3,
+        (1, 5, 2, 6, 3, 7): 4}}
 
 items_connectivity = CONNECTIVITY_SIDES.items()
 SIDES_CONNECTIVITY = {
@@ -2041,7 +2046,7 @@ def put_elem_connectivity(self, id_blk, connectivity) -> bool:
         """Store the nodal connectivity for an element block.
 
         Args:
-            connectivity (numpy.ndarray(float)): A 2D array of node indices.
+            connectivity (numpy.ndarray(int)): A 1D array of node ids.
 
         Returns:
             `True` if successful, otherwise returns `False`.
@@ -2050,8 +2055,11 @@ def put_elem_connectivity(self, id_blk, connectivity) -> bool:
         num_elems = self.get_num_elems_in_blk(id_blk)
         num_nodes = self.get_num_nodes_per_elem(id_blk)
         connect = np.reshape(connectivity, (num_elems, num_nodes))
-        self.dataset.variables[f'connect{idx}'][:] = connect
-        return True
+        try:
+            self.dataset.variables[f'connect{idx}'][:] = connect
+            return True
+        except KeyError:
+            return False
 
     def put_elem_id_map(self, elem_id_map) -> bool:
         """Store an element ID mapping.

exodus_helper/dir_exodus.py

@@ -238,16 +238,16 @@
     'put_coords',
     'put_elem_attr',
     'put_elem_attr_names',
-    'put_elem_attr_values', 
-    'put_elem_blk_info', 
+    'put_elem_attr_values',
+    'put_elem_blk_info',
     'put_elem_blk_name',
     'put_elem_blk_names',
     'put_elem_connectivity',
     'put_elem_face_conn',
     'put_elem_id_map',
     'put_element_attribute_names',
     'put_element_property_value',
-    'put_element_variable_name', 
+    'put_element_variable_name',
     'put_element_variable_values',
     'put_face_count_per_polyhedra',
     'put_face_node_conn',

exodus_helper/reconfigure_mesh.py

@@ -6,6 +6,7 @@
 
 # --------------------------------------------------------------------------- #
 
+from itertools import permutations
 import os
 import sys
 
@@ -18,52 +19,102 @@
 
 # --------------------------------------------------------------------------- #
 
-def scale_mesh(filename, scale=(1., 1., 1.)):
-    """Scale the dimensions of a mesh by multiplying the coordinate values of
-    the mesh by the scaling weights.
-
-    Args:
-        filename (str): Name of a .g file created in association with a mesh.
-        scale (tuple(float), optional): A list of scaling weights.
-            Defaults to (1., 1., 1.).
-
-    Raises:
-        AssertionError: The filename must correspond to an existing .g file.
-    """
-    assert os.path.isfile(filename)
-    root, ext = os.path.splitext(filename)
-
-    mesh = Exodus(filename)
-    dataset_from = mesh.dataset
-
-    dataset_to = netCDF4.Dataset(
-        root + '_scaled' + ext, mode='w', format=dataset_from.data_model)
-
-    for attr in dataset_from.ncattrs():
-        dataset_to.setncattr(attr, dataset_from.getncattr(attr))
+IDXS_EDGES_4 = [[0, 1], [1, 2], [0, 2], [0, 3], [1, 3], [2, 3]]
 
-    dimensions = dataset_from.dimensions
-    for k, v in dimensions.items():
-        dataset_to.createDimension(k, size=v.size)
 
-    variables = dataset_from.variables
-    for k, v in variables.items():
-        attrs = v.ncattrs()
-        fargs = {}
-        if '_FillValue' in attrs:
-            fargs['fill_value'] = v.getncattr('_FillValue')
-            attrs.remove('_FillValue')
-        variable = dataset_to.createVariable(
-            k, v.datatype, dimensions=v.dimensions, **fargs)
-        variable[:] = v[:]
-        for attr in attrs:
-            variable.setncattr(attr, v.getncattr(attr))
+# --------------------------------------------------------------------------- #
 
-    mesh.close()
-    dataset_to.variables['coordx'][:] *= scale[0]
-    dataset_to.variables['coordy'][:] *= scale[1]
-    dataset_to.variables['coordz'][:] *= scale[2]
-    dataset_to.close()
+def convert_tet4_tet10(filename_from, filename_to=None, **kwargs):
+
+    mesh_from = Exodus(filename_from)
+
+    num_elements = mesh_from.get_num_elems()
+    connectivity_to = np.zeros((num_elements, 10), dtype=int)
+    connectivity_from = mesh_from.get_elem_connectivity_full()[:]
+    connectivity_to[:, :4] = connectivity_from
+
+    coords_from = np.stack(mesh_from.get_coords()).T
+    num_nodes_from = len(coords_from)
+
+    dims = mesh_from.dataset.dimensions
+    num_ns = dims['num_ns_global'].size
+    num_ss = dims['num_ss_global'].size
+
+    num_side_sss = [dims[f'num_side_ss{n}'].size for n in range(1, num_ss + 1)]
+
+    pairs_elem_all = {}
+    coords_to = coords_from.tolist()
+    for idx_elem, ids_node in enumerate(connectivity_from):
+        pairs_elem = [tuple(ids_node[i]) for i in IDXS_EDGES_4]
+        for idx_pair, pair in enumerate(pairs_elem):
+            if pair in pairs_elem_all:
+                id_node_new = pairs_elem_all[pair]
+            else:
+                id_node_new = num_nodes_from + len(pairs_elem_all) + 1
+                pairs_elem_all[pair] = id_node_new
+                idxs_old = [p - 1 for p in pair]
+                coords_to.append(np.mean(coords_from[idxs_old, :], axis=0))
+            connectivity_to[idx_elem, 4 + idx_pair] = id_node_new
+
+    ids_node_set = mesh_from.get_node_set_ids()
+    sets_nodes = []
+    for id_node_set in ids_node_set:
+        ids_node = mesh_from.get_node_set_nodes(id_node_set).tolist()
+        pairs = permutations(ids_node, 2)
+        for pair in pairs:
+            id_new = pairs_elem_all.get(pair, pairs_elem_all.get(pair[::-1]))
+            if id_new not in ids_node and id_new is not None:
+                ids_node.append(id_new)
+        sets_nodes.append(ids_node)
+
+    num_nod_nss = [len(s) for s in sets_nodes]
+
+    ids_blk = mesh_from.get_elem_blk_ids()
+    define_maps = True
+    nums_elems = [mesh_from.get_num_elems_in_blk(i) for i in ids_blk]
+    nums_nodes = [10 for i in ids_blk]
+    nums_attrs = [0 for i in ids_blk]
+    elem_types = ['TETRA10' for i in ids_blk]
+
+    dict_attrs, dict_dimensions, dict_variables = mesh_from._get_dicts_netcdf()
+
+    if filename_to is None:
+        filename_to = '_tet10'.join(os.path.splitext(filename_from))
+    mesh_to = Exodus(
+        filename_to,
+        mode='w',
+        num_info=mesh_from.get_num_info_records() + 1,
+        num_dim=mesh_from.numDim,
+        num_nodes=len(coords_to),
+        num_elem=num_elements,
+        num_el_blk=len(ids_blk),
+        num_node_sets=num_ns,
+        num_side_sets=num_ss,
+        num_nod_nss=num_nod_nss,
+        num_side_sss=num_side_sss,
+        title=kwargs.get('title', mesh_from.get_title()),
+        **kwargs)
+
+    mesh_to.put_coords(*np.array(coords_to).T)
+    mesh_to.put_concat_elem_blk(
+        ids_blk, elem_types, nums_elems, nums_nodes, nums_attrs, define_maps)
+    variables = mesh_to.dataset.variables
+    for i in ids_blk:
+        idxs_elem = mesh_from.get_idxs_elem_in_blk(i)
+        variables[f'connect{i}'][:] = connectivity_to[idxs_elem]
+
+    ids_side_set = mesh_from.get_side_set_ids()
+    for id_side_set in ids_side_set:
+        ids_elem, ids_side = mesh_from.get_side_set(id_side_set)
+        num_sides = mesh_from.get_num_faces_in_side_set(id_side_set)
+        mesh_to.put_side_set_params(id_side_set, num_sides)
+        mesh_to.put_side_set(id_side_set, ids_elem, ids_side)
+
+    for id_node_set, ids_node in zip(ids_node_set, sets_nodes):
+        mesh_to.put_node_set_params(id_node_set, len(ids_node))
+        mesh_to.put_node_set(id_node_set, ids_node)
+
+    return mesh_to
 
 
 def create_sets_canonical(filename):
@@ -208,6 +259,54 @@ def create_sets_canonical(filename):
     dataset_to.close()
 
 
+def scale_mesh(filename, scale=(1., 1., 1.)):
+    """Scale the dimensions of a mesh by multiplying the coordinate values of
+    the mesh by the scaling weights.
+
+    Args:
+        filename (str): Name of a .g file created in association with a mesh.
+        scale (tuple(float), optional): A list of scaling weights.
+            Defaults to (1., 1., 1.).
+
+    Raises:
+        AssertionError: The filename must correspond to an existing .g file.
+    """
+    assert os.path.isfile(filename)
+    root, ext = os.path.splitext(filename)
+
+    mesh = Exodus(filename)
+    dataset_from = mesh.dataset
+
+    dataset_to = netCDF4.Dataset(
+        root + '_scaled' + ext, mode='w', format=dataset_from.data_model)
+
+    for attr in dataset_from.ncattrs():
+        dataset_to.setncattr(attr, dataset_from.getncattr(attr))
+
+    dimensions = dataset_from.dimensions
+    for k, v in dimensions.items():
+        dataset_to.createDimension(k, size=v.size)
+
+    variables = dataset_from.variables
+    for k, v in variables.items():
+        attrs = v.ncattrs()
+        fargs = {}
+        if '_FillValue' in attrs:
+            fargs['fill_value'] = v.getncattr('_FillValue')
+            attrs.remove('_FillValue')
+        variable = dataset_to.createVariable(
+            k, v.datatype, dimensions=v.dimensions, **fargs)
+        variable[:] = v[:]
+        for attr in attrs:
+            variable.setncattr(attr, v.getncattr(attr))
+
+    mesh.close()
+    dataset_to.variables['coordx'][:] *= scale[0]
+    dataset_to.variables['coordy'][:] *= scale[1]
+    dataset_to.variables['coordz'][:] *= scale[2]
+    dataset_to.close()
+
+
 # --------------------------------------------------------------------------- #
 
 if __name__ == '__main__':

tests/test_exodus_helper.py

@@ -19,6 +19,7 @@
     from exodus_helper.dir_exodus import _dir_exodus_full
 import exodus_helper
 from exodus_helper.dir_exodus import _attr_inputs
+from exodus_helper.reconfigure_mesh import IDXS_EDGES_4
 
 
 # --------------------------------------------------------------------------- #
@@ -1820,6 +1821,20 @@ def test_calculate_volumes_block(dir_test_file):
 
 # Testing the reconfigure_mesh module --------------------------------------- #
 
+def test_convert_tet4_tet10(dir_test_file):
+    file_path = os.path.join(dir_test_file, 'test_convert_tet.g')
+    mesh_converted = exodus_helper.convert_tet4_tet10(file_path)
+    connectivity = mesh_converted.get_elem_connectivity_full()[:]
+    coords = np.stack(mesh_converted.get_coords()).T
+    for ids_node in connectivity:
+        for idx, edge in enumerate(IDXS_EDGES_4):
+            coords_0 = coords[ids_node[edge[0]] - 1]
+            coords_1 = coords[ids_node[edge[1]] - 1]
+            coords_2 = coords[ids_node[idx + 4] - 1]
+            assert np.allclose(0.5 * (coords_0 + coords_1), coords_2)
+    os.remove(os.path.join(dir_test_file, 'test_convert_tet_tet10.g'))
+
+
 def test_create_sets_canonical(dir_test_file):
     file_path = os.path.join(dir_test_file, 'test_create_sets_canonical.g')
     exodus_helper.create_sets_canonical(file_path)