SeismicAnalysisClass.f90 Source File


Contents


Source Code

module SeismicAnalysisClass
    use fem
    implicit none

    integer(int32) :: WAVE_DISP = 1
    integer(int32) :: WAVE_VELOCITY = 2
    integer(int32) :: WAVE_ACCEL = 3


    type::SeismicAnalysis_
        type(FEMDomain_),pointer :: femdomain
        real(real64),allocatable :: da(:) ! increment of accel.
        real(real64),allocatable :: a(:) ! accel.
        real(real64),allocatable :: a_ext(:) ! External accel.
        real(real64),allocatable :: a_ext_n(:) ! External accel.
        real(real64),allocatable :: v(:) ! velocity
        real(real64),allocatable :: u(:) ! disp.
        real(real64),allocatable :: du(:) ! increment of disp.
        real(real64),allocatable :: wave(:,:)
        real(real64),allocatable :: dwave(:,:)

        real(real64),allocatable :: Density(:)
        real(real64),allocatable :: YoungModulus(:)
        real(real64),allocatable :: PoissonRatio(:)

        real(real64) :: MaxA(3) = 0.0d0
        real(real64) :: MaxV(3) = 0.0d0
        real(real64) :: MaxU(3) = 0.0d0

        integer(int32),allocatable :: WaveNodeList(:)
        integer(int32),allocatable :: FixNodeList_x(:)
        integer(int32),allocatable :: FixNodeList_y(:)
        integer(int32),allocatable :: FixNodeList_z(:)
        real(real64),allocatable :: FixNodeList_Disp_x(:)
        real(real64),allocatable :: FixNodeList_Disp_y(:)
        real(real64),allocatable :: FixNodeList_Disp_z(:)
        character(1) :: wavedirection="z"
        integer(int32) :: wavetype = 0
        real(real64) :: dt=1.0d0
        real(real64) :: error=0.0d0
        real(real64) :: t=0.0d0
        integer(int32) :: step=0
        real(real64) :: alpha = 0.52400d0
        real(real64) :: beta  = 0.00129d0 ! Rayleigh damping parameters, h=1%
        real(real64) :: Newmark_beta  = 0.250d0 ! Nemark-beta method parameters
        real(real64) :: Newmark_gamma  = 0.50d0 ! Nemark-beta method parameters
        logical :: restart=.False.
    contains
        procedure, public :: init => initSeismicAnalysis
        procedure, public :: loadWave => loadWaveSeismicAnalysis
        procedure, public :: fixDisplacement => fixDisplacementSeismicAnalysis 
        procedure, public :: updateWave => updateWaveSeismicAnalysis
        procedure, public :: run => runSeismicAnalysis
        procedure, public :: LinearReyleighNewmark => LinearReyleighNewmarkSeismicAnalysis
        procedure, public :: recordMaxValues => recordMaxValuesSeismicAnalysis
        procedure, public :: save => saveSeismicAnalysis
        procedure, public :: getNewmarkBetaMatrix => getNewmarkBetaMatrixSeismicAnalysis
        procedure, public :: getNewmarkBetaVector => getNewmarkBetaVectorSeismicAnalysis
        procedure, public :: updateVelocityNewmarkBeta => updateVelocityNewmarkBetaSeismicAnalysis
        procedure, public :: updateAccelNewmarkBeta => updateAccelNewmarkBetaSeismicAnalysis
    end type

contains

! ##############################################
subroutine initSeismicAnalysis(obj)
    class(SeismicAnalysis_),intent(inout) :: obj

    obj%U = zeros(obj%femdomain%nn()*obj%femdomain%nd() )
    obj%V = zeros(obj%femdomain%nn()*obj%femdomain%nd() )
    obj%A = zeros(obj%femdomain%nn()*obj%femdomain%nd() )

    if(obj%femdomain%mesh%empty() )then
        print *, "[ERROR] Seismic % init >> obj%femdomain is empty"
        stop
    endif
    obj%Density = zeros(obj%femdomain%ne() )
    obj%YoungModulus = zeros(obj%femdomain%ne() )
    obj%PoissonRatio = zeros(obj%femdomain%ne() )
    obj%A_ext = zeros(obj%femdomain%nn()*obj%femdomain%nd() )
    obj%A_ext_n = zeros(obj%femdomain%nn()*obj%femdomain%nd() )

