FrontISTR  5.9.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  if (associated(heclagmat%D_lagrange)) then
484  spmat%A(m) = heclagmat%D_lagrange(i)
485  else
486  spmat%A(m) = 0.d0
487  endif
488  m=m+1
489  end if
490  enddo
491  endif
492 
493  if (spmat%type == sparse_matrix_type_csr .and. spmat%IRN(ii+1-spmat%offset)/=m) &
494  stop "ERROR: sparse_matrix_contact_set_a"
495  if (m-1 /= spmat%NZ) stop "ERROR: sparse_matrix_contact_set_a"
496  end subroutine sparse_matrix_contact_set_vals
497 
498  subroutine sparse_matrix_para_contact_set_vals(spMAT, hecMAT, hecLagMAT, conMAT)
499  type(sparse_matrix), intent(inout) :: spmat
500  type(hecmwst_matrix), intent(in) :: hecmat
501  type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
502  type(hecmwst_matrix), intent(in) :: conmat
503  integer(kind=kint) :: ndof, ndof2
504  integer(kind=kint) :: m, i, idof, i0, ii, ls, le, l, j, j0, jdof, jdofs
505  integer(kind=kint) :: offset_l, offset_d, offset_u
506 
507  if(spmat%type /= sparse_matrix_type_coo) then
508  write(*,*) 'ERROR: sparse_matrix_para_contact_set_vals on rank ',hecmw_comm_get_rank()
509  write(*,*) 'spMAT%type must be SPARSE_MATRIX_TYPE_COO, only for mumps'
510  call hecmw_abort(hecmw_comm_get_comm())
511  endif
512  ! write(myrank+20,*)'matrix values'
513  ndof=hecmat%NDOF; ndof2=ndof*ndof
514  m=1
515  do i=1,hecmat%NP
516  ! Internal Nodes First
517  if(i <= hecmat%N) then
518  i0=spmat%offset+ndof*(i-1)
519  do idof=1,ndof
520  ii=i0+idof
521  ! Lower
522  if (.not. sparse_matrix_is_sym(spmat)) then
523  ls=hecmat%indexL(i-1)+1
524  le=hecmat%indexL(i)
525  do l=ls,le
526  j=hecmat%itemL(l)
527  if (j > hecmat%N) stop 'j > hecMAT%N'
528  j0=spmat%offset+ndof*(j-1)
529  offset_l=ndof2*(l-1)+ndof*(idof-1)
530  do jdof=1,ndof
531  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
532  stop "ERROR: sparse_matrix_contact_set_a"
533  if (spmat%JCN(m)/=j0+jdof) &
534  stop "ERROR: sparse_matrix_contact_set_a"
535  spmat%A(m)=hecmat%AL(offset_l+jdof) + conmat%AL(offset_l+jdof)
536  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
537  m=m+1
538  enddo
539  enddo
540  endif
541  ! Diag
542  offset_d=ndof2*(i-1)+ndof*(idof-1)
543  if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
544  do jdof=jdofs,ndof
545  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
546  stop "ERROR: sparse_matrix_contact_set_a"
547  if (spmat%JCN(m)/=i0+jdof) &
548  stop "ERROR: sparse_matrix_contact_set_a"
549  ! print *,myrank,'spMAT',size(spMAT%A),m,size(hecMAT%D),size(conMAT%D),offset_d,jdof
550  spmat%A(m)=hecmat%D(offset_d+jdof) + conmat%D(offset_d+jdof)
551  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
552  m=m+1
553  enddo
554  ! Upper
555  ls=hecmat%indexU(i-1)+1
556  le=hecmat%indexU(i)
557  do l=ls,le
558  j=hecmat%itemU(l)
559  if (j <= hecmat%N) then
560  j0=spmat%offset+ndof*(j-1)
561  else
562  j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
563  if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
564  endif
565  offset_u=ndof2*(l-1)+ndof*(idof-1)
566  do jdof=1,ndof
567  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
568  stop "ERROR: sparse_matrix_contact_set_a"
569  if (spmat%JCN(m)/=j0+jdof) &
570  stop "ERROR: sparse_matrix_contact_set_a"
571  spmat%A(m)=hecmat%AU(offset_u+jdof) + conmat%AU(offset_u+jdof)
572  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
573  m=m+1
574  enddo
575  enddo
576  ! Upper Lagrange
577  if (heclagmat%num_lagrange > 0) then
578  j0=spmat%offset+ndof*hecmat%N
579  ls=heclagmat%indexU_lagrange(i-1)+1
580  le=heclagmat%indexU_lagrange(i)
581  do l=ls,le
582  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
583  stop "ERROR: sparse_matrix_contact_set_a"
584  if (spmat%JCN(m)/=j0+heclagmat%itemU_lagrange(l)) &
585  stop "ERROR: sparse_matrix_contact_set_a"
586  spmat%A(m)=heclagmat%AU_lagrange((l-1)*ndof+idof)
587  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
588  m=m+1
589  enddo
590  endif
591  enddo
592  else
593  ! External Nodes
594  i0 = spmat%conv_ext(ndof*(i-hecmat%N))-ndof
595  do idof=1,ndof
596  ii=i0+idof
597  ! Lower
598  ls=hecmat%indexL(i-1)+1
599  le=hecmat%indexL(i)
600  do l=ls,le
601  j=hecmat%itemL(l)
602  if (j <= hecmat%N) then
603  j0=spmat%offset+ndof*(j-1)
604  else
605  j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
606  endif
607  if(sparse_matrix_is_sym(spmat).and.j0 < i0) cycle
608  offset_l=ndof2*(l-1)+ndof*(idof-1)
609  do jdof=1,ndof
610  ! if (spMAT%type==SPARSE_MATRIX_TYPE_COO) spMAT%IRN(m)=ii
611  ! spMAT%JCN(m)=j0+jdof
612  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
613  stop "ERROR: sparse_matrix_contact_set_a"
614  if (spmat%JCN(m)/=j0+jdof) &
615  stop "ERROR: sparse_matrix_contact_set_a"
616  spmat%A(m) = conmat%AL(offset_l+jdof)
617  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
618  m=m+1
619  enddo
620  enddo
621 
622  ! Diag
623  offset_d=ndof2*(i-1)+ndof*(idof-1)
624  if (sparse_matrix_is_sym(spmat)) then; jdofs=idof; else; jdofs=1; endif
625  do jdof=jdofs,ndof
626  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
627  stop "ERROR: sparse_matrix_contact_set_a"
628  if (spmat%JCN(m)/=i0+jdof) &
629  stop "ERROR: sparse_matrix_contact_set_a"
630  ! if (spMAT%type==SPARSE_MATRIX_TYPE_COO) spMAT%IRN(m)=ii
631  ! spMAT%JCN(m)=i0+jdof
632  spmat%A(m) = conmat%D(offset_d+jdof)
633  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
634  m=m+1
635  enddo
636  !
637  ! Upper
638  ls=hecmat%indexU(i-1)+1
639  le=hecmat%indexU(i)
640  do l=ls,le
641  j=hecmat%itemU(l)
642  if (j <= hecmat%N) then
643  j0=spmat%offset+ndof*(j-1)
644  else
645  j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
646  endif
647  if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
648  offset_u=ndof2*(l-1)+ndof*(idof-1)
649  do jdof=1,ndof
650  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
651  stop "ERROR: sparse_matrix_contact_set_a"
652  if (spmat%JCN(m)/=j0+jdof) stop "ERROR: sparse_matrix_contact_set_a"
653  spmat%A(m) = conmat%AU(offset_u+jdof)
654  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
655  m=m+1
656  enddo
657  enddo
658  !
659  ! Upper Lagrange
660  if (heclagmat%num_lagrange > 0) then
661  j0=spmat%offset+ndof*hecmat%N
662  ls=heclagmat%indexU_lagrange(i-1)+1
663  le=heclagmat%indexU_lagrange(i)
664  if (sparse_matrix_is_sym(spmat) .and. j0 < i0) then
665  ! Do nothing
666  else
667  do l=ls,le
668  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
669  stop "ERROR: sparse_matrix_contact_set_a"
670  if (spmat%JCN(m)/=j0+heclagmat%itemU_lagrange(l)) &
671  stop "ERROR: sparse_matrix_contact_set_a"
672  spmat%A(m)=heclagmat%AU_lagrange((l-1)*ndof+idof)
673  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
674  m=m+1
675  enddo
676  endif
677  endif
678  enddo
679 
680  endif
681  enddo
682  ! Lower ROW Lagrange
683  if (heclagmat%num_lagrange > 0) then
684  i0=spmat%offset+ndof*hecmat%N
685  do i=1,heclagmat%num_lagrange
686  ii=i0+i
687  ls=heclagmat%indexL_lagrange(i-1)+1
688  le=heclagmat%indexL_lagrange(i)
689  do l=ls,le
690  j=heclagmat%itemL_lagrange(l)
691  if (j <= hecmat%N) then
692  j0=spmat%offset+ndof*(j-1)
693  else
694  j0=spmat%conv_ext(ndof*(j-hecmat%N))-ndof
695  endif
696  if (sparse_matrix_is_sym(spmat) .and. j0 < i0) cycle
697  do jdof=1,ndof
698  if (spmat%type==sparse_matrix_type_coo .and. spmat%IRN(m)/=ii) &
699  stop "ERROR: sparse_matrix_contact_set_a"
700  if (spmat%JCN(m)/=j0+jdof) &
701  stop "ERROR: sparse_matrix_contact_set_a"
702  spmat%A(m)=heclagmat%AL_lagrange((l-1)*ndof+jdof)
703  ! write(myrank+20,*)m,spMAT%IRN(m),spMAT%JCN(m),spMAT%A(m)
704  m=m+1
705  enddo
706  enddo
707  enddo
708  endif
709  ! print *,myrank,'spMAT%irn',ii,spMAT%offset, m,spMAT%NZ
710  if (spmat%type == sparse_matrix_type_csr) then
711  if(spmat%IRN(ii+1-spmat%offset)/=m) &
712  stop "ERROR: sparse_matrix_contact_set_a"
713  endif
714  if (m-1 /= spmat%NZ) then
715  write(*,*) 'ERROR: sparse_matrix_para_contact_set_vals on rank ',hecmw_comm_get_rank()
716  write(*,*) 'm-1 = ',m-1,', NZ=',spmat%NZ,',num_lagrange=',heclagmat%num_lagrange
717  call hecmw_abort(hecmw_comm_get_comm())
718  endif
719  end subroutine sparse_matrix_para_contact_set_vals
720 
721  subroutine sparse_matrix_contact_set_rhs(spMAT, hecMAT, hecLagMAT)
722  implicit none
723  type (sparse_matrix), intent(inout) :: spmat
724  type (hecmwst_matrix), intent(in) :: hecmat
725  type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
726  integer :: nndof,npndof,ierr,i
727  allocate(spmat%rhs(spmat%N_loc), stat=ierr)
728  if (ierr /= 0) then
729  write(*,*) " Allocation error, spMAT%rhs"
730  call hecmw_abort(hecmw_comm_get_comm())
731  endif
732  nndof=hecmat%N*hecmat%ndof
733  npndof=hecmat%NP*hecmat%ndof
734  do i=1,nndof
735  spmat%rhs(i) = hecmat%b(i)
736  enddo
737  ! skip external DOF in hecMAT%b
738  if (heclagmat%num_lagrange > 0) then
739  do i=1,heclagmat%num_lagrange
740  spmat%rhs(nndof+i) = hecmat%b(npndof+i)
741  enddo
742  endif
743  end subroutine sparse_matrix_contact_set_rhs
744 
745  subroutine sparse_matrix_para_contact_set_rhs(spMAT, hecMAT, hecLagMAT, conMAT)
746  implicit none
747  type (sparse_matrix), intent(inout) :: spmat
748  type (hecmwst_matrix), intent(in) :: hecmat
749  type (hecmwst_matrix_lagrange), intent(in) :: heclagmat
750  type (hecmwst_matrix), intent(in) :: conmat
751  integer :: nndof,npndof,ierr,ndof,i,i0,j
752  real(kind=kreal), allocatable :: rhs_con(:), rhs_con_sum(:)
753 
754  ! assemble contact components of external nodes only
755  allocate(rhs_con(spmat%N),stat=ierr)
756  rhs_con(:) = 0.0d0
757  ndof = hecmat%ndof
758  do i= hecmat%N+1,hecmat%NP
759  i0=spmat%conv_ext(ndof*(i-hecmat%N))-ndof
760  if((i0 < 0.or.i0 > spmat%N)) then
761  do j=1,ndof
762  if(conmat%b((i-1)*ndof+j) /= 0.0d0) then
763  print *,hecmw_comm_get_rank(),'i0',i,spmat%N,'conMAT%b',conmat%b((i-1)*ndof+j)
764  stop
765  endif
766  enddo
767  else
768  if(i0 > spmat%N - ndof) then
769  print *,hecmw_comm_get_rank(),'ext out',hecmat%N,hecmat%NP,i,i0
770  endif
771  do j=1,ndof
772  if(conmat%b((i-1)*ndof+j) /= 0.0d0) then
773  rhs_con(i0+j) = conmat%b((i-1)*ndof+j)
774  endif
775  enddo
776  endif
777  enddo
778  allocate(rhs_con_sum(spmat%N))
779  call hecmw_allreduce_dp(rhs_con, rhs_con_sum, spmat%N, hecmw_sum, hecmw_comm_get_comm())
780  deallocate(rhs_con)
781 
782  allocate(spmat%rhs(spmat%N_loc), stat=ierr)
783  if (ierr /= 0) then
784  write(*,*) " Allocation error, spMAT%rhs"
785  call hecmw_abort(hecmw_comm_get_comm())
786  endif
787  nndof=hecmat%N*hecmat%ndof
788  npndof=hecmat%NP*hecmat%ndof
789  do i=1,nndof
790  spmat%rhs(i) = rhs_con_sum(spmat%offset+i) + hecmat%b(i) + conmat%b(i)
791  end do
792  deallocate(rhs_con_sum)
793  if (heclagmat%num_lagrange > 0) then
794  do i=1,heclagmat%num_lagrange
795  spmat%rhs(nndof+i) = conmat%b(npndof+i)
796  end do
797  endif
798  end subroutine sparse_matrix_para_contact_set_rhs
799 
800  subroutine sparse_matrix_contact_get_rhs(spMAT, hecMAT, hecLagMAT)
801  implicit none
802  type (sparse_matrix), intent(inout) :: spmat
803  type (hecmwst_matrix), intent(inout) :: hecmat
804  type (hecmwst_matrix_lagrange), intent(inout) :: heclagmat
805  integer :: nndof,npndof,i
806  nndof=hecmat%N*hecmat%ndof
807  npndof=hecmat%NP*hecmat%ndof
808  do i=1,nndof
809  hecmat%x(i) = spmat%rhs(i)
810  enddo
811  ! TODO: call update to get external DOF from neighboring domains???
812  if (heclagmat%num_lagrange > 0) then
813  do i=1,heclagmat%num_lagrange
814  hecmat%x(npndof+i) = spmat%rhs(nndof+i)
815  enddo
816  endif
817  deallocate(spmat%rhs)
818  end subroutine sparse_matrix_contact_get_rhs
819 
820 end module m_sparse_matrix_contact
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::hecmw_sum
integer(kind=kint), parameter hecmw_sum
Definition: hecmw_util_f.F90:23
hecmw_util::hecmw_comm_get_size
integer(kind=kint) function hecmw_comm_get_size()
Definition: hecmw_util_f.F90:601
hecmw_util::hecmw_comm_get_comm
integer(kind=kint) function hecmw_comm_get_comm()
Definition: hecmw_util_f.F90:579
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_util::hecmw_comm_get_rank
integer(kind=kint) function hecmw_comm_get_rank()
Definition: hecmw_util_f.F90:590
hecmw_util::hecmw_abort
subroutine hecmw_abort(comm)
Definition: hecmw_util_f.F90:542
m_hecmw_comm_f
Definition: hecmw_comm_f.F90:6
m_hecmw_comm_f::hecmw_allreduce_i1
subroutine hecmw_allreduce_i1(hecMESH, s, ntag)
Definition: hecmw_comm_f.F90:489
m_hecmw_comm_f::hecmw_allreduce_dp
subroutine hecmw_allreduce_dp(val, VALM, n, hec_op, comm)
Definition: hecmw_comm_f.F90:320
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_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:499
m_sparse_matrix_contact::sparse_matrix_contact_get_rhs
subroutine, public sparse_matrix_contact_get_rhs(spMAT, hecMAT, hecLagMAT)
Definition: sparse_matrix_contact.f90:801
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::sparse_matrix_contact_set_rhs
subroutine, public sparse_matrix_contact_set_rhs(spMAT, hecMAT, hecLagMAT)
Definition: sparse_matrix_contact.f90:722
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_para_contact_set_rhs
subroutine, public sparse_matrix_para_contact_set_rhs(spMAT, hecMAT, hecLagMAT, conMAT)
Definition: sparse_matrix_contact.f90:746
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
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
This module provides DOF based sparse matrix data structure (CSR and COO)
Definition: sparse_matrix.f90:6
m_sparse_matrix::sparse_matrix_init
subroutine, public sparse_matrix_init(spMAT, N_loc, NZ)
Definition: sparse_matrix.f90:80
m_sparse_matrix::sparse_matrix_is_sym
logical function, public sparse_matrix_is_sym(spMAT)
Definition: sparse_matrix.f90:172
m_sparse_matrix::sparse_matrix_type_coo
integer(kind=kint), parameter, public sparse_matrix_type_coo
Definition: sparse_matrix.f90:30
m_sparse_matrix::sparse_matrix_type_csr
integer(kind=kint), parameter, public sparse_matrix_type_csr
Definition: sparse_matrix.f90:29
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:236
hecmw_util::hecmwst_matrix_lagrange
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...
Definition: hecmw_util_f.F90:429
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:449