From 4ee85dbf6a9a598cdb5ba46c1054daf19cf1994a Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Tue, 26 Mar 2024 14:08:39 +0100
Subject: [PATCH 1/6] Added GETELEMTABLES action call to LinsolveModule

---
 myjive/implicit/linsolvemodule.py | 13 +++++++++++++
 myjive/names/globnames.py         |  1 +
 2 files changed, 14 insertions(+)

diff --git a/myjive/implicit/linsolvemodule.py b/myjive/implicit/linsolvemodule.py
index ff8f1f6..e7d413f 100644
--- a/myjive/implicit/linsolvemodule.py
+++ b/myjive/implicit/linsolvemodule.py
@@ -30,6 +30,7 @@ class LinsolveModule(SolverModule):
             self, props, "getStrainMatrix", default=False
         )
         self._tnames = optional_argument(self, props, "tables", default=[])
+        self._etnames = optional_argument(self, props, "elemTables", default=[])
 
         solvertype = solverprops[TYPE]
         self._solver = globdat[gn.SOLVERFACTORY].get_solver(solvertype)
@@ -101,6 +102,18 @@ class LinsolveModule(SolverModule):
             table.to_table()
             globdat[gn.TABLES][name] = table
 
+        if gn.ELEMTABLES not in globdat:
+            globdat[gn.ELEMTABLES] = {}
+
+        for name in self._etnames:
+            elemcount = len(globdat[gn.ESET])
+            table = Table(size=elemcount)
+
+            for model in self.get_relevant_models("GETELEMTABLE", globdat[gn.MODELS]):
+                table = model.GETELEMTABLE(name, table, globdat)
+
+            globdat[gn.ELEMTABLES][name] = table
+
         return "ok"
 
     def get_ext_vector(self, globdat):
diff --git a/myjive/names/globnames.py b/myjive/names/globnames.py
index a312e01..444c123 100644
--- a/myjive/names/globnames.py
+++ b/myjive/names/globnames.py
@@ -26,6 +26,7 @@ class GlobNames:
     MATRIX2 = "matrix2"
     MATRIXB = "matrixB"
     TABLES = "tables"
+    ELEMTABLES = "elemTables"
     LAMBDA = "lambda"
     MODULEDATA = "module"
     SLIDERS = "sliders"
-- 
GitLab


From 81726db3e4902662dc9179752d4686d1c6eedf88 Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Tue, 26 Mar 2024 14:09:13 +0100
Subject: [PATCH 2/6] Added two GETELEMTABLE actions to SolidModel

---
 myjivex/models/solidmodel.py | 105 +++++++++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)

diff --git a/myjivex/models/solidmodel.py b/myjivex/models/solidmodel.py
index 6a4f27d..50d14a8 100644
--- a/myjivex/models/solidmodel.py
+++ b/myjivex/models/solidmodel.py
@@ -39,6 +39,13 @@ class SolidModel(Model):
             table, tbwts = self._get_elem_size(table, tbwts, globdat, **kwargs)
         return table, tbwts
 
+    def GETELEMTABLE(self, name, table, globdat, **kwargs):
+        if "strain" in name:
+            table = self._get_strains_by_elem(table, globdat, **kwargs)
+        if "size" in name:
+            table = self._get_elem_size_by_elem(table, globdat, **kwargs)
+        return table
+
     def configure(self, globdat, **props):
 
         # Get props
@@ -438,6 +445,104 @@ class SolidModel(Model):
 
         return table, tbwts
 
+    def _get_strains_by_elem(self, table, globdat, solution=None):
+        if table is None:
+            nodecount = len(globdat[gn.ESET])
+            table = Table(size=nodecount)
+
+        # Convert the table to an XTable and store the original class
+        xtable = to_xtable(table)
+
+        # Get the STATE0 vector if no custom displacement field is provided
+        if solution is None:
+            disp = globdat[gn.STATE0]
+        else:
+            disp = solution
+
+        # Check that only one integration point is used
+        if self._ipcount != 1:
+            raise ValueError(
+                "Strains per element can only be obtained if a single integration point is used"
+            )
+
+        # Add the columns of all strain components to the table
+        if self._rank == 1:
+            jcols = xtable.add_columns(["xx"])
+        elif self._rank == 2:
+            jcols = xtable.add_columns(["xx", "yy", "xy"])
+        elif self._rank == 3:
+            jcols = xtable.add_columns(["xx", "yy", "zz", "xy", "yz", "zx"])
+
+        for ielem in self._ielems:
+            # Get the nodal coordinates of each element
+            inodes = self._elems.get_elem_nodes(ielem)
+            idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
+            coords = self._nodes.get_some_coords(inodes)
+
+            # Get the shape functions, gradients, weights and coordinates of each integration point
+            sfuncs = self._shape.get_shape_functions()
+            grads, weights = self._shape.get_shape_gradients(coords)
+
+            # Get the nodal displacements
+            eldisp = disp[idofs]
+
+            # Reset the element stress matrix and weights
+            eleps = np.zeros((self._shape.node_count(), self._strcount))
+            elwts = np.zeros(self._shape.node_count())
+
+            ip = 0
+
+            # Get the B matrix for each integration point
+            B = self._get_B_matrix(grads[:, :, ip])
+
+            # Get the strain of the element in the integration point
+            strain = np.matmul(B, eldisp)
+
+            # Add the element stresses to the global stresses
+            xtable.set_row_values(ielem, jcols, strain)
+
+        # Convert the table back to the original class
+        table = xtable.to_table()
+
+        return table
+
+    def _get_elem_size_by_elem(self, table, globdat):
+        if table is None:
+            nodecount = len(globdat[gn.ESET])
+            table = Table(size=nodecount)
+
+        # Convert the table to an XTable and store the original class
+        xtable = to_xtable(table)
+
+        # Add the columns of all stress components to the table
+        jcol = xtable.add_column("")
+
+        for ielem in self._ielems:
+            # Get the nodal coordinates of each element
+            inodes = self._elems.get_elem_nodes(ielem)
+            coords = self._nodes.get_some_coords(inodes)
+
+            # Compute the maximum edge length
+            max_edge = 0.0
+            for i, inode in enumerate(inodes):
+                for j, jnode in enumerate(inodes):
+                    icoords = coords[:, i]
+                    jcoords = coords[:, j]
+                    edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
+                    if edge > max_edge:
+                        max_edge = edge
+
+            # Set the element size and weights
+            elsize = max_edge
+
+            # Add the element stresses to the global stresses
+            xtable.set_value(ielem, jcol, elsize)
+
+        # Convert the table back to the original class
+        table = xtable.to_table()
+
+        return table
+
     def _get_N_matrix(self, sfuncs):
         return get_N_matrix_jit(sfuncs, self._dofcount, self._rank)
 
