model new

model large-strain off

block create brick 0 8. 0 1.3 0 8.

block cut joint-set dip 0 dip-dir 0 origin 0. 0. 1 spacing 1. number 100

block group "series=B2"
block group "series=B1" range union pos-z 0 1 pos-z 2 3 pos-z 4 5 pos-z 6 7 

block hide range group "series=B2"
block cut joint-set dip 90 dip-dir 90 origin 2 0. 0 spacing 2. number 100
block hide off

block hide range group "series=B1"
block cut joint-set dip 90 dip-dir 90 origin 1 0. 0 spacing 2. number 100

block hide off

block group 'wall' range pos-z 0.0 8.0

block create brick 0.0 8.0 0.0 1.3 8.0 9.0

block group 'beam' range pos-z 8.0 9.0

block zone generate edge 0.5 

block zone cmodel assign elastic range group 'beam'
block zone property density 2.5 young 3.e7 poisson 0.2 range group 'beam'

block zone cmodel assign elastic range group 'wall'
block zone property density 1.8 young 2.e6 poisson 0.2  range group 'wall'

block contact jmodel assign elastic range group 'wall'
block contact property stiffness-normal 2.5e6 stiffness-shear 1.e6 range group 'wall'                       

block contact jmodel assign elastic range group 'beam'
block contact property stiffness-normal 8.e7 stiffness-shear 3.5e7 range group 'beam'

block contact group "horiz" range orientation normal 0. 0. 1

block gridpoint apply vel 0 0 0 range pos-z -0.01 0.01

block face apply stress 0 0 -500 0 0 0 range pos-z 8.99 9.01

block mech damp local 0.8

model hist mechanical unbal-max

model solve

model save 'initial_8x8'

block contact reset disp

;loading at the place of applied velocity

fish def r_x
    loop foreach gx block.gp.list
        if block.gp.id(gx) = 1856 then ; gridpoint of applied velocity
           reac_vec_gp = block.gp.force.reaction(gx)
        endif
    endloop
    r_x = reac_vec_gp(1)
    r_y = reac_vec_gp(2)
    r_z = reac_vec_gp(3)
end

;normal and shear forces in horizontal wall sections
;
fish def shear_1
   local shear_forces = vector(0,0,0)
   normal_force = 0.0
   loop foreach cx block.subcontact.list
        pos = block.subcontact.pos.z(cx)
        if pos = 2 then ; bracketing horizontal section by height
;          (height 3 maintains equilibrium, height 2 not - see Plot02)
        io.out(pos)
           shear_forces = shear_forces + block.subcontact.force.shear(cx) ; sum of section shear forces
           normal_force = normal_force + block.subcontact.force.norm(cx) ; sum of section normal forces
        endif
    endloop
    shear_1=shear_forces(1)
    shear_2=shear_forces(2)
    shear_3=shear_forces(3)

    shear_total=math.mag(shear_forces)
end


;reactions

fish def reac_z
    local reac_vec = vector(0,0,0)
    loop foreach gx block.gp.list
        if block.gp.pos.z(gx) < 0.01 then
           reac_vec = reac_vec + block.gp.force.reaction(gx)
        endif
    endloop
    reac_x = reac_vec(1)
    reac_z = reac_vec(3)
end

;applied velocity

[freq = 1/20]
fish def sin_
  sin_ = math.sin(2.0*math.pi*freq*mech.time) + 2.0*math.pi*freq*mech.time*math.cos(2.0*math.pi*mech.time*freq)
end

block gridpoint apply vel-x 1.e-4  range id 1856;fish sin_ range id 479 

history interval 1000

;block history displacement-x gridpointid 411 ;position 0.01 0.15 2.1

fish hist name 'Rection_z' reac_z
fish hist name 'Reaction_x' reac_x

fish hist name 'Section_Shear_Force_1' shear_1
fish hist name 'Section_Shear_Force_2' shear_2
fish hist name 'Section_Shear_Force_3' shear_3
fish hist name 'Section_Shear_Total' shear_total


fish hist name 'Section_Normal_Force' normal_force

fish hist name 'Loading_x' r_x
fish hist name 'Loading_y' r_y
fish hist name 'Loading_z' r_z

model cycle 500000

model save 'shear_wall_8x8'

program return