end subroutine
! ##############################################


! ##############################################
subroutine saveSeismicAnalysis(obj,name,ratio)
    class(SeismicAnalysis_),intent(inout) :: obj
    character(*),intent(in) :: name
    real(real64),optional,intent(in) :: ratio
    real(real64) :: rat
    integer(int32) :: i,j

    rat = input(default=1.0d0,option=ratio)

    do i=1,obj%femdomain%nn()
        do j=1,obj%femdomain%nd()
            obj%femdomain%mesh%nodcoord(i,j) = obj%femdomain%mesh%nodcoord(i,j)&
                + rat*obj%U( obj%femdomain%nd()*(i-1) + j )
        enddo
    enddo

    call obj%femdomain%msh(name=name)

    do i=1,obj%femdomain%nn()
        do j=1,obj%femdomain%nd()
            obj%femdomain%mesh%nodcoord(i,j) = obj%femdomain%mesh%nodcoord(i,j)&
                - rat*obj%U( obj%femdomain%nd()*(i-1) + j )
        enddo
    enddo
    

    !obj%femdomain%mesh%nodcoord(:,:) = obj%femdomain%mesh%nodcoord(:,:) &
    !    - rat*reshape(obj%U,obj%femdomain%nn(),obj%femdomain%nd() )
    
end subroutine
! ##############################################


! ##############################################
subroutine loadWaveSeismicAnalysis(obj,x_min,x_max,y_min,y_max,z_min,z_max,direction,wavetype)
    class(SeismicAnalysis_),intent(inout) :: obj
    real(real64),optional,intent(in) :: x_min,x_max,y_min,y_max,z_min,z_max
    integer(int32),intent(in) :: wavetype !
    character(1),optional,intent(in) :: direction ! x, y or z
    obj%wavetype= wavetype
    
    if(obj%wavetype < 0 .or. obj%wavetype >3)then
        print *, "Invalid loadAs :: WAVE_DISP,WAVE_VELOCITY or WAVE_ACCEL"
        stop
    endif

    if(present(direction) )then
        obj%wavedirection = direction
    endif

    obj%WaveNodeList = obj%femdomain%select(&
        x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)
    
end subroutine
! ##############################################


