FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
fstr_contact_elem_alag.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 !-------------------------------------------------------------------------------
6 module m_fstr_contact_elem_alag
7  use hecmw
8  use elementinfo
9  use mcontactdef
10  use msurfelement
11  use m_fstr_contact_geom
12  use m_fstr_contact_elem_common
13  use m_fstr_contact_smoothing
14  implicit none
15 
16  public :: getcontactstiffness_alag
17  public :: getcontactnodalforce_alag
18  public :: gettiedstiffness_alag
19  public :: gettiednodalforce_alag
20  public :: updatecontactmultiplier_alag
21  public :: get_unique_map
22  public :: getintgap
23  public :: build_group_tangent_basis
24  public :: gettangentslip
25  public :: group_return_mapping
26  public :: resolve_lambda_cur
27  public :: getcontactstiffness_alag_surfsurf
28  public :: getcontactnodalforce_alag_surfsurf
29 
30 contains
31 
32  subroutine getcontactstiffness_alag(cstate, tSurf, ele, mu, mut, fcoeff, symm, stiff, force, smoothing_type, edisp, iter, &
33  slvpos)
34 
35  type(tcontactstate), intent(inout) :: cstate
36  type(tsurfelement), intent(in) :: tsurf
37  real(kind=kreal), intent(in) :: ele(:,:)
38  real(kind=kreal), intent(in) :: mu, mut
39  real(kind=kreal), intent(in) :: fcoeff
40  logical, intent(in) :: symm
41  real(kind=kreal), intent(out) :: stiff(:,:)
42  real(kind=kreal), intent(out) :: force(:)
43  integer(kind=kint), optional, intent(in) :: smoothing_type
44  real(kind=kreal), optional, intent(in) :: edisp(:)
45  integer(kind=kint), intent(in) :: iter
46  real(kind=kreal), intent(in) :: slvpos(3)
47 
48  integer :: i, j, nnode
49  real(kind=kreal) :: bn(size(tsurf%nodes)*3+3), ht(2,size(tsurf%nodes)*3+3), gt(2,size(tsurf%nodes)*3+3)
50  real(kind=kreal) :: metric(2,2)
51  real(kind=kreal) :: a(2,2)
52  real(kind=kreal) :: alpha_proj, that_dir(2)
53  real(kind=kreal) :: k_fric(size(tsurf%nodes)*3+3,size(tsurf%nodes)*3+3)
54  real(kind=kreal) :: tmp_vec(2)
55  real(kind=kreal) :: dummy_force(size(tsurf%nodes)*3+3)
56  real(kind=kreal) :: eval_disp(size(tsurf%nodes)*3+3)
57  real(kind=kreal) :: curpos(size(tsurf%nodes)*3+3)
58  real(kind=kreal) :: lam_cone
59  real(kind=kreal) :: htt(size(tsurf%nodes)*3+3)
60  real(kind=kreal) :: tdir(2)
61  real(kind=kreal) :: invmetric(2,2), det, norm_lamt
62 
63  nnode = size(tsurf%nodes)
64 
65  ! Use common mapping routine to compute Bn, metric, Ht, Gt
66  call computecontactmaps_alag(cstate, tsurf, ele, bn, metric, ht, gt, smoothing_type)
67 
68  ! Normal stiffness: stiff = mu * Bn * Bn^T
69  do j = 1, nnode*3+3
70  do i = 1, nnode*3+3
71  stiff(i,j) = mu * bn(i) * bn(j)
72  enddo
73  enddo
74  force(1:nnode*3+3) = bn(:)
75 
76  ! frictional component
77  if( fcoeff /= 0.d0 ) then
78  ! Evaluate trial friction at current displacement for consistent tangent
79  if( present(edisp) ) then
80  eval_disp(1:nnode*3+3) = edisp(1:nnode*3+3)
81  else
82  eval_disp = 0.0d0
83  endif
84  ! Radius of the friction cone. By default it follows the normal force this element
85  ! actually applies, lambda_n + mu*g_n clipped at 0, which is the value
86  ! getContactNodalForce_Alag distributes; the multiplier alone lags that force by the
87  ! penalty term within a substep. With !CONTACT_ALGO, FRICTION_CONE=FROZEN (symm) the
88  ! radius stays at the multiplier of the last augmentation, which keeps the friction
89  ! terms symmetric and leaves the Coulomb condition to the augmentation loop.
90  if( symm ) then
91  lam_cone = cstate%multiplier(1)
92  else
93  curpos(1:3) = slvpos(1:3) + eval_disp(1:3)
94  do j = 1, nnode
95  curpos(j*3+1:j*3+3) = ele(1:3,j) + eval_disp(j*3+1:j*3+3)
96  enddo
97  lam_cone = max( 0.d0, cstate%multiplier(1) + mu*dot_product( bn(1:nnode*3+3), curpos(1:nnode*3+3) ) )
98  endif
99 
100  call computefrictionforce_alag(cstate, fcoeff, lam_cone, metric, &
101  ht, gt, eval_disp, nnode*3+3, dummy_force, &
102  mut, alpha=alpha_proj, that=that_dir)
103 
104  ! Friction tangent operator A in 2D metric space: K_fric = Ht^T * A * Ht
105  if( lam_cone <= 0.0d0 .or. alpha_proj <= 1.0d-20 ) then
106  ! No normal contact force: no friction contribution
107  a = 0.0d0
108  else if( alpha_proj >= 0.999d0 ) then
109  ! Stick: A = mu_t * M (exact)
110  a(1,1) = mut * metric(1,1)
111  a(1,2) = mut * metric(1,2)
112  a(2,1) = mut * metric(2,1)
113  a(2,2) = mut * metric(2,2)
114  else
115  ! Slip: switch tangent by NR iteration count for stability.
116  ! iter <= 2: A = alpha*mu_t*M (stable, no directional correction)
117  ! iter >= 3: A = alpha*mu_t*(M - t_hat x t_hat) (consistent tangent)
118  ! The first two NR steps use M to let t_hat stabilize; after that the consistent tangent is used to regain quadratic convergence.
119  if( iter <= 2 ) then
120  a(1,1) = alpha_proj * mut * metric(1,1)
121  a(1,2) = alpha_proj * mut * metric(1,2)
122  a(2,1) = alpha_proj * mut * metric(2,1)
123  a(2,2) = alpha_proj * mut * metric(2,2)
124  else
125  a(1,1) = alpha_proj * mut * (metric(1,1) - that_dir(1)*that_dir(1))
126  a(1,2) = alpha_proj * mut * (metric(1,2) - that_dir(1)*that_dir(2))
127  a(2,1) = alpha_proj * mut * (metric(2,1) - that_dir(2)*that_dir(1))
128  a(2,2) = alpha_proj * mut * (metric(2,2) - that_dir(2)*that_dir(2))
129  endif
130  endif
131 
132  ! Compute friction stiffness: K_fric = Ht^T * A * Ht (consistent tangent)
133  do j = 1, nnode*3+3
134  tmp_vec = matmul(a, ht(1:2,j))
135  do i = 1, nnode*3+3
136  k_fric(i,j) = dot_product(ht(1:2,i), tmp_vec)
137  enddo
138  enddo
139 
140  ! Coupling block from the radius following the normal force. On the slip branch the
141  ! friction force is f_t = R*tdir with R = fcoeff*lam_cone, and R varies with u through
142  ! g_n, so d(Ht^T f_t)/du gains Ht^T tdir * dR/du = fcoeff*mu * (Ht^T tdir) (x) Bn.
143  ! A stuck node does not use the radius (f_t is the full trial), so the block belongs to
144  ! the slip branch only, and it vanishes where the clip at 0 is active (lam_cone = 0).
145  ! Rows are a slip direction and columns a normal map, so the block is unsymmetric: the
146  ! caller has to set the linear solver up for a general matrix (see fstr_Newton_contactALag).
147  if( .not.symm .and. lam_cone > 0.0d0 .and. alpha_proj > 1.0d-20 .and. alpha_proj < 0.999d0 ) then
148  ! The direction is the tangential multiplier, not the trial direction that_dir, so
149  ! that it stays fixed inside the augmentation step; at the fixed point of the
150  ! augmentation the two coincide and the tangent is still the consistent one.
151  tdir(1:2) = that_dir(1:2)
152  det = metric(1,1)*metric(2,2) - metric(1,2)*metric(2,1)
153  if( abs(det) > 1.0d-20 ) then
154  invmetric(1,1) = metric(2,2)/det
155  invmetric(2,2) = metric(1,1)/det
156  invmetric(1,2) = -metric(1,2)/det
157  invmetric(2,1) = -metric(2,1)/det
158  tmp_vec(1:2) = matmul( invmetric(1:2,1:2), cstate%multiplier(2:3) )
159  norm_lamt = dsqrt( dot_product( cstate%multiplier(2:3), tmp_vec(1:2) ) )
160  if( norm_lamt > 1.0d-20 ) tdir(1:2) = cstate%multiplier(2:3) / norm_lamt
161  endif
162  htt(1:nnode*3+3) = matmul( transpose(ht(1:2,1:nnode*3+3)), tdir(1:2) )
163  do j = 1, nnode*3+3
164  do i = 1, nnode*3+3
165  k_fric(i,j) = k_fric(i,j) + fcoeff * mu * htt(i) * bn(j)
166  enddo
167  enddo
168  endif
169 
170  stiff(1:nnode*3+3,1:nnode*3+3) = stiff(1:nnode*3+3,1:nnode*3+3) + k_fric(1:nnode*3+3,1:nnode*3+3)
171  endif
172 
173  end subroutine getcontactstiffness_alag
174 
175  subroutine get_unique_map(sSurf, maplist, master_idxs, unique_count)
176  type(tcontactsurf) :: ssurf
177  integer(kind=kint), intent(out) :: maplist(:), master_idxs(:)
178  integer(kind=kint), intent(out) :: unique_count
179  integer(kind=kint) :: tmp(max_n_intp)
180  integer(kind=kint) :: i, j, n_intp, ctsurf
181  logical :: found
182 
183  n_intp = ssurf%n_intp
184  maplist = 0
185 
186  unique_count = 0
187 
188  do i = 1, n_intp
189  if( ssurf%states(i)%state == contactfree ) cycle
190  ctsurf = ssurf%states(i)%surface
191  found = .false.
192  ! Search existing groups by master surface index
193  do j = 1, unique_count
194  if (tmp(j) == ctsurf) then
195  maplist(i) = j
196  found = .true.
197  exit
198  endif
199  enddo
200  if (.not. found) then
201  unique_count = unique_count + 1
202  tmp(unique_count) = ctsurf
203  maplist(i) = unique_count
204  endif
205  enddo
206 
207  master_idxs(1:unique_count) = tmp(1:unique_count)
208 
209  end subroutine get_unique_map
210 
222  subroutine getintgap(slave_surf, master, coord, disp, ddisp, &
223  unique_count, maplist, master_idxs, &
224  Snode, Nsnode, gapwnode)
225  type(tcontactsurf) :: slave_surf
226  type(tsurfelement) :: master(:)
227  real(kind=kreal), intent(in) :: coord(:), disp(:), ddisp(:)
228  integer(kind=kint), intent(in) :: unique_count
229  integer(kind=kint), intent(in) :: maplist(:), master_idxs(:)
230  real(kind=kreal), intent(out) :: snode(:,:), nsnode(:,:,:), gapwnode(:,:)
231 
232  integer(kind=kint) :: i, j, g, a, nnode_s, nnode_m, etype, slave, n_intp, ctsurf, nd
233  integer(kind=kint) :: ndlocal(l_max_surface_node+1)
234  real(kind=kreal) :: snode_pos(3,4), weight(max_n_intp)
235  real(kind=kreal) :: ncoord(2), shapefunc_s(4), shapefunc_m(4), direction(3)
236  real(kind=kreal) :: curr_pos(24)
237 
238  nnode_s = size(slave_surf%nodes)
239 
240  snode = 0.d0
241  nsnode = 0.d0
242  gapwnode = 0.d0
243 
244  ! Slave node positions at start of substep (coord + disp), used for IP weights.
245  snode_pos = 0.d0
246  do i = 1, nnode_s
247  slave = slave_surf%nodes(i)
248  snode_pos(:,i) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave)
249  enddo
250  n_intp = slave_surf%n_intp
251  weight = 0.d0
252  call get_intp_weights(slave_surf%etype, nnode_s, n_intp, snode_pos, weight(1:n_intp))
253 
254  ! Accumulate the per-node-within-group constraint accumulator (Nsnode source) and area.
255  do i = 1, n_intp
256  if( slave_surf%states(i)%state == contactfree ) cycle
257  ctsurf = slave_surf%states(i)%surface
258  etype = master(ctsurf)%etype
259  nnode_m = size(master(ctsurf)%nodes)
260  direction = slave_surf%states(i)%direction(1:3)
261  call getintpoint4ss(slave_surf%etype, i, ncoord, n_intp, shapefunc_s)
262  call getshapefunc(etype, slave_surf%states(i)%lpos(1:2), shapefunc_m)
263  g = maplist(i)
264  ! Per-node-within-group: weight the whole IP constraint by the slave node shape function N_s(a).
265  do a = 1, nnode_s
266  snode(g,a) = snode(g,a) + shapefunc_s(a)*weight(i)
267  do j = 1, nnode_s
268  nsnode(g,a,3*j-2:3*j) = nsnode(g,a,3*j-2:3*j) &
269  + shapefunc_s(a)*shapefunc_s(j)*weight(i)*direction(1:3)
270  enddo
271  do j = nnode_s+1, nnode_s+nnode_m
272  nsnode(g,a,3*j-2:3*j) = nsnode(g,a,3*j-2:3*j) &
273  - shapefunc_s(a)*shapefunc_m(j-nnode_s)*weight(i)*direction(1:3)
274  enddo
275  enddo
276  enddo
277 
278  ! Per-node Nsnode (=ANnode/Snode) and weighted gap, at end of substep.
279  do g = 1, unique_count
280  ctsurf = master_idxs(g)
281  nnode_m = size(master(ctsurf)%nodes)
282  ndlocal(1:nnode_s) = slave_surf%nodes(1:nnode_s)
283  ndlocal(nnode_s+1:nnode_s+nnode_m) = master(ctsurf)%nodes(1:nnode_m)
284  do j = 1, nnode_s + nnode_m
285  nd = ndlocal(j)
286  curr_pos(3*j-2:3*j) = coord(3*nd-2:3*nd) + disp(3*nd-2:3*nd) + ddisp(3*nd-2:3*nd)
287  enddo
288  do a = 1, nnode_s
289  ! gapwnode uses the un-normalized accumulator (= ANnode . curr_pos); the residual/aug add mu*gapwnode.
290  gapwnode(g,a) = dot_product(nsnode(g,a,1:(nnode_s+nnode_m)*3), curr_pos(1:(nnode_s+nnode_m)*3))
291  if( snode(g,a) > 0.d0 ) then
292  nsnode(g,a,1:(nnode_s+nnode_m)*3) = nsnode(g,a,1:(nnode_s+nnode_m)*3) / snode(g,a)
293  else
294  nsnode(g,a,1:(nnode_s+nnode_m)*3) = 0.d0
295  endif
296  enddo
297  enddo
298 
299  end subroutine getintgap
300 
308  subroutine build_group_tangent_basis(n_hat, t1, t2)
309  real(kind=kreal), intent(in) :: n_hat(3)
310  real(kind=kreal), intent(out) :: t1(3), t2(3)
311 
312  integer(kind=kint) :: k
313  real(kind=kreal) :: v(3), vn, dotk
314 
315  ! choose the axis least aligned with n_hat for numerical stability
316  k = 1
317  if( abs(n_hat(2)) < abs(n_hat(k)) ) k = 2
318  if( abs(n_hat(3)) < abs(n_hat(k)) ) k = 3
319 
320  ! v = e_k - (e_k . n_hat) n_hat (projection of e_k onto the tangent plane)
321  dotk = n_hat(k)
322  v(1:3) = -dotk * n_hat(1:3)
323  v(k) = v(k) + 1.0d0
324 
325  vn = sqrt( v(1)*v(1) + v(2)*v(2) + v(3)*v(3) )
326  t1(1:3) = v(1:3) / vn
327 
328  ! t2 = n_hat x t1 (completes the right-handed orthonormal frame)
329  t2(1) = n_hat(2)*t1(3) - n_hat(3)*t1(2)
330  t2(2) = n_hat(3)*t1(1) - n_hat(1)*t1(3)
331  t2(3) = n_hat(1)*t1(2) - n_hat(2)*t1(1)
332  end subroutine build_group_tangent_basis
333 
344  subroutine gettangentslip(slave_surf, master, coord, disp, ddisp, &
345  unique_count, maplist, master_idxs, Sigma_node, nacc_node)
346  type(tcontactsurf) :: slave_surf
347  type(tsurfelement) :: master(:)
348  real(kind=kreal), intent(in) :: coord(:), disp(:), ddisp(:)
349  integer(kind=kint), intent(in) :: unique_count
350  integer(kind=kint), intent(in) :: maplist(:), master_idxs(:)
351  real(kind=kreal), intent(out) :: sigma_node(:,:,:)
352  real(kind=kreal), intent(out) :: nacc_node(:,:,:)
353 
354  integer(kind=kint) :: i, j, g, a, nnode_s, nnode_m, etype, slave, n_intp, ctsurf, nd
355  real(kind=kreal) :: snode_pos(3,4), weight(max_n_intp)
356  real(kind=kreal) :: ncoord(2), shapefunc_s(4), shapefunc_m(4), direction(3)
357  real(kind=kreal) :: du_rel(3), du_master(3)
358 
359  nnode_s = size(slave_surf%nodes)
360  sigma_node = 0.d0
361  nacc_node = 0.d0
362 
363  ! Slave node positions at start of substep (coord + disp), used for IP weights
364  ! (identical to getIntGap so the weights match the normal aggregation exactly).
365  snode_pos = 0.d0
366  do i = 1, nnode_s
367  slave = slave_surf%nodes(i)
368  snode_pos(:,i) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave)
369  enddo
370  n_intp = slave_surf%n_intp
371  weight = 0.d0
372  call get_intp_weights(slave_surf%etype, nnode_s, n_intp, snode_pos, weight(1:n_intp))
373 
374  do i = 1, n_intp
375  if( slave_surf%states(i)%state == contactfree ) cycle
376  ctsurf = slave_surf%states(i)%surface
377  etype = master(ctsurf)%etype
378  nnode_m = size(master(ctsurf)%nodes)
379  direction = slave_surf%states(i)%direction(1:3)
380  call getintpoint4ss(slave_surf%etype, i, ncoord, n_intp, shapefunc_s)
381  call getshapefunc(etype, slave_surf%states(i)%lpos(1:2), shapefunc_m)
382  g = maplist(i)
383 
384  ! Mortar relative displacement = Tm_mortar . edisp (edisp = ddisp).
385  ! slave block: + N_s(j) * ddisp(slave_j) ; master block: - N_m(j) * ddisp(master_j).
386  du_rel(1:3) = 0.d0
387  do j = 1, nnode_s
388  nd = slave_surf%nodes(j)
389  du_rel(1:3) = du_rel(1:3) + shapefunc_s(j) * ddisp(3*nd-2:3*nd)
390  enddo
391  du_master(1:3) = 0.d0
392  do j = 1, nnode_m
393  nd = master(ctsurf)%nodes(j)
394  du_master(1:3) = du_master(1:3) + shapefunc_m(j) * ddisp(3*nd-2:3*nd)
395  enddo
396  du_rel(1:3) = du_rel(1:3) - du_master(1:3)
397 
398  ! Per-node: weight the IP slip/normal by the slave node shape function N_s(a).
399  do a = 1, nnode_s
400  sigma_node(g, a, 1:3) = sigma_node(g, a, 1:3) + shapefunc_s(a) * weight(i) * du_rel(1:3)
401  nacc_node(g, a, 1:3) = nacc_node(g, a, 1:3) + shapefunc_s(a) * weight(i) * direction(1:3)
402  enddo
403  enddo
404  end subroutine gettangentslip
405 
413  subroutine group_return_mapping(lam_t_in, rho_t, Dxi, fcoeff, lam_n, eps_fric_band, &
414  lam_t_out, fric_state, alpha, that, update_state)
415  real(kind=kreal), intent(in) :: lam_t_in(2)
416  real(kind=kreal), intent(in) :: rho_t
417  real(kind=kreal), intent(in) :: dxi(2)
418  real(kind=kreal), intent(in) :: fcoeff
419  real(kind=kreal), intent(in) :: lam_n
420  real(kind=kreal), intent(in) :: eps_fric_band
421  real(kind=kreal), intent(out) :: lam_t_out(2)
422  integer(kind=kint), intent(inout) :: fric_state
423  real(kind=kreal), intent(out) :: alpha
424  real(kind=kreal), intent(out) :: that(2)
425  logical, optional, intent(in) :: update_state
426 
427  real(kind=kreal) :: trial(2), norm_trial, radius
428  logical :: is_stick, do_update
429 
430  do_update = .false.
431  if( present(update_state) ) do_update = update_state
432 
433  trial(1:2) = lam_t_in(1:2) + rho_t * dxi(1:2)
434  norm_trial = sqrt( trial(1)*trial(1) + trial(2)*trial(2) ) ! metric = I
435  radius = fcoeff * lam_n
436 
437  if( norm_trial > 1.0d-20 ) then
438  that(1:2) = trial(1:2) / norm_trial
439  else
440  that(1:2) = 0.0d0
441  endif
442 
443  if( lam_n <= 0.0d0 ) then
444  ! No normal force: cone collapses, no friction. Leave fric_state as-is.
445  lam_t_out(1:2) = 0.0d0
446  alpha = 0.0d0
447  else if( do_update ) then
448  ! Augmentation: this is the only place where the stick/slip state is (re)decided.
449  ! Hysteresis band: read the warm-start (previous) state and switch only when
450  ! ||trial|| crosses the asymmetric thresholds. Inside the band the previous
451  ! branch is kept, so marginal flips at the cone boundary cannot drive a limit
452  ! cycle across augmentations.
453  is_stick = ( fric_state == contactstick )
454  if( is_stick ) then
455  if( norm_trial > (1.0d0 + eps_fric_band) * radius ) is_stick = .false. ! genuine slip onset
456  else
457  if( norm_trial <= (1.0d0 - eps_fric_band) * radius ) is_stick = .true. ! genuine stick recovery
458  endif
459 
460  if( is_stick ) then
461  ! Stick branch: keep full trial multiplier.
462  lam_t_out(1:2) = trial(1:2)
463  fric_state = contactstick
464  alpha = 1.0d0
465  else
466  ! Slip branch: project onto cone surface.
467  alpha = radius / norm_trial
468  lam_t_out(1:2) = alpha * trial(1:2)
469  fric_state = contactslip
470  endif
471  else
472  ! Inner Newton-Raphson (read-only): pin the state frozen at the previous
473  ! augmentation, skipping the band re-judgement. fric_state is not written.
474  ! Same isolation as NTS-ALAG (update_multiplier=.false.): the stick/slip
475  ! branch cannot flip on the live trial inside the inner NR, so the residual
476  ! stays smooth and the per-iteration branch-flip limit cycle is broken.
477  if( fric_state == contactstick ) then
478  ! Frozen STICK: keep full trial multiplier (no cone projection).
479  lam_t_out(1:2) = trial(1:2)
480  alpha = 1.0d0
481  else
482  ! Frozen SLIP: project onto the live cone surface with min(1, radius/||trial||), the same
483  ! expression computeFrictionForce_ALag uses on the node-to-surface side. The projection
484  ! never scales a trial up: a node whose trial has come back inside the cone keeps the full
485  ! trial force, which is what the alpha >= 0.999 branch of the tangent is linearised about.
486  ! A node whose state was frozen SLIP at an augmentation where it carried no normal force
487  ! (the lam_n <= 0 branch leaves fric_state untouched) can come back with a zero trial, and
488  ! radius/0 would turn the whole residual into NaN. No trial force means no friction force,
489  ! which is what alpha = 0 gives.
490  if( norm_trial > 1.0d-20 ) then
491  alpha = min( 1.0d0, radius / norm_trial )
492  else
493  alpha = 0.0d0
494  endif
495  lam_t_out(1:2) = alpha * trial(1:2)
496  endif
497  endif
498  end subroutine group_return_mapping
499 
503  subroutine resolve_lambda_cur( surf, master_idxs, unique_count, nnode_s, sorted_idx, &
504  lambda_node, lam_t_cur, fric_state_cur )
505  type(tcontactsurf), intent(in) :: surf
506  integer(kind=kint), intent(in) :: master_idxs(:)
507  integer(kind=kint), intent(in) :: unique_count
508  integer(kind=kint), intent(in) :: nnode_s
509  integer(kind=kint), intent(out) :: sorted_idx(:)
510  real(kind=kreal), intent(out) :: lambda_node(:,:)
511  ! Optional friction warm-start: same reference rule as lambda_n (working -> begin ->
512  ! default), riding the same merge.
513  real(kind=kreal), intent(out), optional :: lam_t_cur(:,:,:)
514  integer(kind=kint), intent(out), optional :: fric_state_cur(:,:)
515  integer(kind=kint) :: r, j, tmp, ib, iw, g, mid
516  logical :: do_fric
517 
518  ! argsort master_idxs ascending (unique_count <= 27, insertion sort)
519  do r = 1, unique_count
520  sorted_idx(r) = r
521  enddo
522  do r = 2, unique_count
523  tmp = sorted_idx(r)
524  j = r - 1
525  do while( j >= 1 )
526  if( master_idxs(sorted_idx(j)) <= master_idxs(tmp) ) exit
527  sorted_idx(j+1) = sorted_idx(j)
528  j = j - 1
529  enddo
530  sorted_idx(j+1) = tmp
531  enddo
532 
533  do_fric = present(lam_t_cur) .and. present(fric_state_cur)
534 
535  ! 2-pointer merge over ascending masters / ascending begin / ascending working
536  ib = 1; iw = 1
537  do r = 1, unique_count
538  mid = master_idxs(sorted_idx(r))
539  do while( ib <= surf%lam_begin_n .and. surf%lam_begin_id(ib) < mid ); ib = ib + 1; enddo
540  do while( iw <= surf%lam_work_n .and. surf%lam_work_id(iw) < mid ); iw = iw + 1; enddo
541  g = sorted_idx(r)
542  if( iw <= surf%lam_work_n .and. surf%lam_work_id(iw) == mid ) then
543  lambda_node(1:nnode_s,g) = surf%lam_work_val(1:nnode_s,iw)
544  if( do_fric ) then
545  lam_t_cur(1:2,1:nnode_s,g) = surf%lam_work_t(1:2,1:nnode_s,iw)
546  fric_state_cur(1:nnode_s,g) = surf%lam_work_fstate(1:nnode_s,iw)
547  endif
548  else if( ib <= surf%lam_begin_n .and. surf%lam_begin_id(ib) == mid ) then
549  lambda_node(1:nnode_s,g) = surf%lam_begin_val(1:nnode_s,ib)
550  if( do_fric ) then
551  lam_t_cur(1:2,1:nnode_s,g) = surf%lam_begin_t(1:2,1:nnode_s,ib)
552  fric_state_cur(1:nnode_s,g) = surf%lam_begin_fstate(1:nnode_s,ib)
553  endif
554  else
555  lambda_node(1:nnode_s,g) = 0.d0
556  if( do_fric ) then
557  lam_t_cur(1:2,1:nnode_s,g) = 0.d0
558  fric_state_cur(1:nnode_s,g) = contactstick
559  endif
560  endif
561  enddo
562  end subroutine resolve_lambda_cur
563 
577  subroutine getcontactstiffness_alag_surfsurf( slave_surf, master, coord, disp, ddisp, &
578  mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, &
579  stiff_n, active_n, stiff_t, active_t )
580  type(tcontactsurf), intent(in) :: slave_surf
581  type(tsurfelement), intent(in) :: master(:)
582  real(kind=kreal), intent(in) :: coord(:)
583  real(kind=kreal), intent(in) :: disp(:)
584  real(kind=kreal), intent(in) :: ddisp(:)
585  real(kind=kreal), intent(in) :: mu, mut
586  real(kind=kreal), intent(in) :: fcoeff
587  logical, intent(in) :: symm
588  real(kind=kreal), intent(in) :: eps_fric_band
589  integer(kind=kint), intent(in) :: unique_count
590  integer(kind=kint), intent(in) :: maplist(:)
591  integer(kind=kint), intent(in) :: master_idxs(:)
592  real(kind=kreal), intent(out) :: stiff_n(:,:,:,:)
593  logical, intent(out) :: active_n(:,:)
594  real(kind=kreal), intent(out), optional :: stiff_t(:,:,:,:)
595  logical, intent(out), optional :: active_t(:,:)
596 
597  integer(kind=kint) :: g, a, j, k, nnode_m, nnode_s
598  integer(kind=kint), allocatable :: sorted_idx(:)
599  real(kind=kreal), allocatable :: snode(:,:), nsnode(:,:,:), gapwnode(:,:), lambda_node(:,:)
600  real(kind=kreal) :: ns(24)
601  ! --- friction consistent tangent ---
602  real(kind=kreal), allocatable :: lam_t_cur(:,:,:), nacc_node(:,:,:), sigma_node(:,:,:)
603  integer(kind=kint), allocatable :: fric_state_cur(:,:)
604  real(kind=kreal) :: nhat(3), t1(3), t2(3), nrm, dxi(2)
605  real(kind=kreal) :: alpha, that(2), lam_t_new(2), amat(2,2), t3d(3,2), m3(3,3)
606  real(kind=kreal) :: wb(l_max_surface_node+1)
607  real(kind=kreal) :: lam_cone, that3d(3)
608  integer(kind=kint) :: na, nb, fstate
609 
610  nnode_s = size(slave_surf%nodes)
611  allocate(snode(unique_count,nnode_s), nsnode(unique_count,nnode_s,24), gapwnode(unique_count,nnode_s))
612  call getintgap(slave_surf, master, coord, disp, ddisp, &
613  unique_count, maplist, master_idxs, &
614  snode, nsnode, gapwnode)
615 
616  allocate(sorted_idx(unique_count), lambda_node(nnode_s,unique_count))
617  stiff_n = 0.d0
618  active_n = .false.
619  if( fcoeff /= 0.d0 ) then
620  allocate(lam_t_cur(2,nnode_s,unique_count), fric_state_cur(nnode_s,unique_count))
621  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
622  lambda_node, lam_t_cur, fric_state_cur)
623  else
624  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
625  lambda_node)
626  endif
627 
628  ! ===== Normal stiffness =====
629  do g = 1, unique_count
630  nnode_m = size(master(master_idxs(g))%nodes)
631  do a = 1, nnode_s
632  ! ALag contact condition per node: augmented force must be positive
633  if( lambda_node(a,g)+mu*gapwnode(g,a) < 0.d0 ) cycle
634  active_n(a,g) = .true.
635  ns = 0.d0
636  ns(1:(nnode_s+nnode_m)*3) = nsnode(g, a, 1:(nnode_s+nnode_m)*3)
637  do j = 1, (nnode_s+nnode_m)*3
638  do k = 1, (nnode_s+nnode_m)*3
639  stiff_n(j,k,a,g) = mu*snode(g,a)*ns(j)*ns(k)
640  enddo
641  enddo
642  enddo
643  enddo
644 
645  if( fcoeff /= 0.d0 ) then
646  ! ===== Friction consistent tangent (per slave node) =====
647  ! Linearization of the per-node friction residual at the same live slip state:
648  ! K_a(b,c) = Snode(g,a) * Wbar(a,b) * Wbar(a,c) * M3_a, M3_a = T3d_a * A_a * T3d_a^T
649  ! with Wbar(a,b) = Nsnode(g,a,b).nhat_a, the same map as the residual back-distribution.
650  stiff_t = 0.d0
651  active_t = .false.
652  allocate(sigma_node(unique_count,nnode_s,3), nacc_node(unique_count,nnode_s,3))
653  call gettangentslip(slave_surf, master, coord, disp, ddisp, &
654  unique_count, maplist, master_idxs, sigma_node, nacc_node)
655  do g = 1, unique_count
656  nnode_m = size(master(master_idxs(g))%nodes)
657  do a = 1, nnode_s
658  if( lambda_node(a,g) <= 0.d0 ) cycle ! no per-node normal force -> no friction
659  ! Radius of the friction cone. With FRICTION_CONE=FROZEN it stays at the multiplier
660  ! of the last augmentation, which keeps the friction terms symmetric and leaves the
661  ! Coulomb condition to the augmentation loop; with !CONTACT_ALGO, FRICTION_CONE=FOLLOW it
662  ! follows the normal force this node actually applies, lambda_node + rho_n*gapwnode,
663  ! the same expression the residual distributes as nrlforce.
664  if( symm ) then
665  lam_cone = lambda_node(a,g)
666  else
667  lam_cone = lambda_node(a,g) + mu*gapwnode(g,a)
668  endif
669  nrm = sqrt( nacc_node(g,a,1)**2 + nacc_node(g,a,2)**2 + nacc_node(g,a,3)**2 )
670  if( nrm < 1.d-30 ) cycle
671  nhat(1:3) = nacc_node(g,a,1:3) / nrm
672  call build_group_tangent_basis(nhat, t1, t2)
673  ! Live per-node slip projection and read-only return mapping (writes only OUT args).
674  dxi(1) = dot_product(t1(1:3), sigma_node(g,a,1:3))
675  dxi(2) = dot_product(t2(1:3), sigma_node(g,a,1:3))
676  fstate = fric_state_cur(a,g)
677  call group_return_mapping(lam_t_cur(1:2,a,g), mut, dxi, fcoeff, lam_cone, &
678  eps_fric_band, lam_t_new, fstate, alpha, that)
679  ! 2D tangent operator A (same construction as getContactStiffness_Alag).
680  if( alpha <= 1.0d-20 ) then
681  amat = 0.d0
682  else if( alpha >= 0.999d0 ) then
683  amat = 0.d0
684  amat(1,1) = mut
685  amat(2,2) = mut
686  else
687  amat(1,1) = alpha * mut * (1.0d0 - that(1)*that(1))
688  amat(1,2) = alpha * mut * (-that(1)*that(2))
689  amat(2,1) = alpha * mut * (-that(2)*that(1))
690  amat(2,2) = alpha * mut * (1.0d0 - that(2)*that(2))
691  endif
692  t3d(1:3,1) = t1(1:3)
693  t3d(1:3,2) = t2(1:3)
694  m3 = matmul( matmul(t3d, amat), transpose(t3d) )
695 
696  ! Per-node averaged mortar weight of each node (= ANnode/Snode, recovered via nhat_a).
697  do na = 1, nnode_s + nnode_m
698  wb(na) = dot_product(nsnode(g,a,3*na-2:3*na), nhat(1:3))
699  enddo
700  active_t(a,g) = .true.
701  do nb = 1, nnode_s + nnode_m
702  do na = 1, nnode_s + nnode_m
703  do k = 1, 3
704  do j = 1, 3
705  stiff_t(3*na-3+j, 3*nb-3+k, a, g) = snode(g,a) * wb(na) * wb(nb) * m3(j,k)
706  enddo
707  enddo
708  enddo
709  enddo
710  ! Coupling block of a cone radius that follows the normal force. On the slip branch
711  ! the friction force is f_t = R*that3d with R = fcoeff*lam_cone, and R varies with u
712  ! through gapwnode: d(gapwnode(g,a))/du = Snode(g,a)*Nsnode(g,a,:), the map the normal
713  ! stiffness uses, so the residual -Wbar(a,b)*f_t gains
714  ! K_a(b,c) += fcoeff*rho_n*Snode(g,a) * Wbar(a,b)*that3d (x) Nsnode(g,a,c).
715  ! Rows are a slip direction and columns a normal map, so the block is unsymmetric and
716  ! the solver is set up for a general matrix (fstr_is_contactALag_symmetric). A stuck
717  ! node does not use the radius (f_t is the full trial), hence the slip-branch window,
718  ! the same one the consistent tangent above uses.
719  if( .not.symm .and. alpha > 1.0d-20 .and. alpha < 0.999d0 ) then
720  that3d(1:3) = that(1)*t1(1:3) + that(2)*t2(1:3)
721  do nb = 1, nnode_s + nnode_m
722  do na = 1, nnode_s + nnode_m
723  do k = 1, 3
724  do j = 1, 3
725  stiff_t(3*na-3+j, 3*nb-3+k, a, g) = stiff_t(3*na-3+j, 3*nb-3+k, a, g) &
726  + fcoeff * mu * snode(g,a) * wb(na) * that3d(j) * nsnode(g,a,3*nb-3+k)
727  enddo
728  enddo
729  enddo
730  enddo
731  endif
732  enddo
733  enddo
734 
735  deallocate(sigma_node, nacc_node, lam_t_cur, fric_state_cur)
736  endif
737 
738  deallocate(sorted_idx)
739  deallocate(snode, nsnode, gapwnode, lambda_node)
740  end subroutine getcontactstiffness_alag_surfsurf
741 
754  subroutine getcontactnodalforce_alag_surfsurf( purpose, slave_surf, master, coord, disp, ddisp, &
755  mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, &
756  ctNForce, active_n, ctTForce, active_t )
757  integer(kind=kint), intent(in) :: purpose
758  type(tcontactsurf), intent(in) :: slave_surf
759  type(tsurfelement), intent(in) :: master(:)
760  real(kind=kreal), intent(in) :: coord(:)
761  real(kind=kreal), intent(in) :: disp(:)
762  real(kind=kreal), intent(in) :: ddisp(:)
763  real(kind=kreal), intent(in) :: mu, mut
764  real(kind=kreal), intent(in) :: fcoeff
765  logical, intent(in) :: symm
766  real(kind=kreal), intent(in) :: eps_fric_band
767  integer(kind=kint), intent(in) :: unique_count
768  integer(kind=kint), intent(in) :: maplist(:)
769  integer(kind=kint), intent(in) :: master_idxs(:)
770  real(kind=kreal), intent(out) :: ctnforce(:,:,:)
771  logical, intent(out) :: active_n(:,:)
772  real(kind=kreal), intent(out), optional :: cttforce(:,:,:)
773  logical, intent(out), optional :: active_t(:,:)
774 
775  integer(kind=kint) :: g, a, j, nnode_m, nnode_s
776  integer(kind=kint), allocatable :: sorted_idx(:)
777  real(kind=kreal), allocatable :: snode(:,:), nsnode(:,:,:), gapwnode(:,:), lambda_node(:,:)
778  real(kind=kreal) :: nrlforce
779  real(kind=kreal) :: ns(24)
780  ! --- friction force back-distribution (live return mapping) ---
781  real(kind=kreal), allocatable :: lam_t_cur(:,:,:), nacc_node(:,:,:), sigma_node(:,:,:)
782  integer(kind=kint), allocatable :: fric_state_cur(:,:)
783  real(kind=kreal) :: nhat(3), t1(3), t2(3), nrm, fvec(3), wbar
784  real(kind=kreal) :: dxi(2), alpha, that(2), lam_t_new(2), lam_cone
785  integer(kind=kint) :: fstate
786 
787  nnode_s = size(slave_surf%nodes)
788  allocate(snode(unique_count,nnode_s), nsnode(unique_count,nnode_s,24), gapwnode(unique_count,nnode_s))
789  call getintgap(slave_surf, master, coord, disp, ddisp, &
790  unique_count, maplist, master_idxs, &
791  snode, nsnode, gapwnode)
792 
793  allocate(sorted_idx(unique_count), lambda_node(nnode_s,unique_count))
794  ctnforce = 0.d0
795  active_n = .false.
796  if( fcoeff /= 0.d0 ) then
797  allocate(lam_t_cur(2,nnode_s,unique_count), fric_state_cur(nnode_s,unique_count))
798  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
799  lambda_node, lam_t_cur, fric_state_cur)
800  else
801  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
802  lambda_node)
803  endif
804 
805  ! ===== Normal force: per-node back-distribution =====
806  ! nrlforce_a = lambda_node(a,g) + mu*gapwnode(g,a) for the residual, lambda_node(a,g) for output.
807  do g = 1, unique_count
808  nnode_m = size(master(master_idxs(g))%nodes)
809  do a = 1, nnode_s
810  nrlforce = lambda_node(a,g) + mu*gapwnode(g,a)
811  ! ALag contact condition per node: augmented force must be positive
812  if( nrlforce < 0.d0 ) cycle
813  active_n(a,g) = .true.
814  ns = 0.d0
815  ns(1:(nnode_s+nnode_m)*3) = nsnode(g, a, 1:(nnode_s+nnode_m)*3)
816  do j = 1, nnode_s + nnode_m
817  if( purpose == kctforresidual ) then
818  ctnforce(3*j-2:3*j,a,g) = -nrlforce*ns(3*j-2:3*j)
819  else if ( purpose == kctforoutput ) then
820  ! Output: multiplier only (converges to true contact force)
821  ctnforce(3*j-2:3*j,a,g) = -lambda_node(a,g)*ns(3*j-2:3*j)
822  end if
823  enddo
824  enddo
825  enddo
826 
827  if( fcoeff /= 0.d0 ) then
828  ! ===== Friction force (per slave node): live return mapping, back-distributed =====
829  ! trial = lam_t_warm(a) + rho_t*Dxi_live(a), projected onto the cone of radius lambda_node(a,g).
830  ! The return mapping is read-only here (the augmentation update is the sole writer of the
831  ! lambda_t / fric_state buffers). The resulting traction is distributed through the per-node
832  ! mortar weight Wbar(a,j) = Nsnode(g,a,j).nhat_a, mirroring the normal back-distribution.
833  ! Output (kctForOutput) keeps the frozen multiplier.
834  cttforce = 0.d0
835  active_t = .false.
836  allocate(sigma_node(unique_count,nnode_s,3), nacc_node(unique_count,nnode_s,3))
837  call gettangentslip(slave_surf, master, coord, disp, ddisp, &
838  unique_count, maplist, master_idxs, sigma_node, nacc_node)
839  do g = 1, unique_count
840  nnode_m = size(master(master_idxs(g))%nodes)
841  do a = 1, nnode_s
842  if( lambda_node(a,g) <= 0.d0 ) cycle ! no per-node normal force -> no friction
843  ! Cone radius: the frozen multiplier with FRICTION_CONE=FROZEN, the normal force this
844  ! node just applied above (nrlforce = lambda_node + rho_n*gapwnode) with FRICTION_CONE=FOLLOW,
845  ! the same radius the tangent uses (see getContactStiffness_Alag_SurfSurf). A negative
846  ! lam_cone reaches group_return_mapping as lam_n <= 0 and gives zero friction, which is
847  ! what the normal back-distribution above does with a negative nrlforce too.
848  if( symm ) then
849  lam_cone = lambda_node(a,g)
850  else
851  lam_cone = lambda_node(a,g) + mu*gapwnode(g,a)
852  endif
853  nrm = sqrt( nacc_node(g,a,1)**2 + nacc_node(g,a,2)**2 + nacc_node(g,a,3)**2 )
854  if( nrm < 1.d-30 ) cycle
855  nhat(1:3) = nacc_node(g,a,1:3) / nrm
856  call build_group_tangent_basis(nhat, t1, t2)
857  if( purpose == kctforresidual ) then
858  ! Live trial: project the live per-node mortar slip onto the per-node frame and return-map.
859  dxi(1) = dot_product(t1(1:3), sigma_node(g,a,1:3))
860  dxi(2) = dot_product(t2(1:3), sigma_node(g,a,1:3))
861  fstate = fric_state_cur(a,g)
862  call group_return_mapping(lam_t_cur(1:2,a,g), mut, dxi, fcoeff, lam_cone, &
863  eps_fric_band, lam_t_new, fstate, alpha, that)
864  fvec(1:3) = lam_t_new(1)*t1(1:3) + lam_t_new(2)*t2(1:3)
865  else
866  ! Output: frozen multiplier only (converges to the true friction force).
867  fvec(1:3) = lam_t_cur(1,a,g)*t1(1:3) + lam_t_cur(2,a,g)*t2(1:3)
868  end if
869 
870  active_t(a,g) = .true.
871  do j = 1, nnode_s + nnode_m
872  ! per-node averaged mortar weight of node j (= ANnode/Snode, recovered via nhat_a)
873  wbar = dot_product(nsnode(g,a,3*j-2:3*j), nhat(1:3))
874  cttforce(3*j-2:3*j,a,g) = -fvec(1:3) * wbar
875  enddo
876  enddo
877  enddo
878 
879  deallocate(sigma_node, nacc_node, lam_t_cur, fric_state_cur)
880  endif
881 
882  deallocate(sorted_idx)
883  deallocate(snode, nsnode, gapwnode, lambda_node)
884  end subroutine getcontactnodalforce_alag_surfsurf
885 
886  subroutine getcontactnodalforce_alag(ctState,tSurf,ndCoord,ndDu,mu,mut,fcoeff,symm,lagrange,ctNForce,ctTForce,cflag, &
887  smoothing_type)
888 
889  use msurfelement
890  type(tcontactstate) :: ctstate
891  type(tsurfelement) :: tsurf
892  integer(kind=kint) :: nnode
893  integer(kind=kint) :: j
894  real(kind=kreal), intent(in) :: mu, mut
895  real(kind=kreal) :: fcoeff
896  logical, intent(in) :: symm
897  real(kind=kreal) :: lagrange
898  real(kind=kreal) :: ndcoord(:), nddu(:)
899  real(kind=kreal) :: ctnforce(:)
900  real(kind=kreal) :: cttforce(:)
901  logical :: cflag
902  integer(kind=kint), optional, intent(in) :: smoothing_type
903 
904  real(kind=kreal) :: normal(3)
905  real(kind=kreal) :: bn(3*l_max_elem_node+3)
906  real(kind=kreal) :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3)
907  real(kind=kreal) :: elemcrd(3, l_max_elem_node)
908  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
909  real(kind=kreal) :: dgn, nrlforce
910  real(kind=kreal) :: lam_cone
911  real(kind=kreal) :: metric(2,2)
912  integer(kind=kint) :: edof
913 
914  nnode = size(tsurf%nodes)
915  edof = nnode*3+3
916 
917  ctnforce = 0.0d0
918  cttforce = 0.0d0
919 
920  normal(1:3) = ctstate%direction(1:3)
921 
922  ! Prepare elemcrd = ndCoord - ndDu (i.e., coord + disp) for computeContactMaps_ALag
923  do j = 1, nnode
924  elemcrd(1:3, j) = ndcoord(j*3+1:j*3+3) - nddu(j*3+1:j*3+3)
925  enddo
926 
927  ! Use common mapping routine to compute Bn, metric, Ht, Gt
928  call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
929  bn, metric, ht, gt, smoothing_type)
930 
931  ! Normal gap: dgn = Bn^T * ndCoord (using normal distribution vector)
932  dgn = dot_product( bn(1:edof), ndcoord(1:edof) )
933 
934  ! Normal force: multiplier + penalty * gap
935  nrlforce = ctstate%multiplier(1) + mu*dgn
936 
937  ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
938  ctnforce(1:edof) = -nrlforce * bn(1:edof)
939 
940  ! Lagrange row (not used in ALagrange, set to 0)
941  ctnforce((nnode+1)*3+1) = 0.d0
942 
943  if( fcoeff == 0.d0 ) return
944 
945  ! --- Tangent component ---
946 
947  ! Prepare edisp from ndDu
948  edisp(1:3) = nddu(1:3) ! slave
949  do j = 1, nnode
950  edisp(j*3+1:j*3+3) = nddu(j*3+1:j*3+3) ! master nodes
951  enddo
952 
953  ! Compute friction force using common routine. With FRICTION_CONE=FOLLOW the cone radius
954  ! is bounded by the normal force just distributed above rather than by the multiplier
955  ! alone, the same radius the tangent uses (see getContactStiffness_Alag).
956  if( symm ) then
957  lam_cone = ctstate%multiplier(1)
958  else
959  lam_cone = max( 0.d0, nrlforce )
960  endif
961  call computefrictionforce_alag(ctstate, fcoeff, lam_cone, metric, &
962  ht, gt, edisp, edof, cttforce, &
963  mut)
964 
965  ! Lagrange row (not used in ALagrange, set to 0)
966  cttforce((nnode+1)*3+1) = 0.d0
967 
968  end subroutine getcontactnodalforce_alag
969 
970  subroutine updatecontactmultiplier_alag(ctState,ndLocal,coord,disp,ddisp,&
971  & mu,mut,fcoeff,tSurf,lgnt,ctchanged,ctNForce,ctTForce,jump_ratio,smoothing_type)
972 
973  type(tcontactstate), intent(inout) :: ctstate
974  integer(kind=kint), intent(in) :: ndlocal(:)
975  real(kind=kreal), intent(in) :: coord(:)
976  real(kind=kreal), intent(in) :: disp(:)
977  real(kind=kreal), intent(in) :: ddisp(:)
978  real(kind=kreal), intent(in) :: mu, mut
979  real(kind=kreal), intent(in) :: fcoeff
980  type(tsurfelement), intent(in) :: tsurf
981  real(kind=kreal), intent(inout) :: lgnt(2)
982  logical, intent(inout) :: ctchanged
983  real(kind=kreal), intent(out) :: ctnforce(:)
984  real(kind=kreal), intent(out) :: cttforce(:)
985  real(kind=kreal), intent(out) :: jump_ratio
986  integer(kind=kint), optional, intent(in) :: smoothing_type
987 
988  integer(kind=kint) :: nnode
989  integer(kind=kint) :: slave, j
990  real(kind=kreal) :: bn(3*l_max_elem_node+3)
991  real(kind=kreal) :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3)
992  real(kind=kreal) :: elemcrd(3,l_max_elem_node)
993  real(kind=kreal) :: curpos(3*l_max_elem_node+3)
994  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
995  real(kind=kreal) :: dgn, nrlforce
996  real(kind=kreal) :: metric(2,2)
997  real(kind=kreal) :: dxy(2)
998  integer(kind=kint) :: edof
999 
1000  nnode = size(ndlocal) - 1
1001  slave = ndlocal(1)
1002  edof = nnode*3+3
1003 
1004  ctnforce = 0.0d0
1005  cttforce = 0.0d0
1006 
1007  ! Prepare elemcrd (coord+disp) and current positions (coord+disp+ddisp)
1008  curpos(1:3) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave) + ddisp(3*slave-2:3*slave)
1009  edisp(1:3) = ddisp(3*slave-2:3*slave)
1010  do j = 1, nnode
1011  elemcrd(1:3,j) = coord(3*ndlocal(j+1)-2:3*ndlocal(j+1)) + disp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1012  curpos(j*3+1:j*3+3) = elemcrd(1:3,j) + ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1013  edisp(j*3+1:j*3+3) = ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1014  enddo
1015 
1016  ! Use common mapping routine to compute Bn, metric, Ht, Gt
1017  call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
1018  bn, metric, ht, gt, smoothing_type)
1019 
1020  ! Normal gap: dgn = Bn^T * curpos (using normal distribution vector)
1021  dgn = dot_product( bn(1:edof), curpos(1:edof) )
1022 
1023  ! Update multiplier and working distance
1024  ctstate%wkdist = -dgn
1025  ctstate%multiplier(1) = ctstate%multiplier(1) - mu*ctstate%wkdist
1026  ctstate%distance = ctstate%wkdist
1027  lgnt(1) = lgnt(1) - ctstate%wkdist
1028 
1029  ! Normal force: use updated multiplier
1030  nrlforce = ctstate%multiplier(1)
1031 
1032  ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
1033  ctnforce(1:edof) = -nrlforce * bn(1:edof)
1034 
1035  if( fcoeff == 0.d0 ) return
1036 
1037  ! --- Tangent component ---
1038 
1039  ! The multiplier has just absorbed mu*g_n above, so it already is the normal force this
1040  ! configuration applies and serves as the cone radius for either FRICTION_CONE setting.
1041  call computefrictionforce_alag(ctstate, fcoeff, ctstate%multiplier(1), metric, &
1042  ht, gt, edisp, edof, cttforce, &
1043  mut, &
1044  update_multiplier=.true., slave_id=slave, ctchanged=ctchanged, &
1045  jump_ratio=jump_ratio)
1046 
1047  ! Tangent displacement for convergence check: use Gt to project curpos directly
1048  dxy = matmul( gt(:,1:edof), curpos(1:edof) )
1049  lgnt(2) = lgnt(2) + dsqrt( dxy(1)*dxy(1) + dxy(2)*dxy(2) )
1050 
1051  end subroutine updatecontactmultiplier_alag
1052 
1053  subroutine gettiedstiffness_alag(cstate, tSurf, mu, stiff, force)
1054 
1055  type(tcontactstate), intent(in) :: cstate
1056  type(tsurfelement), intent(in) :: tsurf
1057  real(kind=kreal), intent(in) :: mu
1058  real(kind=kreal), intent(out) :: stiff(:,:)
1059  real(kind=kreal), intent(out) :: force(:)
1060 
1061  integer :: i, j, nnode, edof
1062  real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
1063  real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
1064 
1065  nnode = size(tsurf%nodes)
1066  edof = nnode*3+3
1067 
1068  stiff = 0.d0
1069 
1070  ! Use common mapping routine to compute Tm
1071  call computetm_tt(cstate, tsurf, 0.0d0, tm, tt)
1072 
1073  ! Tied stiffness: stiff = mu * Tm^T * Tm
1074  do j = 1, edof
1075  do i = 1, edof
1076  stiff(i,j) = mu * dot_product(tm(1:3,i), tm(1:3,j))
1077  enddo
1078  enddo
1079  force(1:edof) = 0.d0 ! not used for tied (3-direction constraint)
1080 
1081  end subroutine gettiedstiffness_alag
1082 
1083  subroutine gettiednodalforce_alag(ctState,tSurf,ndu,mu,ctNForce,ctTForce)
1084 
1085  use msurfelement
1086  type(tcontactstate) :: ctstate
1087  type(tsurfelement) :: tsurf
1088  integer(kind=kint) :: nnode
1089  real(kind=kreal) :: ndu(:)
1090  real(kind=kreal), intent(in) :: mu
1091  real(kind=kreal) :: ctnforce(:)
1092  real(kind=kreal) :: cttforce(:)
1093 
1094  integer(kind=kint) :: edof
1095  real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
1096  real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
1097  real(kind=kreal) :: dg(3)
1098  real(kind=kreal) :: nrlforce(3)
1099 
1100  nnode = size(tsurf%nodes)
1101  edof = nnode*3+3
1102 
1103  ctnforce = 0.0d0
1104  cttforce = 0.0d0
1105 
1106  ! Use common mapping routine to compute Tm
1107  call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt)
1108 
1109  ! Gap vector: dg = Tm * ndu (3-component relative displacement)
1110  dg(1:3) = matmul(tm(1:3, 1:edof), ndu(1:edof))
1111 
1112  ! Force: multiplier + penalty * gap (3 components)
1113  nrlforce(1:3) = ctstate%multiplier(1:3) + mu*dg(1:3)
1114 
1115  ! Distribute force: ctNForce = -Tm^T * nrlforce
1116  ctnforce(1:edof) = -matmul(transpose(tm(1:3, 1:edof)), nrlforce(1:3))
1117 
1118  ! Lagrange row (not used in ALagrange, set to 0)
1119  ctnforce((nnode+1)*3+1) = 0.d0
1120  cttforce((nnode+1)*3+1) = 0.d0
1121 
1122  end subroutine gettiednodalforce_alag
1123 
1124 end module m_fstr_contact_elem_alag
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:754
elementinfo::get_intp_weights
subroutine get_intp_weights(etype, nn, n_intp, elecoord, weight)
Definition: element.f90:999
elementinfo::getintpoint4ss
subroutine getintpoint4ss(fetype, np, pos, n_intp, shapefunc)
Definition: element.f90:624
hecmw
Definition: hecmw.f90:6
m_fstr_contact_elem_alag
Alag method implementations for contact element calculations.
Definition: fstr_contact_elem_alag.f90:6
m_fstr_contact_elem_alag::updatecontactmultiplier_alag
subroutine, public updatecontactmultiplier_alag(ctState, ndLocal, coord, disp, ddisp, mu, mut, fcoeff, tSurf, lgnt, ctchanged, ctNForce, ctTForce, jump_ratio, smoothing_type)
Definition: fstr_contact_elem_alag.f90:972
m_fstr_contact_elem_alag::getcontactstiffness_alag_surfsurf
subroutine, public getcontactstiffness_alag_surfsurf(slave_surf, master, coord, disp, ddisp, mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, stiff_n, active_n, stiff_t, active_t)
Mortar (SURF-SURF) ALag: contact stiffness of one slave segment.
Definition: fstr_contact_elem_alag.f90:580
m_fstr_contact_elem_alag::getcontactstiffness_alag
subroutine, public getcontactstiffness_alag(cstate, tSurf, ele, mu, mut, fcoeff, symm, stiff, force, smoothing_type, edisp, iter, slvpos)
Definition: fstr_contact_elem_alag.f90:34
m_fstr_contact_elem_alag::gettiedstiffness_alag
subroutine, public gettiedstiffness_alag(cstate, tSurf, mu, stiff, force)
Definition: fstr_contact_elem_alag.f90:1054
m_fstr_contact_elem_alag::getcontactnodalforce_alag
subroutine, public getcontactnodalforce_alag(ctState, tSurf, ndCoord, ndDu, mu, mut, fcoeff, symm, lagrange, ctNForce, ctTForce, cflag, smoothing_type)
Definition: fstr_contact_elem_alag.f90:888
m_fstr_contact_elem_alag::gettiednodalforce_alag
subroutine, public gettiednodalforce_alag(ctState, tSurf, ndu, mu, ctNForce, ctTForce)
Definition: fstr_contact_elem_alag.f90:1084
m_fstr_contact_elem_alag::group_return_mapping
subroutine, public group_return_mapping(lam_t_in, rho_t, Dxi, fcoeff, lam_n, eps_fric_band, lam_t_out, fric_state, alpha, that, update_state)
Group-level Coulomb return mapping for the tangent multiplier (metric = I).
Definition: fstr_contact_elem_alag.f90:415
m_fstr_contact_elem_alag::get_unique_map
subroutine, public get_unique_map(sSurf, maplist, master_idxs, unique_count)
Definition: fstr_contact_elem_alag.f90:176
m_fstr_contact_elem_alag::build_group_tangent_basis
subroutine, public build_group_tangent_basis(n_hat, t1, t2)
Build an orthonormal tangent basis (t1,t2) as the orthogonal complement of a group-representative sla...
Definition: fstr_contact_elem_alag.f90:309
m_fstr_contact_elem_alag::gettangentslip
subroutine, public gettangentslip(slave_surf, master, coord, disp, ddisp, unique_count, maplist, master_idxs, Sigma_node, nacc_node)
Aggregate, per master group g and slave-surf node a of one mortar slave segment, the mortar-weighted ...
Definition: fstr_contact_elem_alag.f90:346
m_fstr_contact_elem_alag::getintgap
subroutine, public getintgap(slave_surf, master, coord, disp, ddisp, unique_count, maplist, master_idxs, Snode, Nsnode, gapwnode)
Compute the per-node mortar constraint quantities of one slave segment.
Definition: fstr_contact_elem_alag.f90:225
m_fstr_contact_elem_alag::getcontactnodalforce_alag_surfsurf
subroutine, public getcontactnodalforce_alag_surfsurf(purpose, slave_surf, master, coord, disp, ddisp, mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, ctNForce, active_n, ctTForce, active_t)
Mortar (SURF-SURF) ALag: contact nodal force of one slave segment.
Definition: fstr_contact_elem_alag.f90:757
m_fstr_contact_elem_alag::resolve_lambda_cur
subroutine, public resolve_lambda_cur(surf, master_idxs, unique_count, nnode_s, sorted_idx, lambda_node, lam_t_cur, fric_state_cur)
Mortar: resolve the current per-node lambda of each active group of one slave surf....
Definition: fstr_contact_elem_alag.f90:505
m_fstr_contact_elem_common
Common utilities for contact element calculations.
Definition: fstr_contact_elem_common.f90:6
m_fstr_contact_elem_common::computefrictionforce_alag
subroutine, public computefrictionforce_alag(ctState, fcoeff, lambda_n, metric, Ht, Gt, edisp, edof, ctTForce, mut, update_multiplier, slave_id, ctchanged, norm_trial, alpha, that, jump_ratio)
Compute friction force for ALag method Given tangent maps Ht, Gt and displacement increment,...
Definition: fstr_contact_elem_common.f90:198
m_fstr_contact_elem_common::computecontactmaps_alag
subroutine, public computecontactmaps_alag(ctState, tSurf, ele, Bn, metric, Ht, Gt, smoothing_type)
Compute contact maps for ALag method (normal distribution and tangential displacement maps) This subr...
Definition: fstr_contact_elem_common.f90:125
m_fstr_contact_elem_common::computetm_tt
subroutine, public computetm_tt(ctState, tSurf, fcoeff, Tm, Tt, smoothing_type, Bn)
Compute Tm (relative displacement mapping) and optionally Tt (tangential mapping) This subroutine con...
Definition: fstr_contact_elem_common.f90:27
m_fstr_contact_geom
This module provides geometric calculations for contact.
Definition: fstr_contact_geom.f90:6
m_fstr_contact_smoothing
Contact surface smoothing using Nagata patch interpolation.
Definition: fstr_contact_smoothing.f90:6
mcontactdef
This module manages the data structure for contact calculation.
Definition: fstr_contact_def.F90:12
mcontactdef::contactslip
integer, parameter contactslip
Definition: fstr_contact_def.F90:36
mcontactdef::kctforresidual
integer, parameter kctforresidual
purpose flag for contact force calculation
Definition: fstr_contact_def.F90:70
mcontactdef::max_n_intp
integer, parameter max_n_intp
upper bound of integration points per SURF-SURF slave segment
Definition: fstr_contact_def.F90:63
mcontactdef::kctforoutput
integer, parameter kctforoutput
compute contact force for output (CONT_NFORCE/CONT_FRIC)
Definition: fstr_contact_def.F90:71
mcontactdef::contactfree
integer, parameter contactfree
contact state definition
Definition: fstr_contact_def.F90:33
mcontactdef::contactstick
integer, parameter contactstick
Definition: fstr_contact_def.F90:35
msurfelement
This module manages surface elements in 3D It provides basic definition of surface elements (triangla...
Definition: surf_ele.f90:8
msurfelement::l_max_elem_node
integer(kind=kint), parameter l_max_elem_node
Definition: surf_ele.f90:17
msurfelement::l_max_surface_node
integer(kind=kint), parameter l_max_surface_node
Definition: surf_ele.f90:16
mcontactdef::tcontactstate
This structure records contact status.
Definition: fstr_contact_def.F90:74
mcontactdef::tcontactsurf
Structure to define a slave surface segment of a SURF-SURF (mortar) contact pair.
Definition: fstr_contact_def.F90:98
msurfelement::tsurfelement
Structure to define surface group.
Definition: surf_ele.f90:23