Merge remote-tracking branch 'origin/develop'

femio/config.py

@@ -190,8 +190,15 @@
     'vtk': 'vtk',
 }
 
+# https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html
 DICT_VTK_ID_TO_ELEMENT_TYPE = {
+    1: 'pt',
+    3: 'line',
+    5: 'tri',
+    7: 'polygon',
+    9: 'quad',
     10: 'tet',
+    11: 'voxel',
     12: 'hex',
     13: 'prism',
     14: 'pyr',

femio/fem_data.py

@@ -713,7 +713,9 @@ def flatten(indices, element_type):
             nodes=nodes, elements=elements, nodal_data=nodal_data,
             elemental_data={})
 
-    def to_facets(self, remove_duplicates=True):
+    def to_facets(
+            self, remove_duplicates=True, return_dict_facets=False,
+            dict_facets=None):
         """Convert the FEMData object to the facet data including facets inside
         the solid.
 
@@ -722,19 +724,28 @@ def to_facets(self, remove_duplicates=True):
         remove_duplicates: bool, optional
             If True, remove duplicated faces and remain only one. The default
             is True.
+        return_dict_facets: bool, optional
+            If True, also return dict_facets.
 
         Returns
         -------
         FEMData:
             Facets FEMData object.
         """
-        dict_facets = self.extract_facets(remove_duplicates=remove_duplicates)
-        elements = self.elements.to_surface(dict_facets)
+        if dict_facets is None:
+            dict_facets = self.extract_facets(
+                remove_duplicates=remove_duplicates)
 
-        return FEMData(
+        elements = self.elements.to_surface(dict_facets)
+        facet_fem_data = FEMData(
             nodes=self.nodes, elements=elements, nodal_data=self.nodal_data,
             elemental_data={})
 
+        if return_dict_facets:
+            return facet_fem_data, dict_facets
+        else:
+            return facet_fem_data
+
     @staticmethod
     @njit
     def convert_polyhedron(now_ids, new_ids, poly):

femio/fem_elemental_attribute.py

@@ -94,6 +94,7 @@ def from_meshio(cls, cell_data):
             cls._from_meshio(k, v)
             for k, v in cell_data.items()
             if k in [
+                'line', 'quad', 'triangle',
                 'tetra', 'tetra10', 'hexahedron', 'wedge',
                 'pyramid', 'hexa_prism']})
 
@@ -408,13 +409,33 @@ def _generate_surface_ids_tuple(self, surface_ids_dict):
 
     def _generate_surface(self, tuple_surface_ids):
         ret = {}
+        last_id_end = 0
         for surface_ids in tuple_surface_ids:
             if len(surface_ids) == 0:
                 continue
-            ret.update(self._generate_surface_core(surface_ids))
+            item_surface = self._generate_surface_core(
+                surface_ids, last_id_end)
+            key = list(item_surface.keys())[0]
+            new_value = list(item_surface.values())[0]
+            if key in ret:
+                original_value = ret[key]
+                try:
+                    data = np.concatenate(
+                        [original_value.data, new_value.data])
+                except ValueError:
+                    # Manage if the arrays are not in the same shape
+                    data = np.array(
+                        list(original_value.data) + list(new_value.data))
+                ret[key] = FEMAttribute(
+                    key,
+                    ids=np.concatenate([original_value.ids, new_value.ids]),
+                    data=data)
+            else:
+                ret.update(item_surface)
+            last_id_end = new_value.ids[-1]
         return ret
 
-    def _generate_surface_core(self, surface_ids):
+    def _generate_surface_core(self, surface_ids, last_id_end=0):
         shape = surface_ids.shape
         if len(shape) == 1:
             element_type = 'polygon'
@@ -424,13 +445,16 @@ def _generate_surface_core(self, surface_ids):
                 element_type = 'tri'
             elif n_node_per_element == 4:
                 element_type = 'quad'
+            elif n_node_per_element > 4:
+                element_type = 'polygon'
             else:
                 raise NotImplementedError(
                     f"Unsupported shape of elements: {shape}")
         return {
             element_type:
             FEMAttribute(
-                element_type, np.arange(len(surface_ids)) + 1, surface_ids)}
+                element_type,
+                np.arange(len(surface_ids)) + last_id_end + 1, surface_ids)}
 
     def detect_element_type(self, element_data):
         n_node_per_element = element_data.shape[1]

