Dual_E030 v2.py

import numpy as np
import openseespy.opensees as ops
import pandas as pd
from openseespy.opensees import *
import openseespy.postprocessing.ops_vis as opsv
import openseespy.postprocessing.Get_Rendering as opsplt
import matplotlib.pyplot as plt

#=================================================================
#SISTEMA DE UNIDADES
# Unidades Base
m = 1
kg = 1
s = 1
# Otras Unidades
cm = 0.01*m
N = kg*m/s**2
kN = 1000*N
kgf = 9.81*N
Pa = N/m**2
MPa = 10**6*Pa
inch = 2.54*cm
ft = 12*inch
ksi = 6894757.2932*Pa
kip = ksi*inch**2
psi = 6894.76*Pa
# Constantes Físicas
g = 9.81*m/s**2
#=================================================================
#GEOMETRÍA DE ESTRUCTURA
dxy = 5*m #separación de las grillas en X e Y
dx, dy, dz = dxy, dxy, 3*m

nxy = 4 #número de grillas en X e Y
nx, ny = nxy, nxy

nzmin=5 #cantidad de pisos mínima
nzmax=10 #cantidad de pisos máxima
dnz=1 #intervalo de pisos
nnz=int((nzmax-nzmin)/dnz+1) #número de intervalos para nz (cantidad de pisos)
Nz=np.zeros(nnz) #lista de nz (cantidad de pisos)
nzi=nzmin
for i in range(nnz):
    Nz[i]=nzi
    nzi=nzi+dnz

#==================================================================
#GEOMETRÍA DE ELEMENTOS ESTRUCTURALES
acmin=0.25*m #ancho mínimo de columna
acmax=0.80*m #ancho máximo de columna
dac=0.05*m #intervalo de ancho de columna
nac=int((acmax-acmin)/dac+1) #número de intervalos para ac (ancho de columna)
ac=np.zeros(nac) #lista de ac (ancho de columna)
aci=acmin
for i in range(nac):
    ac[i]=aci
    aci=aci+dac
#print(ac)

bvmin=0.25*m #ancho mínimo de viga
bvmax=0.40*m #ancho máximo de viga
dbv=0.05*m #intervalo de nacho de viga
nbv=int((bvmax-bvmin)/dbv+1) #número de intervalos para bv (ancho de viga)
bv=np.zeros(nbv) #lista de bv (ancho de viga)
bvi=bvmin
for i in range(nbv):
    bv[i]=bvi
    bvi=bvi+dbv
#print(bv)

tmmin=bvmin #espesor mínimo de muro
tmmax=bvmax #espesor máximo de muro
dtm=0.05*m #intervalo de espesor de muro
ntm=int((tmmax-tmmin)/dtm+1) #número de intervalos para tm (espesor de muro)
tm=np.zeros(ntm) #lista de tm (ancho de muro)
tmi=tmmin
for i in range(ntm):
    tm[i]=tmi
    tmi=tmi+dtm
#print(tm)

Lmmin=4*tmmin #Longitud mínima de muro
Lmmax=dxy/2 #longitud máxima de muro
dLm=0.25*m #intervalo de longitud de muro
nLm=int((Lmmax-Lmmin)/dLm+1) #número de intervalos para Lm (longitud de muro)
LM=np.zeros(nLm) #lista de Lm (longitud de muro)
Lmi=Lmmin
for i in range(nLm):
    LM[i]=Lmi
    Lmi=Lmi+dLm
#print(Lm)

hvmin=0.25*m #peralte mínimo de viga
hvmax=2*bvmax #peralte máximo de viga
dhv=0.05*m #intervalo de lperalte de viga
nhv=int((hvmax-hvmin)/dhv+1) #número de intervalos para hv (peralte de viga)
hv=np.zeros(nhv) #lista de hv (peralte de viga)
hvi=hvmin
for i in range(nhv):
    hv[i]=hvi
    hvi=hvi+dhv
#print(hv)

lstf=np.empty((0, 6))

total=nac*nbv*ntm*nLm*nhv
# print(total)
total1=int(total/nLm)
#print(total1)
lst1=np.zeros((total1,6))

lst2=np.zeros((total,6))

cont=0
for j in range(nLm):
    for k in range(nac):
        for ii in range(nbv):
            for jj in range(ntm):
                for kk in range(nhv):
                    a = ac[k]
                    b = bv[ii]
                    t = tm[jj]
                    Lmx, Lmy = LM[j] ,LM[j]
                    h = hv[kk]

                    #Para nz = 1
                    nz = 1
                    
                    VolVx = (b*h*dx*(nx*(ny+1)-4))*nz
                    VolVy = (b*h*dy*(ny*(nx+1)-4))*nz
                    VolCol = (a*a*dz*(nx+1)*(ny+1)+8-12)*nz
                    VolMurosx = t*Lmx*dz*4*nz
                    VolMurosy = t*Lmy*dz*4*nz
                    Volumen = round(VolVx+VolVy+VolCol+VolMurosx+VolMurosy,2)

                    lst1[cont]=[a,b,t,Lmx,h,Volumen]
                    cont=cont+1
    lst1[lst1[:, 5].argsort()]
    #print(lst1)
    #print(cont)
    cont=0
    for i in range(total1):
        lst2[i+j*total1]=lst1[i]
                
for i in range(total):
    #CONDICIONALES
    #tm tiene como minimo bv
    if lst2[i,2] >= lst2[i,1]:
        #Lm tiene como minimo 4tm
        if lst2[i,3] >= 4*lst2[i,2]:
            #hv tiene como minimo 2bv
            if lst2[i,4] >= 2*lst2[i,1]:
                
                b=lst2[i,1]
                h=lst2[i,4]
                Izv = b*h**3/12
                a=lst2[i,0]
                Izc = a**4/12
                
                if Izv <= Izc:
                    
                    lstf=np.append(lstf,[lst2[i]],axis=0)