-- 
GitLab


From 2f6d80a4a0e5e7d44f2f1a89da454f7066a1a8d0 Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Wed, 27 Mar 2024 11:28:21 +0100
Subject: [PATCH 3/6] Major overhaul of gettable actions in SolidModel

---
 myjivex/models/solidmodel.py | 426 ++++++++++++++---------------------
 1 file changed, 175 insertions(+), 251 deletions(-)

diff --git a/myjivex/models/solidmodel.py b/myjivex/models/solidmodel.py
index 50d14a8..b3e108a 100644
--- a/myjivex/models/solidmodel.py
+++ b/myjivex/models/solidmodel.py
@@ -30,18 +30,18 @@ class SolidModel(Model):
 
     def GETTABLE(self, name, table, tbwts, globdat, **kwargs):
         if "stress" in name:
-            table, tbwts = self._get_stresses(table, tbwts, globdat, **kwargs)
+            table, tbwts = self._get_stress_by_node(table, tbwts, globdat, **kwargs)
         elif "strain" in name:
-            table, tbwts = self._get_strains(table, tbwts, globdat, **kwargs)
+            table, tbwts = self._get_strain_by_node(table, tbwts, globdat, **kwargs)
         elif "stiffness" in name:
-            table, tbwts = self._get_stiffness(table, tbwts, globdat, **kwargs)
+            table, tbwts = self._get_stiffness_by_node(table, tbwts, globdat, **kwargs)
         elif "size" in name:
-            table, tbwts = self._get_elem_size(table, tbwts, globdat, **kwargs)
+            table, tbwts = self._get_elem_size_by_node(table, tbwts, globdat, **kwargs)
         return table, tbwts
 
     def GETELEMTABLE(self, name, table, globdat, **kwargs):
         if "strain" in name:
-            table = self._get_strains_by_elem(table, globdat, **kwargs)
+            table = self._get_strain_by_elem(table, globdat, **kwargs)
         if "size" in name:
             table = self._get_elem_size_by_elem(table, globdat, **kwargs)
         return table
@@ -218,18 +218,7 @@ class SolidModel(Model):
 
         return B, wts
 
-    def _get_strains(self, table, tbwts, globdat, solution=None):
-        if table is None:
-            nodecount = len(globdat[gn.NSET])
-            table = Table(size=nodecount)
-
-        if tbwts is None:
-            nodecount = len(globdat[gn.NSET])
-            tbwts = np.zeros(nodecount)
-
-        # Convert the table to an XTable and store the original class
-        xtable = to_xtable(table)
-
+    def _get_strain_by_node(self, table, tbwts, globdat, solution=None):
         # Get the STATE0 vector if no custom displacement field is provided
         if solution is None:
             disp = globdat[gn.STATE0]
@@ -237,311 +226,236 @@ class SolidModel(Model):
             disp = solution
 
         # Add the columns of all strain components to the table
-        if self._rank == 1:
-            jcols = xtable.add_columns(["xx"])
-        elif self._rank == 2:
-            jcols = xtable.add_columns(["xx", "yy", "xy"])
-        elif self._rank == 3:
-            jcols = xtable.add_columns(["xx", "yy", "zz", "xy", "yz", "zx"])
-
-        for ielem in self._ielems:
-            # Get the nodal coordinates of each element
-            inodes = self._elems.get_elem_nodes(ielem)
-            idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
-            coords = self._nodes.get_some_coords(inodes)
-
-            # Get the shape functions, gradients, weights and coordinates of each integration point
-            sfuncs = self._shape.get_shape_functions()
-            grads, weights = self._shape.get_shape_gradients(coords)
-
-            # Get the nodal displacements
-            eldisp = disp[idofs]
-
-            # Reset the element stress matrix and weights
-            eleps = np.zeros((self._shape.node_count(), self._strcount))
-            elwts = np.zeros(self._shape.node_count())
-
-            for ip in range(self._ipcount):
-                # Get the B matrix for each integration point
-                B = self._get_B_matrix(grads[:, :, ip])
-
-                # Get the strain of the element in the integration point
-                strain = np.matmul(B, eldisp)
-
-                # Compute the element strain and weights
-                eleps += np.outer(sfuncs[:, ip], strain)
-                elwts += sfuncs[:, ip].flatten()
-
-            # Add the element weights to the global weights
-            tbwts[inodes] += elwts
-
-            # Add the element stresses to the global stresses
-            xtable.add_block(inodes, jcols, eleps)
+        comps = self._get_gradient_comps()
+        func = self._get_node_strain
+        table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat, disp=disp)
+        return table, tbwts
 
