SEA model examples

Two rooms separated by a wall

# -*- coding: utf-8 -*-
"""
Example script for lecture example
for Ulf Orenius
"""
import numpy as np
import matplotlib.pyplot as plt

# pyva packages
import pyva.models as mds
import pyva.coupling.junctions as con

import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC


import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC

import pyva.useful as uf

mycolors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf', '#17becf', '#17becf']
 
# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels(range='SEA')
# Print option
printsw = True
#printsw = False
figpath = '../source/images/'
plt.close('all')

# wall dimensions
Ly = 4.
Lz = 2.5
S = Lz*Ly

# Additional room dimensions
Lx1 = 3.
Lx2 = 5.
# Absorption area
As1  = 8.
As2  = 10. # As2 = S

# wall thickness
h = 0.05

# junction properties
angle_S = 0
angle_R = 90*np.pi/180

# Create materials
concrete = matC.IsoMat(E=3.8e9,nu=0.33,rho0=1250.)
air = matC.Fluid()

# Create properties
concrete_5cm = stPC.PlateProp(h,concrete)

# Create plate subsystems
wall  = st2Dsys.RectangularPlate(2, Ly,Lz,prop=concrete_5cm, eta = 0.03)
room1 = ac3Dsys.RectangularRoom(1, Lx1, Ly, Lz, air, absorption_area = As1, damping_type= ['surface'] )
# use receiving room with absorption
room2 = ac3Dsys.RectangularRoom(3, Lx2, Ly, Lz, air, absorption_area = As2, damping_type= ['surface'])

J123 = con.AreaJunction((room1,wall,room2))#,area_wall)
#om = 2*np.pi*np.logspace(2,np.log10(10000),3)


# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)
om    = omega.data
freq  = om/2/np.pi 

# define load
# room power
power1mWatt = lC.Load(omega, 0.001*np.ones(omega.shape), dof.DOF(1,0,dof.DOFtype(typestr = 'power')), name = '1Watt')

# wall power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
two_rooms = mds.HybridModel((wall,room1,room2),xdata=omega)
#connect both
two_rooms.add_junction({'areaJ_12':J123})
two_rooms.add_load('1mWatt',power1mWatt) # add 1Watt per band 

#pdof = J123.get_wave_DOF()

# Calculate specific CLF factors
eta_13 = J123.CLF(om,(0,2))
eta_31 = J123.CLF(om,(2,0))
eta_12 = J123.CLF(om,(0,1),)
eta_21 = J123.CLF(om,(1,0))

eta_all = J123.junction_matrix(om)

# check tau
tau_diff = J123.non_res.transmission_diffuse(om,signal = False)



#%% plot the different modal densities
plt.figure(1)
plt.plot(freq,wall.modal_density(om,3),label = 'wall')
plt.plot(freq,room1.modal_density(om),label = 'room1')
plt.plot(freq,room2.modal_density(om),label = 'room2')
# plt.loglog(f,n_1,':',label = 'n1 VA1')
# plt.loglog(f,n_2,':',label = 'n2 VA1')
# plt.loglog(f,n_3,':',label = 'n3 VA1')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$n(\omega)/$Hz$^{-1}$')
plt.legend()
plt.show

if printsw:
     plt.savefig(figpath+'two_rooms_modal_density.png')


#%% plot modal overlap 
plt.figure(2)
plt.plot(freq,wall.modal_overlap(om,3),label = 'wall B')
plt.plot(freq,wall.modal_overlap(om,5),label = 'wall LS')
plt.plot(freq,room1.modal_overlap(om),label = 'room1')
plt.plot(freq,room2.modal_overlap(om),label = 'room2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$M$')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_modal_overlap.png')

#%% plot modes in band 
plt.figure(3)
plt.plot(freq,wall.modes_in_band(om,3),label = 'wall B')
plt.plot(freq,wall.modes_in_band(om,5),label = 'wall LS')
plt.plot(freq,room1.modes_in_band(om),label = 'room1')
plt.plot(freq,room2.modes_in_band(om),label = 'room2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('modes in band')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_modes_in_band.png')


 



#%% PLot radation efficiency
sigma = wall.radiation_efficiency(om,fluid = air)
sigma_simple = wall.radiation_efficiency_simple(om,fluid = air)

