FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_mat_ass.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
11  implicit none
12 
13  private
14 
15  public :: hecmw_mat_ass_elem
16  public :: hecmw_mat_add_node
17  public :: hecmw_array_search_i
18  public :: hecmw_mat_ass_equation
20  public :: hecmw_mat_add_dof
21  public :: hecmw_mat_ass_bc
22  public :: hecmw_mat_ass_contact
23  public :: stf_get_block
24 
25 contains
26 
27  !C
28  !C***
29  !C*** MAT_ASS_ELEM
30  !C***
31  !C
32  subroutine hecmw_mat_ass_elem(hecMAT, nn, nodLOCAL, stiffness)
33  type (hecmwst_matrix) :: hecmat
34  integer(kind=kint) :: nn
35  integer(kind=kint) :: nodlocal(:)
36  real(kind=kreal) :: stiffness(:, :)
37  !** Local variables
38  integer(kind=kint) :: ndof, inod_e, jnod_e, inod, jnod
39  real(kind=kreal) :: a(6,6)
40 
41  ndof = hecmat%NDOF
42 
43  do inod_e = 1, nn
44  inod = nodlocal(inod_e)
45  do jnod_e = 1, nn
46  jnod = nodlocal(jnod_e)
47  !***** Add components
48  call stf_get_block(stiffness, ndof, inod_e, jnod_e, a)
49  call hecmw_mat_add_node(hecmat, inod, jnod, a)
50  enddo
51  enddo
52 
53  end subroutine hecmw_mat_ass_elem
54 
55  subroutine stf_get_block(stiffness, ndof, inod, jnod, a)
56  real(kind=kreal) :: stiffness(:, :), a(:, :)
57  integer(kind=kint) :: ndof, inod, jnod
58  !** Local variables
59  integer(kind=kint) :: row_offset, col_offset, i, j
60 
61  row_offset = ndof*(inod-1)
62  do i = 1, ndof
63 
64  col_offset = ndof*(jnod-1)
65  do j = 1, ndof
66 
67  a(i, j) = stiffness(i + row_offset, j + col_offset)
68  enddo
69  enddo
70  end subroutine stf_get_block
71 
72 
73  subroutine hecmw_mat_add_node(hecMAT, inod, jnod, a)
74  type (hecmwst_matrix) :: hecmat
75  integer(kind=kint) :: inod, jnod
76  real(kind=kreal) :: a(:, :)
77  !** Local variables
78  integer(kind=kint) :: ndof, is, ie, k, idx_base, idx, idof, jdof
79 
80  ndof = hecmat%NDOF
81 
82  if (inod < jnod) then
83  is = hecmat%indexU(inod-1)+1
84  ie = hecmat%indexU(inod)
85  k = hecmw_array_search_i(hecmat%itemU, is, ie, jnod)
86 
87  if (k < is .or. ie < k) then
88  write(*,*) '###ERROR### : cannot find connectivity (1)'
89  write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
91  endif
92 
93  idx_base = ndof**2 * (k-1)
94  do idof = 1, ndof
95  do jdof = 1, ndof
96  idx = idx_base + jdof
97  !$omp atomic
98  hecmat%AU(idx) = hecmat%AU(idx) + a(idof, jdof)
99  enddo
100  idx_base = idx_base + ndof
101  enddo
102 
103  else if (inod > jnod) then
104  is = hecmat%indexL(inod-1)+1
105  ie = hecmat%indexL(inod)
106  k = hecmw_array_search_i(hecmat%itemL, is, ie, jnod)
107 
108  if (k < is .or. ie < k) then
109  write(*,*) '###ERROR### : cannot find connectivity (2)'
110  write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
112  endif
113 
114  idx_base = ndof**2 * (k-1)
115  do idof = 1, ndof
116  do jdof = 1, ndof
117  idx = idx_base + jdof
118  !$omp atomic
119  hecmat%AL(idx) = hecmat%AL(idx) + a(idof, jdof)
120  enddo
121  idx_base = idx_base + ndof
122  enddo
123 
124  else
125  idx_base = ndof**2 * (inod - 1)
126  do idof = 1, ndof
127  do jdof = 1, ndof
128  idx = idx_base + jdof
129  !$omp atomic
130  hecmat%D(idx) = hecmat%D(idx) + a(idof, jdof)
131  enddo
132  idx_base = idx_base + ndof
133  enddo
134  endif
135  end subroutine hecmw_mat_add_node
136 
137 
138  function hecmw_array_search_i(array, is, iE, ival)
139  integer(kind=kint) :: hecmw_array_search_i
140  integer(kind=kint) :: array(:)
141  integer(kind=kint) :: is, ie, ival
142  !** Local variables
143  integer(kind=kint) :: left, right, center, cval
144 
145  left = is
146  right = ie
147  do
148  if (left > right) then
149  center = -1
150  exit
151  endif
152 
153  center = (left + right) / 2
154  cval = array(center)
155 
156  if (ival < cval) then
157  right = center - 1
158  else if (cval < ival) then
159  left = center + 1
160  else
161  exit
162  endif
163  enddo
164 
165  hecmw_array_search_i = center
166 
167  end function hecmw_array_search_i
168 
169 
170  !C
171  !C***
172  !C*** MAT_ASS_EQUATION
173  !C***
174  !C
175  subroutine hecmw_mat_ass_equation ( hecMESH, hecMAT )
176  type (hecmwst_matrix), target :: hecmat
177  type (hecmwst_local_mesh) :: hecmesh
178  !** Local variables
179  real(kind=kreal), pointer :: penalty
180  real(kind=kreal) :: alpha, a1_2inv, ai, aj, factor
181  integer(kind=kint) :: impc, is, ie, i, j, inod, idof, jnod, jdof
182  logical :: is_internal_i, is_internal_j
183 
184  if( hecmw_mat_get_penalized(hecmat) == 1 ) return
185 
186  ! write(*,*) "INFO: imposing MPC by penalty"
187 
188  penalty => hecmat%Rarray(11)
189 
190  if (penalty < 0.0) stop "ERROR: negative penalty"
191  if (penalty < 1.0) write(*,*) "WARNING: penalty ", penalty, " smaller than 1"
192 
193  alpha= hecmw_mat_diag_max(hecmat, hecmesh) * penalty
194  call hecmw_mat_set_penalty_alpha(hecmat, alpha)
195 
196  outer: do impc = 1, hecmesh%mpc%n_mpc
197  is = hecmesh%mpc%mpc_index(impc-1) + 1
198  ie = hecmesh%mpc%mpc_index(impc)
199 
200  do i = is, ie
201  if (hecmesh%mpc%mpc_dof(i) > hecmat%NDOF) cycle outer
202  enddo
203 
204  a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
205 
206 
207  do i = is, ie
208  inod = hecmesh%mpc%mpc_item(i)
209 
210  is_internal_i = (hecmesh%node_ID(2*inod) == hecmw_comm_get_rank())
211 
212  idof = hecmesh%mpc%mpc_dof(i)
213  ai = hecmesh%mpc%mpc_val(i)
214  factor = ai * a1_2inv
215 
216  do j = is, ie
217  jnod = hecmesh%mpc%mpc_item(j)
218 
219  is_internal_j = (hecmesh%node_ID(2*jnod) == hecmw_comm_get_rank())
220  if (.not. (is_internal_i .or. is_internal_j)) cycle
221 
222  jdof = hecmesh%mpc%mpc_dof(j)
223  aj = hecmesh%mpc%mpc_val(j)
224 
225  call hecmw_mat_add_dof(hecmat, inod, idof, jnod, jdof, aj*factor*alpha)
226  enddo
227 
228  enddo
229  enddo outer
230 
231  call hecmw_mat_set_penalized(hecmat, 1)
232 
233  end subroutine hecmw_mat_ass_equation
234 
235 
236  subroutine hecmw_mat_ass_equation_rhs ( hecMESH, hecMAT )
237  type (hecmwst_matrix), target :: hecmat
238  type (hecmwst_local_mesh) :: hecmesh
239  !** Local variables
240  real(kind=kreal) :: alpha, a1_2inv, ai, factor, ci
241  integer(kind=kint) :: ndof, impc, is, ie, i, inod, idof
242 
243  if( hecmw_mat_get_penalized_b(hecmat) == 1) return
244 
245  alpha = hecmw_mat_get_penalty_alpha(hecmat)
246  if (alpha <= 0.0) stop "ERROR: penalty applied on vector before matrix"
247 
248  ndof = hecmat%NDOF
249 
250  outer: do impc = 1, hecmesh%mpc%n_mpc
251  is = hecmesh%mpc%mpc_index(impc-1) + 1
252  ie = hecmesh%mpc%mpc_index(impc)
253 
254  do i = is, ie
255  if (hecmesh%mpc%mpc_dof(i) > ndof) cycle outer
256  enddo
257 
258  a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
259 
260  do i = is, ie
261  inod = hecmesh%mpc%mpc_item(i)
262 
263  idof = hecmesh%mpc%mpc_dof(i)
264  ai = hecmesh%mpc%mpc_val(i)
265  factor = ai * a1_2inv
266 
267  ci = hecmesh%mpc%mpc_const(impc)
268  !$omp atomic
269  hecmat%B(ndof*(inod-1)+idof) = hecmat%B(ndof*(inod-1)+idof) + ci*factor*alpha
270  enddo
271  enddo outer
272 
273  call hecmw_mat_set_penalized_b(hecmat, 1)
274 
275  end subroutine hecmw_mat_ass_equation_rhs
276 
277 
278  subroutine hecmw_mat_add_dof(hecMAT, inod, idof, jnod, jdof, val)
279  type (hecmwst_matrix) :: hecmat
280  integer(kind=kint) :: inod, idof, jnod, jdof
281  real(kind=kreal) :: val
282  !** Local variables
283  integer(kind=kint) :: ndof, is, ie, k, idx
284 
285  ndof = hecmat%NDOF
286  if (inod < jnod) then
287  is = hecmat%indexU(inod-1)+1
288  ie = hecmat%indexU(inod)
289  k = hecmw_array_search_i(hecmat%itemU, is, ie, jnod)
290 
291  if (k < is .or. ie < k) then
292  write(*,*) '###ERROR### : cannot find connectivity (3)'
293  write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
295  return
296  endif
297 
298  idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
299  !$omp atomic
300  hecmat%AU(idx) = hecmat%AU(idx) + val
301 
302  else if (inod > jnod) then
303  is = hecmat%indexL(inod-1)+1
304  ie = hecmat%indexL(inod)
305  k = hecmw_array_search_i(hecmat%itemL, is, ie, jnod)
306 
307  if (k < is .or. ie < k) then
308  write(*,*) '###ERROR### : cannot find connectivity (4)'
309  write(*,*) ' myrank = ', hecmw_comm_get_rank(), ', inod = ', inod, ', jnod = ', jnod
311  return
312  endif
313 
314  idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
315  !$omp atomic
316  hecmat%AL(idx) = hecmat%AL(idx) + val
317 
318  else
319  idx = ndof**2 * (inod - 1) + ndof * (idof - 1) + jdof
320  !$omp atomic
321  hecmat%D(idx) = hecmat%D(idx) + val
322  endif
323 
324  end subroutine hecmw_mat_add_dof
325 
326  !C
327  !C***
328  !C*** MAT_ASS_BC
329  !C***
330  !C
331  subroutine hecmw_mat_ass_bc(hecMAT, inode, idof, RHS, conMAT)
332  type (hecmwst_matrix) :: hecmat
333  integer(kind=kint) :: inode, idof
334  real(kind=kreal) :: rhs, val
335  type (hecmwst_matrix),optional :: conmat
336  integer(kind=kint) :: ndof, in, i, ii, iii, ndof2, k, is, ie, iis, iie, ik, idx
337 
338  ndof = hecmat%NDOF
339  if( ndof < idof ) return
340 
341  !C-- DIAGONAL block
342 
343  hecmat%B(ndof*inode-(ndof-idof)) = rhs
344  if(present(conmat)) conmat%B(ndof*inode-(ndof-idof)) = 0.0d0
345  ndof2 = ndof*ndof
346  ii = ndof2 - idof
347 
348  do i = ndof-1,0,-1
349  if( i .NE. ndof-idof ) then
350  idx = ndof*inode-i
351  val = hecmat%D(ndof2*inode-ii)*rhs
352  !$omp atomic
353  hecmat%B(idx) = hecmat%B(idx) - val
354  if(present(conmat)) then
355  val = conmat%D(ndof2*inode-ii)*rhs
356  !$omp atomic
357  conmat%B(idx) = conmat%B(idx) - val
358  endif
359  endif
360  ii = ii - ndof
361  end do
362 
363  !*Set diagonal row to zero
364  ii = ndof2-1 - (idof-1)*ndof
365 
366  do i = 0, ndof - 1
367  hecmat%D(ndof2*inode-ii+i)=0.d0
368  if(present(conmat)) conmat%D(ndof2*inode-ii+i)=0.d0
369  end do
370 
371  !*Set diagonal column to zero
372  ii = ndof2 - idof
373  do i = 1, ndof
374  if( i.NE.idof ) then
375  hecmat%D(ndof2*inode-ii) = 0.d0
376  if(present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
377  else
378  hecmat%D(ndof2*inode-ii) = 1.d0
379  if(present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
380  endif
381  ii = ii - ndof
382  end do
383 
384  !C-- OFF-DIAGONAL blocks
385 
386  ii = ndof2-1 - (idof-1)*ndof
387  is = hecmat%indexL(inode-1) + 1
388  ie = hecmat%indexL(inode )
389 
390  do k= is, ie
391 
392  !*row (left)
393  do i = 0, ndof - 1
394  hecmat%AL(ndof2*k-ii+i) = 0.d0
395  if(present(conmat)) conmat%AL(ndof2*k-ii+i) = 0.d0
396  end do
397 
398  !*column (upper)
399  in = hecmat%itemL(k)
400  iis = hecmat%indexU(in-1) + 1
401  iie = hecmat%indexU(in )
402  do ik = iis, iie
403  if (hecmat%itemU(ik) .eq. inode) then
404  iii = ndof2 - idof
405  do i = ndof-1,0,-1
406  idx = ndof*in-i
407  val = hecmat%AU(ndof2*ik-iii)*rhs
408  !$omp atomic
409  hecmat%B(idx) = hecmat%B(idx) - val
410  hecmat%AU(ndof2*ik-iii)= 0.d0
411  if(present(conmat)) then
412  val = conmat%AU(ndof2*ik-iii)*rhs
413  !$omp atomic
414  conmat%B(idx) = conmat%B(idx) - val
415  conmat%AU(ndof2*ik-iii)= 0.d0
416  endif
417  iii = iii - ndof
418  end do
419  exit
420  endif
421  enddo
422 
423  enddo
424 
425  ii = ndof2-1 - (idof-1)*ndof
426  is = hecmat%indexU(inode-1) + 1
427  ie = hecmat%indexU(inode )
428 
429  do k= is, ie
430 
431  !*row (right)
432  do i = 0,ndof-1
433  hecmat%AU(ndof2*k-ii+i) = 0.d0
434  if(present(conmat)) conmat%AU(ndof2*k-ii+i) = 0.d0
435  end do
436 
437  !*column (lower)
438  in = hecmat%itemU(k)
439  iis = hecmat%indexL(in-1) + 1
440  iie = hecmat%indexL(in )
441  do ik= iis, iie
442  if (hecmat%itemL(ik) .eq. inode) then
443  iii = ndof2 - idof
444 
445  do i = ndof-1, 0, -1
446  idx = ndof*in-i
447  val = hecmat%AL(ndof2*ik-iii)*rhs
448  !$omp atomic
449  hecmat%B(idx) = hecmat%B(idx) - val
450  hecmat%AL(ndof2*ik-iii) = 0.d0
451  if(present(conmat)) then
452  val = conmat%AL(ndof2*ik-iii)*rhs
453  !$omp atomic
454  conmat%B(idx) = conmat%B(idx) - val
455  conmat%AL(ndof2*ik-iii) = 0.d0
456  endif
457  iii = iii - ndof
458  end do
459  exit
460  endif
461  enddo
462 
463  enddo
464  !*End off - diagonal blocks
465 
466  call hecmw_cmat_ass_bc(hecmat, inode, idof, rhs)
467 
468  end subroutine hecmw_mat_ass_bc
469 
470  !C
471  !C***
472  !C*** MAT_ASS_CONTACT
473  !C***
474  !C
475  subroutine hecmw_mat_ass_contact(hecMAT, nn, nodLOCAL, stiffness)
476  type (hecmwst_matrix) :: hecmat
477  integer(kind=kint) :: nn
478  integer(kind=kint) :: nodlocal(:)
479  real(kind=kreal) :: stiffness(:, :)
480  !** Local variables
481  integer(kind=kint) :: ndof, inod_e, jnod_e, inod, jnod
482  real(kind=kreal) :: a(3,3)
483 
484  ndof = hecmat%NDOF
485  if( ndof .ne. 3 ) then
486  write(*,*) '###ERROR### : ndof=',ndof,'; contact matrix supports only ndof==3'
488  return
489  endif
490 
491  do inod_e = 1, nn
492  inod = nodlocal(inod_e)
493  do jnod_e = 1, nn
494  jnod = nodlocal(jnod_e)
495  !***** Add components
496  call stf_get_block(stiffness, ndof, inod_e, jnod_e, a)
497  call hecmw_cmat_add(hecmat%cmat, inod, jnod, a)
498  enddo
499  enddo
500  call hecmw_cmat_pack(hecmat%cmat)
501 
502  end subroutine hecmw_mat_ass_contact
503 
504 end module hecmw_matrix_ass
subroutine, public hecmw_mat_ass_bc(hecMAT, inode, idof, RHS, conMAT)
integer(kind=kint) function, public hecmw_array_search_i(array, is, iE, ival)
subroutine, public hecmw_mat_add_dof(hecMAT, inod, idof, jnod, jdof, val)
subroutine, public hecmw_mat_ass_equation_rhs(hecMESH, hecMAT)
subroutine, public stf_get_block(stiffness, ndof, inod, jnod, a)
subroutine, public hecmw_mat_ass_contact(hecMAT, nn, nodLOCAL, stiffness)
subroutine, public hecmw_mat_ass_elem(hecMAT, nn, nodLOCAL, stiffness)
subroutine, public hecmw_mat_ass_equation(hecMESH, hecMAT)
subroutine, public hecmw_mat_add_node(hecMAT, inod, jnod, a)
subroutine, public hecmw_cmat_pack(cmat)
subroutine, public hecmw_cmat_ass_bc(hecMAT, inode, idof, RHS)
subroutine, public hecmw_cmat_add(cmat, i, j, val)
real(kind=kreal) function, public hecmw_mat_diag_max(hecMAT, hecMESH)
integer(kind=kint) function, public hecmw_mat_get_penalized(hecMAT)
subroutine, public hecmw_mat_set_penalized_b(hecMAT, penalized_b)
integer(kind=kint) function, public hecmw_mat_get_penalized_b(hecMAT)
real(kind=kreal) function, public hecmw_mat_get_penalty_alpha(hecMAT)
subroutine, public hecmw_mat_set_penalty_alpha(hecMAT, alpha)
subroutine, public hecmw_mat_set_penalized(hecMAT, penalized)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint) function hecmw_comm_get_comm()
integer(kind=4), parameter kreal
integer(kind=kint) function hecmw_comm_get_rank()
subroutine hecmw_abort(comm)