nlstf=len(lstf)

print('Número de modelos a analizar por cantidad de pisos (de 5 a 10 pisos) = %.0f'%(nlstf))
np.savetxt('lista_filtro.txt',lstf)
#==================================================================
#DEFINICIÓN DE FUNCIONES
def GeoModel(dx, dy, h, nx, ny, nz, Lmx, Lmy):
    from numpy import zeros
    Lx, Ly, Lz = dx*nx, dy*ny, h*nz
    NN = (nx+1)*(ny+1)*(nz+1)+8*(nz+1)
    Nodes = zeros((NN,5))
    #Calculando fracciones de masas en los nodos
    mnodes = zeros(5)
    mnodes[0] = (Lmx**2+Lmy**2)/(8*dx*dy) #nudo esquina
    mnodes[1] = 0.25 - (Lmx**2+dy**2-(dx+dy-Lmx)*Lmx)/(8*dx*dy) #nudo extremo de muro en x
    mnodes[2] = 0.50 - (dx+dy-Lmx)*Lmx/(8*dx*dy) #nudo lateral en x
    mnodes[3] = 0.25 - (Lmy**2+dx**2-(dx+dy-Lmy)*Lmy)/(8*dx*dy) #nudo extremo de muro en y
    mnodes[4] = 0.50 - (dx+dy-Lmy)*Lmy/(8*dx*dy) #nudo lateral en y

    # Creando los nodos y asignando coordenadas
    c = 0
    for i in range(nz+1):
        for j in range(ny+1+2):
            for k in range(nx+1+2):
                #Esquinas
                if k == 0 and j == 0:
                    Nodes[c] = [c,0,0,i*h,mnodes[0]]
                    c = c + 1
                elif k == 0 and j == ny+2:
                    Nodes[c] = [c,0,ny*dy,i*h,mnodes[0]]
                    c = c + 1
                elif k == nx+2 and j == 0:
                    Nodes[c] = [c,nx*dx,0,i*h,mnodes[0]]
                    c = c + 1
                elif k == nx+2 and j == ny+2:
                    Nodes[c] = [c,nx*dx,ny*dy,i*h,mnodes[0]]
                    c = c + 1
                #Lm en x
                elif k == 1 and j == 0:
                    Nodes[c] = [c,Lmx,0,i*h,mnodes[1]]
                    c = c + 1
                elif k == 1 and j == ny+2:
                    Nodes[c] = [c,Lmx,ny*dy,i*h,mnodes[1]]
                    c = c + 1
                elif k == nx+1 and j == 0:
                    Nodes[c] = [c,nx*dx-Lmx,0,i*h,mnodes[1]]
                    c = c + 1
                elif k == nx+1 and j == ny+2:
                    Nodes[c] = [c,nx*dx-Lmx,ny*dy,i*h,mnodes[1]]
                    c = c + 1
                #Lm en y
                elif k == 0 and j == 1:
                    Nodes[c] = [c,0,Lmy,i*h,mnodes[3]]
                    c = c + 1
                elif k == nx+2 and j == 1:
                    Nodes[c] = [c,nx*dx,Lmy,i*h,mnodes[3]]
                    c = c + 1
                elif k == 0 and j == ny+1:
                    Nodes[c] = [c,0,ny*dy-Lmy,i*h,mnodes[3]]
                    c = c + 1
                elif k == nx+2 and j == ny+1:
                    Nodes[c] = [c,nx*dx,ny*dy-Lmy,i*h,mnodes[3]]
                    c = c + 1
                #Lm en x que no son puntos
                elif k == 1 and j != 0 and j != ny+2:
                    c = c
                elif k == nx+1 and j != 0 and j != ny+2:
                    c = c
                #Lm en y que no son puntos
                elif j == 1 and k != 0 and k != nx+2:
                    c = c
                elif j == ny+1 and k != 0 and k != nx+2:
                    c = c
                #Medios en zonas de esquina
                elif k == 0 and (j == 2 or j == ny):
                    Nodes[c] = [c,0,(j-1)*dy,i*h,mnodes[4]]
                    c = c + 1
                elif k == nx+2 and (j == 2 or j == ny):
                    Nodes[c] = [c,nx*dx,(j-1)*dy,i*h,mnodes[4]]
                    c = c + 1
                elif j == 0 and (k == 2 or k == nx):
                    Nodes[c] = [c,(k-1)*dx,0,i*h,mnodes[2]]
                    c = c + 1
                elif j == ny+2 and (k == 2 or k == nx):
                    Nodes[c] = [c,(k-1)*dx,ny*dy,i*h,mnodes[2]]
                    c = c + 1
                #Medios centrales
                elif k == 0 and j != 0 and j != 1 and j != 2 and j != ny and j != ny+1 and j != ny+2:
                    Nodes[c] = [c,0,(j-1)*dy,i*h,0.50]
                    c = c + 1
                elif k == nx+2 and j != 0 and j != 1 and j != 2 and j != ny and j != ny+1 and j != ny+2:
                    Nodes[c] = [c,nx*dx,(j-1)*dy,i*h,0.50]
                    c = c + 1
                elif j == 0 and k != 0 and k != 1 and k != 2 and k != nx and k != nx+1 and k != nx+2:
                    Nodes[c] = [c,(k-1)*dx,0,i*h,0.50]
                    c = c + 1
                elif j == ny+2 and k != 0 and k != 1 and k != 2 and k != nx and k != nx+1 and k != nx+2:
                    Nodes[c] = [c,(k-1)*dx,ny*dy,i*h,0.50]
                    c = c + 1
                #Centros
                else:
                    Nodes[c] = [c,(k-1)*dx,(j-1)*dy,i*h,1.00]
                    c = c+1
    Nodes[:((nx+1)*(ny+1)+8),4]=0


    NE = (nx*(ny+1)+ny*(nx+1)+(nx+1)*(ny+1)-4)*nz+8*nz
    Elems = zeros((NE,6))
    # Creando las conexiones de los elementos verticales
    c = 0
    nn = (nx+1)*(ny+1)+8
    #Muros
    for i in range(nz):
        Elems[c] = [c,0+nn*i,1+nn*i,1+nn*(i+1),nn*(i+1),3] #3 muro en x #Asignar los nodos de forma antihoraria
        c = c + 1
        Elems[c] = [c,nx+1+nn*i,nx+2+nn*i,nx+2+nn*(i+1),nx+1+nn*(i+1),3]
        c = c + 1
        Elems[c] = [c,0+nn*i,nx+3+nn*i,nx+3+nn*(i+1),nn*(i+1),4] #4 muro en y
        c = c + 1
        Elems[c] = [c,nx+2+nn*i,nx+4+nn*i,nx+4+nn*(i+1),nx+2+nn*(i+1),4]
        c = c + 1

        Elems[c] = [c,nx+5+(nx+1)*(ny-1)+nn*i,nx+5+(nx+1)*(ny-1)+2+nn*i,nx+5+(nx+1)*(ny-1)+2+nn*(i+1),nx+5+(nx+1)*(ny-1)+nn*(i+1),4]
        c = c + 1
        Elems[c] = [c,nx+5+(nx+1)*(ny-1)+1+nn*i,nx+5+(nx+1)*(ny-1)+1+nx+3+nn*i,nx+5+(nx+1)*(ny-1)+1+nx+3+nn*(i+1),nx+5+(nx+1)*(ny-1)+1+nn*(i+1),4]
        c = c + 1
        Elems[c] = [c,nn-1-(nx+2)+nn*i,nn-(nx+2)+nn*i,nn-(nx+2)+nn*(i+1),nn-1-(nx+2)+nn*(i+1),3]
        c = c + 1
        Elems[c] = [c,nn-2+nn*i,nn-1+nn*i,nn-1+nn*(i+1),nn-2+nn*(i+1),3]
        c = c + 1


    #Columnas
    for i in range(nz):
        for j in range(ny+1):
            for k in range(nx+1):
                if j == 0 and k != 0 and k != nx:
                    Elems[c] = [c,k+1+nn*i,k+1+nn*(i+1),0,0,1]
                    c = c + 1
                elif j != 0 and j != ny:
                    Elems[c] = [c,k+nx+5+(nx+1)*(j-1)+nn*i,k+nx+5+(nx+1)*(j-1)+nn*(i+1),0,0,1]
                    c = c + 1
                elif j == ny and k != 0 and k != nx:
                    Elems[c] = [c,k+nx+5+(nx+1)*(j-1)+3+nn*i,(k+nx+5+(nx+1)*(j-1)+3)+nn*(i+1),0,0,1]
                    c = c + 1
    # Creando las conexiones de los elementos horizontales en x
    for i in range(nz):
        for j in range(ny+1):
            for k in range(nx):
                if j == 0:
                    Elems[c] = [c,k+1+nn*(i+1),k+2+nn*(i+1),0,0,2]
                    c = c + 1
                elif j != 0 and j != ny:
                    Elems[c] = [c,k+nx+5+(nx+1)*(j-1)+nn*(i+1),k+nx+5+1+(nx+1)*(j-1)+nn*(i+1),0,0,2]
                    c = c + 1
                elif j == ny:
                    Elems[c] = [c,k+nx+5+(nx+1)*(j-1)+3+nn*(i+1),k+nx+5+(nx+1)*(j-1)+4+nn*(i+1),0,0,2]
                    c = c + 1

    # Creando las conexiones de los elementos horizontales en y
    for i in range(nz):
        for j in range(nx+1):
            for k in range(ny):
                if j == 0 and k==0:
                    Elems[c] = [c,nx+3+nn*(i+1),nx+5+nn*(i+1),0,0,2]
                    c = c + 1
                elif j == 0 and k!=0:
                    Elems[c] = [c,nx+5+(nx+1)*(k-1)+nn*(i+1),nx+5+(nx+1)*k+nn*(i+1),0,0,2]
                    c = c + 1
                elif j != 0 and j != nx and k==0:
                    Elems[c] = [c,j+1+nn*(i+1),j+1+nx+4+nn*(i+1),0,0,2]
                    c = c + 1
                elif j != 0 and j != nx and k !=0 and k != ny-1:
                    Elems[c] = [c,j+nx+5+(nx+1)*(k-1)+nn*(i+1),j+nx+5+(nx+1)*k+nn*(i+1),0,0,2]
                    c = c + 1
                elif j != 0 and j != nx and k == ny-1:
                    Elems[c] = [c,j+nx+5+(nx+1)*(k-1)+nn*(i+1),j+nx+5+(nx+1)*k+3+nn*(i+1),0,0,2]
                    c = c + 1
                elif j == nx and k==0:
                    Elems[c] = [c,nx+4+nn*(i+1),nx+4+(nx+1)+nn*(i+1),0,0,2]
                    c = c + 1
                elif j == nx and k != 0 and k != ny-1:
                    Elems[c] = [c,nx+4+(nx+1)+(nx+1)*(k-1)+nn*(i+1),nx+4+(nx+1)+(nx+1)*k+nn*(i+1),0,0,2]
                    c = c + 1
                elif j == nx and k == ny-1:
                    Elems[c] = [c,nx+4+(nx+1)+(nx+1)*(k-1)+nn*(i+1),nx+4+(nx+1)+(nx+1)*(k-1)+2+nn*(i+1),0,0,2]
                    c = c + 1

    # Creando centro de diafragmas
    Diap = zeros((nz,4))
    for i in range(nz):
        Diap[i] = [i+1000,Lx/2.0,Ly/2.0,h*(i+1)]
    #
    return Nodes, Elems, Diap