# Show coincidence frequency
f_c = concrete_5cm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

plt.figure(4)
plt.loglog(freq,sigma,label = '[Lep1982]')
plt.loglog(freq,sigma_simple,label = 'ISO EN 12354-1')
#plt.loglog(om_VA1,sigma_VA1,label = 'VAOne')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\sigma$')
#plt.ylim(0,0.005)
#plt.xticks(2*np.pi*fc,fclabels)
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_radiation_efficiency.png')


#%% Plot CLFs
plt.close(5)
plt.figure(5)
eta_13.plot(5,yscale = 'log',fulllegstr = ['$\eta_{13}$'],ls='-')
eta_31.plot(5,yscale = 'log',fulllegstr = ['$\eta_{31}$'],cs=[mycolors[1]],ls = '--')
eta_12.plot(5,yscale = 'log',fulllegstr = ['$\eta_{1,2B}$'],cs=[mycolors[2]],ls = ':',marker = [' '])
eta_21.plot(5,yscale = 'log',fulllegstr = ['$\eta_{2B,1}$'],cs=[mycolors[3]],ls = '-.',marker = [' '])
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\eta$')
plt.ylim(1E-8,1E-2)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])
plt.xscale('log')

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_etas.png')

#%% Solve model
two_rooms.create_SEA_matrix(sym = 1)
two_rooms.solve()

# Derive paths of power input
pow_in_room2 = two_rooms.power_input(3)

#%% Energry result plots
two_rooms.energy.plot(20,xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('wall B','wall LS','room1','room2'))
plt.xlabel('$f_c/$Hz')
plt.ylabel('$E/$ J')
#plt.ylim(bottom=0.1)
plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_energy.png')



#%% Engineering units velocity
plt.figure(30,figsize=(5,5/4*3))
two_rooms.result.plot(30,ID=[2],dof=[3],xscale = 'log',yscale = 'log',
                       fulllegstr = ('wall B',))
#plt.plot(om_VA1,v_wall_B,':',label = 'VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$v_{\\rm rms}/$m s$^{-1}$')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_velocity.png')

#%% engineering units pressure
plt.figure(31,figsize=(5,5/4*3))

two_rooms.result.plot(31,ID=[1,3],xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('room 1','room 2'))
#plt.plot(om_VA1,p_room1,':',label = 'room1 VA1')
#plt.plot(om_VA1,p_room2,':',label = 'room2 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_pressure.png')
    
# Calculate transmission loss from pressure difference
# Note: This is allowed because wall and absorption surface are equal S=A2
p1 = two_rooms.result[2].ydata.flatten()
p2 = two_rooms.result[3].ydata.flatten()
tau = (p2/p1)**2

#%% Plot transmission loss 

plt.figure(32)
plt.semilogx(freq,-10*np.log10(tau),label = 'wall TL')
#plt.semilogx(f,20*np.log10(tau_VA1),':',label = 'pure VA1')
#plt.semilogx(freq,-10*np.log10(mass_law(om)),':',label = 'mass law')
#plt.semilogx(freq,20*np.log10(wall.prop.mass_per_area*freq)-48,':',label = 'mass law diffuse')
plt.semilogx(freq,-10*np.log10(tau_diff),':',label = 'mass law TMM')
plt.xticks(fc[2:],fclabels[2:])
plt.xlabel('$f_c/$Hz')
plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_TL.png')

#%% Plot input power
    
pow_in_room2.plot(33,yscale = 'log',fulllegstr = ['$\Pi_{room1}$','$\Pi_{wall}$'],xscale = 'log')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_power_in.png')






    

Two rooms with floor separated by a wall

# -*- coding: utf-8 -*-
"""
Example script for lecture example
for Ulf Orenius
"""
import numpy as np
import matplotlib.pyplot as plt

# pyva packages
import pyva.models as mds
import pyva.coupling.junctions as con

import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC


import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC

import pyva.useful as uf

mycolors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf', '#17becf', '#17becf']
 
# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels(range='hybrid')
# Print option
printsw = True
printsw = False
figpath = '../source/images/'

plt.close('all')

# wall dimensions
Ly = 4.
Lz = 2.5
S = Lz*Ly

