Page 4 - Demo

Planar Truss Example for Comrel Add-on RCP Consult, 2019-2024 Page 4T = nfe*[None] # transformation matrixfor e in range(nfe): idx = np.r_[dof[IEN[e,0],:], dof[IEN[e,1],:]] # extract the dofs of the element c = np.cos(theta[e]) # cosine inclination angle s = np.sin(theta[e]) # sinus inclination angle T[e] = np.array([[ c, s, 0, 0], # coordinate transformation matrix [-s, c, 0, 0], [ 0, 0, c, s], [ 0, 0, -s, c]]) Kloc = k[e]*np.array([[ 1, 0, -1, 0], # local stiffness matrix [ 0, 0, 0, 0], [-1, 0, 1, 0], [ 0, 0, 0, 0]]) K[np.ix_(idx,idx)] += T[e].T.dot(Kloc).dot(T[e]) # add to K global# Solve systems of equations# f = vector of equivalent nodal forces# q = vector of equilibrium nodal forces# a = displacements#| qd | | Kcc Kcd || ac | | fd |#| | = | || | - | |#| qc | | Kdc Kdd || ad | | fc |Kcc = K[cc,:][:,cc]Kcd = K[cc,:][:,dd]Kdc = K[dd,:][:,cc]Kdd = K[dd,:][:,dd]ac = np.array([[0], [0], [0]]) # displacements for fixed/support nodesad = np.linalg.solve(Kdd, fc - Kdc.dot(ac)) # solutionqd = Kcc.dot(ac) + Kcd.dot(ad)## Assemble vectors of displacements (a) and forces (q)a = np.zeros((ndof,1))q = np.zeros((ndof,1))a[cc] = acq[cc] = qda[dd] = ad # q[dd] = qc = 0## Compute axial loadsN = np.zeros(nfe)for e in range(nfe): idx = np.r_[dof[IEN[e,0],:], dof[IEN[e,1],:]] N[e] = k[e]*np.array([-1, 0, 1, 0]).dot(T[e]).dot(a[idx])## Output vertical deflection at bottom center (N04)uout = a[7]## Visualisationif StrurelPlot : # begin of plot block if StrurelMode == 0 : tit='Planar Truss - Deterministic Solution' # 0 - initial deterministic solution elif StrurelMode == 2 : tit='Planar Truss - Stochastic Solution' # 2 - final stochastic solution if StrurelMode == 0 or StrurelMode == 2 : # Get screen resolution and DPI user32 = ctypes.windll.user32 user32.SetProcessDPIAware() screensize = user32.GetSystemMetrics(0), user32.GetSystemMetrics(1) sysDpi = user32.GetDpiForSystem() px = 1./sysDpi # X and Y coordinates in m XY = np.array([[0, 0], [4, 0], [8, 0], [12, 0], [16, 0], [20, 0], [24, 0], [22, 2], [18, 2], [14, 2], [10, 2], [6, 2], [2, 2]]) ## Create 3 plots for nf in range(3): # Create figure with unique name for each plot plt.figure('%s [%d]' % (tit,nf+1),figsize=[screensize[0]/2*px, screensize[1]/2*px]) if nf == 0: # Set a main subtitle of plot plt.text(12.0, 8.0, 'Initial and deformed states',fontsize=16,ha='center'); bb=np.reshape(a, (nnp, ned)) # deflection vector DEF = XY + bb*2 # deformed coordinates with scale factor 2. fd = np.zeros(ndof) fd[15::2] = np.array(P) fd=np.reshape(fd, (nnp, ned))*0.000025 # forces with scale factor 0.000025 # Plot initial and deformed states of truss with loads for e in range(nfe): plt.plot(XY[IEN[e,:],0],XY[IEN[e,:],1],\# blue color for initial state plt.plot(DEF[IEN[e,:],0],DEF[IEN[e,:],1],\# red color for deformed state x=(XY[IEN[e,0],0]+XY[IEN[e,1],0])/2 y=(XY[IEN[e,0],1]+XY[IEN[e,1],1])/2 plt.text(x, y, 'R'+str(e+1),fontsize=8) # rod number in the middle of rod length for e in range(nnp): plt.text(XY[e,0],XY[e,1],'N'+str(e+1),fontsize=8) # nodal numbers plt.text(DEF[3,0], DEF[3,1]-0.2, 'U = %g m' % bb[3,1],fontsize=8) # value of vertical deflection in node N04 n=0 for e in range(nnp): if(fd[e,1] > 0.1) : n=n+1; plt.arrow(XY[e,0]+fd[e,0],XY[e,1]+fd[e,1],-fd[e,0],-fd[e,1],length_includes_head=True,head_width=0.2,head_length=0.1) # single vector plt.text(XY[e,0]+0.1+fd[e,0],XY[e,1]+fd[e,1],'P%d=%g N' % (n, fd[e,1]*40000),fontsize=8) # force value elif nf == 1: # Set a main subtitle of plot