From 9241c7de91b69c6a4b3567339cbec54f545206a4 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 28 Mar 2007 07:21:47 +0000 Subject: [PATCH] added Mandel notation transformations for sym 3x3x3x3 and 3x3 tensors --- trunk/math.f90 | 179 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 145 insertions(+), 34 deletions(-) diff --git a/trunk/math.f90 b/trunk/math.f90 index 871ec48af..7957705ce 100644 --- a/trunk/math.f90 +++ b/trunk/math.f90 @@ -19,6 +19,24 @@ CONTAINS +! *** Initialize random number generator *** +! *** for later use in mpie_fiber and mpie_disturbOri *** + + SUBROUTINE math_init () + + use prec, only: pReal,pInt + implicit none + + integer (pInt) seed + + call random_seed() + call get_seed(seed) + call halton_seed_set(seed) + call halton_ndim_set(3) + + END SUBROUTINE + + SUBROUTINE math_invert3x3(A, InvA, DetA, err) ! Bestimmung der Determinanten und Inversen einer 3x3-Matrix @@ -219,17 +237,13 @@ ! 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 @@ -257,17 +271,12 @@ ! 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? @@ -275,16 +284,12 @@ 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 @@ -305,27 +310,6 @@ END SUBROUTINE Gauss - - - -! *** Initialize random number generator *** -! *** for later use in mpie_fiber and mpie_disturbOri *** - - SUBROUTINE math_init () - - use prec, only: pReal,pInt - implicit none - - integer (pInt) seed - - call random_seed() - call get_seed(seed) - call halton_seed_set(seed) - call halton_ndim_set(3) - - END SUBROUTINE - - !******************************************************************** ! calculate v6 stress from M6x6 stiffness and v6 strain !******************************************************************** @@ -363,6 +347,133 @@ END FUNCTION +!******************************************************************** +! convert a symmetric 3x3 matrix into Mandel vector 6x1 +!******************************************************************** + FUNCTION math_Mandel33to6(m33) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3) :: m33 + real(pReal), dimension(6) :: math_Mandel33to6 + integer(pInt) i + real(pReal), dimension(6), parameter :: nrm = & + (/1.0_pReal,1.0_pReal,1.0_pReal,dsqrt(2.0_pReal),dsqrt(2.0_pReal),dsqrt(2.0_pReal)/) + integer(pInt), dimension (2,6), parameter :: map = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + /),(/2,6/)) + + forall (i=1:6) math_Mandel33to6(i) = nrm(i)*m33(map(1,i),map(2,i)) + return + + END FUNCTION + + +!******************************************************************** +! convert Mandel 6x1 back to symmetric 3x3 matrix +!******************************************************************** + FUNCTION math_Mandel6to33(v6) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(6) :: v6 + real(pReal), dimension(3,3) :: math_Mandel6to33 + integer(pInt) i,j + real(pReal), dimension(6), parameter :: nrm = & + (/1.0_pReal,1.0_pReal,1.0_pReal,dsqrt(0.5_pReal),dsqrt(0.5_pReal),dsqrt(0.5_pReal)/) + integer(pInt), dimension (2,6), parameter :: map = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + /),(/2,6/)) + + forall (i=1:6) + math_Mandel6to33(map(1,i),map(2,i)) = nrm(i)*v6(i) + math_Mandel6to33(map(2,i),map(1,i)) = nrm(i)*v6(i) + end forall + return + + END FUNCTION + + +!******************************************************************** +! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6 +!******************************************************************** + FUNCTION math_Mandel3333to66(m3333) + + use prec, only: pReal,pInt + implicit none + + real(pReal), dimension(3,3,3,3) :: m3333 + real(pReal), dimension(6,6) :: math_Mandel3333to66 + integer(pInt) i,j + real(pReal), dimension(6), parameter :: nrm = & + (/1.0_pReal,1.0_pReal,1.0_pReal,dsqrt(2.0_pReal),dsqrt(2.0_pReal),dsqrt(2.0_pReal)/) + integer(pInt), dimension (2,6), parameter :: map = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + /),(/2,6/)) + + forall (i=1:6,j=1:6) math_Mandel3333to66(i,j) = & + nrm(i)*nrm(j)*m3333(map(1,i),map(2,i),map(1,j),map(2,j)) + 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 + real(pReal), dimension(6), parameter :: nrm = & + (/1.0_pReal,1.0_pReal,1.0_pReal,dsqrt(0.5_pReal),dsqrt(0.5_pReal),dsqrt(0.5_pReal)/) + integer(pInt), dimension (2,6), parameter :: map = & + reshape((/& + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & + /),(/2,6/)) + + forall (i=1:6,j=1:6) + math_Mandel66to3333(map(1,i),map(2,i),map(1,j),map(2,j)) = nrm(i)*nrm(j)*m66(i,j) + math_Mandel66to3333(map(2,i),map(1,i),map(1,j),map(2,j)) = nrm(i)*nrm(j)*m66(i,j) + math_Mandel66to3333(map(1,i),map(2,i),map(2,j),map(1,j)) = nrm(i)*nrm(j)*m66(i,j) + math_Mandel66to3333(map(2,i),map(1,i),map(2,j),map(1,j)) = nrm(i)*nrm(j)*m66(i,j) + end forall + return + + END FUNCTION + + + + !******************************************************************** ! convert a symmetric 3,3 matrix into an array of 6 !********************************************************************