# Additional room dimensions
Lx1 = 3.
Lx2 = 5.
# Absorption area
As1  = 8.
As2  = 10. # As2 = S

# wall thickness
h   = 0.05
h_f = 0.17 

# junction properties
angle_S = 0
angle_R = 90*np.pi/180

# Create materials
concrete = matC.IsoMat(E=3.8e9,nu=0.33,rho0=1250.)
air = matC.Fluid()

# Create properties
concrete_5cm   = stPC.PlateProp(h,concrete)
concrete_17cm  = stPC.PlateProp(h_f,concrete)

# Create plate subsystems
wall   = st2Dsys.RectangularPlate(2, Ly,Lz,prop=concrete_5cm, eta = 0.03)
floor1 = st2Dsys.RectangularPlate(4, Lx1,Ly,prop=concrete_17cm, eta = 0.03)
floor2 = st2Dsys.RectangularPlate(5, Lx2,Ly,prop=concrete_17cm, eta = 0.03)
room1 = ac3Dsys.RectangularRoom(1, Lx1, Ly, Lz, air, absorption_area = As1, damping_type= ['surface'] )
# use receiving room with absorption
room2 = ac3Dsys.RectangularRoom(3, Lx2, Ly, Lz, air, absorption_area = As2, damping_type= ['surface'])

# Area Junctions
J123 = con.AreaJunction((room1,wall,room2))
J14  = con.AreaJunction((room1,floor1))
J35  = con.AreaJunction((room2,floor2))

# Line Junctions
J425 = con.LineJunction((floor1,wall,floor2),length = Ly, thetas = (0,90,180))


# Frequency range
omega = mC.DataAxis.octave_band(f_min=2*np.pi*50,f_max=2*np.pi*10000)
om    = omega.data
freq  = om/2/np.pi 

# define load
# room power
# power1mWatt = lC.Load(omega, 0.001*np.ones(omega.shape), dof.DOF(1,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
force10Nrms = lC.Load(omega, 10*np.ones(omega.shape), dof.DOF(4,3,dof.DOFtype(typestr = 'force')), name = '10N')

# wall power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
two_rooms = mds.HybridModel((wall,room1,room2,floor1,floor2),xdata=omega)
#connect all
two_rooms.add_junction({'areaJ_123':J123})
two_rooms.add_junction({'areaJ_14':J14})
two_rooms.add_junction({'areaJ_35':J35})
two_rooms.add_junction({'lineJ_425':J425})

#two_rooms.add_load('10Watt',power1mWatt) # add 1Watt per band 
two_rooms.add_load('10N',force10Nrms)# add force excitatio to wave_DOF 3 of system 4
#
# Calculate specific CLF factors
eta_13 = J123.CLF(om,(0,2))
eta_31 = J123.CLF(om,(2,0))
eta_12 = J123.CLF(om,(0,1),)
eta_21 = J123.CLF(om,(1,0))

eta_all = J123.junction_matrix(om)

# check tau
tau_diff = J123.non_res.transmission_diffuse(om,signal = False)

# read reference results (keep if ulf will check with VA1)
#f,n_1,n_2,n_3 = np.loadtxt ('2rooms_m_dens.csv',
#                    unpack = True,
#                    usecols = (0,2,1,5), skiprows = 1,
#                    delimiter = ',')

# f,sigma_VA1 = np.loadtxt ('2rooms_sigma.csv',
#                     unpack = True,
#                     usecols = (0,1), skiprows = 1,
#                     delimiter = ',')

# f,p_room1,p_room2 = np.loadtxt ('2rooms_p_rms.csv',
#                     unpack = True,
#                     usecols = (0,1,2), skiprows = 1,
#                     delimiter = ',')

# f,v_wall_B = np.loadtxt ('2rooms_v_rms.csv',
#                     unpack = True,
#                     usecols = (0,1), skiprows = 1,
#                     delimiter = ',')

# f,pow_wall,pow_room1 = np.loadtxt ('2rooms_pow_in.csv',
#                     unpack = True,
#                     usecols = (0,3,2), skiprows = 1,
#                     delimiter = ',')

# f,tau_VA1,tau_eff_VA1 = np.loadtxt ('2rooms_tau.csv',
#                     unpack = True,
#                     usecols = (0,1,2), skiprows = 1,
#                     delimiter = ',')

