import numpy as np from scipy.sparse.linalg import eigs from computeEquationNumber import computeEquationNumber from assembleGlobalStiffnessMatrix import assembleInitialGlobalMaterialStiffnessMatrix from plotter import plotDeformedStructure2d, plotForceDisplacement, plotInternalForces2dBeam from elemStateDetermin import elemStateDetermin from DOFNumberer import directConstraintEnformcement from fiberSection import meshFiber_rectangle from computeQuadrature import computeQuadrature #Dimension and degrees of freedom of the problem nDim=2 nDofPerNode=3 #nDof per node nNodesPerElement=2 #Frame geometry h=7*1e3 l1=6*1e3 l2=12*1e3 #Node coordinates pts_leftCol = [] pts_middleCol = [] pts_rightCol = [] pts_leftBeam = [] pts_rightBeam = [] nEl_column=5 nEl_beam=5 for i in range(nEl_column+1): y = i * (h / nEl_column ) pts_leftCol.append((0, y)) for i in range(nEl_column+1): y=h-i * (h / nEl_column ) pts_middleCol.append((l1, y)) for i in range(nEl_column+1): y=h-i * (h / nEl_column ) pts_rightCol.append((l1+l2, y)) for i in range(nEl_beam+1): x = i * (l1 / nEl_beam) pts_leftBeam.append((x, h)) for i in range(nEl_beam+1): x = l1 + i * (l2 / nEl_beam) pts_rightBeam.append((x, h)) coordinates_list = pts_leftCol + pts_leftBeam + pts_middleCol + pts_rightBeam + pts_rightCol seen = set() unique_coordinates = [] for coord in coordinates_list: if coord not in seen: unique_coordinates.append(coord) seen.add(coord) coordinates = np.array(unique_coordinates) nNodes=coordinates.shape[0] nDofTot=nNodes*nDofPerNode # Connectivity coord_to_index = {tuple(coord): idx for idx, coord in enumerate(coordinates)} # Create a dictionary to map each coordinate to its index connectivity = [] for i in range(nEl_column): connectivity.append([coord_to_index[pts_leftCol[i]], coord_to_index[pts_leftCol[i + 1]]]) for i in range(nEl_beam): connectivity.append([coord_to_index[pts_leftBeam[i]], coord_to_index[pts_leftBeam[i + 1]]]) for i in range(nEl_column): connectivity.append([coord_to_index[pts_middleCol[i]], coord_to_index[pts_middleCol[i + 1]]]) for i in range(nEl_beam): connectivity.append([coord_to_index[pts_rightBeam[i]], coord_to_index[pts_rightBeam[i + 1]]]) for i in range(nEl_column): connectivity.append([coord_to_index[pts_rightCol[i]], coord_to_index[pts_rightCol[i + 1]]]) connectivity = np.array(connectivity) # print(connectivity) nElements=connectivity.shape[0] #Member properties E=200000 #Young's modulus fy=355 #Yield stress alphaS=0.03 #hardening ratio materialInput=[E,fy,alphaS] bc=300 hc=300 bb=300 hb=700 # Define fiber section nFibersInRow=1 nFibersInCol=10 coordinates_fibers_columns,area_fibers_columns = meshFiber_rectangle(hc,bc,nFibersInRow,nFibersInCol) coordinates_fibers_beams,area_fibers_beams = meshFiber_rectangle(hb,bb,nFibersInRow,nFibersInCol) # Define quadrature rule for numerical integration integrationRule='Gauss-Lobatto' numIntegrPoints=5 integrPoints,integrWeights=computeQuadrature(integrationRule,numIntegrPoints) # Geometric transformation local <-> basic reference frame geoTransfCols='linear' # 'linear' or 'corotational' geoTransfBeams='linear' # Elements properties AllElement_data = [] #contains the properties of each element #2dTruss: '2dTruss', nodes, coordNode1, coordNode2, E, A, geoTransf #2dElasticBeam: '2dElasticBeam', nodes, coordNode1, coordNode2, E, A, I, geoTransf #spring: 'spring', nodes, globalDof, springInput, springOutput where globalDof=0 for X, 1 for Y, 2 for rotation #2dModElasticBeam: '2dModElasticBeam', nodes, coordNode1, coordNode2, E, A, S22,S23,S32,S33, Ie, #2dDispBeamColumn: '2dDispBeamColumn', nodes, coordNode1, coordNode2,integrPoints,integrWeights, materialInput, coordinates_fibers, area_fibers, geoTransf for i in range(nElements): dx=coordinates[connectivity[i,1],0]-coordinates[connectivity[i,0],0] dy=coordinates[connectivity[i,1],1]-coordinates[connectivity[i,0],1] if dx==0:#column elements AllElement_data.append(['2dDispBeamColumn', connectivity[i], coordinates[connectivity[i,0]],coordinates[connectivity[i,1]],integrPoints,integrWeights, materialInput, coordinates_fibers_columns, area_fibers_columns, geoTransfCols]) else:#beam elements AllElement_data.append(['2dDispBeamColumn', connectivity[i], coordinates[connectivity[i,0]],coordinates[connectivity[i,1]],integrPoints,integrWeights, materialInput, coordinates_fibers_beams, area_fibers_beams, geoTransfBeams]) # EqualDof constraints matrix (master node, slave node and its constrained dofs where globalDof=0 for X, 1 for Y, 2 for rotation) allConstraints = [] # Equation number matrix numEquations_noConstraints=computeEquationNumber(nDofPerNode,AllElement_data) #print(numEquations) # Definition of the boundary conditions alpha=0.1 #Apply direct method for constraint enforcement F_ext_nodesDofs=[[nEl_column,[0,alpha],[1,-0.5]],[nEl_column+nEl_beam,[1,-1]],[2*nEl_column+2*nEl_beam,[1,-1]]] #format: [node, [dof1, value1], [dof2, value2], ...] fixedBC_nodesDofs=[[0,0,1,2],[2*nEl_column+nEl_beam,0,1,2],[2*nEl_column+3*nEl_beam,0,1,2]] #format: [node, dof1, dof2, ...] where 0 for X, 1 for Y, 2 for rotation numEquations, F_ext, freeDofs, fixedDofs=directConstraintEnformcement(nDofPerNode,numEquations_noConstraints, F_ext_nodesDofs,fixedBC_nodesDofs,allConstraints) nDofTot=max(max(sublist) for sublist in numEquations)+1 # Assemble initial global material stiffness matrix KMaterial_global=assembleInitialGlobalMaterialStiffnessMatrix(nDofTot,numEquations,AllElement_data) KTot_global=KMaterial_global # Initialize variables v_Array=[] F_int_Array=[] lambda_Vector=[] QLocal_elem_list=[] F_int=np.zeros(nDofTot) lambdaIntegrator=0. v=np.zeros(nDofTot) F_unb=np.zeros(nDofTot) F_ext_tot=np.zeros(nDofTot) DeltaV_u=np.zeros(nDofTot) DeltaV_f=np.zeros(nDofTot) # # Force or displacement control integrator FIXED STEP SIZE # integrator='forceControl' # lambda_max=4803400*0.5 # nIncrements=500 # DeltaLambdaBar=lambda_max/nIncrements # tol=1e-4 # nIterMax=1000 integrator='displacementControl' uMax=1000 nIncrements=500 DeltaVBar=uMax/nIncrements q_ctrDof=nDofPerNode*(nEl_column+1)-3; # Control DOF tol=1e-2 nIterMax=5000 if integrator=='forceControl': aIntegrator=np.zeros(nDofTot) bIntegrator=1 cIntegrator=np.array([DeltaLambdaBar,0]) elif integrator=='displacementControl': aIntegrator=np.zeros(nDofTot) aIntegrator[q_ctrDof]=1 bIntegrator=0 cIntegrator=np.array([DeltaVBar,0]) ## Solve the system for n in range(nIncrements): converged=False i=1 # Iteration while not converged: DeltaV_u[freeDofs] = np.linalg.solve(KTot_global[np.ix_(freeDofs, freeDofs)], -F_unb[freeDofs]) DeltaV_f[freeDofs] = np.linalg.solve(KTot_global[np.ix_(freeDofs, freeDofs)], F_ext[freeDofs]) if i==1: DeltaLambda=(cIntegrator[0]-np.dot(aIntegrator,DeltaV_u))/(np.dot(aIntegrator,DeltaV_f)+bIntegrator) else: DeltaLambda=(cIntegrator[1]-np.dot(aIntegrator,DeltaV_u))/(np.dot(aIntegrator,DeltaV_f)+bIntegrator) lambdaIntegrator += DeltaLambda DeltaV = DeltaV_u + DeltaLambda * DeltaV_f v += DeltaV #Element state determination F_int=np.zeros(nDofTot) KTot_global=np.zeros((nDofTot,nDofTot)) QLocal_elem_step = [] for elemIdx, elem in enumerate(AllElement_data): #Perform element state determination elem, QGlobal_elem, KMaterialGlobal_elem, KGeomGlobal_elem, QLocal_elem = elemStateDetermin(elemIdx,elem,numEquations, v) QLocal_elem_step.append(QLocal_elem.copy()) # Assemble element state determination results dof_elem=np.array(numEquations[elemIdx]).astype(int) F_int[dof_elem]+=QGlobal_elem KTot_global[np.ix_(dof_elem,dof_elem)]+=KMaterialGlobal_elem+KGeomGlobal_elem # Compute unbalanced force vector F_ext_tot+=DeltaLambda*F_ext F_unb=F_int-F_ext_tot # Check convergence # test=np.linalg.norm(F_unb[freeDofs]) test=np.linalg.norm(F_unb[freeDofs])/np.linalg.norm(F_ext_tot[freeDofs]) if test