hecmw_mat_ass.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_matrix_ass
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  implicit none
 
  private
 
  public :: hecmw_mat_ass_elem
  public :: hecmw_mat_add_node
  public :: hecmw_mat_profile_has_node
  public :: hecmw_array_search_i
  public :: hecmw_mat_ass_equation
  public :: hecmw_mat_ass_equation_rhs
  public :: hecmw_mat_add_dof
  public :: hecmw_mat_ass_bc
  public :: hecmw_mat_ass_bc_contactlag
  public :: hecmw_mat_ass_contactlag
  public :: stf_get_block
 
contains
 
  !C
  !C***
  !C*** MAT_ASS_ELEM
  !C***
  !C
  subroutine hecmw_mat_ass_elem(hecMAT, nn, nodLOCAL, stiffness)
    type (hecmwst_matrix)     :: hecmat
    integer(kind=kint) :: nn
    integer(kind=kint) :: nodlocal(:)
    real(kind=kreal) :: stiffness(:, :)
    !** Local variables
    integer(kind=kint) :: ndof, inod_e, jnod_e, inod, jnod
    real(kind=kreal) :: a(6,6)
 
    ndof = hecmat%NDOF
 
    do inod_e = 1, nn
      inod = nodlocal(inod_e)
      do jnod_e = 1, nn
        jnod = nodlocal(jnod_e)
        !***** Add components
        call stf_get_block(stiffness, ndof, inod_e, jnod_e, a)
        call hecmw_mat_add_node(hecmat, inod, jnod, a)
      enddo
    enddo
 
  end subroutine hecmw_mat_ass_elem
 
  subroutine stf_get_block(stiffness, ndof, inod, jnod, a)
    real(kind=kreal) :: stiffness(:, :), a(:, :)
    integer(kind=kint) :: ndof, inod, jnod
    !** Local variables
    integer(kind=kint) :: row_offset, col_offset, i, j
 
    row_offset = ndof*(inod-1)
    do i = 1, ndof
 
      col_offset = ndof*(jnod-1)
      do j = 1, ndof
 
        a(i, j) = stiffness(i + row_offset, j + col_offset)
      enddo
    enddo
  end subroutine stf_get_block
 
 
  subroutine hecmw_mat_add_node(hecMAT, inod, jnod, a)
    type (hecmwst_matrix) :: hecmat
    integer(kind=kint) :: inod, jnod
    real(kind=kreal) :: a(:, :)
    !** Local variables
    integer(kind=kint) :: ndof, is, ie, k, idx_base, idx, idof, jdof
 
    ndof = hecmat%NDOF
 
    if (inod < jnod) then
      is = hecmat%indexU(inod-1)+1
      ie = hecmat%indexU(inod)
      k = hecmw_array_search_i(hecmat%itemU, is, ie, jnod)
 
      if (k < is .or. ie < k) then
        !$omp critical
        write(*,*) '###ERROR### : cannot find connectivity (1)'
        write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
        !$omp end critical
        call hecmw_abort(hecmw_comm_get_comm())
      endif
 
      idx_base = ndof**2 * (k-1)
      do idof = 1, ndof
        do jdof = 1, ndof
          idx = idx_base + jdof
          !$omp atomic
          hecmat%AU(idx) = hecmat%AU(idx) + a(idof, jdof)
        enddo
        idx_base = idx_base + ndof
      enddo
 
    else if (inod > jnod) then
      is = hecmat%indexL(inod-1)+1
      ie = hecmat%indexL(inod)
      k = hecmw_array_search_i(hecmat%itemL, is, ie, jnod)
 
      if (k < is .or. ie < k) then
        !$omp critical
        write(*,*) '###ERROR### : cannot find connectivity (2)'
        write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
        !$omp end critical
        call hecmw_abort(hecmw_comm_get_comm())
      endif
 
      idx_base = ndof**2 * (k-1)
      do idof = 1, ndof
        do jdof = 1, ndof
          idx = idx_base + jdof
          !$omp atomic
          hecmat%AL(idx) = hecmat%AL(idx) + a(idof, jdof)
        enddo
        idx_base = idx_base + ndof
      enddo
 
    else
      idx_base = ndof**2 * (inod - 1)
      do idof = 1, ndof
        do jdof = 1, ndof
          idx = idx_base + jdof
          !$omp atomic
          hecmat%D(idx) = hecmat%D(idx) + a(idof, jdof)
        enddo
        idx_base = idx_base + ndof
      enddo
    endif
  end subroutine hecmw_mat_add_node
 
 
  function hecmw_mat_profile_has_node(hecMAT, inod, jnod)
    logical :: hecmw_mat_profile_has_node
    type (hecmwst_matrix) :: hecmat
    integer(kind=kint) :: inod, jnod
    !** Local variables
    integer(kind=kint) :: is, ie, k
 
    if (inod < jnod) then
      is = hecmat%indexU(inod-1)+1
      ie = hecmat%indexU(inod)
      k = hecmw_array_search_i(hecmat%itemU, is, ie, jnod)
      hecmw_mat_profile_has_node = ( is <= k .and. k <= ie )
    else if (inod > jnod) then
      is = hecmat%indexL(inod-1)+1
      ie = hecmat%indexL(inod)
      k = hecmw_array_search_i(hecmat%itemL, is, ie, jnod)
      hecmw_mat_profile_has_node = ( is <= k .and. k <= ie )
    else
      hecmw_mat_profile_has_node = .true.
    endif
  end function hecmw_mat_profile_has_node
 
 
  function hecmw_array_search_i(array, is, iE, ival)
    integer(kind=kint) :: hecmw_array_search_i
    integer(kind=kint) :: array(:)
    integer(kind=kint) :: is, ie, ival
    !** Local variables
    integer(kind=kint) :: left, right, center, cval
 
    left = is
    right = ie
    do
      if (left > right) then
        center = -1
        exit
      endif
 
      center = (left + right) / 2
      cval = array(center)
 
      if (ival < cval) then
        right = center - 1
      else if (cval < ival) then
        left = center + 1
      else
        exit
      endif
    enddo
 
    hecmw_array_search_i = center
 
  end function hecmw_array_search_i
 
 
  !C
  !C***
  !C*** MAT_ASS_EQUATION
  !C***
  !C
  subroutine hecmw_mat_ass_equation ( hecMESH, hecMAT )
    type (hecmwst_matrix), target :: hecmat
    type (hecmwst_local_mesh)     :: hecmesh
    !** Local variables
    real(kind=kreal), pointer :: penalty
    real(kind=kreal) :: alpha, a1_2inv, ai, aj, factor
    integer(kind=kint) :: impc, is, ie, i, j, inod, idof, jnod, jdof
    logical :: is_internal_i
 
    if( hecmw_mat_get_penalized(hecmat) == 1 ) return
 
    ! write(*,*) "INFO: imposing MPC by penalty"
 
    penalty => hecmat%Rarray(11)
 
    if (penalty < 0.0) stop "ERROR: negative penalty"
    if (penalty < 1.0) write(*,*) "WARNING: penalty ", penalty, " smaller than 1"
 
    alpha= hecmw_mat_diag_max(hecmat, hecmesh) * penalty
    call hecmw_mat_set_penalty_alpha(hecmat, alpha)
 
    outer: do impc = 1, hecmesh%mpc%n_mpc
      is = hecmesh%mpc%mpc_index(impc-1) + 1
      ie = hecmesh%mpc%mpc_index(impc)
 
      do i = is, ie
        if (hecmesh%mpc%mpc_dof(i) > hecmat%NDOF) cycle outer
      enddo
 
      a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
 
 
      do i = is, ie
        inod = hecmesh%mpc%mpc_item(i)
 
        is_internal_i = (hecmesh%node_ID(2*inod) == hecmw_comm_get_rank())
        if (.not. is_internal_i) cycle
 
        idof = hecmesh%mpc%mpc_dof(i)
        ai = hecmesh%mpc%mpc_val(i)
        factor = ai * a1_2inv
 
        do j = is, ie
          jnod = hecmesh%mpc%mpc_item(j)
          jdof = hecmesh%mpc%mpc_dof(j)
          aj = hecmesh%mpc%mpc_val(j)
 
          call hecmw_mat_add_dof(hecmat, inod, idof, jnod, jdof, aj*factor*alpha)
        enddo
 
      enddo
    enddo outer
 
    call hecmw_mat_set_penalized(hecmat, 1)
 
  end subroutine hecmw_mat_ass_equation
 
 
  subroutine hecmw_mat_ass_equation_rhs ( hecMESH, hecMAT )
    type (hecmwst_matrix), target :: hecmat
    type (hecmwst_local_mesh)     :: hecmesh
    !** Local variables
    real(kind=kreal) :: alpha, a1_2inv, ai, factor, ci
    integer(kind=kint) :: ndof, impc, is, ie, i, inod, idof
 
    if( hecmw_mat_get_penalized_b(hecmat) == 1) return
 
    alpha = hecmw_mat_get_penalty_alpha(hecmat)
    if (alpha <= 0.0) stop "ERROR: penalty applied on vector before matrix"
 
    ndof = hecmat%NDOF
 
    outer: do impc = 1, hecmesh%mpc%n_mpc
      is = hecmesh%mpc%mpc_index(impc-1) + 1
      ie = hecmesh%mpc%mpc_index(impc)
 
      do i = is, ie
        if (hecmesh%mpc%mpc_dof(i) > ndof) cycle outer
      enddo
 
      a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
 
      do i = is, ie
        inod = hecmesh%mpc%mpc_item(i)
 
        idof = hecmesh%mpc%mpc_dof(i)
        ai = hecmesh%mpc%mpc_val(i)
        factor = ai * a1_2inv
 
        ci = hecmesh%mpc%mpc_const(impc)
        !$omp atomic
        hecmat%B(ndof*(inod-1)+idof) = hecmat%B(ndof*(inod-1)+idof) + ci*factor*alpha
      enddo
    enddo outer
 
    call hecmw_mat_set_penalized_b(hecmat, 1)
 
  end subroutine hecmw_mat_ass_equation_rhs
 
 
  subroutine hecmw_mat_add_dof(hecMAT, inod, idof, jnod, jdof, val)
    type (hecmwst_matrix) :: hecmat
    integer(kind=kint) :: inod, idof, jnod, jdof
    real(kind=kreal) :: val
    !** Local variables
    integer(kind=kint) :: ndof, is, ie, k, idx
 
    ndof = hecmat%NDOF
    if (inod < jnod) then
      is = hecmat%indexU(inod-1)+1
      ie = hecmat%indexU(inod)
      k = hecmw_array_search_i(hecmat%itemU, is, ie, jnod)
 
      if (k < is .or. ie < k) then
        write(*,*) '###ERROR### : cannot find connectivity (3)'
        write(*,*) '  myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
        call hecmw_abort(hecmw_comm_get_comm())
        return
      endif
 
      idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
      !$omp atomic
      hecmat%AU(idx) = hecmat%AU(idx) + val
 
    else if (inod > jnod) then
      is = hecmat%indexL(inod-1)+1
      ie = hecmat%indexL(inod)
      k = hecmw_array_search_i(hecmat%itemL, is, ie, jnod)
 
      if (k < is .or. ie < k) then
        write(*,*) '###ERROR### : cannot find connectivity (4)'
        write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
        call hecmw_abort(hecmw_comm_get_comm())
        return
      endif
 
      idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
      !$omp atomic
      hecmat%AL(idx) = hecmat%AL(idx) + val
 
    else
      idx = ndof**2 * (inod - 1) + ndof * (idof - 1) + jdof
      !$omp atomic
      hecmat%D(idx) = hecmat%D(idx) + val
    endif
 
  end subroutine hecmw_mat_add_dof
 
  !C
  !C***
  !C*** MAT_ASS_BC
  !C***
  !C
  subroutine hecmw_mat_ass_bc(hecMAT, inode, idof, RHS, conMAT)
    type (hecmwst_matrix)     :: hecmat
    integer(kind=kint) :: inode, idof
    real(kind=kreal) :: rhs, val
    type (hecmwst_matrix),optional     :: conmat
    integer(kind=kint) :: ndof, in, i, ii, iii, ndof2, k, is, ie, iis, iie, ik, idx
 
    ndof = hecmat%NDOF
    if( ndof < idof ) return
 
    !C-- DIAGONAL block
 
    hecmat%B(ndof*inode-(ndof-idof)) = rhs
    if(present(conmat)) conmat%B(ndof*inode-(ndof-idof)) = 0.0d0
    ndof2 = ndof*ndof
    ii  = ndof2 - idof
 
    do i = ndof-1,0,-1
      if( i .NE. ndof-idof ) then
        idx = ndof*inode-i
        val = hecmat%D(ndof2*inode-ii)*rhs
        !$omp atomic
        hecmat%B(idx) = hecmat%B(idx) - val
        if(present(conmat)) then
          val = conmat%D(ndof2*inode-ii)*rhs
          !$omp atomic
          conmat%B(idx) = conmat%B(idx) - val
        endif
      endif
      ii = ii - ndof
    end do
 
    !*Set diagonal row to zero
    ii  = ndof2-1 - (idof-1)*ndof
 
    do i = 0, ndof - 1
      hecmat%D(ndof2*inode-ii+i)=0.d0
      if(present(conmat)) conmat%D(ndof2*inode-ii+i)=0.d0
    end do
 
    !*Set diagonal column to zero
    ii = ndof2 - idof
    do i = 1, ndof
      if( i.NE.idof ) then
        hecmat%D(ndof2*inode-ii) = 0.d0
        if(present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
      else
        hecmat%D(ndof2*inode-ii) = 1.d0
        if(present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
      endif
      ii = ii - ndof
    end do
 
    !C-- OFF-DIAGONAL blocks
 
    ii  = ndof2-1 - (idof-1)*ndof
    is = hecmat%indexL(inode-1) + 1
    ie = hecmat%indexL(inode  )
 
    do k= is, ie
 
      !*row (left)
      do i = 0, ndof - 1
        hecmat%AL(ndof2*k-ii+i) = 0.d0
        if(present(conmat)) conmat%AL(ndof2*k-ii+i) = 0.d0
      end do
 
      !*column (upper)
      in = hecmat%itemL(k)
      iis = hecmat%indexU(in-1) + 1
      iie = hecmat%indexU(in  )
      do ik = iis, iie
        if (hecmat%itemU(ik) .eq. inode) then
          iii = ndof2 - idof
          do i = ndof-1,0,-1
            idx = ndof*in-i
            val = hecmat%AU(ndof2*ik-iii)*rhs
            !$omp atomic
            hecmat%B(idx) = hecmat%B(idx) - val
            hecmat%AU(ndof2*ik-iii)= 0.d0
            if(present(conmat)) then
              val = conmat%AU(ndof2*ik-iii)*rhs
              !$omp atomic
              conmat%B(idx) = conmat%B(idx) - val
              conmat%AU(ndof2*ik-iii)= 0.d0
            endif
            iii = iii - ndof
          end do
          exit
        endif
      enddo
 
    enddo
 
    ii = ndof2-1 - (idof-1)*ndof
    is = hecmat%indexU(inode-1) + 1
    ie = hecmat%indexU(inode  )
 
    do k= is, ie
 
      !*row (right)
      do i = 0,ndof-1
        hecmat%AU(ndof2*k-ii+i) = 0.d0
        if(present(conmat)) conmat%AU(ndof2*k-ii+i) = 0.d0
      end do
 
      !*column (lower)
      in = hecmat%itemU(k)
      iis = hecmat%indexL(in-1) + 1
      iie = hecmat%indexL(in  )
      do ik= iis, iie
        if (hecmat%itemL(ik) .eq. inode) then
          iii  = ndof2 - idof
 
          do i = ndof-1, 0, -1
            idx = ndof*in-i
            val = hecmat%AL(ndof2*ik-iii)*rhs
            !$omp atomic
            hecmat%B(idx) = hecmat%B(idx) - val
            hecmat%AL(ndof2*ik-iii) = 0.d0
            if(present(conmat)) then
              val = conmat%AL(ndof2*ik-iii)*rhs
              !$omp atomic
              conmat%B(idx) = conmat%B(idx) - val
              conmat%AL(ndof2*ik-iii) = 0.d0
            endif
            iii = iii - ndof
          end do
          exit
        endif
      enddo
 
    enddo
    !*End off - diagonal blocks
 
  end subroutine hecmw_mat_ass_bc
 
 
  subroutine hecmw_mat_ass_bc_contactlag(hecMAT,hecLagMAT,inode,idof,RHS)
 
    type(hecmwst_matrix)                 :: hecmat
    type(hecmwst_matrix_lagrange)        :: heclagmat
    integer(kind=kint) :: inode, idof
    integer(kind=kint) :: isu, ieu, isl, iel, i, l, k
    real(kind=kreal)   :: rhs 
 
    isu = heclagmat%indexU_lagrange(inode-1)+1
    ieu = heclagmat%indexU_lagrange(inode)
    do i = isu, ieu
      heclagmat%AU_lagrange((i-1)*3+idof) = 0.0d0
      l = heclagmat%itemU_lagrange(i)
      isl = heclagmat%indexL_lagrange(l-1)+1
      iel = heclagmat%indexL_lagrange(l)
      k = hecmw_array_search_i(heclagmat%itemL_lagrange,isl,iel,inode)
      if(k < isl .or. k > iel) cycle
      hecmat%B(hecmat%NP*hecmat%NDOF+l) = hecmat%B(hecmat%NP*hecmat%NDOF+l) - heclagmat%AL_lagrange((k-1)*3+idof)*rhs
      heclagmat%AL_lagrange((k-1)*3+idof) = 0.0d0
    enddo
 
  end subroutine hecmw_mat_ass_bc_contactlag
 
  subroutine hecmw_mat_ass_contactlag(nnode,ndLocal,id_lagrange,fcoeff,stiffness,hecMAT,hecLagMAT)
 
    type(hecmwst_matrix)                 :: hecmat
    type(hecmwst_matrix_lagrange)        :: heclagmat
    integer(kind=kint) :: nnode, ndlocal(nnode + 1), id_lagrange
    integer(kind=kint) :: i, j, inod, jnod, l
    integer(kind=kint) :: isl, iel, idxl_base, kl, idxl, isu, ieu, idxu_base, ku, idxu
    real(kind=kreal)   :: fcoeff 
    real(kind=kreal)   :: stiffness(:,:) 
    real(kind=kreal)   :: a(3, 3)
 
    i = nnode + 1 + 1
    inod = id_lagrange
    isl = heclagmat%indexL_lagrange(inod-1)+1
    iel = heclagmat%indexL_lagrange(inod)
 
    do j = 1, nnode + 1
      jnod = ndlocal(j)
      isu = heclagmat%indexU_lagrange(jnod-1)+1
      ieu = heclagmat%indexU_lagrange(jnod)
 
      kl = hecmw_array_search_i(heclagmat%itemL_lagrange,isl,iel,jnod)
      if( kl<isl .or. kl>iel ) then
        write(*,*) '###ERROR### : cannot find connectivity (Lagrange1)'
        stop
      endif
      ku = hecmw_array_search_i(heclagmat%itemU_lagrange,isu,ieu,inod)
      if( ku<isu .or. ku>ieu ) then
        write(*,*) '###ERROR### : cannot find connectivity (Lagrange2)'
        stop
      endif
 
      idxl_base = (kl-1)*3
      idxu_base = (ku-1)*3
 
      do l = 1, 3
        idxl = idxl_base + l
        heclagmat%AL_lagrange(idxl) = heclagmat%AL_lagrange(idxl) + stiffness((i-1)*3+1,(j-1)*3+l)
        idxu = idxu_base + l
        heclagmat%AU_lagrange(idxu) = heclagmat%AU_lagrange(idxu) + stiffness((j-1)*3+l,(i-1)*3+1)
      enddo
    enddo
 
 
    if(fcoeff /= 0.0d0)then
 
      do i = 1, nnode + 1
        inod = ndlocal(i)
        do j = 1, nnode + 1
          jnod = ndlocal(j)
          call stf_get_block(stiffness(1:(nnode+1)*3,1:(nnode+1)*3), 3, i, j, a)
          call hecmw_mat_add_node(hecmat, inod, jnod, a)
        enddo
      enddo
 
    endif
 
  end subroutine hecmw_mat_ass_contactlag
 
end module hecmw_matrix_ass