# f,eta_13_VA1,eta_31_VA1,eta_12_VA1,eta_21_VA1,eta_23_VA1,eta_32_VA1 = np.loadtxt ('2rooms_eta.csv',
#                     unpack = True,
#                     usecols = (0,1,2,6,5,3,4), skiprows = 1,
#                     delimiter = ',')


# om_VA1 = f*2*np.pi


#%% plot the different modal densities
plt.figure(1)
plt.plot(freq,wall.modal_density(om,3),label = 'wall')
plt.plot(freq,floor1.modal_density(om,3),label = 'floor1')
plt.plot(freq,floor2.modal_density(om,3),label = 'floor2')
# plt.loglog(f,n_1,':',label = 'n1 VA1')
# plt.loglog(f,n_2,':',label = 'n2 VA1')
# plt.loglog(f,n_3,':',label = 'n3 VA1')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$n(\omega)/$Hz$^{-1}$')
plt.legend()
plt.show

if printsw:
     plt.savefig(figpath+'two_rooms_floor_modal_density.png')


#%% plot modal overlap 
plt.figure(2)
plt.plot(freq,wall.modal_overlap(om,3),label = 'wall B')
#plt.plot(freq,wall.modal_overlap(om,5),label = 'wall LS')
plt.plot(freq,floor1.modal_overlap(om,3),label = 'floor1')
plt.plot(freq,floor2.modal_overlap(om,3),label = 'floor2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$M$')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_modal_overlap.png')