-        # Convert the table back to the original class
-        table = xtable.to_table()
+    def _get_stress_by_node(self, table, tbwts, globdat, solution=None):
+        # Get the STATE0 vector if no custom displacement field is provided
+        if solution is None:
+            disp = globdat[gn.STATE0]
+        else:
+            disp = solution
 
+        # Add the columns of all strain components to the table
+        comps = self._get_gradient_comps()
+        func = self._get_node_stress
+        table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat, disp=disp)
         return table, tbwts
 
-    def _get_stresses(self, table, tbwts, globdat, solution=None):
-        if table is None:
-            nodecount = len(globdat[gn.NSET])
-            table = Table(size=nodecount)
-
-        if tbwts is None:
-            nodecount = len(globdat[gn.NSET])
-            tbwts = np.zeros(nodecount)
+    def _get_stiffness_by_node(self, table, tbwts, globdat):
+        comps = [""]
+        func = self._get_node_stiffness
+        table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat)
+        return table, tbwts
 
-        # Convert the table to an XTable and store the original class
-        xtable = to_xtable(table)
+    def _get_elem_size_by_node(self, table, tbwts, globdat):
+        comps = [""]
+        func = self._get_node_size
+        table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat)
+        return table, tbwts
 
+    def _get_strain_by_elem(self, table, globdat, solution=None):
         # Get the STATE0 vector if no custom displacement field is provided
         if solution is None:
             disp = globdat[gn.STATE0]
         else:
             disp = solution
 
-        # Add the columns of all stress components to the table
-        if self._rank == 1:
-            jcols = xtable.add_columns(["xx"])
-        elif self._rank == 2:
-            jcols = xtable.add_columns(["xx", "yy", "xy"])
-        elif self._rank == 3:
-            jcols = xtable.add_columns(["xx", "yy", "zz", "xy", "yz", "zx"])
-
-        for ielem in self._ielems:
-            # Get the nodal coordinates of each element
-            inodes = self._elems.get_elem_nodes(ielem)
-            idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
-            coords = self._nodes.get_some_coords(inodes)
+        # Check that only one integration point is used
+        if self._ipcount != 1:
+            raise ValueError(
+                "strain per element can only be obtained if a single integration point is used"
+            )
 
-            # Get the shape functions, gradients, weights and coordinates of each integration point
-            sfuncs = self._shape.get_shape_functions()
-            grads, weights = self._shape.get_shape_gradients(coords)
-            ipcoords = self._shape.get_global_integration_points(coords)
+        # Add the columns of all strain components to the table
+        comps = self._get_gradient_comps()
+        func = self._get_elem_strain
 
-            if self._rank == 2:
-                weights *= self._thickness
+        table = self._fill_by_elem(table, comps, func, globdat, disp=disp)
 
-            # Get the nodal displacements
-            eldisp = disp[idofs]
+        return table
 
-            # Reset the element stress matrix and weights
-            elsig = np.zeros((self._shape.node_count(), self._strcount))
-            elwts = np.zeros(self._shape.node_count())
+    def _get_elem_size_by_elem(self, table, globdat, solution=None):
+        # Add the columns of all strain components to the table
+        comps = [""]
+        func = self._get_elem_size
 
-            for ip in range(self._ipcount):
-                # Get the B matrix for each integration point
-                B = self._get_B_matrix(grads[:, :, ip])
+        table = self._fill_by_elem(table, comps, func, globdat)
 
-                # Get the strain and stress of the element in the integration point
-                strain = np.matmul(B, eldisp)
-                stress = self._mat.stress_at_point(strain, ipcoords[:, ip])
+        return table
 
-                # Compute the element stress and weights
-                elsig += np.outer(sfuncs[:, ip], stress)
-                elwts += sfuncs[:, ip].flatten()
+    def _fill_by_elem(self, table, comps, func, globdat, **fargs):
 
-            # Add the element weights to the global weights
-            tbwts[inodes] += elwts
+        xtable = to_xtable(table)
+        jcols = xtable.add_columns(comps)
 
-            # Add the element stresses to the global stresses
-            xtable.add_block(inodes, jcols, elsig)
+        for ielem in self._ielems:
+            strain = func(ielem, globdat, **fargs)
+            xtable.set_row_values(ielem, jcols, strain)
 
-        # Convert the table back to the original class
         table = xtable.to_table()
+        return table
 
-        return table, tbwts
-
-    def _get_stiffness(self, table, tbwts, globdat):
-        if table is None:
-            nodecount = len(globdat[gn.NSET])
-            table = Table(size=nodecount)
-
-        if tbwts is None:
-            nodecount = len(globdat[gn.NSET])
-            tbwts = np.zeros(nodecount)
+    def _fill_by_node(self, table, tbwts, comps, func, globdat, **fargs):
 
-        # Convert the table to an XTable and store the original class
         xtable = to_xtable(table)
-
-        # Add the column of the Young's modulus to the table
-        jcol = xtable.add_column("")
+        jcols = xtable.add_columns(comps)
 
         for ielem in self._ielems:
-            # Get the nodal coordinates of each element
             inodes = self._elems.get_elem_nodes(ielem)