! ##############################################
subroutine fixDisplacementSeismicAnalysis(obj,x_min,x_max,y_min,y_max,z_min,z_max,displacement,direction,release)
    class(SeismicAnalysis_),intent(inout) :: obj
    real(real64),optional,intent(in) :: x_min,x_max,y_min,y_max,z_min,z_max,displacement
    character(*),optional,intent(in) :: direction ! x, y or z
    character(*),optional,intent(in) :: release
    real(real64) :: disp 
    real(real64),allocatable :: buf(:)

    disp = input(default=0.0d0,option=displacement)

    if(present(release) )then
        if(index(release,"x")/=0 )then
            deallocate(obj%FixNodeList_x)
            deallocate(obj%FixNodeList_Disp_x)
        endif
        if(index(release,"X")/=0 )then
            deallocate(obj%FixNodeList_x)
            deallocate(obj%FixNodeList_Disp_x)
        endif
        if(index(release,"y")/=0 )then
            deallocate(obj%FixNodeList_y)
            deallocate(obj%FixNodeList_Disp_y)
        endif
        if(index(release,"Y")/=0 )then
            deallocate(obj%FixNodeList_y)
            deallocate(obj%FixNodeList_Disp_y)
        endif
        if(index(release,"z")/=0 )then
            deallocate(obj%FixNodeList_z)
            deallocate(obj%FixNodeList_Disp_z)
        endif
        if(index(release,"Z")/=0 )then
            deallocate(obj%FixNodeList_z)
            deallocate(obj%FixNodeList_Disp_z)
        endif
        if(index(release,"all")/=0 )then
            deallocate(obj%FixNodeList_x)
            deallocate(obj%FixNodeList_Disp_x)
            deallocate(obj%FixNodeList_y)
            deallocate(obj%FixNodeList_Disp_y)
            deallocate(obj%FixNodeList_z)
            deallocate(obj%FixNodeList_Disp_z)
        endif
        return
    endif

    if(present(direction) )then
        if( trim(direction) == "x" .or.  trim(direction) == "X")then
            obj%FixNodeList_x = hstack(obj%FixNodeList_x , &
             obj%femdomain%select(&
                x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
                )
            buf = zeros(size(obj%FixNodeList_x))
            buf(:) = disp
            if(.not. allocated(obj%FixNodeList_Disp_x) )then
                obj%FixNodeList_Disp_x = buf
            else
                buf(1:size(obj%FixNodeList_Disp_x) ) = obj%FixNodeList_Disp_x(:)
            endif
            obj%FixNodeList_Disp_x = buf

        elseif( trim(direction) == "y" .or.  trim(direction) == "Y")then
            obj%FixNodeList_y = hstack(obj%FixNodeList_y , &
             obj%femdomain%select(&
                x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
                )
                buf = zeros(size(obj%FixNodeList_y))
                buf(:) = disp
                if(.not. allocated(obj%FixNodeList_Disp_y) )then
                else
                    buf(1:size(obj%FixNodeList_Disp_y) ) = obj%FixNodeList_Disp_y(:)
                endif
                obj%FixNodeList_Disp_y = buf

        elseif( trim(direction) == "z" .or.  trim(direction) == "Z")then
            obj%FixNodeList_z = hstack(obj%FixNodeList_z , &
             obj%femdomain%select(&
                x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
                )
                buf = zeros(size(obj%FixNodeList_z))
                buf(:) = disp

            if(.not. allocated(obj%FixNodeList_Disp_z) )then
            else
                buf(1:size(obj%FixNodeList_Disp_z) ) = obj%FixNodeList_Disp_z(:)
            endif
            obj%FixNodeList_Disp_z = buf
    
        elseif( trim(direction) == "all" .or.  trim(direction) == "ALL")then
        
            obj%FixNodeList_x = hstack(obj%FixNodeList_x , &
             obj%femdomain%select(&
            x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
            )
            buf = zeros(size(obj%FixNodeList_x))
            buf(:) = disp

            if(.not. allocated(obj%FixNodeList_Disp_x) )then
            else
                buf(1:size(obj%FixNodeList_Disp_x) ) = obj%FixNodeList_Disp_x(:)
            endif
            obj%FixNodeList_Disp_x = buf
        
            obj%FixNodeList_y = hstack(obj%FixNodeList_y , &
             obj%femdomain%select(&
                x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
                )
                buf = size(zeros(size(obj%FixNodeList_y)))
                buf(:) = disp

                if(.not. allocated(obj%FixNodeList_Disp_y) )then
                else
                    buf(1:size(obj%FixNodeList_Disp_y) ) = obj%FixNodeList_Disp_y(:)
                endif
                obj%FixNodeList_Disp_y = buf    

            obj%FixNodeList_z = hstack(obj%FixNodeList_z , &
             obj%femdomain%select(&
                x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
                )
                buf = zeros(size(obj%FixNodeList_z))
            buf(:) = disp

            if(.not. allocated(obj%FixNodeList_Disp_z) )then
            else
                buf(1:size(obj%FixNodeList_Disp_z) ) = obj%FixNodeList_Disp_z(:)
            endif
            obj%FixNodeList_Disp_z = buf
        else
            print *, "ERROR :: loadWaveSeismicAnalysis >> direction should be x, y or z"
            stop 
        endif
    else
        obj%FixNodeList_x = hstack(obj%FixNodeList_x , &
         obj%femdomain%select(&
            x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
            )
            buf = zeros(size(obj%FixNodeList_x))
            buf(:) = disp
        
            if(.not. allocated(obj%FixNodeList_Disp_x) )then
            else
                buf(1:size(obj%FixNodeList_Disp_x) ) = obj%FixNodeList_Disp_x(:)
            endif
            obj%FixNodeList_Disp_x = buf

        obj%FixNodeList_y = hstack(obj%FixNodeList_y , &
         obj%femdomain%select(&
            x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
            )
            buf = size(zeros(size(obj%FixNodeList_y)))
            buf(:) = disp

            if(.not. allocated(obj%FixNodeList_Disp_y) )then
            else
                buf(1:size(obj%FixNodeList_Disp_y) ) = obj%FixNodeList_Disp_y(:)
            endif
            obj%FixNodeList_Disp_y = buf

        obj%FixNodeList_z = hstack(obj%FixNodeList_z , &
         obj%femdomain%select(&
            x_min=x_min,x_max=x_max,y_min=y_min,y_max=y_max,z_min=z_min,z_max=z_max)&
            )
            buf = zeros(size(obj%FixNodeList_z))
            buf(:) = disp
            
            if(.not. allocated(obj%FixNodeList_Disp_z) )then
            else
                buf(1:size(obj%FixNodeList_Disp_z) ) = obj%FixNodeList_Disp_z(:)
            endif
            obj%FixNodeList_Disp_z = buf

    endif
end subroutine
! ##############################################



subroutine updateWaveSeismicAnalysis(obj,timestep,direction)
    class(SeismicAnalysis_),intent(inout) :: obj
    integer(int32),intent(in) :: timestep
    character(1),optional,intent(in) :: direction ! x, y or z
    integer(int32) :: node_id,i,dir,dim_num

    if(present(direction) )then
        obj%wavedirection = direction
    endif
    if(obj%wavedirection=="x" .or.obj%wavedirection=="X" )then
        dir=1
    endif
    if(obj%wavedirection=="y" .or.obj%wavedirection=="Y" )then
        dir=2
    endif
    if(obj%wavedirection=="z" .or.obj%wavedirection=="Z" )then
        dir=3
    endif

    if(.not. allocated(obj%WaveNodeList) )then
        print *, "Caution >> updateWaveSeismicAnalysis >> no wave"
    endif
    dim_num = obj%femdomain%nd()
    if(obj%wavetype==WAVE_ACCEL)then
        obj%A_ext_n = obj%A_ext
        do i=1,size(obj%WaveNodeList)
            node_id = obj%WaveNodeList(i)
            !obj%A( dim_num*(node_id-1)+dir ) = obj%wave(timestep ,2)
            ! update accel
            obj%A_ext( dim_num*(node_id-1)+dir ) =  obj%wave(timestep ,2)
        enddo
        
    elseif(obj%wavetype==WAVE_VELOCITY)then
        do i=1,size(obj%WaveNodeList)
            node_id = obj%WaveNodeList(i)
            obj%V( dim_num*(node_id-1)+dir ) = obj%wave(timestep ,2)
        enddo
    elseif(obj%wavetype==WAVE_DISP)then
        do i=1,size(obj%WaveNodeList)
            node_id = obj%WaveNodeList(i)
            obj%U( dim_num*(node_id-1)+dir ) = obj%wave(timestep ,2)
        enddo
    endif

end subroutine

! ##############################################
subroutine runSeismicAnalysis(obj,t0,timestep,wave,AccelLimit,disp_magnify_ratio)
    class(SeismicAnalysis_),intent(inout) :: obj
    integer(int32),intent(in) :: timestep(2)

    type(LinearSolver_) :: solver
    type(IO_) :: U, V, A
    real(real64),optional,intent(in) :: t0,disp_magnify_ratio
    real(real64),optional,intent(in) :: wave(:,:),AccelLimit
    integer(int32) :: i,j
    real(real64) :: ratio

    ratio = input(default=1.0d0,option=disp_magnify_ratio)
    if(present(wave) )then
        obj%wave = wave
    endif

    do i=timestep(1),timestep(2)-1
        ! update dt
        obj%dt = abs(obj%wave(i+1,1) - obj%wave(i,1))
        
        ! update time
        obj%step = i
        obj%t = obj%dt*obj%Step
        call obj%updateWave(timestep=obj%step+1)
        ! show info.
        call print("SeismicAnalysis >> "//str(obj%t-obj%dt)//"< t <"//str(obj%t)//" sec.")
        
        ! solve Linear-ElastoDynamic problem with Reyleigh dumping and Newmark Beta
        if(present(AccelLimit) )then
            if(maxval(obj%A)>=AccelLimit)then
                print *, "[Caution] :: runSeismicAnalysis >> exceeds AccelLimit!"
                return    
            endif
        endif

        call obj%LinearReyleighNewmark()

        call obj%recordMaxValues()
        ! Export results
        call obj%save("step_"//str(obj%step),ratio=ratio)

    enddo 

end subroutine
! ##############################################


! ##############################################
subroutine LinearReyleighNewmarkSeismicAnalysis(obj,TOL)
    class(SeismicAnalysis_),intent(inout) :: obj
    type(LinearSolver_) :: solver
    type(IO_) :: f
    real(real64),allocatable :: M_ij(:,:)
    real(real64),allocatable :: C_ij(:,:)
    real(real64),allocatable :: K_ij(:,:)
    real(real64),allocatable :: F_i(:)
    real(real64),allocatable :: dF_i(:)
    real(real64),allocatable :: R_i(:)
    real(real64),allocatable :: U_i(:)
    real(real64),allocatable ::dU_i(:)
    real(real64),allocatable :: V_i(:)
    real(real64),allocatable ::dV_i(:)
    real(real64),allocatable :: A_i(:)
    real(real64),allocatable :: A_ext_i(:)
    real(real64),allocatable :: A_ext_i_n(:)
    real(real64),allocatable ::dA(:)
    real(real64),allocatable ::dV(:)
    real(real64),allocatable ::dU(:)

    real(real64),allocatable ::u_upd(:)
    real(real64),allocatable ::v_upd(:)
    real(real64),allocatable ::a_upd(:)


    real(real64),allocatable ::U_n(:)
    real(real64),allocatable ::V_n(:)
    real(real64),allocatable ::A_n(:)

    real(real64),allocatable :: A_ij(:,:)
    integer(int32) :: i,j,k,l,m,dim_num,n
    integer(int32),allocatable :: FixNodeList(:) ,DomainIDs(:)
    real(real64),allocatable   :: Coordinate(:,:)
    real(real64),optional,intent(in) :: TOL
    real(real64) :: TOL_seismic,center_accel(3),rho,gravity(3),a(0:7)


    gravity(:)=0.0d0
    gravity(3) = -9.81d0
    
    TOL_seismic = input(default=dble(1.0e-14),option=TOL )

    dim_num = size(obj%femdomain%mesh%nodcoord,2)
    n = dim_num * size(obj%femdomain%mesh%nodcoord,1)
    if(.not. allocated(obj%a) )then
        allocate(obj%a(n) )
    endif
    if(.not. allocated(obj%v) )then
        allocate(obj%v(n) )
    endif
    if(.not. allocated(obj%u) )then
        allocate(obj%u(n) )
    endif

        ! Element matrix
        call solver%init(NumberOfNode=[obj%femdomain%nn()],DOF=3)
        obj%dt = abs(obj%dt)
        do i=1,obj%femdomain%ne()
            ! For each element
            ! Ax=b will be installed into solv -obj%A_ext_n)
            rho = obj%Density(i)
            M_ij = obj%femdomain%MassMatrix(ElementID=i,DOF=obj%femdomain%nd(),density=rho )
            K_ij = obj%femdomain%StiffnessMatrix(ElementID=i,E=obj%YoungModulus(i),v=obj%PoissonRatio(i) )
            C_ij = obj%alpha * M_ij + obj%beta * K_ij
            U_i  = obj%femdomain%ElementVector(ElementID=i, GlobalVector=obj%U, DOF=obj%femdomain%nd() )
            V_i  = obj%femdomain%ElementVector(ElementID=i, GlobalVector=obj%V, DOF=obj%femdomain%nd() )
            A_i  = obj%femdomain%ElementVector(ElementID=i, GlobalVector=obj%A, DOF=obj%femdomain%nd() )
            A_ext_i  = obj%femdomain%ElementVector(ElementID=i, GlobalVector=obj%A_ext, DOF=obj%femdomain%nd() )
            A_ext_i_n  = obj%femdomain%ElementVector(ElementID=i, GlobalVector=obj%A_ext_n, DOF=obj%femdomain%nd() )
            
            dim_num = obj%femdomain%nd()
            
            F_i = matmul(M_ij,A_ext_i)
            
            A_ij = obj%getNewmarkBetaMatrix(M=M_ij, C=C_ij,K=K_ij,&
                beta=obj%Newmark_beta,gamma=obj%Newmark_gamma,dt=obj%dt)
            R_i = obj%getNewmarkBetaVector(M=M_ij,C=C_ij,u_n=U_i,v_n=V_i, a_n=A_i, &
                force=F_i,beta=obj%Newmark_beta,gamma=obj%Newmark_gamma,dt=obj%dt)

            !print *, maxval(A_ij),minval(A_ij)
            !print *, maxval(R_i),minval(R_i)
            !stop
            ! Assemble stiffness matrix
            DomainIDs = zeros(obj%femdomain%nne() )
            DomainIDs(:) = 1
            
            call solver%assemble(&
                connectivity=obj%femdomain%connectivity(ElementID=i), &
                DOF = obj%femdomain%nd(), &
                DomainIDs=DomainIDs,&
                eMatrix = A_ij)
            call solver%assemble(&
                connectivity=obj%femdomain%connectivity(ElementID=i), &
                DOF = obj%femdomain%nd(), &
                DomainIDs=DomainIDs,&
                eVector = R_i)
        enddo


        if(allocated(obj%FixNodeList_x) )then
            do i=1,size(obj%FixNodeList_x)
                call solver%fix( NodeID=obj%FixNodeList_x(i)*3-2,&
                entryvalue=obj%FixNodeList_Disp_x(i) ,&
                row_DomainID=1)
            enddo
        endif
        if(allocated(obj%FixNodeList_y) )then
            do i=1,size(obj%FixNodeList_y)
                call solver%fix( NodeID=obj%FixNodeList_y(i)*3-1,&
                entryvalue=obj%FixNodeList_Disp_y(i),&
                row_DomainID=1)
            enddo
        endif
        if(allocated(obj%FixNodeList_z) )then
            do i=1,size(obj%FixNodeList_z)
                call solver%fix( NodeID=obj%FixNodeList_z(i)*3,&
                entryvalue=obj%FixNodeList_Disp_z(i),&
                row_DomainID=1)
            enddo
        endif

        ! Now [A] {du} = {R} is ready
        ! Solve
        
        
        call solver%solve("BiCGSTAB")
        print *, maxval(solver%val),minval(solver%val)
        print *, maxval(solver%x),minval(solver%x)
        
        u_upd = solver%x
        v_upd = obj%updateVelocityNewmarkBeta(u=u_upd,u_n=obj%U,v_n=obj%V,a_n=obj%A,&
            gamma=obj%newmark_gamma,beta=obj%newmark_beta,dt=obj%dt)
        a_upd = obj%updateAccelNewmarkBeta(u=u_upd,u_n=obj%U,v_n=obj%V,a_n=obj%A,&
            gamma=obj%newmark_gamma,beta=obj%newmark_beta,dt=obj%dt)
        
        obj%U = u_upd
        obj%V = v_upd
        obj%A = a_upd
       

        print *, "U"
        print *, minval(obj%U),maxval(obj%U)
        
        print *, "V"
        print *, minval(obj%V),maxval(obj%V)
        
        print *, "A"
        print *, minval(obj%A), maxval(obj%A)
    
end subroutine
! ##############################################


subroutine recordMaxValuesSeismicAnalysis(obj)
    class(SeismicAnalysis_),intent(inout) :: Obj
    real(real64),allocatable :: array(:,:)
    integer(int32) :: i
    
    array = reshape( obj%U,obj%femdomain%nn(),obj%femdomain%nd() )
    array = abs(array)
    do i=1,3
        obj%maxU(i) = maxval( [ abs(obj%maxU(i)) , maxval(array(:,i) ) ]  )
    enddo
    

    array = reshape( obj%V,obj%femdomain%nn(),obj%femdomain%nd() )
    array = abs(array)
    do i=1,3
        obj%maxV(i) = maxval( [ abs(obj%maxV(i)) , maxval(array(:,i) ) ]  )
    enddo


    array = reshape( obj%A,obj%femdomain%nn(),obj%femdomain%nd() )
    array = abs(array)
    do i=1,3
        obj%maxA(i) = maxval( [ abs(obj%maxA(i)) , maxval(array(:,i) ) ]  )
    enddo

end subroutine

function getNewmarkBetaMatrixSeismicAnalysis(obj,M,C,K,beta,gamma,dt) result(ret)
    class(SeismicAnalysis_),intent(in) :: obj
    real(real64),intent(in) :: M(:,:),C(:,:),K(:,:),beta,gamma,dt
    real(real64),allocatable :: ret(:,:)
    integer(int32) :: n

    n = size(M,1)
    ret = zeros(n,n)
    ret(:,:) = 1.0d0/( beta*dt*dt ) * M(:,:) &
        + gamma/( beta*dt ) * C(:,:) &
        + K(:,:)

end function

function getNewmarkBetaVectorSeismicAnalysis(obj,u_n, v_n, a_n, force,beta,gamma,M,C,dt) result(ret)
    class(SeismicAnalysis_),intent(in) :: obj
    real(real64),intent(in) :: u_n(:), v_n(:), a_n(:),force(:),beta,gamma,dt
    real(real64),intent(in) :: M(:,:),C(:,:)
    real(real64),allocatable :: ret(:)
    real(real64) :: a(8)
    integer(int32) :: n

    n = size(U_n,1)
    ret = zeros(n)
    a (1) = gamma/(beta*dt)
    a (2) = - gamma/(beta*dt)
    a (3) = 1.0d0 - gamma/beta
    a (4) = dt*(1.0d0 - gamma/(2.0d0*beta) )
    a (5) = 1.0d0/(beta*dt*dt)
    a (6) = - 1.0d0/(beta*dt*dt)
    a (7) = - 1.0d0/(beta*dt)
    a (8) = 1.0d0-1.0d0/(2.0d0*beta)
    ret(:) = force(:) &
        - a(6)*matmul(M,u_n) - a(7)*matmul(M,v_n) - a(8)*matmul(M,a_n) &
        - a(2)*matmul(C,u_n) - a(3)*matmul(C,v_n) - a(4)*matmul(C,a_n)
end function
! ##########################################################################################
function updateVelocityNewmarkBetaSeismicAnalysis(obj,u,u_n,v_n,a_n,gamma,beta,dt) result(ret)
    class(SeismicAnalysis_),intent(in) :: obj
    real(real64),intent(in) :: u(:), u_n(:), v_n(:), a_n(:),beta,gamma,dt
    real(real64),allocatable :: ret(:)
    integer(int32) :: n

    ret = zeros(size(u) )
    ! usually, gamma=0.5, beta=0.25
    !ret(:) = gamma/(beta*dt)*u(:) &
    !    - gamma/(beta*dt)*u_n(:) &
    !    + (1.0d0 - gamma/beta )*v_n(:) &
    !    + dt*(1.0d0 - gamma/(2.0d0*beta) )*a_n(:)

    ret(:) = gamma/(beta*dt)*u(:) & ! modified
        - gamma/(beta*dt)*u_n(:) &  ! modified
        + (1.0d0 - gamma/beta )*v_n(:) &
        + dt*(1.0d0 - gamma/(2.0d0*beta) )*a_n(:)
    ! confirmed 2021/10/13
end function
! ##########################################################################################


! ##########################################################################################
function updateAccelNewmarkBetaSeismicAnalysis(obj,u,u_n,v_n,a_n,gamma,beta,dt) result(ret)
    class(SeismicAnalysis_),intent(in) :: obj
    real(real64),intent(in) :: u(:), u_n(:), v_n(:), a_n(:),beta,gamma,dt
    real(real64),allocatable :: ret(:)
    integer(int32) :: n

    ret = zeros(size(u) )

    !ret(:) = 1.0d0/(beta*dt*dt)*u(:) &
    !    - 1.0d0/(beta*dt*dt)*u_n(:) &
    !    - 1.0d0/(beta*dt)*v_n(:) &
    !    + (1.0d0 - 1.0d0/(2.0d0*beta) )*a_n(:)
    ! confirmed 2021/10/13
    ret(:) = 1.0d0/(beta*dt*dt)*u(:) &
        - 1.0d0/(beta*dt*dt)*u_n(:) &
        - 1.0d0/(beta*dt)*v_n(:) &
        + (1.0d0 - 1.0d0/(2.0d0*beta) )*a_n(:)
end function
! ##########################################################################################

end module SeismicAnalysisClass