femio/formats/obj/obj.py

@@ -58,12 +58,20 @@ def read_elements(self, string_series):
             for s in str_element_data], dtype=object)
         element_ids = np.arange(len(element_data), dtype=int) + 1
         lengths = np.array([len(e) for e in element_data])
-        unique_lengths = np.unique(lengths)
 
-        elements = {
-            self.DICT_LENGTH2TYPE[length]:
-            self.extract_elements(element_ids, element_data, lengths, length)
-            for length in unique_lengths}
+        keys = list(self.DICT_LENGTH2TYPE.keys())
+        is_polygon = ~np.isin(lengths, keys)
+
+        elements = {}
+        for length in keys:
+            if np.all(lengths != length):
+                continue
+            elements[self.DICT_LENGTH2TYPE[length]] = self.extract_elements(
+                element_ids, element_data, lengths, length)
+        if np.any(is_polygon):
+            elements['polygon'] = FEMAttribute(
+                'polygon', element_ids[is_polygon], element_data[is_polygon]
+            )
         return FEMElementalAttribute('ELEMENT', elements)
 
     def extract_elements(self, element_ids, element_data, lengths, length):

femio/formats/obj/write_obj.py

@@ -1,4 +1,3 @@
-
 import pandas as pd
 
 
@@ -21,22 +20,21 @@ def write(self, file_name=None, *, overwrite=False):
             # Node
             n_node = len(self.fem_data.nodes.ids)
             f.write(pd.DataFrame(
-                index=['v']*n_node, data=self.fem_data.nodes.data
+                index=['v'] * n_node, data=self.fem_data.nodes.data
             ).to_csv(sep=' ', header=False, na_rep='NaN'))
 
             surface_indices, _ \
                 = self.fem_data.extract_surface()
             if isinstance(surface_indices, dict):
                 # Mixed elements
                 for v in surface_indices.values():
-                    n_element = len(v)
-                    f.write(pd.DataFrame(
-                        index=['f']*n_element, data=v+1
-                    ).to_csv(sep=' ', header=False, na_rep='NaN'))
+                    row_strings = [
+                        "f " + " ".join(str(x + 1) for x in row)
+                        for row in v.tolist()]
+                    f.write("\n".join(row_strings) + "\n")
             else:
-                n_element = len(surface_indices)
-                f.write(pd.DataFrame(
-                    index=['f']*n_element, data=surface_indices+1
-                ).to_csv(sep=' ', header=False, na_rep='NaN'))
-
+                row_strings = [
+                    "f " + " ".join(str(x + 1) for x in row)
+                    for row in surface_indices.tolist()]
+                f.write("\n".join(row_strings) + "\n")
         return file_name

femio/functions.py

@@ -71,7 +71,8 @@ def align_nnz(sparses):
         return align_nnz([s.tocsr() for s in sparses])
 
 
-def remove_duplicates(connectivities, return_index=False):
+def remove_duplicates(
+        connectivities, return_index=False, return_inverse=False, end=None):
     """Remove duplicating elements.
 
     Parameters
@@ -80,19 +81,32 @@ def remove_duplicates(connectivities, return_index=False):
         Element connectivities.
     return_index: bool, optional
         If True, return also indices of unique.
+    return_inverse: bool, optional
+        If True, return also the inverse_indice
+    end: int, optional
+        If fed, use only first `end` elements to detect duplication.
 
     Returns
     -------
     connectivities_wo_duplications: numpy.ndarray[int]
         Element connectivities without duplications.
     """
     sorted_connectivities = [
-        np.sort(connectivity) for connectivity in connectivities]
-    _, indices = np.unique(sorted_connectivities, axis=0, return_index=True)
+        np.sort(connectivity)[:end] for connectivity in connectivities]
+    unique = np.unique(
+        sorted_connectivities, axis=0,
+        return_index=True, return_inverse=return_inverse)
+    indices = unique[1]
+    ret = [connectivities[indices]]
+
     if return_index:
-        return connectivities[indices], indices
+        ret.append(indices)
+    if return_inverse:
+        ret.append(unique[-1])
+    if len(ret) == 1:
+        return ret[0]
     else:
-        return connectivities[indices]
+        return tuple(ret)
 
 
 def convert_array2symmetric_matrix(