-            coords = self._nodes.get_some_coords(inodes)
-
-            # Get the shape functions, gradients, weights and coordinates of each integration point
-            sfuncs = self._shape.get_shape_functions()
-            ipcoords = self._shape.get_global_integration_points(coords)
-
-            # Reset the element stress matrix and weights
-            elyoung = np.zeros((self._shape.node_count()))
-            elwts = np.zeros(self._shape.node_count())
-
-            for ip in range(self._ipcount):
-                # Get the stiffness in the integration point
-                E = self._mat._get_E(ipcoords[:, ip])
-
-                # Compute the element stiffness and weights
-                elyoung += E * sfuncs[:, ip]
-                elwts += sfuncs[:, ip].flatten()
+            field, weights = func(ielem, inodes, globdat, **fargs)
 
             # Add the element weights to the global weights
-            tbwts[inodes] += elwts
+            tbwts[inodes] += weights
 
-            # Add the element stiffness to the global stiffness
-            xtable.add_col_values(inodes, jcol, elyoung)
+            # Add the element stress to the global stress
+            xtable.add_block(inodes, jcols, field)
 
         # Convert the table back to the original class
         table = xtable.to_table()
-
         return table, tbwts
 
-    def _get_elem_size(self, table, tbwts, globdat):
-        if table is None:
-            nodecount = len(globdat[gn.NSET])
-            table = Table(size=nodecount)
+    def _get_elem_strain(self, ielem, globdat, disp=None):
 
-        if tbwts is None:
-            nodecount = len(globdat[gn.NSET])
-            tbwts = np.zeros(nodecount)
+        # Get the nodal coordinates of each element
+        inodes = self._elems.get_elem_nodes(ielem)
+        idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
+        coords = self._nodes.get_some_coords(inodes)
 
-        # Convert the table to an XTable and store the original class
-        xtable = to_xtable(table)
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        grads, _ = self._shape.get_shape_gradients(coords)
 
-        # Add the columns of all stress components to the table
-        jcol = xtable.add_column("")
+        # Get the nodal displacements
+        eldisp = disp[idofs]
 
-        for ielem in self._ielems:
-            # Get the nodal coordinates of each element
-            inodes = self._elems.get_elem_nodes(ielem)
-            coords = self._nodes.get_some_coords(inodes)
+        # Get the strain in te first integration point
+        strain = self._get_ip_strain(0, grads, eldisp)
 
-            # Compute the maximum edge length\
-            max_edge = 0.0
-            for i, inode in enumerate(inodes):
-                for j, jnode in enumerate(inodes):
-                    icoords = coords[:, i]
-                    jcoords = coords[:, j]
-                    edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
-                    if edge > max_edge:
-                        max_edge = edge
+        return strain
 
-            # Reset the element size and weights
-            elsize = max_edge * np.ones(self._shape.node_count())
-            elwts = np.ones(self._shape.node_count())
+    def _get_elem_size(self, ielem, globdat):
 
-            # Add the element weights to the global weights
-            tbwts[inodes] += elwts
+        # Get the nodal coordinates of each element
+        inodes = self._elems.get_elem_nodes(ielem)
+        coords = self._nodes.get_some_coords(inodes)
 
-            # Add the element stresses to the global stresses
-            xtable.add_col_values(inodes, jcol, elsize)
+        # Compute the maximum edge length
+        max_edge = 0.0
+        for i, inode in enumerate(inodes):
+            for j, jnode in enumerate(inodes):
+                icoords = coords[:, i]
+                jcoords = coords[:, j]
+                edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
+                if edge > max_edge:
+                    max_edge = edge
 
-        # Convert the table back to the original class
-        table = xtable.to_table()
+        # Set the element size and weights
+        return max_edge
 
-        return table, tbwts
+    def _get_node_strain(self, ielem, inodes, globdat, disp=None):
+        idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
+        coords = self._nodes.get_some_coords(inodes)
 
-    def _get_strains_by_elem(self, table, globdat, solution=None):
-        if table is None:
-            nodecount = len(globdat[gn.ESET])
-            table = Table(size=nodecount)
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        sfuncs = self._shape.get_shape_functions()
+        grads, weights = self._shape.get_shape_gradients(coords)
 
-        # Convert the table to an XTable and store the original class
-        xtable = to_xtable(table)
+        # Get the nodal displacements
+        eldisp = disp[idofs]
 
-        # Get the STATE0 vector if no custom displacement field is provided
-        if solution is None:
-            disp = globdat[gn.STATE0]
-        else:
-            disp = solution
+        # Reset the element stress matrix and weights
+        eleps = np.zeros((self._shape.node_count(), self._strcount))
+        elwts = np.zeros(self._shape.node_count())
 
-        # Check that only one integration point is used
-        if self._ipcount != 1:
-            raise ValueError(
-                "Strains per element can only be obtained if a single integration point is used"
-            )
+        for ip in range(self._ipcount):
+            # Get the strain in each integration point
+            strain = self._get_ip_strain(ip, grads, eldisp)
 
-        # Add the columns of all strain components to the table
-        if self._rank == 1:
-            jcols = xtable.add_columns(["xx"])
-        elif self._rank == 2:
-            jcols = xtable.add_columns(["xx", "yy", "xy"])
-        elif self._rank == 3:
-            jcols = xtable.add_columns(["xx", "yy", "zz", "xy", "yz", "zx"])
+            # Compute the element strain and weights
+            eleps += np.outer(sfuncs[:, ip], strain)
+            elwts += sfuncs[:, ip].flatten()
 
