FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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!-------------------------------------------------------------------------------
7 use hecmw
8 implicit none
9
10 real(kind=kreal), parameter, private :: pi=3.14159265358979d0
11
12contains
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
550end module
Definition: hecmw.f90:6
This module provides aux functions.
Definition: utilities.f90:6
subroutine calinverse(nn, a)
calculate inverse of matrix a
Definition: utilities.f90:258
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