Coverage for src/beamme/four_c/solid_shell_thickness_direction.py: 86%
130 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-28 15:20 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2026 BeamMe Authors
4#
5# Permission is hereby granted, free of charge, to any person obtaining a copy
6# of this software and associated documentation files (the "Software"), to deal
7# in the Software without restriction, including without limitation the rights
8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9# copies of the Software, and to permit persons to whom the Software is
10# furnished to do so, subject to the following conditions:
11#
12# The above copyright notice and this permission notice shall be included in
13# all copies or substantial portions of the Software.
14#
15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21# THE SOFTWARE.
22"""This function allows to reorder the connectivity of solid shell elements such that
23the solid shell direction is correctly represented."""
25import numpy as _np
26import pyvista as _pv
28from beamme.core.conf import bme as _bme
29from beamme.core.element import Element as _Element
30from beamme.core.element_volume import VolumeElement as _VolumeElement
31from beamme.core.element_volume import VolumeHEX8 as _VolumeHEX8
32from beamme.core.mesh import Mesh as _Mesh
33from beamme.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates
36def shape_functions_hex8(xi1, xi2, xi3):
37 """Return the shape functions for a hex8 element."""
38 shape_functions = _np.zeros(8)
39 one_over_eight = 0.125
40 shape_functions[0] = one_over_eight * (1 - xi1) * (1 - xi2) * (1 - xi3)
41 shape_functions[1] = one_over_eight * (1 + xi1) * (1 - xi2) * (1 - xi3)
42 shape_functions[2] = one_over_eight * (1 + xi1) * (1 + xi2) * (1 - xi3)
43 shape_functions[3] = one_over_eight * (1 - xi1) * (1 + xi2) * (1 - xi3)
44 shape_functions[4] = one_over_eight * (1 - xi1) * (1 - xi2) * (1 + xi3)
45 shape_functions[5] = one_over_eight * (1 + xi1) * (1 - xi2) * (1 + xi3)
46 shape_functions[6] = one_over_eight * (1 + xi1) * (1 + xi2) * (1 + xi3)
47 shape_functions[7] = one_over_eight * (1 - xi1) * (1 + xi2) * (1 + xi3)
48 return shape_functions
51def shape_functions_derivative_hex8(xi1, xi2, xi3):
52 """Return the derivative of the shape functions for a hex8 element."""
53 derivatives = _np.zeros((3, 8))
54 one_over_eight = 0.125
56 # Derivatives with respect to xi1
57 derivatives[0, 0] = -one_over_eight * (1 - xi2) * (1 - xi3)
58 derivatives[0, 1] = one_over_eight * (1 - xi2) * (1 - xi3)
59 derivatives[0, 2] = one_over_eight * (1 + xi2) * (1 - xi3)
60 derivatives[0, 3] = -one_over_eight * (1 + xi2) * (1 - xi3)
61 derivatives[0, 4] = -one_over_eight * (1 - xi2) * (1 + xi3)
62 derivatives[0, 5] = one_over_eight * (1 - xi2) * (1 + xi3)
63 derivatives[0, 6] = one_over_eight * (1 + xi2) * (1 + xi3)
64 derivatives[0, 7] = -one_over_eight * (1 + xi2) * (1 + xi3)
66 # Derivatives with respect to xi2
67 derivatives[1, 0] = -one_over_eight * (1 - xi1) * (1 - xi3)
68 derivatives[1, 1] = -one_over_eight * (1 + xi1) * (1 - xi3)
69 derivatives[1, 2] = one_over_eight * (1 + xi1) * (1 - xi3)
70 derivatives[1, 3] = one_over_eight * (1 - xi1) * (1 - xi3)
71 derivatives[1, 4] = -one_over_eight * (1 - xi1) * (1 + xi3)
72 derivatives[1, 5] = -one_over_eight * (1 + xi1) * (1 + xi3)
73 derivatives[1, 6] = one_over_eight * (1 + xi1) * (1 + xi3)
74 derivatives[1, 7] = one_over_eight * (1 - xi1) * (1 + xi3)
76 # Derivatives with respect to xi3
77 derivatives[2, 0] = -one_over_eight * (1 - xi1) * (1 - xi2)
78 derivatives[2, 1] = -one_over_eight * (1 + xi1) * (1 - xi2)
79 derivatives[2, 2] = -one_over_eight * (1 + xi1) * (1 + xi2)
80 derivatives[2, 3] = -one_over_eight * (1 - xi1) * (1 + xi2)
81 derivatives[2, 4] = one_over_eight * (1 - xi1) * (1 - xi2)
82 derivatives[2, 5] = one_over_eight * (1 + xi1) * (1 - xi2)
83 derivatives[2, 6] = one_over_eight * (1 + xi1) * (1 + xi2)
84 derivatives[2, 7] = one_over_eight * (1 - xi1) * (1 + xi2)
86 return derivatives
89def get_hex8_element_center_and_jacobian_mapping(element):
90 """Return the center of a hex8 element and the Jacobian mapping for that point."""
91 nodal_coordinates = _get_nodal_coordinates(element.nodes)
92 if not len(nodal_coordinates) == 8:
93 raise ValueError(f"Expected 8 nodes, got {len(nodal_coordinates)}")
95 N = shape_functions_hex8(0, 0, 0)
96 dN = shape_functions_derivative_hex8(0, 0, 0)
98 reference_position_center = _np.dot(N, nodal_coordinates)
99 jacobian_center = _np.dot(dN, nodal_coordinates)
101 return reference_position_center, jacobian_center
104def get_reordering_index_thickness(jacobian, *, identify_threshold=None):
105 """Return the reordering index from the Jacobian such that the thinnest direction is
106 the 3rd parameter direction.
108 Additionally it is checked, that the thinnest direction is at least
109 identify_threshold times smaller than the next thinnest, to avoid wrongly detected
110 directions.
111 """
112 # The direction with the smallest parameter derivative is the thickness direction
113 parameter_derivative_norms = [
114 _np.linalg.norm(parameter_direction) for parameter_direction in jacobian
115 ]
116 thickness_direction = _np.argmin(parameter_derivative_norms)
118 if identify_threshold is not None:
119 # Check that the minimal parameter direction is at least a given factor
120 # smaller than the other directions. This helps to identify cases where
121 # it is unlikely that a unique direction is found.
122 min_norm = parameter_derivative_norms[thickness_direction]
123 relative_difference = [
124 parameter_derivative_norms[i] / min_norm
125 for i in range(3)
126 if not i == thickness_direction
127 ]
128 if _np.min(relative_difference) < 1.5:
129 raise ValueError("Could not uniquely identify the thickness direction.")
131 return thickness_direction
134def get_reordering_index_director_projection(
135 jacobian, director, *, identify_threshold=None
136):
137 """Return the reordering index from the Jacobian such that the thickness direction
138 is the one that has the largest dot product with the given director."""
139 projections = []
140 for parameter_director in jacobian:
141 parameter_director = parameter_director / _np.linalg.norm(parameter_director)
142 projections.append(_np.abs(_np.dot(director, parameter_director)))
143 thickness_direction = _np.argmax(projections)
145 if identify_threshold is not None:
146 # Check that the maximal dot product is at least a given factor larger than
147 # the other projections. This helps to identify cases where it is unlikely
148 # that a unique direction is found.
149 max_norm = projections[thickness_direction]
150 relative_difference = [
151 projections[i] / max_norm for i in range(3) if not i == thickness_direction
152 ]
153 if _np.max(relative_difference) > 1.0 / identify_threshold:
154 raise ValueError("Could not uniquely identify the thickness direction.")
156 return thickness_direction
159def set_solid_shell_thickness_direction(
160 elements: list[_Element],
161 *,
162 selection_type="thickness",
163 director=None,
164 director_function=None,
165 identify_threshold=2.0,
166):
167 """Set the solid shell directions for all solid shell elements in the element list.
169 Args:
170 ----
171 elements:
172 A list containing all elements that should be checked
173 selection_type:
174 The type of algorithm that shall be used to select the thickness direction
175 "thickness":
176 The "smallest" dimension of the element will be set to the
177 thickness direction
178 "projection_director":
179 The parameter director that aligns most with a given director
180 will be set as the thickness direction
181 "projection_director_function":
182 The parameter director that aligns most with a director obtained
183 by a given function of the element centroid coordinates will be
184 set as the thickness direction
185 identify_threshold: float/None
186 To ensure that the found directions are well-defined, i.e., that not multiple
187 directions are almost equally suited to the thickness direction.
188 """
189 if len(elements) == 0:
190 raise ValueError("Expected a non empty element list")
192 for element in elements:
193 four_c_data = getattr(type(element), "data", None)
194 if four_c_data is None:
195 raise ValueError(
196 "Expected element type to have a data attribute with the 4C element data!"
197 )
198 is_hex8_solid_shell = (
199 four_c_data.four_c_type == "SOLID"
200 and four_c_data.four_c_cell == "HEX8"
201 and four_c_data.element_technology.get("TECH", "") == "shell_eas_ans"
202 )
203 if is_hex8_solid_shell:
204 # Get the element center and the Jacobian at the center
205 (
206 reference_position_center,
207 jacobian_center,
208 ) = get_hex8_element_center_and_jacobian_mapping(element)
210 # Depending on the chosen method, get the thickness direction
211 if selection_type == "thickness":
212 thickness_direction = get_reordering_index_thickness(
213 jacobian_center, identify_threshold=identify_threshold
214 )
215 elif selection_type == "projection_director":
216 thickness_direction = get_reordering_index_director_projection(
217 jacobian_center, director, identify_threshold=identify_threshold
218 )
219 elif selection_type == "projection_director_function":
220 director = director_function(reference_position_center)
221 thickness_direction = get_reordering_index_director_projection(
222 jacobian_center, director, identify_threshold=identify_threshold
223 )
224 else:
225 raise ValueError(
226 f'Got unexpected selection_type of value "{selection_type}"'
227 )
229 n_apply_mapping = 0
230 if thickness_direction == 2:
231 # We already have the orientation we want
232 continue
233 elif thickness_direction == 1:
234 # We need to apply the connectivity mapping once, i.e., the 2nd parameter
235 # direction has to become the 3rd
236 n_apply_mapping = 1
237 elif thickness_direction == 0:
238 # We need to apply the connectivity mapping twice, i.e., the 1nd parameter
239 # direction has to become the 3rd
240 n_apply_mapping = 2
242 # This permutes the parameter coordinate the following way:
243 # [xi,eta,zeta]->[zeta,xi,eta]
244 mapping = [2, 6, 7, 3, 1, 5, 4, 0]
245 for _ in range(n_apply_mapping):
246 element.nodes = [element.nodes[local_index] for local_index in mapping]
249def get_visualization_third_parameter_direction_hex8(mesh: _Mesh):
250 """Return a pyvista mesh with cell data for the third parameter direction for hex8
251 elements."""
252 grid = mesh.get_vtu_representation()
253 grid_solid = grid.extract_cells(
254 grid.cell_data["beamme_type"] == _bme.element_type.solid.value
255 )
257 cell_thickness_direction = []
258 for element in mesh.elements:
259 if isinstance(element, _VolumeHEX8):
260 _, jacobian_center = get_hex8_element_center_and_jacobian_mapping(element)
261 cell_thickness_direction.append(jacobian_center[2])
262 elif isinstance(element, _VolumeElement):
263 cell_thickness_direction.append([0, 0, 0])
265 if not len(cell_thickness_direction) == grid_solid.number_of_cells:
266 raise ValueError(
267 "Expected the same number of cells from the mesh and from the "
268 f"pyvista object. Got {len(cell_thickness_direction)} form the "
269 f"mesh and {grid_solid.number_of_cells} form pyvista"
270 )
272 grid_solid.cell_data["thickness_direction"] = cell_thickness_direction
273 return grid_solid
276def visualize_third_parameter_direction_hex8(mesh: _Mesh):
277 """Visualize the third parameter direction for hex8 elements.
279 This can be used to check the correct definition of the shell thickness for solid
280 shell elements.
281 """
282 grid = get_visualization_third_parameter_direction_hex8(mesh)
283 grid = grid.clean()
284 cell_centers = grid.cell_centers()
285 thickness_direction = cell_centers.glyph(
286 orient="thickness_direction", scale="thickness_direction", factor=5
287 )
289 plotter = _pv.Plotter()
290 plotter.renderer.add_axes()
291 plotter.add_mesh(grid, color="white", show_edges=True, opacity=0.5)
292 plotter.add_mesh(thickness_direction, color="red")
293 plotter.show()