-        for ielem in self._ielems:
-            # Get the nodal coordinates of each element
-            inodes = self._elems.get_elem_nodes(ielem)
-            idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
-            coords = self._nodes.get_some_coords(inodes)
+        return eleps, elwts
 
-            # Get the shape functions, gradients, weights and coordinates of each integration point
-            sfuncs = self._shape.get_shape_functions()
-            grads, weights = self._shape.get_shape_gradients(coords)
+    def _get_node_stress(self, ielem, inodes, globdat, disp=None):
+        idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
+        coords = self._nodes.get_some_coords(inodes)
 
-            # Get the nodal displacements
-            eldisp = disp[idofs]
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        sfuncs = self._shape.get_shape_functions()
+        grads, weights = self._shape.get_shape_gradients(coords)
+        ipcoords = self._shape.get_global_integration_points(coords)
 
-            # Reset the element stress matrix and weights
-            eleps = np.zeros((self._shape.node_count(), self._strcount))
-            elwts = np.zeros(self._shape.node_count())
+        if self._rank == 2:
+            weights *= self._thickness
 
-            ip = 0
+        # Get the nodal displacements
+        eldisp = disp[idofs]
 
-            # Get the B matrix for each integration point
-            B = self._get_B_matrix(grads[:, :, ip])
+        # Reset the element stress matrix and weights
+        elsig = np.zeros((self._shape.node_count(), self._strcount))
+        elwts = np.zeros(self._shape.node_count())
 
-            # Get the strain of the element in the integration point
-            strain = np.matmul(B, eldisp)
+        for ip in range(self._ipcount):
+            # Get the stress in each integration point
+            stress = self._get_ip_stress(ip, ipcoords, grads, eldisp)
 
-            # Add the element stresses to the global stresses
-            xtable.set_row_values(ielem, jcols, strain)
+            # Compute the element stress and weights
+            elsig += np.outer(sfuncs[:, ip], stress)
+            elwts += sfuncs[:, ip].flatten()
 
-        # Convert the table back to the original class
-        table = xtable.to_table()
+        return elsig, elwts
 
-        return table
+    def _get_node_stiffness(self, ielem, inodes, globdat):
+        coords = self._nodes.get_some_coords(inodes)
 
-    def _get_elem_size_by_elem(self, table, globdat):
-        if table is None:
-            nodecount = len(globdat[gn.ESET])
-            table = Table(size=nodecount)
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        sfuncs = self._shape.get_shape_functions()
+        ipcoords = self._shape.get_global_integration_points(coords)
 
-        # Convert the table to an XTable and store the original class
-        xtable = to_xtable(table)
+        # Reset the element stress matrix and weights
+        elyoung = np.zeros((self._shape.node_count(), 1))
+        elwts = np.zeros(self._shape.node_count())
 
-        # Add the columns of all stress components to the table
-        jcol = xtable.add_column("")
+        for ip in range(self._ipcount):
+            # Get the stiffness in the integration point
+            E = self._mat._get_E(ipcoords[:, ip])
 
-        for ielem in self._ielems:
-            # Get the nodal coordinates of each element
-            inodes = self._elems.get_elem_nodes(ielem)
-            coords = self._nodes.get_some_coords(inodes)
+            # Compute the element stiffness and weights
+            elyoung[:,0] += E * sfuncs[:, ip]
+            elwts += sfuncs[:, ip].flatten()
 
-            # Compute the maximum edge length
-            max_edge = 0.0
-            for i, inode in enumerate(inodes):
-                for j, jnode in enumerate(inodes):
-                    icoords = coords[:, i]
-                    jcoords = coords[:, j]
-                    edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
-                    if edge > max_edge:
-                        max_edge = edge
+        return elyoung, elwts
 
-            # Set the element size and weights
-            elsize = max_edge
+    def _get_node_size(self, ielem, inodes, globdat):
+        coords = self._nodes.get_some_coords(inodes)
 
-            # Add the element stresses to the global stresses
-            xtable.set_value(ielem, jcol, elsize)
+        # Compute the maximum edge length\
+        max_edge = 0.0
+        for i, inode in enumerate(inodes):
+            for j, jnode in enumerate(inodes):
+                icoords = coords[:, i]
+                jcoords = coords[:, j]
+                edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
+                if edge > max_edge:
+                    max_edge = edge
 
-        # Convert the table back to the original class
-        table = xtable.to_table()
+        # Reset the element size and weights
+        elsize = max_edge * np.ones((self._shape.node_count(), 1))
+        elwts = np.ones(self._shape.node_count())
 
-        return table
+        return elsize, elwts
+
+    def _get_gradient_comps(self):
+        if self._rank == 1:
+            comps = ["xx"]
+        elif self._rank == 2:
+            comps = ["xx", "yy", "xy"]
+        elif self._rank == 3:
+            comps = ["xx", "yy", "zz", "xy", "yz", "zx"]
+        return comps
 
     def _get_N_matrix(self, sfuncs):
         return get_N_matrix_jit(sfuncs, self._dofcount, self._rank)
@@ -550,3 +464,13 @@ class SolidModel(Model):
         return get_B_matrix_jit(
             grads, self._strcount, self._dofcount, self._shape.node_count(), self._rank
         )
+
+    def _get_ip_strain(self, ip, grads, eldisp):
+        B = self._get_B_matrix(grads[:, :, ip])
+        strain = np.matmul(B, eldisp)
+        return strain
+
+    def _get_ip_stress(self, ip, ipcoords, grads, eldisp):
+        strain = self._get_ip_strain(ip, grads, eldisp)
+        stress = self._mat.stress_at_point(strain, ipcoords[:, ip])
+        return stress
-- 
GitLab