def espectro_E030(T,Z=0.45,U=1.5,S=1.0,Tp=0.4,Tl=2.5,R=1):
    from numpy import zeros
    n = len(T)
    E030 = zeros(n)
    for i in range(n):
        if T[i]>=0 and T[i]<0.2*Tp:
            E030[i]=2.5#1+7.5*T[i]/Tp
        elif T[i]>=0.2*Tp and T[i]<Tp:
            E030[i]=2.5
        elif T[i]>=Tp and T[i]<Tl:
            E030[i] = 2.5*(Tp/T[i])
        elif T[i]>=Tl:
            E030[i] = 2.5*(Tp*Tl/T[i]**2)
        else:
            print("El periodo no puede ser negativo!")
    return E030*Z*U*S/R

def get_static_loads(coef,p,h,T):
    from numpy import zeros
    n = len(h)
    V = coef*sum(p)
    F = zeros(n)
    #
    if T > 0.0 and T <= 0.5:
        k=1.0
    elif T>0.5:
        k = 0.75+0.5*T
    else:
        print('El periodo es negativo!')
    #
    div = 0.
    for i in range(n):
        div = div + p[i]*h[i]**k
    #
    for i in range(n):
        F[i] = p[i]*h[i]**k/div*V
    return F,k

def getCombo(E030,MF,modo,Tmodes,NT,ni):

    # Definimos valores iniciales
    D_ABSx,D_RCSCx = np.zeros(NT),np.zeros(NT)
    Δ_ABSx,Δ_RCSCx = np.zeros(NT),np.zeros(NT)
    V_ABSx,V_RCSCx = np.zeros(NT),np.zeros(NT)
    D_ABSy,D_RCSCy = np.zeros(NT),np.zeros(NT)
    Δ_ABSy,Δ_RCSCy = np.zeros(NT),np.zeros(NT)
    V_ABSy,V_RCSCy = np.zeros(NT),np.zeros(NT)

    # Se realiza la Superpocisión Modal Espectral
    for j in range(1,ni+1):#ni+1
        FPx=modo[j-1].T@MF@Ux
        FPy=modo[j-1].T@MF@Uy
        FPr=modo[j-1].T@MF@Rz
        #
        Sa = E030[j-1]*9.80665
        Sd = Sa/(2*np.pi/Tmodes[j-1])**2
        #
        respDX = Sd*FPx*modo[j-1]
        respAX = Sa*FPx*MF@modo[j-1]
        D_ABSx = D_ABSx + abs(respDX)
        D_RCSCx = D_RCSCx + (respDX)**2
        respDX[3:] = respDX[3:] - respDX[:-3]
        Δ_ABSx = Δ_ABSx + abs(respDX)
        Δ_RCSCx = Δ_RCSCx + (respDX)**2
        V_ABSx = V_ABSx + abs(np.cumsum(respAX[::-1])[::-1])
        V_RCSCx = V_RCSCx + (np.cumsum(respAX[::-1])[::-1])**2

        #
        respDY = Sd*FPy*modo[j-1]
        respAY = Sa*FPy*MF@modo[j-1]
        D_ABSy = D_ABSy + abs(respDY)
        D_RCSCy = D_RCSCy + (respDY)**2
        respDY[3:] = respDY[3:] - respDY[:-3]
        Δ_ABSy = Δ_ABSy + abs(respDY)
        Δ_RCSCy = Δ_RCSCy + (respDY)**2
        V_ABSy = V_ABSy + abs(np.cumsum(respAY[::-1])[::-1])
        V_RCSCy = V_RCSCy + (np.cumsum(respAY[::-1])[::-1])**2
    
    # Se realiza la combinación 25%ABS + 75%RCSC
    D_RCSCx = D_RCSCx**0.5
    Δ_RCSCx = Δ_RCSCx**0.5
    V_RCSCx = V_RCSCx**0.5
    DDx = 0.25*D_ABSx + 0.75*D_RCSCx
    ΔDx = 0.25*Δ_ABSx + 0.75*Δ_RCSCx
    VDx = 0.25*V_ABSx + 0.75*V_RCSCx
    #
    D_RCSCy = D_RCSCy**0.5
    Δ_RCSCy = Δ_RCSCy**0.5
    V_RCSCy = V_RCSCy**0.5
    DDy = 0.25*D_ABSy + 0.75*D_RCSCy
    ΔDy = 0.25*Δ_ABSy + 0.75*Δ_RCSCy
    VDy = 0.25*V_ABSy + 0.75*V_RCSCy
    

    df = pd.DataFrame(columns=['Nivel','VDx(kN)','VDy(kN)','UDx(cm)','UDy(cm)'])
    for i in range(int(NT/3)):
        VDy[1::3][i]=VDx[0::3][i]

    for i in range(int(NT/3)):
        df = df.append({'Nivel':i+1, 'VDx(kN)':VDx[0::3][i]/1000,
        'VDy(kN)':VDy[1::3][i]/1000,'UDx(cm)':DDx[0::3][i]*100,
        'UDy(cm)':DDy[1::3][i]*100}, ignore_index=True)

    return DDx, ΔDx, VDx, DDy, ΔDy, VDy, df