#%% plot modes in band 
plt.figure(3)
plt.plot(freq,wall.modes_in_band(om,3),label = 'wall B')
#plt.plot(freq,wall.modes_in_band(om,5),label = 'wall LS')
plt.plot(freq,floor1.modes_in_band(om,3),label = 'floor1')
plt.plot(freq,floor2.modes_in_band(om,3),label = 'floor2')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.ylabel('modes in band')
plt.ylim(bottom=0.1)
plt.xticks(fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_modes_in_band.png')


 



#%% PLot radation efficiency
sigma1 = floor1.radiation_efficiency(om,fluid = air)
sigma2 = floor2.radiation_efficiency(om,fluid = air)

# Show coincidence frequency
f_cw = concrete_5cm.coincidence_frequency(air.c0)/2/np.pi
f_cf = concrete_17cm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency wall = {0}'.format(f_cw))
print('coincidence frequency floor = {0}'.format(f_cf))

plt.figure(4)
plt.semilogx(freq,sigma1,label = 'floor1')
plt.semilogx(freq,sigma2,label = 'floor2')
#plt.loglog(om_VA1,sigma_VA1,label = 'VAOne')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\sigma$')
#plt.ylim(0,0.005)
#plt.xticks(2*np.pi*fc,fclabels)
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_radiation_efficiency.png')


#%% plot4
def mass_law(omega,theta=0):
    z0 = air.z0
    m  = wall.prop.mass_per_area
    return 1/(1+(m*omega/2/z0)**2)

def eta_f(omega,V,area):
    c0 = air.c0
    return c0*area/4/V/omega*mass_law(omega)

plt.figure(5)
eta_13.plot(5,yscale = 'log',fulllegstr = ['$\eta_{13}$'],ls='-')
eta_31.plot(5,yscale = 'log',fulllegstr = ['$\eta_{31}$'],cs=[mycolors[1]],ls = '--')
#plt.loglog(om_VA1,eta_13_VA1,':',label = '13 VA1')
#plt.loglog(om_VA1,eta_31_VA1,':',label = '31 VA1')
# #%% plot5
#eta_12.plot(4,yscale = 'log',fulllegstr = ['$\eta_{1,2B}$'],cs=[mycolors[2]],ls = ':',marker = [' '])
#eta_21.plot(4,yscale = 'log',fulllegstr = ['$\eta_{2B,1}$'],cs=[mycolors[3]],ls = '-.',marker = [' '])
#plt.loglog(om_VA1,eta_12_VA1,':',label = '13 VA1')
#plt.loglog(om_VA1,eta_21_VA1,':',label = '31 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\eta$')
#plt.ylim(0,0.005)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])
plt.xscale('log')

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show

if printsw:
    plt.savefig(figpath+'two_rooms_floor_etas.png')

#%% Solve model
two_rooms.create_SEA_matrix(sym = 1)
two_rooms.solve()

# Derive paths of power input
pow_in_room1 = two_rooms.power_input(1)
pow_in_room2 = two_rooms.power_input(3)


#%% Enrgry result plots
two_rooms.energy.plot(20,xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                        fulllegstr = ('wall B','wall LS','room1','room2','floor1 B','floor1 LS','floor2 B','floor2 LS'))
plt.xlabel('$f_c/$Hz')
plt.ylabel('$E/$ J')
#plt.ylim(bottom=0.1)
plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_energy.png')



#%% Engineering units velocity
plt.figure(30,figsize=(5,5/4*3))
two_rooms.result.plot(30,ID=[2,4,5],dof=[3,3,3],xscale = 'log',yscale = 'log',
                       fulllegstr = ('wall B','floor1 B','floor2 B'))
#plt.plot(om_VA1,v_wall_B,':',label = 'VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$v_{\\rm rms}/$m s$^{-1}$')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_velocity.png')

#%% engineering units pressure
plt.figure(31,figsize=(5,5/4*3))

two_rooms.result.plot(31,ID=[1,3],xscale = 'log',yscale = 'log',ls = ['-','--',':','-.'],
                       fulllegstr = ('room 1','room 2'))
#plt.plot(om_VA1,p_room1,':',label = 'room1 VA1')
#plt.plot(om_VA1,p_room2,':',label = 'room2 VA1')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
#plt.ylim(bottom=0.1)
ax = plt.gca()
ax.axes.xaxis.set_ticks([])

plt.xticks(2*np.pi*fc[2:],fclabels[2:])
plt.legend()
plt.tight_layout()  	
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_pressure.png')
    
# Calculate transmission loss from pressure difference
# Note: This is allowed because wall and absorption surface are equal S=A2
p1 = two_rooms.result[2].ydata.flatten()
p2 = two_rooms.result[3].ydata.flatten()
tau = (p2/p1)**2


#%% Plot input power
    
pow_in_room1.plot(32,yscale = 'log',xscale = 'log')
#plt.loglog(om_VA1,pow_wall,':')
#plt.loglog(om_VA1,pow_room1,':')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_power_in1.png')

pow_in_room2.plot(33,yscale = 'log',xscale = 'log')
#plt.loglog(om_VA1,pow_wall,':')
#plt.loglog(om_VA1,pow_room1,':')
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.show()

if printsw:
    plt.savefig(figpath+'two_rooms_floor_power_in2.png')





    

Box cover of sound source

# -*- coding: utf-8 -*-
"""
Pyva box example without isolation
"""

import numpy as np
import matplotlib.pyplot as plt

import pyva.models as mds
import pyva.coupling.junctions as jun
import pyva.properties.structuralPropertyClasses as stPC
import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC

import pyva.useful as uf

# my packages

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC
import pyva.properties.materialClasses as matC

plt.close('all')

# x-axis tics for better readability
fc,fclabels = uf.get_3rd_oct_axis_labels()
fc,fclabels = fc[1:],fclabels[1:]

# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)

# Plate dimensions
Lx = 1.2
Ly = 1
Lz = 1.1

# Box dimensions
V = Lx*Ly*Lz
A = 2*(Lx*Ly+Ly*Lz+Lx*Lz)
P = 4*(Lx+Ly+Lz)

# Plates thickness
t = 2.0E-3;

# junction angle
angle_R = (0.,90*np.pi/180)

# Create materials
steel = matC.IsoMat(E=210e9,nu=0.3,rho0=7800, eta = 0.0)
air   = matC.Fluid() 

# Create props
steel2mm = stPC.PlateProp(t,steel)


area_dof = dof.DOF(0,0,dof.DOFtype(typestr='area'))


# %% create models
f_c = steel2mm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

# Create plate subsystems
plate1 = st2Dsys.RectangularPlate(1,Lx,Lz,prop = steel2mm)
plate2 = st2Dsys.RectangularPlate(2,Ly,Lz,prop = steel2mm)
plate3 = st2Dsys.RectangularPlate(3,Lx,Lz,prop = steel2mm)
plate4 = st2Dsys.RectangularPlate(4,Ly,Lz,prop = steel2mm)
plate5 = st2Dsys.RectangularPlate(5,Lx,Ly,prop = steel2mm)

room     = ac3Dsys.Acoustic3DSystem(6, V , A, P, air)

# create semi infinite fluids
sif1 = jun.SemiInfiniteFluid((room,plate1), air)
sif2 = jun.SemiInfiniteFluid((room,plate2), air)
sif3 = jun.SemiInfiniteFluid((room,plate3), air)
sif4 = jun.SemiInfiniteFluid((room,plate4), air)
sif5 = jun.SemiInfiniteFluid((room,plate5), air)

juncs =  { 'j61' : jun.AreaJunction((room,plate1)), 
           'j62' : jun.AreaJunction((room,plate2)),
           'j63' : jun.AreaJunction((room,plate3)),
           'j64' : jun.AreaJunction((room,plate4)),
           'j65' : jun.AreaJunction((room,plate5)),
           'j12' : jun.LineJunction((plate1,plate2), Lz, angle_R),
           'j23' : jun.LineJunction((plate2,plate3), Lz, angle_R),
           'j34' : jun.LineJunction((plate3,plate4), Lz, angle_R),
           'j41' : jun.LineJunction((plate4,plate2), Lz, angle_R),
           'j15' : jun.LineJunction((plate1,plate5), Lx, angle_R),
           'j25' : jun.LineJunction((plate2,plate5), Ly, angle_R),
           'j35' : jun.LineJunction((plate3,plate5), Lx, angle_R),
           'j45' : jun.LineJunction((plate4,plate5), Ly, angle_R)
            }



#om = 2*np.pi*np.logspace(2,np.log10(10000),3)




# define load
# room power
power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(6,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
# plate power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')

#create SEA model
box = mds.HybridModel((plate1,plate2,plate3,plate4,plate5,room),xdata=omega)

#connect both
box.add_junction(juncs)
box.add_SIF({'sif1' : sif1, 
             'sif2' : sif2, 
             'sif3' : sif3, 
             'sif4' : sif4, 
             'sif5' : sif5})

box.add_load('1Watt',power1Watt) # add 1mWatt per band 


#%% solving
box.create_SEA_matrix(sym = 1)
box.solve()

#%% plotting 1

box.result.plot(1,ID = [1,5],xscale = 'log',yscale = 'log',fulllegstr=('a)','a)'))
plt.figure(1)
plt.yscale('log')
plt.ylim(1e-5,0.01)
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
#plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

#%% plotting 2
 
box.result.plot(2,ID = [6],xscale = 'log',yscale = 'linear',fulllegstr=('a)',))   
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
plt.tight_layout()
plt.legend()
plt.show()

# plt.savefig('../source/images/box_pure_pressure.png')

#%% more info
sif1_in = box.power_input('sif1')
sif2_in = box.power_input('sif2')
sif3_in = box.power_input('sif3')
sif4_in = box.power_input('sif4')
sif5_in = box.power_input('sif5')

sif_all = sif1_in.sum()+sif2_in.sum()+sif3_in.sum()+sif4_in.sum()+sif5_in.sum()


#%% plot 3

sif1_in.plot(3,xscale='log',yscale='linear')
plt.figure(3)
#plt.plot(om_VA1,power_in_1_res,':',label = 'VA1 res')
#plt.plot(om_VA1,power_in_1_nonres,':',label = 'VA1 non-res')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_pure_SIF_power_in.png')


#%% plot 4

sif_all.plot(4,xscale='log',res = 'dB',fulllegstr=('a)',))
plt.figure(4)
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
#plt.yticks([40,60,80,100,],[40,60,80,100,])
plt.ylim((60,110))
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_pure_SIF_total.png')

#%% plot 5    
NR = 1/sif_all.ydata[0,:]

plt.figure(5)
plt.semilogx(omega.data,10*np.log10(NR),label = 'a)')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$IL/$dB')
plt.legend()
plt.tight_layout()
plt.show()




    

Box cover with noise control treatment

# -*- coding: utf-8 -*-
"""
Pyva box example with isolation

"""
import numpy as np
import matplotlib.pyplot as plt

import pyva.models as mds
import pyva.coupling.junctions as con
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
import pyva.systems.structure2Dsystems as st2Dsys
import pyva.systems.acoustic3Dsystems as ac3Dsys
import pyva.loads.loadCase as lC

import pyva.useful as uf

# my packages

import pyva.data.dof as dof
import pyva.data.matrixClasses as mC
import pyva.properties.materialClasses as matC

plt.close('all')

# x-axis tics
fc,fclabels = uf.get_3rd_oct_axis_labels()
fc,fclabels = fc[1:],fclabels[1:]

# Frequency range
omega = mC.DataAxis.octave_band(f_max=2*np.pi*10000)
om    = omega.data

# Plate dimensions
Lx = 1.2
Ly = 1
Lz = 1.1

# Room dimensions
V = Lx*Ly*Lz
A = 2*(Lx*Ly+Ly*Lz+Lx*Lz)
P = 4*(Lx+Ly+Lz)


# Plate thickness
t = 2.0E-3;

# junction properties
angle_R = (0.,90*np.pi/180)

# Create materials
steel = matC.IsoMat(E=210e9,nu=0.3,rho0=7800, eta = 0.0)
air   = matC.Fluid()

# Fibre material
rho_bulk = 0.98*1.20 + 30.
fibre1 = matC.EquivalentFluid(porosity = 0.98, \
                               flow_res = 25000.,\
                               tortuosity = 1.02, \
                               length_visc = 90.e-6, \
                               length_therm = 180.e-6,\
                               rho_bulk = rho_bulk , \
                               rho0 = 1.208, \
                               dynamic_viscosity = 1.81e-5 )

# Thickness of fibre layer
h1 = 0.05

# Create props
steel2mm = stPC.PlateProp(t,steel)
# Create Layer

fibre_5cm = iL.FluidLayer(h1,fibre1)
heavy_1kg = iL.MassLayer(0.001, 1000)
# and noise control treatement
nct_mass = mds.TMmodel((fibre_5cm,heavy_1kg)) 
# change order for absorption
nct_mass4abs = mds.TMmodel((heavy_1kg,fibre_5cm))
nct          = mds.TMmodel((fibre_5cm,)) 



# absorption due to trim
alpha_nct = nct.absorption_diffuse(om,in_fluid=air)
alpha_nct_mass = nct_mass4abs.absorption_diffuse(om,in_fluid=air)

# %% plot abs fig 10 
alpha_nct.plot(10,xscale='log',yscale='linear',fulllegstr = ('5cm fibre',))
alpha_nct_mass.plot(10,xscale='log',yscale='linear',fulllegstr = ('5cm fibre + 1kg',), ls = '--')
plt.xlabel('$f_c/$Hz')
plt.ylabel('$\\alpha_s$')
plt.xticks(2*np.pi*fc,fclabels)
plt.tight_layout()
plt.legend()
plt.show()

# %% create absorption
area_dof = dof.DOF(0,0,dof.DOFtype(typestr='area'))
abs_area = alpha_nct*(Lx*Ly)+alpha_nct_mass*(A-Lx*Ly)
# This is a signal those dof must be set to the
abs_area.dof = area_dof

# %% plot 1

abs_area.plot(1,xscale = 'log')

plt.savefig('../source/images/box_isolated_abs_area.png')

# %% plates

f_c = steel2mm.coincidence_frequency(air.c0)/2/np.pi
print('coincidence frequency = {0}'.format(f_c))

# Create plate subsystems
plate1 = st2Dsys.RectangularPlate(1,Lx,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate2 = st2Dsys.RectangularPlate(2,Ly,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate3 = st2Dsys.RectangularPlate(3,Lx,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate4 = st2Dsys.RectangularPlate(4,Ly,Lz,prop = steel2mm,trim=(nct_mass,'none'))
plate5 = st2Dsys.RectangularPlate(5,Lx,Ly,prop = steel2mm,trim=(nct_mass,'none'))

room   = ac3Dsys.Acoustic3DSystem(6, V , A, P, air,\
                        absorption_area = abs_area ,\
                        damping_type= ['eta','surface']    )



# use receiving room with abso
sif1 = con.SemiInfiniteFluid((room,plate1), air)
sif2 = con.SemiInfiniteFluid((room,plate2), air)
sif3 = con.SemiInfiniteFluid((room,plate3), air)
sif4 = con.SemiInfiniteFluid((room,plate4), air)
sif5 = con.SemiInfiniteFluid((room,plate5), air)


juncs =  { 'j61' : con.AreaJunction((room,plate1)), 
           'j62' : con.AreaJunction((room,plate2)),
           'j63' : con.AreaJunction((room,plate3)),
           'j64' : con.AreaJunction((room,plate4)),
           'j65' : con.AreaJunction((room,plate5)),
           'j12' : con.LineJunction((plate1,plate2), Lz, angle_R),
           'j23' : con.LineJunction((plate2,plate3), Lz, angle_R),
           'j34' : con.LineJunction((plate3,plate4), Lz, angle_R),
           'j41' : con.LineJunction((plate4,plate2), Lz, angle_R),
           'j15' : con.LineJunction((plate1,plate5), Lx, angle_R),
           'j25' : con.LineJunction((plate2,plate5), Ly, angle_R),
           'j35' : con.LineJunction((plate3,plate5), Lx, angle_R),
           'j45' : con.LineJunction((plate4,plate5), Ly, angle_R)
           }

# define load
# room power
power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(6,0,dof.DOFtype(typestr = 'power')), name = '1Watt')
# plate power
#power1Watt = lC.Load(omega, np.ones(omega.shape), dof.DOF(2,3,dof.DOFtype(typestr = 'power')), name = '1Watt')


#create SEA model
box = mds.HybridModel((plate1,plate2,plate3,plate4,plate5,room),xdata=omega)

#connect both
box.add_junction(juncs)
box.add_SIF({'sif1' : sif1, 
             'sif2' : sif2, 
             'sif3' : sif3, 
             'sif4' : sif4, 
             'sif5' : sif5})

box.add_load('1Watt',power1Watt) # add 1Watt per band 



#%% Solve
box.create_SEA_matrix(sym = 1)
box.solve()


#%% plotting 3

box.result.plot(2,ID = [1,5],xscale = 'log',yscale = 'log',fulllegstr=('a)','a)'))
plt.figure(2)
plt.yscale('log')
plt.ylim(1e-5,0.01)
plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
#plt.ylabel('$TL/$dB')
plt.legend()
plt.show()

#%% plotting 3
 
box.result.plot(3,ID = [6],xscale = 'log',yscale = 'linear',fulllegstr=('a)',))   

plt.xticks(2*np.pi*fc,fclabels)
plt.xlabel('$f_c/$Hz')
plt.ylabel('$p_{\\rm rms}/$Pa')
plt.tight_layout()
plt.legend()
plt.show()

#%% more info
sif1_in = box.power_input('sif1')
sif2_in = box.power_input('sif2')
sif3_in = box.power_input('sif3')
sif4_in = box.power_input('sif4')
sif5_in = box.power_input('sif5')

sif_all = sif1_in.sum()+sif2_in.sum()+sif3_in.sum()+sif4_in.sum()+sif5_in.sum()


#%% plot 4

sif1_in.plot(4,xscale='log',yscale='linear')
plt.figure(4)
#plt.plot(om_VA1,power_in_1_res,':',label = 'VA1 res')
#plt.plot(om_VA1,power_in_1_nonres,':',label = 'VA1 non-res')
plt.yscale('log')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#%% plot 5

sif_all.plot(5,xscale='log',res = 'dB',fulllegstr=('b)',))
plt.figure(5)
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
#plt.yticks([40,60,80,100,],[40,60,80,100,])
#plt.ylabel('$\Pi_{in}/$W')
plt.legend()
plt.tight_layout()
plt.show()

#plt.savefig('../source/images/box_isolated_SIF_total.png')



#%% plot 6    
NR = 1/sif_all.ydata[0,:]

plt.figure(6)
plt.semilogx(om,10*np.log10(NR),label = 'b)')
plt.xlabel('$f_c/$Hz')
plt.xticks(2*np.pi*fc,fclabels)
plt.ylabel('$IL/$dB')
plt.legend()
plt.tight_layout()
plt.show()