FrontISTR  6.0.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  maplist, 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) :: maplist(:)
350  real(kind=kreal), intent(out) :: sigma_node(:,:,:)
351  real(kind=kreal), intent(out) :: nacc_node(:,:,:)
352 
353  integer(kind=kint) :: i, j, g, a, nnode_s, nnode_m, etype, slave, n_intp, ctsurf, nd
354  real(kind=kreal) :: snode_pos(3,4), weight(max_n_intp)
355  real(kind=kreal) :: ncoord(2), shapefunc_s(4), shapefunc_m(4), direction(3)
356  real(kind=kreal) :: du_rel(3), du_master(3)
357 
358  nnode_s = size(slave_surf%nodes)
359  sigma_node = 0.d0
360  nacc_node = 0.d0
361 
362  ! Slave node positions at start of substep (coord + disp), used for IP weights
363  ! (identical to getIntGap so the weights match the normal aggregation exactly).
364  snode_pos = 0.d0
365  do i = 1, nnode_s
366  slave = slave_surf%nodes(i)
367  snode_pos(:,i) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave)
368  enddo
369  n_intp = slave_surf%n_intp
370  weight = 0.d0
371  call get_intp_weights(slave_surf%etype, nnode_s, n_intp, snode_pos, weight(1:n_intp))
372 
373  do i = 1, n_intp
374  if( slave_surf%states(i)%state == contactfree ) cycle
375  ctsurf = slave_surf%states(i)%surface
376  etype = master(ctsurf)%etype
377  nnode_m = size(master(ctsurf)%nodes)
378  direction = slave_surf%states(i)%direction(1:3)
379  call getintpoint4ss(slave_surf%etype, i, ncoord, n_intp, shapefunc_s)
380  call getshapefunc(etype, slave_surf%states(i)%lpos(1:2), shapefunc_m)
381  g = maplist(i)
382 
383  ! Mortar relative displacement = Tm_mortar . edisp (edisp = ddisp).
384  ! slave block: + N_s(j) * ddisp(slave_j) ; master block: - N_m(j) * ddisp(master_j).
385  du_rel(1:3) = 0.d0
386  do j = 1, nnode_s
387  nd = slave_surf%nodes(j)
388  du_rel(1:3) = du_rel(1:3) + shapefunc_s(j) * ddisp(3*nd-2:3*nd)
389  enddo
390  du_master(1:3) = 0.d0
391  do j = 1, nnode_m
392  nd = master(ctsurf)%nodes(j)
393  du_master(1:3) = du_master(1:3) + shapefunc_m(j) * ddisp(3*nd-2:3*nd)
394  enddo
395  du_rel(1:3) = du_rel(1:3) - du_master(1:3)
396 
397  ! Per-node: weight the IP slip/normal by the slave node shape function N_s(a).
398  do a = 1, nnode_s
399  sigma_node(g, a, 1:3) = sigma_node(g, a, 1:3) + shapefunc_s(a) * weight(i) * du_rel(1:3)
400  nacc_node(g, a, 1:3) = nacc_node(g, a, 1:3) + shapefunc_s(a) * weight(i) * direction(1:3)
401  enddo
402  enddo
403  end subroutine gettangentslip
404 
412  subroutine group_return_mapping(lam_t_in, rho_t, Dxi, fcoeff, lam_n, eps_fric_band, &
413  lam_t_out, fric_state, alpha, that, update_state)
414  real(kind=kreal), intent(in) :: lam_t_in(2)
415  real(kind=kreal), intent(in) :: rho_t
416  real(kind=kreal), intent(in) :: dxi(2)
417  real(kind=kreal), intent(in) :: fcoeff
418  real(kind=kreal), intent(in) :: lam_n
419  real(kind=kreal), intent(in) :: eps_fric_band
420  real(kind=kreal), intent(out) :: lam_t_out(2)
421  integer(kind=kint), intent(inout) :: fric_state
422  real(kind=kreal), intent(out) :: alpha
423  real(kind=kreal), intent(out) :: that(2)
424  logical, optional, intent(in) :: update_state
425 
426  real(kind=kreal) :: trial(2), norm_trial, radius
427  logical :: is_stick, do_update
428 
429  do_update = .false.
430  if( present(update_state) ) do_update = update_state
431 
432  trial(1:2) = lam_t_in(1:2) + rho_t * dxi(1:2)
433  norm_trial = sqrt( trial(1)*trial(1) + trial(2)*trial(2) ) ! metric = I
434  radius = fcoeff * lam_n
435 
436  if( norm_trial > 1.0d-20 ) then
437  that(1:2) = trial(1:2) / norm_trial
438  else
439  that(1:2) = 0.0d0
440  endif
441 
442  if( lam_n <= 0.0d0 ) then
443  ! No normal force: cone collapses, no friction. Leave fric_state as-is.
444  lam_t_out(1:2) = 0.0d0
445  alpha = 0.0d0
446  else if( do_update ) then
447  ! Augmentation: this is the only place where the stick/slip state is (re)decided.
448  ! Hysteresis band: read the warm-start (previous) state and switch only when
449  ! ||trial|| crosses the asymmetric thresholds. Inside the band the previous
450  ! branch is kept, so marginal flips at the cone boundary cannot drive a limit
451  ! cycle across augmentations.
452  is_stick = ( fric_state == contactstick )
453  if( is_stick ) then
454  if( norm_trial > (1.0d0 + eps_fric_band) * radius ) is_stick = .false. ! genuine slip onset
455  else
456  if( norm_trial <= (1.0d0 - eps_fric_band) * radius ) is_stick = .true. ! genuine stick recovery
457  endif
458 
459  if( is_stick ) then
460  ! Stick branch: keep full trial multiplier.
461  lam_t_out(1:2) = trial(1:2)
462  fric_state = contactstick
463  alpha = 1.0d0
464  else
465  ! Slip branch: project onto cone surface.
466  alpha = radius / norm_trial
467  lam_t_out(1:2) = alpha * trial(1:2)
468  fric_state = contactslip
469  endif
470  else
471  ! Inner Newton-Raphson (read-only): pin the state frozen at the previous
472  ! augmentation, skipping the band re-judgement. fric_state is not written.
473  ! Same isolation as NTS-ALAG (update_multiplier=.false.): the stick/slip
474  ! branch cannot flip on the live trial inside the inner NR, so the residual
475  ! stays smooth and the per-iteration branch-flip limit cycle is broken.
476  if( fric_state == contactstick ) then
477  ! Frozen STICK: keep full trial multiplier (no cone projection).
478  lam_t_out(1:2) = trial(1:2)
479  alpha = 1.0d0
480  else
481  ! Frozen SLIP: project onto the live cone surface with min(1, radius/||trial||), the same
482  ! expression computeFrictionForce_ALag uses on the node-to-surface side. The projection
483  ! never scales a trial up: a node whose trial has come back inside the cone keeps the full
484  ! trial force, which is what the alpha >= 0.999 branch of the tangent is linearised about.
485  ! A node whose state was frozen SLIP at an augmentation where it carried no normal force
486  ! (the lam_n <= 0 branch leaves fric_state untouched) can come back with a zero trial, and
487  ! radius/0 would turn the whole residual into NaN. No trial force means no friction force,
488  ! which is what alpha = 0 gives.
489  if( norm_trial > 1.0d-20 ) then
490  alpha = min( 1.0d0, radius / norm_trial )
491  else
492  alpha = 0.0d0
493  endif
494  lam_t_out(1:2) = alpha * trial(1:2)
495  endif
496  endif
497  end subroutine group_return_mapping
498 
502  subroutine resolve_lambda_cur( surf, master_idxs, unique_count, nnode_s, sorted_idx, &
503  lambda_node, lam_t_cur, fric_state_cur )
504  type(tcontactsurf), intent(in) :: surf
505  integer(kind=kint), intent(in) :: master_idxs(:)
506  integer(kind=kint), intent(in) :: unique_count
507  integer(kind=kint), intent(in) :: nnode_s
508  integer(kind=kint), intent(out) :: sorted_idx(:)
509  real(kind=kreal), intent(out) :: lambda_node(:,:)
510  ! Optional friction warm-start: same reference rule as lambda_n (working -> begin ->
511  ! default), riding the same merge.
512  real(kind=kreal), intent(out), optional :: lam_t_cur(:,:,:)
513  integer(kind=kint), intent(out), optional :: fric_state_cur(:,:)
514  integer(kind=kint) :: r, j, tmp, ib, iw, g, mid
515  logical :: do_fric, hit_work, hit_begin
516 
517  ! argsort master_idxs ascending (unique_count <= 27, insertion sort)
518  do r = 1, unique_count
519  sorted_idx(r) = r
520  enddo
521  do r = 2, unique_count
522  tmp = sorted_idx(r)
523  j = r - 1
524  do while( j >= 1 )
525  if( master_idxs(sorted_idx(j)) <= master_idxs(tmp) ) exit
526  sorted_idx(j+1) = sorted_idx(j)
527  j = j - 1
528  enddo
529  sorted_idx(j+1) = tmp
530  enddo
531 
532  do_fric = present(lam_t_cur) .and. present(fric_state_cur)
533 
534  ! 2-pointer merge over ascending masters / ascending begin / ascending working
535  ib = 1; iw = 1
536  do r = 1, unique_count
537  mid = master_idxs(sorted_idx(r))
538  do while( ib <= surf%lam_begin_n )
539  if( surf%lam_begin_id(ib) >= mid ) exit
540  ib = ib + 1
541  enddo
542  do while( iw <= surf%lam_work_n )
543  if( surf%lam_work_id(iw) >= mid ) exit
544  iw = iw + 1
545  enddo
546  g = sorted_idx(r)
547  hit_work = .false.
548  if( iw <= surf%lam_work_n ) hit_work = surf%lam_work_id(iw) == mid
549  hit_begin = .false.
550  if( ib <= surf%lam_begin_n ) hit_begin = surf%lam_begin_id(ib) == mid
551  if( hit_work ) then
552  lambda_node(1:nnode_s,g) = surf%lam_work_val(1:nnode_s,iw)
553  if( do_fric ) then
554  lam_t_cur(1:2,1:nnode_s,g) = surf%lam_work_t(1:2,1:nnode_s,iw)
555  fric_state_cur(1:nnode_s,g) = surf%lam_work_fstate(1:nnode_s,iw)
556  endif
557  else if( hit_begin ) then
558  lambda_node(1:nnode_s,g) = surf%lam_begin_val(1:nnode_s,ib)
559  if( do_fric ) then
560  lam_t_cur(1:2,1:nnode_s,g) = surf%lam_begin_t(1:2,1:nnode_s,ib)
561  fric_state_cur(1:nnode_s,g) = surf%lam_begin_fstate(1:nnode_s,ib)
562  endif
563  else
564  lambda_node(1:nnode_s,g) = 0.d0
565  if( do_fric ) then
566  lam_t_cur(1:2,1:nnode_s,g) = 0.d0
567  fric_state_cur(1:nnode_s,g) = contactstick
568  endif
569  endif
570  enddo
571  end subroutine resolve_lambda_cur
572 
586  subroutine getcontactstiffness_alag_surfsurf( slave_surf, master, coord, disp, ddisp, &
587  mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, &
588  stiff_n, active_n, stiff_t, active_t )
589  type(tcontactsurf), intent(in) :: slave_surf
590  type(tsurfelement), intent(in) :: master(:)
591  real(kind=kreal), intent(in) :: coord(:)
592  real(kind=kreal), intent(in) :: disp(:)
593  real(kind=kreal), intent(in) :: ddisp(:)
594  real(kind=kreal), intent(in) :: mu, mut
595  real(kind=kreal), intent(in) :: fcoeff
596  logical, intent(in) :: symm
597  real(kind=kreal), intent(in) :: eps_fric_band
598  integer(kind=kint), intent(in) :: unique_count
599  integer(kind=kint), intent(in) :: maplist(:)
600  integer(kind=kint), intent(in) :: master_idxs(:)
601  real(kind=kreal), intent(out) :: stiff_n(:,:,:,:)
602  logical, intent(out) :: active_n(:,:)
603  real(kind=kreal), intent(out), optional :: stiff_t(:,:,:,:)
604  logical, intent(out), optional :: active_t(:,:)
605 
606  integer(kind=kint) :: g, a, j, k, nnode_m, nnode_s
607  integer(kind=kint), allocatable :: sorted_idx(:)
608  real(kind=kreal), allocatable :: snode(:,:), nsnode(:,:,:), gapwnode(:,:), lambda_node(:,:)
609  real(kind=kreal) :: ns(24)
610  ! --- friction consistent tangent ---
611  real(kind=kreal), allocatable :: lam_t_cur(:,:,:), nacc_node(:,:,:), sigma_node(:,:,:)
612  integer(kind=kint), allocatable :: fric_state_cur(:,:)
613  real(kind=kreal) :: nhat(3), t1(3), t2(3), nrm, dxi(2)
614  real(kind=kreal) :: alpha, that(2), lam_t_new(2), amat(2,2), t3d(3,2), m3(3,3)
615  real(kind=kreal) :: wb(l_max_surface_node+1)
616  real(kind=kreal) :: lam_cone, that3d(3)
617  integer(kind=kint) :: na, nb, fstate
618 
619  nnode_s = size(slave_surf%nodes)
620  allocate(snode(unique_count,nnode_s), nsnode(unique_count,nnode_s,24), gapwnode(unique_count,nnode_s))
621  call getintgap(slave_surf, master, coord, disp, ddisp, &
622  unique_count, maplist, master_idxs, &
623  snode, nsnode, gapwnode)
624 
625  allocate(sorted_idx(unique_count), lambda_node(nnode_s,unique_count))
626  stiff_n = 0.d0
627  active_n = .false.
628  if( fcoeff /= 0.d0 ) then
629  allocate(lam_t_cur(2,nnode_s,unique_count), fric_state_cur(nnode_s,unique_count))
630  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
631  lambda_node, lam_t_cur, fric_state_cur)
632  else
633  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
634  lambda_node)
635  endif
636 
637  ! ===== Normal stiffness =====
638  do g = 1, unique_count
639  nnode_m = size(master(master_idxs(g))%nodes)
640  do a = 1, nnode_s
641  ! ALag contact condition per node: augmented force must be positive
642  if( lambda_node(a,g)+mu*gapwnode(g,a) < 0.d0 ) cycle
643  active_n(a,g) = .true.
644  ns = 0.d0
645  ns(1:(nnode_s+nnode_m)*3) = nsnode(g, a, 1:(nnode_s+nnode_m)*3)
646  do j = 1, (nnode_s+nnode_m)*3
647  do k = 1, (nnode_s+nnode_m)*3
648  stiff_n(j,k,a,g) = mu*snode(g,a)*ns(j)*ns(k)
649  enddo
650  enddo
651  enddo
652  enddo
653 
654  if( fcoeff /= 0.d0 ) then
655  ! ===== Friction consistent tangent (per slave node) =====
656  ! Linearization of the per-node friction residual at the same live slip state:
657  ! K_a(b,c) = Snode(g,a) * Wbar(a,b) * Wbar(a,c) * M3_a, M3_a = T3d_a * A_a * T3d_a^T
658  ! with Wbar(a,b) = Nsnode(g,a,b).nhat_a, the same map as the residual back-distribution.
659  stiff_t = 0.d0
660  active_t = .false.
661  allocate(sigma_node(unique_count,nnode_s,3), nacc_node(unique_count,nnode_s,3))
662  call gettangentslip(slave_surf, master, coord, disp, ddisp, &
663  maplist, sigma_node, nacc_node)
664  do g = 1, unique_count
665  nnode_m = size(master(master_idxs(g))%nodes)
666  do a = 1, nnode_s
667  if( lambda_node(a,g) <= 0.d0 ) cycle ! no per-node normal force -> no friction
668  ! Radius of the friction cone. With FRICTION_CONE=FROZEN it stays at the multiplier
669  ! of the last augmentation, which keeps the friction terms symmetric and leaves the
670  ! Coulomb condition to the augmentation loop; with !CONTACT_ALGO, FRICTION_CONE=FOLLOW it
671  ! follows the normal force this node actually applies, lambda_node + rho_n*gapwnode,
672  ! the same expression the residual distributes as nrlforce.
673  if( symm ) then
674  lam_cone = lambda_node(a,g)
675  else
676  lam_cone = lambda_node(a,g) + mu*gapwnode(g,a)
677  endif
678  nrm = sqrt( nacc_node(g,a,1)**2 + nacc_node(g,a,2)**2 + nacc_node(g,a,3)**2 )
679  if( nrm < 1.d-30 ) cycle
680  nhat(1:3) = nacc_node(g,a,1:3) / nrm
681  call build_group_tangent_basis(nhat, t1, t2)
682  ! Live per-node slip projection and read-only return mapping (writes only OUT args).
683  dxi(1) = dot_product(t1(1:3), sigma_node(g,a,1:3))
684  dxi(2) = dot_product(t2(1:3), sigma_node(g,a,1:3))
685  fstate = fric_state_cur(a,g)
686  call group_return_mapping(lam_t_cur(1:2,a,g), mut, dxi, fcoeff, lam_cone, &
687  eps_fric_band, lam_t_new, fstate, alpha, that)
688  ! 2D tangent operator A (same construction as getContactStiffness_Alag).
689  if( alpha <= 1.0d-20 ) then
690  amat = 0.d0
691  else if( alpha >= 0.999d0 ) then
692  amat = 0.d0
693  amat(1,1) = mut
694  amat(2,2) = mut
695  else
696  amat(1,1) = alpha * mut * (1.0d0 - that(1)*that(1))
697  amat(1,2) = alpha * mut * (-that(1)*that(2))
698  amat(2,1) = alpha * mut * (-that(2)*that(1))
699  amat(2,2) = alpha * mut * (1.0d0 - that(2)*that(2))
700  endif
701  t3d(1:3,1) = t1(1:3)
702  t3d(1:3,2) = t2(1:3)
703  m3 = matmul( matmul(t3d, amat), transpose(t3d) )
704 
705  ! Per-node averaged mortar weight of each node (= ANnode/Snode, recovered via nhat_a).
706  do na = 1, nnode_s + nnode_m
707  wb(na) = dot_product(nsnode(g,a,3*na-2:3*na), nhat(1:3))
708  enddo
709  active_t(a,g) = .true.
710  do nb = 1, nnode_s + nnode_m
711  do na = 1, nnode_s + nnode_m
712  do k = 1, 3
713  do j = 1, 3
714  stiff_t(3*na-3+j, 3*nb-3+k, a, g) = snode(g,a) * wb(na) * wb(nb) * m3(j,k)
715  enddo
716  enddo
717  enddo
718  enddo
719  ! Coupling block of a cone radius that follows the normal force. On the slip branch
720  ! the friction force is f_t = R*that3d with R = fcoeff*lam_cone, and R varies with u
721  ! through gapwnode: d(gapwnode(g,a))/du = Snode(g,a)*Nsnode(g,a,:), the map the normal
722  ! stiffness uses, so the residual -Wbar(a,b)*f_t gains
723  ! K_a(b,c) += fcoeff*rho_n*Snode(g,a) * Wbar(a,b)*that3d (x) Nsnode(g,a,c).
724  ! Rows are a slip direction and columns a normal map, so the block is unsymmetric and
725  ! the solver is set up for a general matrix (fstr_is_contactALag_symmetric). A stuck
726  ! node does not use the radius (f_t is the full trial), hence the slip-branch window,
727  ! the same one the consistent tangent above uses.
728  if( .not.symm .and. alpha > 1.0d-20 .and. alpha < 0.999d0 ) then
729  that3d(1:3) = that(1)*t1(1:3) + that(2)*t2(1:3)
730  do nb = 1, nnode_s + nnode_m
731  do na = 1, nnode_s + nnode_m
732  do k = 1, 3
733  do j = 1, 3
734  stiff_t(3*na-3+j, 3*nb-3+k, a, g) = stiff_t(3*na-3+j, 3*nb-3+k, a, g) &
735  + fcoeff * mu * snode(g,a) * wb(na) * that3d(j) * nsnode(g,a,3*nb-3+k)
736  enddo
737  enddo
738  enddo
739  enddo
740  endif
741  enddo
742  enddo
743 
744  deallocate(sigma_node, nacc_node, lam_t_cur, fric_state_cur)
745  endif
746 
747  deallocate(sorted_idx)
748  deallocate(snode, nsnode, gapwnode, lambda_node)
749  end subroutine getcontactstiffness_alag_surfsurf
750 
763  subroutine getcontactnodalforce_alag_surfsurf( purpose, slave_surf, master, coord, disp, ddisp, &
764  mu, mut, fcoeff, symm, eps_fric_band, unique_count, maplist, master_idxs, &
765  ctNForce, active_n, ctTForce, active_t )
766  integer(kind=kint), intent(in) :: purpose
767  type(tcontactsurf), intent(in) :: slave_surf
768  type(tsurfelement), intent(in) :: master(:)
769  real(kind=kreal), intent(in) :: coord(:)
770  real(kind=kreal), intent(in) :: disp(:)
771  real(kind=kreal), intent(in) :: ddisp(:)
772  real(kind=kreal), intent(in) :: mu, mut
773  real(kind=kreal), intent(in) :: fcoeff
774  logical, intent(in) :: symm
775  real(kind=kreal), intent(in) :: eps_fric_band
776  integer(kind=kint), intent(in) :: unique_count
777  integer(kind=kint), intent(in) :: maplist(:)
778  integer(kind=kint), intent(in) :: master_idxs(:)
779  real(kind=kreal), intent(out) :: ctnforce(:,:,:)
780  logical, intent(out) :: active_n(:,:)
781  real(kind=kreal), intent(out), optional :: cttforce(:,:,:)
782  logical, intent(out), optional :: active_t(:,:)
783 
784  integer(kind=kint) :: g, a, j, nnode_m, nnode_s
785  integer(kind=kint), allocatable :: sorted_idx(:)
786  real(kind=kreal), allocatable :: snode(:,:), nsnode(:,:,:), gapwnode(:,:), lambda_node(:,:)
787  real(kind=kreal) :: nrlforce
788  real(kind=kreal) :: ns(24)
789  ! --- friction force back-distribution (live return mapping) ---
790  real(kind=kreal), allocatable :: lam_t_cur(:,:,:), nacc_node(:,:,:), sigma_node(:,:,:)
791  integer(kind=kint), allocatable :: fric_state_cur(:,:)
792  real(kind=kreal) :: nhat(3), t1(3), t2(3), nrm, fvec(3), wbar
793  real(kind=kreal) :: dxi(2), alpha, that(2), lam_t_new(2), lam_cone
794  integer(kind=kint) :: fstate
795 
796  nnode_s = size(slave_surf%nodes)
797  allocate(snode(unique_count,nnode_s), nsnode(unique_count,nnode_s,24), gapwnode(unique_count,nnode_s))
798  call getintgap(slave_surf, master, coord, disp, ddisp, &
799  unique_count, maplist, master_idxs, &
800  snode, nsnode, gapwnode)
801 
802  allocate(sorted_idx(unique_count), lambda_node(nnode_s,unique_count))
803  ctnforce = 0.d0
804  active_n = .false.
805  if( fcoeff /= 0.d0 ) then
806  allocate(lam_t_cur(2,nnode_s,unique_count), fric_state_cur(nnode_s,unique_count))
807  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
808  lambda_node, lam_t_cur, fric_state_cur)
809  else
810  call resolve_lambda_cur(slave_surf, master_idxs, unique_count, nnode_s, sorted_idx, &
811  lambda_node)
812  endif
813 
814  ! ===== Normal force: per-node back-distribution =====
815  ! nrlforce_a = lambda_node(a,g) + mu*gapwnode(g,a) for the residual, lambda_node(a,g) for output.
816  do g = 1, unique_count
817  nnode_m = size(master(master_idxs(g))%nodes)
818  do a = 1, nnode_s
819  nrlforce = lambda_node(a,g) + mu*gapwnode(g,a)
820  ! ALag contact condition per node: augmented force must be positive
821  if( nrlforce < 0.d0 ) cycle
822  active_n(a,g) = .true.
823  ns = 0.d0
824  ns(1:(nnode_s+nnode_m)*3) = nsnode(g, a, 1:(nnode_s+nnode_m)*3)
825  do j = 1, nnode_s + nnode_m
826  if( purpose == kctforresidual ) then
827  ctnforce(3*j-2:3*j,a,g) = -nrlforce*ns(3*j-2:3*j)
828  else if ( purpose == kctforoutput ) then
829  ! Output: multiplier only (converges to true contact force)
830  ctnforce(3*j-2:3*j,a,g) = -lambda_node(a,g)*ns(3*j-2:3*j)
831  end if
832  enddo
833  enddo
834  enddo
835 
836  if( fcoeff /= 0.d0 ) then
837  ! ===== Friction force (per slave node): live return mapping, back-distributed =====
838  ! trial = lam_t_warm(a) + rho_t*Dxi_live(a), projected onto the cone of radius lambda_node(a,g).
839  ! The return mapping is read-only here (the augmentation update is the sole writer of the
840  ! lambda_t / fric_state buffers). The resulting traction is distributed through the per-node
841  ! mortar weight Wbar(a,j) = Nsnode(g,a,j).nhat_a, mirroring the normal back-distribution.
842  ! Output (kctForOutput) keeps the frozen multiplier.
843  cttforce = 0.d0
844  active_t = .false.
845  allocate(sigma_node(unique_count,nnode_s,3), nacc_node(unique_count,nnode_s,3))
846  call gettangentslip(slave_surf, master, coord, disp, ddisp, &
847  maplist, sigma_node, nacc_node)
848  do g = 1, unique_count
849  nnode_m = size(master(master_idxs(g))%nodes)
850  do a = 1, nnode_s
851  if( lambda_node(a,g) <= 0.d0 ) cycle ! no per-node normal force -> no friction
852  ! Cone radius: the frozen multiplier with FRICTION_CONE=FROZEN, the normal force this
853  ! node just applied above (nrlforce = lambda_node + rho_n*gapwnode) with FRICTION_CONE=FOLLOW,
854  ! the same radius the tangent uses (see getContactStiffness_Alag_SurfSurf). A negative
855  ! lam_cone reaches group_return_mapping as lam_n <= 0 and gives zero friction, which is
856  ! what the normal back-distribution above does with a negative nrlforce too.
857  if( symm ) then
858  lam_cone = lambda_node(a,g)
859  else
860  lam_cone = lambda_node(a,g) + mu*gapwnode(g,a)
861  endif
862  nrm = sqrt( nacc_node(g,a,1)**2 + nacc_node(g,a,2)**2 + nacc_node(g,a,3)**2 )
863  if( nrm < 1.d-30 ) cycle
864  nhat(1:3) = nacc_node(g,a,1:3) / nrm
865  call build_group_tangent_basis(nhat, t1, t2)
866  if( purpose == kctforresidual ) then
867  ! Live trial: project the live per-node mortar slip onto the per-node frame and return-map.
868  dxi(1) = dot_product(t1(1:3), sigma_node(g,a,1:3))
869  dxi(2) = dot_product(t2(1:3), sigma_node(g,a,1:3))
870  fstate = fric_state_cur(a,g)
871  call group_return_mapping(lam_t_cur(1:2,a,g), mut, dxi, fcoeff, lam_cone, &
872  eps_fric_band, lam_t_new, fstate, alpha, that)
873  fvec(1:3) = lam_t_new(1)*t1(1:3) + lam_t_new(2)*t2(1:3)
874  else
875  ! Output: frozen multiplier only (converges to the true friction force).
876  fvec(1:3) = lam_t_cur(1,a,g)*t1(1:3) + lam_t_cur(2,a,g)*t2(1:3)
877  end if
878 
879  active_t(a,g) = .true.
880  do j = 1, nnode_s + nnode_m
881  ! per-node averaged mortar weight of node j (= ANnode/Snode, recovered via nhat_a)
882  wbar = dot_product(nsnode(g,a,3*j-2:3*j), nhat(1:3))
883  cttforce(3*j-2:3*j,a,g) = -fvec(1:3) * wbar
884  enddo
885  enddo
886  enddo
887 
888  deallocate(sigma_node, nacc_node, lam_t_cur, fric_state_cur)
889  endif
890 
891  deallocate(sorted_idx)
892  deallocate(snode, nsnode, gapwnode, lambda_node)
893  end subroutine getcontactnodalforce_alag_surfsurf
894 
895  subroutine getcontactnodalforce_alag(ctState,tSurf,ndCoord,ndDu,mu,mut,fcoeff,symm,ctNForce,ctTForce, &
896  smoothing_type)
897 
898  use msurfelement
899  type(tcontactstate) :: ctstate
900  type(tsurfelement) :: tsurf
901  integer(kind=kint) :: nnode
902  integer(kind=kint) :: j
903  real(kind=kreal), intent(in) :: mu, mut
904  real(kind=kreal) :: fcoeff
905  logical, intent(in) :: symm
906  real(kind=kreal) :: ndcoord(:), nddu(:)
907  real(kind=kreal) :: ctnforce(:)
908  real(kind=kreal) :: cttforce(:)
909  integer(kind=kint), optional, intent(in) :: smoothing_type
910 
911  real(kind=kreal) :: normal(3)
912  real(kind=kreal) :: bn(3*l_max_elem_node+3)
913  real(kind=kreal) :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3)
914  real(kind=kreal) :: elemcrd(3, l_max_elem_node)
915  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
916  real(kind=kreal) :: dgn, nrlforce
917  real(kind=kreal) :: lam_cone
918  real(kind=kreal) :: metric(2,2)
919  integer(kind=kint) :: edof
920 
921  nnode = size(tsurf%nodes)
922  edof = nnode*3+3
923 
924  ctnforce = 0.0d0
925  cttforce = 0.0d0
926 
927  normal(1:3) = ctstate%direction(1:3)
928 
929  ! Prepare elemcrd = ndCoord - ndDu (i.e., coord + disp) for computeContactMaps_ALag
930  do j = 1, nnode
931  elemcrd(1:3, j) = ndcoord(j*3+1:j*3+3) - nddu(j*3+1:j*3+3)
932  enddo
933 
934  ! Use common mapping routine to compute Bn, metric, Ht, Gt
935  call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
936  bn, metric, ht, gt, smoothing_type)
937 
938  ! Normal gap: dgn = Bn^T * ndCoord (using normal distribution vector)
939  dgn = dot_product( bn(1:edof), ndcoord(1:edof) )
940 
941  ! Normal force: multiplier + penalty * gap
942  nrlforce = ctstate%multiplier(1) + mu*dgn
943 
944  ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
945  ctnforce(1:edof) = -nrlforce * bn(1:edof)
946 
947  ! Lagrange row (not used in ALagrange, set to 0)
948  ctnforce((nnode+1)*3+1) = 0.d0
949 
950  if( fcoeff == 0.d0 ) return
951 
952  ! --- Tangent component ---
953 
954  ! Prepare edisp from ndDu
955  edisp(1:3) = nddu(1:3) ! slave
956  do j = 1, nnode
957  edisp(j*3+1:j*3+3) = nddu(j*3+1:j*3+3) ! master nodes
958  enddo
959 
960  ! Compute friction force using common routine. With FRICTION_CONE=FOLLOW the cone radius
961  ! is bounded by the normal force just distributed above rather than by the multiplier
962  ! alone, the same radius the tangent uses (see getContactStiffness_Alag).
963  if( symm ) then
964  lam_cone = ctstate%multiplier(1)
965  else
966  lam_cone = max( 0.d0, nrlforce )
967  endif
968  call computefrictionforce_alag(ctstate, fcoeff, lam_cone, metric, &
969  ht, gt, edisp, edof, cttforce, &
970  mut)
971 
972  ! Lagrange row (not used in ALagrange, set to 0)
973  cttforce((nnode+1)*3+1) = 0.d0
974 
975  end subroutine getcontactnodalforce_alag
976 
977  subroutine updatecontactmultiplier_alag(ctState,ndLocal,coord,disp,ddisp,&
978  & mu,mut,fcoeff,tSurf,lgnt,ctchanged,ctNForce,ctTForce,jump_ratio,smoothing_type)
979 
980  type(tcontactstate), intent(inout) :: ctstate
981  integer(kind=kint), intent(in) :: ndlocal(:)
982  real(kind=kreal), intent(in) :: coord(:)
983  real(kind=kreal), intent(in) :: disp(:)
984  real(kind=kreal), intent(in) :: ddisp(:)
985  real(kind=kreal), intent(in) :: mu, mut
986  real(kind=kreal), intent(in) :: fcoeff
987  type(tsurfelement), intent(in) :: tsurf
988  real(kind=kreal), intent(inout) :: lgnt(2)
989  logical, intent(inout) :: ctchanged
990  real(kind=kreal), intent(out) :: ctnforce(:)
991  real(kind=kreal), intent(out) :: cttforce(:)
992  real(kind=kreal), intent(out) :: jump_ratio
993  integer(kind=kint), optional, intent(in) :: smoothing_type
994 
995  integer(kind=kint) :: nnode
996  integer(kind=kint) :: slave, j
997  real(kind=kreal) :: bn(3*l_max_elem_node+3)
998  real(kind=kreal) :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3)
999  real(kind=kreal) :: elemcrd(3,l_max_elem_node)
1000  real(kind=kreal) :: curpos(3*l_max_elem_node+3)
1001  real(kind=kreal) :: edisp(3*l_max_elem_node+3)
1002  real(kind=kreal) :: dgn, nrlforce
1003  real(kind=kreal) :: metric(2,2)
1004  real(kind=kreal) :: dxy(2)
1005  integer(kind=kint) :: edof
1006 
1007  nnode = size(ndlocal) - 1
1008  slave = ndlocal(1)
1009  edof = nnode*3+3
1010 
1011  ctnforce = 0.0d0
1012  cttforce = 0.0d0
1013  jump_ratio = 0.0d0
1014 
1015  ! Prepare elemcrd (coord+disp) and current positions (coord+disp+ddisp)
1016  curpos(1:3) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave) + ddisp(3*slave-2:3*slave)
1017  edisp(1:3) = ddisp(3*slave-2:3*slave)
1018  do j = 1, nnode
1019  elemcrd(1:3,j) = coord(3*ndlocal(j+1)-2:3*ndlocal(j+1)) + disp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1020  curpos(j*3+1:j*3+3) = elemcrd(1:3,j) + ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1021  edisp(j*3+1:j*3+3) = ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
1022  enddo
1023 
1024  ! Use common mapping routine to compute Bn, metric, Ht, Gt
1025  call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
1026  bn, metric, ht, gt, smoothing_type)
1027 
1028  ! Normal gap: dgn = Bn^T * curpos (using normal distribution vector)
1029  dgn = dot_product( bn(1:edof), curpos(1:edof) )
1030 
1031  ! Update multiplier and working distance
1032  ctstate%wkdist = -dgn
1033  ctstate%multiplier(1) = ctstate%multiplier(1) - mu*ctstate%wkdist
1034  ctstate%distance = ctstate%wkdist
1035  lgnt(1) = lgnt(1) - ctstate%wkdist
1036 
1037  ! Normal force: use updated multiplier
1038  nrlforce = ctstate%multiplier(1)
1039 
1040  ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
1041  ctnforce(1:edof) = -nrlforce * bn(1:edof)
1042 
1043  if( fcoeff == 0.d0 ) return
1044 
1045  ! --- Tangent component ---
1046 
1047  ! The multiplier has just absorbed mu*g_n above, so it already is the normal force this
1048  ! configuration applies and serves as the cone radius for either FRICTION_CONE setting.
1049  call computefrictionforce_alag(ctstate, fcoeff, ctstate%multiplier(1), metric, &
1050  ht, gt, edisp, edof, cttforce, &
1051  mut, &
1052  update_multiplier=.true., slave_id=slave, ctchanged=ctchanged, &
1053  jump_ratio=jump_ratio)
1054 
1055  ! Tangent displacement for convergence check: use Gt to project curpos directly
1056  dxy = matmul( gt(:,1:edof), curpos(1:edof) )
1057  lgnt(2) = lgnt(2) + dsqrt( dxy(1)*dxy(1) + dxy(2)*dxy(2) )
1058 
1059  end subroutine updatecontactmultiplier_alag
1060 
1061  subroutine gettiedstiffness_alag(cstate, tSurf, mu, stiff, force)
1062 
1063  type(tcontactstate), intent(in) :: cstate
1064  type(tsurfelement), intent(in) :: tsurf
1065  real(kind=kreal), intent(in) :: mu
1066  real(kind=kreal), intent(out) :: stiff(:,:)
1067  real(kind=kreal), intent(out) :: force(:)
1068 
1069  integer :: i, j, nnode, edof
1070  real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
1071  real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
1072 
1073  nnode = size(tsurf%nodes)
1074  edof = nnode*3+3
1075 
1076  stiff = 0.d0
1077 
1078  ! Use common mapping routine to compute Tm
1079  call computetm_tt(cstate, tsurf, 0.0d0, tm, tt)
1080 
1081  ! Tied stiffness: stiff = mu * Tm^T * Tm
1082  do j = 1, edof
1083  do i = 1, edof
1084  stiff(i,j) = mu * dot_product(tm(1:3,i), tm(1:3,j))
1085  enddo
1086  enddo
1087  force(1:edof) = 0.d0 ! not used for tied (3-direction constraint)
1088 
1089  end subroutine gettiedstiffness_alag
1090 
1091  subroutine gettiednodalforce_alag(ctState,tSurf,ndu,mu,ctNForce,ctTForce)
1092 
1093  use msurfelement
1094  type(tcontactstate) :: ctstate
1095  type(tsurfelement) :: tsurf
1096  integer(kind=kint) :: nnode
1097  real(kind=kreal) :: ndu(:)
1098  real(kind=kreal), intent(in) :: mu
1099  real(kind=kreal) :: ctnforce(:)
1100  real(kind=kreal) :: cttforce(:)
1101 
1102  integer(kind=kint) :: edof
1103  real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
1104  real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))
1105  real(kind=kreal) :: dg(3)
1106  real(kind=kreal) :: nrlforce(3)
1107 
1108  nnode = size(tsurf%nodes)
1109  edof = nnode*3+3
1110 
1111  ctnforce = 0.0d0
1112  cttforce = 0.0d0
1113 
1114  ! Use common mapping routine to compute Tm
1115  call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt)
1116 
1117  ! Gap vector: dg = Tm * ndu (3-component relative displacement)
1118  dg(1:3) = matmul(tm(1:3, 1:edof), ndu(1:edof))
1119 
1120  ! Force: multiplier + penalty * gap (3 components)
1121  nrlforce(1:3) = ctstate%multiplier(1:3) + mu*dg(1:3)
1122 
1123  ! Distribute force: ctNForce = -Tm^T * nrlforce
1124  ctnforce(1:edof) = -matmul(transpose(tm(1:3, 1:edof)), nrlforce(1:3))
1125 
1126  ! Lagrange row (not used in ALagrange, set to 0)
1127  ctnforce((nnode+1)*3+1) = 0.d0
1128  cttforce((nnode+1)*3+1) = 0.d0
1129 
1130  end subroutine gettiednodalforce_alag
1131 
1132 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:979
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:589
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:1062
m_fstr_contact_elem_alag::gettangentslip
subroutine, public gettangentslip(slave_surf, master, coord, disp, ddisp, maplist, 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::gettiednodalforce_alag
subroutine, public gettiednodalforce_alag(ctState, tSurf, ndu, mu, ctNForce, ctTForce)
Definition: fstr_contact_elem_alag.f90:1092
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:414
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::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
subroutine, public getcontactnodalforce_alag(ctState, tSurf, ndCoord, ndDu, mu, mut, fcoeff, symm, ctNForce, ctTForce, smoothing_type)
Definition: fstr_contact_elem_alag.f90:897
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:766
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:504
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:99
msurfelement::tsurfelement
Structure to define surface group.
Definition: surf_ele.f90:23