hecmw_mat_con.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_con
  use hecmw_util
  use hecmw_varray_int
  use hecmw_etype
  implicit none
 
  private
 
  public :: hecmw_mat_con
 
contains
  !C***
  !C*** MAT_CON for solver
  !C***
  !C
  subroutine hecmw_mat_con ( hecMESH, hecMAT )
    type (hecmwst_matrix)     :: hecmat
    type (hecmwst_local_mesh) :: hecmesh
 
    type (hecmwst_varray_int), allocatable :: clu(:), cll(:)
    integer(kind=kint)  :: itype, is, ie, iss, ise, ic_type, icel
    integer(kind=kint)  :: i, j, k, kk
    integer(kind=kint)  :: nn, nid(2048), nid1, nid2
 
    !C
    !C +-------+
    !C | INIT. |
    !C +-------+
    !C===
    hecmat%NP= hecmesh%n_node
    hecmat%N = hecmesh%nn_internal
 
    !C===
    !C
    !C +----------------------------------------+
    !C | CONNECTIVITY according to ELEMENT TYPE |
    !C +----------------------------------------+
    !C===
    call hecmw_varray_int_initialize_all( clu, hecmat%NP, 27 )
    call hecmw_varray_int_initialize_all( cll, hecmat%NP, 27 )
    do itype= 1, hecmesh%n_elem_type
      is= hecmesh%elem_type_index(itype-1) + 1
      ie= hecmesh%elem_type_index(itype  )
      ic_type= hecmesh%elem_type_item(itype)
      if( hecmw_is_etype_patch(ic_type) ) cycle
      !C element loop
      do icel= is, ie
        iss = hecmesh%elem_node_index(icel-1)
        ise = hecmesh%elem_node_index(icel  )
        nn = ise-iss
        do j=1,nn
          nid(j) = hecmesh%elem_node_item (iss+j)
        enddo
        do j=1,nn
          nid1 = nid(j)
          do k=1,nn
            nid2 = nid(k)
            if( nid1 < nid2 ) then
              call hecmw_varray_int_insert_if_not_exists( clu(nid1), nid2 )
            else if( nid1 > nid2 ) then
              call hecmw_varray_int_insert_if_not_exists( cll(nid1), nid2 )
            end if
          enddo
        enddo
      enddo
    enddo
 
    !C===
    !C
    !C +---------------------------+
    !C | ALLOCATE hecMAT structure |
    !C +---------------------------+
    !C===
    allocate (hecmat%indexL(0:hecmat%NP), hecmat%indexU(0:hecmat%NP))
 
    hecmat%indexL = 0
    hecmat%indexU = 0
    do i = 1, hecmat%NP
      hecmat%indexL(i) = hecmat%indexL(i-1) + hecmw_varray_int_get_nitem(cll(i))
      hecmat%indexU(i) = hecmat%indexU(i-1) + hecmw_varray_int_get_nitem(clu(i))
    enddo
 
    hecmat%NPL = hecmat%indexL(hecmat%NP)
    hecmat%NPU = hecmat%indexU(hecmat%NP)
 
    allocate (hecmat%itemL(hecmat%NPL), hecmat%itemU(hecmat%NPU))
 
    do i=1,hecmat%NP
      do k=1,hecmw_varray_int_get_nitem(cll(i))
        kk = k + hecmat%indexL(i-1)
        hecmat%itemL(kk) = hecmw_varray_int_get_item(cll(i), k)
      enddo
      do k=1,hecmw_varray_int_get_nitem(clu(i))
        kk = k + hecmat%indexU(i-1)
        hecmat%itemU(kk) = hecmw_varray_int_get_item(clu(i), k)
      enddo
    enddo
 
    call hecmw_varray_int_finalize_all( clu )
    call hecmw_varray_int_finalize_all( cll )
 
  end subroutine
 
end module hecmw_matrix_con