FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
utilities.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_utilities
7  use hecmw
8  implicit none
9 
10  real(kind=kreal), parameter, private :: pi=3.14159265358979d0
11 
12 contains
13 
15  subroutine memget(var,dimn,syze)
16  integer :: var,dimn,syze,bite
17  parameter(bite=1)
18  var = var + dimn*syze*bite
19  end subroutine memget
20 
21 
23  subroutine append_int2name( n, fname, n1 )
24  integer, intent(in) :: n
25  integer, intent(in), optional :: n1
26  character(len=*), intent(inout) :: fname
27  integer :: npos, nlen
28  character(len=128) :: tmpname, tmp
29 
30  npos = scan( fname, '.')
31  nlen = len_trim( fname )
32  if( nlen>128 ) stop "String too long(>128) in append_int2name"
33  if( n>100000 ) stop "Integer too big>100000 in append_int2name"
34  tmpname = fname
35  if( npos==0 ) then
36  write( fname, '(a,i6)') fname(1:nlen),n
37  else
38  write( tmp, '(i6,a)') n,tmpname(npos:nlen)
39  fname = tmpname(1:npos-1) // adjustl(tmp)
40  endif
41  if(present(n1).and.n1/=0)then
42  write(tmp,'(i8)')n1
43  fname = fname(1:len_trim(fname))//'.'//adjustl(tmp)
44  endif
45  end subroutine
46 
48  subroutine insert_int2array( iin, carray )
49  integer, intent(in) :: iin
50  integer, pointer :: carray(:)
51 
52  integer :: i, oldsize
53  integer, pointer :: dumarray(:) => null()
54  if( .not. associated(carray) ) then
55  allocate( carray(1) )
56  carray(1) = iin
57  else
58  oldsize = size( carray )
59  allocate( dumarray(oldsize) )
60  do i=1,oldsize
61  dumarray(i) = carray(i)
62  enddo
63  deallocate( carray )
64  allocate( carray(oldsize+1) )
65  do i=1,oldsize
66  carray(i) = dumarray(i)
67  enddo
68  carray(oldsize+1) = iin
69  endif
70  if( associated(dumarray) ) deallocate( dumarray )
71  end subroutine
72 
74  subroutine tensor_eigen3( tensor, eigval )
75  real(kind=kreal), intent(in) :: tensor(6)
76  real(kind=kreal), intent(out) :: eigval(3)
77 
78  real(kind=kreal) :: i1,i2,i3,r,sita,q, x(3,3), xx(3,3), ii(3,3)
79 
80  ii(:,:)=0.d0
81  ii(1,1)=1.d0; ii(2,2)=1.d0; ii(3,3)=1.d0
82  x(1,1)=tensor(1); x(2,2)=tensor(2); x(3,3)=tensor(3)
83  x(1,2)=tensor(4); x(2,1)=x(1,2)
84  x(2,3)=tensor(5); x(3,2)=x(2,3)
85  x(3,1)=tensor(6); x(1,3)=x(3,1)
86 
87  xx= matmul( x,x )
88  i1= x(1,1)+x(2,2)+x(3,3)
89  i2= 0.5d0*( i1*i1 - (xx(1,1)+xx(2,2)+xx(3,3)) )
90  i3= x(1,1)*x(2,2)*x(3,3)+x(2,1)*x(3,2)*x(1,3)+x(3,1)*x(1,2)*x(2,3) &
91  -x(3,1)*x(2,2)*x(1,3)-x(2,1)*x(1,2)*x(3,3)-x(1,1)*x(3,2)*x(2,3)
92 
93  r=(-2.d0*i1*i1*i1+9.d0*i1*i2-27.d0*i3)/54.d0
94  q=(i1*i1-3.d0*i2)/9.d0
95  sita = acos(r/dsqrt(q*q*q))
96 
97  eigval(1) = -2.d0*q*cos(sita/3.d0)+i1/3.d0
98  eigval(2) = -2.d0*q*cos((sita+2.d0*pi)/3.d0)+i1/3.d0
99  eigval(3) = -2.d0*q*cos((sita-2.d0*pi)/3.d0)+i1/3.d0
100 
101  end subroutine
102 
105  subroutine eigen3 (tensor, eigval, princ)
106  real(kind=kreal), intent(in) :: tensor(6)
107  real(kind=kreal), intent(out) :: eigval(3)
108  real(kind=kreal), intent(out) :: princ(3, 3)
109 
110  integer, parameter :: msweep = 50
111  integer :: i,j, is, ip, iq, ir
112  real(kind=kreal) :: fsum, od, theta, t, c, s, tau, g, h, hd, btens(3,3), factor
113 
114  btens(1,1)=tensor(1); btens(2,2)=tensor(2); btens(3,3)=tensor(3)
115  btens(1,2)=tensor(4); btens(2,1)=btens(1,2)
116  btens(2,3)=tensor(5); btens(3,2)=btens(2,3)
117  btens(3,1)=tensor(6); btens(1,3)=btens(3,1)
118  !
119  ! Initialise princ to the identity
120  !
121  factor = 0.d0
122  do i = 1, 3
123  do j = 1, 3
124  princ(i, j) = 0.d0
125  end do
126  princ(i, i) = 1.d0
127  eigval(i) = btens(i, i)
128  factor = factor + dabs(btens(i,i))
129  end do
130  ! Scaling and iszero/isnan exception
131  if( factor == 0.d0 .or. factor /= factor ) then
132  return
133  else
134  eigval(1:3) = eigval(1:3)/factor
135  btens(1:3,1:3) = btens(1:3,1:3)/factor
136  end if
137 
138  !
139  ! Starts sweeping.
140  !
141  do is = 1, msweep
142  fsum = 0.d0
143  do ip = 1, 2
144  do iq = ip + 1, 3
145  fsum = fsum + abs( btens(ip, iq) )
146  end do
147  end do
148  !
149  ! If the fsum of off-diagonal terms is zero returns
150  !
151  if ( fsum < 1.d-10 ) then
152  eigval(1:3) = eigval(1:3)*factor
153  return
154  endif
155  !
156  ! Performs the sweep in three rotations. One per off diagonal term
157  !
158  do ip = 1, 2
159  do iq = ip + 1, 3
160  od = 100.d0 * abs(btens(ip, iq) )
161  if ( (od+abs(eigval(ip) ) /= abs(eigval(ip) )) &
162  .and. (od+abs(eigval(iq) ) /= abs(eigval(iq) ))) then
163  hd = eigval(iq) - eigval(ip)
164  !
165  ! Evaluates the rotation angle
166  !
167  if ( abs(hd) + od == abs(hd) ) then
168  t = btens(ip, iq) / hd
169  else
170  theta = 0.5d0 * hd / btens(ip, iq)
171  t = 1.d0 / (abs(theta) + sqrt(1.d0 + theta**2) )
172  if ( theta < 0.d0 ) t = - t
173  end if
174  !
175  ! Re-evaluates the diagonal terms
176  !
177  c = 1.d0 / sqrt(1.d0 + t**2)
178  s = t * c
179  tau = s / (1.d0 + c)
180  h = t * btens(ip, iq)
181  eigval(ip) = eigval(ip) - h
182  eigval(iq) = eigval(iq) + h
183  !
184  ! Re-evaluates the remaining off-diagonal terms
185  !
186  ir = 6 - ip - iq
187  g = btens(min(ir, ip), max(ir, ip) )
188  h = btens(min(ir, iq), max(ir, iq) )
189  btens(min(ir, ip), max(ir, ip) ) = g &
190  - s * (h + g * tau)
191  btens(min(ir, iq), max(ir, iq) ) = h &
192  + s * (g - h * tau)
193  !
194  ! Rotates the eigenvectors
195  !
196  do ir = 1, 3
197  g = princ(ir, ip)
198  h = princ(ir, iq)
199  princ(ir, ip) = g - s * (h + g * tau)
200  princ(ir, iq) = h + s * (g - h * tau)
201  end do
202  end if
203  btens(ip, iq) = 0.d0
204  end do
205  end do
206  end do ! over sweeps
207  !
208  ! If convergence is not achieved stops
209  !
210  stop ' Jacobi iteration unable to converge'
211  end subroutine eigen3
212 
214  real(kind=kreal) function determinant( mat )
215  real(kind=kreal) :: mat(6)
216  real(kind=kreal) :: xj(3,3)
217 
218  xj(1,1)=mat(1); xj(2,2)=mat(2); xj(3,3)=mat(3)
219  xj(1,2)=mat(4); xj(2,1)=xj(1,2)
220  xj(2,3)=mat(5); xj(3,2)=xj(2,3)
221  xj(3,1)=mat(6); xj(1,3)=xj(3,1)
222 
223  determinant=xj(1,1)*xj(2,2)*xj(3,3) &
224  +xj(2,1)*xj(3,2)*xj(1,3) &
225  +xj(3,1)*xj(1,2)*xj(2,3) &
226  -xj(3,1)*xj(2,2)*xj(1,3) &
227  -xj(2,1)*xj(1,2)*xj(3,3) &
228  -xj(1,1)*xj(3,2)*xj(2,3)
229  end function determinant
230 
232  real(kind=kreal) function determinant33( XJ )
233  real(kind=kreal) :: xj(3,3)
234 
235  determinant33=xj(1,1)*xj(2,2)*xj(3,3) &
236  +xj(2,1)*xj(3,2)*xj(1,3) &
237  +xj(3,1)*xj(1,2)*xj(2,3) &
238  -xj(3,1)*xj(2,2)*xj(1,3) &
239  -xj(2,1)*xj(1,2)*xj(3,3) &
240  -xj(1,1)*xj(3,2)*xj(2,3)
241  end function determinant33
242 
243  subroutine fstr_chk_alloc( imsg, sub_name, ierr )
244  use hecmw
245  character(*) :: sub_name
246  integer(kind=kint) :: imsg
247  integer(kind=kint) :: ierr
248 
249  if( ierr /= 0 ) then
250  write(imsg,*) 'Memory overflow at ', sub_name
251  write(*,*) 'Memory overflow at ', sub_name
252  call hecmw_abort( hecmw_comm_get_comm( ) )
253  endif
254  end subroutine fstr_chk_alloc
255 
257  subroutine calinverse(NN, A)
258  integer, intent(in) :: NN
259  real(kind=kreal), intent(inout) :: a(nn,nn)
260 
261  integer :: I, J,K,IW,LR,IP(NN)
262  real(kind=kreal) :: w,wmax,pivot,api,eps,det
263  data eps/1.0e-35/
264  det=1.d0
265  lr = 0.0d0
266  do i=1,nn
267  ip(i)=i
268  enddo
269  do k=1,nn
270  wmax=0.d0
271  do i=k,nn
272  w=dabs(a(i,k))
273  if (w.GT.wmax) then
274  wmax=w
275  lr=i
276  endif
277  enddo
278  pivot=a(lr,k)
279  api=abs(pivot)
280  if(api.LE.eps) then
281  write(*,'(''PIVOT ERROR AT'',I5)') k
282  stop
283  end if
284  det=det*pivot
285  if (lr.NE.k) then
286  det=-det
287  iw=ip(k)
288  ip(k)=ip(lr)
289  ip(lr)=iw
290  do j=1,nn
291  w=a(k,j)
292  a(k,j)=a(lr,j)
293  a(lr,j)=w
294  enddo
295  endif
296  do i=1,nn
297  a(k,i)=a(k,i)/pivot
298  enddo
299  do i=1,nn
300  if (i.NE.k) then
301  w=a(i,k)
302  if (w.NE.0.) then
303  do j=1,nn
304  if (j.NE.k) a(i,j)=a(i,j)-w*a(k,j)
305  enddo
306  a(i,k)=-w/pivot
307  endif
308  endif
309  enddo
310  a(k,k)=1.d0/pivot
311  enddo
312 
313  do i=1,nn
314  k=ip(i)
315  if (k.NE.i) then
316  iw=ip(k)
317  ip(k)=ip(i)
318  ip(i)=iw
319  do j=1,nn
320  w=a(j,i)
321  a(j,i)=a(j,k)
322  a(j,k)=w
323  enddo
324  endif
325  enddo
326 
327  end subroutine calinverse
328 
329  subroutine cross_product(v1,v2,vn)
330  real(kind=kreal),intent(in) :: v1(3),v2(3)
331  real(kind=kreal),intent(out) :: vn(3)
332 
333  vn(1) = v1(2)*v2(3) - v1(3)*v2(2)
334  vn(2) = v1(3)*v2(1) - v1(1)*v2(3)
335  vn(3) = v1(1)*v2(2) - v1(2)*v2(1)
336  end subroutine cross_product
337 
338  subroutine transformation(jacob, tm)
339  real(kind=kreal),intent(in) :: jacob(3,3)
340  real(kind=kreal),intent(out) :: tm(6,6)
341 
342  integer :: i,j
343 
344  do i=1,3
345  do j=1,3
346  tm(i,j)= jacob(i,j)*jacob(i,j)
347  enddo
348  tm(i,4) = jacob(i,1)*jacob(i,2)
349  tm(i,5) = jacob(i,2)*jacob(i,3)
350  tm(i,6) = jacob(i,3)*jacob(i,1)
351  enddo
352  tm(4,1) = 2.d0*jacob(1,1)*jacob(2,1)
353  tm(5,1) = 2.d0*jacob(2,1)*jacob(3,1)
354  tm(6,1) = 2.d0*jacob(3,1)*jacob(1,1)
355  tm(4,2) = 2.d0*jacob(1,2)*jacob(2,2)
356  tm(5,2) = 2.d0*jacob(2,2)*jacob(3,2)
357  tm(6,2) = 2.d0*jacob(3,2)*jacob(1,2)
358  tm(4,3) = 2.d0*jacob(1,3)*jacob(2,3)
359  tm(5,3) = 2.d0*jacob(2,3)*jacob(3,3)
360  tm(6,3) = 2.d0*jacob(3,3)*jacob(1,3)
361  tm(4,4) = jacob(1,1)*jacob(2,2) + jacob(1,2)*jacob(2,1)
362  tm(5,4) = jacob(2,1)*jacob(3,2) + jacob(2,2)*jacob(3,1)
363  tm(6,4) = jacob(3,1)*jacob(1,2) + jacob(3,2)*jacob(1,1)
364  tm(4,5) = jacob(1,2)*jacob(2,3) + jacob(1,3)*jacob(2,2)
365  tm(5,5) = jacob(2,2)*jacob(3,3) + jacob(2,3)*jacob(3,2)
366  tm(6,5) = jacob(3,2)*jacob(1,3) + jacob(3,3)*jacob(1,2)
367  tm(4,6) = jacob(1,3)*jacob(2,1) + jacob(1,1)*jacob(2,3)
368  tm(5,6) = jacob(2,3)*jacob(3,1) + jacob(2,1)*jacob(3,3)
369  tm(6,6) = jacob(3,3)*jacob(1,1) + jacob(3,1)*jacob(1,3)
370 
371  end subroutine transformation
372 
373  subroutine get_principal (tensor, eigval, princmatrix)
374 
375  implicit none
376  integer i,j
377  real(kind=kreal) :: tensor(1:6)
378  real(kind=kreal) :: eigval(3)
379  real(kind=kreal) :: princmatrix(3,3)
380  real(kind=kreal) :: princnormal(3,3)
381  real(kind=kreal) :: tempv(3)
382  real(kind=kreal) :: temps
383 
384  call eigen3(tensor,eigval,princnormal)
385 
386  if (eigval(1)<eigval(2)) then
387  temps=eigval(1)
388  eigval(1)=eigval(2)
389  eigval(2)=temps
390  tempv(:)=princnormal(:,1)
391  princnormal(:,1)=princnormal(:,2)
392  princnormal(:,2)=tempv(:)
393  end if
394  if (eigval(1)<eigval(3)) then
395  temps=eigval(1)
396  eigval(1)=eigval(3)
397  eigval(3)=temps
398  tempv(:)=princnormal(:,1)
399  princnormal(:,1)=princnormal(:,3)
400  princnormal(:,3)=tempv(:)
401  end if
402  if (eigval(2)<eigval(3)) then
403  temps=eigval(2)
404  eigval(2)=eigval(3)
405  eigval(3)=temps
406  tempv(:)=princnormal(:,2)
407  princnormal(:,2)=princnormal(:,3)
408  princnormal(:,3)=tempv(:)
409  end if
410 
411  do j=1,3
412  do i=1,3
413  princmatrix(i,j) = princnormal(i,j) * eigval(j)
414  end do
415  end do
416 
417  end subroutine get_principal
418 
419  subroutine eigen3d (tensor, eigval, princ)
420  implicit none
421 
422  real(kind=kreal) :: tensor(6)
423  real(kind=kreal) :: eigval(3)
424  real(kind=kreal) :: princ(3,3)
425 
426  real(kind=kreal) :: s11, s22, s33, s12, s23, s13, j1, j2, j3
427  real(kind=kreal) :: ml,nl
428  complex(kind=kreal):: x1,x2,x3
429  real(kind=kreal):: rtemp
430  integer :: i
431  s11 = tensor(1)
432  s22 = tensor(2)
433  s33 = tensor(3)
434  s12 = tensor(4)
435  s23 = tensor(5)
436  s13 = tensor(6)
437 
438  ! invariants of stress tensor
439  j1 = s11 + s22 + s33
440  j2 = -s11*s22 - s22*s33 - s33*s11 + s12**2 + s23**2 + s13**2
441  j3 = s11*s22*s33 + 2*s12*s23*s13 - s11*s23**2 - s22*s13**2 - s33*s12**2
442  ! Cardano's method
443  ! x^3+ ax^2 + bx +c =0
444  ! s^3 - J1*s^2 -J2s -J3 =0
445  call cardano(-j1, -j2, -j3, x1, x2, x3)
446  eigval(1)= real(x1)
447  eigval(2)= real(x2)
448  eigval(3)= real(x3)
449  if (eigval(1)<eigval(2)) then
450  rtemp=eigval(1)
451  eigval(1)=eigval(2)
452  eigval(2)=rtemp
453  end if
454  if (eigval(1)<eigval(3)) then
455  rtemp=eigval(1)
456  eigval(1)=eigval(3)
457  eigval(3)=rtemp
458  end if
459  if (eigval(2)<eigval(3)) then
460  rtemp=eigval(2)
461  eigval(2)=eigval(3)
462  eigval(3)=rtemp
463  end if
464 
465  do i=1,3
466  if (eigval(i)/(eigval(1)+eigval(2)+eigval(3)) < 1.0d-10 )then
467  eigval(i) = 0.0d0
468  princ(i,:) = 0.0d0
469  exit
470  end if
471  ml = ( s23*s13 - s12*(s33-eigval(i)) ) / ( -s23**2 + (s22-eigval(i))*(s33-eigval(i)) )
472  nl = ( s12**2 - (s22-eigval(i))*(s11-eigval(i)) ) / ( s12*s23 - s13*(s22-eigval(i)) )
473  if (abs(ml) >= huge(ml)) then
474  ml=0.0d0
475  end if
476  if (abs(nl) >= huge(nl)) then
477  nl=0.0d0
478  end if
479  princ(i,1) = eigval(i)/sqrt( 1 + ml**2 + nl**2)
480  princ(i,2) = ml * princ(i,1)
481  princ(i,3) = nl * princ(i,1)
482  end do
483 
484  write(*,*)
485  end subroutine eigen3d
486 
487  subroutine cardano(a,b,c,x1,x2,x3)
488  real(kind=kreal):: a,b,c
489  real(kind=kreal):: p,q,d
490  complex(kind=kreal):: w
491  complex(kind=kreal):: u,v,y
492  complex(kind=kreal):: x1,x2,x3
493  w = (-1.0d0 + sqrt(dcmplx(-3.0d0)))/2.0d0
494  p = -a**2/9.0d0 + b/3.0d0
495  q = 2.0d0/2.7d1*a**3 - a*b/3.0d0 + c
496  d = q**2 + 4.0d0*p**3
497 
498  u = ((-dcmplx(q) + sqrt(dcmplx(d)))/2.0d0)**(1.0d0/3.0d0)
499 
500  if(u.ne.0.0d0) then
501  v = -dcmplx(p)/u
502  x1 = u + v -dcmplx(a)/3.0d0
503  x2 = u*w + v*w**2 -dcmplx(a)/3.0d0
504  x3 = u*w**2 + v*w -dcmplx(a)/3.0d0
505  else
506  y = (-dcmplx(q))**(1.0d0/3.0d0)
507  x1 = y -dcmplx(a)/3.0d0
508  x2 = y*w -dcmplx(a)/3.0d0
509  x3 = y*w**2 -dcmplx(a)/3.0d0
510  end if
511 
512  end subroutine cardano
513 
515  real(kind=kreal),intent(in) :: r(3)
516  real(kind=kreal),intent(inout) :: v(3)
517 
518  real(kind=kreal) :: rotv(3), rv
519  real(kind=kreal) :: cosx, sinc(2)
520  real(kind=kreal) :: x, x2, x4, x6
521  real(kind=kreal), parameter :: c0 = 0.5d0
522  real(kind=kreal), parameter :: c2 = -4.166666666666666d-002
523  real(kind=kreal), parameter :: c4 = 1.388888888888889d-003
524  real(kind=kreal), parameter :: c6 = -2.480158730158730d-005
525 
526  x2 = dot_product(r,r)
527  x = dsqrt(x2)
528  cosx = dcos(x)
529  if( x < 1.d-4 ) then
530  x4 = x2*x2
531  x6 = x4*x2
532  sinc(1) = 1.d0-x2/6.d0+x4/120.d0
533  sinc(2) = c0+c2*x2+c4*x4+c6*x6
534  else
535  sinc(1) = dsin(x)/x
536  sinc(2) = (1.d0-cosx)/x2
537  endif
538 
539  ! calc Rot*v
540  rv = dot_product(r,v)
541  rotv(1:3) = cosx*v(1:3)
542  rotv(1:3) = rotv(1:3)+rv*sinc(2)*r(1:3)
543  rotv(1) = rotv(1) + (-v(2)*r(3)+v(3)*r(2))*sinc(1)
544  rotv(2) = rotv(2) + (-v(3)*r(1)+v(1)*r(3))*sinc(1)
545  rotv(3) = rotv(3) + (-v(1)*r(2)+v(2)*r(1))*sinc(1)
546  v = rotv
547 
548  end subroutine
549 
550 end module
Definition: hecmw.f90:6
This module provides aux functions.
Definition: utilities.f90:6
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
Definition: utilities.f90:106
subroutine cross_product(v1, v2, vn)
Definition: utilities.f90:330
subroutine eigen3d(tensor, eigval, princ)
Definition: utilities.f90:420
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
subroutine insert_int2array(iin, carray)
Insert an integer into a integer array.
Definition: utilities.f90:49
real(kind=kreal) function determinant(mat)
Compute determinant for symmetric 3*3 matrix.
Definition: utilities.f90:215
subroutine append_int2name(n, fname, n1)
Insert an integer at end of a file name.
Definition: utilities.f90:24
subroutine tensor_eigen3(tensor, eigval)
Given symmetric 3x3 matrix M, compute the eigenvalues.
Definition: utilities.f90:75
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:374
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
subroutine fstr_chk_alloc(imsg, sub_name, ierr)
Definition: utilities.f90:244
subroutine cardano(a, b, c, x1, x2, x3)
Definition: utilities.f90:488
subroutine memget(var, dimn, syze)
Record used memeory.
Definition: utilities.f90:16
subroutine rotate_3dvector_by_rodrigues_formula(r, v)
Definition: utilities.f90:515
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:258