FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
static_LIB_1d.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_static_lib_1d
7  use hecmw, only : kint, kreal
8  use elementinfo
9  use mmechgauss
10  use m_matmatrix
11  implicit none
12 
13 contains
14  !
15  !=====================================================================*
16  ! STF_C1
17  !=====================================================================*
19  subroutine stf_c1( etype,nn,ecoord,area,gausses,stiff, u ,temperature )
20  integer(kind=kint), intent(in) :: etype
21  integer(kind=kint), intent(in) :: nn
22  real(kind=kreal), intent(in) :: ecoord(3,nn)
23  real(kind=kreal), intent(in) :: area
24  type(tgaussstatus), intent(in) :: gausses(:)
25  real(kind=kreal), intent(out) :: stiff(:,:)
26  real(kind=kreal), intent(in), optional :: u(:,:)
27  real(kind=kreal), intent(in), optional :: temperature(nn)
28 
29  real(kind=kreal) llen, llen0, elem(3,nn)
30  logical :: ierr
31  real(kind=kreal) ina(1), outa(1), direc(3), direc0(3), coeff, strain
32  integer(kind=kint) :: i,j
33 
34  ! we suppose the same material type in the element
35  if( present(u) ) then
36  elem(:,:) = ecoord(:,:) + u(:,:)
37  else
38  elem(:,:) = ecoord(:,:)
39  endif
40 
41  direc = elem(:,2)-elem(:,1)
42  llen = dsqrt( dot_product(direc, direc) )
43  direc = direc/llen
44  direc0 = ecoord(:,2)-ecoord(:,1)
45  llen0 = dsqrt( dot_product(direc0, direc0) )
46 
47  if( present(temperature) ) then
48  ina(1) = 0.5d0*(temperature(1)+temperature(2))
49  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr, ina )
50  else
51  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr )
52  endif
53  if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
54  coeff = outa(1)*area*llen0/(llen*llen)
55  strain = gausses(1)%strain(1)
56 
57  stiff(:,:) = 0.d0
58  do i=1,3
59  stiff(i,i) = coeff*strain
60  do j=1,3
61  stiff(i,j) = stiff(i,j) + coeff*(1.d0-2.d0*strain)*direc(i)*direc(j)
62  enddo
63  enddo
64 
65  stiff(4:6,1:3) = -stiff(1:3,1:3)
66  stiff(1:3,4:6) = transpose(stiff(4:6,1:3))
67  stiff(4:6,4:6) = stiff(1:3,1:3)
68 
69  end subroutine stf_c1
70 
71  !
73  !---------------------------------------------------------------------*
74  subroutine update_c1( etype, nn, ecoord, area, u, du, qf ,gausses, TT, T0 )
75  !---------------------------------------------------------------------*
76  use m_fstr
77  ! I/F VARIABLES
78  integer(kind=kint), intent(in) :: etype
79  integer(kind=kint), intent(in) :: nn
80  real(kind=kreal), intent(in) :: ecoord(3,nn)
81  real(kind=kreal), intent(in) :: area
82  real(kind=kreal), intent(in) :: u(3,nn)
83  real(kind=kreal), intent(in) :: du(3,nn)
84  real(kind=kreal), intent(out) :: qf(nn*3)
85  type(tgaussstatus), intent(inout) :: gausses(:)
86  real(kind=kreal), intent(in), optional :: tt(nn)
87  real(kind=kreal), intent(in), optional :: t0(nn)
88 
89  ! LOCAL VARIABLES
90  real(kind=kreal) :: direc(3), direc0(3)
91  real(kind=kreal) :: llen, llen0, ina(1), outa(1)
92  real(kind=kreal) :: elem(3,nn)
93  real(kind=kreal) :: young
94  real(kind=kreal) :: ttc, tt0, alp, alp0, epsth
95  logical :: ierr
96 
97  qf(:) = 0.d0
98  ! we suppose the same material type in the element
99  elem(:,:) = ecoord(:,:) + u(:,:) + du(:,:)
100 
101  direc = elem(:,2)-elem(:,1)
102  llen = dsqrt( dot_product(direc, direc) )
103  direc = direc/llen
104  direc0 = ecoord(:,2)-ecoord(:,1)
105  llen0 = dsqrt( dot_product(direc0, direc0) )
106 
107  epsth = 0.d0
108  if( present(tt) .and. present(t0) ) then
109  ttc = 0.5d0*(tt(1)+tt(2))
110  tt0 = 0.5d0*(t0(1)+t0(2))
111 
112  ina(1) = ttc
113  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr, ina )
114  if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
115  young = outa(1)
116 
117  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
118  if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_exapnsion)
119  alp = outa(1)
120 
121  ina(1) = tt0
122  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa(:), ierr, ina )
123  if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_exapnsion)
124  alp0 = outa(1)
125 
126  epsth=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
127  else
128  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa, ierr )
129  if( ierr ) outa(1) = gausses(1)%pMaterial%variables(m_youngs)
130  young = outa(1)
131  endif
132 
133  gausses(1)%strain(1) = dlog(llen/llen0)
134  gausses(1)%stress(1) = young*(gausses(1)%strain(1)-epsth)
135 
136  !set stress and strain for output
137  gausses(1)%strain_out(1) = gausses(1)%strain(1)
138  gausses(1)%stress_out(1) = gausses(1)%stress(1)
139 
140  qf(1) = gausses(1)%stress(1)*area*llen0/llen
141  qf(1:3) = -qf(1)*direc
142  qf(4:6) = -qf(1:3)
143 
144  end subroutine update_c1
145 
146  !----------------------------------------------------------------------*
147  subroutine nodalstress_c1(ETYPE,NN,gausses,ndstrain,ndstress)
148  !----------------------------------------------------------------------*
149  !
150  ! Calculate Strain and Stress increment of solid elements
151  !
152  integer(kind=kint), intent(in) :: etype,nn
153  type(tgaussstatus), intent(in) :: gausses(:)
154  real(kind=kreal), intent(out) :: ndstrain(nn,6)
155  real(kind=kreal), intent(out) :: ndstress(nn,6)
156 
157  ndstrain(1,1:6) = gausses(1)%strain(1:6)
158  ndstress(1,1:6) = gausses(1)%stress(1:6)
159  ndstrain(2,1:6) = gausses(1)%strain(1:6)
160  ndstress(2,1:6) = gausses(1)%stress(1:6)
161 
162  end subroutine
163  !
164  !
165  !
166  !----------------------------------------------------------------------*
167  subroutine elementstress_c1(ETYPE,gausses,strain,stress)
168  !----------------------------------------------------------------------*
169  !
170  ! Calculate Strain and Stress increment of solid elements
171  !
172  integer(kind=kint), intent(in) :: ETYPE
173  type(tgaussstatus), intent(in) :: gausses(:)
174  real(kind=kreal), intent(out) :: strain(6)
175  real(kind=kreal), intent(out) :: stress(6)
176 
177 
178  strain(:) = gausses(1)%strain(1:6)
179  stress(:) = gausses(1)%stress(1:6)
180 
181  end subroutine
182 
183 
184  !----------------------------------------------------------------------*
185  subroutine dl_c1(etype, nn, xx, yy, zz, rho, thick, ltype, params, &
186  vect, nsize)
187  !----------------------------------------------------------------------*
188  !** SET DLOAD
189  ! GRAV LTYPE=4 :GRAVITY FORCE
190  ! I/F VARIABLES
191  integer(kind = kint), intent(in) :: etype, nn
192  real(kind = kreal), intent(in) :: xx(:), yy(:), zz(:)
193  real(kind = kreal), intent(in) :: params(0:6)
194  real(kind = kreal), intent(inout) :: vect(:)
195  real(kind = kreal) :: rho, thick
196  integer(kind = kint) :: ltype, nsize
197  ! LOCAL VARIABLES
198  integer(kind = kint) :: ndof = 3
199  integer(kind = kint) :: ivol, isuf, i
200  real(kind = kreal) :: vx, vy, vz, val, a, aa
201  !--------------------------------------------------------------------
202  val = params(0)
203  !--------------------------------------------------------------
204 
205  ivol = 0
206  isuf = 0
207  if( ltype .LT. 10 ) then
208  ivol = 1
209  else if( ltype .GE. 10 ) then
210  isuf = 1
211  end if
212 
213  !--------------------------------------------------------------------
214  nsize = nn*ndof
215  !--------------------------------------------------------------------
216  vect(1:nsize) = 0.0d0
217 
218  ! Volume force
219  if( ivol .EQ. 1 ) then
220  if( ltype .EQ. 4 ) then
221  aa = dsqrt( ( xx(2)-xx(1) )*( xx(2)-xx(1) ) &
222  +( yy(2)-yy(1) )*( yy(2)-yy(1) ) &
223  +( zz(2)-zz(1) )*( zz(2)-zz(1) ) )
224 
225  a = thick
226  vx = params(1)
227  vy = params(2)
228  vz = params(3)
229  vx = vx/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
230  vy = vy/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
231  vz = vz/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
232 
233  do i = 1, 2
234  vect(3*i-2) = val*rho*a*0.5d0*aa*vx
235  vect(3*i-1) = val*rho*a*0.5d0*aa*vy
236  vect(3*i ) = val*rho*a*0.5d0*aa*vz
237  end do
238 
239  end if
240  end if
241  !--------------------------------------------------------------------
242 
243  return
244  end subroutine
245 
246  !
247  !
248  !----------------------------------------------------------------------*
249  subroutine truss_diag_modify(hecMAT,hecMESH)
250  !----------------------------------------------------------------------*
251  !
252  use hecmw
253  type (hecmwST_matrix) :: hecMAT
254  type (hecmwST_local_mesh) :: hecMESH
255  integer(kind=kint) :: itype, is, iE, ic_type, icel, jS, j, n
256 
257  do itype = 1, hecmesh%n_elem_type
258  ic_type = hecmesh%elem_type_item(itype)
259  if(ic_type == 301)then
260  is = hecmesh%elem_type_index(itype-1) + 1
261  ie = hecmesh%elem_type_index(itype )
262  do icel = is, ie
263  js = hecmesh%elem_node_index(icel-1)
264  do j=1,2
265  n = hecmesh%elem_node_item(js+j)
266  if( hecmat%D(9*n-8) == 0.0d0)then
267  hecmat%D(9*n-8) = 1.0d0
268  !call search_diag_modify(n,1,hecMAT,hecMESH)
269  endif
270  if( hecmat%D(9*n-4) == 0.0d0)then
271  hecmat%D(9*n-4) = 1.0d0
272  !call search_diag_modify(n,2,hecMAT,hecMESH)
273  endif
274  if( hecmat%D(9*n ) == 0.0d0)then
275  hecmat%D(9*n ) = 1.0d0
276  !call search_diag_modify(n,3,hecMAT,hecMESH)
277  endif
278  enddo
279  enddo
280  endif
281  enddo
282 
283  end subroutine
284 
285  !
286  !
287  !----------------------------------------------------------------------*
288  subroutine search_diag_modify(n,nn,hecMAT,hecMESH)
289  !----------------------------------------------------------------------*
290  !
291  use hecmw
292  type (hecmwST_matrix) :: hecMAT
293  type (hecmwST_local_mesh) :: hecMESH
294  integer :: n, nn, is, iE, i, a
295 
296  if(nn == 1)then
297  a = 0
298  is = hecmat%IndexL(n-1)+1
299  ie = hecmat%IndexL(n )
300  do i=is,ie
301  if(hecmat%AL(9*i-8) /= 0.0d0) a = 1
302  if(hecmat%AL(9*i-7) /= 0.0d0) a = 1
303  if(hecmat%AL(9*i-6) /= 0.0d0) a = 1
304  enddo
305  is = hecmat%IndexU(n-1)+1
306  ie = hecmat%IndexU(n )
307  do i=is,ie
308  if(hecmat%AU(9*i-8) /= 0.0d0) a = 1
309  if(hecmat%AU(9*i-7) /= 0.0d0) a = 1
310  if(hecmat%AU(9*i-6) /= 0.0d0) a = 1
311  enddo
312  if(hecmat%D(9*n-7) /= 0.0d0) a = 1
313  if(hecmat%D(9*n-6) /= 0.0d0) a = 1
314  if(a == 0)then
315  hecmat%D(9*n-8) = 1.0d0
316  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:1"
317  endif
318  endif
319  if(nn == 2)then
320  a = 0
321  is = hecmat%IndexL(n-1)+1
322  ie = hecmat%IndexL(n )
323  do i=is,ie
324  if(hecmat%AL(9*i-5) /= 0.0d0) a = 1
325  if(hecmat%AL(9*i-4) /= 0.0d0) a = 1
326  if(hecmat%AL(9*i-3) /= 0.0d0) a = 1
327  enddo
328  is = hecmat%IndexU(n-1)+1
329  ie = hecmat%IndexU(n )
330  do i=is,ie
331  if(hecmat%AU(9*i-5) /= 0.0d0) a = 1
332  if(hecmat%AU(9*i-4) /= 0.0d0) a = 1
333  if(hecmat%AU(9*i-3) /= 0.0d0) a = 1
334  enddo
335  if(hecmat%D(9*n-5) /= 0.0d0) a = 1
336  if(hecmat%D(9*n-3) /= 0.0d0) a = 1
337  if(a == 0)then
338  hecmat%D(9*n-4) = 1.0d0
339  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:2"
340  endif
341  endif
342  if(nn == 3)then
343  a = 0
344  is = hecmat%IndexL(n-1)+1
345  ie = hecmat%IndexL(n )
346  do i=is,ie
347  if(hecmat%AL(9*i-2) /= 0.0d0) a = 1
348  if(hecmat%AL(9*i-1) /= 0.0d0) a = 1
349  if(hecmat%AL(9*i ) /= 0.0d0) a = 1
350  enddo
351  is = hecmat%IndexU(n-1)+1
352  ie = hecmat%IndexU(n )
353  do i=is,ie
354  if(hecmat%AU(9*i-2) /= 0.0d0) a = 1
355  if(hecmat%AU(9*i-1) /= 0.0d0) a = 1
356  if(hecmat%AU(9*i ) /= 0.0d0) a = 1
357  enddo
358  if(hecmat%D(9*n-2) /= 0.0d0) a = 1
359  if(hecmat%D(9*n-1) /= 0.0d0) a = 1
360  if(a == 0)then
361  hecmat%D(9*n ) = 1.0d0
362  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:3"
363  endif
364  endif
365 
366  end subroutine
367 
369  subroutine stf_connector( etype, nn, ecoord, gausses, stiff, u, tincr, temperature )
370  integer(kind=kint), intent(in) :: etype
371  integer(kind=kint), intent(in) :: nn
372  real(kind=kreal), intent(in) :: ecoord(3,nn)
373  type(tgaussstatus), intent(in) :: gausses(:)
374  real(kind=kreal), intent(out) :: stiff(:,:)
375  real(kind=kreal), intent(in) :: u(:,:)
376  real(kind=kreal), intent(in) :: tincr
377  real(kind=kreal), intent(in) :: temperature(nn)
378 
379  integer(kind=kint) :: mtype, ctype
380  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
381 
382  stiff(:,:) = 0.d0
383 
384  ! if spring_d is assigned, add the stiffness to stiff
385  if( getnumofspring_dparam( gausses(1)%pMaterial ) > 0 ) then
386  call stf_spring_d( gausses, stiff )
387  endif
388 
389  ! if spring_a is assigned, add the stiffness to stiff
390  if( getnumofspring_aparam( gausses(1)%pMaterial ) > 0 ) then
391  call stf_spring_a( gausses, ecoord, u, stiff )
392  endif
393 
394  ! if dashpot_d is assigned, add the stiffness to stiff
395  if( getnumofdashpot_dparam( gausses(1)%pMaterial ) > 0 ) then
396  call stf_dashpot_d( gausses, tincr, stiff )
397  endif
398 
399  ! if dashpot_a is assigned, add the stiffness to stiff
400  if( getnumofdashpot_aparam( gausses(1)%pMaterial ) > 0 ) then
401  call stf_dashpot_a( gausses, ecoord, u, tincr, stiff )
402  endif
403 
404  end subroutine
405 
406  !
408  !---------------------------------------------------------------------*
409  subroutine update_connector( etype, nn, ecoord, u, du, qf ,gausses, tincr, TT )
410  !---------------------------------------------------------------------*
411  use m_fstr
412  ! I/F VARIABLES
413  integer(kind=kint), intent(in) :: etype
414  integer(kind=kint), intent(in) :: nn
415  real(kind=kreal), intent(in) :: ecoord(3,nn)
416  real(kind=kreal), intent(in) :: u(3,nn)
417  real(kind=kreal), intent(in) :: du(3,nn)
418  real(kind=kreal), intent(out) :: qf(nn*3)
419  type(tgaussstatus), intent(inout) :: gausses(:)
420  real(kind=kreal), intent(in) :: tincr
421  real(kind=kreal), intent(in), optional :: tt(nn)
422 
423  ! LOCAL VARIABLES
424  integer(kind=kint) :: mtype, ctype
425  real(kind=kreal) llen, llen0, elem(3,nn)
426  real(kind=kreal) direc(3), direc0(3), ratio
427  real(kind=kreal) :: params(4)
428  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
429 
430  qf = 0.d0
431  !clear stress and strain
432  gausses(1)%strain(:) = 0.d0
433  gausses(1)%stress(:) = 0.d0
434 
435  ! if spring_d is assigned, add the equivalent force to qf
436  if( getnumofspring_dparam( gausses(1)%pMaterial ) > 0 ) then
437  call update_spring_d( gausses, ecoord, u, du, qf )
438  endif
439 
440  ! if spring_a is assigned, add the stiffness to stiff
441  if( getnumofspring_aparam( gausses(1)%pMaterial ) > 0 ) then
442  call update_spring_a( gausses, ecoord, u, du, qf )
443  endif
444 
445  ! if dashpot_d is assigned, add the stiffness to stiff
446  if( getnumofdashpot_dparam( gausses(1)%pMaterial ) > 0 ) then
447  call update_dashpot_d( gausses, ecoord, du, tincr, qf )
448  endif
449 
450  ! if dashpot_a is assigned, add the stiffness to stiff
451  if( getnumofdashpot_aparam( gausses(1)%pMaterial ) > 0 ) then
452  call update_dashpot_a( gausses, ecoord, u, du, tincr, qf )
453  endif
454 
455  !set stress and strain for output
456  gausses(1)%strain_out(:) = gausses(1)%strain(:)
457  gausses(1)%stress_out(:) = gausses(1)%stress(:)
458 
459  end subroutine
460 
461  subroutine stf_spring_d( gausses, stiff )
462  type(tgaussstatus), intent(in) :: gausses(:)
463  real(kind=kreal), intent(out) :: stiff(:,:)
464 
465  real(kind=kreal) :: params(1)
466  integer(kind=kint) :: iparams(2)
467  integer(kind=kint) :: i, j, n_ndof, id1, id2
468 
469  n_ndof = getnumofspring_dparam( gausses(1)%pMaterial )
470 
471  do i=1,n_ndof
472  call getconnectorproperty( gausses(1), m_spring_dof, i, params, iparams )
473  id1 = iparams(1)
474  id2 = 3+iparams(2)
475  stiff(id1,id1) = stiff(id1,id1) + params(1)
476  stiff(id1,id2) = stiff(id1,id2) - params(1)
477  stiff(id2,id1) = stiff(id2,id1) - params(1)
478  stiff(id2,id2) = stiff(id2,id2) + params(1)
479  enddo
480  end subroutine
481 
482  subroutine update_spring_d( gausses, ecoord, u, du, qf )
483  type(tgaussstatus), intent(inout) :: gausses(:)
484  real(kind=kreal), intent(in) :: ecoord(:,:)
485  real(kind=kreal), intent(in) :: u(:,:)
486  real(kind=kreal), intent(in) :: du(:,:)
487  real(kind=kreal), intent(out) :: qf(:)
488 
489  ! LOCAL VARIABLES
490  real(kind=kreal) :: params(1)
491  integer(kind=kint) :: iparams(2)
492  real(kind=kreal) :: totaldisp(3,2), stretch, sforce
493  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
494 
495  n_ndof = getnumofspring_dparam( gausses(1)%pMaterial )
496  totaldisp(1:3,1:2) = u(1:3,1:2) + du(1:3,1:2)
497 
498  do i=1,n_ndof
499  call getconnectorproperty( gausses(1), m_spring_dof, i, params, iparams )
500  dof1 = iparams(1)
501  dof2 = iparams(2)
502  id1 = dof1
503  id2 = 3+dof2
504  stretch = totaldisp(dof1,1)-totaldisp(dof2,2)
505  sforce = params(1) * stretch
506  if( dof1 == dof2 ) then
507  gausses(1)%strain(dof1) = gausses(1)%strain(dof1) + stretch
508  gausses(1)%stress(dof1) = gausses(1)%stress(dof1) + sforce
509  endif
510  qf(id1) = qf(id1) + sforce
511  qf(id2) = qf(id2) - sforce
512  enddo
513  end subroutine
514 
515  subroutine stf_spring_a( gausses, ecoord, u, stiff )
516  type(tgaussstatus), intent(in) :: gausses(:)
517  real(kind=kreal), intent(in) :: ecoord(:,:)
518  real(kind=kreal), intent(in) :: u(:,:)
519  real(kind=kreal), intent(out) :: stiff(:,:)
520 
521  real(kind=kreal) :: params(1)
522  integer(kind=kint) :: iparams(2)
523  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
524  real(kind=kreal) :: llen, llen0, elem(3,2)
525  real(kind=kreal) :: direc(3), direc0(3), ratio
526 
527  ! Retrieve connector properties using the first Gauss point
528  call getconnectorproperty( gausses(1), m_spring_axial, 0, params, iparams )
529 
530  ! Calculate the current positions of the elemental nodes after displacement
531  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2)
532 
533  ! Compute the direction vector for the deformed element
534  direc = elem(1:3,1) - elem(1:3,2)
535  llen = dsqrt( dot_product(direc, direc) )
536 
537  ! Compute the direction vector for the undeformed element
538  direc0 = ecoord(1:3,1) - ecoord(1:3,2)
539  llen0 = dsqrt( dot_product(direc0, direc0) )
540 
541  ! Check for near-zero length to avoid division by zero
542  if( llen < 1.d-10 ) then
543  direc(1:3) = 0.d0
544  ratio = 1.d0
545  else
546  direc(1:3) = direc(1:3) / llen
547  ratio = (llen - llen0) / llen
548  endif
549 
550  ! Populate the stiffness matrix
551  do i = 1, 3
552  stiff(i,i) = stiff(i,i) + params(1) * ratio
553  do j = 1, 3
554  stiff(i,j) = stiff(i,j) + params(1) * (1.d0 - ratio) * direc(i) * direc(j)
555  enddo
556  enddo
557 
558  ! Symmetric stiffness matrix assembly for the other degrees of freedom
559  stiff(4:6, 1:3) = -stiff(1:3, 1:3)
560  stiff(1:3, 4:6) = transpose(stiff(4:6, 1:3))
561  stiff(4:6, 4:6) = stiff(1:3, 1:3)
562 
563  end subroutine
564 
565  subroutine update_spring_a( gausses, ecoord, u, du, qf )
566  type(tgaussstatus), intent(inout) :: gausses(:)
567  real(kind=kreal), intent(in) :: ecoord(:,:)
568  real(kind=kreal), intent(in) :: u(:,:)
569  real(kind=kreal), intent(in) :: du(:,:)
570  real(kind=kreal), intent(out) :: qf(:)
571 
572  ! LOCAL VARIABLES
573  real(kind=kreal) :: params(1)
574  integer(kind=kint) :: iparams(2)
575  real(kind=kreal) :: llen, llen0, elem(3,2)
576  real(kind=kreal) :: direc(3), direc0(3)
577 
578  ! Retrieve connector properties for the first Gauss point
579  call getconnectorproperty( gausses(1), m_spring_axial, 0, params, iparams )
580 
581  ! Calculate the current positions of the elemental nodes after applying displacements
582  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2) + du(1:3,1:2)
583 
584  ! Compute the direction vector for the deformed element
585  direc = elem(1:3,1) - elem(1:3,2)
586  llen = dsqrt( dot_product(direc, direc) ) ! Current length of the deformed element
587 
588  ! Compute the direction vector for the undeformed element
589  direc0 = ecoord(1:3,1) - ecoord(1:3,2)
590  llen0 = dsqrt( dot_product(direc0, direc0) ) ! Original length of the element
591 
592  ! Check for near-zero length to avoid division by zero
593  if( llen < 1.d-10 ) then
594  direc(1:3) = 1.d0 ! Set direction to a default value to avoid undefined behavior
595  else
596  direc(1:3) = direc(1:3) / llen ! Normalize the direction vector to unit length
597  endif
598 
599  ! Update strain and stress for the first Gauss point
600  gausses(1)%strain(1) = llen - llen0 ! Calculate strain as the change in length
601  gausses(1)%stress(1) = params(1) * gausses(1)%strain(1) ! Calculate stress using material properties
602 
603  ! Update the internal force vector based on calculated stress
604  qf(1:3) = qf(1:3) + gausses(1)%stress(1)*direc(1:3) ! Add stress contribution to the first three components
605  qf(4:6) = qf(4:6) - gausses(1)%stress(1)*direc(1:3) ! Subtract stress contribution for the next three components
606  end subroutine
607 
608  subroutine stf_dashpot_d( gausses, tincr, stiff )
609  type(tgaussstatus), intent(in) :: gausses(:)
610  real(kind=kreal), intent(in) :: tincr
611  real(kind=kreal), intent(out) :: stiff(:,:)
612 
613  real(kind=kreal) :: params(1)
614  integer(kind=kint) :: iparams(2)
615  integer(kind=kint) :: i, j, n_ndof, id1, id2
616 
617  n_ndof = getnumofdashpot_dparam( gausses(1)%pMaterial )
618 
619  do i=1,n_ndof
620  call getconnectorproperty( gausses(1), m_dashpot_dof, i, params, iparams )
621  id1 = iparams(1)
622  id2 = 3+iparams(2)
623  stiff(id1,id1) = stiff(id1,id1) + params(1)/tincr
624  stiff(id1,id2) = stiff(id1,id2) - params(1)/tincr
625  stiff(id2,id1) = stiff(id2,id1) - params(1)/tincr
626  stiff(id2,id2) = stiff(id2,id2) + params(1)/tincr
627  enddo
628  end subroutine
629 
630  subroutine update_dashpot_d( gausses, ecoord, du, tincr, qf )
631  type(tgaussstatus), intent(inout) :: gausses(:)
632  real(kind=kreal), intent(in) :: ecoord(:,:)
633  real(kind=kreal), intent(in) :: du(:,:)
634  real(kind=kreal), intent(in) :: tincr
635  real(kind=kreal), intent(out) :: qf(:)
636 
637  ! LOCAL VARIABLES
638  real(kind=kreal) :: params(1)
639  integer(kind=kint) :: iparams(2)
640  real(kind=kreal) :: totaldisp(3,2), stretch, dforce
641  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
642 
643  n_ndof = getnumofdashpot_dparam( gausses(1)%pMaterial )
644 
645  do i=1,n_ndof
646  call getconnectorproperty( gausses(1), m_dashpot_dof, i, params, iparams )
647  dof1 = iparams(1)
648  dof2 = iparams(2)
649  id1 = dof1
650  id2 = 3+dof2
651  stretch = du(dof1,1)-du(dof2,2)
652  dforce = params(1) * stretch / tincr
653  if( dof1 == dof2 ) then
654  gausses(1)%strain(dof1) = gausses(1)%strain(dof1) + stretch
655  gausses(1)%stress(dof1) = gausses(1)%stress(dof1) + dforce
656  endif
657  qf(id1) = qf(id1) + dforce
658  qf(id2) = qf(id2) - dforce
659  enddo
660  end subroutine
661 
662  subroutine stf_dashpot_a( gausses, ecoord, u, tincr, stiff )
663  type(tgaussstatus), intent(in) :: gausses(:)
664  real(kind=kreal), intent(in) :: ecoord(:,:)
665  real(kind=kreal), intent(in) :: u(:,:)
666  real(kind=kreal), intent(in) :: tincr
667  real(kind=kreal), intent(out) :: stiff(:,:)
668 
669  real(kind=kreal) :: params(1)
670  integer(kind=kint) :: iparams(2)
671  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
672  real(kind=kreal) :: llen, elem(3,2)
673  real(kind=kreal) :: direc(3), ratio
674 
675  ! Retrieve connector properties using the first Gauss point
676  call getconnectorproperty( gausses(1), m_dashpot_axial, 0, params, iparams )
677 
678  ! Calculate the current positions of the elemental nodes after displacement
679  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2)
680 
681  ! Compute the direction vector for the deformed element
682  direc = elem(1:3,1) - elem(1:3,2)
683  llen = dsqrt( dot_product(direc, direc) )
684 
685  ! Check for near-zero length to avoid division by zero
686  if( llen < 1.d-10 ) then
687  direc(1:3) = 0.d0
688  ratio = 1.d0
689  else
690  direc(1:3) = direc(1:3) / llen
691  ratio = gausses(1)%strain(1) / llen
692  endif
693 
694  ! Populate the stiffness matrix
695  do i = 1, 3
696  stiff(i,i) = stiff(i,i) + params(1) * ratio / tincr
697  do j = 1, 3
698  stiff(i,j) = stiff(i,j) + params(1) * (1.d0 - ratio) * direc(i) * direc(j) / tincr
699  enddo
700  enddo
701 
702  ! Symmetric stiffness matrix assembly for the other degrees of freedom
703  stiff(4:6, 1:3) = -stiff(1:3, 1:3)
704  stiff(1:3, 4:6) = transpose(stiff(4:6, 1:3))
705  stiff(4:6, 4:6) = stiff(1:3, 1:3)
706 
707  end subroutine
708 
709  subroutine update_dashpot_a( gausses, ecoord, u, du, tincr, qf )
710  type(tgaussstatus), intent(inout) :: gausses(:)
711  real(kind=kreal), intent(in) :: ecoord(:,:)
712  real(kind=kreal), intent(in) :: u(:,:)
713  real(kind=kreal), intent(in) :: du(:,:)
714  real(kind=kreal), intent(in) :: tincr
715  real(kind=kreal), intent(out) :: qf(:)
716 
717  ! LOCAL VARIABLES
718  real(kind=kreal) :: params(1)
719  integer(kind=kint) :: iparams(2)
720  real(kind=kreal) :: llen, llen0, elem(3,2)
721  real(kind=kreal) :: direc(3), direc0(3)
722 
723  ! Retrieve connector properties for the first Gauss point
724  call getconnectorproperty( gausses(1), m_dashpot_axial, 0, params, iparams )
725 
726  ! Calculate the current positions of the elemental nodes after applying displacements
727  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2) + du(1:3,1:2)
728 
729  ! Compute the direction vector for the deformed element
730  direc = elem(1:3,1) - elem(1:3,2)
731  llen = dsqrt( dot_product(direc, direc) ) ! Current length of the deformed element
732 
733  ! Compute the direction vector for the deformed element at the start of increment
734  direc0 = ecoord(1:3,1) + u(1:3,1) - ecoord(1:3,2) - u(1:3,2)
735  llen0 = dsqrt( dot_product(direc0, direc0) ) ! Original length of the element
736 
737  ! Check for near-zero length to avoid division by zero
738  if( llen < 1.d-10 ) then
739  direc(1:3) = 1.d0 ! Set direction to a default value to avoid undefined behavior
740  else
741  direc(1:3) = direc(1:3) / llen ! Normalize the direction vector to unit length
742  endif
743 
744  ! Update strain and stress for the first Gauss point
745  gausses(1)%strain(1) = llen - llen0 ! Calculate strain as the change in length
746  gausses(1)%stress(1) = params(1) * gausses(1)%strain(1) / tincr ! Calculate stress using material properties
747 
748  ! Update the internal force vector based on calculated stress
749  qf(1:3) = qf(1:3) + gausses(1)%stress(1)*direc(1:3) ! Add stress contribution to the first three components
750  qf(4:6) = qf(4:6) - gausses(1)%stress(1)*direc(1:3) ! Subtract stress contribution for the next three components
751  end subroutine
752 
753 
754 end module m_static_lib_1d
m_static_lib_1d::update_spring_d
subroutine update_spring_d(gausses, ecoord, u, du, qf)
Definition: static_LIB_1d.f90:483
m_matmatrix::getconnectorproperty
subroutine getconnectorproperty(gauss, ctype, dofid, rparams, iparams, stretch)
Definition: calMatMatrix.f90:339
m_static_lib_1d::dl_c1
subroutine dl_c1(etype, nn, xx, yy, zz, rho, thick, ltype, params, vect, nsize)
Definition: static_LIB_1d.f90:187
m_static_lib_1d::update_spring_a
subroutine update_spring_a(gausses, ecoord, u, du, qf)
Definition: static_LIB_1d.f90:566
m_static_lib_1d::update_dashpot_d
subroutine update_dashpot_d(gausses, ecoord, du, tincr, qf)
Definition: static_LIB_1d.f90:631
m_static_lib_1d::stf_dashpot_a
subroutine stf_dashpot_a(gausses, ecoord, u, tincr, stiff)
Definition: static_LIB_1d.f90:663
m_static_lib_1d::truss_diag_modify
subroutine truss_diag_modify(hecMAT, hecMESH)
Definition: static_LIB_1d.f90:250
m_static_lib_1d::stf_c1
subroutine stf_c1(etype, nn, ecoord, area, gausses, stiff, u, temperature)
This subroutine calculate stiff matrix of 2-nodes truss element.
Definition: static_LIB_1d.f90:20
m_static_lib_1d::elementstress_c1
subroutine elementstress_c1(ETYPE, gausses, strain, stress)
Definition: static_LIB_1d.f90:168
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
m_static_lib_1d::stf_spring_d
subroutine stf_spring_d(gausses, stiff)
Definition: static_LIB_1d.f90:462
m_static_lib_1d::update_c1
subroutine update_c1(etype, nn, ecoord, area, u, du, qf, gausses, TT, T0)
Update strain and stress inside element.
Definition: static_LIB_1d.f90:75
m_static_lib_1d::nodalstress_c1
subroutine nodalstress_c1(ETYPE, NN, gausses, ndstrain, ndstress)
Definition: static_LIB_1d.f90:148
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_static_lib_1d::update_dashpot_a
subroutine update_dashpot_a(gausses, ecoord, u, du, tincr, qf)
Definition: static_LIB_1d.f90:710
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:136
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
hecmw
Definition: hecmw.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
m_static_lib_1d::search_diag_modify
subroutine search_diag_modify(n, nn, hecMAT, hecMESH)
Definition: static_LIB_1d.f90:289
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_static_lib_1d::stf_dashpot_d
subroutine stf_dashpot_d(gausses, tincr, stiff)
Definition: static_LIB_1d.f90:609
m_static_lib_1d::stf_spring_a
subroutine stf_spring_a(gausses, ecoord, u, stiff)
Definition: static_LIB_1d.f90:516
m_static_lib_1d
This module provide common functions of 3D truss elements.
Definition: static_LIB_1d.f90:6
m_static_lib_1d::stf_connector
subroutine stf_connector(etype, nn, ecoord, gausses, stiff, u, tincr, temperature)
This subroutine calculate stiff matrix of 2-nodes connector element.
Definition: static_LIB_1d.f90:370
m_static_lib_1d::update_connector
subroutine update_connector(etype, nn, ecoord, u, du, qf, gausses, tincr, TT)
Update strain and stress inside spring element.
Definition: static_LIB_1d.f90:410