Merge branch '235-implement-stresses' into 'developer'

Resolve "Implement stresses"

See merge request AM/AMfe!234

amfe/component/mesh_component.py

@@ -36,6 +36,7 @@ def __init__(self, mesh=Mesh()):
         self._neumann = NeumannManager()
         self._assembly = Assembly()
         self._constraints = ConstraintManager()
+        self._elements_on_nodes = np.zeros(mesh.no_of_nodes)
 
     def __str__(self):
         """
@@ -140,6 +141,18 @@ def _assign_material_by_eleids(self, materialobj, eleids, physics):
         self._C_csr = self._assembly.preallocate(self._mapping.no_of_dofs, self._mapping.elements2global)
         self._M_csr = self._C_csr.copy()
         self._f_glob_int = np.zeros(self._C_csr.shape[1])
+        self._update_elements_on_nodes()
+
+    def _update_elements_on_nodes(self):
+        logger = logging.getLogger(__name__)
+        logger.debug('Update number of elements on nodes...')
+        self._elements_on_nodes = np.zeros(self._mesh.no_of_nodes)
+        for nodeidx in range(self._mesh.no_of_nodes):
+            nodeids = self._mesh.get_nodeids_by_nodeidxs([nodeidx])
+            element_ids_all = self._mesh.get_elementids_by_nodeids(nodeids)
+            element_ids = np.intersect1d(self._ele_obj_df['fk_mesh'].values, element_ids_all)
+            self._elements_on_nodes[nodeidx] = element_ids.size
+        logger.debug('Update number of elements on nodes finished.')
 
     # -- ASSIGN NEUMANN CONDITION METHODS -----------------------------------------------------------------
     def assign_neumann(self, name, condition, tag_values, tag='_groups', ignore_nonexistent=False):

amfe/component/structural_component.py

@@ -4,7 +4,8 @@
 #
 
 from scipy.sparse import csc_matrix
-
+import numpy as np
+import logging
 from .mesh_component import MeshComponent
 from amfe.assembly.structural_assembly import StructuralAssembly
 from amfe.component.constants import ELEPROTOTYPEHELPERLIST, SHELLELEPROTOTYPEHELPERLIST
@@ -31,6 +32,8 @@ def __init__(self, mesh=Mesh()):
         self._M_csr = None
         self._f_glob_int = None
         self._f_glob_ext = None
+        self._stresses = None
+        self._strains = None
 
     def g_holo(self, q, t):
         """
@@ -277,3 +280,36 @@ def f_ext(self, q, dq, t):
                                                          neumann_connectivities, neumann_dofs, q, t,
                                                          f_glob=self._f_glob_ext)
         return self._f_glob_ext
+
+    def strains_and_stresses(self, q, dq, t):
+        """
+        Update strain- and stress-fields
+
+        Parameters
+        ----------
+        q : ndarray
+            Displacement field in voigt notation.
+        dq : ndarray
+            Velocity field in voigt notation.
+        t : float
+            time
+
+        Returns
+        -------
+        strains : ndarray
+            nodal strains for each node
+
+        stresses : ndarray
+            nodal stresses for each node
+        """
+        logger = logging.getLogger(__name__)
+        logger.info('assembling strains and stresses...')
+
+        K_csr, f_glob, self._stresses, self._strains = self._assembly.assemble_k_f_S_E(self._mesh.nodes, self.ele_obj,
+                                                      self._mesh.get_iconnectivity_by_elementids(
+                                                          self._ele_obj_df['fk_mesh'].values),
+                                                      self._mapping.get_dofs_by_ids(
+                                                          self._ele_obj_df['fk_mapping'].values), self._elements_on_nodes, q, t,
+                                                      self._C_csr, self._f_glob_int)
+
+        return self._strains, self._stresses

amfe/io/postprocessing/reader/amfe_solution_reader.py

