From a6d1e5c911eb774d207715e72e255643684438f6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Feb 2011 12:49:18 +0000 Subject: [PATCH] now deleted old _single files --- code/math_single.f90 | 2888 ------------------------------------ code/mesh.f90 | 4 +- code/mesh_single.f90 | 3360 ------------------------------------------ 3 files changed, 2 insertions(+), 6250 deletions(-) delete mode 100644 code/math_single.f90 delete mode 100644 code/mesh_single.f90 diff --git a/code/math_single.f90 b/code/math_single.f90 deleted file mode 100644 index dfc228b36..000000000 --- a/code/math_single.f90 +++ /dev/null @@ -1,2888 +0,0 @@ -!* $Id: math.f90 748 2011-02-04 15:41:32Z MPIE\c.kords $ -!############################################################## - MODULE math -!############################################################## - - - use prec, only: pReal,pInt - implicit none - - real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal), parameter :: inDeg = 180.0_pReal/pi - real(pReal), parameter :: inRad = pi/180.0_pReal - real(pReal), parameter :: NaN = 0.0_pReal/0.0_pReal ! Not a number -! *** 3x3 Identity *** - real(pReal), dimension(3,3), parameter :: math_I3 = & - reshape( (/ & - 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/) - -! *** 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/)) - - - -! Symmetry operations as quaternions -! 24 for cubic, 12 for hexagonal = 36 -integer(pInt), dimension(2), parameter :: math_NsymOperations = (/24,12/) -real(pReal), dimension(4,36), parameter :: math_symOperations = & - reshape((/& - 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations - 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, & ! 2-fold symmetry - 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, & - 0.0_pReal, -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, 0.0_pReal, & - 0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & ! 3-fold symmetry - -0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & - 0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & - -0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & - 0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & - -0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & - 0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & - -0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & - 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & ! 4-fold symmetry - 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & - -0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & - 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & - -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & - 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! hexagonal symmetry operations - 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & ! 2-fold symmetry - 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.5_pReal, 0.866025403784439_pReal, 0.0_pReal, & - 0.0_pReal, -0.5_pReal, 0.866025403784439_pReal, 0.0_pReal, & - 0.0_pReal, 0.866025403784439_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, -0.866025403784439_pReal, 0.5_pReal, 0.0_pReal, & - 0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & ! 6-fold symmetry - -0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & - 0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & - -0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & - 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & - /),(/4,36/)) - - - - CONTAINS - -!************************************************************************** -! initialization of module -!************************************************************************** - SUBROUTINE math_init () - - use prec, only: pReal,pInt - use numerics, only: fixedSeed - use IO, only: IO_error - implicit none - - real(pReal), dimension(3,3) :: R,R2 - real(pReal), dimension(3) :: Eulers - real(pReal), dimension(4) :: q,q2,axisangle - real(pReal), dimension(2) :: rnd - integer(pInt), dimension(1) :: randInit - - - write(6,*) - write(6,*) '<<<+- math init -+>>>' - write(6,*) '$Id: math.f90 748 2011-02-04 15:41:32Z MPIE\c.kords $' - write(6,*) - - if (fixedSeed > 0_pInt) then - randInit = fixedSeed - call random_seed(put=randInit) - else - call random_seed() - endif - - call random_seed(get=randInit) - write(6,*) 'random seed: ',randInit(1) - write(6,*) - - call halton_seed_set(randInit(1)) - call halton_ndim_set(3) - - ! --- check rotation dictionary --- - - ! +++ q -> a -> q +++ - q = math_qRnd(); - axisangle = math_QuaternionToAxisAngle(q); - q2 = math_AxisAngleToQuaternion(axisangle(1:3),axisangle(4)) - if ( any(abs( q-q2) > 1.0e-8 ) .and. & - any(abs(-q-q2) > 1.0e-8 ) ) & - call IO_error(670) - - ! +++ q -> R -> q +++ - R = math_QuaternionToR(q); - q2 = math_RToQuaternion(R) - if ( any(abs( q-q2) > 1.0e-8 ) .and. & - any(abs(-q-q2) > 1.0e-8 ) ) & - call IO_error(671) - - ! +++ q -> euler -> q +++ - Eulers = math_QuaternionToEuler(q); - q2 = math_EulerToQuaternion(Eulers) - if ( any(abs( q-q2) > 1.0e-8 ) .and. & - any(abs(-q-q2) > 1.0e-8 ) ) & - call IO_error(672) - - ! +++ R -> euler -> R +++ - Eulers = math_RToEuler(R); - R2 = math_EulerToR(Eulers) - if ( any(abs( R-R2) > 1.0e-8 ) ) & - call IO_error(673) - - - ENDSUBROUTINE - - - -!************************************************************************** -! Quicksort algorithm for two-dimensional integer arrays -! -! Sorting is done with respect to array(1,:) -! and keeps array(2:N,:) linked to it. -!************************************************************************** - RECURSIVE SUBROUTINE qsort(a, istart, iend) - - implicit none - integer(pInt), dimension(:,:) :: a - integer(pInt) :: istart,iend,ipivot - - if (istart < iend) then - ipivot = math_partition(a,istart, iend) - call qsort(a, istart, ipivot-1) - call qsort(a, ipivot+1, iend) - endif - return - - ENDSUBROUTINE - -!************************************************************************** -! Partitioning required for quicksort -!************************************************************************** - integer(pInt) function math_partition(a, istart, iend) - - implicit none - integer(pInt), dimension(:,:) :: a - 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 - - j = iend - x = a(1,istart) - do -! find the first element on the right side less than or equal to the pivot point - do j = j, istart, -1 - if (a(1,j) <= x) exit - enddo -! find the first element on the left side greater than the pivot point - do i = i, iend - if (a(1,i) > x) exit - enddo - if (i < j ) then ! if the indexes do not cross, exchange values - do k = 1,d - tmp = a(k,i) - a(k,i) = a(k,j) - a(k,j) = tmp - enddo - else ! if they do cross, exchange left value with pivot and return with the partition index - do k = 1,d - tmp = a(k,istart) - a(k,istart) = a(k,j) - a(k,j) = tmp - enddo - math_partition = j - return - endif - enddo - - endfunction - - -!************************************************************************** -! range of integers starting at one -!************************************************************************** - pure function math_range(N) - - use prec, only: pInt - implicit none - - integer(pInt), intent(in) :: N - integer(pInt) i - integer(pInt), dimension(N) :: math_range - - forall (i=1:N) math_range(i) = i - return - - endfunction - -!************************************************************************** -! second rank identity tensor of specified dimension -!************************************************************************** - pure function math_identity2nd(dimen) - - use prec, only: pReal, pInt - implicit none - - integer(pInt), intent(in) :: dimen - integer(pInt) i - real(pReal), dimension(dimen,dimen) :: math_identity2nd - - math_identity2nd = 0.0_pReal - forall (i=1:dimen) math_identity2nd(i,i) = 1.0_pReal - return - - endfunction - -!************************************************************************** -! permutation tensor e_ijk used for computing cross product of two tensors -! e_ijk = 1 if even permutation of ijk -! e_ijk = -1 if odd permutation of ijk -! e_ijk = 0 otherwise -!************************************************************************** - pure function math_civita(i,j,k) ! change its name from math_permut - ! to math_civita <<>> - - use prec, only: pReal, pInt - implicit none - - integer(pInt), intent(in) :: i,j,k - real(pReal) math_civita - - math_civita = 0.0_pReal - if (((i == 1).and.(j == 2).and.(k == 3)) .or. & - ((i == 2).and.(j == 3).and.(k == 1)) .or. & - ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal - if (((i == 1).and.(j == 3).and.(k == 2)) .or. & - ((i == 2).and.(j == 1).and.(k == 3)) .or. & - ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal - return - - endfunction - -!************************************************************************** -! kronecker delta function d_ij -! d_ij = 1 if i = j -! d_ij = 0 otherwise -!************************************************************************** - pure function math_delta(i,j) - - use prec, only: pReal, pInt - implicit none - - integer(pInt), intent (in) :: i,j - real(pReal) math_delta - - math_delta = 0.0_pReal - if (i == j) math_delta = 1.0_pReal - - return - - endfunction - -!************************************************************************** -! fourth rank identity tensor of specified dimension -!************************************************************************** - pure function math_identity4th(dimen) - - use prec, only: pReal, pInt - implicit none - - integer(pInt), intent(in) :: dimen - integer(pInt) i,j,k,l - real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th - - forall (i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) math_identity4th(i,j,k,l) = & - 0.5_pReal*(math_I3(i,k)*math_I3(j,k)+math_I3(i,l)*math_I3(j,k)) - return - - endfunction - - -!************************************************************************** -! vector product a x b -!************************************************************************** - pure function math_vectorproduct(A,B) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: A,B - real(pReal), dimension(3) :: math_vectorproduct - - math_vectorproduct(1) = A(2)*B(3)-A(3)*B(2) - math_vectorproduct(2) = A(3)*B(1)-A(1)*B(3) - math_vectorproduct(3) = A(1)*B(2)-A(2)*B(1) - - return - - endfunction - - -!************************************************************************** -! tensor product a \otimes b -!************************************************************************** - pure function math_tensorproduct(A,B) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: A,B - real(pReal), dimension(3,3) :: math_tensorproduct - integer(pInt) i,j - - forall (i=1:3,j=1:3) math_tensorproduct(i,j) = A(i)*B(j) - - return - - endfunction - - -!************************************************************************** -! matrix multiplication 3x3 = 1 -!************************************************************************** - pure function math_mul3x3(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - real(pReal), dimension(3), intent(in) :: A,B - real(pReal), dimension(3) :: C - real(pReal) math_mul3x3 - - forall (i=1:3) C(i) = A(i)*B(i) - math_mul3x3 = sum(C) - - return - - endfunction - - -!************************************************************************** -! matrix multiplication 6x6 = 1 -!************************************************************************** - pure function math_mul6x6(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - real(pReal), dimension(6), intent(in) :: A,B - real(pReal), dimension(6) :: C - real(pReal) math_mul6x6 - - forall (i=1:6) C(i) = A(i)*B(i) - math_mul6x6 = sum(C) - - return - - endfunction - - -!************************************************************************** -! matrix multiplication 33x33 = 1 (double contraction --> ij * ij) -!************************************************************************** - pure function math_mul33xx33(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(3,3), intent(in) :: A,B - real(pReal), dimension(3,3) :: C - real(pReal) math_mul33xx33 - - forall (i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j) - math_mul33xx33 = sum(C) - return - - endfunction - -!************************************************************************** -! matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij) -!************************************************************************** - pure function math_mul3333xx33(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(3,3,3,3), intent(in) :: A - real(pReal), dimension(3,3), intent(in) :: B - real(pReal), dimension(3,3) :: C,math_mul3333xx33 - - do i = 1,3 - do j = 1,3 - math_mul3333xx33(i,j) = sum(A(i,j,:,:)*B(:,:)) - enddo; enddo - return - - endfunction - - -!************************************************************************** -! matrix multiplication 33x33 = 3x3 -!************************************************************************** - pure function math_mul33x33(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(3,3), intent(in) :: A,B - real(pReal), dimension(3,3) :: math_mul33x33 - - forall (i=1:3,j=1:3) math_mul33x33(i,j) = & - A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) - return - - endfunction - - -!************************************************************************** -! matrix multiplication 66x66 = 6x6 -!************************************************************************** - pure function math_mul66x66(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(6,6), intent(in) :: A,B - real(pReal), dimension(6,6) :: math_mul66x66 - - forall (i=1:6,j=1:6) math_mul66x66(i,j) = & - A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & - A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) - return - - endfunction - - -!************************************************************************** -! matrix multiplication 99x99 = 9x9 -!************************************************************************** - pure function math_mul99x99(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(9,9), intent(in) :: A,B - - real(pReal), dimension(9,9) :: math_mul99x99 - - - forall (i=1:9,j=1:9) math_mul99x99(i,j) = & - A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & - A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) + & - A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j) - return - - endfunction - - -!************************************************************************** -! matrix multiplication 33x3 = 3 -!************************************************************************** - pure function math_mul33x3(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - real(pReal), dimension(3,3), intent(in) :: A - real(pReal), dimension(3), intent(in) :: B - real(pReal), dimension(3) :: math_mul33x3 - - forall (i=1:3) math_mul33x3(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) - return - - endfunction - - !************************************************************************** -! matrix multiplication complex(33) x real(3) = complex(3) -!************************************************************************** - pure function math_mul33x3_complex(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - complex(pReal), dimension(3,3), intent(in) :: A - real(pReal), dimension(3), intent(in) :: B - complex(pReal), dimension(3) :: math_mul33x3_complex - - forall (i=1:3) math_mul33x3_complex(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) - return - - endfunction - - -!************************************************************************** -! matrix multiplication 66x6 = 6 -!************************************************************************** - pure function math_mul66x6(A,B) - - use prec, only: pReal, pInt - implicit none - - integer(pInt) i - real(pReal), dimension(6,6), intent(in) :: A - real(pReal), dimension(6), intent(in) :: B - real(pReal), dimension(6) :: math_mul66x6 - - forall (i=1:6) math_mul66x6(i) = & - A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & - A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) - return - - endfunction - - -!************************************************************************** -! random quaternion -!************************************************************************** - function math_qRnd() - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4) :: math_qRnd - real(pReal), dimension(3) :: rnd - - call halton(3,rnd) - math_qRnd(1) = cos(2.0_pReal*pi*rnd(1))*sqrt(rnd(3)) - math_qRnd(2) = sin(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3)) - math_qRnd(3) = cos(2.0_pReal*pi*rnd(2))*sqrt(1.0_pReal-rnd(3)) - math_qRnd(4) = sin(2.0_pReal*pi*rnd(1))*sqrt(rnd(3)) - - endfunction - - -!************************************************************************** -! quaternion multiplication q1xq2 = q12 -!************************************************************************** - pure function math_qMul(A,B) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: A, B - real(pReal), dimension(4) :: math_qMul - - math_qMul(1) = A(1)*B(1) - A(2)*B(2) - A(3)*B(3) - A(4)*B(4) - math_qMul(2) = A(1)*B(2) + A(2)*B(1) + A(3)*B(4) - A(4)*B(3) - math_qMul(3) = A(1)*B(3) - A(2)*B(4) + A(3)*B(1) + A(4)*B(2) - math_qMul(4) = A(1)*B(4) + A(2)*B(3) - A(3)*B(2) + A(4)*B(1) - - endfunction - - -!************************************************************************** -! quaternion dotproduct -!************************************************************************** - pure function math_qDot(A,B) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: A, B - real(pReal) math_qDot - - math_qDot = A(1)*B(1) + A(2)*B(2) + A(3)*B(3) + A(4)*B(4) - - endfunction - - -!************************************************************************** -! quaternion conjugation -!************************************************************************** - pure function math_qConj(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(4) :: math_qConj - - math_qConj(1) = Q(1) - math_qConj(2:4) = -Q(2:4) - - endfunction - - -!************************************************************************** -! quaternion norm -!************************************************************************** - pure function math_qNorm(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal) math_qNorm - - math_qNorm = sqrt(max(0.0_pReal, Q(1)*Q(1) + Q(2)*Q(2) + Q(3)*Q(3) + Q(4)*Q(4))) - - endfunction - - -!************************************************************************** -! quaternion inversion -!************************************************************************** - pure function math_qInv(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(4) :: math_qInv - real(pReal) squareNorm - - math_qInv = 0.0_pReal - - squareNorm = math_qDot(Q,Q) - if (squareNorm > tiny(squareNorm)) & - math_qInv = math_qConj(Q) / squareNorm - - endfunction - - -!************************************************************************** -! action of a quaternion on a vector (rotate vector v with Q) -!************************************************************************** - pure function math_qRot(Q,v) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3), intent(in) :: v - real(pReal), dimension(3) :: math_qRot - real(pReal), dimension(4,4) :: T - integer(pInt) i, j - - do i = 1,4 - do j = 1,i - T(i,j) = Q(i) * Q(j) - enddo - enddo - - math_qRot(1) = -v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)) - math_qRot(2) = v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)) - math_qRot(3) = v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3)) - - math_qRot = 2.0_pReal * math_qRot + v - - endfunction - - -!************************************************************************** -! transposition of a 3x3 matrix -!************************************************************************** -pure function math_transpose3x3(A) - - use prec, only: pReal,pInt - implicit none - - real(pReal),dimension(3,3),intent(in) :: A - real(pReal),dimension(3,3) :: math_transpose3x3 - integer(pInt) i,j - - forall(i=1:3, j=1:3) math_transpose3x3(i,j) = A(j,i) - return - - endfunction - - -!************************************************************************** -! Cramer inversion of 3x3 matrix (function) -!************************************************************************** - pure function math_inv3x3(A) - -! direct Cramer inversion of matrix A. -! returns all zeroes if not possible, i.e. if det close to zero - - use prec, only: pReal,pInt - implicit none - - real(pReal),dimension(3,3),intent(in) :: A - real(pReal) DetA - - real(pReal),dimension(3,3) :: math_inv3x3 - - math_inv3x3 = 0.0_pReal - - 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 - math_inv3x3(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) / DetA - math_inv3x3(2,1) = ( -A(2,1) * A(3,3) + A(2,3) * A(3,1) ) / DetA - math_inv3x3(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) / DetA - - math_inv3x3(1,2) = ( -A(1,2) * A(3,3) + A(1,3) * A(3,2) ) / DetA - math_inv3x3(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1) ) / DetA - math_inv3x3(3,2) = ( -A(1,1) * A(3,2) + A(1,2) * A(3,1) ) / DetA - - math_inv3x3(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2) ) / DetA - math_inv3x3(2,3) = ( -A(1,1) * A(2,3) + A(1,3) * A(2,1) ) / DetA - math_inv3x3(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA - endif - return - - endfunction - - - -!************************************************************************** -! Cramer inversion of 3x3 matrix (subroutine) -!************************************************************************** - PURE 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 - - 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 <= 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(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1) ) / DetA - - error = .false. - endif - return - - ENDSUBROUTINE - - - -!************************************************************************** -! Gauss elimination to invert 6x6 matrix -!************************************************************************** - PURE 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 - real(pReal) LogAbsDetA - real(pReal),dimension(dimen,dimen) :: B - - InvA = math_identity2nd(dimen) - B = A - CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error) - RETURN - - ENDSUBROUTINE math_invert - - - -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - PURE 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 - ENDDO - -! 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 - ENDIF - ENDDO - ENDDO - - 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. - ENDIF -! 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. - ENDIF -! 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) - ENDDO - DO K = 1, dimen - B(J,K) = B(J,K) - Quote * B(I,K) - ENDDO - ENDDO -! 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 - ENDIF - ENDDO - ENDDO - - IF (PivotWert .LT. EpsAbs) RETURN ! Pivotelement = 0? - - ENDDO - -! 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) - ENDDO - B(I,L) = B(I,L) / A(I,I) - ENDDO - ENDDO - -! 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) - ENDDO - ENDDO - ENDIF - -! 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) - ENDDO - - error = .false. - - RETURN - - ENDSUBROUTINE 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) = 0.5_pReal * (m(i,j) + m(j,i)) - - endfunction - - -!******************************************************************** -! symmetrize a 6x6 matrix -!******************************************************************** - pure function math_symmetric6x6(m) - - use prec, only: pReal,pInt - implicit none - - integer(pInt) i,j - real(pReal), dimension(6,6), intent(in) :: m - real(pReal), dimension(6,6) :: math_symmetric6x6 - - forall (i=1:6,j=1:6) math_symmetric6x6(i,j) = 0.5_pReal * (m(i,j) + m(j,i)) - - endfunction - - -!******************************************************************** -! equivalent scalar quantity of a full strain tensor -!******************************************************************** - pure function math_equivStrain33(m) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3), intent(in) :: m - real(pReal) math_equivStrain33,e11,e22,e33,s12,s23,s31 - - e11 = (2.0_pReal*m(1,1)-m(2,2)-m(3,3))/3.0_pReal - e22 = (2.0_pReal*m(2,2)-m(3,3)-m(1,1))/3.0_pReal - e33 = (2.0_pReal*m(3,3)-m(1,1)-m(2,2))/3.0_pReal - s12 = 2.0_pReal*m(1,2) - s23 = 2.0_pReal*m(2,3) - s31 = 2.0_pReal*m(3,1) - - math_equivStrain33 = 2.0_pReal*(1.50_pReal*(e11**2.0_pReal+e22**2.0_pReal+e33**2.0_pReal) + & - 0.75_pReal*(s12**2.0_pReal+s23**2.0_pReal+s31**2.0_pReal))**(0.5_pReal)/3.0_pReal - return - - endfunction - -!******************************************************************** -! determinant of a 3x3 matrix -!******************************************************************** - pure function math_det3x3(m) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3), intent(in) :: m - real(pReal) math_det3x3 - - math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) & - -m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) & - +m(1,3)*(m(2,1)*m(3,2)-m(2,2)*m(3,1)) - return - - endfunction - - -!******************************************************************** -! euclidic norm of a 3x1 vector -!******************************************************************** - pure function math_norm3(v) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3), intent(in) :: v - real(pReal) math_norm3 - - math_norm3 = sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) - return - - endfunction - - -!******************************************************************** -! convert 3x3 matrix into vector 9x1 -!******************************************************************** - pure function math_Plain33to9(m33) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3), intent(in) :: 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 - - endfunction - - -!******************************************************************** -! convert Plain 9x1 back to 3x3 matrix -!******************************************************************** - pure function math_Plain9to33(v9) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(9), intent(in) :: 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 - - endfunction - - -!******************************************************************** -! convert symmetric 3x3 matrix into Mandel vector 6x1 -!******************************************************************** - pure function math_Mandel33to6(m33) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3), intent(in) :: m33 - real(pReal), dimension(6) :: math_Mandel33to6 - integer(pInt) i - - forall (i=1:6) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i)) - return - - endfunction - - -!******************************************************************** -! convert Mandel 6x1 back to symmetric 3x3 matrix -!******************************************************************** - pure function math_Mandel6to33(v6) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6), intent(in) :: v6 - real(pReal), dimension(3,3) :: math_Mandel6to33 - integer(pInt) i - - forall (i=1:6) - math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i) - math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i) - end forall - return - - endfunction - - -!******************************************************************** -! convert 3x3x3x3 tensor into plain matrix 9x9 -!******************************************************************** - pure function math_Plain3333to99(m3333) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3,3,3), intent(in) :: 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 - - endfunction - -!******************************************************************** -! plain matrix 9x9 into 3x3x3x3 tensor -!******************************************************************** - pure function math_Plain99to3333(m99) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(9,9), intent(in) :: m99 - real(pReal), dimension(3,3,3,3) :: math_Plain99to3333 - integer(pInt) i,j - - forall (i=1:9,j=1:9) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),& - mapPlain(1,j),mapPlain(2,j)) = m99(i,j) - return - - endfunction - -!******************************************************************** -! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6 -!******************************************************************** - pure function math_Mandel3333to66(m3333) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 - real(pReal), dimension(6,6) :: math_Mandel3333to66 - integer(pInt) i,j - - forall (i=1:6,j=1:6) math_Mandel3333to66(i,j) = & - nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) - return - - endfunction - - - -!******************************************************************** -! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - pure function math_Mandel66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6), intent(in) :: 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 - - endfunction - - - -!******************************************************************** -! convert Voigt matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - pure function math_Voigt66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6), intent(in) :: 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 - - endfunction - - - -!******************************************************************** -! Euler angles (in radians) from rotation matrix -!******************************************************************** - pure function math_RtoEuler(R) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension (3,3), intent(in) :: R - real(pReal), dimension(3) :: math_RtoEuler - real(pReal) sqhkl, squvw, sqhk, val - - sqhkl=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)+R(3,3)*R(3,3)) - squvw=sqrt(R(1,1)*R(1,1)+R(2,1)*R(2,1)+R(3,1)*R(3,1)) - sqhk=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)) -! calculate PHI - val=R(3,3)/sqhkl - - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(2) = acos(val) - - if(math_RtoEuler(2) < 1.0e-8_pReal) then -! calculate phi2 - math_RtoEuler(3) = 0.0_pReal -! calculate phi1 - val=R(1,1)/squvw - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(1) = acos(val) - if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - else -! calculate phi2 - val=R(2,3)/sqhk - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(3) = acos(val) - if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) -! calculate phi1 - val=-R(3,2)/sin(math_RtoEuler(2)) - if(val > 1.0_pReal) val = 1.0_pReal - if(val < -1.0_pReal) val = -1.0_pReal - - math_RtoEuler(1) = acos(val) - if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) - end if - return - - endfunction - - -!******************************************************************** -! quaternion (w+ix+jy+kz) from orientation matrix -!******************************************************************** - pure function math_RtoQuaternion(R) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension (3,3), intent(in) :: R - real(pReal), dimension(4) :: absQ,math_RtoQuaternion - integer(pInt), dimension(1) :: largest - - ! math adopted from http://code.google.com/p/mtex/source/browse/trunk/geometry/geometry_tools/mat2quat.m - - math_RtoQuaternion = 0.0_pReal - - absQ(1) = 0.5_pReal * sqrt(1.0_pReal+R(1,1)+R(2,2)+R(3,3)) - absQ(2) = 0.5_pReal * sqrt(1.0_pReal+R(1,1)-R(2,2)-R(3,3)) - absQ(3) = 0.5_pReal * sqrt(1.0_pReal-R(1,1)+R(2,2)-R(3,3)) - absQ(4) = 0.5_pReal * sqrt(1.0_pReal-R(1,1)-R(2,2)+R(3,3)) - - largest = maxloc(absQ) - select case(largest(1)) - case (1) - - math_RtoQuaternion(2) = R(2,3)-R(3,2) - math_RtoQuaternion(3) = R(3,1)-R(1,3) - math_RtoQuaternion(4) = R(1,2)-R(2,1) - - case (2) - math_RtoQuaternion(1) = R(2,3)-R(3,2) - - math_RtoQuaternion(3) = R(1,2)+R(2,1) - math_RtoQuaternion(4) = R(3,1)+R(1,3) - - case (3) - math_RtoQuaternion(1) = R(3,1)-R(1,3) - math_RtoQuaternion(2) = R(1,2)+R(2,1) - - math_RtoQuaternion(4) = R(2,3)+R(3,2) - - case (4) - math_RtoQuaternion (1) = R(1,2)-R(2,1) - math_RtoQuaternion (2) = R(3,1)+R(1,3) - math_RtoQuaternion (3) = R(3,2)+R(2,3) - - end select - - math_RtoQuaternion = math_RtoQuaternion*0.25_pReal/maxval(absQ) - math_RtoQuaternion(largest(1)) = maxval(absQ) - - return - - endfunction - - -!**************************************************************** -! rotation matrix from Euler angles (in radians) -!**************************************************************** - pure function math_EulerToR (Euler) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_EulerToR - real(pReal) c1, c, c2, s1, s, s2 - - C1 = cos(Euler(1)) - C = cos(Euler(2)) - C2 = cos(Euler(3)) - S1 = sin(Euler(1)) - S = sin(Euler(2)) - S2 = sin(Euler(3)) - math_EulerToR(1,1)=C1*C2-S1*S2*C - math_EulerToR(1,2)=S1*C2+C1*S2*C - math_EulerToR(1,3)=S2*S - math_EulerToR(2,1)=-C1*S2-S1*C2*C - math_EulerToR(2,2)=-S1*S2+C1*C2*C - math_EulerToR(2,3)=C2*S - math_EulerToR(3,1)=S1*S - math_EulerToR(3,2)=-C1*S - math_EulerToR(3,3)=C - return - - endfunction - - -!******************************************************************** -! quaternion (w+ix+jy+kz) from 3-1-3 Euler angles (in radians) -!******************************************************************** - pure function math_EulerToQuaternion(eulerangles) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: eulerangles - real(pReal), dimension(4) :: math_EulerToQuaternion - real(pReal), dimension(3) :: halfangles - real(pReal) c, s - - halfangles = 0.5_pReal * eulerangles - - c = cos(halfangles(2)) - s = sin(halfangles(2)) - - math_EulerToQuaternion(1) = cos(halfangles(1)+halfangles(3)) * c - math_EulerToQuaternion(2) = cos(halfangles(1)-halfangles(3)) * s - math_EulerToQuaternion(3) = sin(halfangles(1)-halfangles(3)) * s - math_EulerToQuaternion(4) = sin(halfangles(1)+halfangles(3)) * c - - endfunction - - -!**************************************************************** -! rotation matrix from axis and angle (in radians) -!**************************************************************** - pure function math_AxisAngleToR(axis,omega) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - real(pReal), dimension(3) :: axisNrm - real(pReal), dimension(3,3) :: math_AxisAngleToR - real(pReal) norm,s,c,c1 - integer(pInt) i - - norm = sqrt(math_mul3x3(axis,axis)) - if (norm > 1.0e-8_pReal) then ! non-zero rotation - forall (i=1:3) axisNrm(i) = axis(i)/norm ! normalize axis to be sure - - s = sin(omega) - c = cos(omega) - c1 = 1.0_pReal - c - - ! formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html - ! below is transposed form to get passive rotation - - math_AxisAngleToR(1,1) = c + c1*axisNrm(1)**2 - math_AxisAngleToR(2,1) = -s*axisNrm(3) + c1*axisNrm(1)*axisNrm(2) - math_AxisAngleToR(3,1) = s*axisNrm(2) + c1*axisNrm(1)*axisNrm(3) - - math_AxisAngleToR(1,2) = s*axisNrm(3) + c1*axisNrm(2)*axisNrm(1) - math_AxisAngleToR(2,2) = c + c1*axisNrm(2)**2 - math_AxisAngleToR(3,2) = -s*axisNrm(1) + c1*axisNrm(2)*axisNrm(3) - - math_AxisAngleToR(1,3) = -s*axisNrm(2) + c1*axisNrm(3)*axisNrm(1) - math_AxisAngleToR(2,3) = s*axisNrm(1) + c1*axisNrm(3)*axisNrm(2) - math_AxisAngleToR(3,3) = c + c1*axisNrm(3)**2 - else - math_AxisAngleToR = math_I3 - endif - - return - - endfunction - - -!**************************************************************** -! quaternion (w+ix+jy+kz) from axis and angle (in radians) -!**************************************************************** - pure function math_AxisAngleToQuaternion(axis,omega) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - real(pReal), dimension(3) :: axisNrm - real(pReal), dimension(4) :: math_AxisAngleToQuaternion - real(pReal) s,c,norm - integer(pInt) i - - norm = sqrt(math_mul3x3(axis,axis)) - if (norm > 1.0e-8_pReal) then ! non-zero rotation - forall (i=1:3) axisNrm(i) = axis(i)/norm ! normalize axis to be sure - ! formula taken from http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters - s = sin(omega/2.0_pReal) - c = cos(omega/2.0_pReal) - math_AxisAngleToQuaternion(1) = c - math_AxisAngleToQuaternion(2:4) = s * axisNrm(1:3) - else - math_AxisAngleToQuaternion = (/1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal/) ! no rotation - endif - - return - - endfunction - - -!******************************************************************** -! orientation matrix from quaternion (w+ix+jy+kz) -!******************************************************************** - pure function math_QuaternionToR(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3,3) :: math_QuaternionToR, T,S - real(pReal) w2 - integer(pInt) i, j - - forall (i = 1:3, j = 1:3) & - T(i,j) = Q(i+1) * Q(j+1) - S = reshape( (/0.0_pReal, Q(4), -Q(3), & - -Q(4),0.0_pReal, +Q(2), & - Q(3), -Q(2),0.0_pReal/),(/3,3/)) ! notation is transposed! - - math_QuaternionToR = (2.0_pReal * Q(1)*Q(1) - 1.0_pReal) * math_I3 + & - 2.0_pReal * T - & - 2.0_pReal * Q(1) * S - - return - - endfunction - - -!******************************************************************** -! 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz) -!******************************************************************** - pure function math_QuaternionToEuler(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3) :: math_QuaternionToEuler - - math_QuaternionToEuler(2) = acos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3))) - - if (abs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then - math_QuaternionToEuler(1) = 2.0_pReal*acos(Q(1)) - math_QuaternionToEuler(3) = 0.0_pReal - else - math_QuaternionToEuler(1) = atan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4)) - if (math_QuaternionToEuler(1) < 0.0_pReal) & - math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi - - math_QuaternionToEuler(3) = atan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4)) - if (math_QuaternionToEuler(3) < 0.0_pReal) & - math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi - endif - - if (math_QuaternionToEuler(2) < 0.0_pReal) & - math_QuaternionToEuler(2) = math_QuaternionToEuler(2) + pi - - return - - endfunction - - -!******************************************************************** -! axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) -!******************************************************************** - pure function math_QuaternionToAxisAngle(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal) halfAngle, sinHalfAngle - real(pReal), dimension(4) :: math_QuaternionToAxisAngle - - halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg - sinHalfAngle = sin(halfAngle) - - if (sinHalfAngle <= 1.0e-4_pReal) then ! very small rotation angle? - math_QuaternionToAxisAngle = 0.0_pReal - else - math_QuaternionToAxisAngle(1:3) = Q(2:4)/sinHalfAngle - math_QuaternionToAxisAngle(4) = halfAngle*2.0_pReal - endif - - return - - endfunction - - -!******************************************************************** -! Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) -!******************************************************************** - pure function math_QuaternionToRodrig(Q) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3) :: math_QuaternionToRodrig - - if (Q(1) /= 0.0_pReal) then ! unless rotation by 180 deg - math_QuaternionToRodrig = Q(2:4)/Q(1) - else - math_QuaternionToRodrig = NaN ! Rodrig is unbound for 180 deg... - endif - - return - - endfunction - - -!************************************************************************** -! misorientation angle between two sets of Euler angles -!************************************************************************** - pure function math_EulerMisorientation(EulerA,EulerB) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: EulerA,EulerB - real(pReal), dimension(3,3) :: r - real(pReal) math_EulerMisorientation, tr - - r = math_mul33x33(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA))) - - tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal - math_EulerMisorientation = abs(0.5_pReal*pi-asin(tr)) - return - - endfunction - - -!************************************************************************** -! figures whether unit quat falls into stereographic standard triangle -!************************************************************************** -pure function math_QuaternionInSST(Q, symmetryType) - - use prec, only: pReal, pInt - implicit none - - !*** input variables - real(pReal), dimension(4), intent(in) :: Q ! orientation - integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal - - !*** output variables - logical math_QuaternionInSST - - !*** local variables - real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q - - Rodrig = math_QuaternionToRodrig(Q) - select case (symmetryType) - case (1) - math_QuaternionInSST = Rodrig(1) > Rodrig(2) .and. & - Rodrig(2) > Rodrig(3) .and. & - Rodrig(3) > 0.0_pReal - case (2) - math_QuaternionInSST = Rodrig(1) > sqrt(3.0_pReal)*Rodrig(2) .and. & - Rodrig(2) > 0.0_pReal .and. & - Rodrig(3) > 0.0_pReal - case default - math_QuaternionInSST = .true. - end select - -endfunction - - -!************************************************************************** -! calculates the disorientation for 2 unit quaternions -!************************************************************************** -function math_QuaternionDisorientation(Q1, Q2, symmetryType) - - use prec, only: pReal, pInt - use IO, only: IO_error - implicit none - - !*** input variables - real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation - Q2 ! 2nd orientation - integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal - - !*** output variables - real(pReal), dimension(4) :: math_QuaternionDisorientation ! disorientation - - !*** local variables - real(pReal), dimension(4) :: dQ,dQsymA,mis - integer(pInt) i,j,k,s - - dQ = math_qMul(math_qConj(Q1),Q2) - math_QuaternionDisorientation = dQ - - select case (symmetryType) - case (0) - if (math_QuaternionDisorientation(1) < 0.0_pReal) & - math_QuaternionDisorientation = -math_QuaternionDisorientation ! keep omega within 0 to 180 deg - - case (1,2) - s = sum(math_NsymOperations(1:symmetryType-1)) - do i = 1,2 - dQ = math_qConj(dQ) ! switch order of "from -- to" - do j = 1,math_NsymOperations(symmetryType) ! run through first crystal's symmetries - dQsymA = math_qMul(math_symOperations(:,s+j),dQ) ! apply sym - do k = 1,math_NsymOperations(symmetryType) ! run through 2nd crystal's symmetries - mis = math_qMul(dQsymA,math_symOperations(:,s+k)) ! apply sym - if (mis(1) < 0.0_pReal) & ! want positive angle - mis = -mis - if (mis(1)-math_QuaternionDisorientation(1) > -1e-8_pReal .and. & - math_QuaternionInSST(mis,symmetryType)) & - math_QuaternionDisorientation = mis ! found better one - enddo; enddo; enddo - - case default - call IO_error(550,symmetryType) ! complain about unknown symmetry - end select - - -endfunction - - -!******************************************************************** -! draw a random sample from Euler space -!******************************************************************** - function math_sampleRandomOri() - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleRandomOri, rnd - - call halton(3,rnd) - math_sampleRandomOri(1) = rnd(1)*2.0_pReal*pi - math_sampleRandomOri(2) = acos(2.0_pReal*rnd(2)-1.0_pReal) - math_sampleRandomOri(3) = rnd(3)*2.0_pReal*pi - return - - endfunction - - -!******************************************************************** -! draw a random sample from Gauss component -! with noise (in radians) half-width -!******************************************************************** - function math_sampleGaussOri(center,noise) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleGaussOri, center, disturb - real(pReal), dimension(3), parameter :: origin = (/0.0_pReal,0.0_pReal,0.0_pReal/) - real(pReal), dimension(5) :: rnd - real(pReal) noise,scatter,cosScatter - integer(pInt) i - -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 - scatter = 0.95_pReal * noise - cosScatter = cos(scatter) - - do - call halton(5,rnd) - forall (i=1:3) rnd(i) = 2.0_pReal*rnd(i)-1.0_pReal ! expand 1:3 to range [-1,+1] - disturb(1) = scatter * rnd(1) ! phi1 - disturb(2) = sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)) ! Phi - disturb(3) = scatter * rnd(2) ! phi2 - if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(origin,disturb)/scatter)**2)) exit - enddo - - math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) - - return - - endfunction - - -!******************************************************************** -! draw a random sample from Fiber component -! with noise (in radians) -!******************************************************************** - function math_sampleFiberOri(alpha,beta,noise) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3) :: math_sampleFiberOri, fiberInC,fiberInS,axis - real(pReal), dimension(2) :: alpha,beta, rnd - real(pReal), dimension(3,3) :: oRot,fRot,pRot - real(pReal) noise, scatter, cos2Scatter, angle - integer(pInt), dimension(2,3), parameter :: rotMap = reshape((/2,3, 3,1, 1,2/),(/2,3/)) - integer(pInt) i - -! Helming uses different distribution with Bessel functions -! therefore the gauss scatter width has to be scaled differently - scatter = 0.95_pReal * noise - cos2Scatter = cos(2.0_pReal*scatter) - -! fiber axis in crystal coordinate system - fiberInC(1)=sin(alpha(1))*cos(alpha(2)) - fiberInC(2)=sin(alpha(1))*sin(alpha(2)) - fiberInC(3)=cos(alpha(1)) -! fiber axis in sample coordinate system - fiberInS(1)=sin(beta(1))*cos(beta(2)) - fiberInS(2)=sin(beta(1))*sin(beta(2)) - fiberInS(3)=cos(beta(1)) - -! ---# rotation matrix from sample to crystal system #--- - angle = -acos(dot_product(fiberInC,fiberInS)) - if(angle /= 0.0_pReal) then -! rotation axis between sample and crystal system (cross product) - forall(i=1:3) axis(i) = fiberInC(rotMap(1,i))*fiberInS(rotMap(2,i))-fiberInC(rotMap(2,i))*fiberInS(rotMap(1,i)) - oRot = math_AxisAngleToR(math_vectorproduct(fiberInC,fiberInS),angle) - else - oRot = math_I3 - end if - -! ---# rotation matrix about fiber axis (random angle) #--- - call halton(1,rnd) - fRot = math_AxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi) - -! ---# rotation about random axis perpend to fiber #--- -! random axis pependicular to fiber axis - call halton(2,axis) - if (fiberInS(3) /= 0.0_pReal) then - axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2))/fiberInS(3) - else if(fiberInS(2) /= 0.0_pReal) then - axis(3)=axis(2) - axis(2)=-(axis(1)*fiberInS(1)+axis(3)*fiberInS(3))/fiberInS(2) - else if(fiberInS(1) /= 0.0_pReal) then - axis(3)=axis(1) - axis(1)=-(axis(2)*fiberInS(2)+axis(3)*fiberInS(3))/fiberInS(1) - end if - -! scattered rotation angle - do - call halton(2,rnd) - angle = acos(cos2Scatter+(1.0_pReal-cos2Scatter)*rnd(1)) - if (rnd(2) <= exp(-1.0_pReal*(angle/scatter)**2)) exit - enddo - call halton(1,rnd) - if (rnd(1) <= 0.5) angle = -angle - pRot = math_AxisAngleToR(axis,angle) - -! ---# apply the three rotations #--- - math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))) - - return - - endfunction - - - -!******************************************************************** -! symmetric Euler angles for given symmetry string -! 'triclinic' or '', 'monoclinic', 'orthotropic' -!******************************************************************** - pure function math_symmetricEulers(sym,Euler) - - use prec, only: pReal, pInt - implicit none - - integer(pInt), 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 (4) ! all done - - case (2) ! return only first - math_symmetricEulers(:,2:3) = 0.0_pReal - - case default ! return blank - math_symmetricEulers = 0.0_pReal - end select - - return - - endfunction - - - -!******************************************************************** -! draw a random sample from Gauss variable -!******************************************************************** -function math_sampleGaussVar(meanvalue, stddev, width) - -use prec, only: pReal, pInt -implicit none - -!*** input variables -real(pReal), intent(in) :: meanvalue, & ! meanvalue of gauss distribution - stddev ! standard deviation of gauss distribution -real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation - -!*** output variables -real(pReal) math_sampleGaussVar - -!*** local variables -real(pReal), dimension(2) :: rnd ! random numbers -real(pReal) scatter, & ! normalized scatter around meanvalue - myWidth - -if (stddev == 0.0) then - math_sampleGaussVar = meanvalue - return -endif - -if (present(width)) then - myWidth = width -else - myWidth = 3.0_pReal ! use +-3*sigma as default value for scatter -endif - -do - call halton(2, rnd) - scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) - if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) & ! test if scattered value is drawn - exit -enddo - -math_sampleGaussVar = scatter * stddev - -endfunction - - - -!**************************************************************** - pure subroutine math_pDecomposition(FE,U,R,error) -!-----FE = R.U -!**************************************************************** - use prec, only: pReal, pInt - implicit none - - real(pReal), intent(in) :: FE(3,3) - real(pReal), intent(out) :: R(3,3), U(3,3) - logical, intent(out) :: error - real(pReal) CE(3,3),EW1,EW2,EW3,EB1(3,3),EB2(3,3),EB3(3,3),UI(3,3),det - - error = .false. - ce = math_mul33x33(transpose(FE),FE) - - CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) - U=sqrt(EW1)*EB1+sqrt(EW2)*EB2+sqrt(EW3)*EB3 - call math_invert3x3(U,UI,det,error) - if (.not. error) R = math_mul33x33(FE,UI) - - return - - ENDSUBROUTINE - - -!********************************************************************** - pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) -!**** EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M - - use prec, only: pReal, pInt - implicit none - - real(pReal), intent(in) :: M(3,3) - real(pReal), intent(out) :: EB1(3,3),EB2(3,3),EB3(3,3),EW1,EW2,EW3 - real(pReal) HI1M,HI2M,HI3M,TOL,R,S,T,P,Q,RHO,PHI,Y1,Y2,Y3,D1,D2,D3 - real(pReal) C1,C2,C3,M1(3,3),M2(3,3),M3(3,3),arg - TOL=1.e-14_pReal - CALL math_hi(M,HI1M,HI2M,HI3M) - R=-HI1M - S= HI2M - T=-HI3M - P=S-R**2.0_pReal/3.0_pReal - Q=2.0_pReal/27.0_pReal*R**3.0_pReal-R*S/3.0_pReal+T - EB1=0.0_pReal - EB2=0.0_pReal - EB3=0.0_pReal - IF((ABS(P).LT.TOL).AND.(ABS(Q).LT.TOL))THEN -! DREI GLEICHE EIGENWERTE - EW1=HI1M/3.0_pReal - EW2=EW1 - EW3=EW1 -! this is not really correct, but this way U is calculated -! correctly in PDECOMPOSITION (correct is EB?=I) - EB1(1,1)=1.0_pReal - EB2(2,2)=1.0_pReal - EB3(3,3)=1.0_pReal - ELSE - RHO=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal - arg=-Q/RHO/2.0_pReal - if(arg.GT.1) arg=1 - if(arg.LT.-1) arg=-1 - PHI=acos(arg) - Y1=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal) - Y2=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) - Y3=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) - EW1=Y1-R/3.0_pReal - EW2=Y2-R/3.0_pReal - EW3=Y3-R/3.0_pReal - C1=ABS(EW1-EW2) - C2=ABS(EW2-EW3) - C3=ABS(EW3-EW1) - - IF(C1.LT.TOL) THEN -! EW1 is equal to EW2 - D3=1.0_pReal/(EW3-EW1)/(EW3-EW2) - M1=M-EW1*math_I3 - M2=M-EW2*math_I3 - EB3=math_mul33x33(M1,M2)*D3 - - EB1=math_I3-EB3 -! both EB2 and EW2 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW2=0.0_pReal - ELSE IF(C2.LT.TOL) THEN -! EW2 is equal to EW3 - D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) - M2=M-math_I3*EW2 - M3=M-math_I3*EW3 - EB1=math_mul33x33(M2,M3)*D1 - EB2=math_I3-EB1 -! both EB3 and EW3 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW3=0.0_pReal - ELSE IF(C3.LT.TOL) THEN -! EW1 is equal to EW3 - D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) - M1=M-math_I3*EW1 - M3=M-math_I3*EW3 - EB2=math_mul33x33(M1,M3)*D2 - EB1=math_I3-EB2 -! both EB3 and EW3 are set to zero so that they do not -! contribute to U in PDECOMPOSITION - EW3=0.0_pReal - ELSE -! all three eigenvectors are different - D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) - D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) - D3=1.0_pReal/(EW3-EW1)/(EW3-EW2) - M1=M-EW1*math_I3 - M2=M-EW2*math_I3 - M3=M-EW3*math_I3 - EB1=math_mul33x33(M2,M3)*D1 - EB2=math_mul33x33(M1,M3)*D2 - EB3=math_mul33x33(M1,M2)*D3 - - END IF - END IF - RETURN - ENDSUBROUTINE - - -!********************************************************************** -!**** HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M - - PURE SUBROUTINE math_hi(M,HI1M,HI2M,HI3M) - use prec, only: pReal, pInt - implicit none - - real(pReal), intent(in) :: M(3,3) - real(pReal), intent(out) :: HI1M, HI2M, HI3M - - HI1M=M(1,1)+M(2,2)+M(3,3) - HI2M=HI1M**2/2.0_pReal-(M(1,1)**2+M(2,2)**2+M(3,3)**2)/2.0_pReal-M(1,2)*M(2,1)-M(1,3)*M(3,1)-M(2,3)*M(3,2) - HI3M=math_det3x3(M) -! QUESTION: is 3rd equiv det(M) ?? if yes, use function math_det !agreed on YES - return - - ENDSUBROUTINE - - - SUBROUTINE get_seed(seed) -! -!******************************************************************************* -! -!! GET_SEED returns a seed for the random number generator. -! -! -! Discussion: -! -! The seed depends on the current time, and ought to be (slightly) -! different every millisecond. Once the seed is obtained, a random -! number generator should be called a few times to further process -! the seed. -! -! Modified: -! -! 27 June 2000 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Output, integer SEED, a pseudorandom seed value. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt) seed - real(pReal) temp - character ( len = 10 ) time - character ( len = 8 ) today - integer(pInt) values(8) - character ( len = 5 ) zone - - call date_and_time ( today, time, zone, values ) - - temp = 0.0D+00 - - temp = temp + dble ( values(2) - 1 ) / 11.0D+00 - temp = temp + dble ( values(3) - 1 ) / 30.0D+00 - temp = temp + dble ( values(5) ) / 23.0D+00 - temp = temp + dble ( values(6) ) / 59.0D+00 - temp = temp + dble ( values(7) ) / 59.0D+00 - temp = temp + dble ( values(8) ) / 999.0D+00 - temp = temp / 6.0D+00 - - if ( temp <= 0.0D+00 ) then - temp = 1.0D+00 / 3.0D+00 - else if ( 1.0D+00 <= temp ) then - temp = 2.0D+00 / 3.0D+00 - end if - - seed = int ( dble ( huge ( 1 ) ) * temp , pInt) -! -! Never use a seed of 0 or maximum integer. -! - if ( seed == 0 ) then - seed = 1 - end if - - if ( seed == huge ( 1 ) ) then - seed = seed - 1 - end if - - return - - ENDSUBROUTINE - - - subroutine halton ( ndim, r ) -! -!******************************************************************************* -! -!! HALTON computes the next element in the Halton sequence. -! -! -! Modified: -! -! 09 March 2003 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer NDIM, the dimension of the element. -! -! Output, real R(NDIM), the next element of the current Halton -! sequence. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, ONLY: pReal, pInt - implicit none - - integer(pInt) ndim - - integer(pInt) base(ndim) - real(pReal) r(ndim) - integer(pInt) seed - integer(pInt) value(1) - - call halton_memory ( 'GET', 'SEED', 1, value ) - seed = value(1) - - call halton_memory ( 'GET', 'BASE', ndim, base ) - - call i_to_halton ( seed, base, ndim, r ) - - value(1) = 1 - call halton_memory ( 'INC', 'SEED', 1, value ) - - return - - ENDSUBROUTINE - - - subroutine halton_memory ( action, name, ndim, value ) -! -!******************************************************************************* -! -!! HALTON_MEMORY sets or returns quantities associated with the Halton sequence. -! -! -! Modified: -! -! 09 March 2003 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, character ( len = * ) ACTION, the desired action. -! 'GET' means get the value of a particular quantity. -! 'SET' means set the value of a particular quantity. -! 'INC' means increment the value of a particular quantity. -! (Only the SEED can be incremented.) -! -! Input, character ( len = * ) NAME, the name of the quantity. -! 'BASE' means the Halton base or bases. -! 'NDIM' means the spatial dimension. -! 'SEED' means the current Halton seed. -! -! Input/output, integer NDIM, the dimension of the quantity. -! If ACTION is 'SET' and NAME is 'BASE', then NDIM is input, and -! is the number of entries in VALUE to be put into BASE. -! -! Input/output, integer VALUE(NDIM), contains a value. -! If ACTION is 'SET', then on input, VALUE contains values to be assigned -! to the internal variable. -! If ACTION is 'GET', then on output, VALUE contains the values of -! the specified internal variable. -! If ACTION is 'INC', then on input, VALUE contains the increment to -! be added to the specified internal variable. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - character ( len = * ) action - integer(pInt), allocatable, save :: base(:) - logical, save :: first_call = .true. - integer(pInt) i - character ( len = * ) name - integer(pInt) ndim - integer(pInt), save :: ndim_save = 0 - integer(pInt), save :: seed = 1 - integer(pInt) value(*) - - if ( first_call ) then - ndim_save = 1 - allocate ( base(ndim_save) ) - base(1) = 2 - first_call = .false. - end if -! -! Set -! - if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then - - if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then - - if ( ndim_save /= ndim ) then - deallocate ( base ) - ndim_save = ndim - allocate ( base(ndim_save) ) - end if - - base(1:ndim) = value(1:ndim) - - else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then - - if ( ndim_save /= value(1) ) then - deallocate ( base ) - ndim_save = value(1) - allocate ( base(ndim_save) ) - do i = 1, ndim_save - base(i) = prime ( i ) - enddo - else - ndim_save = value(1) - end if - else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - seed = value(1) - end if -! -! Get -! - else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then - if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then - if ( ndim /= ndim_save ) then - deallocate ( base ) - ndim_save = ndim - allocate ( base(ndim_save) ) - do i = 1, ndim_save - base(i) = prime(i) - enddo - end if - value(1:ndim_save) = base(1:ndim_save) - else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then - value(1) = ndim_save - else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - value(1) = seed - end if -! -! Increment -! - else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then - if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then - seed = seed + value(1) - end if - end if - - return - - ENDSUBROUTINE - - - subroutine halton_ndim_set ( ndim ) -! -!******************************************************************************* -! -!! HALTON_NDIM_SET sets the dimension for a Halton sequence. -! -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer NDIM, the dimension of the Halton vectors. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt) ndim - integer(pInt) value(1) - - value(1) = ndim - call halton_memory ( 'SET', 'NDIM', 1, value ) - - return - - ENDSUBROUTINE - - - subroutine halton_seed_set ( seed ) -! -!******************************************************************************* -! -!! HALTON_SEED_SET sets the "seed" for the Halton sequence. -! -! -! Discussion: -! -! Calling HALTON repeatedly returns the elements of the -! Halton sequence in order, starting with element number 1. -! An internal counter, called SEED, keeps track of the next element -! to return. Each time the routine is called, the SEED-th element -! is computed, and then SEED is incremented by 1. -! -! To restart the Halton sequence, it is only necessary to reset -! SEED to 1. It might also be desirable to reset SEED to some other value. -! This routine allows the user to specify any value of SEED. -! -! The default value of SEED is 1, which restarts the Halton sequence. -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer SEED, the seed for the Halton sequence. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt), parameter :: ndim = 1 - - integer(pInt) seed - integer(pInt) value(ndim) - - value(1) = seed - call halton_memory ( 'SET', 'SEED', ndim, value ) - - return - - ENDSUBROUTINE - - - subroutine i_to_halton ( seed, base, ndim, r ) -! -!******************************************************************************* -! -!! I_TO_HALTON computes an element of a Halton sequence. -! -! -! Reference: -! -! J H Halton, -! On the efficiency of certain quasi-random sequences of points -! in evaluating multi-dimensional integrals, -! Numerische Mathematik, -! Volume 2, pages 84-90, 1960. -! -! Modified: -! -! 26 February 2001 -! -! Author: -! -! John Burkardt -! -! Parameters: -! -! Input, integer SEED, the index of the desired element. -! Only the absolute value of SEED is considered. SEED = 0 is allowed, -! and returns R = 0. -! -! Input, integer BASE(NDIM), the Halton bases, which should be -! distinct prime numbers. This routine only checks that each base -! is greater than 1. -! -! Input, integer NDIM, the dimension of the sequence. -! -! Output, real R(NDIM), the SEED-th element of the Halton sequence -! for the given bases. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, ONLY: pReal, pInt - implicit none - - integer(pInt) ndim - - integer(pInt) base(ndim) - real(pReal) base_inv(ndim) - integer(pInt) digit(ndim) - integer(pInt) i - real(pReal) r(ndim) - integer(pInt) seed - integer(pInt) seed2(ndim) - - seed2(1:ndim) = abs ( seed ) - - r(1:ndim) = 0.0_pReal - - if ( any ( base(1:ndim) <= 1 ) ) then -!$OMP CRITICAL (write2out) - write ( *, '(a)' ) ' ' - write ( *, '(a)' ) 'I_TO_HALTON - Fatal error!' - write ( *, '(a)' ) ' An input base BASE is <= 1!' - do i = 1, ndim - write ( *, '(i6,i6)' ) i, base(i) - enddo - call flush(6) -!$OMP END CRITICAL (write2out) - stop - end if - - base_inv(1:ndim) = 1.0_pReal / real ( base(1:ndim), pReal ) - - do while ( any ( seed2(1:ndim) /= 0 ) ) - digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim) ) - r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal ) * base_inv(1:ndim) - base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal ) - seed2(1:ndim) = seed2(1:ndim) / base(1:ndim) - enddo - - return - - ENDSUBROUTINE - - - function prime ( n ) -! -!******************************************************************************* -! -!! PRIME returns any of the first PRIME_MAX prime numbers. -! -! -! Note: -! -! PRIME_MAX is 1500, and the largest prime stored is 12553. -! -! Modified: -! -! 21 June 2002 -! -! Author: -! -! John Burkardt -! -! Reference: -! -! Milton Abramowitz and Irene Stegun, -! Handbook of Mathematical Functions, -! US Department of Commerce, 1964, pages 870-873. -! -! Daniel Zwillinger, -! CRC Standard Mathematical Tables and Formulae, -! 30th Edition, -! CRC Press, 1996, pages 95-98. -! -! Parameters: -! -! Input, integer N, the index of the desired prime number. -! N = -1 returns PRIME_MAX, the index of the largest prime available. -! N = 0 is legal, returning PRIME = 1. -! It should generally be true that 0 <= N <= PRIME_MAX. -! -! Output, integer PRIME, the N-th prime. If N is out of range, PRIME -! is returned as 0. -! -! Modified: -! -! 29 April 2005 -! -! Author: -! -! Franz Roters -! - use prec, only: pReal, pInt - implicit none - - integer(pInt), parameter :: prime_max = 1500 - - integer(pInt), save :: icall = 0 - integer(pInt) n - integer(pInt), save, dimension ( prime_max ) :: npvec - integer(pInt) prime - - if ( icall == 0 ) then - - icall = 1 - - npvec(1:100) = (/& - 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & - 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & - 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & - 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & - 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & - 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & - 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & - 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & - 467, 479, 487, 491, 499, 503, 509, 521, 523, 541 /) - - npvec(101:200) = (/ & - 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & - 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & - 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & - 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & - 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & - 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & - 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & - 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & - 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & - 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223 /) - - npvec(201:300) = (/ & - 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & - 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & - 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & - 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & - 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & - 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & - 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & - 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & - 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & - 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987 /) - - npvec(301:400) = (/ & - 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & - 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & - 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & - 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & - 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & - 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & - 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & - 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & - 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & - 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741 /) - - npvec(401:500) = (/ & - 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & - 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & - 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & - 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & - 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & - 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & - 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & - 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & - 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & - 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571 /) - - npvec(501:600) = (/ & - 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & - 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & - 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & - 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & - 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & - 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & - 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & - 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & - 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & - 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409 /) - - npvec(601:700) = (/ & - 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & - 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & - 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & - 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & - 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & - 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & - 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & - 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & - 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & - 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279 /) - - npvec(701:800) = (/ & - 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & - 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & - 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & - 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & - 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & - 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & - 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & - 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & - 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & - 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133 /) - - npvec(801:900) = (/ & - 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & - 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & - 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & - 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & - 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & - 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & - 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & - 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & - 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & - 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997 /) - - npvec(901:1000) = (/ & - 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & - 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & - 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & - 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & - 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & - 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & - 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & - 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & - 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & - 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919 /) - - npvec(1001:1100) = (/ & - 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & - 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & - 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & - 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & - 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & - 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & - 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & - 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & - 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & - 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831 /) - - npvec(1101:1200) = (/ & - 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & - 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & - 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & - 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & - 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & - 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & - 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & - 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & - 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & - 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733 /) - - npvec(1201:1300) = (/ & - 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & - 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & - 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, & - 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, & - 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, & - 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, & - 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, & - 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, & - 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, & - 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657 /) - - npvec(1301:1400) = (/ & - 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, & - 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, & - 10861,10867,10883,10889,10891,10903,10909,19037,10939,10949, & - 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, & - 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, & - 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, & - 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, & - 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, & - 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, & - 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657 /) - - npvec(1401:1500) = (/ & - 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, & - 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, & - 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, & - 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, & - 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, & - 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, & - 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, & - 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, & - 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, & - 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553 /) - - end if - - if ( n == -1 ) then - prime = prime_max - else if ( n == 0 ) then - prime = 1 - else if ( n <= prime_max ) then - prime = npvec(n) - else - prime = 0 -!$OMP CRITICAL (write2out) - - write ( 6, '(a)' ) ' ' - write ( 6, '(a)' ) 'PRIME - Fatal error!' - write ( 6, '(a,i6)' ) ' Illegal prime index N = ', n - write ( 6, '(a,i6)' ) ' N must be between 0 and PRIME_MAX =',prime_max - call flush(6) -!$OMP END CRITICAL (write2out) - - stop - end if - - return - - endfunction - -!************************************************************************** -! volume of tetrahedron given by four vertices -!************************************************************************** - pure function math_volTetrahedron(v1,v2,v3,v4) - - use prec, only: pReal - implicit none - - real(pReal) math_volTetrahedron - real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4 - real(pReal), dimension (3,3) :: m - - m(:,1) = v1-v2 - m(:,2) = v2-v3 - m(:,3) = v3-v4 - - math_volTetrahedron = math_det3x3(m)/6.0_pReal - return - - endfunction - - - - END MODULE math diff --git a/code/mesh.f90 b/code/mesh.f90 index 97f23716e..eee8c7db6 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -3040,7 +3040,7 @@ endsubroutine !*********************************************************** subroutine mesh_build_ipVolumes() - use prec, only: pInt + use prec, only: pInt, tol_gravityNodePos use math, only: math_volTetrahedron implicit none @@ -3070,7 +3070,7 @@ endsubroutine do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last if (gravityNode(j)) then ! valid node index do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list - if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < 1.0e-100_pReal)) then ! found duplicate + if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate gravityNode(j) = .false. ! delete first instance gravityNodePos(:,j) = 0.0_pReal exit ! continue with next suspect diff --git a/code/mesh_single.f90 b/code/mesh_single.f90 deleted file mode 100644 index ce4d7311f..000000000 --- a/code/mesh_single.f90 +++ /dev/null @@ -1,3360 +0,0 @@ -!* $Id: mesh.f90 765 2011-02-16 16:23:08Z MPIE\c.kords $ -!############################################################## - MODULE mesh -!############################################################## - - use prec, only: pReal,pInt - implicit none - -! --------------------------- -! _Nelems : total number of elements in mesh -! _NcpElems : total number of CP elements in mesh -! _Nnodes : total number of nodes in mesh -! _maxNnodes : max number of nodes in any CP element -! _maxNips : max number of IPs in any CP element -! _maxNipNeighbors : max number of IP neighbors in any CP element -! _maxNsharedElems : max number of CP elements sharing a node -! -! _element : FEid, type(internal representation), material, texture, node indices -! _node : x,y,z coordinates (initially!) -! _sharedElem : entryCount and list of elements containing node -! -! _mapFEtoCPelem : [sorted FEid, corresponding CPid] -! _mapFEtoCPnode : [sorted FEid, corresponding CPid] -! -! MISSING: these definitions should actually reside in the -! FE-solver specific part (different for MARC/ABAQUS)..! -! Hence, I suggest to prefix with "FE_" -! -! _Nnodes : # nodes in a specific type of element (how we use it) -! _NoriginalNodes : # nodes in a specific type of element (how it is originally defined by marc) -! _Nips : # IPs in a specific type of element -! _NipNeighbors : # IP neighbors in a specific type of element -! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and -! (negative) neighbor faces per own IP in a specific type of element -! _NfaceNodes : # nodes per face in a specific type of element - -! _nodeOnFace : list of node indices on each face of a specific type of element -! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP -! _nodesAtIP : map IP index to two node indices in a specific type of element -! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] -! _NsubNodes : # subnodes required to fully define all IP volumes - -! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype -! --------------------------- - integer(pInt) mesh_Nelems, mesh_NcpElems, mesh_NelemSets, mesh_maxNelemInSet - integer(pInt) mesh_Nmaterials - integer(pInt) mesh_Nnodes, mesh_maxNnodes, mesh_maxNips, mesh_maxNipNeighbors, mesh_maxNsharedElems, mesh_maxNsubNodes - integer(pInt), dimension(2) :: mesh_maxValStateVar = 0_pInt - character(len=64), dimension(:), allocatable :: mesh_nameElemSet, & ! names of elementSet - mesh_nameMaterial, & ! names of material in solid section - mesh_mapMaterial ! name of elementSet for material - integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet - integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode - integer(pInt), dimension(:,:), allocatable :: mesh_element, & - mesh_sharedElem, & - mesh_nodeTwins ! node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) - integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood - - real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element - real(pReal), dimension(:,:), allocatable :: mesh_node, & - mesh_ipVolume ! volume associated with IP - real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP - mesh_ipCenterOfGravity ! center of gravity of IP - real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP - - integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP - integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor - integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent - integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace - - logical :: noPart ! for cases where the ABAQUS input file does not use part/assembly information - logical, dimension(3) :: mesh_periodicSurface ! flag indicating periodic outer surfaces (used for fluxes) - - integer(pInt) :: hypoelasticTableStyle - integer(pInt) :: initialcondTableStyle - integer(pInt), parameter :: FE_Nelemtypes = 9 - integer(pInt), parameter :: FE_maxNnodes = 8 - integer(pInt), parameter :: FE_maxNsubNodes = 56 - integer(pInt), parameter :: FE_maxNips = 27 - integer(pInt), parameter :: FE_maxNipNeighbors = 6 - integer(pInt), parameter :: FE_maxmaxNnodesAtIP = 8 - integer(pInt), parameter :: FE_NipFaceNodes = 4 - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & - (/8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 8, & ! element 21 - 8, & ! element 117 - 8 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NoriginalNodes = & - (/8, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 8, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 20,& ! element 21 - 8, & ! element 117 - 20 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & - (/8, & ! element 7 - 1, & ! element 134 - 4, & ! element 11 - 9, & ! element 27 - 4, & ! element 157 - 6, & ! element 136 - 27,& ! element 21 - 1, & ! element 117 - 8 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & - (/6, & ! element 7 - 4, & ! element 134 - 4, & ! element 11 - 4, & ! element 27 - 6, & ! element 157 - 6, & ! element 136 - 6, & ! element 21 - 6, & ! element 117 - 6 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NsubNodes = & - (/19,& ! element 7 - 0, & ! element 134 - 5, & ! element 11 - 12,& ! element 27 - 0, & ! element 157 - 15,& ! element 136 - 56,& ! element 21 - 0, & ! element 117 - 19 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_NfaceNodes = & - reshape((/& - 4,4,4,4,4,4, & ! element 7 - 3,3,3,3,0,0, & ! element 134 - 2,2,2,2,0,0, & ! element 11 - 2,2,2,2,0,0, & ! element 27 - 3,3,3,3,0,0, & ! element 157 - 3,4,4,4,3,0, & ! element 136 - 4,4,4,4,4,4, & ! element 21 - 4,4,4,4,4,4, & ! element 117 - 4,4,4,4,4,4 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /),(/FE_maxNipNeighbors,FE_Nelemtypes/)) - integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_maxNnodesAtIP = & - (/1, & ! element 7 - 4, & ! element 134 - 1, & ! element 11 - 2, & ! element 27 - 1, & ! element 157 - 1, & ! element 136 - 4, & ! element 21 - 8, & ! element 117 - 1 & ! element 57 (c3d20r == c3d8 --> copy of 7) - /) - integer(pInt), dimension(FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes), parameter :: FE_nodeOnFace = & - reshape((/& - 1,2,3,4 , & ! element 7 - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,0 , & ! element 134 - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,0,0 , & ! element 11 - 2,3,0,0 , & - 3,4,0,0 , & - 4,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,0,0 , & ! element 27 - 2,3,0,0 , & - 3,4,0,0 , & - 4,1,0,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,3,0 , & ! element 157 - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 , & - 1,2,3,0 , & ! element 136 - 1,4,5,2 , & - 2,5,6,3 , & - 1,3,6,4 , & - 4,6,5,0 , & - 0,0,0,0 , & - 1,2,3,4 , & ! element 21 - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,4 , & ! element 117 - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,4 , & ! element 57 (c3d20r == c3d8 --> copy of 7) - 2,1,5,6 , & - 3,2,6,7 , & - 4,3,7,8 , & - 4,1,5,8 , & - 8,7,6,5 & - /),(/FE_NipFaceNodes,FE_maxNipNeighbors,FE_Nelemtypes/)) - - CONTAINS -! --------------------------- -! subroutine mesh_init() -! function mesh_FEtoCPelement(FEid) -! function mesh_build_ipNeighorhood() -! --------------------------- - - -!*********************************************************** -! initialization -!*********************************************************** - subroutine mesh_init (ip,element) - - use mpie_interface - use prec, only: pInt - use IO, only: IO_error,IO_open_InputFile,IO_abaqus_hasNoPart - use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode - - implicit none - - integer(pInt), parameter :: fileUnit = 222 - integer(pInt) e,element,ip - - write(6,*) - write(6,*) '<<<+- mesh init -+>>>' - write(6,*) '$Id: mesh.f90 765 2011-02-16 16:23:08Z MPIE\c.kords $' - write(6,*) - - call mesh_build_FEdata() ! --- get properties of the different types of elements - - if (IO_open_inputFile(fileUnit)) then ! --- parse info from input file... - - select case (FEsolver) - case ('Spectral') - call mesh_spectral_count_nodesAndElements(fileUnit) - call mesh_spectral_count_cpElements() - call mesh_spectral_map_elements() - call mesh_spectral_map_nodes() - call mesh_spectral_count_cpSizes() - call mesh_spectral_build_nodes(fileUnit) - call mesh_spectral_build_elements(fileUnit) - - case ('Marc') - call mesh_marc_get_tableStyles(fileUnit) - call mesh_marc_count_nodesAndElements(fileUnit) - call mesh_marc_count_elementSets(fileUnit) - call mesh_marc_map_elementSets(fileUnit) - call mesh_marc_count_cpElements(fileUnit) - call mesh_marc_map_elements(fileUnit) - call mesh_marc_map_nodes(fileUnit) - call mesh_marc_build_nodes(fileUnit) - call mesh_marc_count_cpSizes(fileunit) - call mesh_marc_build_elements(fileUnit) - call mesh_marc_get_mpieOptions(fileUnit) - case ('Abaqus') - noPart = IO_abaqus_hasNoPart(fileUnit) - call mesh_abaqus_count_nodesAndElements(fileUnit) - call mesh_abaqus_count_elementSets(fileUnit) - call mesh_abaqus_count_materials(fileUnit) - call mesh_abaqus_map_elementSets(fileUnit) - call mesh_abaqus_map_materials(fileUnit) - call mesh_abaqus_count_cpElements(fileUnit) - call mesh_abaqus_map_elements(fileUnit) - call mesh_abaqus_map_nodes(fileUnit) - call mesh_abaqus_build_nodes(fileUnit) - call mesh_abaqus_count_cpSizes(fileunit) - call mesh_abaqus_build_elements(fileUnit) - end select - close (fileUnit) - - call mesh_build_subNodeCoords() - call mesh_build_ipVolumes() - call mesh_build_ipAreas() - call mesh_build_nodeTwins() - call mesh_build_sharedElems() - call mesh_build_ipNeighborhood() - call mesh_tell_statistics() - - parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive - else - call IO_error(101) ! cannot open input file - endif - - FEsolving_execElem = (/1,mesh_NcpElems/) - allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt - forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e)) - - allocate(calcMode(mesh_maxNips,mesh_NcpElems)) - write(6,*) '<<<+- mesh init done -+>>>' - calcMode = .false. ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc" - lastMode = .true. ! and its mode is already known... - endsubroutine - - - -!*********************************************************** -! mapping of FE element types to internal representation -!*********************************************************** - function FE_mapElemtype(what) - - use IO, only: IO_lc - - implicit none - - character(len=*), intent(in) :: what - integer(pInt) FE_mapElemtype - - select case (IO_lc(what)) - case ( '7', & - 'c3d8') - FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick - case ('134', & - 'c3d4') - FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron - case ( '11', & - 'cpe4') - FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain - case ( '27', & - 'cpe8') - FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral - case ('157') - FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations - case ('136', & - 'c3d6') - FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral - case ( '21', & - 'c3d20') - FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted quadratic hexahedral - case ( '117', & - '123', & - 'c3d8r') - FE_mapElemtype = 8 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration - case ( '57', & - 'c3d20r') - FE_mapElemtype = 9 ! Three-dimensional Arbitrarily Distorted quad hexahedral with reduced integration - case default - FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..! - endselect - - endfunction - - - -!*********************************************************** -! FE to CP id mapping by binary search thru lookup array -! -! valid questions are 'elem', 'node' -!*********************************************************** - function mesh_FEasCP(what,id) - - use prec, only: pInt - use IO, only: IO_lc - implicit none - - character(len=*), intent(in) :: what - integer(pInt), intent(in) :: id - integer(pInt), dimension(:,:), pointer :: lookupMap - integer(pInt) mesh_FEasCP, lower,upper,center - - mesh_FEasCP = 0_pInt - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - endselect - - lower = 1_pInt - upper = size(lookupMap,2) - - ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - if (lookupMap(1,lower) == id) then - mesh_FEasCP = lookupMap(2,lower) - return - elseif (lookupMap(1,upper) == id) then - mesh_FEasCP = lookupMap(2,upper) - return - endif - - ! binary search in between bounds - do while (upper-lower > 1) - center = (lower+upper)/2 - if (lookupMap(1,center) < id) then - lower = center - elseif (lookupMap(1,center) > id) then - upper = center - else - mesh_FEasCP = lookupMap(2,center) - exit - endif - enddo - return - - endfunction - - -!*********************************************************** -! find face-matching element of same type -!*********************************************************** -subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace) - -use prec, only: pInt -implicit none - -!*** output variables -integer(pInt), intent(out) :: matchingElem, & ! matching CP element ID - matchingFace ! matching FE face ID - -!*** input variables -integer(pInt), intent(in) :: face, & ! FE face ID - elem ! FE elem ID - -!*** local variables -integer(pInt), dimension(FE_NfaceNodes(face,mesh_element(2,elem))) :: & - myFaceNodes ! global node ids on my face -integer(pInt) myType, & - candidateType, & - candidateElem, & - candidateFace, & - candidateFaceNode, & - minNsharedElems, & - NsharedElems, & - lonelyNode, & - i, & - n, & - dir ! periodicity direction -integer(pInt), dimension(:), allocatable :: element_seen -logical checkTwins - - -minNsharedElems = mesh_maxNsharedElems + 1_pInt ! init to worst case -matchingFace = 0_pInt -matchingElem = 0_pInt ! intialize to "no match found" -myType = mesh_element(2,elem) ! figure elemType - -do n = 1,FE_NfaceNodes(face,myType) ! loop over nodes on face - myFaceNodes(n) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(n,face,myType),elem)) ! CP id of face node - NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node - if (NsharedElems < minNsharedElems) then - minNsharedElems = NsharedElems ! remember min # shared elems - lonelyNode = n ! remember most lonely node - endif -enddo - -allocate(element_seen(minNsharedElems)) -element_seen = 0_pInt - -checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements - candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem - if (all(element_seen /= candidateElem)) then ! element seen for the first time? - element_seen(i) = candidateElem - candidateType = mesh_element(2,candidateElem) ! figure elemType of candidate -checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate - if (FE_NfaceNodes(candidateFace,candidateType) /= FE_NfaceNodes(face,myType) & ! incompatible face - .or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face - cycle checkCandidateFace - endif - checkTwins = .false. - do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face - candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) - if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes - checkTwins = .true. ! perhaps the twin nodes do match - exit - endif - enddo - if(checkTwins) then -checkCandidateFaceTwins: do dir = 1,3 - do n = 1,FE_NfaceNodes(candidateFace,candidateType) ! loop through nodes on face - candidateFaceNode = mesh_FEasCP('node', mesh_element(4+FE_nodeOnFace(n,candidateFace,candidateType),candidateElem)) - if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either - if (dir == 3) then - cycle checkCandidateFace - else - cycle checkCandidateFaceTwins ! try twins in next dimension - endif - endif - enddo - exit checkCandidateFaceTwins - enddo checkCandidateFaceTwins - endif - matchingFace = candidateFace - matchingElem = candidateElem - exit checkCandidate ! found my matching candidate - enddo checkCandidateFace - endif -enddo checkCandidate - -deallocate(element_seen) - -endsubroutine - - -!******************************************************************** -! get properties of different types of finite elements -! -! assign globals: -! FE_nodesAtIP, FE_ipNeighbor, FE_subNodeParent, FE_subNodeOnIPFace -!******************************************************************** - subroutine mesh_build_FEdata () - - use prec, only: pInt - implicit none - - allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Nelemtypes)) ; FE_nodesAtIP = 0_pInt - allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_ipNeighbor = 0_pInt - allocate(FE_subNodeParent(FE_maxNips,FE_maxNsubNodes,FE_Nelemtypes)) ; FE_subNodeParent = 0_pInt - allocate(FE_subNodeOnIPFace(FE_NipFaceNodes,FE_maxNipNeighbors,FE_maxNips,FE_Nelemtypes)) ; FE_subNodeOnIPFace = 0_pInt - - ! fill FE_nodesAtIP with data - FE_nodesAtIP(:,:FE_Nips(1),1) = & ! element 7 - reshape((/& - 1, & - 2, & - 4, & - 3, & - 5, & - 6, & - 8, & - 7 & - /),(/FE_maxNnodesAtIP(1),FE_Nips(1)/)) - FE_nodesAtIP(:,:FE_Nips(2),2) = & ! element 134 - reshape((/& - 1,2,3,4 & - /),(/FE_maxNnodesAtIP(2),FE_Nips(2)/)) - FE_nodesAtIP(:,:FE_Nips(3),3) = & ! element 11 - reshape((/& - 1, & - 2, & - 4, & - 3 & - /),(/FE_maxNnodesAtIP(3),FE_Nips(3)/)) - FE_nodesAtIP(:,:FE_Nips(4),4) = & ! element 27 - reshape((/& - 1,0, & - 1,2, & - 2,0, & - 1,4, & - 0,0, & - 2,3, & - 4,0, & - 3,4, & - 3,0 & - /),(/FE_maxNnodesAtIP(4),FE_Nips(4)/)) - FE_nodesAtIP(:,:FE_Nips(5),5) = & ! element 157 - reshape((/& - 1, & - 2, & - 3, & - 4 & - /),(/FE_maxNnodesAtIP(5),FE_Nips(5)/)) - FE_nodesAtIP(:,:FE_Nips(6),6) = & ! element 136 - reshape((/& - 1, & - 2, & - 3, & - 4, & - 5, & - 6 & - /),(/FE_maxNnodesAtIP(6),FE_Nips(6)/)) - FE_nodesAtIP(:,:FE_Nips(7),7) = & ! element 21 - reshape((/& - 1,0, 0,0, & - 1,2, 0,0, & - 2,0, 0,0, & - 1,4, 0,0, & - 1,3, 2,4, & - 2,3, 0,0, & - 4,0, 0,0, & - 3,4, 0,0, & - 3,0, 0,0, & - 1,5, 0,0, & - 1,6, 2,5, & - 2,6, 0,0, & - 1,8, 4,5, & - 0,0, 0,0, & - 2,7, 3,6, & - 4,8, 0,0, & - 3,8, 4,7, & - 3,7, 0,0, & - 5,0, 0,0, & - 5,6, 0,0, & - 6,0, 0,0, & - 5,8, 0,0, & - 5,7, 6,8, & - 6,7, 0,0, & - 8,0, 0,0, & - 7,8, 0,0, & - 7,0, 0,0 & - /),(/FE_maxNnodesAtIP(7),FE_Nips(7)/)) -! FE_nodesAtIP(:,:FE_Nips(8),8) = & ! element 117 (c3d8r --> single IP per element, so no need for this mapping) -! reshape((/& -! 1,2,3,4,5,6,7,8 & -! /),(/FE_maxNnodesAtIP(8),FE_Nips(8)/)) - FE_nodesAtIP(:,:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& - 1, & - 2, & - 4, & - 3, & - 5, & - 6, & - 8, & - 7 & - /),(/FE_maxNnodesAtIP(9),FE_Nips(9)/)) - - ! fill FE_ipNeighbor with data - FE_ipNeighbor(:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 - reshape((/& - 2,-5, 3,-2, 5,-1, & - -3, 1, 4,-2, 6,-1, & - 4,-5,-4, 1, 7,-1, & - -3, 3,-4, 2, 8,-1, & - 6,-5, 7,-2,-6, 1, & - -3, 5, 8,-2,-6, 2, & - 8,-5,-4, 5,-6, 3, & - -3, 7,-4, 6,-6, 4 & - /),(/FE_NipNeighbors(1),FE_Nips(1)/)) - FE_ipNeighbor(:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 - reshape((/& - -1,-2,-3,-4 & - /),(/FE_NipNeighbors(2),FE_Nips(2)/)) - FE_ipNeighbor(:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 - reshape((/& - 2,-4, 3,-1, & - -2, 1, 4,-1, & - 4,-4,-3, 1, & - -2, 3,-3, 2 & - /),(/FE_NipNeighbors(3),FE_Nips(3)/)) - FE_ipNeighbor(:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 - reshape((/& - 2,-4, 4,-1, & - 3, 1, 5,-1, & - -2, 2, 6,-1, & - 5,-4, 7, 1, & - 6, 4, 8, 2, & - -2, 5, 9, 3, & - 8,-4,-3, 4, & - 9, 7,-3, 5, & - -2, 8,-3, 6 & - /),(/FE_NipNeighbors(4),FE_Nips(4)/)) - FE_ipNeighbor(:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 - reshape((/& - 2,-4, 3,-2, 4,-1, & - 3,-2, 1,-3, 4,-1, & - 1,-3, 2,-4, 4,-1, & - 1,-3, 2,-4, 3,-2 & - /),(/FE_NipNeighbors(5),FE_Nips(5)/)) - FE_ipNeighbor(:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 - reshape((/& - 2,-4, 3,-2, 4,-1, & - -3, 1, 3,-2, 5,-1, & - 2,-4,-3, 1, 6,-1, & - 5,-4, 6,-2,-5, 1, & - -3, 4, 6,-2,-5, 2, & - 5,-4,-3, 4,-5, 3 & - /),(/FE_NipNeighbors(6),FE_Nips(6)/)) - FE_ipNeighbor(:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 - reshape((/& - 2,-5, 4,-2,10,-1, & - 3, 1, 5,-2,11,-1, & - -3, 2, 6,-2,12,-1, & - 5,-5, 7, 1,13,-1, & - 6, 4, 8, 2,14,-1, & - -3, 5, 9, 3,15,-1, & - 8,-5,-4, 4,16,-1, & - 9, 7,-4, 5,17,-1, & - -3, 8,-4, 6,18,-1, & - 11,-5,13,-2,19, 1, & - 12,10,14,-2,20, 2, & - -3,11,15,-2,21, 3, & - 14,-5,16,10,22, 4, & - 15,13,17,11,23, 5, & - -3,14,18,12,24, 6, & - 17,-5,-4,13,25, 7, & - 18,16,-4,14,26, 8, & - -3,17,-4,15,27, 9, & - 20,-5,22,-2,-6,10, & - 21,19,23,-2,-6,11, & - -3,20,24,-2,-6,12, & - 23,-5,25,19,-6,13, & - 24,22,26,20,-6,14, & - -3,23,27,21,-6,15, & - 26,-5,-4,22,-6,16, & - 27,25,-4,23,-6,17, & - -3,26,-4,24,-6,18 & - /),(/FE_NipNeighbors(7),FE_Nips(7)/)) -FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 - reshape((/& - -3,-5,-4,-2,-6,-1 & - /),(/FE_NipNeighbors(8),FE_Nips(8)/)) - FE_ipNeighbor(:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& - 2,-5, 3,-2, 5,-1, & - -3, 1, 4,-2, 6,-1, & - 4,-5,-4, 1, 7,-1, & - -3, 3,-4, 2, 8,-1, & - 6,-5, 7,-2,-6, 1, & - -3, 5, 8,-2,-6, 2, & - 8,-5,-4, 5,-6, 3, & - -3, 7,-4, 6,-6, 4 & - /),(/FE_NipNeighbors(9),FE_Nips(9)/)) - - ! fill FE_subNodeParent with data - FE_subNodeParent(:FE_Nips(1),:FE_NsubNodes(1),1) = & ! element 7 - reshape((/& - 1, 2, 0, 0, 0, 0, 0, 0, & - 2, 3, 0, 0, 0, 0, 0, 0, & - 3, 4, 0, 0, 0, 0, 0, 0, & - 4, 1, 0, 0, 0, 0, 0, 0, & - 1, 5, 0, 0, 0, 0, 0, 0, & - 2, 6, 0, 0, 0, 0, 0, 0, & - 3, 7, 0, 0, 0, 0, 0, 0, & - 4, 8, 0, 0, 0, 0, 0, 0, & - 5, 6, 0, 0, 0, 0, 0, 0, & - 6, 7, 0, 0, 0, 0, 0, 0, & - 7, 8, 0, 0, 0, 0, 0, 0, & - 8, 5, 0, 0, 0, 0, 0, 0, & - 1, 2, 3, 4, 0, 0, 0, 0, & - 1, 2, 6, 5, 0, 0, 0, 0, & - 2, 3, 7, 6, 0, 0, 0, 0, & - 3, 4, 8, 7, 0, 0, 0, 0, & - 1, 4, 8, 5, 0, 0, 0, 0, & - 5, 6, 7, 8, 0, 0, 0, 0, & - 1, 2, 3, 4, 5, 6, 7, 8 & - /),(/FE_Nips(1),FE_NsubNodes(1)/)) -!FE_subNodeParent(:FE_Nips(2),:FE_NsubNodes(2),2) ! element 134 has no subnodes - FE_subNodeParent(:FE_Nips(3),:FE_NsubNodes(3),3) = & ! element 11 - reshape((/& - 1, 2, 0, 0, & - 2, 3, 0, 0, & - 3, 4, 0, 0, & - 4, 1, 0, 0, & - 1, 2, 3, 4 & - /),(/FE_Nips(3),FE_NsubNodes(3)/)) - FE_subNodeParent(:FE_Nips(4),:FE_NsubNodes(4),4) = & ! element 27 - reshape((/& - 1, 1, 2, 0, 0, 0, 0, 0, 0, & - 1, 2, 2, 0, 0, 0, 0, 0, 0, & - 2, 2, 3, 0, 0, 0, 0, 0, 0, & - 2, 3, 3, 0, 0, 0, 0, 0, 0, & - 3, 3, 4, 0, 0, 0, 0, 0, 0, & - 3, 4, 4, 0, 0, 0, 0, 0, 0, & - 4, 4, 1, 0, 0, 0, 0, 0, 0, & - 4, 1, 1, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 2, 2, 4, 4, 3, & - 2, 2, 2, 2, 1, 1, 3, 3, 4, & - 3, 3, 3, 3, 2, 2, 4, 4, 1, & - 4, 4, 4, 4, 1, 1, 3, 3, 2 & - /),(/FE_Nips(4),FE_NsubNodes(4)/)) - !FE_subNodeParent(:FE_Nips(5),:FE_NsubNodes(5),5) = & ! element 157 - ! reshape((/& - ! *still to be defined* - ! /),(/FE_Nips(5),FE_NsubNodes(5)/)) - FE_subNodeParent(:FE_Nips(6),:FE_NsubNodes(6),6) = & ! element 136 - reshape((/& - 1, 2, 0, 0, 0, 0, & - 2, 3, 0, 0, 0, 0, & - 3, 1, 0, 0, 0, 0, & - 1, 4, 0, 0, 0, 0, & - 2, 5, 0, 0, 0, 0, & - 3, 6, 0, 0, 0, 0, & - 4, 5, 0, 0, 0, 0, & - 5, 6, 0, 0, 0, 0, & - 6, 4, 0, 0, 0, 0, & - 1, 2, 3, 0, 0, 0, & - 1, 2, 4, 5, 0, 0, & - 2, 3, 5, 6, 0, 0, & - 1, 3, 4, 6, 0, 0, & - 4, 5, 6, 0, 0, 0, & - 1, 2, 3, 4, 5, 6 & - /),(/FE_Nips(6),FE_NsubNodes(6)/)) - FE_subNodeParent(:FE_Nips(7),:FE_NsubNodes(7),7) = & ! element 21 - reshape((/& - 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 1, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 3, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 5, 6, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 6, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 7, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 8, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 2, 2, 4, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 2, 2, 2, 1, 1, 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 3, 3, 3, 2, 2, 4, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 4, 4, 4, 1, 1, 3, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 2, 2, 5, 5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 2, 2, 2, 1, 1, 6, 6, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 2, 2, 2, 2, 3, 3, 6, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 3, 3, 3, 2, 2, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 3, 3, 3, 3, 4, 4, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 4, 4, 4, 3, 3, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 4, 4, 4, 4, 1, 1, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 4, 4, 5, 5, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 5, 5, 5, 5, 1, 1, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 6, 6, 6, 6, 2, 2, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 6, 6, 6, 6, 2, 2, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 7, 7, 7, 7, 3, 3, 6, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 7, 7, 7, 7, 3, 3, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 8, 8, 8, 8, 4, 4, 7, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 8, 8, 8, 8, 4, 4, 5, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 5, 5, 5, 5, 1, 1, 8, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 5, 5, 5, 5, 6, 6, 8, 8, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 6, 6, 6, 6, 5, 5, 7, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 7, 7, 7, 7, 6, 6, 8, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 8, 8, 8, 8, 5, 5, 7, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 5, 5, 5, 5, 3, 3, 6, 6, 8, 8, 7, & - 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 3, 3, 3, 3, 6, 6, 6, 6, 4, 4, 5, 5, 7, 7, 8, & - 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 4, 4, 4, 4, 7, 7, 7, 7, 1, 1, 6, 6, 8, 8, 5, & - 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1, 3, 3, 3, 3, 8, 8, 8, 8, 2, 2, 5, 5, 7, 7, 6, & - 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 7, 7, 3, & - 6, 6, 6, 6, 6, 6, 6, 6, 2, 2, 2, 2, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 8, 8, 4, & - 7, 7, 7, 7, 7, 7, 7, 7, 3, 3, 3, 3, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 4, 4, 5, 5, 1, & - 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 7, 7, 1, 1, 3, 3, 6, 6, 2 & - /),(/FE_Nips(7),FE_NsubNodes(7)/)) -!FE_subNodeParent(:FE_Nips(8),:FE_NsubNodes(8),8) ! element 117 has no subnodes - FE_subNodeParent(:FE_Nips(9),:FE_NsubNodes(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& - 1, 2, 0, 0, 0, 0, 0, 0, & - 2, 3, 0, 0, 0, 0, 0, 0, & - 3, 4, 0, 0, 0, 0, 0, 0, & - 4, 1, 0, 0, 0, 0, 0, 0, & - 1, 5, 0, 0, 0, 0, 0, 0, & - 2, 6, 0, 0, 0, 0, 0, 0, & - 3, 7, 0, 0, 0, 0, 0, 0, & - 4, 8, 0, 0, 0, 0, 0, 0, & - 5, 6, 0, 0, 0, 0, 0, 0, & - 6, 7, 0, 0, 0, 0, 0, 0, & - 7, 8, 0, 0, 0, 0, 0, 0, & - 8, 5, 0, 0, 0, 0, 0, 0, & - 1, 2, 3, 4, 0, 0, 0, 0, & - 1, 2, 6, 5, 0, 0, 0, 0, & - 2, 3, 7, 6, 0, 0, 0, 0, & - 3, 4, 8, 7, 0, 0, 0, 0, & - 1, 4, 8, 5, 0, 0, 0, 0, & - 5, 6, 7, 8, 0, 0, 0, 0, & - 1, 2, 3, 4, 5, 6, 7, 8 & - /),(/FE_Nips(9),FE_NsubNodes(9)/)) - - ! fill FE_subNodeOnIPFace with data - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(1),:FE_Nips(1),1) = & ! element 7 - reshape((/& - 9,21,27,22, & ! 1 - 1,13,25,12, & - 12,25,27,21, & - 1, 9,22,13, & - 13,22,27,25, & - 1,12,21, 9, & - 2,10,23,14, & ! 2 - 9,22,27,21, & - 10,21,27,23, & - 2,14,22, 9, & - 14,23,27,22, & - 2, 9,21,10, & - 11,24,27,21, & ! 3 - 4,12,25,16, & - 4,16,24,11, & - 12,21,27,25, & - 16,25,27,24, & - 4,11,21,12, & - 3,15,23,10, & ! 4 - 11,21,27,24, & - 3,11,24,15, & - 10,23,27,21, & - 15,24,27,23, & - 3,10,21,11, & - 17,22,27,26, & ! 5 - 5,20,25,13, & - 20,26,27,25, & - 5,13,22,17, & - 5,17,26,20, & - 13,25,27,22, & - 6,14,23,18, & ! 6 - 17,26,27,22, & - 18,23,27,26, & - 6,17,22,14, & - 6,18,26,17, & - 14,22,27,23, & - 19,26,27,24, & ! 7 - 8,16,25,20, & - 8,19,24,16, & - 20,25,27,26, & - 8,20,26,19, & - 16,24,27,25, & - 7,18,23,15, & ! 8 - 19,24,27,26, & - 7,15,24,19, & - 18,26,27,23, & - 7,19,26,18, & - 15,23,27,24 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(1),FE_Nips(1)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(2),:FE_Nips(2),2) = & ! element 134 - reshape((/& - 1, 1, 3, 2, & ! 1 - 1, 1, 2, 4, & - 2, 2, 3, 4, & - 1, 1, 4, 3 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(2),FE_Nips(2)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(3),:FE_Nips(3),3) = & ! element 11 - reshape((/& - 5, 9, 0, 0, & ! 1 - 1, 8, 0, 0, & - 8, 9, 0, 0, & - 1, 5, 0, 0, & - 2, 6, 0, 0, & ! 2 - 5, 9, 0, 0, & - 6, 9, 0, 0, & - 2, 5, 0, 0, & - 3, 6, 0, 0, & ! 3 - 7, 9, 0, 0, & - 3, 7, 0, 0, & - 6, 9, 0, 0, & - 7, 9, 0, 0, & ! 4 - 4, 8, 0, 0, & - 4, 7, 0, 0, & - 8, 9, 0, 0 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(3),FE_Nips(3)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(4),:FE_Nips(4),4) = & ! element 27 - reshape((/& - 9,17, 0, 0, & ! 1 - 1,16, 0, 0, & - 16,17, 0, 0, & - 1, 9, 0, 0, & - 10,18, 0, 0, & ! 2 - 9,17, 0, 0, & - 17,18, 0, 0, & - 9,10, 0, 0, & - 2,11, 0, 0, & ! 3 - 10,18, 0, 0, & - 11,18, 0, 0, & - 2,10, 0, 0, & - 17,20, 0, 0, & ! 4 - 15,16, 0, 0, & - 15,20, 0, 0, & - 16,17, 0, 0, & - 18,19, 0, 0, & ! 5 - 17,20, 0, 0, & - 19,20, 0, 0, & - 17,18, 0, 0, & - 11,12, 0, 0, & ! 6 - 18,19, 0, 0, & - 12,19, 0, 0, & - 11,18, 0, 0, & - 14,20, 0, 0, & ! 7 - 4,15, 0, 0, & - 4,14, 0, 0, & - 15,20, 0, 0, & - 13,19, 0, 0, & ! 8 - 14,20, 0, 0, & - 13,14, 0, 0, & - 19,20, 0, 0, & - 3,12, 0, 0, & ! 9 - 13,19, 0, 0, & - 3,13, 0, 0, & - 12,19, 0, 0 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(4),FE_Nips(4)/)) - !FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(5),:FE_Nips(5),5) = & ! element 157 - ! reshape((/& - ! *still to be defined* - ! /),(/FE_NipFaceNodes,FE_NipNeighbors(5),FE_Nips(5)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(6),:FE_Nips(6),6) = & ! element 136 - reshape((/& - 7,16,21,17, & ! 1 - 1,10,19, 9, & - 9,19,21,16, & - 1, 7,17,10, & - 10,17,21,19, & - 1, 9,16, 7, & - 2, 8,18,11, & ! 2 - 7,17,21,16, & - 8,16,21,18, & - 2,11,17, 7, & - 11,18,21,17, & - 2, 7,16, 8, & - 8,18,21,16, & ! 3 - 3, 9,19,12, & - 3,12,18, 8, & - 9,16,21,19, & - 12,19,21,18, & - 3, 8,16, 9, & - 13,17,21,20, & ! 4 - 4,15,19,10, & - 15,20,21,19, & - 4,10,17,13, & - 4,13,20,15, & - 10,19,21,17, & - 5,11,18,14, & ! 5 - 13,20,21,17, & - 14,18,21,20, & - 5,13,17,11, & - 5,14,20,13, & - 11,17,21,18, & - 14,20,21,18, & ! 6 - 6,12,19,15, & - 6,14,18,12, & - 15,19,21,20, & - 6,15,20,14, & - 12,18,21,19 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(6),FE_Nips(6)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(7),:FE_Nips(7),7) = & ! element 21 - reshape((/& - 9,33,57,37, & ! 1 - 1,17,44,16, & - 33,16,44,57, & - 1, 9,37,17, & - 17,37,57,44, & - 1,16,33, 9, & - 10,34,58,38, & ! 2 - 9,37,57,33, & - 34,33,57,58, & - 9,10,38,37, & - 37,38,58,57, & - 9,33,34,10, & - 2,11,39,18, & ! 3 - 10,38,58,34, & - 11,34,58,39, & - 10, 2,18,38, & - 38,18,39,58, & - 10,34,11, 2, & - 33,36,60,57, & ! 4 - 16,44,43,15, & - 36,15,43,60, & - 16,33,57,44, & - 44,57,60,43, & - 16,15,36,33, & - 34,35,59,58, & ! 5 - 33,57,60,36, & - 35,36,60,59, & - 33,34,58,57, & - 57,58,59,60, & - 33,36,35,34, & - 11,12,40,39, & ! 6 - 34,58,59,35, & - 12,35,59,40, & - 34,11,39,58, & - 58,39,40,59, & - 34,35,12,11, & - 36,14,42,60, & ! 7 - 15,43,20, 4, & - 14, 4,20,42, & - 15,36,60,43, & - 43,60,42,20, & - 15, 4,14,36, & - 35,13,41,59, & ! 8 - 36,60,42,14, & - 13,14,42,41, & - 36,35,59,60, & - 60,59,41,42, & - 36,14,13,35, & - 12, 3,19,40, & ! 9 - 35,59,41,13, & - 3,13,41,19, & - 35,12,40,59, & - 59,40,19,41, & - 35,13, 3,12, & - 37,57,61,45, & ! 10 - 17,21,52,44, & - 57,44,52,61, & - 17,37,45,21, & - 21,45,61,52, & - 17,44,57,37, & - 38,58,62,46, & ! 11 - 37,45,61,57, & - 58,57,61,62, & - 37,38,46,45, & - 45,46,62,61, & - 37,57,58,38, & - 18,39,47,22, & ! 12 - 38,46,62,58, & - 39,58,62,47, & - 38,18,22,46, & - 46,22,47,62, & - 38,58,39,18, & - 57,60,64,61, & ! 13 - 44,52,51,43, & - 60,43,51,64, & - 44,57,61,52, & - 52,61,64,51, & - 44,43,60,57, & - 58,59,63,62, & ! 14 - 57,61,64,60, & - 59,60,64,63, & - 57,58,62,61, & - 61,62,63,64, & - 57,60,59,58, & - 39,40,48,47, & ! 15 - 58,62,63,59, & - 40,59,63,48, & - 58,39,47,62, & - 62,47,48,63, & - 58,59,40,39, & - 60,42,50,64, & ! 16 - 43,51,24,20, & - 42,20,24,50, & - 43,60,64,51, & - 51,64,50,24, & - 43,20,42,60, & - 59,41,49,63, & ! 17 - 60,64,50,42, & - 41,42,50,49, & - 60,59,63,64, & - 64,63,49,50, & - 60,42,41,59, & - 40,19,23,48, & ! 18 - 59,63,49,41, & - 19,41,49,23, & - 59,40,48,63, & - 63,48,23,49, & - 59,41,19,40, & - 45,61,53,25, & ! 19 - 21, 5,32,52, & - 61,52,32,53, & - 21,45,25, 5, & - 5,25,53,32, & - 21,52,61,45, & - 46,62,54,26, & ! 20 - 45,25,53,61, & - 62,61,53,54, & - 45,46,26,25, & - 25,26,54,53, & - 45,61,62,46, & - 22,47,27, 6, & ! 21 - 46,26,54,62, & - 47,62,54,27, & - 46,22, 6,26, & - 26, 6,27,54, & - 46,62,47,22, & - 61,64,56,53, & ! 22 - 52,32,31,51, & - 64,51,31,56, & - 52,61,53,32, & - 32,53,56,31, & - 52,51,64,61, & - 62,63,55,54, & ! 23 - 61,53,56,64, & - 63,64,56,55, & - 61,62,54,53, & - 53,54,55,56, & - 61,64,63,62, & - 47,48,28,27, & ! 24 - 62,54,55,63, & - 48,63,55,28, & - 62,47,27,54, & - 54,27,28,55, & - 62,63,48,47, & - 64,50,30,56, & ! 25 - 51,31, 8,24, & - 50,24, 8,30, & - 51,64,56,31, & - 31,56,30, 8, & - 51,24,50,64, & - 63,49,29,55, & ! 26 - 64,56,30,50, & - 49,50,30,29, & - 64,63,55,56, & - 56,55,29,30, & - 64,50,49,63, & - 48,23, 7,28, & ! 27 - 63,55,29,49, & - 23,49,29, 7, & - 63,48,28,55, & - 55,28, 7,29, & - 63,49,23,48 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(7),FE_Nips(7)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117 - reshape((/& - 2, 3, 7, 6, & ! 1 - 1, 5, 8, 4, & - 3, 4, 8, 7, & - 1, 2, 6, 5, & - 5, 6, 7, 8, & - 1, 4, 3, 2 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(8),FE_Nips(8)/)) - FE_subNodeOnIPFace(:FE_NipFaceNodes,:FE_NipNeighbors(9),:FE_Nips(9),9) = & ! element 57 (c3d20r == c3d8 --> copy of 7) - reshape((/& - 9,21,27,22, & ! 1 - 1,13,25,12, & - 12,25,27,21, & - 1, 9,22,13, & - 13,22,27,25, & - 1,12,21, 9, & - 2,10,23,14, & ! 2 - 9,22,27,21, & - 10,21,27,23, & - 2,14,22, 9, & - 14,23,27,22, & - 2, 9,21,10, & - 11,24,27,21, & ! 3 - 4,12,25,16, & - 4,16,24,11, & - 12,21,27,25, & - 16,25,27,24, & - 4,11,21,12, & - 3,15,23,10, & ! 4 - 11,21,27,24, & - 3,11,24,15, & - 10,23,27,21, & - 15,24,27,23, & - 3,10,21,11, & - 17,22,27,26, & ! 5 - 5,20,25,13, & - 20,26,27,25, & - 5,13,22,17, & - 5,17,26,20, & - 13,25,27,22, & - 6,14,23,18, & ! 6 - 17,26,27,22, & - 18,23,27,26, & - 6,17,22,14, & - 6,18,26,17, & - 14,22,27,23, & - 19,26,27,24, & ! 7 - 8,16,25,20, & - 8,19,24,16, & - 20,25,27,26, & - 8,20,26,19, & - 16,24,27,25, & - 7,18,23,15, & ! 8 - 19,24,27,26, & - 7,15,24,19, & - 18,26,27,23, & - 7,19,26,18, & - 15,23,27,24 & - /),(/FE_NipFaceNodes,FE_NipNeighbors(9),FE_Nips(9)/)) - - return - - endsubroutine - - -!******************************************************************** -! figure out table styles (Marc only) -! -! initialcondTableStyle, hypoelasticTableStyle -!******************************************************************** - subroutine mesh_marc_get_tableStyles (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 6 - integer(pInt), dimension (1+2*maxNchunks) :: pos - - integer(pInt) unit - character(len=300) line - - initialcondTableStyle = 0_pInt - hypoelasticTableStyle = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_stringValue(line,pos,1)) == 'table' .and. pos(1) .GT. 5) then - initialcondTableStyle = IO_intValue(line,pos,4) - hypoelasticTableStyle = IO_intValue(line,pos,5) - exit - endif - enddo - -620 return - - endsubroutine - - -!******************************************************************** -! get any additional mpie options from input file (Marc only) -! -! mesh_periodicSurface -!******************************************************************** -subroutine mesh_marc_get_mpieOptions(unit) - -use prec, only: pInt -use IO -implicit none - -integer(pInt), intent(in) :: unit - -integer(pInt), parameter :: maxNchunks = 5 -integer(pInt), dimension (1+2*maxNchunks) :: pos -integer(pInt) chunk -character(len=300) line - -mesh_periodicSurface = (/.false., .false., .false./) - -610 FORMAT(A300) - -rewind(unit) -do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - - if (IO_lc(IO_stringValue(line,pos,1)) == '$mpie') then ! found keyword for user defined input - if (IO_lc(IO_stringValue(line,pos,2)) == 'periodic' & ! found keyword 'periodic' - .and. pos(1) > 3) then ! and there is at least one more chunk to read - do chunk = 2,pos(1) ! loop through chunks (skipping the keyword) - select case(IO_stringValue(line,pos,chunk)) ! chunk matches keyvalues x,y or z? - case('x') - mesh_periodicSurface(1) = .true. - case('y') - mesh_periodicSurface(2) = .true. - case('z') - mesh_periodicSurface(3) = .true. - end select - enddo - endif - endif -enddo - -620 return - -endsubroutine - - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** - subroutine mesh_spectral_count_nodesAndElements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension (1+2*maxNchunks) :: pos - integer(pInt) a,b,c,i - - integer(pInt) unit - character(len=1024) line - - mesh_Nnodes = 0_pInt - mesh_Nelems = 0_pInt - - rewind(unit) - do - read(unit,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! skip empty lines - pos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_StringValue(line,pos,1)) == 'resolution') then - do i = 2,6,2 - select case (IO_lc(IO_stringValue(line,pos,i))) - case('a') - a = IO_intValue(line,pos,i+1) - case('b') - b = IO_intValue(line,pos,i+1) - case('c') - c = IO_intValue(line,pos,i+1) - end select - enddo - mesh_Nelems = a * b * c - mesh_Nnodes = (1 + a)*(1 + b)*(1 + c) - exit - endif - enddo - -100 return - - endsubroutine - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** - subroutine mesh_marc_count_nodesAndElements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension (1+2*maxNchunks) :: pos - - integer(pInt) unit - character(len=300) line - - mesh_Nnodes = 0_pInt - mesh_Nelems = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_StringValue(line,pos,1)) == 'sizing') then - mesh_Nelems = IO_IntValue (line,pos,3) - mesh_Nnodes = IO_IntValue (line,pos,4) - exit - endif - enddo - -620 return - - endsubroutine - -!******************************************************************** -! count overall number of nodes and elements in mesh -! -! mesh_Nelems, mesh_Nnodes -!******************************************************************** - subroutine mesh_abaqus_count_nodesAndElements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit - logical inPart - - mesh_Nnodes = 0 - mesh_Nelems = 0 - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if (inPart .or. noPart) then - select case ( IO_lc(IO_stringValue(line,pos,1))) - case('*node') - if( & - IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' & - ) & - mesh_Nnodes = mesh_Nnodes + IO_countDataLines(unit) - case('*element') - if( & - IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' & - ) then - mesh_Nelems = mesh_Nelems + IO_countDataLines(unit) - endif - endselect - endif - enddo - -620 if (mesh_Nnodes < 2) call IO_error(900) - if (mesh_Nelems == 0) call IO_error(901) - - return - - endsubroutine - - -!******************************************************************** -! count overall number of element sets in mesh -! -! mesh_NelemSets, mesh_maxNelemInSet -!******************************************************************** - subroutine mesh_marc_count_elementSets (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - - integer(pInt) unit - character(len=300) line - - mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_StringValue(line,pos,1)) == 'define' .and. & - IO_lc(IO_StringValue(line,pos,2)) == 'element' ) then - mesh_NelemSets = mesh_NelemSets + 1_pInt - mesh_maxNelemInSet = max(mesh_maxNelemInSet, & - IO_countContinousIntValues(unit)) - endif - enddo - -620 return - - endsubroutine - - -!******************************************************************** -! count overall number of element sets in mesh -! -! mesh_NelemSets, mesh_maxNelemInSet -!******************************************************************** - subroutine mesh_abaqus_count_elementSets (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit - logical inPart - - mesh_NelemSets = 0_pInt - mesh_maxNelemInSet = mesh_Nelems ! have to be conservative, since Abaqus allows for recursive definitons - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) & - mesh_NelemSets = mesh_NelemSets + 1_pInt - enddo - -620 continue - if (mesh_NelemSets == 0) call IO_error(902) - - return - endsubroutine - - -!******************************************************************** -! count overall number of solid sections sets in mesh (Abaqus only) -! -! mesh_Nmaterials -!******************************************************************** - subroutine mesh_abaqus_count_materials (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit - logical inPart - - mesh_Nmaterials = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & - IO_lc(IO_StringValue(line,pos,2)) == 'section' ) & - mesh_Nmaterials = mesh_Nmaterials + 1_pInt - enddo - -620 if (mesh_Nmaterials == 0) call IO_error(903) - - return - endsubroutine - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** - subroutine mesh_spectral_count_cpElements () - - use prec, only: pInt - implicit none - - mesh_NcpElems = mesh_Nelems - return - - endsubroutine - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** - subroutine mesh_marc_count_cpElements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: pos - - integer(pInt) unit,i - character(len=300) line - - mesh_NcpElems = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - - if ( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic') then - do i=1,3+hypoelasticTableStyle ! Skip 3 or 4 lines - read (unit,610,END=620) line - enddo - mesh_NcpElems = mesh_NcpElems + IO_countContinousIntValues(unit) - exit - endif - enddo - -620 return - - endsubroutine - - -!******************************************************************** -! count overall number of cpElements in mesh -! -! mesh_NcpElems -!******************************************************************** - subroutine mesh_abaqus_count_cpElements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i - logical materialFound - character (len=64) materialName,elemSetName - - mesh_NcpElems = 0 - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_stringValue(line,pos,1)) ) - case('*material') - materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) exit ! found one - enddo - if (i <= mesh_Nmaterials) then ! found one? - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do i = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(i)) & ! matched? - mesh_NcpElems = mesh_NcpElems + mesh_mapElemSet(1,i) ! add those elem count - enddo - endif - materialFound = .false. - endif - endselect - enddo - -620 if (mesh_NcpElems == 0) call IO_error(906) - - return - endsubroutine - - -!******************************************************************** -! map element sets -! -! allocate globals: mesh_nameElemSet, mesh_mapElemSet -!******************************************************************** - subroutine mesh_marc_map_elementSets (unit) - - use prec, only: pInt - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,elemSet - - allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' - allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt - -610 FORMAT(A300) - - elemSet = 0_pInt - rewind(unit) - do - read (unit,610,END=640) line - pos = IO_stringPos(line,maxNchunks) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'define' ) .and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'element' ) ) then - elemSet = elemSet+1 - mesh_nameElemSet(elemSet) = IO_stringValue(line,pos,4) - mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_maxNelemInSet,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - endif - enddo - -640 return - - endsubroutine - - -!******************************************************************** -! Build element set mapping -! -! allocate globals: mesh_nameElemSet, mesh_mapElemSet -!******************************************************************** - subroutine mesh_abaqus_map_elementSets (unit) - - use prec, only: pInt - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,elemSet,i - logical inPart - - allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' - allocate (mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt - -610 FORMAT(A300) - - elemSet = 0_pInt - inPart = .false. - rewind(unit) - do - read (unit,610,END=640) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) then - elemSet = elemSet + 1_pInt - mesh_nameElemSet(elemSet) = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'elset') - mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,elemSet-1) - endif - enddo - -640 do i = 1,elemSet -! write(6,*)'elemSetName: ',mesh_nameElemSet(i) -! write(6,*)'elems in Elset',mesh_mapElemSet(:,i) - if (mesh_mapElemSet(1,i) == 0) call IO_error(ID=904,ext_msg=mesh_nameElemSet(i)) - enddo - - return - endsubroutine - - -!******************************************************************** -! map solid section (Abaqus only) -! -! allocate globals: mesh_nameMaterial, mesh_mapMaterial -!******************************************************************** - subroutine mesh_abaqus_map_materials (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 20 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,count - logical inPart,materialFound - character(len=64) elemSetName,materialName - - allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' - allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' - -610 FORMAT(A300) - - count = 0 - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if ( (inPart .or. noPart) .and. & - IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. & - IO_lc(IO_StringValue(line,pos,2)) == 'section' ) then - - elemSetName = '' - materialName = '' - - do i = 3,pos(1) - if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') /= '') & - elemSetName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'elset') - if (IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') /= '') & - materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,i)),'material') - enddo - - if (elemSetName /= '' .and. materialName /= '') then - count = count + 1_pInt - mesh_nameMaterial(count) = materialName ! name of material used for this section - mesh_mapMaterial(count) = elemSetName ! mapped to respective element set - endif - endif - enddo - -620 if (count==0) call IO_error(905) - do i=1,count -! write(6,*)'name of materials: ',i,mesh_nameMaterial(i) -! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i) - if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(905) - enddo - - return - endsubroutine - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** - subroutine mesh_spectral_map_nodes () - - use prec, only: pInt - - implicit none - integer(pInt) i - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - - forall (i = 1:mesh_Nnodes) & - mesh_mapFEtoCPnode(:,i) = i - - return - - endsubroutine - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** - subroutine mesh_marc_map_nodes (unit) - - use prec, only: pInt - use math, only: qsort - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt) unit,i - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - -610 FORMAT(A300) - - node_count = 0_pInt - - rewind(unit) - do - read (unit,610,END=650) line - pos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=650) line ! skip crap line - do i = 1,mesh_Nnodes - read (unit,610,END=650) line - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) - mesh_mapFEtoCPnode(2,i) = i - enddo - exit - endif - enddo - -650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) - - return - - endsubroutine - - - -!******************************************************************** -! map nodes from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPnode -!******************************************************************** - subroutine mesh_abaqus_map_nodes (unit) - - use prec, only: pInt - use math, only: qsort - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,count,cpNode - logical inPart - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - -610 FORMAT(A300) - - cpNode = 0_pInt - inPart = .false. - rewind(unit) - do - read (unit,610,END=650) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & - ) then - count = IO_countDataLines(unit) - do i = 1,count - backspace(unit) - enddo - do i = 1,count - read (unit,610,END=650) line - pos = IO_stringPos(line,maxNchunks) - cpNode = cpNode + 1_pInt - mesh_mapFEtoCPnode(1,cpNode) = IO_intValue(line,pos,1) - mesh_mapFEtoCPnode(2,cpNode) = cpNode - enddo - endif - enddo - -650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) - - if (size(mesh_mapFEtoCPnode) == 0) call IO_error(908) - return - - endsubroutine - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** - subroutine mesh_spectral_map_elements () - - use prec, only: pInt - - implicit none - integer(pInt) i - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - - forall (i = 1:mesh_NcpElems) & - mesh_mapFEtoCPelem(:,i) = i - - return - - endsubroutine - - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** - subroutine mesh_marc_map_elements (unit) - - use prec, only: pInt - use math, only: qsort - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt), dimension (1+mesh_NcpElems) :: contInts - integer(pInt) unit,i,cpElem - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - -610 FORMAT(A300) - - cpElem = 0_pInt - rewind(unit) - do - read (unit,610,END=660) line - pos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then - do i=1,3+hypoelasticTableStyle ! skip three (or four if new table style!) lines - read (unit,610,END=660) line - enddo - contInts = IO_continousIntValues(unit,mesh_NcpElems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) - do i = 1,contInts(1) - cpElem = cpElem+1 - mesh_mapFEtoCPelem(1,cpElem) = contInts(1+i) - mesh_mapFEtoCPelem(2,cpElem) = cpElem - enddo - endif - enddo - -660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems - - return - endsubroutine - - -!******************************************************************** -! map elements from FE id to internal (consecutive) representation -! -! allocate globals: mesh_mapFEtoCPelem -!******************************************************************** - subroutine mesh_abaqus_map_elements (unit) - - use prec, only: pInt - use math, only: qsort - use IO - - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,j,cpElem - logical materialFound - character (len=64) materialName,elemSetName - - allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt - -610 FORMAT(A300) - - cpElem = 0_pInt - rewind(unit) - do - read (unit,610,END=660) line - pos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_stringValue(line,pos,1)) ) - case('*material') - materialName = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'name') ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if (IO_lc(IO_stringValue(line,pos,2)) == 'material' .and. materialFound) then - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) exit ! found one - enddo - if (i <= mesh_Nmaterials) then ! found one? - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - do i = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(i)) then ! matched? - do j = 1,mesh_mapElemSet(1,i) - cpElem = cpElem + 1_pInt - mesh_mapFEtoCPelem(1,cpElem) = mesh_mapElemSet(1+j,i) ! store FE id - mesh_mapFEtoCPelem(2,cpElem) = cpElem ! store our id - enddo - endif - enddo - endif - materialFound = .false. - endif - endselect - enddo - -660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems - - if (size(mesh_mapFEtoCPelem) < 2) call IO_error(907) - - return - endsubroutine - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** -subroutine mesh_spectral_count_cpSizes () - - use prec, only: pInt - implicit none - - integer(pInt) t - - t = FE_mapElemtype('C3D8R') ! fake 3D hexahedral 8 node 1 IP element - - mesh_maxNnodes = FE_Nnodes(t) - mesh_maxNips = FE_Nips(t) - mesh_maxNipNeighbors = FE_NipNeighbors(t) - mesh_maxNsubNodes = FE_NsubNodes(t) - - endsubroutine - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** -subroutine mesh_marc_count_cpSizes (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,t,e - - mesh_maxNnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNipNeighbors = 0_pInt - mesh_maxNsubNodes = 0_pInt - -610 FORMAT(A300) - rewind(unit) - do - read (unit,610,END=630) line - pos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=630) line ! Garbage line - do i=1,mesh_Nelems ! read all elements - read (unit,610,END=630) line - pos = IO_stringPos(line,maxNchunks) ! limit to id and type - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then - t = FE_mapElemtype(IO_stringValue(line,pos,2)) - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - endif - enddo - exit - endif - enddo - -630 return - endsubroutine - - -!******************************************************************** -! get maximum count of nodes, IPs, IP neighbors, and subNodes -! among cpElements -! -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsubNodes -!******************************************************************** - subroutine mesh_abaqus_count_cpSizes (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,count,t - logical inPart - - mesh_maxNnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNipNeighbors = 0_pInt - mesh_maxNsubNodes = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type - if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_count_cpSizes') - count = IO_countDataLines(unit) - do i = 1,count - backspace(unit) - enddo - do i = 1,count - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - if (mesh_FEasCP('elem',IO_intValue(line,pos,1)) /= 0) then ! disregard non CP elems - mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) - mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) - mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) - mesh_maxNsubNodes = max(mesh_maxNsubNodes,FE_NsubNodes(t)) - endif - enddo - endif - enddo - -620 return - - endsubroutine - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** - subroutine mesh_spectral_build_nodes (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension (1+2*maxNchunks) :: pos - integer(pInt) a,b,c,n,i - real(pReal) x,y,z - logical gotResolution,gotDimension - - integer(pInt) unit - character(len=64) tag - character(len=1024) line - - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt - - a = 1_pInt - b = 1_pInt - c = 1_pInt - x = 1.0_pReal - y = 1.0_pReal - z = 1.0_pReal - - gotResolution = .false. - gotDimension = .false. - - rewind(unit) - do - read(unit,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! skip empty lines - pos = IO_stringPos(line,maxNchunks) - - select case ( IO_lc(IO_StringValue(line,pos,1)) ) - case ('resolution') - gotResolution = .true. - do i = 2,6,2 - tag = IO_lc(IO_stringValue(line,pos,i)) - select case (tag) - case('a') - a = 1 + IO_intValue(line,pos,i+1) - case('b') - b = 1 + IO_intValue(line,pos,i+1) - case('c') - c = 1 + IO_intValue(line,pos,i+1) - end select - enddo - case ('dimension') - gotDimension = .true. - do i = 2,6,2 - tag = IO_lc(IO_stringValue(line,pos,i)) - select case (tag) - case('x') - x = IO_floatValue(line,pos,i+1) - case('y') - y = IO_floatValue(line,pos,i+1) - case('z') - z = IO_floatValue(line,pos,i+1) - end select - enddo - end select - if (gotDimension .and. gotResolution) exit - enddo - -! --- sanity checks --- - - if (.not. gotDimension .or. .not. gotResolution) call IO_error(42) - if (a < 2 .or. b < 2 .or. c < 2) call IO_error(43) - if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(44) - - forall (n = 0:mesh_Nnodes-1) - mesh_node(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal)) - mesh_node(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal)) - mesh_node(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal)) - end forall - -100 return - - endsubroutine - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** - subroutine mesh_marc_build_nodes (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) - integer(pInt), parameter :: maxNchunks = 1 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,j,m - - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=670) line - pos = IO_stringPos(line,maxNchunks) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=670) line ! skip crap line - do i=1,mesh_Nnodes - read (unit,610,END=670) line - m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1)) - forall (j = 1:3) mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) - enddo - exit - endif - enddo - -670 return - - endsubroutine - - -!******************************************************************** -! store x,y,z coordinates of all nodes in mesh -! -! allocate globals: -! _node -!******************************************************************** - subroutine mesh_abaqus_build_nodes (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt) unit,i,j,m,count - logical inPart - - allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=670) line - pos = IO_stringPos(line,maxNchunks) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. & - ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'file' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & - ) then - count = IO_countDataLines(unit) ! how many nodes are defined here? - do i = 1,count - backspace(unit) ! rewind to first entry - enddo - do i = 1,count - read (unit,610,END=670) line - pos = IO_stringPos(line,maxNchunks) - m = mesh_FEasCP('node',IO_intValue(line,pos,1)) - forall (j=1:3) mesh_node(j,m) = IO_floatValue(line,pos,j+1) - enddo - endif - enddo - -670 if (size(mesh_node,2) /= mesh_Nnodes) call IO_error(909) - return - - endsubroutine - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** - subroutine mesh_spectral_build_elements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension (1+2*maxNchunks) :: pos - integer(pInt) a,b,c,e,i,homog - logical gotResolution,gotDimension,gotHomogenization - - integer(pInt) unit - character(len=1024) line - - a = 1_pInt - b = 1_pInt - c = 1_pInt - - gotResolution = .false. - gotDimension = .false. - gotHomogenization = .false. - - rewind(unit) - do - read(unit,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! skip empty lines - pos = IO_stringPos(line,maxNchunks) - - select case ( IO_lc(IO_StringValue(line,pos,1)) ) - case ('dimension') - gotDimension = .true. - - case ('homogenization') - gotHomogenization = .true. - homog = IO_intValue(line,pos,2) - - case ('resolution') - gotResolution = .true. - do i = 2,6,2 - select case (IO_lc(IO_stringValue(line,pos,i))) - case('a') - a = IO_intValue(line,pos,i+1) - case('b') - b = IO_intValue(line,pos,i+1) - case('c') - c = IO_intValue(line,pos,i+1) - end select - enddo - end select - if (gotDimension .and. gotHomogenization .and. gotResolution) exit - enddo - -100 allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - - e = 0_pInt - do while (e < mesh_NcpElems) - read(unit,'(a1024)',END=110) line - if (IO_isBlank(line)) cycle ! skip empty lines - pos(1:1+2*1) = IO_stringPos(line,1) - - e = e+1 ! valid element entry - mesh_element ( 1,e) = e ! FE id - mesh_element ( 2,e) = FE_mapElemtype('C3D8R') ! elem type - mesh_element ( 3,e) = homog ! homogenization - mesh_element ( 4,e) = IO_IntValue(line,pos,1) ! microstructure - mesh_element ( 5,e) = e + (e-1)/a + (e-1)/a/b*(a+1) ! base node - mesh_element ( 6,e) = mesh_element ( 5,e) + 1 - mesh_element ( 7,e) = mesh_element ( 5,e) + (a+1) + 1 - mesh_element ( 8,e) = mesh_element ( 5,e) + (a+1) - mesh_element ( 9,e) = mesh_element ( 5,e) + (a+1)*(b+1) ! second floor base node - mesh_element (10,e) = mesh_element ( 9,e) + 1 - mesh_element (11,e) = mesh_element ( 9,e) + (a+1) + 1 - mesh_element (12,e) = mesh_element ( 9,e) + (a+1) - mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) !needed for statistics - mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e)) - enddo - -110 return - - endsubroutine - - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** - subroutine mesh_marc_build_elements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 66 - integer(pInt), dimension (1+2*maxNchunks) :: pos - character(len=300) line - - integer(pInt), dimension(1+mesh_NcpElems) :: contInts - integer(pInt) unit,i,j,sv,val,e - - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos(1:1+2*1) = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=620) line ! Garbage line - do i = 1,mesh_Nelems - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type) - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then ! disregard non CP elems - mesh_element(1,e) = IO_IntValue (line,pos,1) ! FE id - mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type - forall (j = 1:FE_Nnodes(mesh_element(2,e))) & - mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes - call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line - endif - enddo - exit - endif - enddo - -620 rewind(unit) ! just in case "initial state" apears before "connectivity" - read (unit,610,END=620) line - do - pos(1:1+2*2) = IO_stringPos(line,2) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then - if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style - read (unit,610,END=630) line ! read line with index of state var - pos(1:1+2*1) = IO_stringPos(line,1) - sv = IO_IntValue(line,pos,1) ! figure state variable index - if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest - read (unit,610,END=620) line ! read line with value of state var - pos(1:1+2*1) = IO_stringPos(line,1) - do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? - val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value - mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index - if (initialcondTableStyle == 2) then - read (unit,610,END=630) line ! read extra line - read (unit,610,END=630) line ! read extra line - endif - contInts = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,mesh_NelemSets) ! get affected elements - do i = 1,contInts(1) - e = mesh_FEasCP('elem',contInts(1+i)) - mesh_element(1+sv,e) = val - enddo - if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style - read (unit,610,END=630) line - pos(1:1+2*1) = IO_stringPos(line,1) - enddo - endif - else - read (unit,610,END=630) line - endif - enddo - -630 return - - endsubroutine - - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** - subroutine mesh_abaqus_build_elements (unit) - - use prec, only: pInt - use IO - implicit none - - integer(pInt), parameter :: maxNchunks = 65 - integer(pInt), dimension (1+2*maxNchunks) :: pos - - integer(pInt) unit,i,j,count,e,t,homog,micro - logical inPart,materialFound - character (len=64) materialName,elemSetName - character(len=300) line - - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - - inPart = .false. - rewind(unit) - do - read (unit,610,END=620) line - pos(1:1+2*2) = IO_stringPos(line,2) - if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. - if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & - IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. - - if( (inPart .or. noPart) .and. & - IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. & - ( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. & - IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) & - ) then - t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type - if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_build_elements') - count = IO_countDataLines(unit) - do i = 1,count - backspace(unit) - enddo - do i = 1,count - read (unit,610,END=620) line - pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max - e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (e /= 0) then ! disregard non CP elems - mesh_element(1,e) = IO_intValue(line,pos,1) ! FE id - mesh_element(2,e) = t ! elem type - forall (j=1:FE_Nnodes(t)) & - mesh_element(4+j,e) = IO_intValue(line,pos,1+j) ! copy FE ids of nodes to position 5: - call IO_skipChunks(unit,FE_NoriginalNodes(t)-(pos(1)-1)) ! read on (even multiple lines) if FE_NoriginalNodes exceeds required node count - endif - enddo - endif - enddo - - -620 rewind(unit) ! just in case "*material" definitions apear before "*element" - - materialFound = .false. - do - read (unit,610,END=630) line - pos = IO_stringPos(line,maxNchunks) - select case ( IO_lc(IO_StringValue(line,pos,1))) - case('*material') - materialName = IO_extractValue(IO_lc(IO_StringValue(line,pos,2)),'name') ! extract name=value - materialFound = materialName /= '' ! valid name? - case('*user') - if ( IO_lc(IO_StringValue(line,pos,2)) == 'material' .and. & - materialFound ) then - do i = 1,mesh_Nmaterials ! look thru material names - if (materialName == mesh_nameMaterial(i)) exit ! found one - enddo - if (i <= mesh_Nmaterials) then ! found one? - elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet - read (unit,610,END=630) line ! read homogenization and microstructure - pos(1:1+2*2) = IO_stringPos(line,2) - homog = NINT(IO_floatValue(line,pos,1)) - micro = NINT(IO_floatValue(line,pos,2)) - do i = 1,mesh_NelemSets ! look thru all elemSet definitions - if (elemSetName == mesh_nameElemSet(i)) then ! matched? - do j = 1,mesh_mapElemSet(1,i) - e = mesh_FEasCP('elem',mesh_mapElemSet(1+j,i)) - mesh_element(3,e) = homog ! store homogenization - mesh_element(4,e) = micro ! store microstructure - mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),homog) - mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),micro) - enddo - endif - enddo - endif - materialFound = .false. - endif - endselect - enddo - -630 return - - endsubroutine - - -!******************************************************************** -! get maximum count of shared elements among cpElements and -! build list of elements shared by each node in mesh -! -! _maxNsharedElems -! _sharedElem -!******************************************************************** -subroutine mesh_build_sharedElems() - -use prec, only: pInt -implicit none - -integer(pint) e, & ! element index - t, & ! element type - node, & ! CP node index - j, & ! node index per element - dim, & ! dimension index - nodeTwin ! node twin in the specified dimension -integer(pInt), dimension (mesh_Nnodes) :: node_count -integer(pInt), dimension (:), allocatable :: node_seen - -allocate(node_seen(maxval(FE_Nnodes))) - -node_count = 0_pInt - -do e = 1,mesh_NcpElems - t = mesh_element(2,e) - - node_seen = 0_pInt ! reset node duplicates - do j = 1,FE_Nnodes(t) ! check each node of element - node = mesh_FEasCP('node',mesh_element(4+j,e)) - if (all(node_seen /= node)) then - node_count(node) = node_count(node) + 1_pInt ! if FE node not yet encountered -> count it - do dim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(dim,node) - if (nodeTwin > 0) & ! if i am a twin of some node... - node_count(nodeTwin) = node_count(nodeTwin) + 1_pInt ! -> count me again for the twin node - enddo - endif - node_seen(j) = node ! remember this node to be counted already - enddo -enddo - -mesh_maxNsharedElems = maxval(node_count) ! most shared node - -allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes)) -mesh_sharedElem = 0_pInt - -do e = 1,mesh_NcpElems - t = mesh_element(2,e) - node_seen = 0_pInt - do j = 1,FE_Nnodes(t) - node = mesh_FEasCP('node',mesh_element(4+j,e)) - if (all(node_seen /= node)) then - mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1_pInt ! count for each node the connected elements - mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id - do dim = 1,3 ! check in each dimension... - nodeTwin = mesh_nodeTwins(dim,node) - if (nodeTwin > 0) then ! if i am a twin of some node... - mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1_pInt ! ...count me again for the twin - mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id - endif - enddo - endif - node_seen(j) = node - enddo -enddo - -deallocate(node_seen) - -endsubroutine - - - -!*********************************************************** -! build up of IP neighborhood -! -! allocate globals -! _ipNeighborhood -!*********************************************************** -subroutine mesh_build_ipNeighborhood() - -use prec, only: pInt -implicit none - -integer(pInt) myElem, & ! my CP element index - myIP, & - myType, & ! my element type - myFace, & - neighbor, & ! neighor index - neighboringIPkey, & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element) - candidateIP, & - neighboringType, & ! element type of neighbor - NlinkedNodes, & ! number of linked nodes - twin_of_linkedNode, & ! node twin of a specific linkedNode - NmatchingNodes, & ! number of matching nodes - dir, & ! direction of periodicity - matchingElem, & ! CP elem number of matching element - matchingFace, & ! face ID of matching element - k, a, anchor -integer(pInt), dimension(FE_maxmaxNnodesAtIP) :: & - linkedNodes, & - matchingNodes -logical checkTwins - -allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) -mesh_ipNeighborhood = 0_pInt - -linkedNodes = 0_pInt - -do myElem = 1,mesh_NcpElems ! loop over cpElems - myType = mesh_element(2,myElem) ! get elemType - do myIP = 1,FE_Nips(myType) ! loop over IPs of elem - do neighbor = 1,FE_NipNeighbors(myType) ! loop over neighbors of IP - neighboringIPkey = FE_ipNeighbor(neighbor,myIP,myType) - - !*** if the key is positive, the neighbor is inside the element - !*** that means, we have already found our neighboring IP - - if (neighboringIPkey > 0) then - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey - - - !*** if the key is negative, the neighbor resides in a neighboring element - !*** that means, we have to look through the face indicated by the key and see which element is behind that face - - elseif (neighboringIPkey < 0) then ! neighboring element's IP - myFace = -neighboringIPkey - call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match - if (matchingElem > 0_pInt) then ! found match? - neighboringType = mesh_element(2,matchingElem) - - !*** trivial solution if neighbor has only one IP - - if (FE_Nips(neighboringType) == 1_pInt) then - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1_pInt - cycle - endif - - !*** find those nodes which build the link to the neighbor - - NlinkedNodes = 0_pInt - linkedNodes = 0_pInt - do a = 1,FE_maxNnodesAtIP(myType) ! figure my anchor nodes on connecting face - anchor = FE_nodesAtIP(a,myIP,myType) - if (anchor /= 0_pInt) then ! valid anchor node - if (any(FE_nodeOnFace(:,myFace,myType) == anchor)) then ! ip anchor sits on face? - NlinkedNodes = NlinkedNodes + 1_pInt - linkedNodes(NlinkedNodes) = mesh_element(4+anchor,myElem) ! FE id of anchor node - else ! something went wrong with the linkage, since not all anchors sit on my face - NlinkedNodes = 0_pInt - linkedNodes = 0_pInt - exit - endif - endif - enddo - - !*** loop through the ips of my neighbor - !*** and try to find an ip with matching nodes - !*** also try to match with node twins - -checkCandidateIP: do candidateIP = 1,FE_Nips(neighboringType) - NmatchingNodes = 0_pInt - matchingNodes = 0_pInt - do a = 1,FE_maxNnodesAtIP(neighboringType) ! check each anchor node of that ip - anchor = FE_nodesAtIP(a,candidateIP,neighboringType) - if (anchor /= 0_pInt) then ! valid anchor node - if (any(FE_nodeOnFace(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face? - NmatchingNodes = NmatchingNodes + 1_pInt - matchingNodes(NmatchingNodes) = mesh_element(4+anchor,matchingElem) ! FE id of neighbor's anchor node - else ! no matching, because not all nodes sit on the matching face - NmatchingNodes = 0_pInt - matchingNodes = 0_pInt - exit - endif - endif - enddo - - if (NmatchingNodes /= NlinkedNodes) & ! this ip has wrong count of anchors on face - cycle checkCandidateIP - - !*** check "normal" nodes whether they match or not - - checkTwins = .false. - do a = 1,NlinkedNodes - if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode - checkTwins = .true. - exit ! no need to search further - endif - enddo - - !*** if no match found, then also check node twins - - if(checkTwins) then -periodicityDirection: do dir = 1,3 - do a = 1,NlinkedNodes - twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a)) - if (twin_of_linkedNode == 0_pInt & ! twin of linkedNode does not exist... - .or. all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode - if (dir < 3) then ! no match in this direction... - cycle periodicityDirection ! ... so try in different direction - else ! no matching in any direction... - cycle checkCandidateIP ! ... so check next candidateIP - endif - endif - enddo - exit periodicityDirection - enddo periodicityDirection - endif - - !*** we found a match !!! - - mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem - mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP - exit checkCandidateIP - - enddo checkCandidateIP - endif ! end of valid external matching - endif ! end of internal/external matching - enddo - enddo -enddo - -endsubroutine - - - -!*********************************************************** -! assignment of coordinates for subnodes in each cp element -! -! allocate globals -! _subNodeCoord -!*********************************************************** - subroutine mesh_build_subNodeCoords() - - use prec, only: pInt,pReal - implicit none - - integer(pInt) e,t,n,p - - allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) ; mesh_subNodeCoord = 0.0_pReal - - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do n = 1,FE_Nnodes(t) - mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type - enddo - do n = 1,FE_NsubNodes(t) ! now for the true subnodes - do p = 1,FE_Nips(t) ! loop through possible parent nodes - if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = & - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + & - mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents - enddo - mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0) - enddo - enddo - - return - - endsubroutine - - -!*********************************************************** -! calculation of IP volume -! -! allocate globals -! _ipVolume -!*********************************************************** - subroutine mesh_build_ipVolumes() - - use prec, only: pInt - use math, only: math_volTetrahedron - implicit none - - integer(pInt) e,f,t,i,j,k,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles - logical(pInt), dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav - real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav - real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face - real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra - real(pReal), dimension(3) :: centerOfGravity - - allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal - allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal - - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - gravityNode = .false. ! reset flagList - gravityNodePos = 0.0_pReal ! reset coordinates - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1,FE_NipFaceNodes ! loop over nodes on interface - gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true. - gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - enddo - enddo - - do j = 1,mesh_maxNnodes+mesh_maxNsubNodes-1 ! walk through entire flagList except last - if (gravityNode(j)) then ! valid node index - do k = j+1,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list - if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < 1.0e-36_pReal)) then ! found duplicate - gravityNode(j) = .false. ! delete first instance - gravityNodePos(:,j) = 0.0_pReal - exit ! continue with next suspect - endif - enddo - endif - enddo - centerOfGravity = sum(gravityNodePos,2)/count(gravityNode) - - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG - forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface - volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG - nPos(:,1+mod(n+j-1,FE_NipFaceNodes)), & - nPos(:,1+mod(n+j-0,FE_NipFaceNodes)), & - centerOfGravity) - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface - enddo - mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them - mesh_ipCenterOfGravity(:,i,e) = centerOfGravity - enddo - enddo - return - - endsubroutine - - -!*********************************************************** -! calculation of IP interface areas -! -! allocate globals -! _ipArea, _ipAreaNormal -!*********************************************************** - subroutine mesh_build_ipAreas() - - use prec, only: pInt,pReal - use math - implicit none - - integer(pInt) e,f,t,i,j,n - integer(pInt), parameter :: Ntriangles = FE_NipFaceNodes-2 ! each interface is made up of this many triangles - real(pReal), dimension (3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face - real(pReal), dimension(3,Ntriangles,FE_NipFaceNodes) :: normal - real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: area - - allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal - allocate(mesh_ipAreaNormal(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal - - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e) - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) ! start at each interface node and build valid triangles to cover interface - normal(:,j,n) = math_vectorproduct(nPos(:,1+mod(n+j-1,FE_NipFaceNodes)) - nPos(:,n), & ! calc their normal vectors - nPos(:,1+mod(n+j-0,FE_NipFaceNodes)) - nPos(:,n)) - area(j,n) = sqrt(sum(normal(:,j,n)*normal(:,j,n))) ! and area - end forall - forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles, area(j,n) > 0.0_pReal) & - normal(:,j,n) = normal(:,j,n) / area(j,n) ! make unit normal - - mesh_ipArea(f,i,e) = sum(area) / (FE_NipFaceNodes*2.0_pReal) ! area of parallelograms instead of triangles - mesh_ipAreaNormal(:,f,i,e) = sum(sum(normal,3),2) / count(area > 0.0_pReal) ! average of all valid normals - enddo - enddo - enddo - return - - endsubroutine - - -!*********************************************************** -! assignment of twin nodes for each cp node -! -! allocate globals -! _nodeTwins -!*********************************************************** -subroutine mesh_build_nodeTwins() - -use prec, only: pInt, pReal -implicit none - -integer(pInt) dir, & ! direction of periodicity - node, & - minimumNode, & - maximumNode, & - n1, & - n2 -integer(pInt), dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes -real(pReal) minCoord, maxCoord, & ! extreme positions in one dimension - tolerance ! tolerance below which positions are assumed identical -real(pReal), dimension(3) :: distance ! distance between two nodes in all three coordinates -logical, dimension(mesh_Nnodes) :: node_seen - -allocate(mesh_nodeTwins(3,mesh_Nnodes)) -mesh_nodeTwins = 0_pInt - -tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal - -do dir = 1,3 ! check periodicity in directions of x,y,z - if (mesh_periodicSurface(dir)) then ! only if periodicity is requested - - - !*** find out which nodes sit on the surface - !*** and have a minimum or maximum position in this dimension - - minimumNodes = 0_pInt - maximumNodes = 0_pInt - minCoord = minval(mesh_node(dir,:)) - maxCoord = maxval(mesh_node(dir,:)) - do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes - if (abs(mesh_node(dir,node) - minCoord) <= tolerance) then - minimumNodes(1) = minimumNodes(1) + 1_pInt - minimumNodes(minimumNodes(1)+1) = node - elseif (abs(mesh_node(dir,node) - maxCoord) <= tolerance) then - maximumNodes(1) = maximumNodes(1) + 1_pInt - maximumNodes(maximumNodes(1)+1) = node - endif - enddo - - - !*** find the corresponding node on the other side with the same position in this dimension - - node_seen = .false. - do n1 = 1,minimumNodes(1) - minimumNode = minimumNodes(n1+1) - if (node_seen(minimumNode)) & - cycle - do n2 = 1,maximumNodes(1) - maximumNode = maximumNodes(n2+1) - distance = abs(mesh_node(:,minimumNode) - mesh_node(:,maximumNode)) - if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance) - mesh_nodeTwins(dir,minimumNode) = maximumNode - mesh_nodeTwins(dir,maximumNode) = minimumNode - node_seen(maximumNode) = .true. ! remember this node, we don't have to look for his partner again - exit - endif - enddo - enddo - - endif -enddo - -endsubroutine - - - -!*********************************************************** -! write statistics regarding input file parsing -! to the output file -! -!*********************************************************** -subroutine mesh_tell_statistics() - -use prec, only: pInt -use math, only: math_range -use IO, only: IO_error -use debug, only: verboseDebugger - -implicit none - -integer(pInt), dimension (:,:), allocatable :: mesh_HomogMicro -character(len=64) fmt - -integer(pInt) i,e,n,f,t - -if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(110) ! no homogenization specified -if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(120) ! no microstructure specified - -allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt -do e = 1,mesh_NcpElems - if (mesh_element(3,e) < 1_pInt) call IO_error(110,e) ! no homogenization specified - if (mesh_element(4,e) < 1_pInt) call IO_error(120,e) ! no homogenization specified - mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = & - mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure -enddo - -if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write (6,*) - write (6,*) "Input Parser: SUBNODE COORDINATES" - write (6,*) - write(6,'(a5,x,a5,x,a15,x,a15,x,a20,3(x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - do n = 1,FE_NipFaceNodes ! loop over nodes on interface - write(6,'(i5,x,i5,x,i15,x,i15,x,i20,3(x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& - mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& - mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& - mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) - enddo - enddo - enddo - enddo - write(6,*) - write(6,*) 'Input Parser: IP COORDINATES' - write(6,'(a5,x,a5,3(x,a12))') 'elem','IP','x','y','z' - do e = 1,mesh_NcpElems - do i = 1,FE_Nips(mesh_element(2,e)) - write (6,'(i5,x,i5,3(x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e) - enddo - enddo - write (6,*) - write (6,*) "Input Parser: ELEMENT VOLUME" - write (6,*) - write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume) - write (6,*) - write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --" - do e = 1,mesh_NcpElems - do i = 1,FE_Nips(mesh_element(2,e)) - write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) - do f = 1,FE_NipNeighbors(mesh_element(2,e)) - write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) - enddo - enddo - enddo - write (6,*) - write (6,*) "Input Parser: NODE TWINS" - write (6,*) - write(6,'(a6,3(3(x),a6))') ' node','twin_x','twin_y','twin_z' - do n = 1,mesh_Nnodes ! loop over cpNodes - write(6,'(i6,3(3(x),i6))') n, mesh_nodeTwins(1:3,n) - enddo - write(6,*) - write(6,*) "Input Parser: IP NEIGHBORHOOD" - write(6,*) - write(6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor" - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP - write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) - enddo - enddo - enddo - !$OMP END CRITICAL (write2out) -endif - - -!$OMP CRITICAL (write2out) - write (6,*) - write (6,*) "Input Parser: STATISTICS" - write (6,*) - write (6,*) mesh_Nelems, " : total number of elements in mesh" - write (6,*) mesh_NcpElems, " : total number of CP elements in mesh" - write (6,*) mesh_Nnodes, " : total number of nodes in mesh" - write (6,*) mesh_maxNnodes, " : max number of nodes in any CP element" - write (6,*) mesh_maxNips, " : max number of IPs in any CP element" - write (6,*) mesh_maxNipNeighbors, " : max number of IP neighbors in any CP element" - write (6,*) mesh_maxNsubNodes, " : max number of (additional) subnodes in any CP element" - write (6,*) mesh_maxNsharedElems, " : max number of CP elements sharing a node" - write (6,*) - write (6,*) "Input Parser: HOMOGENIZATION/MICROSTRUCTURE" - write (6,*) - write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index" - write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index" - write (6,*) - write (fmt,"(a,i5,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))" - write (6,fmt) "+-",math_range(mesh_maxValStateVar(2)) - write (fmt,"(a,i5,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))" - do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations - write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstrcutures - enddo - write(6,*) - write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS" - write(6,*) - write(6,*) "periodic surface : ", mesh_periodicSurface - write(6,*) - call flush(6) -!$OMP END CRITICAL (write2out) - -deallocate(mesh_HomogMicro) - -endsubroutine - - -END MODULE mesh -