FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_matrix_contact_lagrange.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 !-------------------------------------------------------------------------------
5 
7  use hecmw_util
8  implicit none
9 
10  private
11  public :: noderelated
18 
20  type noderelated
21  integer(kind=kint) :: num_node = 0, num_lagrange = 0
22  integer(kind=kint), pointer :: id_node(:) => null()
23  integer(kind=kint), pointer :: id_lagrange(:) => null()
24  end type
25 
26 contains
27 
29  subroutine hecmw_construct_noderelated_from_hecmat(hecMAT, NPL, NPU, list_nodeRelated)
30  type(hecmwst_matrix), intent(in) :: hecmat
31  integer(kind=kint), intent(out) :: npl
32  integer(kind=kint), intent(out) :: npu
33  type(noderelated), pointer, intent(inout) :: list_noderelated(:)
34 
35  integer(kind=kint) :: numl, numu, num_noderelated
36  integer(kind=kint) :: i, j
37  integer(kind=kint) :: ierr
38 
39  npl = hecmat%NPL; npu = hecmat%NPU
40 
41  if( associated(list_noderelated) ) return
42  allocate(list_noderelated(hecmat%NP), stat=ierr)
43  if( ierr /= 0) stop " Allocation error, list_nodeRelated "
44 
45  do i = 1, hecmat%NP
46 
47  numl = hecmat%indexL(i) - hecmat%indexL(i-1)
48  numu = hecmat%indexU(i) - hecmat%indexU(i-1)
49 
50  num_noderelated = numl + numu + 1
51 
52  allocate(list_noderelated(i)%id_node(num_noderelated), stat=ierr)
53  if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_node "
54 
55  list_noderelated(i)%num_node = num_noderelated
56 
57  do j = 1, numl
58  list_noderelated(i)%id_node(j) = hecmat%itemL(hecmat%indexL(i-1)+j)
59  enddo
60  list_noderelated(i)%id_node(numl+1) = i
61  do j = 1, numu
62  list_noderelated(i)%id_node(numl+1+j) = hecmat%itemU(hecmat%indexU(i-1)+j)
63  enddo
64 
65  enddo
66 
68 
70  subroutine hecmw_init_noderelated_from_org(np,num_lagrange,is_contact_active,list_nodeRelated_org,list_nodeRelated)
71  integer(kind=kint), intent(in) :: np, num_lagrange
72  logical, intent(in) :: is_contact_active
73  type(noderelated), pointer, intent(in) :: list_noderelated_org(:)
74  type(noderelated), pointer, intent(inout) :: list_noderelated(:)
75 
76  integer(kind=kint) :: num_noderelated_org
77  integer(kind=kint) :: i, ierr
78 
79  allocate(list_noderelated(np+num_lagrange), stat=ierr)
80  if( ierr /= 0) stop " Allocation error, list_nodeRelated "
81 
82  do i = 1, np !hecMAT%NP
83  num_noderelated_org = list_noderelated_org(i)%num_node
84  allocate(list_noderelated(i)%id_node(num_noderelated_org), stat=ierr)
85  if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_node "
86  list_noderelated(i)%num_node = num_noderelated_org
87  list_noderelated(i)%id_node(1:num_noderelated_org) = list_noderelated_org(i)%id_node(1:num_noderelated_org)
88  enddo
89 
90  if( is_contact_active ) then
91  do i = np+1, np+num_lagrange !hecMAT%NP+1, hecMAT%NP+num_lagrange
92  allocate(list_noderelated(i)%id_lagrange(5), stat=ierr)
93  if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_lagrange "
94  list_noderelated(i)%num_lagrange = 0
95  list_noderelated(i)%id_lagrange = 0
96  enddo
97  endif
98 
99  end subroutine hecmw_init_noderelated_from_org
100 
101  subroutine hecmw_ass_noderelated_from_contact_pair(np, nnode, ndLocal, count_lagrange, permission, &
102  & necessary_to_insert_node, list_nodeRelated_org, list_nodeRelated, countNon0LU_node, countNon0LU_lagrange )
103  integer(kind=kint), intent(in) :: np
104  integer(kind=kint), intent(in) :: nnode
105  integer(kind=kint), intent(in) :: ndlocal(:)
106  integer(kind=kint), intent(in) :: count_lagrange
107  logical, intent(inout) :: permission
108  logical, intent(in) :: necessary_to_insert_node
109  type(noderelated), pointer, intent(in) :: list_noderelated_org(:)
110  type(noderelated), pointer, intent(inout) :: list_noderelated(:)
111  integer(kind=kint), intent(inout) :: countnon0lu_node, countnon0lu_lagrange
113 
114  integer(kind=kint) :: k, l, num, num_noderelated_org, ierr
115 
116  do k = 1, nnode+1
117 
118  if( .not. associated(list_noderelated(ndlocal(k))%id_lagrange) )then
119  num = 10
120  ! if( k == 1 ) num = 1
121  allocate(list_noderelated(ndlocal(k))%id_lagrange(num),stat=ierr)
122  if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_lagrange "
123  list_noderelated(ndlocal(k))%num_lagrange = 0
124  list_noderelated(ndlocal(k))%id_lagrange = 0
125  endif
126 
127  if( necessary_to_insert_node ) then
128  num_noderelated_org = list_noderelated_org(ndlocal(k))%num_node
129  if( list_noderelated(ndlocal(k))%num_node == num_noderelated_org )then
130  num = 10
131  if(k==1) num = 4
132  call reallocate_memory(num,list_noderelated(ndlocal(k)))
133  endif
134  endif
135 
136  if( count_lagrange > 0 ) call insert_lagrange(k,count_lagrange,list_noderelated(ndlocal(k)),countnon0lu_lagrange,permission)
137 
138  do l = k, nnode+1
139  if( necessary_to_insert_node ) then
140  if( k /= l) then
141  num_noderelated_org = list_noderelated_org(ndlocal(k))%num_node
142  call insert_node(ndlocal(l),list_noderelated(ndlocal(k)),countnon0lu_node)
143  num_noderelated_org = list_noderelated_org(ndlocal(l))%num_node
144  if( list_noderelated(ndlocal(l))%num_node == num_noderelated_org )then
145  num = 10
146  call reallocate_memory(num,list_noderelated(ndlocal(l)))
147  endif
148  call insert_node(ndlocal(k),list_noderelated(ndlocal(l)),countnon0lu_node)
149  endif
150  endif
151 
152  if(k == 1 .and. count_lagrange > 0) &
153  call insert_lagrange(0,ndlocal(l),list_noderelated(np+count_lagrange),countnon0lu_lagrange,permission)
154 
155  enddo
156 
157  enddo
159 
161  subroutine insert_lagrange(i,id_lagrange,list_node,countNon0_lagrange,permission)
162 
163  type(noderelated), intent(inout) :: list_node
164  integer(kind=kint), intent(in) :: i, id_lagrange
166  integer(kind=kint), intent(inout) :: countnon0_lagrange
168  logical, intent(inout) :: permission
169 
170  integer(kind=kint) :: ierr, num_lagrange, location
171  integer(kind=kint), allocatable :: id_lagrange_save(:)
172 
173  character(len=1) :: answer
174 
175  ierr = 0
176 
177  num_lagrange = count(list_node%id_lagrange /= 0 )
178 
179  ! if( i == 1 .and. num_lagrange /= 0) return
180  if( i == 1 .and. num_lagrange /= 0 .and. .not. permission) then
181  1 write(*,*) '##Error: node is both slave and master node simultaneously !'
182  write(*,*) ' Please check contact surface definition !'
183  write(*,'('' Do you want to continue(y/n)):'',$)')
184  read(*,'(A1)',err=1) answer
185  if(answer == 'Y' .OR. answer == 'y')then
186  permission = .true.
187  else
188  stop
189  endif
190  endif
191 
192  if (num_lagrange == 0)then
193  list_node%num_lagrange = 1
194  list_node%id_lagrange(1) = id_lagrange
195  countnon0_lagrange = countnon0_lagrange + 1
196  else
197  allocate(id_lagrange_save(num_lagrange))
198  id_lagrange_save(1:num_lagrange) = list_node%id_lagrange(1:num_lagrange)
199  location = find_locationinarray(id_lagrange,num_lagrange,list_node%id_lagrange)
200  if(location /= 0)then
201  num_lagrange = num_lagrange + 1
202  if( num_lagrange > size(list_node%id_lagrange)) then
203  deallocate(list_node%id_lagrange)
204  allocate(list_node%id_lagrange(num_lagrange),stat=ierr)
205  if( ierr /= 0 ) stop " Allocation error, list_nodeRelated%id_lagrange "
206  endif
207  list_node%num_lagrange = num_lagrange
208  list_node%id_lagrange(location) = id_lagrange
209  if(location /= 1) list_node%id_lagrange(1:location-1) = id_lagrange_save(1:location-1)
210  if(location /= num_lagrange) list_node%id_lagrange(location+1:num_lagrange) = id_lagrange_save(location:num_lagrange-1)
211  countnon0_lagrange = countnon0_lagrange + 1
212  endif
213  deallocate(id_lagrange_save)
214  endif
215 
216  end subroutine insert_lagrange
217 
219  subroutine insert_node(id_node,list_node,countNon0_node)
220 
221  type(noderelated) :: list_node
222  integer(kind=kint) :: id_node
224  integer(kind=kint) :: countnon0_node
225  integer(kind=kint) :: ierr, num_node, location
226  integer(kind=kint),allocatable :: id_node_save(:)
227 
228  ierr = 0
229 
230  num_node = list_node%num_node
231  allocate(id_node_save(num_node))
232  id_node_save(1:num_node) = list_node%id_node(1:num_node)
233  location = find_locationinarray(id_node,num_node,list_node%id_node)
234  if(location /= 0)then
235  num_node = num_node + 1
236  if( num_node > size(list_node%id_node)) then
237  deallocate(list_node%id_node)
238  allocate(list_node%id_node(num_node),stat=ierr)
239  if( ierr /= 0) stop " Allocation error, list_nodeRelated%id_node "
240  endif
241  list_node%num_node = num_node
242  list_node%id_node(location) = id_node
243  if(location /= 1) list_node%id_node(1:location-1) = id_node_save(1:location-1)
244  if(location /= num_node) list_node%id_node(location+1:num_node) = id_node_save(location:num_node-1)
245  countnon0_node = countnon0_node + 1
246  endif
247  deallocate(id_node_save)
248 
249  end subroutine insert_node
250 
252  integer function find_locationinarray(item,n,a)
253  integer(kind=kint), intent(in) :: item, n
254  integer(kind=kint), pointer, intent(inout) :: a(:)
255 
256  integer(kind=kint) :: l, r, m
257 
258  find_locationinarray = 0
259 
260  l = 1 ; r = n ; m = (l+r)/2
261  if( item == a(l) .or. item == a(r) )then
262  return
263  elseif( item < a(l) )then
264  find_locationinarray = 1
265  return
266  elseif( item > a(r) )then
267  find_locationinarray = n + 1
268  return
269  endif
270 
271  do while ( l <= r)
272  if( item > a(m) ) then
273  l = m + 1
274  m = (l + r)/2
275  elseif( item < a(m) ) then
276  r = m - 1
277  m = (l + r)/2
278  elseif( item == a(m) )then
279  return
280  endif
281  enddo
282 
283  find_locationinarray = m + 1
284 
285  end function find_locationinarray
286 
288  subroutine reallocate_memory(num,list_node)
289  integer(kind=kint), intent(in) :: num
290  type(noderelated),intent(inout) :: list_node
291 
292  integer(kind=kint) :: num_node_org
294  integer(kind=kint) :: id_save(1000)
295  integer(kind=kint) :: ierr
296 
297  num_node_org = size(list_node%id_node)
298  id_save(1:num_node_org) = list_node%id_node(1:num_node_org)
299  deallocate(list_node%id_node)
300  allocate(list_node%id_node(num_node_org+num),stat=ierr)
301  if( ierr /= 0) stop " reAllocation error, list_nodeRelated%id_node "
302  list_node%id_node = 0
303  list_node%id_node(1:num_node_org) = id_save(1:num_node_org)
304 
305  end subroutine reallocate_memory
306 
308  subroutine hecmw_construct_hecmat_from_noderelated(n, np, ndof, nnz, num_lagrange, list_nodeRelated, hecMAT)
309  integer(kind=kint), intent(in) :: n
310  integer(kind=kint), intent(in) :: np
311  integer(kind=kint), intent(in) :: ndof
312  integer(kind=kint), intent(in) :: nnz
313  integer(kind=kint), intent(in) :: num_lagrange
314  type(noderelated), pointer, intent(in) :: list_noderelated(:)
315  type(hecmwst_matrix), intent(inout) :: hecmat
316 
317  integer(kind=kint) :: countnon0l_node, countnon0u_node
318  integer(kind=kint) :: numi_node
319  integer(kind=kint) :: i, j, ierr
320  integer(kind=kint) :: ndof2
321 
322  ndof2 = ndof*ndof
323 
324  hecmat%N = n
325  hecmat%NP = np
326  hecmat%NDOF = ndof
327  hecmat%NPL = nnz
328  hecmat%NPU = nnz
329 
330  if(associated(hecmat%indexL)) deallocate(hecmat%indexL)
331  if(associated(hecmat%indexU)) deallocate(hecmat%indexU)
332  if(associated(hecmat%itemL)) deallocate(hecmat%itemL)
333  if(associated(hecmat%itemU)) deallocate(hecmat%itemU)
334  if(associated(hecmat%AL)) deallocate(hecmat%AL)
335  if(associated(hecmat%AU)) deallocate(hecmat%AU)
336  if(associated(hecmat%B)) deallocate(hecmat%B)
337  if(associated(hecmat%X)) deallocate(hecmat%X)
338  if(associated(hecmat%D)) deallocate(hecmat%D)
339 
340  allocate(hecmat%indexL(0:np), stat=ierr)
341  if ( ierr /= 0) stop " Allocation error, hecMAT%indexL "
342  hecmat%indexL = 0
343 
344  allocate(hecmat%indexU(0:np), stat=ierr)
345  if ( ierr /= 0) stop " Allocation error, hecMAT%indexU "
346  hecmat%indexU = 0
347 
348  allocate(hecmat%itemL(nnz), stat=ierr)
349  if ( ierr /= 0) stop " Allocation error, hecMAT%itemL "
350  hecmat%itemL = 0
351 
352  allocate(hecmat%itemU(nnz), stat=ierr)
353  if ( ierr /= 0) stop " Allocation error, hecMAT%itemU "
354  hecmat%itemU = 0
355 
356  countnon0l_node = 0
357  countnon0u_node = 0
358  do i = 1, np
359  numi_node = count(list_noderelated(i)%id_node /= 0)
360 
361  do j = 1, numi_node
362  if( list_noderelated(i)%id_node(j) < i ) then
363  countnon0l_node = countnon0l_node + 1
364  hecmat%itemL(countnon0l_node) = list_noderelated(i)%id_node(j)
365  elseif( list_noderelated(i)%id_node(j) > i ) then
366  countnon0u_node = countnon0u_node + 1
367  hecmat%itemU(countnon0u_node) = list_noderelated(i)%id_node(j)
368  endif
369  enddo
370  hecmat%indexL(i) = countnon0l_node
371  hecmat%indexU(i) = countnon0u_node
372  end do
373 
374  allocate(hecmat%AL(nnz*ndof2), stat=ierr)
375  if ( ierr /= 0 ) stop " Allocation error, hecMAT%AL "
376  hecmat%AL = 0.0d0
377 
378  allocate(hecmat%AU(nnz*ndof2), stat=ierr)
379  if ( ierr /= 0 ) stop " Allocation error, hecMAT%AU "
380  hecmat%AU = 0.0d0
381 
382  allocate(hecmat%D(np*ndof2+num_lagrange))
383  hecmat%D = 0.0d0
384 
385  allocate(hecmat%B(np*ndof+num_lagrange))
386  hecmat%B = 0.0d0
387 
388  allocate(hecmat%X(np*ndof+num_lagrange))
389  hecmat%X = 0.0d0
390 
391  end subroutine
392 
395  & np,ndof,num_lagrange,numNon0_lagrange,is_contact_active,list_nodeRelated,hecLagMAT)
396  integer(kind=kint), intent(in) :: np
397  integer(kind=kint), intent(in) :: ndof
398  integer(kind=kint), intent(in) :: num_lagrange
399  integer(kind=kint), intent(in) :: numnon0_lagrange
400  logical, intent(in) :: is_contact_active
401  type(noderelated), pointer, intent(in) :: list_noderelated(:)
402  type(hecmwst_matrix_lagrange), intent(inout) :: heclagmat
403 
404  integer(kind=kint) :: countnon0u_lagrange, countnon0l_lagrange
405  integer(kind=kint) :: numi_lagrange
406  integer(kind=kint) :: i, j, ierr
407 
408  heclagmat%num_lagrange = num_lagrange
409  heclagmat%numL_lagrange = numnon0_lagrange
410  heclagmat%numU_lagrange = numnon0_lagrange
411 
412  if(associated(heclagmat%indexL_lagrange)) deallocate(heclagmat%indexL_lagrange)
413  if(associated(heclagmat%indexU_lagrange)) deallocate(heclagmat%indexU_lagrange)
414  if(associated(heclagmat%itemL_lagrange)) deallocate(heclagmat%itemL_lagrange)
415  if(associated(heclagmat%itemU_lagrange)) deallocate(heclagmat%itemU_lagrange)
416  if(associated(heclagmat%AL_lagrange)) deallocate(heclagmat%AL_lagrange)
417  if(associated(heclagmat%AU_lagrange)) deallocate(heclagmat%AU_lagrange)
418  if(associated(heclagmat%Lagrange)) deallocate(heclagmat%Lagrange)
419 
420  if( is_contact_active ) then
421  ! init indexU_lagrange
422  allocate(heclagmat%indexU_lagrange(0:np), stat=ierr)
423  if ( ierr /= 0) stop " Allocation error, hecLagMAT%indexU_lagrange "
424  heclagmat%indexU_lagrange = 0
425 
426  ! init itemU_lagrange
427  allocate(heclagmat%itemU_lagrange(numnon0_lagrange), stat=ierr)
428  if ( ierr /= 0) stop " Allocation error, hecLagMAT%itemU_lagrange "
429  heclagmat%itemU_lagrange = 0
430 
431  ! setup upper lagrange CRS matrix
432  countnon0u_lagrange = 0
433  do i = 1, np
434  numi_lagrange = list_noderelated(i)%num_lagrange
435  do j = 1, numi_lagrange
436  countnon0u_lagrange = countnon0u_lagrange + 1
437  heclagmat%itemU_lagrange(countnon0u_lagrange) = list_noderelated(i)%id_lagrange(j)
438  enddo
439  heclagmat%indexU_lagrange(i) = countnon0u_lagrange
440  end do
441 
442  ! init indexL_lagrange
443  allocate(heclagmat%indexL_lagrange(0:num_lagrange), stat=ierr)
444  if ( ierr /= 0) stop " Allocation error, hecLagMAT%indexL_lagrange "
445  heclagmat%indexL_lagrange = 0
446 
447  ! init itemL_lagrange
448  allocate(heclagmat%itemL_lagrange(numnon0_lagrange), stat=ierr)
449  if ( ierr /= 0) stop " Allocation error, hecLagMAT%itemL_lagrange "
450  heclagmat%itemL_lagrange = 0
451 
452  ! setup lower lagrange CRS matrix
453  countnon0l_lagrange = 0
454  do i = 1, num_lagrange
455  numi_lagrange = list_noderelated(np+i)%num_lagrange
456  do j = 1, numi_lagrange
457  countnon0l_lagrange = countnon0l_lagrange + 1
458  heclagmat%itemL_lagrange(countnon0l_lagrange) = list_noderelated(np+i)%id_lagrange(j)
459  enddo
460  heclagmat%indexL_lagrange(i) = countnon0l_lagrange
461  enddo
462 
463  ! init AU_lagrange
464  allocate(heclagmat%AU_lagrange(ndof*numnon0_lagrange), stat=ierr)
465  if ( ierr /= 0 ) stop " Allocation error, hecLagMAT%AU_lagrange "
466  heclagmat%AU_lagrange = 0.0d0
467 
468  ! init AL_lagrange
469  allocate(heclagmat%AL_lagrange(ndof*numnon0_lagrange), stat=ierr)
470  if ( ierr /= 0 ) stop " Allocation error, hecLagMAT%AL_lagrange "
471  heclagmat%AL_lagrange = 0.0d0
472 
473  ! init Lagrange
474  allocate(heclagmat%Lagrange(num_lagrange))
475  heclagmat%Lagrange = 0.0d0
476  endif
477 
478  end subroutine
479 
481  subroutine hecmw_finalize_noderelated(list_nodeRelated)
482  type(noderelated), pointer, intent(inout) :: list_noderelated(:)
483 
484  integer(kind=kint) :: i
485 
486  do i = 1, size(list_noderelated)
487  list_noderelated(i)%num_node = 0
488  list_noderelated(i)%num_lagrange = 0
489  if(associated(list_noderelated(i)%id_node)) deallocate(list_noderelated(i)%id_node)
490  if(associated(list_noderelated(i)%id_lagrange)) deallocate(list_noderelated(i)%id_lagrange)
491  end do
492 
493  deallocate(list_noderelated)
494  end subroutine
495 
hecmw_matrix_contact_lagrange
Definition: hecmw_matrix_contact_lagrange.f90:6
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_matrix_contact_lagrange::hecmw_finalize_noderelated
subroutine, public hecmw_finalize_noderelated(list_nodeRelated)
This subroutine finalize list_nodeRelated.
Definition: hecmw_matrix_contact_lagrange.f90:482
hecmw_matrix_contact_lagrange::hecmw_init_noderelated_from_org
subroutine, public hecmw_init_noderelated_from_org(np, num_lagrange, is_contact_active, list_nodeRelated_org, list_nodeRelated)
Get original list of related nodes.
Definition: hecmw_matrix_contact_lagrange.f90:71
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_matrix_contact_lagrange::hecmw_ass_noderelated_from_contact_pair
subroutine, public hecmw_ass_noderelated_from_contact_pair(np, nnode, ndLocal, count_lagrange, permission, necessary_to_insert_node, list_nodeRelated_org, list_nodeRelated, countNon0LU_node, countNon0LU_lagrange)
Definition: hecmw_matrix_contact_lagrange.f90:103
hecmw_matrix_contact_lagrange::hecmw_construct_hecmat_from_noderelated
subroutine, public hecmw_construct_hecmat_from_noderelated(n, np, ndof, nnz, num_lagrange, list_nodeRelated, hecMAT)
This subroutine constructs hecMAT structure from list_nodeRelated.
Definition: hecmw_matrix_contact_lagrange.f90:309
hecmw_matrix_contact_lagrange::hecmw_construct_noderelated_from_hecmat
subroutine, public hecmw_construct_noderelated_from_hecmat(hecMAT, NPL, NPU, list_nodeRelated)
This subroutine saves original matrix structure constructed originally by hecMW_matrix.
Definition: hecmw_matrix_contact_lagrange.f90:30
hecmw_matrix_contact_lagrange::hecmw_construct_heclagmat_from_noderelated
subroutine, public hecmw_construct_heclagmat_from_noderelated(np, ndof, num_lagrange, numNon0_lagrange, is_contact_active, list_nodeRelated, hecLagMAT)
Construct hecLagMAT structure.
Definition: hecmw_matrix_contact_lagrange.f90:396
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444