@@ -66,6 +66,34 @@ def parse(self, builder):
             builder.write_field('acceleration', PostProcessDataType.VECTOR, t, acceleration_field, index,
                                 MeshEntityType.NODE)
 
+        if self._amfesolution.strain[0] is not None:
+            strains = np.array(self._amfesolution.strain).T
+            strains_field = self._convert_data_2_normal(strains)
+            index = self._meshcomponent.mesh.nodes_df.index.values
+            builder.write_field('strains_normal', PostProcessDataType.VECTOR, t, strains_field, index,
+                                MeshEntityType.NODE)
+
+        if self._amfesolution.strain[0] is not None:
+            strains = np.array(self._amfesolution.strain).T
+            strains_field = self._convert_data_2_shear(strains)
+            index = self._meshcomponent.mesh.nodes_df.index.values
+            builder.write_field('strains_shear', PostProcessDataType.VECTOR, t, strains_field, index,
+                                MeshEntityType.NODE)
+
+        if self._amfesolution.stress[0] is not None:
+            stresses = np.array(self._amfesolution.stress).T
+            stresses_field = self._convert_data_2_normal(stresses)
+            index = self._meshcomponent.mesh.nodes_df.index.values
+            builder.write_field('stresses_normal', PostProcessDataType.VECTOR, t, stresses_field, index,
+                                MeshEntityType.NODE)
+
+        if self._amfesolution.stress[0] is not None:
+            stresses = np.array(self._amfesolution.stress).T
+            stresses_field = self._convert_data_2_shear(stresses)
+            index = self._meshcomponent.mesh.nodes_df.index.values
+            builder.write_field('stresses_shear', PostProcessDataType.VECTOR, t, stresses_field, index,
+                                MeshEntityType.NODE)
+
     def _convert_data_2_field(self, data_array):
         """
         Converts a global 1-dimensional array that contains vector field data to a 2-dimensional array.
@@ -82,6 +110,59 @@ def _convert_data_2_field(self, data_array):
         data : ndarray
             2-d-array, reordered solution-data, columns contain x, y and z components
         """
+        data, nodeidxs, nodes2dofs = self._preallocate_dataarray()
+
+        data[nodeidxs, 0, :] = data_array[nodes2dofs[:, 0], :]
+        data[nodeidxs, 1, :] = data_array[nodes2dofs[:, 1], :]
+        if nodes2dofs.shape[1] > 2:
+            data[nodeidxs, 2, :] = data_array[nodes2dofs[:, 2], :]
+        else:
+            data[nodeidxs, 2, :] = np.zeros((len(nodeidxs), len(self._amfesolution.t)))
+        return data
+
+    def _convert_data_2_normal(self, data_array):
+        """
+        Selects only the first three entries of the data_array, assuming, these are the normal-parts of a strain-/
+        stress-tensor and reorders a data-array to match the ordering, which is expected by the Postprocessor-writer.
+
+        Parameters
+        ----------
+        data_array : ndarray
+            array of strains/stresses, that shall be read in.
+
+        Returns
+        -------
+        data : ndarray
+            reordered solution-data
+        """
+        data, nodeidxs, _ = self._preallocate_dataarray()
+
+        for ientry in range(3):
+            data[nodeidxs, ientry, :] = data_array[ientry, nodeidxs, :]
+        return data
+
+    def _convert_data_2_shear(self, data_array):
+        """
+        Selects only the fourth to sixth entries of the data_array, assuming, these are the shear-parts of a strain-/
+        stress-tensor and reorders a data-array to match the ordering, which is expected by the Postprocessor-writer.
+
+        Parameters
+        ----------
+        data_array : ndarray
+            array of strains/stresses, that shall be read in.
+
+        Returns
+        -------
+        data : ndarray
+            reordered solution-data
+        """
+        data, nodeidxs, _ = self._preallocate_dataarray()
+
+        for ientry in range(3, 6):
+            data[nodeidxs, ientry-3, :] = data_array[ientry, nodeidxs, :]
+        return data
+
+    def _preallocate_dataarray(self):
         no_of_timesteps = len(self._amfesolution.t)
 
         mapping = self._meshcomponent.mapping
