beef.general
General purpose functions.
1''' 2General purpose functions. 3''' 4 5import numpy as np 6from scipy.linalg import null_space as null 7 8 9def n2d_ix(node_ixs, n_dofs=6): 10 r''' 11 Convert node indices to dof indices. 12 13 Arguments 14 ---------- 15 node_ixs : int 16 list/array of indices of nodes 17 n_dofs : 6, optional 18 number of dofs per node 19 20 Returns 21 ---------- 22 dof_ixs : int 23 array of DOF indices corresponding to input node indices 24 ''' 25 return np.hstack([node_ix*n_dofs+np.arange(n_dofs) for node_ix in node_ixs]) 26 27 28def extract_dofs(mat, dof_ix=[0,1,2], n_dofs=6, axis=0): 29 r''' 30 Extract selected dofs from matrix. 31 32 Arguments 33 ------------ 34 mat : float 35 numpy array (matrix) to extract DOFs from 36 dof_ix : int, optional 37 list of DOF indices (local numbering of DOFs in node) to extract 38 n_dofs : 6, optional 39 number of DOFs per node 40 axis : 0, optional 41 axis to pick along 42 43 Returns 44 ---------- 45 mat_sub : float 46 subselection of matrix 47 ''' 48 get_dofs = np.vstack([np.arange(dof, mat.shape[axis],n_dofs) for dof in dof_ix]).T.flatten() 49 return mat.take(get_dofs, axis=axis) 50 51def convert_dofs(dofs_in, node_type='beam3d', sort_output=True): 52 r''' 53 Convert string DOFs to indices. 54 55 Arguments 56 ---------- 57 dofs_in : {'all', 'trans', 'rot'} or int 58 string or list of indices 59 node_type : {'beam2d', 'beam3d'}, optional 60 type of nodes (defining number of DOFs per node) 61 sort_output : True, optional 62 whether or not to sort the DOF indices output 63 64 Returns 65 -------- 66 dofs_out : int 67 array of DOF indices 68 ''' 69 dof_dict = dict(beam2d=dict(all=np.arange(0,3), trans=np.arange(0,2), rot=np.arange(2,3)), 70 beam3d=dict(all=np.arange(0,6), trans=np.arange(0,3), rot=np.arange(3,6))) 71 72 if dofs_in in ['all', 'trans', 'rot']: 73 dofs_out = dof_dict[node_type][dofs_in] 74 else: 75 dofs_out = dofs_in 76 77 if sort_output: 78 dofs_out = np.sort(dofs_out) 79 80 return dofs_out 81 82 83def convert_dofs_list(dofs_list_in, n_nodes, node_type='beam3d', sort_output=True): 84 r''' 85 Convert DOFs list to correct format. 86 87 Arguments 88 ----------- 89 dofs_list_in : str,int 90 list of strings and indices corresponding to DOFs to 91 n_nodes : int 92 number of nodes the DOFs are related to 93 node_type : {'beam3d', 'beam2d'} 94 type of nodes (defining number of DOFs per node) 95 sort_output : True 96 whether or not to sort output 97 98 Returns 99 ----------- 100 dofs_list_out : int 101 array of DOF indices 102 ''' 103 104 contains_strings = np.any([type(d) is str for d in dofs_list_in]) 105 106 if type(dofs_list_in) is not list: # single input (all, rot or trans) 107 dofs_list_out = [dofs_list_in]*n_nodes 108 elif ~contains_strings and (len(dofs_list_in)<=6) and (np.max(dofs_list_in)<6): # 109 dofs_list_out = [dofs_list_in]*n_nodes 110 elif len(dofs_list_in)!=n_nodes: 111 raise TypeError('Wrong input format of "dofs"') 112 else: 113 dofs_list_out = dofs_list_in 114 115 for ix, dofs_list_out_i in enumerate(dofs_list_out): 116 dofs_list_out[ix] = convert_dofs(dofs_list_out_i, node_type=node_type, sort_output=sort_output) 117 118 return dofs_list_out 119 120 121def transform_unit(e1, e2=None, e3=None, warnings=False): 122 r''' 123 Establish transformation matrix from e1 and temporary e2 or e3 vectors. 124 125 Arguments 126 ----------- 127 e1 : float 128 unit vector describing element longitudinal direction 129 e2 : float, optional 130 temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate y-direction) 131 e3 : float, optional 132 temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate z-direction) 133 if both e2 and e3 are different from None, e2 is used (e3 disregarded) 134 135 Returns 136 ----------- 137 T : float 138 transformation matrix 139 140 141 Notes 142 ----------- 143 Input vectors \(\{e_1\}\) and \(\{e_{2p}\}\) are used to establish a third unit vector, as the cross product of the two, as follows: 144 145 $$ 146 \{e_3\} = \{e_1\} \\times \{e_{2p}\} 147 $$ 148 149 Then, a final vector is created based on a second cross product: 150 $$ 151 \{e_2\} = \{e_3\} \\times \{e_1\} 152 $$ 153 154 Finally, the transformation matrix is established as follows: 155 156 $$ 157 [T] = \\begin{bmatrix} 158 \{e_1\}^T \\\\ 159 \{e_2\}^T \\\\ 160 \{e_3\}^T 161 \\end{bmatrix} 162 $$ 163 164 where the unit vectors established are the rows of the transformation matrix. 165 166 ''' 167 168 e1 = np.array(e1).flatten() 169 170 if (e2 is not None) and (e3 is not None): 171 e3 = None 172 173 if e2 is not None: 174 e2 = np.array(e2).flatten() 175 e3_tmp = np.cross(e1, e2) # Direction of the third unit vector 176 if np.all(e3==0): 177 if warnings: 178 print('Warning: e1 and e2 identical. Check orientations carefully!') 179 180 e2 = e2 + 0.1 181 e3 = np.cross(e1, e2) 182 183 e2 = np.cross(e3_tmp, e1) # Direction of the second unit vector 184 185 e1 = e1/np.linalg.norm(e1) # Normalize the direction vectors to become unit vectors 186 e2 = e2/np.linalg.norm(e2) 187 e3 = np.cross(e1, e2) 188 189 elif e3 is not None: 190 e3 = np.array(e3).flatten() 191 e2_tmp = np.cross(e3, e1) # Direction of the third unit vector 192 e3 = np.cross(e1, e2_tmp) 193 e1 = e1/np.linalg.norm(e1) # Normalize the direction vectors to become unit vectors 194 e3 = e3/np.linalg.norm(e3) 195 e2 = np.cross(e3, e1) 196 197 if e2 is None and e3 is None: 198 raise ValueError('Specify either e2 or e3') 199 200 T = np.vstack([e1,e2,e3]) 201 202 return T 203 204 205def gdof_ix_from_nodelabels(all_node_labels, node_labels, dof_ix=[0,1,2]): # from nlfe - not debugged 206 ''' 207 Get global DOF indices from node labels. 208 209 Arguments 210 ------------ 211 all_node_labels : int 212 list or array of all node labels 213 node_labels : int 214 list or array of the node labels to get the indices of 215 dof_ix : int 216 local DOF indices to use 217 218 Returns 219 --------- 220 gdof_ix : int 221 global indices of DOFs of requested nodes 222 ''' 223 if type(node_labels) is not list: 224 node_labels = [node_labels] 225 226 node_ix = [np.where(nl==all_node_labels)[0] for nl in node_labels] 227 gdof_ix = gdof_from_nodedof(node_ix, dof_ix) 228 229 return gdof_ix 230 231 232def gdof_from_nodedof(node_ix, dof_ixs, n_dofs=3, merge=True): 233 r''' 234 Get global DOF from node DOF. 235 236 Arguments 237 ---------- 238 node_ix : int 239 node indices to establish global DOF indices of 240 dof_ixs : int 241 local DOF indices to include 242 n_dofs : 3, optional 243 number of DOFs (all) 244 merge : True, optional 245 whether or not to merge the global DOF indices from all nodes, to a single list 246 247 Returns 248 --------- 249 gdof_ix : int 250 global indices of DOFs of requested nodes 251 252 ''' 253 gdof_ix = [] 254 255 if type(node_ix) is not list: 256 node_ix = [node_ix] 257 258 if type(dof_ixs) is not list: 259 dof_ixs = [dof_ixs]*len(node_ix) 260 261 elif len(dof_ixs) != len(node_ix): 262 dof_ixs = [dof_ixs]*len(node_ix) 263 264 for ix, n in enumerate(node_ix): 265 gdof_ix.append(n*n_dofs + np.array(dof_ixs[ix])) 266 267 if merge: 268 gdof_ix = np.array(gdof_ix).flatten() 269 270 return gdof_ix 271 272def B_to_dofpairs(B, master_val=1): 273 r''' 274 Establish pairs of indices of DOFs to couple from compatibility matrix B. 275 276 Arguments 277 ---------- 278 B : float 279 Lagrange compatiblity \([B]\) matrix 280 master_val : float 281 value used to identify master DOF 282 283 Returns 284 ---------- 285 dof_pairs : int 286 indices of DOFs paired through input \([B]\) matrix 287 288 ''' 289 n_constr, n_dofs = B.shape 290 dof_pairs = [None]*n_constr 291 292 for icon in range(0, n_constr): 293 master = np.where(B[icon,:] == master_val)[0][0] 294 slave = np.where(B[icon,:] == -master_val)[0] 295 if len(slave) != 0: 296 slave = slave[0] 297 else: 298 slave = None 299 300 dof_pairs[icon] = [master, slave] 301 dof_pairs = np.array(dof_pairs).T 302 303 return dof_pairs 304 305 306def dof_pairs_to_Linv(dof_pairs, n_dofs): 307 r''' 308 Establish quasi-inverse of L from dof pairs describing constraints. 309 310 Arguments 311 ---------- 312 dof_pairs : int 313 list of lists of indices of DOFs paired 314 n_dofs : int 315 number of DOFs 316 317 Returns 318 ---------- 319 Linv : int 320 321 Notes 322 ---------- 323 \([L_{inv}]\) is a constructed quantity, conventient as it can construct the constrainted u (unique, free DOFs) from the full as follows: 324 325 $$ 326 \{u_c\} = [L_{inv}] \{u\} 327 $$ 328 329 It is worth noting that it only works for fully hard-constrained systems (no weak connections, or similar). 330 331 ''' 332 333 slave_nodes = dof_pairs[dof_pairs[:,1]!=None, 1] 334 grounded_nodes = dof_pairs[dof_pairs[:,1]==None, 0] 335 336 remove = np.unique(np.hstack([grounded_nodes, slave_nodes])) 337 all_dofs = np.arange(0,n_dofs) 338 primal_dofs = np.setdiff1d(all_dofs, remove) 339 340 Linv = np.zeros([len(primal_dofs), n_dofs]) 341 342 for primal_dof_ix, primal_dof in enumerate(primal_dofs): 343 Linv[primal_dof_ix, primal_dof] = 1 344 345 return Linv 346 347 348def compatibility_matrix(dof_pairs, n_dofs): 349 r''' 350 Establish compatibility matrix from specified pairs of DOFs. 351 352 Arguments 353 ---------- 354 dof_pairs : int 355 list of lists of indices of DOFs paired 356 n_dofs : int 357 number of DOFs in full system; defines the number of rows of \([B]\) 358 359 Returns 360 ---------- 361 B : int 362 numpy array describing compatibility matrix \([B]\) 363 ''' 364 n_constraints = dof_pairs.shape[0] 365 B = np.zeros([n_constraints, n_dofs]) 366 367 for constraint_ix, dof_pair in enumerate(dof_pairs): 368 if dof_pair[1] is None: # connected dof is None --> fixed to ground 369 B[constraint_ix, dof_pair[0]] = 1 370 else: 371 B[constraint_ix, dof_pair[0]] = 1 372 B[constraint_ix, dof_pair[1]] = -1 373 return B 374 375 376def lagrange_constrain(mat, dof_pairs, null=False): 377 r''' 378 Lagrange constrain matrix from specified DOF pairs. 379 380 Arguments 381 ----------------- 382 mat : float 383 matrix to constrain 384 dof_pairs : int 385 list of lists of indices of DOFs paired 386 null : False, optional 387 to create a 0 matrix of the correct size, this can be set to True 388 389 Returns 390 ---------- 391 mat_fixed : float 392 constrained matrix 393 ''' 394 # Lagrange multiplier method - expand matrix 395 n_constraints = len(dof_pairs) 396 n_dofs = mat.shape[0] 397 398 if not null: 399 B = compatibility_matrix(dof_pairs, n_dofs) 400 else: 401 B = np.zeros([n_constraints, n_dofs]) 402 403 O = np.zeros([n_constraints, n_constraints]) 404 mat_fixed = np.vstack([np.hstack([mat, B.T]), np.hstack([B, O])]) 405 406 return mat_fixed 407 408def lagrange_constrain_from_B(mat, B, null=False): 409 r''' 410 Lagrange constrain matrix from specified compatibility matrix. 411 412 Arguments 413 ----------------- 414 mat : float 415 matrix to constrain 416 B : int 417 numpy array describing compatibility matrix \([B]\) 418 null : False, optional 419 to create a 0 matrix of the correct size, this can be set to True 420 421 Returns 422 ---------- 423 mat_fixed : float 424 constrained matrix 425 ''' 426 # Lagrange multiplier method - expand matrix 427 n_constraints = B.shape[0] 428 O = np.zeros([n_constraints, n_constraints]) 429 mat_fixed = np.vstack([np.hstack([mat, B.T]), np.hstack([B, O])]) 430 431 return mat_fixed 432 433 434def expand_vector_with_zeros(vec, n_constraints): 435 r''' 436 Append vector with zeros based on number of constraints (n_constraints). 437 ''' 438 vec_expanded = np.hstack([vec.flatten(), np.zeros(n_constraints)])[np.newaxis, :].T 439 return vec_expanded 440 441 442def blkdiag(mat, n): 443 return np.kron(np.eye(n), mat) 444 445 446def nodes_to_beam_element_matrix(node_labels, first_element_label=1): 447 r''' 448 Establish element matrix definition assuming ordered node labels. 449 450 Arguments 451 ----------- 452 node_labels : int 453 list of integer node labels 454 first_element_label : int 455 first integer used in element matrix 456 457 Returns 458 ----------- 459 element_matrix : int 460 element definition where each row represents an element as [element_label_i, node_label_1, node_label_2] 461 462 ''' 463 n_nodes = len(node_labels) 464 n_elements = n_nodes-1 465 466 element_matrix = np.empty([n_elements, 3]) 467 element_matrix[:, 0] = np.arange(first_element_label,first_element_label+n_elements) 468 element_matrix[:, 1] = node_labels[0:-1] 469 element_matrix[:, 2] = node_labels[1:] 470 471 return element_matrix 472 473 474def get_phys_modes(phi, B, lambd=None, return_as_ix=False): 475 r''' 476 Get physical modes. 477 478 Arguments 479 ----------- 480 phi : float 481 B : int 482 numpy array describing compatibility matrix \([B]\) 483 lambd : float, optional 484 standard value None does not 485 return_as_ix : False, optional 486 whether or not to output as indices instead of eigenvectors and eigenvalues 487 488 Returns 489 ---------- 490 lambd_phys : float 491 physical eigenvalue array 492 phi_phys : float 493 physical eigenvector array (modal transformation matrix) 494 phys_ix : int 495 indices of selected DOFs (if `return_as_ix=True` is passed) 496 ''' 497 n_lagr = B.shape[0] 498 u_ = phi[0:-n_lagr, :] #physical dofs 499 l_ = phi[-n_lagr:,:] #lagr dofs 500 501 # Compatibility at interface 502 mask1 = np.all((B @ u_)==0, axis=0) 503 504 # Equilibrium at interface 505 L = null(B) 506 g = -B.T @ l_ 507 mask2 = np.all((L.T @ g)==0, axis=0) 508 509 phys_ix = np.logical_and(mask1, mask1) 510 511 if return_as_ix: 512 return phys_ix 513 else: 514 if lambd is None: 515 lambd_phys = None 516 else: 517 lambd_phys = lambd[phys_ix] 518 519 phi_phys = u_[:, phys_ix] 520 521 return lambd_phys, phi_phys 522 523 524def ensure_list(list_in, levels=1, increase_only=True): 525 r''' 526 Ensure input variable is list. 527 528 Arguments 529 ---------- 530 list_in : float, int 531 list or scalar 532 levels : int 533 number of levels of list (nest-level) 534 increase_only : True, optional 535 if True, only increase amount of list-wrapping 536 537 Returns 538 ---------- 539 list_out : float, int 540 list 541 ''' 542 543 dimensions = np.array(list_in).ndim 544 545 if type(list_in) is not list: 546 list_out = [list_in] 547 dimensions = dimensions+1 548 else: 549 list_out = list_in * 1 550 551 552 if not increase_only: 553 while dimensions>levels: 554 list_out = list_out[0] 555 dimensions = dimensions - 1 556 557 while dimensions<levels: 558 list_out = [list_out] 559 dimensions = dimensions + 1 560 561 return list_out 562 563 564def feature_matrix(master_dofs, values, slave_dofs=None, ndofs=None, return_small=False): 565 r''' 566 Arguments 567 --------- 568 master_dofs : int 569 list of indices of master DOFs 570 values : float 571 list of amplitudes/values for features 572 slave_dofs : int, optional 573 list of indices of slave DOFs (standard value None, makes grounded connections from master_dofs) 574 ndofs : int, optional 575 number of modes used to construct full matrix (not relevant if `return_small=True`) 576 return_small : False, optional 577 whether or not to return small or full matrix - if `return_small=True`, returns small ndofs-by-ndofs matrix and corresponding indices in global matrix; if `return_small=False` returns big ndofs_global-by-ndofs_global matrix 578 579 Returns 580 --------- 581 mat : float 582 feature matrix 583 dof_ixs : int 584 list of indices of DOFs for relevant feature connectivity (only returned if `return_small=True`) 585 586 ''' 587 588 ndofs_small = len(master_dofs) 589 if type(values) is float or type(values) is int: 590 values = [float(values)]*ndofs_small 591 elif len(values) != len(master_dofs): 592 raise ValueError('Length of master_dofs and values must match') 593 594 if slave_dofs is None: 595 small = np.diag(values) 596 slave_dofs = [] 597 else: 598 if len(slave_dofs) != len(master_dofs): 599 raise ValueError('Length of master_dofs and slave_dofs must match') 600 small = np.diag(np.hstack([values,values])) 601 for ix in range(ndofs_small): 602 small[ix, ix+ndofs_small] = small[ix+ndofs_small, ix] = -values[ix] 603 604 dof_ixs = np.hstack([master_dofs, slave_dofs]).astype(int) 605 606 if return_small: 607 return small, dof_ixs 608 else: 609 if len(set(dof_ixs)) != len(dof_ixs): 610 raise ValueError('Non-unique dof indices are provided') 611 if ndofs is None: 612 ndofs = np.max(dof_ixs)+1 613 614 large = np.zeros([ndofs, ndofs]) 615 large[np.ix_(dof_ixs, dof_ixs)] = small 616 617 return large 618 619def basic_coupled(): 620 return np.array([[1, -1], [-1, 1]]) 621 622def generic_beam_mat(L, yy=0, yz=0, yt=0, zy=0, zz=0, zt=0, ty=0, tz=0, tt=0): 623 mat = np.zeros([12,12]) 624 625 mat[0:6, 0:6] = np.array([ 626 [0, 0, 0, 0, 0, 0 ], 627 [0, 156*yy, 156*yz, 147*yt, -22*L*yz, 22*L*yy ], 628 [0, 156*zy, 156*zz, 147*zt, -22*L*zz, 22*L*zy ], 629 [0, 147*ty, 147*tz, 140*tt, -21*L*tz, 21*L*ty ], 630 [0, -22*L*zy, -22*L*zz, -21*L*zt, 4*L**2*zz, -4*L**2*zy ], 631 [0, 22*L*yy, 22*L*yz, 21*L*yt, -4*L**2*yz, 4*L**2*yy ], 632 ]) 633 634 mat[0:6, 6:12] = np.array([ 635 [0, 0, 0, 0, 0, 0 ], 636 [0, 54*yy, 54*yz, 63*yt, 13*L*yz, -13*L*yy ], 637 [0, 54*zy, 54*zz, 63*zt, 13*L*zz, -13*L*zy ], 638 [0, 63*ty, 63*tz, 70*tt, 14*L*tz, -14*L*ty ], 639 [0, -13*L*zy, -13*L*zz, -14*L*zt, -3*L**2*zz, 3*L**2*zy ], 640 [0, 13*L*yy, 13*L*yz, 14*L*yt, 3*L**2*yz, -3*L**2*yy ], 641 ]) 642 643 mat[6:12, 0:6] = np.array([ 644 [0, 0, 0, 0, 0, 0 ], 645 [0, 54*yy, 54*yz, 63*yt, -13*L*yz, 13*L*yy ], 646 [0, 54*zy, 54*zz, 63*zt, -13*L*zz, 13*L*zy ], 647 [0, 63*ty, 63*tz, 70*tt, -14*L*tz, 14*L*ty ], 648 [0, 13*L*zy, 13*L*zz, 14*L*zt, -3*L**2*zz, 3*L**2*zy ], 649 [0, -13*L*yy, -13*L*yz, -14*L*yt, 3*L**2*yz, -3*L**2*yy ], 650 ]) 651 652 mat[6:12,6:12] = np.array([ 653 [0, 0, 0, 0, 0, 0 ], 654 [0, 156*yy, 156*yz, 147*yt, 22*L*yz, -22*L*yy ], 655 [0, 156*zy, 156*zz, 147*zt, 22*L*zz, -22*L*zy ], 656 [0, 147*ty, 147*tz, 140*tt, 21*L*tz, -21*L*ty ], 657 [0, 22*L*zy, 22*L*zz, 21*L*zt, 4*L**2*zz, -4*L**2*zy ], 658 [0, -22*L*yy, -22*L*yz, -21*L*yt, -4*L**2*yz, 4*L**2*yy ], 659 ]) 660 661 return mat * L / 420 662 663 664def bar_foundation_stiffness(L, kx, ky, kz): #only z and y currently, will be extended! 665 mat = np.vstack([ 666 [kx*1/4, 0, 0, 0, 0, 0, kx*1/4, 0, 0, 0, 0, 0], #ux1 667 [0, ky*1/3, 0, 0, 0, 0, 0, ky*1/6, 0, 0, 0, 0], #uy1 668 [0, 0, kz*1/3, 0, 0, 0, 0, 0, kz*1/6, 0, 0, 0], #uz1 669 670 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rx1 671 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #ry1 672 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rz1 673 674 [kx*1/4, 0, 0, 0, 0, 0, kx*1/4, 0, 0, 0, 0, 0], #ux2 675 [0, ky*1/6, 0, 0, 0, 0, 0, ky*1/3, 0, 0, 0, 0], #uy2 676 [0, 0, kz*1/6, 0, 0, 0, 0, 0, kz*1/3, 0, 0, 0], #uz2 677 678 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rx2 679 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #ry2 680 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) * L #rz2 681 return mat
10def n2d_ix(node_ixs, n_dofs=6): 11 r''' 12 Convert node indices to dof indices. 13 14 Arguments 15 ---------- 16 node_ixs : int 17 list/array of indices of nodes 18 n_dofs : 6, optional 19 number of dofs per node 20 21 Returns 22 ---------- 23 dof_ixs : int 24 array of DOF indices corresponding to input node indices 25 ''' 26 return np.hstack([node_ix*n_dofs+np.arange(n_dofs) for node_ix in node_ixs])
Convert node indices to dof indices.
Arguments
node_ixs : int
list/array of indices of nodes
n_dofs : 6, optional
number of dofs per node
Returns
- dof_ixs (int): array of DOF indices corresponding to input node indices
29def extract_dofs(mat, dof_ix=[0,1,2], n_dofs=6, axis=0): 30 r''' 31 Extract selected dofs from matrix. 32 33 Arguments 34 ------------ 35 mat : float 36 numpy array (matrix) to extract DOFs from 37 dof_ix : int, optional 38 list of DOF indices (local numbering of DOFs in node) to extract 39 n_dofs : 6, optional 40 number of DOFs per node 41 axis : 0, optional 42 axis to pick along 43 44 Returns 45 ---------- 46 mat_sub : float 47 subselection of matrix 48 ''' 49 get_dofs = np.vstack([np.arange(dof, mat.shape[axis],n_dofs) for dof in dof_ix]).T.flatten() 50 return mat.take(get_dofs, axis=axis)
Extract selected dofs from matrix.
Arguments
mat : float numpy array (matrix) to extract DOFs from dof_ix : int, optional list of DOF indices (local numbering of DOFs in node) to extract n_dofs : 6, optional number of DOFs per node axis : 0, optional axis to pick along
Returns
- mat_sub (float): subselection of matrix
52def convert_dofs(dofs_in, node_type='beam3d', sort_output=True): 53 r''' 54 Convert string DOFs to indices. 55 56 Arguments 57 ---------- 58 dofs_in : {'all', 'trans', 'rot'} or int 59 string or list of indices 60 node_type : {'beam2d', 'beam3d'}, optional 61 type of nodes (defining number of DOFs per node) 62 sort_output : True, optional 63 whether or not to sort the DOF indices output 64 65 Returns 66 -------- 67 dofs_out : int 68 array of DOF indices 69 ''' 70 dof_dict = dict(beam2d=dict(all=np.arange(0,3), trans=np.arange(0,2), rot=np.arange(2,3)), 71 beam3d=dict(all=np.arange(0,6), trans=np.arange(0,3), rot=np.arange(3,6))) 72 73 if dofs_in in ['all', 'trans', 'rot']: 74 dofs_out = dof_dict[node_type][dofs_in] 75 else: 76 dofs_out = dofs_in 77 78 if sort_output: 79 dofs_out = np.sort(dofs_out) 80 81 return dofs_out
Convert string DOFs to indices.
Arguments
dofs_in : {'all', 'trans', 'rot'} or int string or list of indices node_type : {'beam2d', 'beam3d'}, optional type of nodes (defining number of DOFs per node) sort_output : True, optional whether or not to sort the DOF indices output
Returns
- dofs_out (int): array of DOF indices
84def convert_dofs_list(dofs_list_in, n_nodes, node_type='beam3d', sort_output=True): 85 r''' 86 Convert DOFs list to correct format. 87 88 Arguments 89 ----------- 90 dofs_list_in : str,int 91 list of strings and indices corresponding to DOFs to 92 n_nodes : int 93 number of nodes the DOFs are related to 94 node_type : {'beam3d', 'beam2d'} 95 type of nodes (defining number of DOFs per node) 96 sort_output : True 97 whether or not to sort output 98 99 Returns 100 ----------- 101 dofs_list_out : int 102 array of DOF indices 103 ''' 104 105 contains_strings = np.any([type(d) is str for d in dofs_list_in]) 106 107 if type(dofs_list_in) is not list: # single input (all, rot or trans) 108 dofs_list_out = [dofs_list_in]*n_nodes 109 elif ~contains_strings and (len(dofs_list_in)<=6) and (np.max(dofs_list_in)<6): # 110 dofs_list_out = [dofs_list_in]*n_nodes 111 elif len(dofs_list_in)!=n_nodes: 112 raise TypeError('Wrong input format of "dofs"') 113 else: 114 dofs_list_out = dofs_list_in 115 116 for ix, dofs_list_out_i in enumerate(dofs_list_out): 117 dofs_list_out[ix] = convert_dofs(dofs_list_out_i, node_type=node_type, sort_output=sort_output) 118 119 return dofs_list_out
Convert DOFs list to correct format.
Arguments
dofs_list_in : str,int list of strings and indices corresponding to DOFs to n_nodes : int number of nodes the DOFs are related to node_type : {'beam3d', 'beam2d'} type of nodes (defining number of DOFs per node) sort_output : True whether or not to sort output
Returns
- dofs_list_out (int): array of DOF indices
122def transform_unit(e1, e2=None, e3=None, warnings=False): 123 r''' 124 Establish transformation matrix from e1 and temporary e2 or e3 vectors. 125 126 Arguments 127 ----------- 128 e1 : float 129 unit vector describing element longitudinal direction 130 e2 : float, optional 131 temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate y-direction) 132 e3 : float, optional 133 temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate z-direction) 134 if both e2 and e3 are different from None, e2 is used (e3 disregarded) 135 136 Returns 137 ----------- 138 T : float 139 transformation matrix 140 141 142 Notes 143 ----------- 144 Input vectors \(\{e_1\}\) and \(\{e_{2p}\}\) are used to establish a third unit vector, as the cross product of the two, as follows: 145 146 $$ 147 \{e_3\} = \{e_1\} \\times \{e_{2p}\} 148 $$ 149 150 Then, a final vector is created based on a second cross product: 151 $$ 152 \{e_2\} = \{e_3\} \\times \{e_1\} 153 $$ 154 155 Finally, the transformation matrix is established as follows: 156 157 $$ 158 [T] = \\begin{bmatrix} 159 \{e_1\}^T \\\\ 160 \{e_2\}^T \\\\ 161 \{e_3\}^T 162 \\end{bmatrix} 163 $$ 164 165 where the unit vectors established are the rows of the transformation matrix. 166 167 ''' 168 169 e1 = np.array(e1).flatten() 170 171 if (e2 is not None) and (e3 is not None): 172 e3 = None 173 174 if e2 is not None: 175 e2 = np.array(e2).flatten() 176 e3_tmp = np.cross(e1, e2) # Direction of the third unit vector 177 if np.all(e3==0): 178 if warnings: 179 print('Warning: e1 and e2 identical. Check orientations carefully!') 180 181 e2 = e2 + 0.1 182 e3 = np.cross(e1, e2) 183 184 e2 = np.cross(e3_tmp, e1) # Direction of the second unit vector 185 186 e1 = e1/np.linalg.norm(e1) # Normalize the direction vectors to become unit vectors 187 e2 = e2/np.linalg.norm(e2) 188 e3 = np.cross(e1, e2) 189 190 elif e3 is not None: 191 e3 = np.array(e3).flatten() 192 e2_tmp = np.cross(e3, e1) # Direction of the third unit vector 193 e3 = np.cross(e1, e2_tmp) 194 e1 = e1/np.linalg.norm(e1) # Normalize the direction vectors to become unit vectors 195 e3 = e3/np.linalg.norm(e3) 196 e2 = np.cross(e3, e1) 197 198 if e2 is None and e3 is None: 199 raise ValueError('Specify either e2 or e3') 200 201 T = np.vstack([e1,e2,e3]) 202 203 return T
Establish transformation matrix from e1 and temporary e2 or e3 vectors.
Arguments
e1 : float unit vector describing element longitudinal direction e2 : float, optional temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate y-direction) e3 : float, optional temporary unit vector describing a chosen vector that's perpendicular to the longitudinal direction (approximate z-direction) if both e2 and e3 are different from None, e2 is used (e3 disregarded)
Returns
- T (float): transformation matrix
Notes
Input vectors ({e_1}) and ({e_{2p}}) are used to establish a third unit vector, as the cross product of the two, as follows:
$$ {e_3} = {e_1} \times {e_{2p}} $$
Then, a final vector is created based on a second cross product: $$ {e_2} = {e_3} \times {e_1} $$
Finally, the transformation matrix is established as follows:
$$ [T] = \begin{bmatrix} {e_1}^T \\ {e_2}^T \\ {e_3}^T \end{bmatrix} $$
where the unit vectors established are the rows of the transformation matrix.
206def gdof_ix_from_nodelabels(all_node_labels, node_labels, dof_ix=[0,1,2]): # from nlfe - not debugged 207 ''' 208 Get global DOF indices from node labels. 209 210 Arguments 211 ------------ 212 all_node_labels : int 213 list or array of all node labels 214 node_labels : int 215 list or array of the node labels to get the indices of 216 dof_ix : int 217 local DOF indices to use 218 219 Returns 220 --------- 221 gdof_ix : int 222 global indices of DOFs of requested nodes 223 ''' 224 if type(node_labels) is not list: 225 node_labels = [node_labels] 226 227 node_ix = [np.where(nl==all_node_labels)[0] for nl in node_labels] 228 gdof_ix = gdof_from_nodedof(node_ix, dof_ix) 229 230 return gdof_ix
Get global DOF indices from node labels.
Arguments
all_node_labels : int list or array of all node labels node_labels : int list or array of the node labels to get the indices of dof_ix : int local DOF indices to use
Returns
- gdof_ix (int): global indices of DOFs of requested nodes
233def gdof_from_nodedof(node_ix, dof_ixs, n_dofs=3, merge=True): 234 r''' 235 Get global DOF from node DOF. 236 237 Arguments 238 ---------- 239 node_ix : int 240 node indices to establish global DOF indices of 241 dof_ixs : int 242 local DOF indices to include 243 n_dofs : 3, optional 244 number of DOFs (all) 245 merge : True, optional 246 whether or not to merge the global DOF indices from all nodes, to a single list 247 248 Returns 249 --------- 250 gdof_ix : int 251 global indices of DOFs of requested nodes 252 253 ''' 254 gdof_ix = [] 255 256 if type(node_ix) is not list: 257 node_ix = [node_ix] 258 259 if type(dof_ixs) is not list: 260 dof_ixs = [dof_ixs]*len(node_ix) 261 262 elif len(dof_ixs) != len(node_ix): 263 dof_ixs = [dof_ixs]*len(node_ix) 264 265 for ix, n in enumerate(node_ix): 266 gdof_ix.append(n*n_dofs + np.array(dof_ixs[ix])) 267 268 if merge: 269 gdof_ix = np.array(gdof_ix).flatten() 270 271 return gdof_ix
Get global DOF from node DOF.
Arguments
node_ix : int node indices to establish global DOF indices of dof_ixs : int local DOF indices to include n_dofs : 3, optional number of DOFs (all) merge : True, optional whether or not to merge the global DOF indices from all nodes, to a single list
Returns
- gdof_ix (int): global indices of DOFs of requested nodes
273def B_to_dofpairs(B, master_val=1): 274 r''' 275 Establish pairs of indices of DOFs to couple from compatibility matrix B. 276 277 Arguments 278 ---------- 279 B : float 280 Lagrange compatiblity \([B]\) matrix 281 master_val : float 282 value used to identify master DOF 283 284 Returns 285 ---------- 286 dof_pairs : int 287 indices of DOFs paired through input \([B]\) matrix 288 289 ''' 290 n_constr, n_dofs = B.shape 291 dof_pairs = [None]*n_constr 292 293 for icon in range(0, n_constr): 294 master = np.where(B[icon,:] == master_val)[0][0] 295 slave = np.where(B[icon,:] == -master_val)[0] 296 if len(slave) != 0: 297 slave = slave[0] 298 else: 299 slave = None 300 301 dof_pairs[icon] = [master, slave] 302 dof_pairs = np.array(dof_pairs).T 303 304 return dof_pairs
Establish pairs of indices of DOFs to couple from compatibility matrix B.
Arguments
B : float
Lagrange compatiblity ([B]) matrix
master_val : float
value used to identify master DOF
Returns
- dof_pairs (int): indices of DOFs paired through input ([B]) matrix
307def dof_pairs_to_Linv(dof_pairs, n_dofs): 308 r''' 309 Establish quasi-inverse of L from dof pairs describing constraints. 310 311 Arguments 312 ---------- 313 dof_pairs : int 314 list of lists of indices of DOFs paired 315 n_dofs : int 316 number of DOFs 317 318 Returns 319 ---------- 320 Linv : int 321 322 Notes 323 ---------- 324 \([L_{inv}]\) is a constructed quantity, conventient as it can construct the constrainted u (unique, free DOFs) from the full as follows: 325 326 $$ 327 \{u_c\} = [L_{inv}] \{u\} 328 $$ 329 330 It is worth noting that it only works for fully hard-constrained systems (no weak connections, or similar). 331 332 ''' 333 334 slave_nodes = dof_pairs[dof_pairs[:,1]!=None, 1] 335 grounded_nodes = dof_pairs[dof_pairs[:,1]==None, 0] 336 337 remove = np.unique(np.hstack([grounded_nodes, slave_nodes])) 338 all_dofs = np.arange(0,n_dofs) 339 primal_dofs = np.setdiff1d(all_dofs, remove) 340 341 Linv = np.zeros([len(primal_dofs), n_dofs]) 342 343 for primal_dof_ix, primal_dof in enumerate(primal_dofs): 344 Linv[primal_dof_ix, primal_dof] = 1 345 346 return Linv
Establish quasi-inverse of L from dof pairs describing constraints.
Arguments
dof_pairs : int
list of lists of indices of DOFs paired
n_dofs : int
number of DOFs
Returns
- Linv (int):
Notes
([L_{inv}]) is a constructed quantity, conventient as it can construct the constrainted u (unique, free DOFs) from the full as follows:
$$ {u_c} = [L_{inv}] {u} $$
It is worth noting that it only works for fully hard-constrained systems (no weak connections, or similar).
349def compatibility_matrix(dof_pairs, n_dofs): 350 r''' 351 Establish compatibility matrix from specified pairs of DOFs. 352 353 Arguments 354 ---------- 355 dof_pairs : int 356 list of lists of indices of DOFs paired 357 n_dofs : int 358 number of DOFs in full system; defines the number of rows of \([B]\) 359 360 Returns 361 ---------- 362 B : int 363 numpy array describing compatibility matrix \([B]\) 364 ''' 365 n_constraints = dof_pairs.shape[0] 366 B = np.zeros([n_constraints, n_dofs]) 367 368 for constraint_ix, dof_pair in enumerate(dof_pairs): 369 if dof_pair[1] is None: # connected dof is None --> fixed to ground 370 B[constraint_ix, dof_pair[0]] = 1 371 else: 372 B[constraint_ix, dof_pair[0]] = 1 373 B[constraint_ix, dof_pair[1]] = -1 374 return B
Establish compatibility matrix from specified pairs of DOFs.
Arguments
dof_pairs : int
list of lists of indices of DOFs paired
n_dofs : int
number of DOFs in full system; defines the number of rows of ([B])
Returns
- B (int): numpy array describing compatibility matrix ([B])
377def lagrange_constrain(mat, dof_pairs, null=False): 378 r''' 379 Lagrange constrain matrix from specified DOF pairs. 380 381 Arguments 382 ----------------- 383 mat : float 384 matrix to constrain 385 dof_pairs : int 386 list of lists of indices of DOFs paired 387 null : False, optional 388 to create a 0 matrix of the correct size, this can be set to True 389 390 Returns 391 ---------- 392 mat_fixed : float 393 constrained matrix 394 ''' 395 # Lagrange multiplier method - expand matrix 396 n_constraints = len(dof_pairs) 397 n_dofs = mat.shape[0] 398 399 if not null: 400 B = compatibility_matrix(dof_pairs, n_dofs) 401 else: 402 B = np.zeros([n_constraints, n_dofs]) 403 404 O = np.zeros([n_constraints, n_constraints]) 405 mat_fixed = np.vstack([np.hstack([mat, B.T]), np.hstack([B, O])]) 406 407 return mat_fixed
Lagrange constrain matrix from specified DOF pairs.
Arguments
mat : float matrix to constrain dof_pairs : int list of lists of indices of DOFs paired null : False, optional to create a 0 matrix of the correct size, this can be set to True
Returns
- mat_fixed (float): constrained matrix
409def lagrange_constrain_from_B(mat, B, null=False): 410 r''' 411 Lagrange constrain matrix from specified compatibility matrix. 412 413 Arguments 414 ----------------- 415 mat : float 416 matrix to constrain 417 B : int 418 numpy array describing compatibility matrix \([B]\) 419 null : False, optional 420 to create a 0 matrix of the correct size, this can be set to True 421 422 Returns 423 ---------- 424 mat_fixed : float 425 constrained matrix 426 ''' 427 # Lagrange multiplier method - expand matrix 428 n_constraints = B.shape[0] 429 O = np.zeros([n_constraints, n_constraints]) 430 mat_fixed = np.vstack([np.hstack([mat, B.T]), np.hstack([B, O])]) 431 432 return mat_fixed
Lagrange constrain matrix from specified compatibility matrix.
Arguments
mat : float matrix to constrain B : int numpy array describing compatibility matrix ([B]) null : False, optional to create a 0 matrix of the correct size, this can be set to True
Returns
- mat_fixed (float): constrained matrix
435def expand_vector_with_zeros(vec, n_constraints): 436 r''' 437 Append vector with zeros based on number of constraints (n_constraints). 438 ''' 439 vec_expanded = np.hstack([vec.flatten(), np.zeros(n_constraints)])[np.newaxis, :].T 440 return vec_expanded
Append vector with zeros based on number of constraints (n_constraints).
447def nodes_to_beam_element_matrix(node_labels, first_element_label=1): 448 r''' 449 Establish element matrix definition assuming ordered node labels. 450 451 Arguments 452 ----------- 453 node_labels : int 454 list of integer node labels 455 first_element_label : int 456 first integer used in element matrix 457 458 Returns 459 ----------- 460 element_matrix : int 461 element definition where each row represents an element as [element_label_i, node_label_1, node_label_2] 462 463 ''' 464 n_nodes = len(node_labels) 465 n_elements = n_nodes-1 466 467 element_matrix = np.empty([n_elements, 3]) 468 element_matrix[:, 0] = np.arange(first_element_label,first_element_label+n_elements) 469 element_matrix[:, 1] = node_labels[0:-1] 470 element_matrix[:, 2] = node_labels[1:] 471 472 return element_matrix
Establish element matrix definition assuming ordered node labels.
Arguments
node_labels : int list of integer node labels first_element_label : int first integer used in element matrix
Returns
- element_matrix (int): element definition where each row represents an element as [element_label_i, node_label_1, node_label_2]
475def get_phys_modes(phi, B, lambd=None, return_as_ix=False): 476 r''' 477 Get physical modes. 478 479 Arguments 480 ----------- 481 phi : float 482 B : int 483 numpy array describing compatibility matrix \([B]\) 484 lambd : float, optional 485 standard value None does not 486 return_as_ix : False, optional 487 whether or not to output as indices instead of eigenvectors and eigenvalues 488 489 Returns 490 ---------- 491 lambd_phys : float 492 physical eigenvalue array 493 phi_phys : float 494 physical eigenvector array (modal transformation matrix) 495 phys_ix : int 496 indices of selected DOFs (if `return_as_ix=True` is passed) 497 ''' 498 n_lagr = B.shape[0] 499 u_ = phi[0:-n_lagr, :] #physical dofs 500 l_ = phi[-n_lagr:,:] #lagr dofs 501 502 # Compatibility at interface 503 mask1 = np.all((B @ u_)==0, axis=0) 504 505 # Equilibrium at interface 506 L = null(B) 507 g = -B.T @ l_ 508 mask2 = np.all((L.T @ g)==0, axis=0) 509 510 phys_ix = np.logical_and(mask1, mask1) 511 512 if return_as_ix: 513 return phys_ix 514 else: 515 if lambd is None: 516 lambd_phys = None 517 else: 518 lambd_phys = lambd[phys_ix] 519 520 phi_phys = u_[:, phys_ix] 521 522 return lambd_phys, phi_phys
Get physical modes.
Arguments
phi : float B : int numpy array describing compatibility matrix ([B]) lambd : float, optional standard value None does not return_as_ix : False, optional whether or not to output as indices instead of eigenvectors and eigenvalues
Returns
- lambd_phys (float): physical eigenvalue array
- phi_phys (float): physical eigenvector array (modal transformation matrix)
- phys_ix (int):
indices of selected DOFs (if
return_as_ix=True
is passed)
525def ensure_list(list_in, levels=1, increase_only=True): 526 r''' 527 Ensure input variable is list. 528 529 Arguments 530 ---------- 531 list_in : float, int 532 list or scalar 533 levels : int 534 number of levels of list (nest-level) 535 increase_only : True, optional 536 if True, only increase amount of list-wrapping 537 538 Returns 539 ---------- 540 list_out : float, int 541 list 542 ''' 543 544 dimensions = np.array(list_in).ndim 545 546 if type(list_in) is not list: 547 list_out = [list_in] 548 dimensions = dimensions+1 549 else: 550 list_out = list_in * 1 551 552 553 if not increase_only: 554 while dimensions>levels: 555 list_out = list_out[0] 556 dimensions = dimensions - 1 557 558 while dimensions<levels: 559 list_out = [list_out] 560 dimensions = dimensions + 1 561 562 return list_out
Ensure input variable is list.
Arguments
list_in : float, int list or scalar levels : int number of levels of list (nest-level) increase_only : True, optional if True, only increase amount of list-wrapping
Returns
- list_out (float, int): list
565def feature_matrix(master_dofs, values, slave_dofs=None, ndofs=None, return_small=False): 566 r''' 567 Arguments 568 --------- 569 master_dofs : int 570 list of indices of master DOFs 571 values : float 572 list of amplitudes/values for features 573 slave_dofs : int, optional 574 list of indices of slave DOFs (standard value None, makes grounded connections from master_dofs) 575 ndofs : int, optional 576 number of modes used to construct full matrix (not relevant if `return_small=True`) 577 return_small : False, optional 578 whether or not to return small or full matrix - if `return_small=True`, returns small ndofs-by-ndofs matrix and corresponding indices in global matrix; if `return_small=False` returns big ndofs_global-by-ndofs_global matrix 579 580 Returns 581 --------- 582 mat : float 583 feature matrix 584 dof_ixs : int 585 list of indices of DOFs for relevant feature connectivity (only returned if `return_small=True`) 586 587 ''' 588 589 ndofs_small = len(master_dofs) 590 if type(values) is float or type(values) is int: 591 values = [float(values)]*ndofs_small 592 elif len(values) != len(master_dofs): 593 raise ValueError('Length of master_dofs and values must match') 594 595 if slave_dofs is None: 596 small = np.diag(values) 597 slave_dofs = [] 598 else: 599 if len(slave_dofs) != len(master_dofs): 600 raise ValueError('Length of master_dofs and slave_dofs must match') 601 small = np.diag(np.hstack([values,values])) 602 for ix in range(ndofs_small): 603 small[ix, ix+ndofs_small] = small[ix+ndofs_small, ix] = -values[ix] 604 605 dof_ixs = np.hstack([master_dofs, slave_dofs]).astype(int) 606 607 if return_small: 608 return small, dof_ixs 609 else: 610 if len(set(dof_ixs)) != len(dof_ixs): 611 raise ValueError('Non-unique dof indices are provided') 612 if ndofs is None: 613 ndofs = np.max(dof_ixs)+1 614 615 large = np.zeros([ndofs, ndofs]) 616 large[np.ix_(dof_ixs, dof_ixs)] = small 617 618 return large
Arguments
master_dofs : int
list of indices of master DOFs
values : float
list of amplitudes/values for features
slave_dofs : int, optional
list of indices of slave DOFs (standard value None, makes grounded connections from master_dofs)
ndofs : int, optional
number of modes used to construct full matrix (not relevant if return_small=True
)
return_small : False, optional
whether or not to return small or full matrix - if return_small=True
, returns small ndofs-by-ndofs matrix and corresponding indices in global matrix; if return_small=False
returns big ndofs_global-by-ndofs_global matrix
Returns
- mat (float): feature matrix
- dof_ixs (int):
list of indices of DOFs for relevant feature connectivity (only returned if
return_small=True
)
623def generic_beam_mat(L, yy=0, yz=0, yt=0, zy=0, zz=0, zt=0, ty=0, tz=0, tt=0): 624 mat = np.zeros([12,12]) 625 626 mat[0:6, 0:6] = np.array([ 627 [0, 0, 0, 0, 0, 0 ], 628 [0, 156*yy, 156*yz, 147*yt, -22*L*yz, 22*L*yy ], 629 [0, 156*zy, 156*zz, 147*zt, -22*L*zz, 22*L*zy ], 630 [0, 147*ty, 147*tz, 140*tt, -21*L*tz, 21*L*ty ], 631 [0, -22*L*zy, -22*L*zz, -21*L*zt, 4*L**2*zz, -4*L**2*zy ], 632 [0, 22*L*yy, 22*L*yz, 21*L*yt, -4*L**2*yz, 4*L**2*yy ], 633 ]) 634 635 mat[0:6, 6:12] = np.array([ 636 [0, 0, 0, 0, 0, 0 ], 637 [0, 54*yy, 54*yz, 63*yt, 13*L*yz, -13*L*yy ], 638 [0, 54*zy, 54*zz, 63*zt, 13*L*zz, -13*L*zy ], 639 [0, 63*ty, 63*tz, 70*tt, 14*L*tz, -14*L*ty ], 640 [0, -13*L*zy, -13*L*zz, -14*L*zt, -3*L**2*zz, 3*L**2*zy ], 641 [0, 13*L*yy, 13*L*yz, 14*L*yt, 3*L**2*yz, -3*L**2*yy ], 642 ]) 643 644 mat[6:12, 0:6] = np.array([ 645 [0, 0, 0, 0, 0, 0 ], 646 [0, 54*yy, 54*yz, 63*yt, -13*L*yz, 13*L*yy ], 647 [0, 54*zy, 54*zz, 63*zt, -13*L*zz, 13*L*zy ], 648 [0, 63*ty, 63*tz, 70*tt, -14*L*tz, 14*L*ty ], 649 [0, 13*L*zy, 13*L*zz, 14*L*zt, -3*L**2*zz, 3*L**2*zy ], 650 [0, -13*L*yy, -13*L*yz, -14*L*yt, 3*L**2*yz, -3*L**2*yy ], 651 ]) 652 653 mat[6:12,6:12] = np.array([ 654 [0, 0, 0, 0, 0, 0 ], 655 [0, 156*yy, 156*yz, 147*yt, 22*L*yz, -22*L*yy ], 656 [0, 156*zy, 156*zz, 147*zt, 22*L*zz, -22*L*zy ], 657 [0, 147*ty, 147*tz, 140*tt, 21*L*tz, -21*L*ty ], 658 [0, 22*L*zy, 22*L*zz, 21*L*zt, 4*L**2*zz, -4*L**2*zy ], 659 [0, -22*L*yy, -22*L*yz, -21*L*yt, -4*L**2*yz, 4*L**2*yy ], 660 ]) 661 662 return mat * L / 420
665def bar_foundation_stiffness(L, kx, ky, kz): #only z and y currently, will be extended! 666 mat = np.vstack([ 667 [kx*1/4, 0, 0, 0, 0, 0, kx*1/4, 0, 0, 0, 0, 0], #ux1 668 [0, ky*1/3, 0, 0, 0, 0, 0, ky*1/6, 0, 0, 0, 0], #uy1 669 [0, 0, kz*1/3, 0, 0, 0, 0, 0, kz*1/6, 0, 0, 0], #uz1 670 671 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rx1 672 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #ry1 673 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rz1 674 675 [kx*1/4, 0, 0, 0, 0, 0, kx*1/4, 0, 0, 0, 0, 0], #ux2 676 [0, ky*1/6, 0, 0, 0, 0, 0, ky*1/3, 0, 0, 0, 0], #uy2 677 [0, 0, kz*1/6, 0, 0, 0, 0, 0, kz*1/3, 0, 0, 0], #uz2 678 679 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #rx2 680 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], #ry2 681 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) * L #rz2 682 return mat