FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
sparse_matrix_contact.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
8 use m_fstr
12 implicit none
13
14 private
19 !
22
23contains
24
25 !!$#define NEED_DIAG
26
27 subroutine sparse_matrix_contact_init_prof(spMAT, hecMAT, fstrMAT, hecMESH)
28 type (sparse_matrix), intent(inout) :: spmat
29 type (hecmwst_matrix), intent(in) :: hecmat
30 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
31 type (hecmwst_local_mesh), intent(in) :: hecmesh
32 integer(kind=kint) :: ndof, ndof2, n_loc, nl, nu, nz, nn,nlag
33 ndof=hecmat%NDOF; ndof2=ndof*ndof
34 ! ---- For Parallel Contact with Multi-Partition Domains
35 if(paracontactflag) then
36 nn = hecmat%NP
37 else
38 nn = hecmat%N
39 endif
40 n_loc=hecmat%N*ndof+fstrmat%num_lagrange
41 if (sparse_matrix_is_sym(spmat)) then
42 nu=hecmat%indexU(nn)
43 nz=nn*(ndof2+ndof)/2+nu*ndof2 &
44 +fstrmat%numU_lagrange*ndof
45 else
46 nl=hecmat%indexL(nn)
47 nu=hecmat%indexU(nn)
48 nz=(nn+nu+nl)*ndof2 &
49 +(fstrmat%numL_lagrange+fstrmat%numU_lagrange)*ndof
50 endif
51 ! diagonal elements must be allocated for CRS format even if they are all zero
52 if (spmat%type==sparse_matrix_type_csr) nz=nz+fstrmat%numL_lagrange
53 call sparse_matrix_init(spmat, n_loc, nz)
54 call sparse_matrix_hec_set_conv_ext(spmat, hecmesh, ndof)
55 if(fstrmat%num_lagrange > 0) &
56 print '(I3,A,4I10,A,2I10)',myrank,' sparse_matrix_init ',hecmat%N,hecmat%NP,n_loc,nz,' LAG',fstrmat%num_lagrange,spmat%offset
57 nlag = fstrmat%num_lagrange
58 call hecmw_allreduce_i1(hecmesh, nlag, hecmw_sum)
59 if(myrank == 0) print *,'total number of contact nodes:',nlag
60 spmat%timelog = hecmat%Iarray(22)
61 if(paracontactflag) then
62 call sparse_matrix_para_contact_set_prof(spmat, hecmat, fstrmat)
63 else
64 call sparse_matrix_contact_set_prof(spmat, hecmat, fstrmat)
65 endif
67
68 subroutine sparse_matrix_contact_set_prof(spMAT, hecMAT, fstrMAT)
69 type(sparse_matrix), intent(inout) :: spmat
70 type(hecmwst_matrix), intent(in) :: hecmat
71 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
72 integer(kind=kint) :: ndof, ndof2
73 integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
74 !integer(kind=kint) :: offset_l, offset_d, offset_u
75 ! CONVERT TO CSR or COO STYLE
76 ndof=hecmat%NDOF; ndof2=ndof*ndof
77 m=1
78 ii=0
79 do i=1,hecmat%N
80 do idof=1,ndof
81 i0=spmat%offset+ndof*(i-1)
82 ii=i0+idof
83 if (spmat%type==sparse_matrix_type_csr) spmat%IRN(ii-spmat%offset)=m
84 ! Lower
85 if (.not. sparse_matrix_is_sym(spmat)) then
86 ls=hecmat%indexL(i-1)+1
87 le=hecmat%indexL(i)
88 do l=ls,le
89 j=hecmat%itemL(l)
90 !if (j <= hecMAT%N) then
91 j0=spmat%offset+ndof*(j-1)
92 !else
93 ! j0=spMAT%conv_ext(ndof*(j-hecMAT%N))-ndof
94 !endif
95 !offset_l=ndof2*(l-1)+ndof*(idof-1)
96 do jdof=1,ndof
97 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
98 spmat%JCN(m)=j0+jdof
99 !spMAT%A(m)=hecMAT%AL(offset_l+jdof)
100 m=m+1
101 enddo
102 enddo
103 endif
104 ! Diag
105 !offset_d=ndof2*(i-1)+ndof*(idof-1)
106 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
107 do jdof=jdofs,ndof
108 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
109 spmat%JCN(m)=i0+jdof
110 !spMAT%A(m)=hecMAT%D(offset_d+jdof)
111 m=m+1
112 enddo
113 ! Upper
114 ls=hecmat%indexU(i-1)+1
115 le=hecmat%indexU(i)
116 do l=ls,le
117 j=hecmat%itemU(l)
118 if (j <= hecmat%N) then
119 j0=spmat%offset+ndof*(j-1)
120 else
121 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
122 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
123 endif
124 !offset_u=ndof2*(l-1)+ndof*(idof-1)
125 do jdof=1,ndof
126 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
127 spmat%JCN(m)=j0+jdof
128 !spMAT%A(m)=hecMAT%AU(offset_u+jdof)
129 m=m+1
130 enddo
131 enddo
132 ! Upper Lagrange
133 if (fstrmat%num_lagrange > 0) then
134 j0=spmat%offset+ndof*hecmat%N
135 ls=fstrmat%indexU_lagrange(i-1)+1
136 le=fstrmat%indexU_lagrange(i)
137 do l=ls,le
138 j=fstrmat%itemU_lagrange(l)
139 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
140 spmat%JCN(m)=j0+j
141 m=m+1
142 enddo
143 endif
144 enddo
145 enddo
146 ! Lower Lagrange
147 if (fstrmat%num_lagrange > 0) then
148 i0=spmat%offset+ndof*hecmat%N
149 do i=1,fstrmat%num_lagrange
150 ii=i0+i
151 if (spmat%type==sparse_matrix_type_csr) spmat%IRN(ii-spmat%offset)=m
152 if (.not. sparse_matrix_is_sym(spmat)) then
153 ls=fstrmat%indexL_lagrange(i-1)+1
154 le=fstrmat%indexL_lagrange(i)
155 do l=ls,le
156 j=fstrmat%itemL_lagrange(l)
157 j0=spmat%offset+ndof*(j-1)
158 do jdof=1,ndof
159 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
160 spmat%JCN(m)=j0+jdof
161 m=m+1
162 enddo
163 enddo
164 endif
165 ! Diag( Only for CSR )
166 if (spmat%type==sparse_matrix_type_csr) then
167 spmat%JCN(m)=ii
168 m=m+1
169 end if
170 enddo
171 endif
172 if (spmat%type == sparse_matrix_type_csr) spmat%IRN(ii+1-spmat%offset)=m
173 if (m-1 < spmat%NZ) spmat%NZ=m-1
174 if (m-1 /= spmat%NZ) then
175 write(*,*) 'ERROR: sparse_matrix_contact_set_prof on rank ',myrank
176 write(*,*) 'm-1 = ',m-1,', NZ=',spmat%NZ,',num_lagrange=',fstrmat%num_lagrange
177 call hecmw_abort(hecmw_comm_get_comm())
178 endif
179 end subroutine sparse_matrix_contact_set_prof
180
181 subroutine sparse_matrix_para_contact_set_prof(spMAT, hecMAT, fstrMAT)
182 type(sparse_matrix), intent(inout) :: spmat
183 type(hecmwst_matrix), intent(in) :: hecmat
184 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
185 integer(kind=kint) :: ndof, ndof2
186 integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
187 !integer(kind=kint) :: offset_l, offset_d, offset_u
188 ! CONVERT TO CSR or COO STYLE
189 if(spmat%type /= sparse_matrix_type_coo) then
190 write(*,*) 'ERROR: sparse_matrix_para_contact_set_prof on rank ',myrank
191 write(*,*) 'spMAT%type must be SPARSE_MATRIX_TYPE_COO, only for mumps'
192 call hecmw_abort(hecmw_comm_get_comm())
193 endif
194 ! write(myrank+20,*)'matrix profile'
195 ndof=hecmat%NDOF; ndof2=ndof*ndof
196 m=1
197 do i=1,hecmat%NP
198 ! Internal Nodes First
199 if(i <= hecmat%N) then
200 do idof=1,ndof
201 i0=spmat%offset+ndof*(i-1)
202 ii=i0+idof
203 ! Exact Lower Triangle (No External Nodes)
204 if(.not. sparse_matrix_is_sym(spmat)) then
205 ls=hecmat%indexL(i-1)+1
206 le=hecmat%indexL(i)
207 do l=ls,le
208 j=hecmat%itemL(l)
209 j0=spmat%offset+ndof*(j-1)
210 !offset_l=ndof2*(l-1)+ndof*(idof-1)
211 do jdof=1,ndof
212 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
213 spmat%JCN(m)=j0+jdof
214 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
215 !spMAT%A(m)=hecMAT%AL(offset_l+jdof)
216 m=m+1
217 enddo
218 enddo
219 endif
220 ! Diag Part
221 !offset_d=ndof2*(i-1)+ndof*(idof-1)
222 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
223 do jdof=jdofs,ndof
224 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
225 spmat%JCN(m)=i0+jdof
226 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
227 !spMAT%A(m)=hecMAT%D(offset_d+jdof)
228 m=m+1
229 enddo
230 ! Upper Triangle (Possible External Nodes)
231 ls=hecmat%indexU(i-1)+1
232 le=hecmat%indexU(i)
233 do l=ls,le
234 j=hecmat%itemU(l)
235 if (j <= hecmat%N) then
236 j0=spmat%offset+ndof*(j-1)
237 else
238 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
239 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
240 endif
241 !offset_u=ndof2*(l-1)+ndof*(idof-1)
242 do jdof=1,ndof
243 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
244 spmat%JCN(m)=j0+jdof
245 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
246 !spMAT%A(m)=hecMAT%AU(offset_u+jdof)
247 m=m+1
248 enddo
249 enddo
250 ! Upper COL Lagrange
251 if (fstrmat%num_lagrange > 0) then
252 j0=spmat%offset+ndof*hecmat%N
253 ls=fstrmat%indexU_lagrange(i-1)+1
254 le=fstrmat%indexU_lagrange(i)
255 do l=ls,le
256 j=fstrmat%itemU_lagrange(l)
257 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
258 spmat%JCN(m)=j0+j
259 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
260 m=m+1
261 enddo
262 endif
263 enddo
264 else
265 ! External Nodes
266 i0 = spmat%conv_ext(ndof*(i-hecmat%N))-ndof
267 do idof=1,ndof
268 ii=i0+idof
269 ! Lower
270 ls=hecmat%indexL(i-1)+1
271 le=hecmat%indexL(i)
272 do l=ls,le
273 j=hecmat%itemL(l)
274 if (j <= hecmat%N) then
275 j0=spmat%offset+ndof*(j-1)
276 else
277 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
278 endif
279 if(sparse_matrix_is_sym(spmat).and.j0 < i0) cycle
280 do jdof=1,ndof
281 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
282 spmat%JCN(m)=j0+jdof
283 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
284 !spMAT%A(m)=hecMAT%AL(offset_l+jdof)
285 m=m+1
286 enddo
287 enddo
288
289 ! Diag
290 !offset_d=ndof2*(i-1)+ndof*(idof-1)
291 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
292 do jdof=jdofs,ndof
293 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
294 spmat%JCN(m)=i0+jdof
295 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
296 !spMAT%A(m)=hecMAT%D(offset_d+jdof)
297 m=m+1
298 enddo
299 !
300 ! Upper
301 ls=hecmat%indexU(i-1)+1
302 le=hecmat%indexU(i)
303 do l=ls,le
304 j=hecmat%itemU(l)
305 if (j <= hecmat%N) then
306 j0=spmat%offset+ndof*(j-1)
307 else
308 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
309 endif
310 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
311 !offset_u=ndof2*(l-1)+ndof*(idof-1)
312 do jdof=1,ndof
313 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
314 spmat%JCN(m)=j0+jdof
315 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
316 !spMAT%A(m)=hecMAT%AU(offset_u+jdof)
317 m=m+1
318 enddo
319 enddo
320 !
321 ! Upper COL Lagrange
322 if (fstrmat%num_lagrange > 0) then
323 j0=spmat%offset+ndof*hecmat%N
324 ls=fstrmat%indexU_lagrange(i-1)+1
325 le=fstrmat%indexU_lagrange(i)
326 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) then
327
328 else
329 do l=ls,le
330 j=fstrmat%itemU_lagrange(l)
331 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
332 spmat%JCN(m)=j0+j
333 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
334 m=m+1
335 enddo
336 endif
337 endif
338 enddo
339
340 endif
341 enddo
342 ! Lower ROW Lagrange
343 if (fstrmat%num_lagrange > 0) then
344 i0=spmat%offset+ndof*hecmat%N
345 do i=1,fstrmat%num_lagrange
346 ii=i0+i
347 ls=fstrmat%indexL_lagrange(i-1)+1
348 le=fstrmat%indexL_lagrange(i)
349 do l=ls,le
350 j=fstrmat%itemL_lagrange(l)
351 if (j <= hecmat%N) then
352 j0=spmat%offset+ndof*(j-1)
353 else
354 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
355 endif
356 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
357 ! j0=spMAT%OFFSET+ndof*(j-1)
358 do jdof=1,ndof
359 if (spmat%type==sparse_matrix_type_coo) spmat%IRN(m)=ii
360 spmat%JCN(m)=j0+jdof
361 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m)
362 m=m+1
363 enddo
364 enddo
365 enddo
366 endif
367
368 ! if (sparse_matrix_is_sym(spMAT) .and. m-1 < spMAT%NZ) spMAT%NZ=m-1
369 if (m-1 /= spmat%NZ) then
370 write(*,*) 'ERROR: sparse_matrix_para_contact_set_prof on rank ',myrank
371 write(*,*) 'm-1 = ',m-1,', NZ=',spmat%NZ,',num_lagrange=',fstrmat%num_lagrange
372 call hecmw_abort(hecmw_comm_get_comm())
373 endif
374 end subroutine sparse_matrix_para_contact_set_prof
375
376 subroutine sparse_matrix_contact_set_vals(spMAT, hecMAT, fstrMAT)
377 type(sparse_matrix), intent(inout) :: spmat
378 type(hecmwst_matrix), intent(in) :: hecmat
379 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
380 integer(kind=kint) :: ndof, ndof2
381 integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
382 integer(kind=kint) :: offset_l, offset_d, offset_u
383 ndof=hecmat%NDOF; ndof2=ndof*ndof
384 m=1
385 ii=0
386 do i=1,hecmat%N
387 do idof=1,ndof
388 i0=spmat%offset+ndof*(i-1)
389 ii=i0+idof
390 if (spmat%type==sparse_matrix_type_csr .and. spmat%IRN(ii-spmat%offset)/=m) &
391 stop "ERROR: sparse_matrix_contact_set_a"
392 ! Lower
393 if (.not. sparse_matrix_is_sym(spmat)) then
394 ls=hecmat%indexL(i-1)+1
395 le=hecmat%indexL(i)
396 do l=ls,le
397 j=hecmat%itemL(l)
398 !if (j <= hecMAT%N) then
399 j0=spmat%offset+ndof*(j-1)
400 !else
401 ! j0=spMAT%conv_ext(ndof*(j-hecMAT%N))-ndof
402 !endif
403 offset_l=ndof2*(l-1)+ndof*(idof-1)
404 do jdof=1,ndof
405 if ( spmat%type==sparse_matrix_type_coo ) then
406 if( spmat%IRN(m)/=ii ) stop "ERROR: sparse_matrix_contact_set_a"
407 end if
408 if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_contact_set_a"
409 spmat%A(m)=hecmat%AL(offset_l+jdof)
410 m=m+1
411 enddo
412 enddo
413 endif
414 ! Diag
415 offset_d=ndof2*(i-1)+ndof*(idof-1)
416 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
417 do jdof=jdofs,ndof
418 if ( spmat%type==sparse_matrix_type_coo ) then
419 if( spmat%IRN(m)/=ii ) stop "ERROR: sparse_matrix_contact_set_a"
420 end if
421 if (spmat%JCN(m)/=i0+jdof) stop "ERROR: sparse_matrix_contact_set_a"
422 spmat%A(m)=hecmat%D(offset_d+jdof)
423 m=m+1
424 enddo
425 ! Upper
426 ls=hecmat%indexU(i-1)+1
427 le=hecmat%indexU(i)
428 do l=ls,le
429 j=hecmat%itemU(l)
430 if (j <= hecmat%N) then
431 j0=spmat%offset+ndof*(j-1)
432 else
433 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
434 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
435 endif
436 offset_u=ndof2*(l-1)+ndof*(idof-1)
437 do jdof=1,ndof
438 if ( spmat%type==sparse_matrix_type_coo ) then
439 if( spmat%IRN(m)/=ii ) stop "ERROR: sparse_matrix_contact_set_a"
440 end if
441 if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_contact_set_a"
442 spmat%A(m)=hecmat%AU(offset_u+jdof)
443 m=m+1
444 enddo
445 enddo
446 ! Upper Lagrange
447 if (fstrmat%num_lagrange > 0) then
448 j0=spmat%offset+ndof*hecmat%N
449 ls=fstrmat%indexU_lagrange(i-1)+1
450 le=fstrmat%indexU_lagrange(i)
451 do l=ls,le
452 if ( spmat%type==sparse_matrix_type_coo ) then
453 if( spmat%IRN(m)/=ii ) stop "ERROR: sparse_matrix_contact_set_a"
454 end if
455 if (spmat%JCN(m)/=j0+fstrmat%itemU_lagrange(l)) &
456 stop "ERROR: sparse_matrix_contact_set_a"
457 spmat%A(m)=fstrmat%AU_lagrange((l-1)*ndof+idof)
458 m=m+1
459 enddo
460 endif
461 enddo
462 enddo
463 ! Lower Lagrange
464 if (fstrmat%num_lagrange > 0) then
465 i0=spmat%offset+ndof*hecmat%N
466 do i=1,fstrmat%num_lagrange
467 ii=i0+i
468 if (spmat%type==sparse_matrix_type_csr) spmat%IRN(ii-spmat%offset)=m
469 if (.not. sparse_matrix_is_sym(spmat)) then
470 ls=fstrmat%indexL_lagrange(i-1)+1
471 le=fstrmat%indexL_lagrange(i)
472 do l=ls,le
473 j=fstrmat%itemL_lagrange(l)
474 j0=spmat%offset+ndof*(j-1)
475 do jdof=1,ndof
476 if ( spmat%type==sparse_matrix_type_coo ) then
477 if( spmat%IRN(m)/=ii ) stop "ERROR: sparse_matrix_contact_set_a"
478 end if
479 if (spmat%JCN(m)/=j0+jdof) &
480 stop "ERROR: sparse_matrix_contact_set_a"
481 spmat%A(m)=fstrmat%AL_lagrange((l-1)*ndof+jdof)
482 m=m+1
483 enddo
484 enddo
485 endif
486 ! Diag( Only for CSR )
487 if (spmat%type==sparse_matrix_type_csr) then
488 spmat%A(m)=0.d0
489 m=m+1
490 end if
491 enddo
492 endif
493
494 if (spmat%type == sparse_matrix_type_csr .and. spmat%IRN(ii+1-spmat%offset)/=m) &
495 stop "ERROR: sparse_matrix_contact_set_a"
496 if (m-1 /= spmat%NZ) stop "ERROR: sparse_matrix_contact_set_a"
497 end subroutine sparse_matrix_contact_set_vals
498
499 subroutine sparse_matrix_para_contact_set_vals(spMAT, hecMAT, fstrMAT, conMAT)
500 type(sparse_matrix), intent(inout) :: spmat
501 type(hecmwst_matrix), intent(in) :: hecmat
502 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
503 type(hecmwst_matrix), intent(in) :: conmat
504 integer(kind=kint) :: ndof, ndof2
505 integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
506 integer(kind=kint) :: offset_l, offset_d, offset_u
507
508 if(spmat%type /= sparse_matrix_type_coo) then
509 write(*,*) 'ERROR: sparse_matrix_para_contact_set_vals on rank ',myrank
510 write(*,*) 'spMAT%type must be SPARSE_MATRIX_TYPE_COO, only for mumps'
511 call hecmw_abort(hecmw_comm_get_comm())
512 endif
513 ! write(myrank+20,*)'matrix values'
514 ndof=hecmat%NDOF; ndof2=ndof*ndof
515 m=1
516 do i=1,hecmat%NP
517 ! Internal Nodes First
518 if(i <= hecmat%N) then
519 i0=spmat%offset+ndof*(i-1)
520 do idof=1,ndof
521 ii=i0+idof
522 ! Lower
523 if (.not. sparse_matrix_is_sym(spmat)) then
524 ls=hecmat%indexL(i-1)+1
525 le=hecmat%indexL(i)
526 do l=ls,le
527 j=hecmat%itemL(l)
528 if (j > hecmat%N) stop 'j > hecMAT%N'
529 j0=spmat%offset+ndof*(j-1)
530 offset_l=ndof2*(l-1)+ndof*(idof-1)
531 do jdof=1,ndof
532 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
533 stop "ERROR: sparse_matrix_contact_set_a"
534 if (spmat%JCN(m)/=j0+jdof) &
535 stop "ERROR: sparse_matrix_contact_set_a"
536 spmat%A(m)=hecmat%AL(offset_l+jdof) + conmat%AL(offset_l+jdof)
537 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
538 m=m+1
539 enddo
540 enddo
541 endif
542 ! Diag
543 offset_d=ndof2*(i-1)+ndof*(idof-1)
544 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
545 do jdof=jdofs,ndof
546 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
547 stop "ERROR: sparse_matrix_contact_set_a"
548 if (spmat%JCN(m)/=i0+jdof) &
549 stop "ERROR: sparse_matrix_contact_set_a"
550 ! print *,myrank,'spMAT',size(spMAT%A),m,size(hecMAT%D),size(conMAT%D),offset_d,jdof
551 spmat%A(m)=hecmat%D(offset_d+jdof) + conmat%D(offset_d+jdof)
552 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
553 m=m+1
554 enddo
555 ! Upper
556 ls=hecmat%indexU(i-1)+1
557 le=hecmat%indexU(i)
558 do l=ls,le
559 j=hecmat%itemU(l)
560 if (j <= hecmat%N) then
561 j0=spmat%offset+ndof*(j-1)
562 else
563 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
564 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
565 endif
566 offset_u=ndof2*(l-1)+ndof*(idof-1)
567 do jdof=1,ndof
568 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
569 stop "ERROR: sparse_matrix_contact_set_a"
570 if (spmat%JCN(m)/=j0+jdof) &
571 stop "ERROR: sparse_matrix_contact_set_a"
572 spmat%A(m)=hecmat%AU(offset_u+jdof) + conmat%AU(offset_u+jdof)
573 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
574 m=m+1
575 enddo
576 enddo
577 ! Upper Lagrange
578 if (fstrmat%num_lagrange > 0) then
579 j0=spmat%offset+ndof*hecmat%N
580 ls=fstrmat%indexU_lagrange(i-1)+1
581 le=fstrmat%indexU_lagrange(i)
582 do l=ls,le
583 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
584 stop "ERROR: sparse_matrix_contact_set_a"
585 if (spmat%JCN(m)/=j0+fstrmat%itemU_lagrange(l)) &
586 stop "ERROR: sparse_matrix_contact_set_a"
587 spmat%A(m)=fstrmat%AU_lagrange((l-1)*ndof+idof)
588 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
589 m=m+1
590 enddo
591 endif
592 enddo
593 else
594 ! External Nodes
595 i0 = spmat%conv_ext(ndof*(i-hecmat%N))-ndof
596 do idof=1,ndof
597 ii=i0+idof
598 ! Lower
599 ls=hecmat%indexL(i-1)+1
600 le=hecmat%indexL(i)
601 do l=ls,le
602 j=hecmat%itemL(l)
603 if (j <= hecmat%N) then
604 j0=spmat%offset+ndof*(j-1)
605 else
606 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
607 endif
608 if(sparse_matrix_is_sym(spmat).and.j0 < i0) cycle
609 offset_l=ndof2*(l-1)+ndof*(idof-1)
610 do jdof=1,ndof
611 ! if (spMAT%type==SPARSE_MATRIX_TYPE_COO) spMAT%IRN(m)=ii
612 ! spMAT%JCN(m)=j0+jdof
613 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
614 stop "ERROR: sparse_matrix_contact_set_a"
615 if (spmat%JCN(m)/=j0+jdof) &
616 stop "ERROR: sparse_matrix_contact_set_a"
617 spmat%A(m) = conmat%AL(offset_l+jdof)
618 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
619 m=m+1
620 enddo
621 enddo
622
623 ! Diag
624 offset_d=ndof2*(i-1)+ndof*(idof-1)
625 if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
626 do jdof=jdofs,ndof
627 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
628 stop "ERROR: sparse_matrix_contact_set_a"
629 if (spmat%JCN(m)/=i0+jdof) &
630 stop "ERROR: sparse_matrix_contact_set_a"
631 ! if (spMAT%type==SPARSE_MATRIX_TYPE_COO) spMAT%IRN(m)=ii
632 ! spMAT%JCN(m)=i0+jdof
633 spmat%A(m) = conmat%D(offset_d+jdof)
634 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
635 m=m+1
636 enddo
637 !
638 ! Upper
639 ls=hecmat%indexU(i-1)+1
640 le=hecmat%indexU(i)
641 do l=ls,le
642 j=hecmat%itemU(l)
643 if (j <= hecmat%N) then
644 j0=spmat%offset+ndof*(j-1)
645 else
646 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
647 endif
648 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
649 offset_u=ndof2*(l-1)+ndof*(idof-1)
650 do jdof=1,ndof
651 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
652 stop "ERROR: sparse_matrix_contact_set_a"
653 if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_contact_set_a"
654 spmat%A(m) = conmat%AU(offset_u+jdof)
655 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
656 m=m+1
657 enddo
658 enddo
659 !
660 ! Upper Lagrange
661 if (fstrmat%num_lagrange > 0) then
662 j0=spmat%offset+ndof*hecmat%N
663 ls=fstrmat%indexU_lagrange(i-1)+1
664 le=fstrmat%indexU_lagrange(i)
665 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) then
666 ! Do nothing
667 else
668 do l=ls,le
669 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
670 stop "ERROR: sparse_matrix_contact_set_a"
671 if (spmat%JCN(m)/=j0+fstrmat%itemU_lagrange(l)) &
672 stop "ERROR: sparse_matrix_contact_set_a"
673 spmat%A(m)=fstrmat%AU_lagrange((l-1)*ndof+idof)
674 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
675 m=m+1
676 enddo
677 endif
678 endif
679 enddo
680
681 endif
682 enddo
683 ! Lower ROW Lagrange
684 if (fstrmat%num_lagrange > 0) then
685 i0=spmat%offset+ndof*hecmat%N
686 do i=1,fstrmat%num_lagrange
687 ii=i0+i
688 ls=fstrmat%indexL_lagrange(i-1)+1
689 le=fstrmat%indexL_lagrange(i)
690 do l=ls,le
691 j=fstrmat%itemL_lagrange(l)
692 if (j <= hecmat%N) then
693 j0=spmat%offset+ndof*(j-1)
694 else
695 j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
696 endif
697 if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
698 do jdof=1,ndof
699 if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
700 stop "ERROR: sparse_matrix_contact_set_a"
701 if (spmat%JCN(m)/=j0+jdof) &
702 stop "ERROR: sparse_matrix_contact_set_a"
703 spmat%A(m)=fstrmat%AL_lagrange((l-1)*ndof+jdof)
704 ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
705 m=m+1
706 enddo
707 enddo
708 enddo
709 endif
710 ! print *,myrank,'spMAT%irn',ii,spMAT%offset, m,spMAT%NZ
711 if (spmat%type == sparse_matrix_type_csr) then
712 if(spmat%IRN(ii+1-spmat%offset)/=m) &
713 stop "ERROR: sparse_matrix_contact_set_a"
714 endif
715 if (m-1 /= spmat%NZ) then
716 write(*,*) 'ERROR: sparse_matrix_para_contact_set_vals on rank ',myrank
717 write(*,*) 'm-1 = ',m-1,', NZ=',spmat%NZ,',num_lagrange=',fstrmat%num_lagrange
718 call hecmw_abort(hecmw_comm_get_comm())
719 endif
721
722 subroutine sparse_matrix_contact_set_rhs(spMAT, hecMAT, fstrMAT)
723 implicit none
724 type (sparse_matrix), intent(inout) :: spmat
725 type (hecmwst_matrix), intent(in) :: hecmat
726 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
727 integer :: nndof,npndof,ierr,i
728 allocate(spmat%rhs(spmat%N_loc), stat=ierr)
729 if (ierr /= 0) then
730 write(*,*) " Allocation error, spMAT%rhs"
731 call hecmw_abort(hecmw_comm_get_comm())
732 endif
733 nndof=hecmat%N*hecmat%ndof
734 npndof=hecmat%NP*hecmat%ndof
735 do i=1,nndof
736 spmat%rhs(i) = hecmat%b(i)
737 enddo
738 ! skip external DOF in hecMAT%b
739 if (fstrmat%num_lagrange > 0) then
740 do i=1,fstrmat%num_lagrange
741 spmat%rhs(nndof+i) = hecmat%b(npndof+i)
742 enddo
743 endif
744 end subroutine sparse_matrix_contact_set_rhs
745
746 subroutine sparse_matrix_para_contact_set_rhs(spMAT, hecMAT, fstrMAT, conMAT)
747 implicit none
748 type (sparse_matrix), intent(inout) :: spmat
749 type (hecmwst_matrix), intent(in) :: hecmat
750 type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
751 type (hecmwst_matrix), intent(in) :: conmat
752 integer :: nndof,npndof,ierr,ndof,i,i0,j
753 real(kind=kreal), allocatable :: rhs_con(:), rhs_con_sum(:)
754
755 ! assemble contact components of external nodes only
756 allocate(rhs_con(spmat%N),stat=ierr)
757 rhs_con(:) = 0.0d0
758 ndof = hecmat%ndof
759 do i= hecmat%N+1,hecmat%NP
760 i0=spmat%conv_ext(ndof*(i-hecmat%N))-ndof
761 if((i0 < 0.or.i0 > spmat%N)) then
762 do j=1,ndof
763 if(conmat%b((i-1)*ndof+j) /= 0.0d0) then
764 print *,myrank,'i0',i,spmat%N,'conMAT%b',conmat%b((i-1)*ndof+j)
765 stop
766 endif
767 enddo
768 else
769 if(i0 > spmat%N - ndof) then
770 print *,myrank,'ext out',hecmat%N,hecmat%NP,i,i0
771 endif
772 do j=1,ndof
773 if(conmat%b((i-1)*ndof+j) /= 0.0d0) then
774 rhs_con(i0+j) = conmat%b((i-1)*ndof+j)
775 endif
776 enddo
777 endif
778 enddo
779 allocate(rhs_con_sum(spmat%N))
780 call hecmw_allreduce_dp(rhs_con, rhs_con_sum, spmat%N, hecmw_sum, hecmw_comm_get_comm())
781 deallocate(rhs_con)
782
783 allocate(spmat%rhs(spmat%N_loc), stat=ierr)
784 if (ierr /= 0) then
785 write(*,*) " Allocation error, spMAT%rhs"
786 call hecmw_abort(hecmw_comm_get_comm())
787 endif
788 nndof=hecmat%N*hecmat%ndof
789 npndof=hecmat%NP*hecmat%ndof
790 do i=1,nndof
791 spmat%rhs(i) = rhs_con_sum(spmat%offset+i) + hecmat%b(i) + conmat%b(i)
792 end do
793 deallocate(rhs_con_sum)
794 if (fstrmat%num_lagrange > 0) then
795 do i=1,fstrmat%num_lagrange
796 spmat%rhs(nndof+i) = conmat%b(npndof+i)
797 end do
798 endif
800
801 subroutine sparse_matrix_contact_get_rhs(spMAT, hecMAT, fstrMAT)
802 implicit none
803 type (sparse_matrix), intent(inout) :: spmat
804 type (hecmwst_matrix), intent(inout) :: hecmat
805 type (fstrst_matrix_contact_lagrange), intent(inout) :: fstrmat
806 integer :: nndof,npndof,i
807 nndof=hecmat%N*hecmat%ndof
808 npndof=hecmat%NP*hecmat%ndof
809 do i=1,nndof
810 hecmat%x(i) = spmat%rhs(i)
811 enddo
812 ! TODO: call update to get external DOF from neighboring domains???
813 if (fstrmat%num_lagrange > 0) then
814 do i=1,fstrmat%num_lagrange
815 hecmat%x(npndof+i) = spmat%rhs(nndof+i)
816 enddo
817 endif
818 deallocate(spmat%rhs)
819 end subroutine sparse_matrix_contact_get_rhs
820
This module provides functions of reconstructing.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
logical paracontactflag
PARALLEL CONTACT FLAG.
Definition: m_fstr.f90:84
This module provides conversion routines between HEC and FISTR data structures and DOF based sparse m...
subroutine, public sparse_matrix_contact_init_prof(spmat, hecmat, fstrmat, hecmesh)
subroutine, public sparse_matrix_para_contact_set_rhs(spmat, hecmat, fstrmat, conmat)
subroutine, public sparse_matrix_contact_set_vals(spmat, hecmat, fstrmat)
subroutine, public sparse_matrix_contact_get_rhs(spmat, hecmat, fstrmat)
subroutine, public sparse_matrix_contact_set_rhs(spmat, hecmat, fstrmat)
subroutine, public sparse_matrix_para_contact_set_vals(spmat, hecmat, fstrmat, conmat)
This module provides conversion routines between HEC data structure and DOF based sparse matrix struc...
subroutine, public sparse_matrix_hec_set_conv_ext(spmat, hecmesh, ndof)
This module provides DOF based sparse matrix data structure (CSR and COO)
integer(kind=kint), parameter, public sparse_matrix_type_coo
integer(kind=kint), parameter, public sparse_matrix_type_csr
subroutine, public sparse_matrix_init(spmat, n_loc, nz)
logical function, public sparse_matrix_is_sym(spmat)
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...