From c74f974974454d1249174c137f7267d715bb6998 Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Wed, 27 Mar 2024 11:42:30 +0100
Subject: [PATCH 4/6] Added get_stress_by_elem and get_stiffness_by_elem
 functions to SolidModel

---
 myjivex/models/solidmodel.py | 76 ++++++++++++++++++++++++++++++++++--
 1 file changed, 72 insertions(+), 4 deletions(-)

diff --git a/myjivex/models/solidmodel.py b/myjivex/models/solidmodel.py
index b3e108a..946bdb5 100644
--- a/myjivex/models/solidmodel.py
+++ b/myjivex/models/solidmodel.py
@@ -29,10 +29,10 @@ class SolidModel(Model):
         return B, wts
 
     def GETTABLE(self, name, table, tbwts, globdat, **kwargs):
-        if "stress" in name:
-            table, tbwts = self._get_stress_by_node(table, tbwts, globdat, **kwargs)
-        elif "strain" in name:
+        if "strain" in name:
             table, tbwts = self._get_strain_by_node(table, tbwts, globdat, **kwargs)
+        elif "stress" in name:
+            table, tbwts = self._get_stress_by_node(table, tbwts, globdat, **kwargs)
         elif "stiffness" in name:
             table, tbwts = self._get_stiffness_by_node(table, tbwts, globdat, **kwargs)
         elif "size" in name:
@@ -42,7 +42,11 @@ class SolidModel(Model):
     def GETELEMTABLE(self, name, table, globdat, **kwargs):
         if "strain" in name:
             table = self._get_strain_by_elem(table, globdat, **kwargs)
-        if "size" in name:
+        elif "stress" in name:
+            table = self._get_stress_by_elem(table, globdat, **kwargs)
+        elif "stiffness" in name:
+            table = self._get_stiffness_by_elem(table, globdat, **kwargs)
+        elif "size" in name:
             table = self._get_elem_size_by_elem(table, globdat, **kwargs)
         return table
 
@@ -277,6 +281,36 @@ class SolidModel(Model):
 
         return table
 
+    def _get_stress_by_elem(self, table, globdat, solution=None):
+        # Get the STATE0 vector if no custom displacement field is provided
+        if solution is None:
+            disp = globdat[gn.STATE0]
+        else:
+            disp = solution
+
+        # Check that only one integration point is used
+        if self._ipcount != 1:
+            raise ValueError(
+                "strain per element can only be obtained if a single integration point is used"
+            )
+
+        # Add the columns of all strain components to the table
+        comps = self._get_gradient_comps()
+        func = self._get_elem_stress
+
+        table = self._fill_by_elem(table, comps, func, globdat, disp=disp)
+
+        return table
+
+    def _get_stiffness_by_elem(self, table, globdat, solution=None):
+        # Add the columns of all strain components to the table
+        comps = [""]
+        func = self._get_elem_stiffness
+
+        table = self._fill_by_elem(table, comps, func, globdat)
+
+        return table
+
     def _get_elem_size_by_elem(self, table, globdat, solution=None):
         # Add the columns of all strain components to the table
         comps = [""]
@@ -335,6 +369,40 @@ class SolidModel(Model):
 
         return strain
 
+    def _get_elem_stress(self, ielem, globdat, disp=None):
+
+        # Get the nodal coordinates of each element
+        inodes = self._elems.get_elem_nodes(ielem)
+        idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
+        coords = self._nodes.get_some_coords(inodes)
+
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        grads, _ = self._shape.get_shape_gradients(coords)
+        ipcoords = self._shape.get_global_integration_points(coords)
+
+        # Get the nodal displacements
+        eldisp = disp[idofs]
+
+        # Get the strain in te first integration point
+        stress = self._get_ip_stress(0, ipcoords, grads, eldisp)
+
+        return stress
+
+    def _get_elem_stiffness(self, ielem, globdat):
+
+        # Get the nodal coordinates of each element
+        inodes = self._elems.get_elem_nodes(ielem)
+        coords = self._nodes.get_some_coords(inodes)
+
+        # Get the shape functions, gradients, weights and coordinates of each integration point
+        grads, _ = self._shape.get_shape_gradients(coords)
+        ipcoords = self._shape.get_global_integration_points(coords)
+
+        # Get the stiffness in the first integration point
+        stiffness = self._mat._get_E(ipcoords[:, 0])
+
+        return stiffness
+
     def _get_elem_size(self, ielem, globdat):
 
         # Get the nodal coordinates of each element
-- 
GitLab


From 31c0331689931cd76b5de39cd2a1aa5c3e1f9d88 Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Wed, 27 Mar 2024 11:43:05 +0100
Subject: [PATCH 5/6] Made get_node_size depend on get_elem_size in solidmodel

---
 myjivex/models/solidmodel.py | 16 ++--------------
 1 file changed, 2 insertions(+), 14 deletions(-)

diff --git a/myjivex/models/solidmodel.py b/myjivex/models/solidmodel.py
index 946bdb5..8c93321 100644
--- a/myjivex/models/solidmodel.py
+++ b/myjivex/models/solidmodel.py
@@ -492,25 +492,13 @@ class SolidModel(Model):
             E = self._mat._get_E(ipcoords[:, ip])
 
             # Compute the element stiffness and weights
-            elyoung[:,0] += E * sfuncs[:, ip]
+            elyoung[:, 0] += E * sfuncs[:, ip]
             elwts += sfuncs[:, ip].flatten()
 
         return elyoung, elwts
 
     def _get_node_size(self, ielem, inodes, globdat):