#====================================================================
#PROPIEDADES DE MATERIALES Y SECCIONES
fc = 210*kgf/cm**2
E = 15100*(fc/(kgf/cm**2))**0.5*kgf/cm**2
G = 0.5*E/(1+0.2)
# Densidad del concreto
ρ = 2400*kg/m**3
#===================================================================
#ITERACIONES
contador=0
df6 = pd.DataFrame(columns=['#Niveles','Nmodel','a(m)','b(m)','h(m)','t(m)','Lm(m)','V(kN)','Δmax(‰)','VolConc(m3)']) #Para guardar los datos
for zz in range(nnz):
    nz = int(Nz[zz])
    Lmopt = 0
    for nmodel in range(nlstf) :
        a = lstf[nmodel,0]
        b = lstf[nmodel,1]
        t = lstf[nmodel,2]
        Lm = lstf[nmodel,3]
        h = lstf[nmodel,4]
        Volumen = lstf[nmodel,5]*nz

        print('Modelo N° %.0f para edificio de %.0f niveles'%(nmodel+1,nz))

        #CONDICIONAL PARA PASAR AL SIGUIENTE LM DEBES DE ENCONTRAR EL MODELO ÓPTIMO PARA UN LM
        if Lm == Lmopt :
            continue

        # Viga
        Av = b*h
        Izv = b*h**3/12
        Iyv = b**3*h/12
        aa, bb = max(b,h),min(b,h)
        β= 1/3-0.21*bb/aa*(1-(bb/aa)**4/12)
        Jxxv = β*bb**3*aa
        # Columna
        Ac = a**2
        Izc = a**4/12
        Iyc = a**4/12
        β= 1/3-0.21*1.*(1-(1.)**4/12)
        Jxxc = β*a**4

        #CONDICIONALES
        #if t > b :
        #    continue
        #elif Lm < 4*t :
        #    continue
        #elif Izv > Izc :
        #    continue
        
        #CREACIÓN DEL MODELO
        ops.wipe()
        ops.model('basic', '-ndm', 3, '-ndf', 6)
        RigidDiaphragm = 'ON'

        #Muro
        Lmx, Lmy = Lm , Lm
        ops.uniaxialMaterial('Elastic', 1, E) #Concreto Axial
        ops.uniaxialMaterial('Elastic', 2, 2*10**6*kgf/cm**2) #Acero
        ops.uniaxialMaterial('Elastic', 3, G) #Concreto Cortante
        #Muros en x
        ancho = 20*cm
        mufx = int(round(Lmx/ancho))
        ttx = np.zeros(mufx)
        ttx[:] = t
        wwx = np.zeros(mufx)
        wwx[:] = Lmx/(mufx)
        ρρx = np.zeros(mufx)
        ρρx[:] = 0 #Cuantía vertical en muros
        concx = np.zeros(mufx)
        concx[:] = int(1)
        acerox = np.zeros(mufx)
        acerox[:] = int(2)
        #Muros en y
        ancho = 20*cm
        mufy = int(round(Lmy/ancho))
        tty = np.zeros(mufy)
        tty[:] = t
        wwy = np.zeros(mufy)
        wwy[:] = Lmy/(mufy)
        ρρy = np.zeros(mufy)
        ρρy[:] = 0 #Cuantía vertical en muros
        concy = np.zeros(mufy)
        concy[:] = int(1)
        aceroy = np.zeros(mufy)
        aceroy[:] = int(2)

        # Nodos del Modelo
        Nodes, Elems, Diap = GeoModel(dx,dy,dz,nx,ny,nz,Lmx,Lmy)

        #CREAMOS NODOS DEL MODELO
        # Creamos los nodos
        for Ni in Nodes:
            ops.node(int(Ni[0]), *Ni[1:4])

        # Definimos diafragmas rígidos
        if RigidDiaphragm == 'ON':
            dirDia = 3 # perpendicular al plano del diafragma
            for Nd in Diap:
                ops.node(int(Nd[0]), *Nd[1:4])
                ops.fix(int(Nd[0]),*[0,0,1,1,1,0])
                NodesDi = []
                for Ni in Nodes:
                    if Ni[3]==Nd[3]:
                        NodesDi.append(int(Ni[0]))
                ops.rigidDiaphragm(dirDia,int(Nd[0]),*NodesDi)
        
        #ASIGNAMOS RESTRICCIONES EN LA BASE
        # Restricciones
        ops.fixZ(0.0, *[1,1,1,1,1,1], '-tol', 1e-6)

        #EJES LOCALES
        #Establecemos transformación geométrica
        ops.geomTransf('PDelta', int(1), *[1, 0, 0])
        ops.geomTransf('Linear', int(2), *[1,-1, 0])

        #DEFINIMOS ELEMENTOS CON SUS PROPIEDADES
        # Creamos los elementos
        for Ele in Elems:
            if int(Ele[5]) == 1: # 1 Columna
                ops.element('elasticBeamColumn', int(Ele[0]), int(Ele[1]), int(Ele[2]), Ac, E, G, Jxxc, Iyc, Izc, int(Ele[5]),'-mass', ρ*Ac)
            elif int(Ele[5]) == 2: # 2 Viga
                ops.element('elasticBeamColumn', int(Ele[0]), int(Ele[1]), int(Ele[2]), Av, E, G, Jxxv, Iyv, Izv, int(Ele[5]),'-mass', ρ*Av)#*(dx-a)/dx)
            elif int(Ele[5]) == 3: # 3 Muro en x
                ops.element('MVLEM_3D', int(Ele[0]), int(Ele[1]), int(Ele[2]), int(Ele[3]), int(Ele[4]), mufx, '-thick', *ttx[:], '-width', *wwx[:], '-rho', *ρρx[:], '-matConcrete', *concx[:], '-matSteel', *acerox[:], '-matShear', int(3), '-Poisson', 0.2, '-Density', ρ)
            elif int(Ele[5]) == 4: # 4 Muro en y
                ops.element('MVLEM_3D', int(Ele[0]), int(Ele[1]), int(Ele[2]), int(Ele[3]), int(Ele[4]), mufy, '-thick', *tty[:], '-width', *wwy[:], '-rho', *ρρy[:], '-matConcrete', *concy[:], '-matSteel', *aceroy[:], '-matShear', int(3), '-Poisson', 0.2, '-Density', ρ)


        #ASIGNACIÓN DE MASAS Y MODOS DE VIBRACIÓN
        # Aplicando Cargas vivas y muertas
        wLive = 250*kg/m**2
        wLosa = 300*kg/m**2
        wAcab = 100*kg/m**2
        wTabi = 150*kg/m**2
        wTotal = 1.0*(wLosa+wAcab+wTabi)+0.25*wLive

        Carga = wTotal*dx*dy*m**2
        for Ni in Nodes:
            ops.mass(int(Ni[0]),Ni[4]*Carga,Ni[4]*Carga,0.0)

        # Obtenemos los modos
        Nmodes = 3*nz

        vals = ops.eigen(Nmodes)
        Tmodes = np.zeros(len(vals))
        for i in range(Nmodes):
            Tmodes[i] = 2*np.pi/vals[i]**0.5
            #print("T[%i]: %.5f"%(i+1,Tmodes[i]))

        #ANÁLISIS PARA OBTENER LA MATRIZ DE MASAS
        # Realizamos un análisis para obtener la matriz de Masas
        ops.wipeAnalysis()
        ops.system('FullGeneral')
        ops.numberer("Plain")
        ops.constraints('Transformation') 
        ops.algorithm('Linear')
        ops.analysis('Transient')
        ops.integrator('GimmeMCK',1.0,0.0,0.0)
        ops.analyze(1,0.0) 

        # Obtenemos la matriz de Masas
        N = ops.systemSize()         # Número de Grados de Libertad
        Mmatrix = ops.printA('-ret')
        Mmatrix = np.array(Mmatrix).reshape((N,N))
        MF = Mmatrix[-3*nz:,-3*nz:]
        #print(MF)

        #ANÁLISIS ESTÁTICO
        #Análisis Estático en X
        np.set_printoptions(precision=3,linewidth=300,suppress=True)
        H = np.arange(1,nz+1)*dz
        P = sum(MF[0::3,0::3])*9.80665 # Peso por nivel
        #print(H,P)
        Ro = 7.0
        E030 = espectro_E030(Tmodes,Z=0.45,U=1.0,S=1.0,Tp=0.4,Tl=2.5,R=Ro)
        F, k = get_static_loads(E030[0],P,H,Tmodes[0])
        CR = E030[0]/(0.45*1.*1.)
        #print('C/R=',CR)
        #print(E030[0],k)

        #CONDICIONALES
        if CR < 0.11 :
            #print('C/R=',CR)
            continue

        ##Aplicamos fuerzas nodales
        ops.timeSeries('Linear',1)
        ops.pattern('Plain',1,1)
        Le = ny*dy*0.05
        for i in range(nz):
            #print(int(Diap[i][0]))
            ops.load(int(Diap[i][0]),F[i],0.,0.,0.,0.,F[i]*Le)
        
        #Realizamos el análisis
        ops.wipeAnalysis()
        ops.constraints('Transformation')
        ops.numberer('Plain')
        ops.system('FullGeneral')
        ops.algorithm('Linear')
        ops.integrator('LoadControl',1)
        ops.analysis('Static')
        ops.analyze(1)

        #Calculando cortantes
        VSx = np.cumsum(F[::-1])[::-1]

        #Resultados del análisis estático en X
        # Desplazamientos
        df1 = pd.DataFrame(columns=['Nivel','Vx(kN)','UxMax(cm)','UyMax(cm)','DriftX(‰)','DriftY(‰)'])
        tempX, tempY = 0., 0.
        for i in range(nz):
            desX = ops.nodeDisp(int(Diap[i][0]),1)
            desY = ops.nodeDisp(int(Diap[i][0]),2)
            rotZ = ops.nodeDisp(int(Diap[i][0]),6)
            desX = desX + abs(rotZ*ny*dy/2)
            desY = desY + abs(rotZ*nx*dx/2)
            desX, desY = desX*0.75*Ro, desY*0.75*Ro
            driftX = 1000.*(desX-tempX)/dz
            driftY = 1000.*(desY-tempY)/dz 
            tempX, tempY = desX, desY
            df1 = df1.append({'Nivel':i+1,'Vx(kN)':VSx[i]/1000,'UxMax(cm)':desX*100,'UyMax(cm)':desY*100,
                            'DriftX(‰)':driftX,'DriftY(‰)':driftY}, ignore_index=True)
        #print('\nANÁLISIS ESTÁTICO EN X')
        #print(df1.round(4))

        ops.reactions()
        Vmurox=0
        Vcolumx=0
        for i in range((nx+1)*(ny+1)+8):
            if i == 0 or i == 1 or i == nx+1 or i == nx+2 or i == nx+3 or i == nx+4 or i == ((nx+1)*(ny+1)+8)-1 or i == ((nx+1)*(ny+1)+8)-2 or i == ((nx+1)*(ny+1)+8)-2-nx or i == ((nx+1)*(ny+1)+8)-2-nx-1 or i == ((nx+1)*(ny+1)+8)-2-nx-2 or i == ((nx+1)*(ny+1)+8)-2-nx-3:
                Vmurox = Vmurox + ops.nodeReaction(i,1)
            else: Vcolumx = Vcolumx + ops.nodeReaction(i,1)
        PVmurox = abs(Vmurox)/VSx[0]*100
        PVcolumx = abs(Vcolumx)/VSx[0]*100

        #print('Vmurox = %.4f kN = %.2f%% de VSx'%(abs(Vmurox/kN),PVmurox))
        #print('Vcolumx = %.4f kN = %.2f%% de VSx'%(abs(Vcolumx/kN),PVcolumx))

        #CONDICIONALES
        if PVmurox <= 20 or PVmurox >= 80 :
            #print('PVmurox=',PVmurox)
            continue

        #Análisis Estático en Y
        np.set_printoptions(precision=3,linewidth=300,suppress=True)
        H = np.arange(1,nz+1)*dz
        P = sum(MF[0::3,0::3])*9.80665 # Peso por nivel
        #print(H,P)
        Ro = 7.0
        E030 = espectro_E030(Tmodes,Z=0.45,U=1.0,S=1.0,Tp=0.4,Tl=2.5,R=Ro)
        F, k = get_static_loads(E030[0],P,H,Tmodes[0])
        CR = E030[0]/(0.45*1.*1.)
        #print('C/R=',CR)
        #print(E030[0],k)

        #CONDICIONALES
        if CR < 0.11 :
            #print('C/R=',CR)
            continue

        ##Aplicamos fuerzas nodales
        ops.loadConst('-time', 0.0)
        ops.remove('timeSeries',1)
        ops.remove('loadPattern',1)
        ops.timeSeries('Linear',1)
        ops.pattern('Plain',1,1)
        Le = nx*dx*0.05
        for i in range(nz):
            #print(int(Diap[i][0]))
            ops.load(int(Diap[i][0]),0.,F[i],0.,0.,0.,F[i]*Le)


        #Realizamos el análisis
        ops.wipeAnalysis()
        ops.constraints('Transformation')
        ops.numberer('Plain')
        ops.system('FullGeneral')
        ops.algorithm('Linear')
        ops.integrator('LoadControl',1)
        ops.analysis('Static')
        ops.analyze(1)

        #Calculando cortantes
        VSy = np.cumsum(F[::-1])[::-1]

        #Resultados del análisis estático en Y
        # Desplazamientos
        df2 = pd.DataFrame(columns=['Nivel','Vy(kN)','UxMax(cm)','UyMax(cm)','DriftX(‰)','DriftY(‰)'])
        tempX, tempY = 0., 0.
        for i in range(nz):
            desX = ops.nodeDisp(int(Diap[i][0]),1)
            desY = ops.nodeDisp(int(Diap[i][0]),2)
            rotZ = ops.nodeDisp(int(Diap[i][0]),6)
            desX = desX + abs(rotZ*ny*dy/2)
            desY = desY + abs(rotZ*nx*dx/2)
            desX, desY = desX*0.75*Ro, desY*0.75*Ro
            driftX = 1000.*(desX-tempX)/dz
            driftY = 1000.*(desY-tempY)/dz 
            tempX, tempY = desX, desY
            df2 = df2.append({'Nivel':i+1,'Vy(kN)':VSy[i]/1000,'UxMax(cm)':desX*100,'UyMax(cm)':desY*100,
                            'DriftX(‰)':driftX,'DriftY(‰)':driftY}, ignore_index=True)
        #print('\nANÁLISIS ESTÁTICO EN Y')
        #print(df2.round(4))

        ops.reactions()
        Vmuroy=0
        Vcolumy=0
        for i in range((nx+1)*(ny+1)+8):
            if i == 0 or i == 1 or i == nx+1 or i == nx+2 or i == nx+3 or i == nx+4 or i == ((nx+1)*(ny+1)+8)-1 or i == ((nx+1)*(ny+1)+8)-2 or i == ((nx+1)*(ny+1)+8)-2-nx or i == ((nx+1)*(ny+1)+8)-2-nx-1 or i == ((nx+1)*(ny+1)+8)-2-nx-2 or i == ((nx+1)*(ny+1)+8)-2-nx-3:
                Vmuroy = Vmuroy + ops.nodeReaction(i,2)
            else: Vcolumy = Vcolumy + ops.nodeReaction(i,2)
        PVmuroy = abs(Vmuroy)/VSy[0]*100
        PVcolumy = abs(Vcolumy)/VSy[0]*100

        #print('Vmuroy = %.4f kN = %.2f%% de VSy'%(abs(Vmuroy/kN),PVmuroy))
        #print('Vcolumy = %.4f kN = %.2f%% de VSy'%(abs(Vcolumy/kN),PVcolumy))

        #CONDICIONALES
        if PVmuroy <= 20 or PVmuroy >= 80 :
            #print('PVmuroy=',PVmuroy)
            continue

        #MASAS EFECTIVAS
        Tags = ops.getNodeTags()
        # print(Tags)
        modo = np.zeros((Nmodes,3*nz))
        for j in range(1,Nmodes+1):
            ind = 0
            for i in Tags[-nz:]:
                temp = ops.nodeEigenvector(i,j)
                modo[j-1,[ind,ind+1,ind+2]] = temp[0],temp[1],temp[-1]
                ind = ind + 3

        # Definimos valores iniciales
        Ux,Uy,Rz = np.zeros(3*nz),np.zeros(3*nz),np.zeros(3*nz)
        Ux[0::3]=1
        Uy[1::3]=1
        Rz[2::3]=1
        SUMx, SUMy, SUMr = 0., 0., 0.
        ni = 0

        Mx = sum(sum(MF[0::3,0::3]))
        My = sum(sum(MF[1::3,1::3]))
        Mr = sum(sum(MF[2::3,2::3]))

        df3 = pd.DataFrame(columns=['Modo','T(s)','FPRx','FPRy','FPRr','SumUx','SumUy','SumRz'])
        for j in range(1,Nmodes+1):
            FPx = modo[j-1].T@MF@Ux
            FPy = modo[j-1].T@MF@Uy
            FPr = modo[j-1].T@MF@Rz
            FPRx = FPx**2/Mx
            FPRy = FPy**2/My
            FPRr = FPr**2/Mr
            SUMx = SUMx + FPRx
            SUMy = SUMy + FPRy
            SUMr = SUMr + FPRr
            #
            if min(SUMx,SUMy,SUMr)>=0.90 and ni==0:
                ni = j
            df3 = df3.append({'Modo':j, 'T(s)':Tmodes[j-1],'FPRx':FPRx,'FPRy':FPRy,'FPRr':FPRr,'SumUx':SUMx,
                            'SumUy':SUMy,'SumRz':SUMr}, ignore_index=True)
        #print(df3.round(5))
        #print('N° mínimo de Modos a considerar:',ni)

        #CONDICIONALES
        if ni == 0 :
            continue

        DDx, ΔDx, VDx, DDy, ΔDy, VDy, df4 = getCombo(E030,MF,modo,Tmodes,3*nz,ni)
        #print('\nANÁLISIS DINÁMICO SIN ESCALAR')
        df4 = df4.astype({'Nivel':int})
        #print(df4.round(4))

        # Escalamiento de los resultados del análisis dinámico
        if VDx[0::3][0]<0.80*VSx[0]:
            FSx  = 0.80*VSx[0]/VDx[0::3][0]
            msjx = 'SI es necesario aplicar un Factor de Escala en X: %.4f'%FSx
        else:
            FSx = 1.
            msjx = 'NO es necesario escalar en X'

        if VDy[1::3][0]<0.80*VSy[0]:
            FSy  = 0.80*VSy[0]/VDy[1::3][0]
            msjy = 'SI es necesario aplicar un Factor de Escala en Y: %.4f'%FSy
        else:
            FSy = 1.
            msjy = 'NO es necesario escalar en Y'

        texto1 = '\nAl comparar la cortante basal obtenida en el análisis dinámico en X \n\
        (%.2f kN) y el 80%% de la cortante basal del análisis estático en X (%.2f kN), \n\
        se obtiene que %s. '%(VDx[0::3][0]/1000,0.80*VSx[0]/1000,msjx)
        texto1 = texto1 + '\nEn la dirección Y, la cortante basal obtenida en el análisis \n\
        dinámico es %.2f kN y el 80%% de la cortante basal del análisis estático es %.2f kN. \n\
        Por lo que %s.'%(VDy[1::3][0]/1000,0.80*VSy[0]/1000,msjy)
        #print(texto1)

        # Se aplican los Factores de Escala
        #print('\nANÁLISIS DINÁMICO FINAL')
        df5 = pd.DataFrame(columns=['Nivel','Vx(kN)','Vy(kN)','Ux(cm)','Uy(cm)','Δx(‰)','Δy(‰)'])
        for i in range(nz):
            Δx = 0.75*Ro*ΔDx[0::3][i]/dz
            Δy = 0.75*Ro*ΔDy[1::3][i]/dz
            #
            df5 = df5.append({'Nivel':i+1, 'Vx(kN)':FSx*VDx[0::3][i]/1000,
                'Vy(kN)':FSy*VDy[1::3][i]/1000,'Ux(cm)':0.75*Ro*DDx[0::3][i]*100,
                'Uy(cm)':0.75*Ro*DDy[1::3][i]*100,'Δx(‰)':Δx*1000,'Δy(‰)':Δy*1000}, ignore_index=True)
        df5 = df5.astype({'Nivel':int})
        #print(df5.round(4))

        # Distorsiones máximas
        vecX = np.array(df5.loc[:,'Δx(‰)'])
        vecY = np.array(df5.loc[:,'Δy(‰)'])

        maxdriftx = vecX.max()
        maxdrifty = vecY.max()

        #CONDICIONALES
        if maxdriftx >= 7 or maxdrifty >= 7 :
            continue

        #Se guarda el Lm del modelo óptimo
        Lmopt = Lm

        #GUARDAMOS LOS DATOS
        print('OK')
        contador = contador + 1
        
        df6 = df6.append({'#Niveles':'%.0f'%nz,'Nmodel':'%.0f'%(nmodel+1),'a(m)':a,'b(m)':b,'h(m)':h,'t(m)':t,'Lm(m)':Lmx,'V(kN)':df5['Vx(kN)'][0],'Δmax(‰)': maxdriftx,'VolConc(m3)':Volumen}, ignore_index=True)
    
print('Número de Modelos correctos %.0f:'%(contador))

print(df6.round(4))
df6.to_csv('Modelos_analizados_v2.csv')