From cfaa0e696da95e6a84303abc8ff34f225604967a Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 15 Feb 2008 12:42:27 +0000 Subject: [PATCH] exchanged TINY (from prec) with intrinsic "tiny" function (Fortran90) --- trunk/math.f90 | 858 ++++++++++++++++++++++++++++--------------------- 1 file changed, 484 insertions(+), 374 deletions(-) diff --git a/trunk/math.f90 b/trunk/math.f90 index 5c9f81054..0700b1d03 100644 --- a/trunk/math.f90 +++ b/trunk/math.f90 @@ -2,7 +2,8 @@ !############################################################## MODULE math !############################################################## - + + use prec, only: pReal,pInt implicit none @@ -15,34 +16,53 @@ 1.0_pReal,0.0_pReal,0.0_pReal, & 0.0_pReal,1.0_pReal,0.0_pReal, & 0.0_pReal,0.0_pReal,1.0_pReal /),(/3,3/)) -! *** Mandel notation *** - integer(pInt), dimension (2,6), parameter :: mapMandel = & - reshape((/& - 1,1, & - 2,2, & - 3,3, & - 1,2, & - 2,3, & - 1,3 & - /),(/2,6/)) - real(pReal), dimension(6), parameter :: nrmMandel = & - (/1.0_pReal,1.0_pReal,1.0_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal/) - real(pReal), dimension(6), parameter :: invnrmMandel = & - (/1.0_pReal,1.0_pReal,1.0_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal/) -! *** Voigt notation *** - integer(pInt), dimension (2,6), parameter :: mapVoigt = & - reshape((/& - 1,1, & - 2,2, & - 3,3, & - 2,3, & - 1,3, & - 1,2 & - /),(/2,6/)) - real(pReal), dimension(6), parameter :: nrmVoigt = & - (/1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal/) - real(pReal), dimension(6), parameter :: invnrmVoigt = & - (/1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal/) + +! *** Mandel notation *** + integer(pInt), dimension (2,6), parameter :: mapMandel = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + /),(/2,6/)) + + real(pReal), dimension(6), parameter :: nrmMandel = & + (/1.0_pReal,1.0_pReal,1.0_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal/) + real(pReal), dimension(6), parameter :: invnrmMandel = & + (/1.0_pReal,1.0_pReal,1.0_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal,0.7071067811865476_pReal/) + +! *** Voigt notation *** + integer(pInt), dimension (2,6), parameter :: mapVoigt = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 2,3, & + 1,3, & + 1,2 & + /),(/2,6/)) + + real(pReal), dimension(6), parameter :: nrmVoigt = & + (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/) + real(pReal), dimension(6), parameter :: invnrmVoigt = & + (/1.0_pReal,1.0_pReal,1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal/) + +! *** Plain notation *** + integer(pInt), dimension (2,9), parameter :: mapPlain = & + reshape((/& + 1,1, & + 1,2, & + 1,3, & + 2,1, & + 2,2, & + 2,3, & + 3,1, & + 3,2, & + 3,3 & + /),(/2,9/)) + CONTAINS @@ -95,8 +115,10 @@ integer(pInt) :: istart,iend,d,i,j,k,x,tmp d = size(a,1) ! number of linked data -! set the starting and ending points, and the pivot point - i = istart +! set the starting and ending points, and the pivot point + + i = istart + j = iend x = a(1,istart) do @@ -163,278 +185,295 @@ END FUNCTION + !************************************************************************** ! Cramer inversion of 3x3 matrix !************************************************************************** - SUBROUTINE math_invert3x3(A, InvA, DetA, error) - -! Bestimmung der Determinanten und Inversen einer 3x3-Matrix - -! A = Matrix A -! InvA = Inverse of A -! DetA = Determinant of A + SUBROUTINE math_invert3x3(A, InvA, DetA, error) + +! Bestimmung der Determinanten und Inversen einer 3x3-Matrix +! A = Matrix A +! InvA = Inverse of A +! DetA = Determinant of A ! error = logical - - use prec, only: pReal,pInt - implicit none - - logical, intent(out) :: error + + use prec, only: pReal,pInt + implicit none + + logical, intent(out) :: error + real(pReal),dimension(3,3),intent(in) :: A real(pReal),dimension(3,3),intent(out) :: InvA real(pReal), intent(out) :: DetA - - DetA = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) )& - - A(1,2) * ( A(2,1) * A(3,3) - A(2,3) * A(3,1) )& - + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) - - if (DetA <= 0.0000001) then - error = .true. + + DetA = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) )& + - A(1,2) * ( A(2,1) * A(3,3) - A(2,3) * A(3,1) )& + + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) + + if (DetA <= tiny(DetA)) then + error = .true. else - InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA - InvA(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA - InvA(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA - - InvA(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA - InvA(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA - InvA(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA - - InvA(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA - InvA(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA + InvA(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA + InvA(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA + InvA(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA + + InvA(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA + InvA(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA + InvA(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA + + InvA(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA + InvA(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA InvA(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA - error = .false. + error = .false. endif return - END SUBROUTINE + END SUBROUTINE + !************************************************************************** ! Gauss elimination to invert 6x6 matrix !************************************************************************** - SUBROUTINE math_invert6x6(A, InvA, AnzNegEW, error) - -! Invertieren einer 6x6-Matrix - -! A = Matrix A -! InvA = Inverse von A -! AnzNegEW = Anzahl der negativen Eigenwerte von A -! error = logical -! = false: Inversion wurde durchgefuehrt. -! = true: Die Inversion in SymGauss wurde wegen eines verschwindenen -! Pivotelement abgebrochen. - - use prec, only: pReal,pInt - implicit none - - real(pReal),dimension(6,6), intent(in) :: A - real(pReal),dimension(6,6), intent(out) :: InvA - integer(pInt), intent(out) :: AnzNegEW - logical, intent(out) :: error - - integer(pInt) i - real(pReal) LogAbsDetA - real(pReal),dimension(6,6) :: B - - InvA = 0.0_pReal - - forall(i = 1:6) InvA(i,i) = 1.0_pReal - B = A - CALL Gauss (B,InvA,LogAbsDetA,AnzNegEW,error) - - RETURN - - END SUBROUTINE math_invert6x6 - -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - SUBROUTINE Gauss (A,B,LogAbsDetA,NegHDK,error) - -! Loesung eines linearen Gleichungsssystem A * X = B mit Hilfe des -! GAUSS-Algorithmusses - -! Zur numerischen Stabilisierung wird eine Zeilen- und Spaltenpivotsuche -! durchgefuehrt. - -! Eingabeparameter: -! -! A(6,6) = Koeffizientenmatrix A -! B(6,6) = rechte Seiten B -! -! Ausgabeparameter: -! -! B(6,6) = Matrix der Unbekanntenvektoren X -! LogAbsDetA = 10-Logarithmus des Betrages der Determinanten von A -! NegHDK = Anzahl der negativen Hauptdiagonalkoeffizienten nach der -! Vorwaertszerlegung -! error = logical -! = false: Das Gleichungssystem wurde geloest. -! = true : Matrix A ist singulaer. - -! A und B werden veraendert! - - use prec, only: pReal,pInt - implicit none + SUBROUTINE math_invert(dimen,A, InvA, AnzNegEW, error) + +! Invertieren einer dimen x dimen - Matrix +! A = Matrix A +! InvA = Inverse von A +! AnzNegEW = Anzahl der negativen Eigenwerte von A +! error = logical +! = false: Inversion wurde durchgefuehrt. +! = true: Die Inversion in SymGauss wurde wegen eines verschwindenen +! Pivotelement abgebrochen. + + use prec, only: pReal,pInt + implicit none + + integer(pInt), intent(in) :: dimen + real(pReal),dimension(dimen,dimen), intent(in) :: A + real(pReal),dimension(dimen,dimen), intent(out) :: InvA + integer(pInt), intent(out) :: AnzNegEW + logical, intent(out) :: error + integer(pInt) i + real(pReal) LogAbsDetA + real(pReal),dimension(dimen,dimen) :: B + + InvA = math_identity2nd(dimen) + B = A + CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error) + + RETURN + + END SUBROUTINE math_invert + + + +! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + SUBROUTINE Gauss (dimen,A,B,LogAbsDetA,NegHDK,error) + +! Loesung eines linearen Gleichungsssystem A * X = B mit Hilfe des +! GAUSS-Algorithmus +! Zur numerischen Stabilisierung wird eine Zeilen- und Spaltenpivotsuche +! durchgefuehrt. + +! Eingabeparameter: +! +! A(dimen,dimen) = Koeffizientenmatrix A +! B(dimen,dimen) = rechte Seiten B +! +! Ausgabeparameter: +! +! B(dimen,dimen) = Matrix der Unbekanntenvektoren X +! LogAbsDetA = 10-Logarithmus des Betrages der Determinanten von A +! NegHDK = Anzahl der negativen Hauptdiagonalkoeffizienten nach der +! Vorwaertszerlegung +! error = logical +! = false: Das Gleichungssystem wurde geloest. +! = true : Matrix A ist singulaer. + +! A und B werden veraendert! + + use prec, only: pReal,pInt + implicit none + + logical error + integer (pInt) dimen,NegHDK + real(pReal) LogAbsDetA + real(pReal) A(dimen,dimen), B(dimen,dimen) + + INTENT (IN) dimen + INTENT (OUT) LogAbsDetA, NegHDK, error + INTENT (INOUT) A, B + + LOGICAL SortX + integer (pInt) PivotZeile, PivotSpalte, StoreI, I, IP1, J, K, L + integer (pInt) XNr(dimen) + real(pReal) AbsA, PivotWert, EpsAbs, Quote + real(pReal) StoreA(dimen), StoreB(dimen) + + error = .true. + NegHDK = 1 + SortX = .FALSE. + +! Unbekanntennumerierung + + DO I = 1, dimen + XNr(I) = I + END DO + +! Genauigkeitsschranke und Bestimmung des groessten Pivotelementes + + PivotWert = ABS(A(1,1)) + PivotZeile = 1 + PivotSpalte = 1 + + DO I = 1, dimen + DO J = 1, dimen + AbsA = ABS(A(I,J)) + IF (AbsA .GT. PivotWert) THEN + PivotWert = AbsA + PivotZeile = I + PivotSpalte = J + END IF + END DO + END DO + + IF (PivotWert .LT. 0.0000001) RETURN ! Pivotelement = 0? + + EpsAbs = PivotWert * 0.1_pReal ** PRECISION(1.0_pReal) + +! V O R W A E R T S T R I A N G U L A T I O N + + DO I = 1, dimen - 1 +! Zeilentausch? + IF (PivotZeile .NE. I) THEN + StoreA(I:dimen) = A(I,I:dimen) + A(I,I:dimen) = A(PivotZeile,I:dimen) + A(PivotZeile,I:dimen) = StoreA(I:dimen) + StoreB(1:dimen) = B(I,1:dimen) + B(I,1:dimen) = B(PivotZeile,1:dimen) + B(PivotZeile,1:dimen) = StoreB(1:dimen) + SortX = .TRUE. + END IF +! Spaltentausch? + IF (PivotSpalte .NE. I) THEN + StoreA(1:dimen) = A(1:dimen,I) + A(1:dimen,I) = A(1:dimen,PivotSpalte) + A(1:dimen,PivotSpalte) = StoreA(1:dimen) + StoreI = XNr(I) + XNr(I) = XNr(PivotSpalte) + XNr(PivotSpalte) = StoreI + SortX = .TRUE. + END IF +! Triangulation + DO J = I + 1, dimen + Quote = A(J,I) / A(I,I) + DO K = I + 1, dimen + A(J,K) = A(J,K) - Quote * A(I,K) + END DO + DO K = 1, dimen + B(J,K) = B(J,K) - Quote * B(I,K) + END DO + END DO +! Bestimmung des groessten Pivotelementes + IP1 = I + 1 + PivotWert = ABS(A(IP1,IP1)) + PivotZeile = IP1 + PivotSpalte = IP1 + DO J = IP1, dimen + DO K = IP1, dimen + AbsA = ABS(A(J,K)) + IF (AbsA .GT. PivotWert) THEN + PivotWert = AbsA + PivotZeile = J + PivotSpalte = K + END IF + END DO + END DO + + IF (PivotWert .LT. EpsAbs) RETURN ! Pivotelement = 0? + + END DO + +! R U E C K W A E R T S A U F L O E S U N G + + DO I = dimen, 1, -1 + DO L = 1, dimen + DO J = I + 1, dimen + B(I,L) = B(I,L) - A(I,J) * B(J,L) + END DO + B(I,L) = B(I,L) / A(I,I) + END DO + END DO + +! Sortieren der Unbekanntenvektoren? + + IF (SortX) THEN + DO L = 1, dimen + StoreA(1:dimen) = B(1:dimen,L) + DO I = 1, dimen + J = XNr(I) + B(J,L) = StoreA(I) + END DO + END DO + END IF + +! Determinante + + LogAbsDetA = 0.0_pReal + NegHDK = 0 + + DO I = 1, dimen + IF (A(I,I) .LT. 0.0_pReal) NegHDK = NegHDK + 1 + AbsA = ABS(A(I,I)) + LogAbsDetA = LogAbsDetA + LOG10(AbsA) + END DO + + error = .false. + + RETURN + + END SUBROUTINE Gauss + + + +!******************************************************************** +! symmetrize a 3x3 matrix +!******************************************************************** + FUNCTION math_symmetric3x3(m) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3) :: math_symmetric3x3,m + integer(pInt) i,j + + forall (i=1:3,j=1:3) math_symmetric3x3(i,j) = 1.0_pReal/2.0_pReal * & + (m(i,j) + m(j,i)) + + END FUNCTION + + +!******************************************************************** +! symmetrize a 6x6 matrix +!******************************************************************** + FUNCTION math_symmetric6x6(m) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(6,6) :: math_symmetric6x6,m + integer(pInt) i,j + + forall (i=1:6,j=1:6) math_symmetric6x6(i,j) = 1.0_pReal/2.0_pReal * & + (m(i,j) + m(j,i)) + + END FUNCTION + - logical error - integer (pInt) NegHDK - real(pReal) LogAbsDetA - real(pReal) A(6,6), B(6,6) - - INTENT (OUT) LogAbsDetA, NegHDK, error - INTENT (INOUT) A, B - - LOGICAL SortX - integer (pInt) PivotZeile, PivotSpalte, StoreI, I, IP1, J, K, L - integer (pInt) XNr(6) - real(pReal) AbsA, PivotWert, EpsAbs, Quote - real(pReal) StoreA(6), StoreB(6) - - error = .true. - NegHDK = 1 - SortX = .FALSE. - -! Unbekanntennumerierung - - DO I = 1, 6 - XNr(I) = I - END DO - -! Genauigkeitsschranke und Bestimmung des groessten Pivotelementes - - PivotWert = ABS(A(1,1)) - PivotZeile = 1 - PivotSpalte = 1 - - DO I = 1, 6 - DO J = 1, 6 - AbsA = ABS(A(I,J)) - IF (AbsA .GT. PivotWert) THEN - PivotWert = AbsA - PivotZeile = I - PivotSpalte = J - END IF - END DO - END DO - - IF (PivotWert .LT. 0.0000001) RETURN ! Pivotelement = 0? - - EpsAbs = PivotWert * 0.1_pReal ** PRECISION(1.0_pReal) - -! V O R W A E R T S T R I A N G U L A T I O N - - DO I = 1, 6 - 1 - -! Zeilentausch? - - IF (PivotZeile .NE. I) THEN - - StoreA(I:6) = A(I,I:6) - A(I,I:6) = A(PivotZeile,I:6) - A(PivotZeile,I:6) = StoreA(I:6) - - StoreB(1:6) = B(I,1:6) - B(I,1:6) = B(PivotZeile,1:6) - B(PivotZeile,1:6) = StoreB(1:6) - - SortX = .TRUE. - - END IF - -! Spaltentausch? - - IF (PivotSpalte .NE. I) THEN - - StoreA(1:6) = A(1:6,I) - A(1:6,I) = A(1:6,PivotSpalte) - A(1:6,PivotSpalte) = StoreA(1:6) - - StoreI = XNr(I) - XNr(I) = XNr(PivotSpalte) - XNr(PivotSpalte) = StoreI - - SortX = .TRUE. - - END IF - -! Triangulation - - DO J = I + 1, 6 - Quote = A(J,I) / A(I,I) - DO K = I + 1, 6 - A(J,K) = A(J,K) - Quote * A(I,K) - END DO - DO K = 1, 6 - B(J,K) = B(J,K) - Quote * B(I,K) - END DO - END DO - -! Bestimmung des groessten Pivotelementes - - IP1 = I + 1 - PivotWert = ABS(A(IP1,IP1)) - PivotZeile = IP1 - PivotSpalte = IP1 - - DO J = IP1, 6 - DO K = IP1, 6 - AbsA = ABS(A(J,K)) - IF (AbsA .GT. PivotWert) THEN - PivotWert = AbsA - PivotZeile = J - PivotSpalte = K - END IF - END DO - END DO - - IF (PivotWert .LT. EpsAbs) RETURN ! Pivotelement = 0? - - END DO - -! R U E C K W A E R T S A U F L O E S U N G - - DO I = 6, 1, -1 - DO L = 1, 6 - DO J = I + 1, 6 - B(I,L) = B(I,L) - A(I,J) * B(J,L) - END DO - B(I,L) = B(I,L) / A(I,I) - END DO - END DO - -! Sortieren der Unbekanntenvektoren? - - IF (SortX) THEN - - DO L = 1, 6 - StoreA(1:6) = B(1:6,L) - DO I = 1, 6 - J = XNr(I) - B(J,L) = StoreA(I) - END DO - END DO - END IF - -! Determinante - - LogAbsDetA = 0.0_pReal - NegHDK = 0 - - DO I = 1, 6 - IF (A(I,I) .LT. 0.0_pReal) NegHDK = NegHDK + 1 - AbsA = ABS(A(I,I)) - LogAbsDetA = LogAbsDetA + LOG10(AbsA) - END DO - - error = .false. - - RETURN - - END SUBROUTINE Gauss - - !******************************************************************** ! determinant of a 3x3 matrix !******************************************************************** @@ -453,6 +492,42 @@ END FUNCTION +!******************************************************************** +! convert 3x3 matrix into vector 9x1 +!******************************************************************** + FUNCTION math_Plain33to9(m33) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3) :: m33 + real(pReal), dimension(9) :: math_Plain33to9 + integer(pInt) i + + forall (i=1:9) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) + return + + END FUNCTION + + +!******************************************************************** +! convert Plain 9x1 back to 3x3 matrix +!******************************************************************** + FUNCTION math_Plain9to33(v9) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(9) :: v9 + real(pReal), dimension(3,3) :: math_Plain9to33 + integer(pInt) i + + forall (i=1:9) math_Plain9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) + return + + END FUNCTION + + !******************************************************************** ! convert symmetric 3x3 matrix into Mandel vector 6x1 !******************************************************************** @@ -492,6 +567,25 @@ END FUNCTION +!******************************************************************** +! convert 3x3x3x3 tensor into plain matrix 9x9 +!******************************************************************** + FUNCTION math_Plain3333to99(m3333) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3,3,3) :: m3333 + real(pReal), dimension(9,9) :: math_Plain3333to99 + integer(pInt) i,j + + forall (i=1:9,j=1:9) math_Plain3333to99(i,j) = & + m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) + return + + END FUNCTION + + !******************************************************************** ! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6 !******************************************************************** @@ -511,51 +605,53 @@ END FUNCTION -!******************************************************************** -! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - FUNCTION math_Mandel66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6) :: m66 - real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 - integer(pInt) i,j - - forall (i=1:6,j=1:6) - math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) - end forall - return - - END FUNCTION - - -!******************************************************************** -! convert Voigt matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - FUNCTION math_Voigt66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6) :: m66 - real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 - integer(pInt) i,j - - forall (i=1:6,j=1:6) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) - end forall - return - - END FUNCTION - +!******************************************************************** +! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor +!******************************************************************** + FUNCTION math_Mandel66to3333(m66) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(6,6) :: m66 + real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 + integer(pInt) i,j + + forall (i=1:6,j=1:6) + math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) + math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) + math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) + math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = invnrmMandel(i)*invnrmMandel(j)*m66(i,j) + end forall + return + + END FUNCTION + + + +!******************************************************************** +! convert Voigt matrix 6x6 back to symmetric 3x3x3x3 tensor +!******************************************************************** + FUNCTION math_Voigt66to3333(m66) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(6,6) :: m66 + real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 + integer(pInt) i,j + + forall (i=1:6,j=1:6) + math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) + math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) + math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) + math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) + end forall + return + + END FUNCTION + + !******************************************************************** ! Euler angles from orientation matrix @@ -725,11 +821,16 @@ real(pReal), dimension(5) :: rnd real(pReal) noise,scatter,cosScatter integer(pInt) i - -if (noise==0.0) then - math_sampleGaussOri = center - return -endif + + +if (noise==0.0) then + + math_sampleGaussOri = center + + return + +endif + ! Helming uses different distribution with Bessel functions ! therefore the gauss scatter width has to be scaled differently @@ -822,48 +923,57 @@ endif return END FUNCTION - - -!******************************************************************** -! symmetric Euler angles for given symmetry string -! 'triclinic' or '', 'monoclinic', 'orthotropic' -!******************************************************************** - PURE FUNCTION math_symmetricEulers(sym,Euler) - - use prec, only: pReal, pInt - implicit none - - character(len=80), intent(in) :: sym - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_symmetricEulers - integer(pInt) i,j - - math_symmetricEulers(1,1) = pi+Euler(1) - math_symmetricEulers(2,1) = Euler(2) - math_symmetricEulers(3,1) = Euler(3) - - math_symmetricEulers(1,2) = pi-Euler(1) - math_symmetricEulers(2,2) = pi-Euler(2) - math_symmetricEulers(3,2) = pi+Euler(3) - - math_symmetricEulers(1,3) = 2.0_pReal*pi-Euler(1) - math_symmetricEulers(2,3) = pi-Euler(2) - math_symmetricEulers(3,3) = pi+Euler(3) - - forall (i=1:3,j=1:3) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi) - - select case (sym) - case ('orthotropic') ! all done - - case ('monoclinic') ! return only first - math_symmetricEulers(:,2:3) = 0.0_pReal - case default ! return blank - math_symmetricEulers = 0.0_pReal - end select - return - + + + + +!******************************************************************** + +! symmetric Euler angles for given symmetry string + +! 'triclinic' or '', 'monoclinic', 'orthotropic' + +!******************************************************************** + + PURE FUNCTION math_symmetricEulers(sym,Euler) + + use prec, only: pReal, pInt + implicit none + + character(len=80), intent(in) :: sym + real(pReal), dimension(3), intent(in) :: Euler + real(pReal), dimension(3,3) :: math_symmetricEulers + integer(pInt) i,j + + math_symmetricEulers(1,1) = pi+Euler(1) + math_symmetricEulers(2,1) = Euler(2) + math_symmetricEulers(3,1) = Euler(3) + + math_symmetricEulers(1,2) = pi-Euler(1) + math_symmetricEulers(2,2) = pi-Euler(2) + math_symmetricEulers(3,2) = pi+Euler(3) + + math_symmetricEulers(1,3) = 2.0_pReal*pi-Euler(1) + math_symmetricEulers(2,3) = pi-Euler(2) + math_symmetricEulers(3,3) = pi+Euler(3) + + forall (i=1:3,j=1:3) math_symmetricEulers(j,i) = modulo(math_symmetricEulers(j,i),2.0_pReal*pi) + + select case (sym) + case ('orthotropic') ! all done + + case ('monoclinic') ! return only first + math_symmetricEulers(:,2:3) = 0.0_pReal + + case default ! return blank + math_symmetricEulers = 0.0_pReal + end select + + return + END FUNCTION - + + !**************************************************************** subroutine math_pDecomposition(FE,U,R,error)