FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
20 public :: hecmw_mat_add_dof
21 public :: hecmw_mat_ass_bc
22 public :: hecmw_mat_ass_contact
23 public :: stf_get_block
24
25contains
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
504end module hecmw_matrix_ass
subroutine, public hecmw_mat_ass_contact(hecmat, nn, nodlocal, stiffness)
subroutine, public hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiffness)
subroutine, public hecmw_mat_add_node(hecmat, inod, jnod, a)
subroutine, public hecmw_mat_ass_bc(hecmat, inode, idof, rhs, conmat)
subroutine, public hecmw_mat_ass_equation_rhs(hecmesh, hecmat)
subroutine, public hecmw_mat_add_dof(hecmat, inod, idof, jnod, jdof, val)
integer(kind=kint) function, public hecmw_array_search_i(array, is, ie, ival)
subroutine, public stf_get_block(stiffness, ndof, inod, jnod, a)
subroutine, public hecmw_mat_ass_equation(hecmesh, hecmat)
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)
subroutine, public hecmw_mat_set_penalized(hecmat, penalized)
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_diag_max(hecmat, hecmesh)
real(kind=kreal) function, public hecmw_mat_get_penalty_alpha(hecmat)
integer(kind=kint) function, public hecmw_mat_get_penalized(hecmat)
subroutine, public hecmw_mat_set_penalty_alpha(hecmat, alpha)
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)