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