-        coords = self._nodes.get_some_coords(inodes)
-
-        # Compute the maximum edge length\
-        max_edge = 0.0
-        for i, inode in enumerate(inodes):
-            for j, jnode in enumerate(inodes):
-                icoords = coords[:, i]
-                jcoords = coords[:, j]
-                edge = np.sqrt(np.sum((icoords - jcoords) ** 2))
-                if edge > max_edge:
-                    max_edge = edge
-
-        # Reset the element size and weights
+        max_edge = self._get_elem_size(ielem, globdat)
         elsize = max_edge * np.ones((self._shape.node_count(), 1))
         elwts = np.ones(self._shape.node_count())
 
-- 
GitLab


From c2a7f2e4face768a67f91a71535f731683debecb Mon Sep 17 00:00:00 2001
From: Anne Poot <a.poot-1@tudelft.nl>
Date: Wed, 27 Mar 2024 11:51:27 +0100
Subject: [PATCH 6/6] Compacted code in SolidModel

---
 myjivex/models/solidmodel.py | 110 ++++-------------------------------
 1 file changed, 10 insertions(+), 100 deletions(-)

diff --git a/myjivex/models/solidmodel.py b/myjivex/models/solidmodel.py
index 8c93321..0d34c14 100644
--- a/myjivex/models/solidmodel.py
+++ b/myjivex/models/solidmodel.py
@@ -4,7 +4,7 @@ from myjive.names import GlobNames as gn
 from myjive.model.model import Model
 from ..materials import new_material
 import myjive.util.proputils as pu
-from myjive.util import Table, to_xtable
+from myjive.util import to_xtable
 from myjive.util.proputils import mandatory_argument, mandatory_dict, optional_argument
 from .jit.solidmodel import get_N_matrix_jit, get_B_matrix_jit
 
@@ -223,26 +223,16 @@ class SolidModel(Model):
         return B, wts
 
     def _get_strain_by_node(self, table, tbwts, globdat, solution=None):
-        # Get the STATE0 vector if no custom displacement field is provided
-        if solution is None:
-            disp = globdat[gn.STATE0]
-        else:
-            disp = solution
+        disp = globdat[gn.STATE0] if solution is None else solution
 
-        # Add the columns of all strain components to the table
         comps = self._get_gradient_comps()
         func = self._get_node_strain
         table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat, disp=disp)
         return table, tbwts
 
     def _get_stress_by_node(self, table, tbwts, globdat, solution=None):
-        # Get the STATE0 vector if no custom displacement field is provided
-        if solution is None:
-            disp = globdat[gn.STATE0]
-        else:
-            disp = solution
+        disp = globdat[gn.STATE0] if solution is None else solution
 
-        # Add the columns of all strain components to the table
         comps = self._get_gradient_comps()
         func = self._get_node_stress
         table, tbwts = self._fill_by_node(table, tbwts, comps, func, globdat, disp=disp)
@@ -261,67 +251,41 @@ class SolidModel(Model):
         return table, tbwts
 
     def _get_strain_by_elem(self, table, globdat, solution=None):
-        # Get the STATE0 vector if no custom displacement field is provided
-        if solution is None:
-            disp = globdat[gn.STATE0]
-        else:
-            disp = solution
-
-        # Check that only one integration point is used
+        disp = globdat[gn.STATE0] if solution is None else solution
         if self._ipcount != 1:
-            raise ValueError(
-                "strain per element can only be obtained if a single integration point is used"
-            )
+            raise ValueError("element strain only works with a single Gauss point")
 
-        # Add the columns of all strain components to the table
         comps = self._get_gradient_comps()
         func = self._get_elem_strain
-
         table = self._fill_by_elem(table, comps, func, globdat, disp=disp)
-
         return table
 
     def _get_stress_by_elem(self, table, globdat, solution=None):
-        # Get the STATE0 vector if no custom displacement field is provided
-        if solution is None:
-            disp = globdat[gn.STATE0]
-        else:
-            disp = solution
-
-        # Check that only one integration point is used
+        disp = globdat[gn.STATE0] if solution is None else solution
         if self._ipcount != 1:
-            raise ValueError(
-                "strain per element can only be obtained if a single integration point is used"
-            )
+            raise ValueError("element stress only works with a single Gauss point")
 
-        # Add the columns of all strain components to the table
         comps = self._get_gradient_comps()
         func = self._get_elem_stress
-
         table = self._fill_by_elem(table, comps, func, globdat, disp=disp)
-
         return table
 
     def _get_stiffness_by_elem(self, table, globdat, solution=None):
-        # Add the columns of all strain components to the table
+        if self._ipcount != 1:
+            raise ValueError("element stiffness only works with a single Gauss point")
+
         comps = [""]
         func = self._get_elem_stiffness
-
         table = self._fill_by_elem(table, comps, func, globdat)
-
         return table
 
     def _get_elem_size_by_elem(self, table, globdat, solution=None):
-        # Add the columns of all strain components to the table
         comps = [""]
         func = self._get_elem_size
-
         table = self._fill_by_elem(table, comps, func, globdat)
-
         return table
 
     def _fill_by_elem(self, table, comps, func, globdat, **fargs):
-
         xtable = to_xtable(table)
         jcols = xtable.add_columns(comps)
 
@@ -333,7 +297,6 @@ class SolidModel(Model):
         return table
 
     def _fill_by_node(self, table, tbwts, comps, func, globdat, **fargs):
-
         xtable = to_xtable(table)
         jcols = xtable.add_columns(comps)
 
