DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
math.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/math.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/math.f90"
5 !--------------------------------------------------------------------------------------------------
11 !--------------------------------------------------------------------------------------------------
12 module math
13  use prec
14  use io
15  use numerics
16 
17  implicit none
18  public
19 
20 
21 
22 
23 
24 
25 
26 
27  real(preal), parameter :: pi = acos(-1.0_preal)
28  real(preal), parameter :: indeg = 180.0_preal/pi
29  real(preal), parameter :: inrad = pi/180.0_preal
30  complex(pReal), parameter :: twopiimg = cmplx(0.0_preal,2.0_preal*pi)
31 
32  real(preal), dimension(3,3), parameter :: &
33  math_i3 = reshape([&
34  1.0_preal,0.0_preal,0.0_preal, &
35  0.0_preal,1.0_preal,0.0_preal, &
36  0.0_preal,0.0_preal,1.0_preal &
37  ],[3,3])
38 
39  real(preal), dimension(6), parameter, private :: &
40  nrmmandel = [&
41  1.0_preal, 1.0_preal, 1.0_preal, &
42  sqrt(2.0_preal), sqrt(2.0_preal), sqrt(2.0_preal) ]
43 
44  real(preal), dimension(6), parameter, private :: &
45  invnrmmandel = 1.0_preal/nrmmandel
46 
47  integer, dimension (2,6), parameter, private :: &
48  mapnye = reshape([&
49  1,1, &
50  2,2, &
51  3,3, &
52  1,2, &
53  2,3, &
54  1,3 &
55  ],[2,6])
56 
57  integer, dimension (2,6), parameter, private :: &
58  mapvoigt = reshape([&
59  1,1, &
60  2,2, &
61  3,3, &
62  2,3, &
63  1,3, &
64  1,2 &
65  ],[2,6])
66 
67  integer, dimension (2,9), parameter, private :: &
68  mapplain = reshape([&
69  1,1, &
70  1,2, &
71  1,3, &
72  2,1, &
73  2,2, &
74  2,3, &
75  3,1, &
76  3,2, &
77  3,3 &
78  ],[2,9])
79 
80  interface math_eye
81  module procedure math_identity2nd
82  end interface math_eye
83 
84 
85 !---------------------------------------------------------------------------------------------------
86  private :: &
87  unittest
88 
89 contains
90 
91 !--------------------------------------------------------------------------------------------------
93 !--------------------------------------------------------------------------------------------------
94 subroutine math_init
95 
96  real(pReal), dimension(4) :: randTest
97  integer :: randSize
98  integer, dimension(:), allocatable :: randInit
99 
100  write(6,'(/,a)') ' <<<+- math init -+>>>'; flush(6)
101 
102  call random_seed(size=randsize)
103  allocate(randinit(randsize))
104  if (randomseed > 0) then
105  randinit = randomseed
106  else
107  call random_seed()
108  call random_seed(get = randinit)
109  randinit(2:randsize) = randinit(1)
110  endif
111 
112  call random_seed(put = randinit)
113  call random_number(randtest)
114 
115  write(6,'(a,i2)') ' size of random seed: ', randsize
116  write(6,'(a,i0)') ' value of random seed: ', randinit(1)
117  write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randtest
118 
119  call random_seed(put = randinit)
120 
121  call unittest
122 
123 end subroutine math_init
124 
125 
126 !--------------------------------------------------------------------------------------------------
128 ! Sorting is done with respect to array(sort,:) and keeps array(/=sort,:) linked to it.
129 ! default: sort=1
130 !--------------------------------------------------------------------------------------------------
131 recursive subroutine math_sort(a, istart, iend, sortDim)
132 
133  integer, dimension(:,:), intent(inout) :: a
134  integer, intent(in),optional :: istart,iend, sortdim
135  integer :: ipivot,s,e,d
136 
137  if(present(istart)) then
138  s = istart
139  else
140  s = lbound(a,2)
141  endif
142 
143  if(present(iend)) then
144  e = iend
145  else
146  e = ubound(a,2)
147  endif
148 
149  if(present(sortdim)) then
150  d = sortdim
151  else
152  d = 1
153  endif
154 
155  if (s < e) then
156  ipivot = qsort_partition(a,s, e, d)
157  call math_sort(a, s, ipivot-1, d)
158  call math_sort(a, ipivot+1, e, d)
159  endif
160 
161 
162  contains
163 
164  !-------------------------------------------------------------------------------------------------
166  !-------------------------------------------------------------------------------------------------
167  integer function qsort_partition(a, istart, iend, sort)
168 
169  integer, dimension(:,:), intent(inout) :: a
170  integer, intent(in) :: istart,iend,sort
171  integer, dimension(size(a,1)) :: tmp
172  integer :: i,j
173 
174  do
175  ! find the first element on the right side less than or equal to the pivot point
176  do j = iend, istart, -1
177  if (a(sort,j) <= a(sort,istart)) exit
178  enddo
179  ! find the first element on the left side greater than the pivot point
180  do i = istart, iend
181  if (a(sort,i) > a(sort,istart)) exit
182  enddo
183  cross: if (i >= j) then ! exchange left value with pivot and return with the partition index
184  tmp = a(:,istart)
185  a(:,istart) = a(:,j)
186  a(:,j) = tmp
187  qsort_partition = j
188  return
189  else cross ! exchange values
190  tmp = a(:,i)
191  a(:,i) = a(:,j)
192  a(:,j) = tmp
193  endif cross
194  enddo
195 
196  end function qsort_partition
197 
198 end subroutine math_sort
199 
200 
201 !--------------------------------------------------------------------------------------------------
205 !--------------------------------------------------------------------------------------------------
206 pure function math_expand(what,how)
207 
208  real(preal), dimension(:), intent(in) :: what
209  integer, dimension(:), intent(in) :: how
210  real(preal), dimension(sum(how)) :: math_expand
211  integer :: i
212 
213  if (sum(how) == 0) return
214 
215  do i = 1, size(how)
216  math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
217  enddo
218 
219 end function math_expand
220 
221 
222 !--------------------------------------------------------------------------------------------------
224 !--------------------------------------------------------------------------------------------------
225 pure function math_range(N)
226 
227  integer, intent(in) :: n
228  integer :: i
229  integer, dimension(N) :: math_range
230 
231  math_range = [(i,i=1,n)]
232 
233 end function math_range
234 
235 
236 !--------------------------------------------------------------------------------------------------
238 !--------------------------------------------------------------------------------------------------
239 pure function math_identity2nd(d)
240 
241  integer, intent(in) :: d
242  integer :: i
243  real(preal), dimension(d,d) :: math_identity2nd
244 
245  math_identity2nd = 0.0_preal
246  do i=1,d
247  math_identity2nd(i,i) = 1.0_preal
248  enddo
249 
250 end function math_identity2nd
251 
252 
253 !--------------------------------------------------------------------------------------------------
255 ! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
256 !--------------------------------------------------------------------------------------------------
257 pure function math_identity4th(d)
258 
259  integer, intent(in) :: d
260  integer :: i,j,k,l
261  real(preal), dimension(d,d,d,d) :: math_identity4th
262  real(preal), dimension(d,d) :: identity2nd
263 
264  identity2nd = math_identity2nd(d)
265  do i=1,d; do j=1,d; do k=1,d; do l=1,d
266  math_identity4th(i,j,k,l) = 0.5_preal &
267  *(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
268  enddo; enddo; enddo; enddo
269 
270 end function math_identity4th
271 
272 
273 !--------------------------------------------------------------------------------------------------
275 ! e_ijk = 1 if even permutation of ijk
276 ! e_ijk = -1 if odd permutation of ijk
277 ! e_ijk = 0 otherwise
278 !--------------------------------------------------------------------------------------------------
279 real(preal) pure function math_levicivita(i,j,k)
280 
281  integer, intent(in) :: i,j,k
282 
283  if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
284  math_levicivita = +1.0_preal
285  elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
286  math_levicivita = -1.0_preal
287  else
288  math_levicivita = 0.0_preal
289  endif
290 
291 end function math_levicivita
292 
293 
294 !--------------------------------------------------------------------------------------------------
296 ! d_ij = 1 if i = j
297 ! d_ij = 0 otherwise
298 !--------------------------------------------------------------------------------------------------
299 real(preal) pure function math_delta(i,j)
300 
301  integer, intent (in) :: i,j
302 
303  math_delta = merge(0.0_preal, 1.0_preal, i /= j)
304 
305 end function math_delta
306 
307 
308 !--------------------------------------------------------------------------------------------------
310 !--------------------------------------------------------------------------------------------------
311 pure function math_cross(A,B)
312 
313  real(preal), dimension(3), intent(in) :: a,b
314  real(preal), dimension(3) :: math_cross
315 
316  math_cross = [ a(2)*b(3) -a(3)*b(2), &
317  a(3)*b(1) -a(1)*b(3), &
318  a(1)*b(2) -a(2)*b(1) ]
319 
320 end function math_cross
321 
322 
323 !--------------------------------------------------------------------------------------------------
325 !--------------------------------------------------------------------------------------------------
326 pure function math_outer(A,B)
327 
328  real(preal), dimension(:), intent(in) :: a,b
329  real(preal), dimension(size(A,1),size(B,1)) :: math_outer
330  integer :: i,j
331 
332  do i=1,size(a,1); do j=1,size(b,1)
333  math_outer(i,j) = a(i)*b(j)
334  enddo; enddo
335 
336 end function math_outer
337 
338 
339 !--------------------------------------------------------------------------------------------------
341 !--------------------------------------------------------------------------------------------------
342 real(preal) pure function math_inner(a,b)
343 
344  real(preal), dimension(:), intent(in) :: a
345  real(preal), dimension(size(A,1)), intent(in) :: b
346 
347  math_inner = sum(a*b)
348 
349 end function math_inner
350 
351 
352 !--------------------------------------------------------------------------------------------------
354 !--------------------------------------------------------------------------------------------------
355 real(preal) pure function math_tensordot(a,b)
356 
357  real(preal), dimension(3,3), intent(in) :: a,b
358 
359  math_tensordot = sum(a*b)
360 
361 end function math_tensordot
362 
363 
364 !--------------------------------------------------------------------------------------------------
366 !--------------------------------------------------------------------------------------------------
367 pure function math_mul3333xx33(A,B)
368 
369  real(preal), dimension(3,3,3,3), intent(in) :: a
370  real(preal), dimension(3,3), intent(in) :: b
371  real(preal), dimension(3,3) :: math_mul3333xx33
372  integer :: i,j
373 
374  do i=1,3; do j=1,3
375  math_mul3333xx33(i,j) = sum(a(i,j,1:3,1:3)*b(1:3,1:3))
376  enddo; enddo
377 
378 end function math_mul3333xx33
379 
380 
381 !--------------------------------------------------------------------------------------------------
383 !--------------------------------------------------------------------------------------------------
384 pure function math_mul3333xx3333(A,B)
385 
386  integer :: i,j,k,l
387  real(preal), dimension(3,3,3,3), intent(in) :: a
388  real(preal), dimension(3,3,3,3), intent(in) :: b
389  real(preal), dimension(3,3,3,3) :: math_mul3333xx3333
390 
391  do i=1,3; do j=1,3; do k=1,3; do l=1,3
392  math_mul3333xx3333(i,j,k,l) = sum(a(i,j,1:3,1:3)*b(1:3,1:3,k,l))
393  enddo; enddo; enddo; enddo
394 
395 end function math_mul3333xx3333
396 
397 
398 !--------------------------------------------------------------------------------------------------
400 !--------------------------------------------------------------------------------------------------
401 pure function math_exp33(A,n)
402 
403  real(preal), dimension(3,3), intent(in) :: a
404  integer, intent(in), optional :: n
405  real(preal), dimension(3,3) :: b, math_exp33
406 
407  real(preal) :: invfac
408  integer :: n_,i
409 
410  if (present(n)) then
411  n_ = n
412  else
413  n_ = 5
414  endif
415 
416  invfac = 1.0_preal ! 0!
417  b = math_i3
418  math_exp33 = math_i3 ! A^0 = I
419 
420  do i = 1, n_
421  invfac = invfac/real(i,preal) ! invfac = 1/(i!)
422  b = matmul(b,a)
423  math_exp33 = math_exp33 + invfac*b ! exp = SUM (A^i)/(i!)
424  enddo
425 
426 end function math_exp33
427 
428 
429 !--------------------------------------------------------------------------------------------------
432 ! if determinant is close to zero
433 !--------------------------------------------------------------------------------------------------
434 pure function math_inv33(A)
435 
436  real(preal), dimension(3,3), intent(in) :: a
437  real(preal), dimension(3,3) :: math_inv33
438 
439  real(preal) :: deta
440  logical :: error
441 
442  call math_invert33(math_inv33,deta,error,a)
443  if(error) math_inv33 = 0.0_preal
444 
445 end function math_inv33
446 
447 
448 !--------------------------------------------------------------------------------------------------
451 ! Returns an error if not possible, i.e. if determinant is close to zero
452 !--------------------------------------------------------------------------------------------------
453 pure subroutine math_invert33(InvA, DetA, error, A)
454 
455  real(preal), dimension(3,3), intent(out) :: inva
456  real(preal), intent(out) :: deta
457  logical, intent(out) :: error
458  real(preal), dimension(3,3), intent(in) :: a
459 
460  inva(1,1) = a(2,2) * a(3,3) - a(2,3) * a(3,2)
461  inva(2,1) = -a(2,1) * a(3,3) + a(2,3) * a(3,1)
462  inva(3,1) = a(2,1) * a(3,2) - a(2,2) * a(3,1)
463 
464  deta = a(1,1) * inva(1,1) + a(1,2) * inva(2,1) + a(1,3) * inva(3,1)
465 
466  if (deq0(deta)) then
467  inva = 0.0_preal
468  error = .true.
469  else
470  inva(1,2) = -a(1,2) * a(3,3) + a(1,3) * a(3,2)
471  inva(2,2) = a(1,1) * a(3,3) - a(1,3) * a(3,1)
472  inva(3,2) = -a(1,1) * a(3,2) + a(1,2) * a(3,1)
473 
474  inva(1,3) = a(1,2) * a(2,3) - a(1,3) * a(2,2)
475  inva(2,3) = -a(1,1) * a(2,3) + a(1,3) * a(2,1)
476  inva(3,3) = a(1,1) * a(2,2) - a(1,2) * a(2,1)
477 
478  inva = inva/deta
479  error = .false.
480  endif
481 
482 end subroutine math_invert33
483 
484 
485 !--------------------------------------------------------------------------------------------------
487 !--------------------------------------------------------------------------------------------------
488 function math_invsym3333(A)
489 
490  real(preal),dimension(3,3,3,3) :: math_invsym3333
491 
492  real(preal),dimension(3,3,3,3),intent(in) :: a
493 
494  integer :: ierr
495  integer, dimension(6) :: ipiv6
496  real(preal), dimension(6,6) :: temp66
497  real(preal), dimension(6*(64+2)) :: work
498  logical :: error
499  external :: &
500  dgetrf, &
501  dgetri
502 
503  temp66 = math_sym3333to66(a)
504  call dgetrf(6,6,temp66,6,ipiv6,ierr)
505  error = (ierr /= 0)
506  call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr)
507  error = error .or. (ierr /= 0)
508  if (error) then
509  call io_error(400, ext_msg = 'math_invSym3333')
510  else
511  math_invsym3333 = math_66tosym3333(temp66)
512  endif
513 
514 end function math_invsym3333
515 
516 
517 !--------------------------------------------------------------------------------------------------
519 !--------------------------------------------------------------------------------------------------
520 subroutine math_invert(InvA, error, A)
521 
522  real(pReal), dimension(:,:), intent(in) :: A
523  real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
524  logical, intent(out) :: error
525 
526  integer, dimension(size(A,1)) :: ipiv
527  real(pReal), dimension(size(A,1)*(64+2)) :: work
528  integer :: ierr
529  external :: &
530  dgetrf, &
531  dgetri
532 
533  inva = a
534  call dgetrf(size(a,1),size(a,1),inva,size(a,1),ipiv,ierr)
535  error = (ierr /= 0)
536  call dgetri(size(a,1),inva,size(a,1),ipiv,work,size(work,1),ierr)
537  error = error .or. (ierr /= 0)
538 
539 end subroutine math_invert
540 
541 
542 !--------------------------------------------------------------------------------------------------
544 !--------------------------------------------------------------------------------------------------
545 pure function math_symmetric33(m)
546 
547  real(preal), dimension(3,3) :: math_symmetric33
548  real(preal), dimension(3,3), intent(in) :: m
549 
550  math_symmetric33 = 0.5_preal * (m + transpose(m))
551 
552 end function math_symmetric33
553 
554 
555 !--------------------------------------------------------------------------------------------------
557 !--------------------------------------------------------------------------------------------------
558 pure function math_symmetric66(m)
559 
560  real(preal), dimension(6,6) :: math_symmetric66
561  real(preal), dimension(6,6), intent(in) :: m
562 
563  math_symmetric66 = 0.5_preal * (m + transpose(m))
564 
565 end function math_symmetric66
566 
567 
568 !--------------------------------------------------------------------------------------------------
570 !--------------------------------------------------------------------------------------------------
571 pure function math_skew33(m)
572 
573  real(preal), dimension(3,3) :: math_skew33
574  real(preal), dimension(3,3), intent(in) :: m
575 
576  math_skew33 = m - math_symmetric33(m)
577 
578 end function math_skew33
579 
580 
581 !--------------------------------------------------------------------------------------------------
583 !--------------------------------------------------------------------------------------------------
584 pure function math_spherical33(m)
585 
586  real(preal), dimension(3,3) :: math_spherical33
587  real(preal), dimension(3,3), intent(in) :: m
588 
589  math_spherical33 = math_i3 * math_trace33(m)/3.0_preal
590 
591 end function math_spherical33
592 
593 
594 !--------------------------------------------------------------------------------------------------
596 !--------------------------------------------------------------------------------------------------
597 pure function math_deviatoric33(m)
598 
599  real(preal), dimension(3,3) :: math_deviatoric33
600  real(preal), dimension(3,3), intent(in) :: m
601 
602  math_deviatoric33 = m - math_spherical33(m)
603 
604 end function math_deviatoric33
605 
606 
607 !--------------------------------------------------------------------------------------------------
609 !--------------------------------------------------------------------------------------------------
610 real(preal) pure function math_trace33(m)
611 
612  real(preal), dimension(3,3), intent(in) :: m
613 
614  math_trace33 = m(1,1) + m(2,2) + m(3,3)
615 
616 end function math_trace33
617 
618 
619 !--------------------------------------------------------------------------------------------------
621 !--------------------------------------------------------------------------------------------------
622 real(preal) pure function math_det33(m)
623 
624  real(preal), dimension(3,3), intent(in) :: m
625 
626  math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
627  - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
628  + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1))
629 
630 end function math_det33
631 
632 
633 !--------------------------------------------------------------------------------------------------
635 !--------------------------------------------------------------------------------------------------
636 real(preal) pure function math_detsym33(m)
637 
638  real(preal), dimension(3,3), intent(in) :: m
639 
640  math_detsym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) &
641  + m(1,1)*m(2,2)*m(3,3) + 2.0_preal * m(1,2)*m(1,3)*m(2,3)
642 
643 end function math_detsym33
644 
645 
646 !--------------------------------------------------------------------------------------------------
648 !--------------------------------------------------------------------------------------------------
649 pure function math_33to9(m33)
650 
651  real(preal), dimension(9) :: math_33to9
652  real(preal), dimension(3,3), intent(in) :: m33
653 
654  integer :: i
655 
656  do i = 1, 9
657  math_33to9(i) = m33(mapplain(1,i),mapplain(2,i))
658  enddo
659 
660 end function math_33to9
661 
662 
663 !--------------------------------------------------------------------------------------------------
665 !--------------------------------------------------------------------------------------------------
666 pure function math_9to33(v9)
667 
668  real(preal), dimension(3,3) :: math_9to33
669  real(preal), dimension(9), intent(in) :: v9
670 
671  integer :: i
672 
673  do i = 1, 9
674  math_9to33(mapplain(1,i),mapplain(2,i)) = v9(i)
675  enddo
676 
677 end function math_9to33
678 
679 
680 !--------------------------------------------------------------------------------------------------
683 ! components according to Mandel. Advisable for matrix operations.
684 ! Unweighted conversion only changes order according to Nye
685 !--------------------------------------------------------------------------------------------------
686 pure function math_sym33to6(m33,weighted)
687 
688  real(preal), dimension(6) :: math_sym33to6
689  real(preal), dimension(3,3), intent(in) :: m33
690  logical, optional, intent(in) :: weighted
691 
692  real(preal), dimension(6) :: w
693  integer :: i
694 
695  if(present(weighted)) then
696  w = merge(nrmmandel,1.0_preal,weighted)
697  else
698  w = nrmmandel
699  endif
700 
701  do i = 1, 6
702  math_sym33to6(i) = w(i)*m33(mapnye(1,i),mapnye(2,i))
703  enddo
704 
705 end function math_sym33to6
706 
707 
708 !--------------------------------------------------------------------------------------------------
711 ! components according to Mandel. Advisable for matrix operations.
712 ! Unweighted conversion only changes order according to Nye
713 !--------------------------------------------------------------------------------------------------
714 pure function math_6tosym33(v6,weighted)
715 
716  real(preal), dimension(3,3) :: math_6tosym33
717  real(preal), dimension(6), intent(in) :: v6
718  logical, optional, intent(in) :: weighted
719 
720  real(preal), dimension(6) :: w
721  integer :: i
722 
723  if(present(weighted)) then
724  w = merge(invnrmmandel,1.0_preal,weighted)
725  else
726  w = invnrmmandel
727  endif
728 
729  do i=1,6
730  math_6tosym33(mapnye(1,i),mapnye(2,i)) = w(i)*v6(i)
731  math_6tosym33(mapnye(2,i),mapnye(1,i)) = w(i)*v6(i)
732  enddo
733 
734 end function math_6tosym33
735 
736 
737 !--------------------------------------------------------------------------------------------------
739 !--------------------------------------------------------------------------------------------------
740 pure function math_3333to99(m3333)
741 
742  real(preal), dimension(9,9) :: math_3333to99
743  real(preal), dimension(3,3,3,3), intent(in) :: m3333
744 
745  integer :: i,j
746 
747  do i=1,9; do j=1,9
748  math_3333to99(i,j) = m3333(mapplain(1,i),mapplain(2,i),mapplain(1,j),mapplain(2,j))
749  enddo; enddo
750 
751 end function math_3333to99
752 
753 
754 !--------------------------------------------------------------------------------------------------
756 !--------------------------------------------------------------------------------------------------
757 pure function math_99to3333(m99)
758 
759  real(preal), dimension(3,3,3,3) :: math_99to3333
760  real(preal), dimension(9,9), intent(in) :: m99
761 
762  integer :: i,j
763 
764  do i=1,9; do j=1,9
765  math_99to3333(mapplain(1,i),mapplain(2,i),mapplain(1,j),mapplain(2,j)) = m99(i,j)
766  enddo; enddo
767 
768 end function math_99to3333
769 
770 
771 !--------------------------------------------------------------------------------------------------
774 ! components according to Mandel. Advisable for matrix operations.
775 ! Unweighted conversion only rearranges order according to Nye
776 !--------------------------------------------------------------------------------------------------
777 pure function math_sym3333to66(m3333,weighted)
778 
779  real(preal), dimension(6,6) :: math_sym3333to66
780  real(preal), dimension(3,3,3,3), intent(in) :: m3333
781  logical, optional, intent(in) :: weighted
782 
783  real(preal), dimension(6) :: w
784  integer :: i,j
785 
786  if(present(weighted)) then
787  w = merge(nrmmandel,1.0_preal,weighted)
788  else
789  w = nrmmandel
790  endif
791 
792  do i=1,6; do j=1,6
793  math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapnye(1,i),mapnye(2,i),mapnye(1,j),mapnye(2,j))
794  enddo; enddo
795 
796 end function math_sym3333to66
797 
798 
799 !--------------------------------------------------------------------------------------------------
802 ! components according to Mandel. Advisable for matrix operations.
803 ! Unweighted conversion only rearranges order according to Nye
804 !--------------------------------------------------------------------------------------------------
805 pure function math_66tosym3333(m66,weighted)
806 
807  real(preal), dimension(3,3,3,3) :: math_66tosym3333
808  real(preal), dimension(6,6), intent(in) :: m66
809  logical, optional, intent(in) :: weighted
810 
811  real(preal), dimension(6) :: w
812  integer :: i,j
813 
814  if(present(weighted)) then
815  w = merge(invnrmmandel,1.0_preal,weighted)
816  else
817  w = invnrmmandel
818  endif
819 
820  do i=1,6; do j=1,6
821  math_66tosym3333(mapnye(1,i),mapnye(2,i),mapnye(1,j),mapnye(2,j)) = w(i)*w(j)*m66(i,j)
822  math_66tosym3333(mapnye(2,i),mapnye(1,i),mapnye(1,j),mapnye(2,j)) = w(i)*w(j)*m66(i,j)
823  math_66tosym3333(mapnye(1,i),mapnye(2,i),mapnye(2,j),mapnye(1,j)) = w(i)*w(j)*m66(i,j)
824  math_66tosym3333(mapnye(2,i),mapnye(1,i),mapnye(2,j),mapnye(1,j)) = w(i)*w(j)*m66(i,j)
825  enddo; enddo
826 
827 end function math_66tosym3333
828 
829 
830 !--------------------------------------------------------------------------------------------------
832 !--------------------------------------------------------------------------------------------------
833 pure function math_voigt66to3333(m66)
834 
835  real(preal), dimension(3,3,3,3) :: math_voigt66to3333
836  real(preal), dimension(6,6), intent(in) :: m66
837  integer :: i,j
838 
839  do i=1,6; do j=1, 6
840  math_voigt66to3333(mapvoigt(1,i),mapvoigt(2,i),mapvoigt(1,j),mapvoigt(2,j)) = m66(i,j)
841  math_voigt66to3333(mapvoigt(2,i),mapvoigt(1,i),mapvoigt(1,j),mapvoigt(2,j)) = m66(i,j)
842  math_voigt66to3333(mapvoigt(1,i),mapvoigt(2,i),mapvoigt(2,j),mapvoigt(1,j)) = m66(i,j)
843  math_voigt66to3333(mapvoigt(2,i),mapvoigt(1,i),mapvoigt(2,j),mapvoigt(1,j)) = m66(i,j)
844  enddo; enddo
845 
846 end function math_voigt66to3333
847 
848 
849 !--------------------------------------------------------------------------------------------------
851 !--------------------------------------------------------------------------------------------------
852 real(pReal) function math_samplegaussvar(meanvalue, stddev, width)
853 
854  real(preal), intent(in) :: meanvalue, & !< meanvalue of gauss distribution
855  stddev
856  real(preal), intent(in), optional :: width
857 
858  real(preal), dimension(2) :: rnd ! random numbers
859  real(preal) :: scatter, & ! normalized scatter around meanvalue
860  width_
861 
862  if (abs(stddev) < tol_math_check) then
863  math_samplegaussvar = meanvalue
864  else
865  if (present(width)) then
866  width_ = width
867  else
868  width_ = 3.0_preal ! use +-3*sigma as default scatter
869  endif
870 
871  do
872  call random_number(rnd)
873  scatter = width_ * (2.0_preal * rnd(1) - 1.0_preal)
874  if (rnd(2) <= exp(-0.5_preal * scatter ** 2.0_preal)) exit ! test if scattered value is drawn
875  enddo
876 
877  math_samplegaussvar = scatter * stddev
878  endif
879 
880 end function math_samplegaussvar
881 
882 
883 !--------------------------------------------------------------------------------------------------
885 ! ToDo: has wrong oder of arguments
886 !--------------------------------------------------------------------------------------------------
887 subroutine math_eigh(m,w,v,error)
888 
889  real(pReal), dimension(:,:), intent(in) :: m
890  real(pReal), dimension(size(m,1)), intent(out) :: w
891  real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: v
892 
893  logical, intent(out) :: error
894  integer :: ierr
895  real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
896  external :: &
897  dsyev
898 
899  v = m ! copy matrix to input (doubles as output) array
900  call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
901  error = (ierr /= 0)
902 
903 end subroutine math_eigh
904 
905 
906 !--------------------------------------------------------------------------------------------------
912 ! ToDo: has wrong oder of arguments
913 !--------------------------------------------------------------------------------------------------
914 subroutine math_eigh33(m,w,v)
915 
916  real(pReal), dimension(3,3),intent(in) :: m
917  real(pReal), dimension(3), intent(out) :: w
918  real(pReal), dimension(3,3),intent(out) :: v
919 
920  real(pReal) :: T, U, norm, threshold
921  logical :: error
922 
923  w = math_eigvalsh33(m)
924 
925  v(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
926  m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), &
927  m(1, 2)**2]
928 
929  t = maxval(abs(w))
930  u = max(t, t**2)
931  threshold = sqrt(5.68e-14_preal * u**2)
932 
933  v(1:3,1) = [ v(1,2) + m(1, 3) * w(1), &
934  v(2,2) + m(2, 3) * w(1), &
935  (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)]
936  norm = norm2(v(1:3, 1))
937  fallback1: if(norm < threshold) then
938  call math_eigh(m,w,v,error)
939  else fallback1
940  v(1:3,1) = v(1:3, 1) / norm
941  v(1:3,2) = [ v(1,2) + m(1, 3) * w(2), &
942  v(2,2) + m(2, 3) * w(2), &
943  (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)]
944  norm = norm2(v(1:3, 2))
945  fallback2: if(norm < threshold) then
946  call math_eigh(m,w,v,error)
947  else fallback2
948  v(1:3,2) = v(1:3, 2) / norm
949  v(1:3,3) = math_cross(v(1:3,1),v(1:3,2))
950  endif fallback2
951  endif fallback1
952 
953 end subroutine math_eigh33
954 
955 
956 
957 
958 !--------------------------------------------------------------------------------------------------
960 !--------------------------------------------------------------------------------------------------
961 function math_rotationalpart(m)
962 
963  real(preal), intent(in), dimension(3,3) :: m
964  real(preal), dimension(3,3) :: math_rotationalpart
965  real(preal), dimension(3,3) :: u , uinv
966 
967  u = eigenvectorbasis(matmul(transpose(m),m))
968  uinv = math_inv33(u)
969 
970  inversionfailed: if (all(deq0(uinv))) then
971  math_rotationalpart = math_i3
972  call io_warning(650)
973  else inversionfailed
974  math_rotationalpart = matmul(m,uinv)
975  endif inversionfailed
976 
977 contains
978  !--------------------------------------------------------------------------------------------------
980  !--------------------------------------------------------------------------------------------------
981  pure function eigenvectorbasis(m)
982 
983  real(preal), dimension(3,3) :: eigenvectorbasis
984  real(preal), dimension(3,3), intent(in) :: m
985 
986  real(preal), dimension(3) :: i, v
987  real(preal) :: p, q, rho, phi
988  real(preal), parameter :: tol=1.e-14_preal
989  real(preal), dimension(3,3,3) :: n, eb
990 
991  i = math_invariantssym33(m)
992 
993  p = i(2)-i(1)**2.0_preal/3.0_preal
994  q = -2.0_preal/27.0_preal*i(1)**3.0_preal+product(i(1:2))/3.0_preal-i(3)
995 
996  threesimilareigvals: if(all(abs([p,q]) < tol)) then
997  v = i(1)/3.0_preal
998  ! this is not really correct, but at least the basis is correct
999  eb = 0.0_preal
1000  eb(1,1,1)=1.0_preal
1001  eb(2,2,2)=1.0_preal
1002  eb(3,3,3)=1.0_preal
1003  else threesimilareigvals
1004  rho=sqrt(-3.0_preal*p**3.0_preal)/9.0_preal
1005  phi=acos(math_clip(-q/rho*0.5_preal,-1.0_preal,1.0_preal))
1006  v = 2.0_preal*rho**(1.0_preal/3.0_preal)* [cos((phi )/3.0_preal), &
1007  cos((phi+2.0_preal*pi)/3.0_preal), &
1008  cos((phi+4.0_preal*pi)/3.0_preal) &
1009  ] + i(1)/3.0_preal
1010  n(1:3,1:3,1) = m-v(1)*math_i3
1011  n(1:3,1:3,2) = m-v(2)*math_i3
1012  n(1:3,1:3,3) = m-v(3)*math_i3
1013  twosimilareigvals: if(abs(v(1)-v(2)) < tol) then
1014  eb(1:3,1:3,3) = matmul(n(1:3,1:3,1),n(1:3,1:3,2))/((v(3)-v(1))*(v(3)-v(2)))
1015  eb(1:3,1:3,1) = math_i3-eb(1:3,1:3,3)
1016  eb(1:3,1:3,2) = 0.0_preal
1017  elseif (abs(v(2)-v(3)) < tol) then twosimilareigvals
1018  eb(1:3,1:3,1) = matmul(n(1:3,1:3,2),n(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
1019  eb(1:3,1:3,2) = math_i3-eb(1:3,1:3,1)
1020  eb(1:3,1:3,3) = 0.0_preal
1021  elseif (abs(v(3)-v(1)) < tol) then twosimilareigvals
1022  eb(1:3,1:3,2) = matmul(n(1:3,1:3,3),n(1:3,1:3,1))/((v(2)-v(3))*(v(2)-v(1)))
1023  eb(1:3,1:3,3) = math_i3-eb(1:3,1:3,2)
1024  eb(1:3,1:3,1) = 0.0_preal
1025  else twosimilareigvals
1026  eb(1:3,1:3,1) = matmul(n(1:3,1:3,2),n(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
1027  eb(1:3,1:3,2) = matmul(n(1:3,1:3,3),n(1:3,1:3,1))/((v(2)-v(3))*(v(2)-v(1)))
1028  eb(1:3,1:3,3) = matmul(n(1:3,1:3,1),n(1:3,1:3,2))/((v(3)-v(1))*(v(3)-v(2)))
1029  endif twosimilareigvals
1030  endif threesimilareigvals
1031 
1032  eigenvectorbasis = sqrt(v(1)) * eb(1:3,1:3,1) &
1033  + sqrt(v(2)) * eb(1:3,1:3,2) &
1034  + sqrt(v(3)) * eb(1:3,1:3,3)
1035 
1036  end function eigenvectorbasis
1037 
1038 end function math_rotationalpart
1039 
1040 
1041 !--------------------------------------------------------------------------------------------------
1043 ! will return NaN on error
1044 !--------------------------------------------------------------------------------------------------
1045 function math_eigvalsh(m)
1046 
1047  real(preal), dimension(:,:), intent(in) :: m
1048  real(preal), dimension(size(m,1)) :: math_eigvalsh
1049 
1050  real(preal), dimension(size(m,1),size(m,1)) :: m_
1051  integer :: ierr
1052  real(preal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
1053  external :: &
1054  dsyev
1055 
1056  m_= m ! copy matrix to input (will be destroyed)
1057  call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr)
1058  if (ierr /= 0) math_eigvalsh = ieee_value(1.0_preal,ieee_quiet_nan)
1059 
1060 end function math_eigvalsh
1061 
1062 
1063 !--------------------------------------------------------------------------------------------------
1069 !--------------------------------------------------------------------------------------------------
1070 function math_eigvalsh33(m)
1071 
1072  real(preal), intent(in), dimension(3,3) :: m
1073  real(preal), dimension(3) :: math_eigvalsh33,i
1074  real(preal) :: p, q, rho, phi
1075  real(preal), parameter :: tol=1.e-14_preal
1076 
1077  i = math_invariantssym33(m) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206
1078 
1079  p = i(2)-i(1)**2.0_preal/3.0_preal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
1080  q = product(i(1:2))/3.0_preal &
1081  - 2.0_preal/27.0_preal*i(1)**3.0_preal &
1082  - i(3) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
1083 
1084  if(all(abs([p,q]) < tol)) then
1085  math_eigvalsh33 = math_eigvalsh(m)
1086  else
1087  rho=sqrt(-3.0_preal*p**3.0_preal)/9.0_preal
1088  phi=acos(math_clip(-q/rho*0.5_preal,-1.0_preal,1.0_preal))
1089  math_eigvalsh33 = 2.0_preal*rho**(1.0_preal/3.0_preal)* &
1090  [cos(phi/3.0_preal), &
1091  cos((phi+2.0_preal*pi)/3.0_preal), &
1092  cos((phi+4.0_preal*pi)/3.0_preal) &
1093  ] &
1094  + i(1)/3.0_preal
1095  endif
1096 
1097 end function math_eigvalsh33
1098 
1099 
1100 !--------------------------------------------------------------------------------------------------
1102 !--------------------------------------------------------------------------------------------------
1103 pure function math_invariantssym33(m)
1104 
1105  real(preal), dimension(3,3), intent(in) :: m
1106  real(preal), dimension(3) :: math_invariantssym33
1107 
1108  math_invariantssym33(1) = math_trace33(m)
1109  math_invariantssym33(2) = m(1,1)*m(2,2) + m(1,1)*m(3,3) + m(2,2)*m(3,3) &
1110  -(m(1,2)**2 + m(1,3)**2 + m(2,3)**2)
1111  math_invariantssym33(3) = math_detsym33(m)
1112 
1113 end function math_invariantssym33
1114 
1115 
1116 !--------------------------------------------------------------------------------------------------
1118 !--------------------------------------------------------------------------------------------------
1119 integer pure function math_factorial(n)
1120 
1121  integer, intent(in) :: n
1122 
1123  math_factorial = product(math_range(n))
1124 
1125 end function math_factorial
1126 
1127 
1128 !--------------------------------------------------------------------------------------------------
1130 !--------------------------------------------------------------------------------------------------
1131 integer pure function math_binomial(n,k)
1132 
1133  integer, intent(in) :: n, k
1134  integer :: i, k_, n_
1135 
1136  k_ = min(k,n-k)
1137  n_ = n
1138  math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n
1139  do i = 1, k_
1140  math_binomial = (math_binomial * n_)/i
1141  n_ = n_ -1
1142  enddo
1143 
1144 end function math_binomial
1145 
1146 
1147 !--------------------------------------------------------------------------------------------------
1149 !--------------------------------------------------------------------------------------------------
1150 integer pure function math_multinomial(alpha)
1151 
1152  integer, intent(in), dimension(:) :: alpha
1153  integer :: i
1154 
1155  math_multinomial = 1
1156  do i = 1, size(alpha)
1157  math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
1158  enddo
1159 
1160 end function math_multinomial
1161 
1162 
1163 !--------------------------------------------------------------------------------------------------
1165 !--------------------------------------------------------------------------------------------------
1166 real(preal) pure function math_voltetrahedron(v1,v2,v3,v4)
1167 
1168  real(preal), dimension (3), intent(in) :: v1,v2,v3,v4
1169  real(preal), dimension (3,3) :: m
1170 
1171  m(1:3,1) = v1-v2
1172  m(1:3,2) = v1-v3
1173  m(1:3,3) = v1-v4
1174 
1175  math_voltetrahedron = abs(math_det33(m))/6.0_preal
1176 
1177 end function math_voltetrahedron
1178 
1179 
1180 !--------------------------------------------------------------------------------------------------
1182 !--------------------------------------------------------------------------------------------------
1183 real(preal) pure function math_areatriangle(v1,v2,v3)
1184 
1185  real(preal), dimension (3), intent(in) :: v1,v2,v3
1186 
1187  math_areatriangle = 0.5_preal * norm2(math_cross(v1-v2,v1-v3))
1188 
1189 end function math_areatriangle
1190 
1191 
1192 !--------------------------------------------------------------------------------------------------
1194 ! Will return NaN if left > right
1195 !--------------------------------------------------------------------------------------------------
1196 real(preal) pure elemental function math_clip(a, left, right)
1197 
1198  real(preal), intent(in) :: a
1199  real(preal), intent(in), optional :: left, right
1200 
1201  math_clip = a
1202  if (present(left)) math_clip = max(left,math_clip)
1203  if (present(right)) math_clip = min(right,math_clip)
1204  if (present(left) .and. present(right)) &
1205  math_clip = merge(ieee_value(1.0_preal,ieee_quiet_nan),math_clip, left>right)
1206 
1207 end function math_clip
1208 
1209 
1210 !--------------------------------------------------------------------------------------------------
1212 !--------------------------------------------------------------------------------------------------
1213 subroutine unittest
1214 
1215  integer, dimension(2,4) :: &
1216  sort_in_ = reshape([+1,+5, +5,+6, -1,-1, +3,-2],[2,4])
1217  integer, dimension(2,4), parameter :: &
1218  sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4])
1219 
1220  integer, dimension(5) :: range_out_ = [1,2,3,4,5]
1221  integer, dimension(3) :: ijk
1222 
1223  real(preal) :: det
1224  real(preal), dimension(3) :: v3_1,v3_2,v3_3,v3_4
1225  real(preal), dimension(6) :: v6
1226  real(preal), dimension(9) :: v9
1227  real(preal), dimension(3,3) :: t33,t33_2
1228  real(preal), dimension(6,6) :: t66
1229  real(preal), dimension(9,9) :: t99,t99_2
1230  real(preal), dimension(:,:), &
1231  allocatable :: txx,txx_2
1232  real(preal) :: r
1233  integer :: d
1234  logical :: e
1235 
1236  if (any(abs([1.0_preal,2.0_preal,2.0_preal,3.0_preal,3.0_preal,3.0_preal] - &
1237  math_expand([1.0_preal,2.0_preal,3.0_preal],[1,2,3,0])) > tol_math_check)) &
1238  call io_error(0,ext_msg='math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]')
1239 
1240  if (any(abs([1.0_preal,2.0_preal,2.0_preal] - &
1241  math_expand([1.0_preal,2.0_preal,3.0_preal],[1,2])) > tol_math_check)) &
1242  call io_error(0,ext_msg='math_expand [1,2,3] by [1,2] => [1,2,2]')
1243 
1244  if (any(abs([1.0_preal,2.0_preal,2.0_preal,1.0_preal,1.0_preal,1.0_preal] - &
1245  math_expand([1.0_preal,2.0_preal],[1,2,3])) > tol_math_check)) &
1246  call io_error(0,ext_msg='math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]')
1247 
1248  call math_sort(sort_in_,1,3,2)
1249  if(any(sort_in_ /= sort_out_)) &
1250  call io_error(0,ext_msg='math_sort')
1251 
1252  if(any(math_range(5) /= range_out_)) &
1253  call io_error(0,ext_msg='math_range')
1254 
1255  if(any(dneq(math_exp33(math_i3,0),math_i3))) &
1256  call io_error(0,ext_msg='math_exp33(math_I3,1)')
1257  if(any(dneq(math_exp33(math_i3,256),exp(1.0_preal)*math_i3))) &
1258  call io_error(0,ext_msg='math_exp33(math_I3,256)')
1259 
1260  call random_number(v9)
1261  if(any(dneq(math_33to9(math_9to33(v9)),v9))) &
1262  call io_error(0,ext_msg='math_33to9/math_9to33')
1263 
1264  call random_number(t99)
1265  if(any(dneq(math_3333to99(math_99to3333(t99)),t99))) &
1266  call io_error(0,ext_msg='math_3333to99/math_99to3333')
1267 
1268  call random_number(v6)
1269  if(any(dneq(math_sym33to6(math_6tosym33(v6)),v6))) &
1270  call io_error(0,ext_msg='math_sym33to6/math_6toSym33')
1271 
1272  call random_number(t66)
1273  if(any(dneq(math_sym3333to66(math_66tosym3333(t66)),t66))) &
1274  call io_error(0,ext_msg='math_sym3333to66/math_66toSym3333')
1275 
1276  call random_number(v6)
1277  if(any(dneq0(math_6tosym33(v6) - math_symmetric33(math_6tosym33(v6))))) &
1278  call io_error(0,ext_msg='math_symmetric33')
1279 
1280  call random_number(v3_1)
1281  call random_number(v3_2)
1282  call random_number(v3_3)
1283  call random_number(v3_4)
1284 
1285  if(dneq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
1286  math_voltetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_preal)) &
1287  call io_error(0,ext_msg='math_volTetrahedron')
1288 
1289  call random_number(t33)
1290  if(dneq(math_det33(math_symmetric33(t33)),math_detsym33(math_symmetric33(t33)),tol=1.0e-12_preal)) &
1291  call io_error(0,ext_msg='math_det33/math_detSym33')
1292 
1293  if(any(dneq0(math_identity2nd(3),math_inv33(math_i3)))) &
1294  call io_error(0,ext_msg='math_inv33(math_I3)')
1295 
1296  do while(abs(math_det33(t33))<1.0e-9_preal)
1297  call random_number(t33)
1298  enddo
1299  if(any(dneq0(matmul(t33,math_inv33(t33)) - math_identity2nd(3),tol=1.0e-9_preal))) &
1300  call io_error(0,ext_msg='math_inv33')
1301 
1302  call math_invert33(t33_2,det,e,t33)
1303  if(any(dneq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_preal)) .or. e) &
1304  call io_error(0,ext_msg='math_invert33: T:T^-1 != I')
1305  if(dneq(det,math_det33(t33),tol=1.0e-12_preal)) &
1306  call io_error(0,ext_msg='math_invert33 (determinant)')
1307 
1308  call math_invert(t33_2,e,t33)
1309  if(any(dneq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_preal)) .or. e) &
1310  call io_error(0,ext_msg='math_invert t33')
1311 
1312  do while(math_det33(t33)<1.0e-2_preal) ! O(det(F)) = 1
1313  call random_number(t33)
1314  enddo
1315  t33_2 = math_rotationalpart(transpose(t33))
1316  t33 = math_rotationalpart(t33)
1317  if(any(dneq0(matmul(t33_2,t33) - math_i3,tol=1.0e-10_preal))) &
1318  call io_error(0,ext_msg='math_rotationalPart')
1319 
1320  call random_number(r)
1321  d = int(r*5.0_preal) + 1
1322  txx = math_identity2nd(d)
1323  allocate(txx_2(d,d))
1324  call math_invert(txx_2,e,txx)
1325  if(any(dneq0(txx_2,txx) .or. e)) &
1326  call io_error(0,ext_msg='math_invert(txx)/math_identity2nd')
1327 
1328  call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix
1329  if(any(dneq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_preal)) .or. e) &
1330  call io_error(0,ext_msg='math_invert(t99)')
1331 
1332  if(any(dneq(math_clip([4.0_preal,9.0_preal],5.0_preal,6.5_preal),[5.0_preal,6.5_preal]))) &
1333  call io_error(0,ext_msg='math_clip')
1334 
1335  if(math_factorial(10) /= 3628800) &
1336  call io_error(0,ext_msg='math_factorial')
1337 
1338  if(math_binomial(49,6) /= 13983816) &
1339  call io_error(0,ext_msg='math_binomial')
1340 
1341  ijk = cshift([1,2,3],int(r*1.0e2_preal))
1342  if(dneq(math_levicivita(ijk(1),ijk(2),ijk(3)),+1.0_preal)) &
1343  call io_error(0,ext_msg='math_LeviCivita(even)')
1344  ijk = cshift([3,2,1],int(r*2.0e2_preal))
1345  if(dneq(math_levicivita(ijk(1),ijk(2),ijk(3)),-1.0_preal)) &
1346  call io_error(0,ext_msg='math_LeviCivita(odd)')
1347  ijk = cshift([2,2,1],int(r*2.0e2_preal))
1348  if(dneq0(math_levicivita(ijk(1),ijk(2),ijk(3))))&
1349  call io_error(0,ext_msg='math_LeviCivita')
1350 
1351 end subroutine unittest
1352 
1353 end module math
math::math_expand
pure real(preal) function, dimension(sum(how)) math_expand(what, how)
vector expansion
Definition: math.f90:207
math::math_symmetric33
pure real(preal) function, dimension(3, 3) math_symmetric33(m)
symmetrize a 3x3 matrix
Definition: math.f90:546
math::invnrmmandel
real(preal), dimension(6), parameter, private invnrmmandel
backward weighting for Mandel notation
Definition: math.f90:44
math::math_binomial
integer pure function math_binomial(n, k)
binomial coefficient
Definition: math.f90:1132
eigenvectorbasis
pure real(preal) function, dimension(3, 3) eigenvectorbasis(m)
eigenvector basis of positive-definite 3x3 matrix
Definition: math.f90:982
math::math_voigt66to3333
pure real(preal) function, dimension(3, 3, 3, 3) math_voigt66to3333(m66)
convert 66 Voigt matrix into symmetric 3x3x3x3 matrix
Definition: math.f90:834
math::math_outer
pure real(preal) function, dimension(size(a, 1), size(b, 1)) math_outer(A, B)
outer product of arbitrary sized vectors (A ⊗ B / i,j)
Definition: math.f90:327
math::math_99to3333
pure real(preal) function, dimension(3, 3, 3, 3) math_99to3333(m99)
convert 9x9 matrix into 3x3x3x3 matrix
Definition: math.f90:758
math::math_3333to99
pure real(preal) function, dimension(9, 9) math_3333to99(m3333)
convert 3x3x3x3 matrix into 9x9 matrix
Definition: math.f90:741
math::math_identity2nd
pure real(preal) function, dimension(d, d) math_identity2nd(d)
second rank identity tensor of specified dimension
Definition: math.f90:240
math::indeg
real(preal), parameter indeg
conversion from radian into degree
Definition: math.f90:28
numerics::randomseed
integer, public, protected randomseed
fixed seeding for pseudo-random number generator, Default 0: use random seed
Definition: numerics.f90:1470
math::math_sym3333to66
pure real(preal) function, dimension(6, 6) math_sym3333to66(m3333, weighted)
convert symmetric 3x3x3x3 matrix into 6x6 matrix
Definition: math.f90:778
math::math_spherical33
pure real(preal) function, dimension(3, 3) math_spherical33(m)
hydrostatic part of a 3x3 matrix
Definition: math.f90:585
math::math_33to9
pure real(preal) function, dimension(9) math_33to9(m33)
convert 3x3 matrix into vector 9
Definition: math.f90:650
math::math_eigvalsh
real(preal) function, dimension(size(m, 1)) math_eigvalsh(m)
Eigenvalues of symmetric matrix.
Definition: math.f90:1046
math::inrad
real(preal), parameter inrad
conversion from degree into radian
Definition: math.f90:29
math::math_range
pure integer function, dimension(n) math_range(N)
range of integers starting at one
Definition: math.f90:226
math::math_voltetrahedron
real(preal) pure function math_voltetrahedron(v1, v2, v3, v4)
volume of tetrahedron given by four vertices
Definition: math.f90:1167
math::math_eye
Definition: math.f90:80
math::math_eigh
subroutine math_eigh(m, w, v, error)
eigenvalues and eigenvectors of symmetric matrix
Definition: math.f90:888
math::math_symmetric66
pure real(preal) function, dimension(6, 6) math_symmetric66(m)
symmetrize a 6x6 matrix
Definition: math.f90:559
prec
setting precision for real and int type
Definition: prec.f90:13
math::unittest
subroutine, private unittest
check correctness of some math functions
Definition: math.f90:1214
math::mapnye
integer, dimension(2, 6), parameter, private mapnye
arrangement in Nye notation.
Definition: math.f90:47
math::math_areatriangle
real(preal) pure function math_areatriangle(v1, v2, v3)
area of triangle given by three vertices
Definition: math.f90:1184
math::math_det33
real(preal) pure function math_det33(m)
determinant of a 3x3 matrix
Definition: math.f90:623
math::math_66tosym3333
pure real(preal) function, dimension(3, 3, 3, 3) math_66tosym3333(m66, weighted)
convert 66 matrix into symmetric 3x3x3x3 matrix
Definition: math.f90:806
math::math_mul3333xx33
pure real(preal) function, dimension(3, 3) math_mul3333xx33(A, B)
matrix double contraction 3333x33 = 33 (ijkl,kl)
Definition: math.f90:368
math::twopiimg
complex(preal), parameter twopiimg
Re(0.0), Im(2xPi)
Definition: math.f90:30
math::math_invariantssym33
pure real(preal) function, dimension(3) math_invariantssym33(m)
invariants of symmetrix 3x3 matrix
Definition: math.f90:1104
math::mapplain
integer, dimension(2, 9), parameter, private mapplain
arrangement in Plain notation
Definition: math.f90:67
math::math_mul3333xx3333
pure real(preal) function, dimension(3, 3, 3, 3) math_mul3333xx3333(A, B)
matrix multiplication 3333x3333 = 3333 (ijkl,klmn)
Definition: math.f90:385
io
input/output functions, partly depending on chosen solver
Definition: IO.f90:12
math::nrmmandel
real(preal), dimension(6), parameter, private nrmmandel
forward weighting for Mandel notation
Definition: math.f90:39
prec::preal
integer, parameter preal
number with 15 significant digits, up to 1e+-307 (typically 64 bit)
Definition: prec.f90:20
math::math_invsym3333
real(preal) function, dimension(3, 3, 3, 3) math_invsym3333(A)
Inversion of symmetriced 3x3x3x3 matrix.
Definition: math.f90:489
math::math_levicivita
real(preal) pure function math_levicivita(i, j, k)
permutation tensor e_ijk
Definition: math.f90:280
math::math_multinomial
integer pure function math_multinomial(alpha)
multinomial coefficient
Definition: math.f90:1151
math::math_inner
real(preal) pure function math_inner(A, B)
inner product of arbitrary sized vectors (A · B / i,i)
Definition: math.f90:343
math::math_samplegaussvar
real(preal) function math_samplegaussvar(meanvalue, stddev, width)
draw a random sample from Gauss variable
Definition: math.f90:853
math::math_invert33
pure subroutine math_invert33(InvA, DetA, error, A)
Cramer inversion of 3x3 matrix (subroutine)
Definition: math.f90:454
math::math_eigvalsh33
real(preal) function, dimension(3) math_eigvalsh33(m)
eigenvalues of symmetric 3x3 matrix using an analytical expression
Definition: math.f90:1071
math::math_clip
real(preal) pure elemental function math_clip(a, left, right)
limits a scalar value to a certain range (either one or two sided)
Definition: math.f90:1197
qsort_partition
integer function qsort_partition(a, istart, iend, sort)
Partitioning required for quicksort.
Definition: math.f90:168
math::math_invert
subroutine math_invert(InvA, error, A)
invert quadratic matrix of arbitrary dimension
Definition: math.f90:521
math::math_trace33
real(preal) pure function math_trace33(m)
trace of a 3x3 matrix
Definition: math.f90:611
math::math_delta
real(preal) pure function math_delta(i, j)
kronecker delta function d_ij
Definition: math.f90:300
math::math_tensordot
real(preal) pure function math_tensordot(A, B)
double contraction of 3x3 matrices (A : B / ij,ij)
Definition: math.f90:356
math::math_sym33to6
pure real(preal) function, dimension(6) math_sym33to6(m33, weighted)
convert symmetric 3x3 matrix into 6 vector
Definition: math.f90:687
math
Mathematical library, including random number generation and tensor representations.
Definition: math.f90:12
math::math_eigh33
subroutine math_eigh33(m, w, v)
eigenvalues and eigenvectors of symmetric 3x3 matrix using an analytical expression and the general L...
Definition: math.f90:915
math::math_rotationalpart
real(preal) function, dimension(3, 3) math_rotationalpart(m)
rotational part from polar decomposition of 3x3 tensor
Definition: math.f90:962
math::math_deviatoric33
pure real(preal) function, dimension(3, 3) math_deviatoric33(m)
deviatoric part of a 3x3 matrix
Definition: math.f90:598
math::math_9to33
pure real(preal) function, dimension(3, 3) math_9to33(v9)
convert 9 vector into 3x3 matrix
Definition: math.f90:667
io::unittest
subroutine unittest
check correctness of some IO functions
Definition: IO.f90:662
math::math_factorial
integer pure function math_factorial(n)
factorial
Definition: math.f90:1120
math::math_sort
recursive subroutine math_sort(a, istart, iend, sortDim)
Quicksort algorithm for two-dimensional integer arrays.
Definition: math.f90:132
math::math_init
subroutine math_init
initialization of random seed generator and internal checks
Definition: math.f90:95
math::pi
real(preal), parameter pi
ratio of a circle's circumference to its diameter
Definition: math.f90:27
math::math_cross
pure real(preal) function, dimension(3) math_cross(A, B)
cross product a x b
Definition: math.f90:312
math::math_6tosym33
pure real(preal) function, dimension(3, 3) math_6tosym33(v6, weighted)
convert 6 vector into symmetric 3x3 matrix
Definition: math.f90:715
math::math_inv33
pure real(preal) function, dimension(3, 3) math_inv33(A)
Cramer inversion of 3x3 matrix (function)
Definition: math.f90:435
numerics
Managing of parameters related to numerics.
Definition: numerics.f90:10
math::math_detsym33
real(preal) pure function math_detsym33(m)
determinant of a symmetric 3x3 matrix
Definition: math.f90:637
math::math_skew33
pure real(preal) function, dimension(3, 3) math_skew33(m)
skew part of a 3x3 matrix
Definition: math.f90:572
math::mapvoigt
integer, dimension(2, 6), parameter, private mapvoigt
arrangement in Voigt notation
Definition: math.f90:57
math::math_exp33
pure real(preal) function, dimension(3, 3) math_exp33(A, n)
3x3 matrix exponential up to series approximation order n (default 5)
Definition: math.f90:402
math::math_i3
real(preal), dimension(3, 3), parameter math_i3
3x3 Identity
Definition: math.f90:32
math::math_identity4th
pure real(preal) function, dimension(d, d, d, d) math_identity4th(d)
symmetric fourth rank identity tensor of specified dimension
Definition: math.f90:258