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