FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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
10  use m_sparse_matrix
12  implicit none
13 
14  private
19  !
22 
23 contains
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
66  end subroutine sparse_matrix_contact_init_prof
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 
821 end module m_sparse_matrix_contact
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_get_rhs(spMAT, hecMAT, fstrMAT)
subroutine, public sparse_matrix_contact_set_rhs(spMAT, hecMAT, fstrMAT)
subroutine, public sparse_matrix_contact_set_vals(spMAT, hecMAT, fstrMAT)
subroutine, public sparse_matrix_para_contact_set_vals(spMAT, hecMAT, fstrMAT, conMAT)
subroutine, public sparse_matrix_para_contact_set_rhs(spMAT, hecMAT, fstrMAT, conMAT)
subroutine, public sparse_matrix_contact_init_prof(spMAT, hecMAT, fstrMAT, hecMESH)
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)
subroutine, public sparse_matrix_init(spMAT, N_loc, NZ)
logical function, public sparse_matrix_is_sym(spMAT)
integer(kind=kint), parameter, public sparse_matrix_type_coo
integer(kind=kint), parameter, public sparse_matrix_type_csr
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...