SDOF_SFS_MODEL.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
subroutine sdof_sfs_model(sID, gr_acc, direction)
 
  use speed_sci
 
  implicit none
 
  integer*4 :: sID, direction
  real*8, dimension(3) :: gr_acc
  real*8, dimension(4) :: tn, tn1, tn2, tn3, tn4, temp, du, df, vec, a_old, v_old, u_old, f_old
 
  tn = 0; tn1 = 0; tn2 = 0; tn3 = 0; tn4 = 0; temp = 0; du = 0; df = 0; vec = 0
  a_old = sys(sid)%a(:,direction); v_old = sys(sid)%v(:,direction); u_old = sys(sid)%u(:,direction); f_old = sys(sid)%f(:,direction)
 
  tn1(1) = sys(sid)%Ms(1,1)*gr_acc(direction)
  tn1(2) = sys(sid)%Mf*gr_acc(direction)
  tn1(4) = (sys(sid)%Ms(1,1) + sys(sid)%Mf)*gr_acc(3)
 
  temp = (1 - 2*sys(sid)%beta_newmark)/(2*sys(sid)%beta_newmark)*a_old + 1/(sys(sid)%beta_newmark*sys(sid)%dt)*v_old
  tn2 = matmul(sys(sid)%MAT_M,temp)
 
  temp = (sys(sid)%gamma_newmark - 2*sys(sid)%beta_newmark)/(2*sys(sid)%beta_newmark)*sys(sid)%dt*a_old + &
          (sys(sid)%gamma_newmark - sys(sid)%beta_newmark)/(sys(sid)%beta_newmark)*v_old
  tn3 = matmul(sys(sid)%MAT_C,temp)
 
  tn4 = matmul(sys(sid)%MAT_KS,u_old)
 
  tn = tn1 + tn2 + tn3 - tn4 - f_old
 
  du = matmul(sys(sid)%MAT_MCinv,tn)
 
  df = matmul(sys(sid)%MAT_KI,du)
 
  sys(sid)%u(:,direction) = u_old + du
 
  sys(sid)%f(:,direction) = f_old + df
 
  vec = matmul(sys(sid)%MAT_KS,sys(sid)%u(:,direction))
  sys(sid)%fs(direction) = vec(1)
  sys(sid)%fb(direction) = sys(sid)%f(2,direction) + sys(sid)%Mf*(gr_acc(direction) - a_old(2))
 
  sys(sid)%v(:,direction) = v_old + (2*sys(sid)%beta_newmark - sys(sid)%gamma_newmark)/(2*sys(sid)%beta_newmark)*sys(sid)%dt*a_old + &
    sys(sid)%gamma_newmark/(sys(sid)%beta_newmark*sys(sid)%dt)*(du - v_old*sys(sid)%dt)
 
  sys(sid)%a(:,direction) = (2*sys(sid)%beta_newmark - 1)/(2*sys(sid)%beta_newmark)*a_old + &
    1/(sys(sid)%beta_newmark*sys(sid)%dt2)*(du - v_old*sys(sid)%dt)
 
  sys(sid)%IntForce(1,direction) = sys(sid)%f(2,direction)
 
  ! UnComment this after verifying
  !sys(sID)%SDOFItF(3) = sys(sID)%f(4,direction)
 
  return
end subroutine sdof_sfs_model