-
Notifications
You must be signed in to change notification settings - Fork 0
/
StructuralLoss.py
307 lines (258 loc) · 19.6 KB
/
StructuralLoss.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
import torch
from models.layers.mesh import Mesh
from utils import save_mesh, export_vector
DOF = 6 # degrees of freedom per vertex
class StructuralLoss:
def __init__(self, file=None, mesh=None, beam_properties=None, beam_have_load=False, device='cpu'):
self.device = torch.device(device)
self.set_beam_properties(beam_properties)
self.beam_have_load = beam_have_load
if file is not None:
self.initial_mesh = Mesh(file, device=device)
elif mesh is not None:
self.initial_mesh = mesh
else:
raise ValueError('No reference mesh specified.')
self.initialize_containers()
self.beam_model_solve(self.initial_mesh)
def __call__(self, mesh, loss_type):
self.beam_model_solve(mesh)
# Giving right loss corresponding with loss_type.
if loss_type == 'norm_vertex_deformations':
return torch.sum(torch.norm(self.vertex_deformations[:, :3], p=2, dim=1))
elif loss_type == 'mean_beam_energy':
return torch.mean(self.beam_energy)
elif loss_type == 'no_axial_mean_beam_energy':
return torch.mean(self.beam_energy_without_axial)
# Store beam properties involved in the task.
# Custom properties are passed through a list whose elements follow this order:
#
# -- [0] Poisson's ratio, default is 0.3;
# -- [1] Young's modulus, default is 2.1*10^8 kN/m^2;
# -- [2] beam section area, default is 1*10^-3 m^2;
# -- [3] moment of inertia2 Ixx = Iyy, default is 4.189828*10^-8, round section;
# -- [4] moment of inertia3 Ixx = Iyy, default is 4.189828*10^-8, round section;
# -- [5] polar moment, default is 8.379656e-8;
# -- [6] shear section factor, default is 1.2;
# -- [7] weight per surface unit, default is 3 kN/m^2;
# -- [8] beam density, default is 78.5 kN/m^3.
#
def set_beam_properties(self, beam_properties):
if beam_properties is not None:
self.properties = beam_properties
else:
self.properties = [0.3, 21e7, 1e-3, 4.189828e-8, 4.189828e-8, 8.379656e-8, 1.2, 3, 78.5]
# Computing G := young/(2 * (1 + poisson)).
self.properties.append(self.properties[1] / (2 * (1+self.properties[0])))
# Converting to torch.tensor on current device.
self.properties = torch.tensor(self.properties, device=self.device)
# Re-usable tensors are initialized just once at the beginning.
def initialize_containers(self):
# Beam local frames container: (#edges, 3, 3) torch.tensor.
self.beam_frames = torch.zeros(self.initial_mesh.edges.shape[0], 3, 3, device=self.device)
# Beam stiff matrices container: (#edges, 2*DOF, 2*DOF) torch.tensor.
self.beam_stiff_matrices = torch.zeros(self.initial_mesh.edges.shape[0], 2*DOF, 2*DOF, device=self.device)
# Per-vertex loads vector.
self.load = torch.zeros(self.initial_mesh.vertices.shape[0] * DOF, device=self.device)
###########################################################################################################
# Storing augmented stiffmatrix non-zero indices.
self.augmented_stiff_idx = self.make_augmented_stiffmatrix_nnz_indices()
###########################################################################################################
# Bulding bool tensor masks related to constraints.
# Non-constrained vertex mask:
self.non_constrained_vertices = torch.logical_or(self.initial_mesh.vertex_is_red, self.initial_mesh.vertex_is_blue).logical_not()
# Red vertices refer to all dofs constrainess, blue vertices to just translation constrainess.
blue_dofs = torch.kron(self.initial_mesh.vertex_is_blue, torch.tensor([True] * int(DOF/2) + [False] * int(DOF/2), device=self.device))
red_dofs = torch.kron(self.initial_mesh.vertex_is_red, torch.tensor([True] * DOF, device=self.device))
# Non-constrained vertex dofs mask.
self.dofs_non_constrained_mask = torch.logical_and(blue_dofs.logical_not(), red_dofs.logical_not())
# Please note, + operator in this case makes row-wise arranged torch.arange(DOF) copies each time shiftedy by
# corresponding 6 * self.mesh.edges[:, 0].unsqueeze(-1) row element.
# (...).unsqueeze(-1) expand (#edges, ) tensor to (#edges, 1) column tensor.
def make_edge_endpoints_dofs_matrix(self):
endpts_0 = 6 * self.initial_mesh.edges[:, 0].unsqueeze(-1) + torch.arange(DOF, device=self.device)
endpts_1 = 6 * self.initial_mesh.edges[:, 1].unsqueeze(-1) + torch.arange(DOF, device=self.device)
return torch.cat((endpts_0, endpts_1), dim=1)
# Builds augmented stiffmatrix non-zero element tensor according to vertex dofs.
def make_augmented_stiffmatrix_nnz_indices(self):
self.endpoints_dofs_matrix = self.make_edge_endpoints_dofs_matrix()
dim1_idx = self.endpoints_dofs_matrix.view(-1, 1).expand(-1, 2*DOF).flatten()
dim2_idx = self.endpoints_dofs_matrix.expand(2*DOF, -1, -1).transpose(0, 1).flatten()
return torch.stack([dim1_idx, dim2_idx], dim=0)
# Stores load matrix and builds beam frames.
# Load matrix (no. vertices x DOF) has row-structure (Fx, Fy, Fz, Mx, My, Mz) whose values are referred to global ref system.
def set_loads_and_beam_frames(self, mesh):
#####################################################################################################################################
# Load computation.
# Computing beam loads on each vertex: -1/2 * <beam volume> * <beam density> (beam load is equally parted to endpoints) along z-axis.
# Some details: vec.scatter_add(0, idx, src) with vec, idx, src 1d tensors, add at vec positions specified
# by idx corresponding src values.
if self.beam_have_load:
beam_loads = torch.mul(mesh.edge_lengths, -1/2 * self.properties[2] * self.properties[8])
on_vertices_beam_loads = torch.zeros(mesh.vertices.shape[0], device=self.device)
on_vertices_beam_loads.scatter_add_(0, mesh.edges.flatten(), torch.stack([beam_loads] * 2, dim=1).flatten())
# Computing face loads on each vertex: -1/3 * <face areas> * <weight per surface unit> (face load is equally parted to vertices) along z.
# Some details: vec.scatter_add(0, idx, src) with vec, idx, src 1d tensors, add at vec positions specified
# by idx corresponding src values.
face_loads = torch.mul(mesh.face_areas, -1/3 * self.properties[7])
on_vertices_face_loads = torch.zeros(mesh.vertices.shape[0], device=self.device)
on_vertices_face_loads.scatter_add_(0, mesh.faces.flatten(), torch.stack([face_loads] * 3, dim=1).flatten())
# Summing beam and face components to compute per-vertex loads.
if self.beam_have_load:
self.load[DOF*torch.arange(mesh.vertices.shape[0], device=self.device) + 2] = on_vertices_beam_loads + on_vertices_face_loads
else:
self.load[DOF*torch.arange(mesh.vertices.shape[0], device=self.device) + 2] = on_vertices_face_loads
#######################################################################################################################################
# Beam frames computation.
self.beam_frames[:, 0, :] = mesh.edge_directions
self.beam_frames[:, 1, :] = mesh.edge_normals
self.beam_frames[:, 2, :] = mesh.cross_dirs
# Execute all stiffness and resistence computations.
def beam_model_solve(self, mesh):
self.set_loads_and_beam_frames(mesh)
self.build_stiff_matrix(mesh)
self.compute_stiff_deformation(mesh)
# Stiffness matrices in beam reference systems are computed and then aggregated to compound a global stiff matrix.
def build_stiff_matrix(self, mesh):
###########################################################################################################
# Assembling beam-local stiff matrices whose have the following structure:
# k1 = properties.young * properties.cross_area / self.mesh.edge_lengths[edge_id]
# k2 = 12 * properties.young * properties.inertia2 / (self.mesh.edge_lengths[edge_id]**3)
# k3 = 12 * properties.young * properties.inertia3 / (self.mesh.edge_lengths[edge_id]**3)
# k4 = 6 * properties.young * properties.inertia3 / (self.mesh.edge_lengths[edge_id]**2)
# k5 = 6 * properties.young * properties.inertia2 /(self.mesh.edge_lengths[edge_id]**2)
# k6 = self.G * properties.polar / self.mesh.edge_lengths[edge_id]
# k7 = properties.young * properties.inertia2 / self.mesh.edge_lengths[edge_id]
# k8 = properties.young * properties.inertia3 / self.mesh.edge_lengths[edge_id]
# k_e = [[k1, 0, 0, 0, 0, 0, -k1, 0, 0, 0, 0, 0],
# [0, k3, 0, 0, 0, k4, 0, -k3, 0, 0, 0, k4],
# [0, 0, k2, 0, -k5, 0, 0, 0, -k2, 0, -k5, 0],
# [0, 0, 0, k6, 0, 0, 0, 0, 0, -k6, 0, 0],
# [0, 0, -k5, 0, 4*k7, 0, 0, 0, k5, 0, 2*k7, 0],
# [0, k4, 0, 0, 0, 4*k8, 0, -k4, 0, 0, 0, 2*k8],
# [-k1, 0, 0, 0, 0, 0, k1, 0, 0, 0, 0, 0],
# [0, -k3, 0, 0, 0, -k4, 0, k3, 0, 0, 0, -k4],
# [0, 0, -k2, 0, k5, 0, 0, 0, k2, 0, k4, 0],
# [0, 0, 0, -k6, 0, 0, 0, 0, 0, k6, 0, 0],
# [0, 0, -k5, 0, 2*k7, 0, 0, 0, k5, 0, 4*k7, 0],
# [0, k4, 0, 0, 0, 2*k8, 0, -k4, 0, 0, 0, 4*k8]]
#
# Computing squares and cubes of beam_lenghts tensor.
squared_beam_lenghts = torch.pow(mesh.edge_lengths, 2)
cubed_beam_lenghts = torch.pow(mesh.edge_lengths, 3)
# Filling non empty entries in self.beam_stiff_matrices.
# k1 and -k1
self.beam_stiff_matrices[:, 0, 0] = self.properties[1] * self.properties[2] / mesh.edge_lengths
self.beam_stiff_matrices[:, 6, 6] = self.beam_stiff_matrices[:, 0, 0]
self.beam_stiff_matrices[:, 6, 0] = self.beam_stiff_matrices[:, 0, 6] = -self.beam_stiff_matrices[:, 0, 0]
# k2 and -k2
self.beam_stiff_matrices[:, 2, 2] = 12 * self.properties[1] * self.properties[3] / cubed_beam_lenghts
self.beam_stiff_matrices[:, 8, 8] = self.beam_stiff_matrices[:, 2, 2]
self.beam_stiff_matrices[:, 8, 2] = self.beam_stiff_matrices[:, 2, 8] = -self.beam_stiff_matrices[:, 2, 2]
# k3 and -k3
self.beam_stiff_matrices[:, 1, 1] = 12 * self.properties[1] * self.properties[4] / cubed_beam_lenghts
self.beam_stiff_matrices[:, 7, 7] = self.beam_stiff_matrices[:, 1, 1]
self.beam_stiff_matrices[:, 7, 1] = self.beam_stiff_matrices[:, 1, 7] = -self.beam_stiff_matrices[:, 1, 1]
# k4 and -k4
self.beam_stiff_matrices[:, 5, 1] = 6 * self.properties[1] * self.properties[4] / squared_beam_lenghts
self.beam_stiff_matrices[:, 1, 5] = self.beam_stiff_matrices[:, 11, 1] = self.beam_stiff_matrices[:, 5, 1]
self.beam_stiff_matrices[:, 1, 11] = self.beam_stiff_matrices[:, 8, 10] = self.beam_stiff_matrices[:, 5, 1]
self.beam_stiff_matrices[:, 7, 5] = self.beam_stiff_matrices[:, 11, 7] = - self.beam_stiff_matrices[:, 5, 1]
self.beam_stiff_matrices[:, 5, 7] = self.beam_stiff_matrices[:, 7, 11] = - self.beam_stiff_matrices[:, 5, 1]
# k5 and -k5
self.beam_stiff_matrices[:, 8, 4] = 6 * self.properties[1] * self.properties[3] / squared_beam_lenghts
self.beam_stiff_matrices[:, 4, 8] = self.beam_stiff_matrices[:, 10, 8] = self.beam_stiff_matrices[:, 8, 4]
self.beam_stiff_matrices[:, 4, 2] = self.beam_stiff_matrices[:, 10, 2] = -self.beam_stiff_matrices[:, 8, 4]
self.beam_stiff_matrices[:, 2, 4] = self.beam_stiff_matrices[:, 2, 10] = -self.beam_stiff_matrices[:, 8, 4]
# k6 and -k6
self.beam_stiff_matrices[:, 3, 3] = self.properties[9] * self.properties[5] / mesh.edge_lengths
self.beam_stiff_matrices[:, 9, 9] = self.beam_stiff_matrices[:, 3, 3]
self.beam_stiff_matrices[:, 9, 3] = self.beam_stiff_matrices[:, 3, 9] = -self.beam_stiff_matrices[:, 3, 3]
# k7 and multiples
k7 = self.properties[1] * self.properties[3] / mesh.edge_lengths
self.beam_stiff_matrices[:, 4, 4] = self.beam_stiff_matrices[:, 10, 10] = torch.mul(4, k7)
self.beam_stiff_matrices[:, 10, 4] = self.beam_stiff_matrices[:, 4, 10] = torch.mul(2, k7)
# k8 and multiples
k8 = self.properties[1] * self.properties[4] / mesh.edge_lengths
self.beam_stiff_matrices[:, 5, 5] = self.beam_stiff_matrices[:, 11, 11] = torch.mul(4, k8)
self.beam_stiff_matrices[:, 11, 5] = self.beam_stiff_matrices[:, 5, 11] = torch.mul(2, k8)
###########################################################################################################
# Assembling beam-local to global transition matrices via Kronecker product: container is again a
# 3-dimensional torch.tensor.
transition_matrices = torch.kron(torch.eye(4, 4, device=self.device), self.beam_frames)
###########################################################################################################
# Building beam contributions to global stiff matrix: container (beam_contributions) is a 3d torch.tensor,
# another contanier (self.beam_forces_contributions) for products beam_stiff_matrix @ transition_matrix
# is saved in order not to repeat steps in node forces computation phase.
# Please note: @ operator for 3d tensors produces a 'batched' 2d matrix multiplication along sections.
self.beam_forces_contributions = self.beam_stiff_matrices @ transition_matrices
beam_contributions = torch.transpose(transition_matrices, 1, 2) @ self.beam_forces_contributions
#torch.transpose(_, 1, 2) transpose along dimensions 1 and 2
###########################################################################################################
# Building global stiff matrix by adding all beam contributions.
# Global stiff matrix: (DOF*#vertices, DOF*#vertices), computed by densing an augmented sparse
# (DOF*#vertices, DOF*#vertices), admitting several entries for the same ij coordinate couple.
size = (DOF * mesh.vertices.shape[0], DOF * mesh.vertices.shape[0])
augmented_stiff_matrix = torch.sparse_coo_tensor(self.augmented_stiff_idx, beam_contributions.flatten(), size=size, device=self.device)
self.stiff_matrix = augmented_stiff_matrix.to_dense()
# Freeing memory space.
del transition_matrices, beam_contributions, augmented_stiff_matrix, squared_beam_lenghts, cubed_beam_lenghts
# Compute vertex deformations by solving a stiff-matrix-based linear system.
def compute_stiff_deformation(self, mesh):
# Solving reduced linear system and adding zeros in constrained dofs.
self.vertex_deformations = torch.zeros(mesh.vertices.shape[0] * DOF, device=self.device)
sys_sol = torch.linalg.solve(self.stiff_matrix[self.dofs_non_constrained_mask][:, self.dofs_non_constrained_mask], self.load[self.dofs_non_constrained_mask])
self.vertex_deformations[self.dofs_non_constrained_mask] = sys_sol
# Freeing memory space.
del self.stiff_matrix, sys_sol
# Computing beam resulting forces and energies.
self.compute_beam_force_and_energy(mesh)
# Making deformation tensor by reshaping self.vertex_deformations.
self.vertex_deformations = self.vertex_deformations.view(mesh.vertices.shape[0], DOF)
# Computing beam resulting forces and energies.
def compute_beam_force_and_energy(self, mesh):
# edge_dofs_deformations aggregates vertex_deformation in a edge-endpoints-wise manner: (#edges, 2*DOF) torch.tensor
edge_dofs_deformations = self.vertex_deformations[self.endpoints_dofs_matrix]
#############################################################################################################################
# Computing resulting forces at nodes via 'batched' matrix multiplication @.
# Some details:
# edge_dofs_deformations.unsqueeze(2) expands (#edges, 2*DOF) to #edge batches of (2*DOF, 1), i.e. (#edges, 2*DOF, 1) tensor;
# vice-versa, (...).squeeze(2) contracts (#edges, 2*DOF, 1) -> (#edges, 2*DOF)
node_forces = (self.beam_forces_contributions @ edge_dofs_deformations.unsqueeze(2)).squeeze(2)
# Averaging force components 0:DOF with D0F:2*DOF of opposite sign.
mean_forces = torch.mul(0.5, node_forces[:, :DOF] - node_forces[:, DOF:2*DOF])
# Computing beam energy according this coordinate system:
# axes: 1=elementAxis; 2=EdgeNormal; 3=InPlaneAxis
# output rows: [Axial_startNode; Shear2_startNode; Shear3_startNode; Torque_startNode; Bending3_startNode; Bending2_startNode;
# ...Axial_endNode; Shear2_endNode; Shear3_endNode; Torque_endNode; Bending3_endNode; Bending2_endNode]
self.beam_energy_without_axial = mesh.edge_lengths/2 * (self.properties[6] * mean_forces[:, 1]**2 / (self.properties[9] * self.properties[2]) +
self.properties[6] * mean_forces[:, 2]**2 / (self.properties[9] * self.properties[2]) +
mean_forces[:, 3]**2 / (self.properties[9] * self.properties[5]) +
mean_forces[:, 4]**2 / (self.properties[1] * self.properties[4]) +
mean_forces[:, 5]**2 / (self.properties[1] * self.properties[3]) )
self.beam_energy = self.beam_energy_without_axial + mesh.edge_lengths/2 * mean_forces[:, 0]**2/(self.properties[1]*self.properties[2])
# Freeing memory space.
del edge_dofs_deformations, self.beam_forces_contributions
def clean_attributes(self):
self.load.detach_()
self.beam_frames.detach_()
self.beam_stiff_matrices.detach_()
# Displace initial mesh with self.beam_model_solve() computed translations.
def stress_mesh(self):
# REQUIRES: self.beam_model_solve() has to be called before executing this.
if not hasattr(self, 'vertex_deformations'):
raise RuntimeError("self.beam_model_solve() method not called yet.")
# Updating mesh vertices.
stressed_mesh = self.initial_mesh.update_verts(self.vertex_deformations[:, :int(DOF/2)])
return stressed_mesh
def save_grid_shell(self, mesh):
quality = torch.norm(self.vertex_deformations[:, :3], p=2, dim=1)
save_mesh(mesh, 'stressed_mesh.ply', v_quality=quality.unsqueeze(1))
#############################################################################################################
# TEST StructuralLoss.py ###
if __name__ == '__main__':
lc = StructuralLoss(file='meshes/casestudy_compr.ply', device='cpu')
lc.save_grid_shell(lc.initial_mesh)
export_vector(torch.norm(lc.vertex_deformations[:, :3], dim=1), 'load.csv')
export_vector(lc.initial_mesh.edges, 'edges.csv')
export_vector(lc.beam_energy, 'energy.csv')