@@ -93,10 +174,5 @@ def _convert_data_2_field(self, data_array):
         # Write data
         nodeids = mapping._nodal2global.index.values
         nodeidxs = self._meshcomponent.mesh.get_nodeidxs_by_nodeids(nodeids)
-        data[nodeidxs, 0, :] = data_array[nodal2global[:, 0], :]
-        data[nodeidxs, 1, :] = data_array[nodal2global[:, 1], :]
-        if nodal2global.shape[1] > 2:
-            data[nodeidxs, 2, :] = data_array[nodal2global[:, 2], :]
-        else:
-            data[nodeidxs, 2, :] = np.zeros((len(nodeidxs), no_of_timesteps))
-        return data
+
+        return data, nodeidxs, nodal2global

amfe/io/postprocessing/writer/hdf5_post_processor_writer.py

@@ -168,6 +168,7 @@ def _write_field(self, fp, name, field_type, t, data, index, mesh_entity_type):
                 data[nodesiloc, :] = data[:, :]
             else:
                 raise NotImplementedError('Field Data Type {} is not supported for converting'.format(field_type.name))
+
             dataset = fp.create_array(resultsroot, name, data)
             dataset.attrs.data_type = field_type.name
             dataset.attrs.mesh_entitiy_type = mesh_entity_type.name

amfe/material.py

@@ -428,7 +428,6 @@ def plane_stress(self, plane_stress):
             self._plane_stress = plane_stress
             self.notify()
 
-
     def S_Sv_and_C(self, E):
         ''' '''
 

amfe/mesh.py

@@ -107,6 +107,13 @@ def __str__(self):
                "No of boundary elements\t{2:,}\nNo of nodes: \t\t\t{3:,}\nAdress in RAM\t\t\t{4}\n" \
                .format(self.dimension, self.no_of_elements, self.no_of_boundary_elements, self.no_of_nodes, id(self))
 
+    @property
+    def element_ids(self):
+        """
+        Returns all element-ids, that are part of the mesh
+        """
+        return self._el_df.index.values
+
     @property
     def el_df(self):
         return self._el_df
@@ -318,7 +325,7 @@ def get_elementidxs_by_elementids(self, elementids):
 
     def get_elementids_by_elementidxs(self, elementidxs):
         """
-        Returns elementids belonging to elements with elementidxs in connectivity array
+        Returns elementids of elements at the given dataframe-indices
 
         Parameters
         ----------
@@ -331,6 +338,38 @@ def get_elementids_by_elementidxs(self, elementidxs):
         """
         return self._el_df.iloc[elementidxs].index.values
 
+    def get_elementids_by_nodeids(self, nodeids):
+        """
+        Returns elementids of elements, that have at least on of the given nodeids in their connectivity
+
+        Parameters
+        ----------
+        nodeids : iterable
+            nodeids, the connectivity has to be searched for
+
+        Returns
+        -------
+        eleids : ndarray
+            elementids, that belong to the node
+        """
+        #eleids = np.array([], dtype=int)
+        #for element_id in self.element_ids:
+#            connectivity = self.get_connectivity_by_elementids(np.array([element_id]))
+#            for nodeid in nodeids:
+#                if nodeid in connectivity[0]:
+#                    eleids = np.append(eleids, np.array([element_id]))
+#                    break
+
+        def df_fun(arr):
+            for nodeid in nodeids:
+                if nodeid in arr:
+                    return True
+            return False
+
+        eleids = self._el_df.index[self._el_df['connectivity'].apply(df_fun)].values
+
+        return eleids
+
     def get_nodeid_by_coordinates(self, x, y, z=None, epsilon=1e-12):
         """