FrontISTR  5.9.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  subroutine elementstress_c1(ETYPE,gausses,strain,stress)
166  !----------------------------------------------------------------------*
167  !
168  ! Calculate Strain and Stress increment of solid elements
169  !
170  integer(kind=kint), intent(in) :: ETYPE
171  type(tgaussstatus), intent(in) :: gausses(:)
172  real(kind=kreal), intent(out) :: strain(6)
173  real(kind=kreal), intent(out) :: stress(6)
174 
175 
176  strain(:) = gausses(1)%strain(1:6)
177  stress(:) = gausses(1)%stress(1:6)
178 
179  end subroutine
180 
181  !----------------------------------------------------------------------*
182  subroutine dl_c1(etype, nn, xx, yy, zz, rho, thick, ltype, params, &
183  vect, nsize)
184  !----------------------------------------------------------------------*
185  !** SET DLOAD
186  ! GRAV LTYPE=4 :GRAVITY FORCE
187  ! I/F VARIABLES
188  integer(kind = kint), intent(in) :: etype, nn
189  real(kind = kreal), intent(in) :: xx(:), yy(:), zz(:)
190  real(kind = kreal), intent(in) :: params(0:6)
191  real(kind = kreal), intent(inout) :: vect(:)
192  real(kind = kreal) :: rho, thick
193  integer(kind = kint) :: ltype, nsize
194  ! LOCAL VARIABLES
195  integer(kind = kint) :: ndof = 3
196  integer(kind = kint) :: ivol, isuf, i
197  real(kind = kreal) :: vx, vy, vz, val, a, aa
198  !--------------------------------------------------------------------
199  val = params(0)
200  !--------------------------------------------------------------
201 
202  ivol = 0
203  isuf = 0
204  if( ltype .LT. 10 ) then
205  ivol = 1
206  else if( ltype .GE. 10 ) then
207  isuf = 1
208  end if
209 
210  !--------------------------------------------------------------------
211  nsize = nn*ndof
212  !--------------------------------------------------------------------
213  vect(1:nsize) = 0.0d0
214 
215  ! Volume force
216  if( ivol .EQ. 1 ) then
217  if( ltype .EQ. 4 ) then
218  aa = dsqrt( ( xx(2)-xx(1) )*( xx(2)-xx(1) ) &
219  +( yy(2)-yy(1) )*( yy(2)-yy(1) ) &
220  +( zz(2)-zz(1) )*( zz(2)-zz(1) ) )
221 
222  a = thick
223  vx = params(1)
224  vy = params(2)
225  vz = params(3)
226  vx = vx/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
227  vy = vy/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
228  vz = vz/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
229 
230  do i = 1, 2
231  vect(3*i-2) = val*rho*a*0.5d0*aa*vx
232  vect(3*i-1) = val*rho*a*0.5d0*aa*vy
233  vect(3*i ) = val*rho*a*0.5d0*aa*vz
234  end do
235 
236  end if
237  end if
238  !--------------------------------------------------------------------
239 
240  return
241  end subroutine
242 
243  !----------------------------------------------------------------------*
244  subroutine truss_diag_modify(hecMAT,hecMESH)
245  !----------------------------------------------------------------------*
246  !
247  use hecmw
248  type (hecmwST_matrix) :: hecMAT
249  type (hecmwST_local_mesh) :: hecMESH
250  integer(kind=kint) :: itype, is, iE, ic_type, icel, jS, j, n
251 
252  do itype = 1, hecmesh%n_elem_type
253  ic_type = hecmesh%elem_type_item(itype)
254  if(ic_type == 301)then
255  is = hecmesh%elem_type_index(itype-1) + 1
256  ie = hecmesh%elem_type_index(itype )
257  do icel = is, ie
258  js = hecmesh%elem_node_index(icel-1)
259  do j=1,2
260  n = hecmesh%elem_node_item(js+j)
261  if( hecmat%D(9*n-8) == 0.0d0)then
262  hecmat%D(9*n-8) = 1.0d0
263  !call search_diag_modify(n,1,hecMAT,hecMESH)
264  endif
265  if( hecmat%D(9*n-4) == 0.0d0)then
266  hecmat%D(9*n-4) = 1.0d0
267  !call search_diag_modify(n,2,hecMAT,hecMESH)
268  endif
269  if( hecmat%D(9*n ) == 0.0d0)then
270  hecmat%D(9*n ) = 1.0d0
271  !call search_diag_modify(n,3,hecMAT,hecMESH)
272  endif
273  enddo
274  enddo
275  endif
276  enddo
277 
278  end subroutine
279 
280  !----------------------------------------------------------------------*
281  subroutine search_diag_modify(n,nn,hecMAT,hecMESH)
282  !----------------------------------------------------------------------*
283  !
284  use hecmw
285  type (hecmwST_matrix) :: hecMAT
286  type (hecmwST_local_mesh) :: hecMESH
287  integer :: n, nn, is, iE, i, a
288 
289  if(nn == 1)then
290  a = 0
291  is = hecmat%IndexL(n-1)+1
292  ie = hecmat%IndexL(n )
293  do i=is,ie
294  if(hecmat%AL(9*i-8) /= 0.0d0) a = 1
295  if(hecmat%AL(9*i-7) /= 0.0d0) a = 1
296  if(hecmat%AL(9*i-6) /= 0.0d0) a = 1
297  enddo
298  is = hecmat%IndexU(n-1)+1
299  ie = hecmat%IndexU(n )
300  do i=is,ie
301  if(hecmat%AU(9*i-8) /= 0.0d0) a = 1
302  if(hecmat%AU(9*i-7) /= 0.0d0) a = 1
303  if(hecmat%AU(9*i-6) /= 0.0d0) a = 1
304  enddo
305  if(hecmat%D(9*n-7) /= 0.0d0) a = 1
306  if(hecmat%D(9*n-6) /= 0.0d0) a = 1
307  if(a == 0)then
308  hecmat%D(9*n-8) = 1.0d0
309  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:1"
310  endif
311  endif
312  if(nn == 2)then
313  a = 0
314  is = hecmat%IndexL(n-1)+1
315  ie = hecmat%IndexL(n )
316  do i=is,ie
317  if(hecmat%AL(9*i-5) /= 0.0d0) a = 1
318  if(hecmat%AL(9*i-4) /= 0.0d0) a = 1
319  if(hecmat%AL(9*i-3) /= 0.0d0) a = 1
320  enddo
321  is = hecmat%IndexU(n-1)+1
322  ie = hecmat%IndexU(n )
323  do i=is,ie
324  if(hecmat%AU(9*i-5) /= 0.0d0) a = 1
325  if(hecmat%AU(9*i-4) /= 0.0d0) a = 1
326  if(hecmat%AU(9*i-3) /= 0.0d0) a = 1
327  enddo
328  if(hecmat%D(9*n-5) /= 0.0d0) a = 1
329  if(hecmat%D(9*n-3) /= 0.0d0) a = 1
330  if(a == 0)then
331  hecmat%D(9*n-4) = 1.0d0
332  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:2"
333  endif
334  endif
335  if(nn == 3)then
336  a = 0
337  is = hecmat%IndexL(n-1)+1
338  ie = hecmat%IndexL(n )
339  do i=is,ie
340  if(hecmat%AL(9*i-2) /= 0.0d0) a = 1
341  if(hecmat%AL(9*i-1) /= 0.0d0) a = 1
342  if(hecmat%AL(9*i ) /= 0.0d0) a = 1
343  enddo
344  is = hecmat%IndexU(n-1)+1
345  ie = hecmat%IndexU(n )
346  do i=is,ie
347  if(hecmat%AU(9*i-2) /= 0.0d0) a = 1
348  if(hecmat%AU(9*i-1) /= 0.0d0) a = 1
349  if(hecmat%AU(9*i ) /= 0.0d0) a = 1
350  enddo
351  if(hecmat%D(9*n-2) /= 0.0d0) a = 1
352  if(hecmat%D(9*n-1) /= 0.0d0) a = 1
353  if(a == 0)then
354  hecmat%D(9*n ) = 1.0d0
355  !write(*,"(a,i,a,i,a)")"### FIX DIAGONAL n:",n,", ID:",hecMESH%global_node_ID(n),", dof:3"
356  endif
357  endif
358 
359  end subroutine
360 
362  subroutine stf_connector( etype, nn, ecoord, gausses, stiff, u, temperature )
363  integer(kind=kint), intent(in) :: etype
364  integer(kind=kint), intent(in) :: nn
365  real(kind=kreal), intent(in) :: ecoord(3,nn)
366  type(tgaussstatus), intent(in) :: gausses(:)
367  real(kind=kreal), intent(out) :: stiff(:,:)
368  real(kind=kreal), intent(in) :: u(:,:)
369  real(kind=kreal), intent(in) :: temperature(nn)
370 
371  integer(kind=kint) :: mtype, ctype
372  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
373 
374  stiff(:,:) = 0.d0
375 
376  ! if spring_d is assigned, add the stiffness to stiff
377  if( getnumofspring_dparam( gausses(1)%pMaterial ) > 0 ) then
378  call stf_spring_d( gausses, stiff )
379  endif
380 
381  ! if spring_a is assigned, add the stiffness to stiff
382  if( getnumofspring_aparam( gausses(1)%pMaterial ) > 0 ) then
383  call stf_spring_a( gausses, ecoord, u, stiff )
384  endif
385 
386  end subroutine
387 
389  subroutine dmp_connector( etype, nn, ecoord, gausses, damping, u, temperature )
390  integer(kind=kint), intent(in) :: etype
391  integer(kind=kint), intent(in) :: nn
392  real(kind=kreal), intent(in) :: ecoord(3,nn)
393  type(tgaussstatus), intent(in) :: gausses(:)
394  real(kind=kreal), intent(out) :: damping(:,:)
395  real(kind=kreal), intent(in) :: u(:,:)
396  real(kind=kreal), intent(in) :: temperature(nn)
397 
398  integer(kind=kint) :: mtype, ctype
399  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
400 
401  damping(:,:) = 0.d0
402 
403  ! if dashpot_d is assigned, calcurate the damping matrix
404  if( getnumofdashpot_dparam( gausses(1)%pMaterial ) > 0 ) then
405  call stf_dashpot_d( gausses, damping )
406  endif
407 
408  ! if dashpot_a is assigned, calcurate the damping matrix
409  if( getnumofdashpot_aparam( gausses(1)%pMaterial ) > 0 ) then
410  call stf_dashpot_a( gausses, ecoord, u, damping )
411  endif
412 
413  end subroutine
414 
416  !---------------------------------------------------------------------*
417  subroutine update_connector( etype, nn, ecoord, u, du, qf ,gausses, TT )
418  !---------------------------------------------------------------------*
419  use m_fstr
420  ! I/F VARIABLES
421  integer(kind=kint), intent(in) :: etype
422  integer(kind=kint), intent(in) :: nn
423  real(kind=kreal), intent(in) :: ecoord(3,nn)
424  real(kind=kreal), intent(in) :: u(3,nn)
425  real(kind=kreal), intent(in) :: du(3,nn)
426  real(kind=kreal), intent(out) :: qf(nn*3)
427  type(tgaussstatus), intent(inout) :: gausses(:)
428  real(kind=kreal), intent(in), optional :: tt(nn)
429 
430  ! LOCAL VARIABLES
431  integer(kind=kint) :: mtype, ctype
432  real(kind=kreal) llen, llen0, elem(3,nn)
433  real(kind=kreal) direc(3), direc0(3), ratio
434  real(kind=kreal) :: params(4)
435  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
436 
437  qf = 0.d0
438  !clear stress and strain
439  gausses(1)%strain(:) = 0.d0
440  gausses(1)%stress(:) = 0.d0
441 
442  ! if spring_d is assigned, add the equivalent force to qf
443  if( getnumofspring_dparam( gausses(1)%pMaterial ) > 0 ) then
444  call update_spring_d( gausses, ecoord, u, du, qf )
445  endif
446 
447  ! if spring_a is assigned, add the stiffness to stiff
448  if( getnumofspring_aparam( gausses(1)%pMaterial ) > 0 ) then
449  call update_spring_a( gausses, ecoord, u, du, qf )
450  endif
451 
452  !set stress and strain for output
453  gausses(1)%strain_out(:) = gausses(1)%strain(:)
454  gausses(1)%stress_out(:) = gausses(1)%stress(:)
455 
456  end subroutine
457 
458  subroutine stf_spring_d( gausses, stiff )
459  type(tgaussstatus), intent(in) :: gausses(:)
460  real(kind=kreal), intent(out) :: stiff(:,:)
461 
462  real(kind=kreal) :: params(1)
463  integer(kind=kint) :: iparams(2)
464  integer(kind=kint) :: i, j, n_ndof, id1, id2
465 
466  n_ndof = getnumofspring_dparam( gausses(1)%pMaterial )
467 
468  do i=1,n_ndof
469  call getconnectorproperty( gausses(1), m_spring_dof, i, params, iparams )
470  id1 = iparams(1)
471  id2 = 3+iparams(2)
472  stiff(id1,id1) = stiff(id1,id1) + params(1)
473  stiff(id1,id2) = stiff(id1,id2) - params(1)
474  stiff(id2,id1) = stiff(id2,id1) - params(1)
475  stiff(id2,id2) = stiff(id2,id2) + params(1)
476  enddo
477  end subroutine
478 
479  subroutine update_spring_d( gausses, ecoord, u, du, qf )
480  type(tgaussstatus), intent(inout) :: gausses(:)
481  real(kind=kreal), intent(in) :: ecoord(:,:)
482  real(kind=kreal), intent(in) :: u(:,:)
483  real(kind=kreal), intent(in) :: du(:,:)
484  real(kind=kreal), intent(out) :: qf(:)
485 
486  ! LOCAL VARIABLES
487  real(kind=kreal) :: params(1)
488  integer(kind=kint) :: iparams(2)
489  real(kind=kreal) :: totaldisp(3,2), stretch, sforce
490  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
491 
492  n_ndof = getnumofspring_dparam( gausses(1)%pMaterial )
493  totaldisp(1:3,1:2) = u(1:3,1:2) + du(1:3,1:2)
494 
495  do i=1,n_ndof
496  call getconnectorproperty( gausses(1), m_spring_dof, i, params, iparams )
497  dof1 = iparams(1)
498  dof2 = iparams(2)
499  id1 = dof1
500  id2 = 3+dof2
501  stretch = totaldisp(dof1,1)-totaldisp(dof2,2)
502  sforce = params(1) * stretch
503  if( dof1 == dof2 ) then
504  gausses(1)%strain(dof1) = gausses(1)%strain(dof1) + stretch
505  gausses(1)%stress(dof1) = gausses(1)%stress(dof1) + sforce
506  endif
507  qf(id1) = qf(id1) + sforce
508  qf(id2) = qf(id2) - sforce
509  enddo
510  end subroutine
511 
512  subroutine stf_spring_a( gausses, ecoord, u, stiff )
513  type(tgaussstatus), intent(in) :: gausses(:)
514  real(kind=kreal), intent(in) :: ecoord(:,:)
515  real(kind=kreal), intent(in) :: u(:,:)
516  real(kind=kreal), intent(out) :: stiff(:,:)
517 
518  real(kind=kreal) :: params(1)
519  integer(kind=kint) :: iparams(2)
520  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
521  real(kind=kreal) :: llen, llen0, elem(3,2)
522  real(kind=kreal) :: direc(3), direc0(3), ratio
523 
524  ! Retrieve connector properties using the first Gauss point
525  call getconnectorproperty( gausses(1), m_spring_axial, 0, params, iparams )
526 
527  ! Calculate the current positions of the elemental nodes after displacement
528  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2)
529 
530  ! Compute the direction vector for the deformed element
531  direc = elem(1:3,1) - elem(1:3,2)
532  llen = dsqrt( dot_product(direc, direc) )
533 
534  ! Compute the direction vector for the undeformed element
535  direc0 = ecoord(1:3,1) - ecoord(1:3,2)
536  llen0 = dsqrt( dot_product(direc0, direc0) )
537 
538  ! Check for near-zero length to avoid division by zero
539  if( llen < 1.d-10 ) then
540  direc(1:3) = 0.d0
541  ratio = 1.d0
542  else
543  direc(1:3) = direc(1:3) / llen
544  ratio = (llen - llen0) / llen
545  endif
546 
547  ! Populate the stiffness matrix
548  do i = 1, 3
549  stiff(i,i) = stiff(i,i) + params(1) * ratio
550  do j = 1, 3
551  stiff(i,j) = stiff(i,j) + params(1) * (1.d0 - ratio) * direc(i) * direc(j)
552  enddo
553  enddo
554 
555  ! Symmetric stiffness matrix assembly for the other degrees of freedom
556  stiff(4:6, 1:3) = -stiff(1:3, 1:3)
557  stiff(1:3, 4:6) = transpose(stiff(4:6, 1:3))
558  stiff(4:6, 4:6) = stiff(1:3, 1:3)
559 
560  end subroutine
561 
562  subroutine update_spring_a( gausses, ecoord, u, du, qf )
563  type(tgaussstatus), intent(inout) :: gausses(:)
564  real(kind=kreal), intent(in) :: ecoord(:,:)
565  real(kind=kreal), intent(in) :: u(:,:)
566  real(kind=kreal), intent(in) :: du(:,:)
567  real(kind=kreal), intent(out) :: qf(:)
568 
569  ! LOCAL VARIABLES
570  real(kind=kreal) :: params(1)
571  integer(kind=kint) :: iparams(2)
572  real(kind=kreal) :: llen, llen0, elem(3,2)
573  real(kind=kreal) :: direc(3), direc0(3)
574 
575  ! Retrieve connector properties for the first Gauss point
576  call getconnectorproperty( gausses(1), m_spring_axial, 0, params, iparams )
577 
578  ! Calculate the current positions of the elemental nodes after applying displacements
579  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2) + du(1:3,1:2)
580 
581  ! Compute the direction vector for the deformed element
582  direc = elem(1:3,1) - elem(1:3,2)
583  llen = dsqrt( dot_product(direc, direc) ) ! Current length of the deformed element
584 
585  ! Compute the direction vector for the undeformed element
586  direc0 = ecoord(1:3,1) - ecoord(1:3,2)
587  llen0 = dsqrt( dot_product(direc0, direc0) ) ! Original length of the element
588 
589  ! Check for near-zero length to avoid division by zero
590  if( llen < 1.d-10 ) then
591  direc(1:3) = 1.d0 ! Set direction to a default value to avoid undefined behavior
592  else
593  direc(1:3) = direc(1:3) / llen ! Normalize the direction vector to unit length
594  endif
595 
596  ! Update strain and stress for the first Gauss point
597  gausses(1)%strain(1) = llen - llen0 ! Calculate strain as the change in length
598  gausses(1)%stress(1) = params(1) * gausses(1)%strain(1) ! Calculate stress using material properties
599 
600  ! Update the internal force vector based on calculated stress
601  qf(1:3) = qf(1:3) + gausses(1)%stress(1)*direc(1:3) ! Add stress contribution to the first three components
602  qf(4:6) = qf(4:6) - gausses(1)%stress(1)*direc(1:3) ! Subtract stress contribution for the next three components
603  end subroutine
604 
605  subroutine stf_dashpot_d( gausses, stiff )
606  type(tgaussstatus), intent(in) :: gausses(:)
607  real(kind=kreal), intent(out) :: stiff(:,:)
608 
609  real(kind=kreal) :: params(1)
610  integer(kind=kint) :: iparams(2)
611  integer(kind=kint) :: i, j, n_ndof, id1, id2
612 
613  n_ndof = getnumofdashpot_dparam( gausses(1)%pMaterial )
614 
615  do i=1,n_ndof
616  call getconnectorproperty( gausses(1), m_dashpot_dof, i, params, iparams )
617  id1 = iparams(1)
618  id2 = 3+iparams(2)
619  stiff(id1,id1) = stiff(id1,id1) + params(1)
620  stiff(id1,id2) = stiff(id1,id2) - params(1)
621  stiff(id2,id1) = stiff(id2,id1) - params(1)
622  stiff(id2,id2) = stiff(id2,id2) + params(1)
623  enddo
624  end subroutine
625 
626  subroutine stf_dashpot_a( gausses, ecoord, u, stiff )
627  type(tgaussstatus), intent(in) :: gausses(:)
628  real(kind=kreal), intent(in) :: ecoord(:,:)
629  real(kind=kreal), intent(in) :: u(:,:)
630  real(kind=kreal), intent(out) :: stiff(:,:)
631 
632  real(kind=kreal) :: params(1)
633  integer(kind=kint) :: iparams(2)
634  integer(kind=kint) :: i, j, n_ndof, dof1, dof2, id1, id2
635  real(kind=kreal) :: llen, elem(3,2)
636  real(kind=kreal) :: direc(3), ratio
637 
638  ! Retrieve connector properties using the first Gauss point
639  call getconnectorproperty( gausses(1), m_dashpot_axial, 0, params, iparams )
640 
641  ! Calculate the current positions of the elemental nodes after displacement
642  elem(1:3,1:2) = ecoord(1:3,1:2) + u(1:3,1:2)
643 
644  ! Compute the direction vector for the deformed element
645  direc = elem(1:3,1) - elem(1:3,2)
646  llen = dsqrt( dot_product(direc, direc) )
647 
648  ! Check for near-zero length to avoid division by zero
649  if( llen < 1.d-10 ) then
650  direc(1:3) = 0.d0
651  ratio = 1.d0
652  else
653  direc(1:3) = direc(1:3) / llen
654  ratio = gausses(1)%strain(1) / llen
655  endif
656 
657  ! Populate the stiffness matrix
658  do i = 1, 3
659  stiff(i,i) = stiff(i,i) + params(1) * ratio
660  do j = 1, 3
661  stiff(i,j) = stiff(i,j) + params(1) * (1.d0 - ratio) * direc(i) * direc(j)
662  enddo
663  enddo
664 
665  ! Symmetric stiffness matrix assembly for the other degrees of freedom
666  stiff(4:6, 1:3) = -stiff(1:3, 1:3)
667  stiff(1:3, 4:6) = transpose(stiff(4:6, 1:3))
668  stiff(4:6, 4:6) = stiff(1:3, 1:3)
669 
670  end subroutine
671 
672 end module m_static_lib_1d
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
hecmw
Definition: hecmw.f90:6
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.F90:138
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::getconnectorproperty
subroutine getconnectorproperty(gauss, ctype, dofid, rparams, iparams, stretch)
Definition: calMatMatrix.f90:339
m_static_lib_1d
This module provide common functions of 3D truss elements.
Definition: static_LIB_1d.f90:6
m_static_lib_1d::stf_dashpot_a
subroutine stf_dashpot_a(gausses, ecoord, u, stiff)
Definition: static_LIB_1d.f90:627
m_static_lib_1d::update_spring_d
subroutine update_spring_d(gausses, ecoord, u, du, qf)
Definition: static_LIB_1d.f90:480
m_static_lib_1d::update_spring_a
subroutine update_spring_a(gausses, ecoord, u, du, qf)
Definition: static_LIB_1d.f90:563
m_static_lib_1d::stf_spring_d
subroutine stf_spring_d(gausses, stiff)
Definition: static_LIB_1d.f90:459
m_static_lib_1d::elementstress_c1
subroutine elementstress_c1(ETYPE, gausses, strain, stress)
Definition: static_LIB_1d.f90:166
m_static_lib_1d::nodalstress_c1
subroutine nodalstress_c1(ETYPE, NN, gausses, ndstrain, ndstress)
Definition: static_LIB_1d.f90:148
m_static_lib_1d::search_diag_modify
subroutine search_diag_modify(n, nn, hecMAT, hecMESH)
Definition: static_LIB_1d.f90:282
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::stf_connector
subroutine stf_connector(etype, nn, ecoord, gausses, stiff, u, temperature)
This subroutine calculate stiff matrix of 2-nodes connector element.
Definition: static_LIB_1d.f90:363
m_static_lib_1d::stf_spring_a
subroutine stf_spring_a(gausses, ecoord, u, stiff)
Definition: static_LIB_1d.f90:513
m_static_lib_1d::dmp_connector
subroutine dmp_connector(etype, nn, ecoord, gausses, damping, u, temperature)
This subroutine calculate damping matrix of 2-nodes connector element.
Definition: static_LIB_1d.f90:390
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::stf_dashpot_d
subroutine stf_dashpot_d(gausses, stiff)
Definition: static_LIB_1d.f90:606
m_static_lib_1d::truss_diag_modify
subroutine truss_diag_modify(hecMAT, hecMESH)
Definition: static_LIB_1d.f90:245
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:184
m_static_lib_1d::update_connector
subroutine update_connector(etype, nn, ecoord, u, du, qf, gausses, TT)
Update strain and stress inside spring element.
Definition: static_LIB_1d.f90:418
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13