1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6module hecmw_matrix_con
7  use hecmw_util
8  implicit none
9
10  private
11
12  public :: hecmw_mat_con
13
14  integer(kind=kint) :: NU, NL
15  integer(kind=kint), pointer :: INL(:), INU(:)
16  integer(kind=kint), pointer :: IAL(:,:), IAU(:,:)
17
18contains
19  !C***
20  !C*** MAT_CON for solver
21  !C***
22  !C
23  subroutine hecmw_mat_con ( hecMESH, hecMAT )
24
25    use hecmw_util
26    use hecmw_matrix_contact
27
28    implicit none
29    type (hecmwST_matrix)     :: hecMAT
30    type (hecmwST_local_mesh) :: hecMESH
31
32    call hecmw_mat_con0  (hecMESH, hecMAT)
33    call hecmw_mat_con1  (         hecMAT)
34    call hecmw_cmat_init (hecMAT%cmat)
35
36  end subroutine hecmw_mat_con
37  !C
38  !C***
39  !C*** MAT_CON0 for solver
40  !C***
41  !C
42  subroutine hecmw_mat_con0 ( hecMESH, hecMAT )
43
44    use hecmw_util
45    use hecmw_etype
46
47    implicit none
48    integer(kind=kint) ierr,itype,is,iE,ic_type,nn,icel,j,k,inod
49    type (hecmwST_matrix)     :: hecMAT
50    type (hecmwST_local_mesh) :: hecMESH
51    integer(kind=kint) nid(20)
52
53    integer(kind=kint), dimension(2048) :: NCOL1, NCOL2
54    !C
55    !C +-------+
56    !C | INIT. |
57    !C +-------+
58    !C===
59    hecMAT%NP= hecMESH%n_node
60    hecMAT%N = hecMESH%nn_internal
61
62    NU= 10
63    NL= 10
64
65    allocate (INL(hecMAT%NP), IAL(hecMAT%NP,NL))
66    allocate (INU(hecMAT%NP), IAU(hecMAT%NP,NU))
67
68    INL= 0
69    IAL= 0
70    INU= 0
71    IAU= 0
72    !C===
73    !C
74    !C +----------------------------------------+
75    !C | CONNECTIVITY according to ELEMENT TYPE |
76    !C +----------------------------------------+
77    !C===
78    do
79      ierr = 0
80      do itype= 1, hecMESH%n_elem_type
81        is= hecMESH%elem_type_index(itype-1) + 1
82        iE= hecMESH%elem_type_index(itype  )
83        ic_type= hecMESH%elem_type_item(itype)
84        if( hecmw_is_etype_patch(ic_type) ) cycle
85        !C Set number of nodes
86        nn = hecmw_get_max_node(ic_type)
87        !C element loop
88        do icel= is, iE
89          is= hecMESH%elem_node_index(icel-1)
90          do j=1,nn
91            nid(j)= hecMESH%elem_node_item (is+j)
92          enddo
93          do j=1,nn
94            do k=1,nn
95              if( k .ne. j ) then
96                call hecmw_FIND_NODE( hecMAT,nid(j),nid(k), ierr  )
97                if( ierr.ne.0 ) then
98                  call hecmw_mat_con0_clear (ierr)
99                  exit
100                endif
101              endif
102            enddo
103            if( ierr.ne.0 ) exit
104          enddo
105          if( ierr.ne.0 ) exit
106        enddo
107        if( ierr.ne.0 ) exit
108      enddo
109      if( ierr.eq.0 ) exit
110    enddo
111    !C===
112    !C
113    !C +---------+
114    !C | SORTING |
115    !C +---------+
116    !C===
117    do inod= 1, hecMAT%NP
118      NN= INL(inod)
119      do k= 1, NN
120        NCOL1(k)= IAL(inod,k)
121      enddo
122      call hecmw_mSORT (NCOL1, NCOL2, NN)
123      do k= NN, 1, -1
124        IAL(inod,NN-k+1)= NCOL1(NCOL2(k))
125      enddo
126      NN= INU(inod)
127      do k= 1, NN
128        NCOL1(k)= IAU(inod,k)
129      enddo
130      call hecmw_mSORT (NCOL1, NCOL2, NN)
131      do k= NN, 1, -1
132        IAU(inod,NN-k+1)= NCOL1(NCOL2(k))
133      enddo
134    enddo
135    !C===
136  contains
137    !C
138    !C*** MAT_CON0_CLEAR
139    !C
140    subroutine hecmw_mat_con0_clear (IERR)
141
142      implicit none
143      integer(kind=kint) IERR
144
145      deallocate (INL, IAL, INU, IAU)
146
147      if (IERR.eq.1) NL= NL + 5
148      if (IERR.eq.2) NU= NU + 5
149      allocate (INL(hecMAT%NP),IAL(hecMAT%NP,NL))
150      allocate (INU(hecMAT%NP),IAU(hecMAT%NP,NU))
151
152      INL= 0
153      IAL= 0
154      INU= 0
155      IAU= 0
156
157    end subroutine hecmw_mat_con0_clear
158  end subroutine hecmw_mat_con0
159  !C
160  !C***
161  !C*** FIND_TS_NODE
162  !C***
163  !C
164  subroutine hecmw_FIND_NODE ( hecMAT, ip1,ip2, IERR )
165
166    use hecmw_util
167
168    implicit none
169    integer(kind=kint) ip1,ip2,IERR
170    integer(kind=kint) kk,icou
171    type (hecmwST_matrix)     :: hecMAT
172
173    if (ip1.gt.ip2) then
174      do kk= 1, INL(ip1)
175        if (ip2.eq.IAL(ip1,kk)) return
176      enddo
177      icou= INL(ip1) + 1
178      if (icou.gt.NL) then
179        IERR= 1
180        return
181      endif
182      IAL(ip1,icou)= ip2
183      INL(ip1     )= icou
184      return
185    endif
186
187    if (ip2.gt.ip1) then
188      do kk= 1, INU(ip1)
189        if (ip2.eq.IAU(ip1,kk)) return
190      enddo
191      icou= INU(ip1) + 1
192      if (icou.gt.NU) then
193        IERR= 2
194        return
195      endif
196      IAU(ip1,icou)= ip2
197      INU(ip1     )= icou
198      return
199    endif
200
201  end subroutine hecmw_FIND_NODE
202  !C
203  !C***
204  !C*** fstr_mSORT
205  !C***
206  !C
207  subroutine hecmw_mSORT (STEM,INUM,NN)
208    use hecmw_util
209
210    implicit none
211    integer(kind=kint) NN
212    integer(kind=kint) STEM(NN), INUM(NN)
213    integer(kind=kint) ii,jj,ITEM
214    do ii = 1,NN
215      INUM(ii)= ii
216    enddo
217    do ii= 1,NN-1
218      !CDIR NOVECTOR
219      do jj= 1,NN-ii
220        if (STEM(INUM(jj)) .lt. STEM(INUM(jj+1))) then
221          ITEM      = INUM(jj+1)
222          INUM(jj+1)= INUM(jj)
223          INUM(jj)  = ITEM
224        endif
225      enddo
226    enddo
227    return
228  end subroutine hecmw_mSORT
229  !C
230  !C***
231  !C*** MAT_CON1 for solver
232  !C***
233  !C
234  subroutine hecmw_mat_con1 (hecMAT)
235
236    use hecmw_util
237
238    implicit none
239    integer(kind=kint) i,k,kk
240    type (hecmwST_matrix    ) :: hecMAT
241
242    allocate (hecMAT%indexL(0:hecMAT%NP), hecMAT%indexU(0:hecMAT%NP))
243
244    hecMAT%indexL = 0
245    hecMAT%indexU = 0
246    do i = 1, hecMAT%NP
247      hecMAT%indexL(i) = hecMAT%indexL(i-1) + INL(i)
248      hecMAT%indexU(i) = hecMAT%indexU(i-1) + INU(i)
249    enddo
250
251    hecMAT%NPL = hecMAT%indexL(hecMAT%NP)
252    hecMAT%NPU = hecMAT%indexU(hecMAT%NP)
253
254    allocate (hecMAT%itemL(hecMAT%NPL), hecMAT%itemU(hecMAT%NPU))
255
256    do i = 1, hecMAT%NP
257      do k = 1, INL(i)
258        kk  = k + hecMAT%indexL(i-1)
259        hecMAT%itemL(kk) =     IAL(i,k)
260      enddo
261      do k= 1, INU(i)
262        kk  = k + hecMAT%indexU(i-1)
263        hecMAT%itemU(kk) =     IAU(i,k)
264      enddo
265    enddo
266
267    deallocate (INL, INU, IAL, IAU)
268
269  end subroutine hecmw_mat_con1
270end module hecmw_matrix_con
271