Coverage for src/meshpy/core/coupling.py: 91%

1# The MIT License (MIT) 

2# 

3# Copyright (c) 2018-2025 MeshPy 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 module implements a class to couple geometry together.""" 

23 

24from typing import List as _List 

25from typing import Union as _Union 

26 

27import numpy as _np 

28 

29import meshpy.core.conf as _conf 

30from meshpy.core.boundary_condition import ( 

31 BoundaryConditionBase as _BoundaryConditionBase, 

32) 

33from meshpy.core.conf import mpy as _mpy 

34from meshpy.core.geometry_set import GeometrySet as _GeometrySet 

35from meshpy.core.geometry_set import GeometrySetBase as _GeometrySetBase 

36from meshpy.core.node import Node as _Node 

37 

38 

39class Coupling(_BoundaryConditionBase): 

40 """Represents a coupling between geometries in 4C.""" 

41 

42 def __init__( 

43 self, 

44 geometry: _Union[_GeometrySetBase, _List[_Node]], 

45 coupling_type: _Union[_conf.BoundaryCondition, str], 

46 coupling_dof_type, 

47 *, 

48 check_overlapping_nodes: bool = True, 

49 **kwargs, 

50 ): 

51 """Initialize this object. 

52 

53 Args: 

54 geometry: Geometry set or nodes that should be coupled. 

55 coupling_type: Type of coupling. 

56 coupling_dof_type: mpy.coupling_dof, str 

57 If this is a string it is the string that will be used in the input 

58 file, otherwise it has to be of type mpy.coupling_dof. 

59 check_overlapping_nodes: If all nodes of this coupling condition 

60 have to be at the same physical position. 

61 """ 

62 

63 if isinstance(geometry, _GeometrySetBase): 

64 pass 

65 elif isinstance(geometry, list): 

66 geometry = _GeometrySet(geometry) 

67 else: 

68 raise TypeError( 

69 f"Coupling expects a GeometrySetBase item, got {type(geometry)}" 

70 ) 

71 

72 # Couplings only work for point sets 

73 if ( 

74 isinstance(geometry, _GeometrySetBase) 

75 and geometry.geometry_type is not _mpy.geo.point 

76 ): 

77 raise TypeError("Couplings are only implemented for point sets.") 

78 

79 super().__init__(geometry, bc_type=coupling_type, **kwargs) 

80 self.coupling_dof_type = coupling_dof_type 

81 self.check_overlapping_nodes = check_overlapping_nodes 

82 

83 # Perform sanity checks for this boundary condition 

84 self.check() 

85 

86 def check(self): 

87 """Check that all nodes that are coupled have the same position 

88 (depending on the check_overlapping_nodes parameter).""" 

89 

90 if not self.check_overlapping_nodes: 

91 return 

92 

93 nodes = self.geometry_set.get_points() 

94 diff = _np.zeros([len(nodes), 3]) 

95 for i, node in enumerate(nodes): 

96 # Get the difference to the first node 

97 diff[i, :] = node.coordinates - nodes[0].coordinates 

98 if _np.max(_np.linalg.norm(diff, axis=1)) > _mpy.eps_pos: 

99 raise ValueError( 

100 "The nodes given to Coupling do not have the same position." 

101 ) 

102 

103 def dump_to_list(self): 

104 """Return a list with a single item representing this coupling 

105 condition.""" 

106 

107 if isinstance(self.coupling_dof_type, dict): 

108 data = self.coupling_dof_type 

109 else: 

110 # In this case we have to check which beams are connected to the node. 

111 # TODO: Coupling also makes sense for different beam types, this can 

112 # be implemented at some point. 

113 nodes = self.geometry_set.get_points() 

114 beam_type = nodes[0].element_link[0].beam_type 

115 for node in nodes: 

116 for element in node.element_link: 

117 if beam_type is not element.beam_type: 

118 raise ValueError( 

119 f'The first element in this coupling is of the type "{beam_type}" ' 

120 f'another one is of type "{element.beam_type}"! They have to be ' 

121 "of the same kind." 

122 ) 

123 if beam_type is _mpy.beam.kirchhoff and element.rotvec is False: 

124 raise ValueError( 

125 "Couplings for Kirchhoff beams and rotvec==False not yet implemented." 

126 ) 

127 

128 # In 4C it is not possible to couple beams of the same type, but 

129 # with different centerline discretizations, e.g. Beam3rHerm2Line3 

130 # and Beam3rLine2Line2, therefore, we check that all beams are 

131 # exactly the same type and discretization. 

132 # TODO: Remove this check once it is possible to couple beams, but 

133 # then also the syntax in the next few lines has to be adapted. 

134 beam_four_c_type = type(nodes[0].element_link[0]) 

135 for node in nodes: 

136 for element in node.element_link: 

137 if beam_four_c_type is not type(element): 

138 raise ValueError( 

139 "Coupling beams of different types is not yet possible!" 

140 ) 

141 

142 data = beam_four_c_type.get_coupling_dict(self.coupling_dof_type) 

143 

144 return [{"E": self.geometry_set.i_global, **data}] 

145 

146 

147def coupling_factory(geometry, coupling_type, coupling_dof_type, **kwargs): 

148 """Create coupling conditions for the nodes in geometry. 

149 

150 Some solvers only allow coupling conditions containing two points at 

151 once, in that case we have to create multiple coupling conditions 

152 between the individual points to ensure the correct representation 

153 of the coupling. 

154 """ 

155 

156 if coupling_type is _mpy.bc.point_coupling_penalty: 

157 # Penalty point couplings in 4C can only contain two nodes. In this case 

158 # we expect the given geometry to be a list of nodes. 

159 main_node = geometry[0] 

160 return [ 

161 Coupling([main_node, node], coupling_type, coupling_dof_type, **kwargs) 

162 for node in geometry[1:] 

163 ] 

164 else: 

165 return [Coupling(geometry, coupling_type, coupling_dof_type, **kwargs)]