@@ -341,75 +304,46 @@ class SolidModel(Model):
             inodes = self._elems.get_elem_nodes(ielem)
             field, weights = func(ielem, inodes, globdat, **fargs)
 
-            # Add the element weights to the global weights
             tbwts[inodes] += weights
-
-            # Add the element stress to the global stress
             xtable.add_block(inodes, jcols, field)
 
-        # Convert the table back to the original class
         table = xtable.to_table()
         return table, tbwts
 
     def _get_elem_strain(self, ielem, globdat, disp=None):
-
-        # Get the nodal coordinates of each element
         inodes = self._elems.get_elem_nodes(ielem)
         idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         grads, _ = self._shape.get_shape_gradients(coords)
 
-        # Get the nodal displacements
         eldisp = disp[idofs]
-
-        # Get the strain in te first integration point
         strain = self._get_ip_strain(0, grads, eldisp)
-
         return strain
 
     def _get_elem_stress(self, ielem, globdat, disp=None):
-
-        # Get the nodal coordinates of each element
         inodes = self._elems.get_elem_nodes(ielem)
         idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         grads, _ = self._shape.get_shape_gradients(coords)
         ipcoords = self._shape.get_global_integration_points(coords)
 
-        # Get the nodal displacements
         eldisp = disp[idofs]
-
-        # Get the strain in te first integration point
         stress = self._get_ip_stress(0, ipcoords, grads, eldisp)
-
         return stress
 
     def _get_elem_stiffness(self, ielem, globdat):
-
-        # Get the nodal coordinates of each element
         inodes = self._elems.get_elem_nodes(ielem)
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         grads, _ = self._shape.get_shape_gradients(coords)
         ipcoords = self._shape.get_global_integration_points(coords)
 
-        # Get the stiffness in the first integration point
         stiffness = self._mat._get_E(ipcoords[:, 0])
-
         return stiffness
 
     def _get_elem_size(self, ielem, globdat):
-
-        # Get the nodal coordinates of each element
         inodes = self._elems.get_elem_nodes(ielem)
         coords = self._nodes.get_some_coords(inodes)
 
-        # Compute the maximum edge length
         max_edge = 0.0
         for i, inode in enumerate(inodes):
             for j, jnode in enumerate(inodes):
@@ -419,29 +353,20 @@ class SolidModel(Model):
                 if edge > max_edge:
                     max_edge = edge
 
-        # Set the element size and weights
         return max_edge
 
     def _get_node_strain(self, ielem, inodes, globdat, disp=None):
         idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         sfuncs = self._shape.get_shape_functions()
         grads, weights = self._shape.get_shape_gradients(coords)
 
-        # Get the nodal displacements
         eldisp = disp[idofs]
-
-        # Reset the element stress matrix and weights
         eleps = np.zeros((self._shape.node_count(), self._strcount))
         elwts = np.zeros(self._shape.node_count())
 
         for ip in range(self._ipcount):
-            # Get the strain in each integration point
             strain = self._get_ip_strain(ip, grads, eldisp)
-
-            # Compute the element strain and weights
             eleps += np.outer(sfuncs[:, ip], strain)
             elwts += sfuncs[:, ip].flatten()
 
@@ -450,8 +375,6 @@ class SolidModel(Model):
     def _get_node_stress(self, ielem, inodes, globdat, disp=None):
         idofs = globdat[gn.DOFSPACE].get_dofs(inodes, DOFTYPES[0 : self._rank])
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         sfuncs = self._shape.get_shape_functions()
         grads, weights = self._shape.get_shape_gradients(coords)
         ipcoords = self._shape.get_global_integration_points(coords)
@@ -459,18 +382,12 @@ class SolidModel(Model):
         if self._rank == 2:
             weights *= self._thickness
 
-        # Get the nodal displacements
         eldisp = disp[idofs]
-
-        # Reset the element stress matrix and weights
         elsig = np.zeros((self._shape.node_count(), self._strcount))
         elwts = np.zeros(self._shape.node_count())
 
         for ip in range(self._ipcount):
-            # Get the stress in each integration point
             stress = self._get_ip_stress(ip, ipcoords, grads, eldisp)
-
-            # Compute the element stress and weights
             elsig += np.outer(sfuncs[:, ip], stress)
             elwts += sfuncs[:, ip].flatten()
 
@@ -478,20 +395,14 @@ class SolidModel(Model):
 
     def _get_node_stiffness(self, ielem, inodes, globdat):
         coords = self._nodes.get_some_coords(inodes)
-
-        # Get the shape functions, gradients, weights and coordinates of each integration point
         sfuncs = self._shape.get_shape_functions()
         ipcoords = self._shape.get_global_integration_points(coords)
 
-        # Reset the element stress matrix and weights
         elyoung = np.zeros((self._shape.node_count(), 1))
         elwts = np.zeros(self._shape.node_count())
 
         for ip in range(self._ipcount):
-            # Get the stiffness in the integration point
             E = self._mat._get_E(ipcoords[:, ip])
-
-            # Compute the element stiffness and weights
             elyoung[:, 0] += E * sfuncs[:, ip]
             elwts += sfuncs[:, ip].flatten()
 
@@ -501,7 +412,6 @@ class SolidModel(Model):
         max_edge = self._get_elem_size(ielem, globdat)
         elsize = max_edge * np.ones((self._shape.node_count(), 1))
         elwts = np.ones(self._shape.node_count())
-
         return elsize, elwts
 
     def _get_gradient_comps(self):
-- 
GitLab