34 integer(kind=kint) :: nn
35 integer(kind=kint) :: nodlocal(:)
36 real(kind=
kreal) :: stiffness(:, :)
38 integer(kind=kint) :: ndof, inod_e, jnod_e, inod, jnod
39 real(kind=
kreal) :: a(6,6)
44 inod = nodlocal(inod_e)
46 jnod = nodlocal(jnod_e)
56 real(kind=
kreal) :: stiffness(:, :), a(:, :)
57 integer(kind=kint) :: ndof, inod, jnod
59 integer(kind=kint) :: row_offset, col_offset, i, j
61 row_offset = ndof*(inod-1)
64 col_offset = ndof*(jnod-1)
67 a(i, j) = stiffness(i + row_offset, j + col_offset)
75 integer(kind=kint) :: inod, jnod
76 real(kind=
kreal) :: a(:, :)
78 integer(kind=kint) :: ndof, is, ie, k, idx_base, idx, idof, jdof
83 is = hecmat%indexU(inod-1)+1
84 ie = hecmat%indexU(inod)
87 if (k < is .or. ie < k)
then
88 write(*,*)
'###ERROR### : cannot find connectivity (1)'
93 idx_base = ndof**2 * (k-1)
98 hecmat%AU(idx) = hecmat%AU(idx) + a(idof, jdof)
100 idx_base = idx_base + ndof
103 else if (inod > jnod)
then
104 is = hecmat%indexL(inod-1)+1
105 ie = hecmat%indexL(inod)
108 if (k < is .or. ie < k)
then
109 write(*,*)
'###ERROR### : cannot find connectivity (2)'
114 idx_base = ndof**2 * (k-1)
117 idx = idx_base + jdof
119 hecmat%AL(idx) = hecmat%AL(idx) + a(idof, jdof)
121 idx_base = idx_base + ndof
125 idx_base = ndof**2 * (inod - 1)
128 idx = idx_base + jdof
130 hecmat%D(idx) = hecmat%D(idx) + a(idof, jdof)
132 idx_base = idx_base + ndof
140 integer(kind=kint) :: array(:)
141 integer(kind=kint) :: is, ie, ival
143 integer(kind=kint) :: left, right, center, cval
148 if (left > right)
then
153 center = (left + right) / 2
156 if (ival < cval)
then
158 else if (cval < ival)
then
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
188 penalty => hecmat%Rarray(11)
190 if (penalty < 0.0) stop
"ERROR: negative penalty"
191 if (penalty < 1.0)
write(*,*)
"WARNING: penalty ", penalty,
" smaller than 1"
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)
201 if (hecmesh%mpc%mpc_dof(i) > hecmat%NDOF) cycle outer
204 a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
208 inod = hecmesh%mpc%mpc_item(i)
212 idof = hecmesh%mpc%mpc_dof(i)
213 ai = hecmesh%mpc%mpc_val(i)
214 factor = ai * a1_2inv
217 jnod = hecmesh%mpc%mpc_item(j)
220 if (.not. (is_internal_i .or. is_internal_j)) cycle
222 jdof = hecmesh%mpc%mpc_dof(j)
223 aj = hecmesh%mpc%mpc_val(j)
240 real(kind=
kreal) :: alpha, a1_2inv, ai, factor, ci
241 integer(kind=kint) :: ndof, impc, is, ie, i, inod, idof
246 if (alpha <= 0.0) stop
"ERROR: penalty applied on vector before matrix"
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)
255 if (hecmesh%mpc%mpc_dof(i) > ndof) cycle outer
258 a1_2inv = 1.0 / hecmesh%mpc%mpc_val(is)**2
261 inod = hecmesh%mpc%mpc_item(i)
263 idof = hecmesh%mpc%mpc_dof(i)
264 ai = hecmesh%mpc%mpc_val(i)
265 factor = ai * a1_2inv
267 ci = hecmesh%mpc%mpc_const(impc)
269 hecmat%B(ndof*(inod-1)+idof) = hecmat%B(ndof*(inod-1)+idof) + ci*factor*alpha
280 integer(kind=kint) :: inod, idof, jnod, jdof
281 real(kind=
kreal) :: val
283 integer(kind=kint) :: ndof, is, ie, k, idx
286 if (inod < jnod)
then
287 is = hecmat%indexU(inod-1)+1
288 ie = hecmat%indexU(inod)
291 if (k < is .or. ie < k)
then
292 write(*,*)
'###ERROR### : cannot find connectivity (3)'
298 idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
300 hecmat%AU(idx) = hecmat%AU(idx) + val
302 else if (inod > jnod)
then
303 is = hecmat%indexL(inod-1)+1
304 ie = hecmat%indexL(inod)
307 if (k < is .or. ie < k)
then
308 write(*,*)
'###ERROR### : cannot find connectivity (4)'
314 idx = ndof**2 * (k-1) + ndof * (idof-1) + jdof
316 hecmat%AL(idx) = hecmat%AL(idx) + val
319 idx = ndof**2 * (inod - 1) + ndof * (idof - 1) + jdof
321 hecmat%D(idx) = hecmat%D(idx) + val
333 integer(kind=kint) :: inode, idof
334 real(kind=
kreal) :: rhs, val
336 integer(kind=kint) :: ndof, in, i, ii, iii, ndof2, k, is, ie, iis, iie, ik, idx
339 if( ndof < idof )
return
343 hecmat%B(ndof*inode-(ndof-idof)) = rhs
344 if(
present(conmat)) conmat%B(ndof*inode-(ndof-idof)) = 0.0d0
349 if( i .NE. ndof-idof )
then
351 val = hecmat%D(ndof2*inode-ii)*rhs
353 hecmat%B(idx) = hecmat%B(idx) - val
354 if(
present(conmat))
then
355 val = conmat%D(ndof2*inode-ii)*rhs
357 conmat%B(idx) = conmat%B(idx) - val
364 ii = ndof2-1 - (idof-1)*ndof
367 hecmat%D(ndof2*inode-ii+i)=0.d0
368 if(
present(conmat)) conmat%D(ndof2*inode-ii+i)=0.d0
375 hecmat%D(ndof2*inode-ii) = 0.d0
376 if(
present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
378 hecmat%D(ndof2*inode-ii) = 1.d0
379 if(
present(conmat)) conmat%D(ndof2*inode-ii) = 0.d0
386 ii = ndof2-1 - (idof-1)*ndof
387 is = hecmat%indexL(inode-1) + 1
388 ie = hecmat%indexL(inode )
394 hecmat%AL(ndof2*k-ii+i) = 0.d0
395 if(
present(conmat)) conmat%AL(ndof2*k-ii+i) = 0.d0
400 iis = hecmat%indexU(in-1) + 1
401 iie = hecmat%indexU(in )
403 if (hecmat%itemU(ik) .eq. inode)
then
407 val = hecmat%AU(ndof2*ik-iii)*rhs
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
414 conmat%B(idx) = conmat%B(idx) - val
415 conmat%AU(ndof2*ik-iii)= 0.d0
425 ii = ndof2-1 - (idof-1)*ndof
426 is = hecmat%indexU(inode-1) + 1
427 ie = hecmat%indexU(inode )
433 hecmat%AU(ndof2*k-ii+i) = 0.d0
434 if(
present(conmat)) conmat%AU(ndof2*k-ii+i) = 0.d0
439 iis = hecmat%indexL(in-1) + 1
440 iie = hecmat%indexL(in )
442 if (hecmat%itemL(ik) .eq. inode)
then
447 val = hecmat%AL(ndof2*ik-iii)*rhs
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
454 conmat%B(idx) = conmat%B(idx) - val
455 conmat%AL(ndof2*ik-iii) = 0.d0
477 integer(kind=kint) :: nn
478 integer(kind=kint) :: nodlocal(:)
479 real(kind=
kreal) :: stiffness(:, :)
481 integer(kind=kint) :: ndof, inod_e, jnod_e, inod, jnod
482 real(kind=
kreal) :: a(3,3)
485 if( ndof .ne. 3 )
then
486 write(*,*)
'###ERROR### : ndof=',ndof,
'; contact matrix supports only ndof==3'
492 inod = nodlocal(inod_e)
494 jnod = nodlocal(jnod_e)
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)
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)
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)