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