diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index e81da9537..6ea3d0845 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -16,11 +16,14 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## -!* $Id$ -!############################################################## +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief reading in of data when doing a restart +!-------------------------------------------------------------------------------------------------- module FEsolving -!############################################################## use prec, only: & pInt, & pReal @@ -37,15 +40,15 @@ module FEsolving theDelta = 0.0_pReal logical, public :: & - outdatedFFN1 = .false., & - symmetricSolver = .false., & - restartWrite = .false., & - restartRead = .false., & - terminallyIll = .false., & - parallelExecution = .true., & - lastMode = .true., & - lastIncConverged = .false., & - outdatedByNewInc = .false., & + outdatedFFN1 = .false., & !< toDo + symmetricSolver = .false., & !< use a symmetric solver (FEM) + restartWrite = .false., & !< write current state to enable restart + restartRead = .false., & !< restart information to continue calculation from saved state + terminallyIll = .false., & !< at least one material point is terminally ill + parallelExecution = .true., & !< OpenMP multicore calculation + lastMode = .true., & !< toDo + lastIncConverged = .false., & !< toDo + outdatedByNewInc = .false., & !< toDo cutBack = .false. integer(pInt), dimension(:,:), allocatable, public :: & @@ -64,10 +67,9 @@ module FEsolving contains -!*********************************************************** -! determine whether a symmetric solver is used -! and whether restart is requested -!*********************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief determine whether a symmetric solver is used and whether restart is requested +!-------------------------------------------------------------------------------------------------- subroutine FE_init use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) @@ -114,21 +116,21 @@ subroutine FE_init call IO_warning(warning_ID=34_pInt) restartInc = 1_pInt endif - restartRead = restartInc > 1_pInt ! only read in if "true" restart requested + restartRead = restartInc > 1_pInt ! only read in if "true" restart requested #else call IO_open_inputFile(fileunit,modelName) rewind(fileunit) do read (fileunit,'(a1024)',END=100) line positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('solver') - read (fileunit,'(a1024)',END=100) line ! next line + read (fileunit,'(a1024)',END=100) line ! next line positions = IO_stringPos(line,maxNchunks) symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt) case ('restart') - read (fileunit,'(a1024)',END=100) line ! next line + read (fileunit,'(a1024)',END=100) line ! next line positions = IO_stringPos(line,maxNchunks) restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt diff --git a/code/debug.f90 b/code/debug.f90 index f8f55bb3a..8dfcec6be 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -1,7 +1,7 @@ -! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH +! Copyright 2011,2012 Max-Planck-Institut für Eisenforschung GmbH ! ! This file is part of DAMASK, -! the Düsseldorf Advanced MAterial Simulation Kit. +! the Düsseldorf Advanced Material Simulation Kit. ! ! DAMASK is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by @@ -16,11 +16,16 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## +!-------------------------------------------------------------------------------------------------- !* $Id$ -!############################################################## +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Reading in and interpretating the debugging settings for the various modules +!-------------------------------------------------------------------------------------------------- module debug -!############################################################## use prec, only: & pInt, & pReal, & @@ -44,10 +49,10 @@ module debug debug_debug = 1_pInt, & debug_math = 2_pInt, & debug_FEsolving = 3_pInt, & - debug_mesh = 4_pInt, & ! stores debug level for mesh part of DAMASK - debug_material = 5_pInt, & ! stores debug level for material part of DAMASK - debug_lattice = 6_pInt, & ! stores debug level for lattice part of DAMASK - debug_constitutive = 7_pInt, & ! stores debug level for constitutive part of DAMASK + debug_mesh = 4_pInt, & !< stores debug level for mesh part of DAMASK bitwise coded + debug_material = 5_pInt, & !< stores debug level for material part of DAMASK bitwise coded + debug_lattice = 6_pInt, & !< stores debug level for lattice part of DAMASK bitwise coded + debug_constitutive = 7_pInt, & !< stores debug level for constitutive part of DAMASK bitwise coded debug_crystallite = 8_pInt, & debug_homogenization = 9_pInt, & debug_CPFEM = 10_pInt, & diff --git a/code/math.f90 b/code/math.f90 index ae4a9be63..dbdccde6a 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -16,6 +16,9 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! +#ifdef Spectral +#include "kdtree2.f90" +#endif !-------------------------------------------------------------------------------------------------- !* $Id$ !-------------------------------------------------------------------------------------------------- @@ -24,14 +27,11 @@ !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @brief Mathematical library, including random number generation and tensor represenations !-------------------------------------------------------------------------------------------------- -#ifdef Spectral -#include "kdtree2.f90" -#endif -module math +module math use, intrinsic :: iso_c_binding use prec, only: pReal,pInt - + implicit none real(pReal), parameter, public :: PI = 3.14159265358979323846264338327950288419716939937510_pReal real(pReal), parameter, public :: INDEG = 180.0_pReal/pi @@ -53,17 +53,17 @@ module math 1_pInt,2_pInt, & 2_pInt,3_pInt, & 1_pInt,3_pInt & - ],[2,6]) !< Mandel notation + ],[2,6]) !< arrangement in Mandel notation real(pReal), dimension(6), parameter, private :: & nrmMandel = [& 1.0_pReal, 1.0_pReal, 1.0_pReal,& - 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal] - + 1.414213562373095_pReal, 1.414213562373095_pReal, 1.414213562373095_pReal] !< weighting for Mandel notation (forward) + real(pReal), dimension(6), parameter , public :: & invnrmMandel = [& 1.0_pReal, 1.0_pReal, 1.0_pReal,& - 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal] + 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal] !< weighting for Mandel notation (backward) integer(pInt), dimension (2,6), parameter, private :: & mapVoigt = reshape([& @@ -73,13 +73,12 @@ module math 2_pInt,3_pInt, & 1_pInt,3_pInt, & 1_pInt,2_pInt & - ],[2,6]) !< Voigt notation + ],[2,6]) !< arrangement in Voigt notation - real(pReal), dimension(6), parameter, private :: & - nrmVoigt = 1.0_pReal, & - invnrmVoigt = 1.0_pReal + real(pReal), dimension(6), parameter, private :: & + nrmVoigt = 1.0_pReal, & !< weighting for Voigt notation (forward) + invnrmVoigt = 1.0_pReal !< weighting for Voigt notation (backward) -! *** Plain notation *** integer(pInt), dimension (2,9), parameter, private :: & mapPlain = reshape([& 1_pInt,1_pInt, & @@ -91,13 +90,11 @@ module math 3_pInt,1_pInt, & 3_pInt,2_pInt, & 3_pInt,3_pInt & - ],[2,9]) + ],[2,9]) !< arrangement in Plain notation -! Symmetry operations as quaternions -! 24 for cubic, 12 for hexagonal = 36 integer(pInt), dimension(2), parameter, private :: & - math_NsymOperations = [24_pInt,12_pInt] - + math_NsymOperations = [24_pInt,12_pInt] !< Symmetry operations as quaternions 24 for cubic, 12 for hexagonal = 36 + real(pReal), dimension(4,36), parameter, private :: & math_symOperations = reshape([& 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations @@ -153,7 +150,7 @@ real(pReal), dimension(4,36), parameter, private :: & Gauss contains - + !-------------------------------------------------------------------------------------------------- !> @brief initialization of random seed generator !-------------------------------------------------------------------------------------------------- @@ -163,26 +160,27 @@ subroutine math_init use prec, only: tol_math_check use numerics, only: fixedSeed use IO, only: IO_error - + implicit none integer(pInt) :: i real(pReal), dimension(3,3) :: R,R2 real(pReal), dimension(3) :: Eulers real(pReal), dimension(4) :: q,q2,axisangle,randTest ! the following variables are system dependend and shound NOT be pInt - integer :: randSize ! gfortran requires a variable length to compile + integer :: randSize ! gfortran requires a variable length to compile integer, dimension(:), allocatable :: randInit ! if recalculations of former randomness (with given seed) is necessary ! comment the first random_seed call out, set randSize to 1, and use ifort character(len=64) :: error_msg - + !$OMP CRITICAL (write2out) write(6,*) '' write(6,*) '<<<+- math init -+>>>' write(6,*) '$Id$' #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - + call random_seed(size=randSize) + if (allocated(randInit)) deallocate(randInit) allocate(randInit(randSize)) if (fixedSeed > 0_pInt) then randInit(1:randSize) = int(fixedSeed) ! fixedSeed is of type pInt, randInit not @@ -198,13 +196,12 @@ subroutine math_init enddo !$OMP CRITICAL (write2out) - ! this critical block did cause trouble at IWM write(6,*) 'value of random seed: ', randInit(1) write(6,*) 'size of random seed: ', randSize write(6,'(a,4(/,26x,f17.14))') ' start of random sequence: ', randTest write(6,*) '' !$OMP END CRITICAL (write2out) - + call random_seed(put=randInit) call random_seed(get=randInit) @@ -221,8 +218,8 @@ subroutine math_init any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) call IO_error(401_pInt,ext_msg=error_msg) - endif - + endif + ! +++ q -> R -> q +++ R = math_QuaternionToR(q); q2 = math_RToQuaternion(R) @@ -230,8 +227,8 @@ subroutine math_init any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) call IO_error(402_pInt,ext_msg=error_msg) - endif - + endif + ! +++ q -> euler -> q +++ Eulers = math_QuaternionToEuler(q); q2 = math_EulerToQuaternion(Eulers) @@ -239,7 +236,7 @@ subroutine math_init any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) call IO_error(403_pInt,ext_msg=error_msg) - endif + endif ! +++ R -> euler -> R +++ Eulers = math_RToEuler(R); @@ -247,17 +244,17 @@ subroutine math_init if ( any(abs( R-R2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2)) call IO_error(404_pInt,ext_msg=error_msg) - endif - + endif + end subroutine math_init -!************************************************************************** -! Quicksort algorithm for two-dimensional integer arrays -! + +!-------------------------------------------------------------------------------------------------- +!> @brief 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 @@ -270,13 +267,13 @@ recursive subroutine qsort(a, istart, iend) call qsort(a, istart, ipivot-1_pInt) call qsort(a, ipivot+1_pInt, iend) endif - + end subroutine qsort -!************************************************************************** -! Partitioning required for quicksort -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Partitioning required for quicksort +!-------------------------------------------------------------------------------------------------- integer(pInt) function math_partition(a, istart, iend) implicit none @@ -318,12 +315,12 @@ integer(pInt) function math_partition(a, istart, iend) enddo end function math_partition - -!************************************************************************** -! range of integers starting at one -!************************************************************************** -pure function math_range(N) + +!-------------------------------------------------------------------------------------------------- +!> @brief range of integers starting at one +!-------------------------------------------------------------------------------------------------- +pure function math_range(N) implicit none integer(pInt), intent(in) :: N @@ -335,28 +332,28 @@ pure function math_range(N) end function math_range -!************************************************************************** -! second rank identity tensor of specified dimension -!************************************************************************** -pure function math_identity2nd(dimen) +!-------------------------------------------------------------------------------------------------- +!> @brief second rank identity tensor of specified dimension +!-------------------------------------------------------------------------------------------------- +pure function math_identity2nd(dimen) 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_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal + math_identity2nd = 0.0_pReal + forall (i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal end function math_identity2nd -!************************************************************************** -! permutation tensor e_ijk used for computing cross product of two tensors +!-------------------------------------------------------------------------------------------------- +!> @brief 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) implicit none @@ -374,11 +371,11 @@ pure function math_civita(i,j,k) end function math_civita -!************************************************************************** -! kronecker delta function d_ij +!-------------------------------------------------------------------------------------------------- +!> @brief kronecker delta function d_ij ! d_ij = 1 if i = j ! d_ij = 0 otherwise -!************************************************************************** +!-------------------------------------------------------------------------------------------------- pure function math_delta(i,j) implicit none @@ -391,10 +388,10 @@ pure function math_delta(i,j) end function math_delta -!************************************************************************** -! fourth rank identity tensor of specified dimension -!************************************************************************** -pure function math_identity4th(dimen) +!-------------------------------------------------------------------------------------------------- +!> @brief fourth rank identity tensor of specified dimension +!-------------------------------------------------------------------------------------------------- +pure function math_identity4th(dimen) implicit none integer(pInt), intent(in) :: dimen @@ -402,15 +399,15 @@ pure function math_identity4th(dimen) real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th forall (i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt: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)) + 0.5_pReal*(math_I3(i,k)*math_I3(j,k)+math_I3(i,l)*math_I3(j,k)) end function math_identity4th - -!************************************************************************** -! vector product a x b -!************************************************************************** -pure function math_vectorproduct(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief vector product a x b +!-------------------------------------------------------------------------------------------------- +pure function math_vectorproduct(A,B) implicit none real(pReal), dimension(3), intent(in) :: A,B @@ -423,26 +420,26 @@ pure function math_vectorproduct(A,B) end function math_vectorproduct -!************************************************************************** -! tensor product a \otimes b -!************************************************************************** -pure function math_tensorproduct(A,B) +!-------------------------------------------------------------------------------------------------- +!> @brief tensor product a \otimes b +!-------------------------------------------------------------------------------------------------- +pure function math_tensorproduct(A,B) implicit none real(pReal), dimension(3), intent(in) :: A,B real(pReal), dimension(3,3) :: math_tensorproduct integer(pInt) :: i,j - + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_tensorproduct(i,j) = A(i)*B(j) end function math_tensorproduct -!************************************************************************** -! matrix multiplication 3x3 = 1 -!************************************************************************** -pure function math_mul3x3(A,B) +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 3x3 = 1 +!-------------------------------------------------------------------------------------------------- +pure function math_mul3x3(A,B) implicit none @@ -457,10 +454,10 @@ pure function math_mul3x3(A,B) end function math_mul3x3 -!************************************************************************** -! matrix multiplication 6x6 = 1 -!************************************************************************** -pure function math_mul6x6(A,B) +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 6x6 = 1 +!-------------------------------------------------------------------------------------------------- +pure function math_mul6x6(A,B) implicit none @@ -474,11 +471,11 @@ pure function math_mul6x6(A,B) end function math_mul6x6 - -!************************************************************************** -! matrix multiplication 33x33 = 1 (double contraction --> ij * ij) -!************************************************************************** -pure function math_mul33xx33(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 33x33 = 1 (double contraction --> ij * ij) +!-------------------------------------------------------------------------------------------------- +pure function math_mul33xx33(A,B) implicit none @@ -492,11 +489,11 @@ pure function math_mul33xx33(A,B) end function math_mul33xx33 - -!************************************************************************** -! matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij) -!************************************************************************** -pure function math_mul3333xx33(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij) +!-------------------------------------------------------------------------------------------------- +pure function math_mul3333xx33(A,B) implicit none @@ -511,10 +508,10 @@ pure function math_mul3333xx33(A,B) end function math_mul3333xx33 -!************************************************************************** -! matrix multiplication 3333x3333 = 3333 (ijkl *klmn = ijmn) -!************************************************************************** -pure function math_mul3333xx3333(A,B) +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 3333x3333 = 3333 (ijkl *klmn = ijmn) +!-------------------------------------------------------------------------------------------------- +pure function math_mul3333xx3333(A,B) implicit none integer(pInt) :: i,j,k,l @@ -530,12 +527,12 @@ pure function math_mul3333xx3333(A,B) enddo; enddo; enddo; enddo end function math_mul3333xx3333 - -!************************************************************************** -! matrix multiplication 33x33 = 33 -!************************************************************************** -pure function math_mul33x33(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 33x33 = 33 +!-------------------------------------------------------------------------------------------------- +pure function math_mul33x33(A,B) implicit none integer(pInt) :: i,j @@ -548,10 +545,10 @@ pure function math_mul33x33(A,B) end function math_mul33x33 -!************************************************************************** -! matrix multiplication 66x66 = 66 -!************************************************************************** -pure function math_mul66x66(A,B) +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 66x66 = 66 +!-------------------------------------------------------------------------------------------------- +pure function math_mul66x66(A,B) implicit none integer(pInt) :: i,j @@ -564,11 +561,11 @@ pure function math_mul66x66(A,B) end function math_mul66x66 - -!************************************************************************** -! matrix multiplication 99x99 = 99 -!************************************************************************** -pure function math_mul99x99(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 99x99 = 99 +!-------------------------------------------------------------------------------------------------- +pure function math_mul99x99(A,B) use prec, only: pReal, pInt @@ -586,11 +583,11 @@ pure function math_mul99x99(A,B) end function math_mul99x99 - -!************************************************************************** -! matrix multiplication 33x3 = 3 -!************************************************************************** -pure function math_mul33x3(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 33x3 = 3 +!-------------------------------------------------------------------------------------------------- +pure function math_mul33x3(A,B) implicit none integer(pInt) :: i @@ -601,11 +598,12 @@ pure function math_mul33x3(A,B) forall (i=1_pInt:3_pInt) math_mul33x3(i) = sum(A(i,1:3)*B) end function math_mul33x3 - - !************************************************************************** -! matrix multiplication complex(33) x real(3) = complex(3) -!************************************************************************** -pure function math_mul33x3_complex(A,B) + + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication complex(33) x real(3) = complex(3) +!-------------------------------------------------------------------------------------------------- +pure function math_mul33x3_complex(A,B) implicit none integer(pInt) :: i @@ -617,11 +615,11 @@ pure function math_mul33x3_complex(A,B) end function math_mul33x3_complex - -!************************************************************************** -! matrix multiplication 66x6 = 6 -!************************************************************************** -pure function math_mul66x6(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief matrix multiplication 66x6 = 6 +!-------------------------------------------------------------------------------------------------- +pure function math_mul66x6(A,B) implicit none @@ -636,16 +634,16 @@ pure function math_mul66x6(A,B) end function math_mul66x6 - -!************************************************************************** -! random quaternion -!************************************************************************** -function math_qRnd() + +!-------------------------------------------------------------------------------------------------- +!> @brief random quaternion +!-------------------------------------------------------------------------------------------------- +function math_qRnd() implicit none real(pReal), dimension(4) :: math_qRnd real(pReal), dimension(3) :: rnd - + call halton(3_pInt,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)) @@ -654,11 +652,11 @@ function math_qRnd() end function math_qRnd - -!************************************************************************** -! quaternion multiplication q1xq2 = q12 -!************************************************************************** -pure function math_qMul(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion multiplication q1xq2 = q12 +!-------------------------------------------------------------------------------------------------- +pure function math_qMul(A,B) implicit none real(pReal), dimension(4), intent(in) :: A, B @@ -671,11 +669,11 @@ pure function math_qMul(A,B) end function math_qMul - -!************************************************************************** -! quaternion dotproduct -!************************************************************************** -pure function math_qDot(A,B) + +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion dotproduct +!-------------------------------------------------------------------------------------------------- +pure function math_qDot(A,B) implicit none real(pReal), dimension(4), intent(in) :: A, B @@ -685,11 +683,11 @@ pure function math_qDot(A,B) end function math_qDot - -!************************************************************************** -! quaternion conjugation -!************************************************************************** -pure function math_qConj(Q) + +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion conjugation +!-------------------------------------------------------------------------------------------------- +pure function math_qConj(Q) implicit none real(pReal), dimension(4), intent(in) :: Q @@ -700,44 +698,44 @@ pure function math_qConj(Q) end function math_qConj - -!************************************************************************** -! quaternion norm -!************************************************************************** -pure function math_qNorm(Q) + +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion norm +!-------------------------------------------------------------------------------------------------- +pure function math_qNorm(Q) 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))) end function math_qNorm -!************************************************************************** -! quaternion inversion -!************************************************************************** -pure function math_qInv(Q) +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion inversion +!-------------------------------------------------------------------------------------------------- +pure function math_qInv(Q) 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 - + end function math_qInv - -!************************************************************************** -! action of a quaternion on a vector (rotate vector v with Q) -!************************************************************************** -pure function math_qRot(Q,v) + +!-------------------------------------------------------------------------------------------------- +!> @brief action of a quaternion on a vector (rotate vector v with Q) +!-------------------------------------------------------------------------------------------------- +pure function math_qRot(Q,v) implicit none real(pReal), dimension(4), intent(in) :: Q @@ -745,51 +743,50 @@ pure function math_qRot(Q,v) real(pReal), dimension(3) :: math_qRot real(pReal), dimension(4,4) :: T integer(pInt) :: i, j - + do i = 1_pInt,4_pInt do j = 1_pInt,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 end function math_qRot - -!************************************************************************** -! transposition of a 33 matrix -!************************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief transposition of a 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_transpose33(A) implicit none real(pReal),dimension(3,3),intent(in) :: A real(pReal),dimension(3,3) :: math_transpose33 integer(pInt) :: i,j - + forall(i=1_pInt:3_pInt, j=1_pInt:3_pInt) math_transpose33(i,j) = A(j,i) end function math_transpose33 - -!************************************************************************** -! Cramer inversion of 33 matrix (function) -!************************************************************************** -pure function math_inv33(A) + +!-------------------------------------------------------------------------------------------------- +!> @brief Cramer inversion of 33 matrix (function) ! direct Cramer inversion of matrix A. ! returns all zeroes if not possible, i.e. if det close to zero +!-------------------------------------------------------------------------------------------------- +pure function math_inv33(A) implicit none - real(pReal),dimension(3,3),intent(in) :: A real(pReal) :: DetA real(pReal),dimension(3,3) :: math_inv33 - + math_inv33 = 0.0_pReal DetA = A(1,1) * (A(2,2) * A(3,3) - A(2,3) * A(3,2))& @@ -813,9 +810,9 @@ pure function math_inv33(A) end function math_inv33 -!************************************************************************** -! Cramer inversion of 33 matrix (subroutine) -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Cramer inversion of 33 matrix (subroutine) +!-------------------------------------------------------------------------------------------------- pure subroutine math_invert33(A, InvA, DetA, error) ! Bestimmung der Determinanten und Inversen einer 33-Matrix @@ -849,36 +846,36 @@ pure subroutine math_invert33(A, InvA, DetA, error) 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 end subroutine math_invert33 -!************************************************************************** -! Inversion of symmetriced 3x3x3x3 tensor. -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Inversion of symmetriced 3x3x3x3 tensor. +!-------------------------------------------------------------------------------------------------- function math_invSym3333(A) use IO, only: IO_error - + implicit none real(pReal),dimension(3,3,3,3) :: math_invSym3333 - + real(pReal),dimension(3,3,3,3),intent(in) :: A integer(pInt) :: ierr1, ierr2 integer(pInt), dimension(6) :: ipiv6 real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6) :: work6 - + temp66_real = math_Mandel3333to66(A) call dgetrf(6,6,temp66_real,6,ipiv6,ierr1) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr2) if (ierr1*ierr2 == 0_pInt) then math_invSym3333 = math_Mandel66to3333(temp66_real) - else + else call IO_error(400_pInt, ext_msg = 'math_invSym3333') endif @@ -1069,9 +1066,9 @@ pure subroutine Gauss (dimen,A,B,LogAbsDetA,NegHDK,error) end subroutine Gauss -!******************************************************************** -! symmetrize a 33 matrix -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief symmetrize a 33 matrix +!-------------------------------------------------------------------------------------------------- function math_symmetric33(m) implicit none @@ -1083,11 +1080,11 @@ function math_symmetric33(m) forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_symmetric33(i,j) = 0.5_pReal * (m(i,j) + m(j,i)) end function math_symmetric33 - -!******************************************************************** -! symmetrize a 66 matrix -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief symmetrize a 66 matrix +!-------------------------------------------------------------------------------------------------- pure function math_symmetric66(m) implicit none @@ -1095,15 +1092,15 @@ pure function math_symmetric66(m) integer(pInt) :: i,j real(pReal), dimension(6,6), intent(in) :: m real(pReal), dimension(6,6) :: math_symmetric66 - + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_symmetric66(i,j) = 0.5_pReal * (m(i,j) + m(j,i)) end function math_symmetric66 - -!******************************************************************** -! skew part of a 33 matrix -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief skew part of a 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_skew33(m) implicit none @@ -1111,15 +1108,15 @@ pure function math_skew33(m) real(pReal), dimension(3,3) :: math_skew33 real(pReal), dimension(3,3), intent(in) :: m integer(pInt) :: i,j - + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_skew33(i,j) = m(i,j) - 0.5_pReal * (m(i,j) + m(j,i)) end function math_skew33 - -!******************************************************************** -! deviatoric part of a 33 matrix -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief deviatoric part of a 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_deviatoric33(m) implicit none @@ -1136,9 +1133,9 @@ pure function math_deviatoric33(m) end function math_deviatoric33 -!******************************************************************** -! equivalent scalar quantity of a full strain tensor -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief equivalent scalar quantity of a full strain tensor +!-------------------------------------------------------------------------------------------------- pure function math_equivStrain33(m) implicit none @@ -1207,10 +1204,10 @@ pure function math_det33(m) end function math_det33 - -!******************************************************************** -! norm of a 33 matrix -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief norm of a 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_norm33(m) implicit none @@ -1222,10 +1219,10 @@ pure function math_norm33(m) end function - -!******************************************************************** -! euclidic norm of a 3 vector -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief euclidic norm of a 3 vector +!-------------------------------------------------------------------------------------------------- pure function math_norm3(v) implicit none @@ -1234,13 +1231,13 @@ pure function math_norm3(v) real(pReal) :: math_norm3 math_norm3 = sqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) - + end function math_norm3 - -!******************************************************************** -! convert 33 matrix into vector 9 -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief convert 33 matrix into vector 9 +!-------------------------------------------------------------------------------------------------- pure function math_Plain33to9(m33) implicit none @@ -1248,15 +1245,15 @@ pure function math_Plain33to9(m33) real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(9) :: math_Plain33to9 integer(pInt) :: i - + forall (i=1_pInt:9_pInt) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) end function math_Plain33to9 - - -!******************************************************************** -! convert Plain 9 back to 33 matrix -!******************************************************************** + + +!-------------------------------------------------------------------------------------------------- +!> @brief convert Plain 9 back to 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_Plain9to33(v9) implicit none @@ -1264,15 +1261,16 @@ pure function math_Plain9to33(v9) real(pReal), dimension(9), intent(in) :: v9 real(pReal), dimension(3,3) :: math_Plain9to33 integer(pInt) :: i - + forall (i=1_pInt:9_pInt) math_Plain9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) end function math_Plain9to33 - -!******************************************************************** -! convert symmetric 33 matrix into Mandel vector 6 -!******************************************************************** + + +!-------------------------------------------------------------------------------------------------- +!> @brief convert symmetric 33 matrix into Mandel vector 6 +!-------------------------------------------------------------------------------------------------- pure function math_Mandel33to6(m33) implicit none @@ -1280,15 +1278,15 @@ pure function math_Mandel33to6(m33) real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(6) :: math_Mandel33to6 integer(pInt) :: i - + forall (i=1_pInt:6_pInt) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i)) end function math_Mandel33to6 -!******************************************************************** -! convert Mandel 6 back to symmetric 33 matrix -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert Mandel 6 back to symmetric 33 matrix +!-------------------------------------------------------------------------------------------------- pure function math_Mandel6to33(v6) implicit none @@ -1296,7 +1294,7 @@ pure function math_Mandel6to33(v6) real(pReal), dimension(6), intent(in) :: v6 real(pReal), dimension(3,3) :: math_Mandel6to33 integer(pInt) :: i - + forall (i=1_pInt:6_pInt) math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i) math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i) @@ -1305,9 +1303,9 @@ pure function math_Mandel6to33(v6) end function math_Mandel6to33 -!******************************************************************** -! convert 3333 tensor into plain matrix 99 -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert 3333 tensor into plain matrix 99 +!-------------------------------------------------------------------------------------------------- pure function math_Plain3333to99(m3333) implicit none @@ -1315,15 +1313,15 @@ pure function math_Plain3333to99(m3333) real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(9,9) :: math_Plain3333to99 integer(pInt) :: i,j - + forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain3333to99(i,j) = & m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) end function math_Plain3333to99 - -!******************************************************************** -! plain matrix 99 into 3333 tensor -!******************************************************************** + +!-------------------------------------------------------------------------------------------------- +!> @brief plain matrix 99 into 3333 tensor +!-------------------------------------------------------------------------------------------------- pure function math_Plain99to3333(m99) implicit none @@ -1331,16 +1329,16 @@ pure function math_Plain99to3333(m99) real(pReal), dimension(9,9), intent(in) :: m99 real(pReal), dimension(3,3,3,3) :: math_Plain99to3333 integer(pInt) :: i,j - + forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),& mapPlain(1,j),mapPlain(2,j)) = m99(i,j) end function math_Plain99to3333 -!******************************************************************** -! convert Mandel matrix 66 into Plain matrix 66 -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert Mandel matrix 66 into Plain matrix 66 +!-------------------------------------------------------------------------------------------------- pure function math_Mandel66toPlain66(m66) implicit none @@ -1348,7 +1346,7 @@ pure function math_Mandel66toPlain66(m66) real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6) :: math_Mandel66toPlain66 integer(pInt) :: i,j - + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) & math_Mandel66toPlain66(i,j) = invnrmMandel(i) * invnrmMandel(j) * m66(i,j) return @@ -1356,9 +1354,9 @@ pure function math_Mandel66toPlain66(m66) end function -!******************************************************************** -! convert Plain matrix 66 into Mandel matrix 66 -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert Plain matrix 66 into Mandel matrix 66 +!-------------------------------------------------------------------------------------------------- pure function math_Plain66toMandel66(m66) implicit none @@ -1366,7 +1364,7 @@ pure function math_Plain66toMandel66(m66) real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6) :: math_Plain66toMandel66 integer(pInt) i,j - + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) & math_Plain66toMandel66(i,j) = nrmMandel(i) * nrmMandel(j) * m66(i,j) return @@ -1374,9 +1372,9 @@ pure function math_Plain66toMandel66(m66) end function -!******************************************************************** -! convert symmetric 3333 tensor into Mandel matrix 66 -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert symmetric 3333 tensor into Mandel matrix 66 +!-------------------------------------------------------------------------------------------------- pure function math_Mandel3333to66(m3333) implicit none @@ -1384,16 +1382,16 @@ pure function math_Mandel3333to66(m3333) real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(6,6) :: math_Mandel3333to66 integer(pInt) :: i,j - + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_Mandel3333to66(i,j) = & nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) end function math_Mandel3333to66 -!******************************************************************** -! convert Mandel matrix 66 back to symmetric 3333 tensor -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert Mandel matrix 66 back to symmetric 3333 tensor +!-------------------------------------------------------------------------------------------------- pure function math_Mandel66to3333(m66) implicit none @@ -1401,8 +1399,8 @@ pure function math_Mandel66to3333(m66) real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 integer(pInt) :: i,j - - forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) + + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) 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) @@ -1412,9 +1410,9 @@ pure function math_Mandel66to3333(m66) end function math_Mandel66to3333 -!******************************************************************** -! convert Voigt matrix 66 back to symmetric 3333 tensor -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief convert Voigt matrix 66 back to symmetric 3333 tensor +!-------------------------------------------------------------------------------------------------- pure function math_Voigt66to3333(m66) implicit none @@ -1422,8 +1420,8 @@ pure function math_Voigt66to3333(m66) real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 integer(pInt) :: i,j - - forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) + + forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) 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) @@ -1433,9 +1431,9 @@ pure function math_Voigt66to3333(m66) end function math_Voigt66to3333 -!******************************************************************** -! Euler angles (in radians) from rotation matrix -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Euler angles (in radians) from rotation matrix +!-------------------------------------------------------------------------------------------------- pure function math_RtoEuler(R) implicit none @@ -1449,10 +1447,10 @@ pure function math_RtoEuler(R) sqhk=sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3)) ! calculate PHI myVal=R(3,3)/sqhkl - + if(myVal > 1.0_pReal) myVal = 1.0_pReal if(myVal < -1.0_pReal) myVal = -1.0_pReal - + math_RtoEuler(2) = acos(myVal) if(math_RtoEuler(2) < 1.0e-8_pReal) then @@ -1462,7 +1460,7 @@ pure function math_RtoEuler(R) myVal=R(1,1)/squvw if(myVal > 1.0_pReal) myVal = 1.0_pReal if(myVal < -1.0_pReal) myVal = -1.0_pReal - + math_RtoEuler(1) = acos(myVal) if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) else @@ -1470,24 +1468,24 @@ pure function math_RtoEuler(R) myVal=R(2,3)/sqhk if(myVal > 1.0_pReal) myVal = 1.0_pReal if(myVal < -1.0_pReal) myVal = -1.0_pReal - + math_RtoEuler(3) = acos(myVal) if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3) ! calculate phi1 myVal=-R(3,2)/sin(math_RtoEuler(2)) if(myVal > 1.0_pReal) myVal = 1.0_pReal if(myVal < -1.0_pReal) myVal = -1.0_pReal - + math_RtoEuler(1) = acos(myVal) if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1) end if - + end function math_RtoEuler -!******************************************************************** -! quaternion (w+ix+jy+kz) from orientation matrix -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion (w+ix+jy+kz) from orientation matrix +!-------------------------------------------------------------------------------------------------- ! math adopted from http://code.google.com/p/mtex/source/browse/trunk/geometry/geometry_tools/mat2quat.m pure function math_RtoQuaternion(R) @@ -1506,7 +1504,7 @@ pure function math_RtoQuaternion(R) largest = maxloc(absQ) - max_absQ=0.5_pReal * sqrt(absQ(largest(1))) + max_absQ=0.5_pReal * sqrt(absQ(largest(1))) select case(largest(1)) case (1_pInt) @@ -1514,19 +1512,19 @@ pure function math_RtoQuaternion(R) 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_pInt) math_RtoQuaternion(1) = R(2,3)-R(3,2) !2---------------------------------- math_RtoQuaternion(3) = R(1,2)+R(2,1) math_RtoQuaternion(4) = R(3,1)+R(1,3) - + case (3_pInt) math_RtoQuaternion(1) = R(3,1)-R(1,3) math_RtoQuaternion(2) = R(1,2)+R(2,1) !3---------------------------------- math_RtoQuaternion(4) = R(2,3)+R(3,2) - + case (4_pInt) math_RtoQuaternion (1) = R(1,2)-R(2,1) math_RtoQuaternion (2) = R(3,1)+R(1,3) @@ -1536,13 +1534,13 @@ pure function math_RtoQuaternion(R) math_RtoQuaternion = math_RtoQuaternion*0.25_pReal/max_absQ math_RtoQuaternion(largest(1)) = max_absQ - + end function math_RtoQuaternion -!**************************************************************** -! rotation matrix from Euler angles (in radians) -!**************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief rotation matrix from Euler angles (in radians) +!-------------------------------------------------------------------------------------------------- pure function math_EulerToR(Euler) implicit none @@ -1567,13 +1565,13 @@ pure function math_EulerToR(Euler) math_EulerToR(3,1)=S1*S math_EulerToR(3,2)=-C1*S math_EulerToR(3,3)=C - -end function math_EulerToR - -!******************************************************************** -! quaternion (w+ix+jy+kz) from 3-1-3 Euler angles (in radians) -!******************************************************************** +end function math_EulerToR + + +!-------------------------------------------------------------------------------------------------- +!> @brief quaternion (w+ix+jy+kz) from 3-1-3 Euler angles (in radians) +!-------------------------------------------------------------------------------------------------- pure function math_EulerToQuaternion(eulerangles) implicit none @@ -1582,23 +1580,23 @@ pure function math_EulerToQuaternion(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 - + end function math_EulerToQuaternion -!**************************************************************** -! rotation matrix from axis and angle (in radians) -!**************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief rotation matrix from axis and angle (in radians) +!-------------------------------------------------------------------------------------------------- pure function math_AxisAngleToR(axis,omega) implicit none @@ -1617,18 +1615,18 @@ pure function math_AxisAngleToR(axis,omega) 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.0_pReal - math_AxisAngleToR(2,1) = -s*axisNrm(3) + c1*axisNrm(1)*axisNrm(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.0_pReal 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.0_pReal @@ -1669,9 +1667,9 @@ pure function math_AxisAngleToQuaternion(axis,omega) end function math_AxisAngleToQuaternion -!******************************************************************** -! orientation matrix from quaternion (w+ix+jy+kz) -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief orientation matrix from quaternion (w+ix+jy+kz) +!-------------------------------------------------------------------------------------------------- pure function math_QuaternionToR(Q) implicit none @@ -1679,7 +1677,7 @@ pure function math_QuaternionToR(Q) real(pReal), dimension(4), intent(in) :: Q real(pReal), dimension(3,3) :: math_QuaternionToR, T,S integer(pInt) :: i, j - + forall (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & T(i,j) = Q(i+1_pInt) * Q(j+1_pInt) S = reshape( (/0.0_pReal, Q(4), -Q(3), & @@ -1693,9 +1691,9 @@ pure function math_QuaternionToR(Q) end function math_QuaternionToR -!******************************************************************** -! 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz) -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz) +!-------------------------------------------------------------------------------------------------- pure function math_QuaternionToEuler(Q) implicit none @@ -1708,8 +1706,8 @@ pure function math_QuaternionToEuler(Q) if (abs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then acos_arg=Q(1) - if(acos_arg > 1.0_pReal)acos_arg = 1.0_pReal - if(acos_arg < -1.0_pReal)acos_arg = -1.0_pReal + if(acos_arg > 1.0_pReal)acos_arg = 1.0_pReal + if(acos_arg < -1.0_pReal)acos_arg = -1.0_pReal math_QuaternionToEuler(1) = 2.0_pReal*acos(acos_arg) math_QuaternionToEuler(3) = 0.0_pReal else @@ -1728,20 +1726,20 @@ pure function math_QuaternionToEuler(Q) end function math_QuaternionToEuler -!******************************************************************** -! axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) +!-------------------------------------------------------------------------------------------------- pure function math_QuaternionToAxisAngle(Q) implicit none real(pReal), dimension(4), intent(in) :: Q real(pReal) :: halfAngle, sinHalfAngle - real(pReal), dimension(4) :: math_QuaternionToAxisAngle + 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 @@ -1772,9 +1770,9 @@ pure function math_QuaternionToRodrig(Q) end function math_QuaternionToRodrig -!************************************************************************** -! misorientation angle between two sets of Euler angles -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief misorientation angle between two sets of Euler angles +!-------------------------------------------------------------------------------------------------- pure function math_EulerMisorientation(EulerA,EulerB) implicit none @@ -1791,23 +1789,23 @@ pure function math_EulerMisorientation(EulerA,EulerB) end function math_EulerMisorientation -!************************************************************************** -! figures whether unit quat falls into stereographic standard triangle -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief figures whether unit quat falls into stereographic standard triangle +!-------------------------------------------------------------------------------------------------- pure function math_QuaternionInSST(Q, symmetryType) implicit none - !*** input variables + !*** 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) if (any(Rodrig/=Rodrig)) then math_QuaternionInSST = .false. @@ -1825,38 +1823,38 @@ pure function math_QuaternionInSST(Q, symmetryType) math_QuaternionInSST = .true. end select endif - + end function math_QuaternionInSST -!************************************************************************** -! calculates the disorientation for 2 unit quaternions -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the disorientation for 2 unit quaternions +!-------------------------------------------------------------------------------------------------- function math_QuaternionDisorientation(Q1, Q2, symmetryType) use IO, only: IO_error implicit none - - !*** input variables + + !*** 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_pInt) if (math_QuaternionDisorientation(1) < 0.0_pReal) & math_QuaternionDisorientation = -math_QuaternionDisorientation ! keep omega within 0 to 180 deg - + case (1_pInt,2_pInt) s = sum(math_NsymOperations(1:symmetryType-1_pInt)) do i = 1_pInt,2_pInt @@ -1871,17 +1869,17 @@ function math_QuaternionDisorientation(Q1, Q2, symmetryType) math_QuaternionInSST(mis,symmetryType)) & math_QuaternionDisorientation = mis ! found better one enddo; enddo; enddo - + case default call IO_error(450_pInt,symmetryType) ! complain about unknown symmetry end select - + end function math_QuaternionDisorientation -!******************************************************************** -! draw a random sample from Euler space -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief draw a random sample from Euler space +!-------------------------------------------------------------------------------------------------- function math_sampleRandomOri() implicit none @@ -1896,10 +1894,9 @@ function math_sampleRandomOri() end function math_sampleRandomOri -!******************************************************************** -! draw a random sample from Gauss component -! with noise (in radians) half-width -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief draw a random sample from Gauss component with noise (in radians) half-width +!-------------------------------------------------------------------------------------------------- function math_sampleGaussOri(center,noise) implicit none @@ -1922,17 +1919,17 @@ endif do call halton(5_pInt,rnd) - forall (i=1_pInt:3_pInt) rnd(i) = 2.0_pReal*rnd(i)-1.0_pReal ! expand 1:3 to range [-1,+1] + forall (i=1_pInt:3_pInt) 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_pReal)) exit + if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(origin,disturb)/scatter)**2_pReal)) exit enddo math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) - + end function math_sampleGaussOri - + !-------------------------------------------------------------------------------------------------- !> @brief draw a random sample from Fiber component with noise (in radians) @@ -1982,7 +1979,6 @@ function math_sampleFiberOri(alpha,beta,noise) ! ---# rotation about random axis perpend to fiber #--- ! random axis pependicular to fiber axis - axis(1:2) = rnd(2:3) if (fiberInS(3) /= 0.0_pReal) then axis(3)=-(axis(1)*fiberInS(1)+axis(2)*fiberInS(2))/fiberInS(3) @@ -2003,7 +1999,6 @@ function math_sampleFiberOri(alpha,beta,noise) exit end if enddo - if (rnd(6) <= 0.5) angle = -angle pRot = math_AxisAngleToR(axis,angle) @@ -2014,10 +2009,9 @@ function math_sampleFiberOri(alpha,beta,noise) end function math_sampleFiberOri -!******************************************************************** -! symmetric Euler angles for given symmetry string -! 'triclinic' or '', 'monoclinic', 'orthotropic' -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief symmetric Euler angles for given symmetry string 'triclinic' or '', 'monoclinic', 'orthotropic' +!-------------------------------------------------------------------------------------------------- pure function math_symmetricEulers(sym,Euler) implicit none @@ -2026,7 +2020,7 @@ pure function math_symmetricEulers(sym,Euler) 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) @@ -2096,12 +2090,12 @@ math_sampleGaussVar = scatter * stddev end function math_sampleGaussVar - -!**************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief not yet done +!-------------------------------------------------------------------------------------------------- subroutine math_spectralDecompositionSym33(M,values,vectors,error) -!**************************************************************** - implicit none + implicit none real(pReal), dimension(3,3), intent(in) :: M real(pReal), dimension(3), intent(out) :: values real(pReal), dimension(3,3), intent(out) :: vectors @@ -2109,18 +2103,19 @@ subroutine math_spectralDecompositionSym33(M,values,vectors,error) integer(pInt) :: info real(pReal), dimension((64+2)*3) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f - + vectors = M ! copy matrix to input (doubles as output) array call DSYEV('V','U',3,vectors,3,values,work,(64+2)*3,info) error = (info == 0_pInt) - + end subroutine -!**************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief FE = R.U +!-------------------------------------------------------------------------------------------------- pure subroutine math_pDecomposition(FE,U,R,error) -!-----FE = R.U -!**************************************************************** + implicit none real(pReal), intent(in), dimension(3,3) :: FE @@ -2140,9 +2135,10 @@ pure subroutine math_pDecomposition(FE,U,R,error) end subroutine math_pDecomposition -!********************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M +!-------------------------------------------------------------------------------------------------- pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) -!**** EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M implicit none @@ -2186,7 +2182,7 @@ pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) EW2=Y2-R/3.0_pReal EW3=Y3-R/3.0_pReal C1=ABS(EW1-EW2) - C2=ABS(EW2-EW3) + C2=ABS(EW2-EW3) C3=ABS(EW3-EW1) IF(C1.LT.TOL) THEN @@ -2212,7 +2208,7 @@ pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) EW3=0.0_pReal ELSE IF(C3.LT.TOL) THEN ! EW1 is equal to EW3 - D2=1.0_pReal/(EW2-EW1)/(EW2-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 @@ -2223,7 +2219,7 @@ pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) ELSE ! all three eigenvectors are different D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) - D2=1.0_pReal/(EW2-EW1)/(EW2-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 @@ -2238,9 +2234,10 @@ pure subroutine math_spectral1(M,EW1,EW2,EW3,EB1,EB2,EB3) end subroutine math_spectral1 -!********************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief Eigenvalues of symmetric 3X3 matrix M +!-------------------------------------------------------------------------------------------------- function math_eigenvalues33(M) -!**** Eigenvalues of symmetric 3X3 matrix M implicit none @@ -2256,7 +2253,7 @@ function math_eigenvalues33(M) 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 - + if((abs(P) < TOL) .and. (abs(Q) < TOL)) THEN ! three equivalent eigenvalues math_eigenvalues33(1) = HI1M/3.0_pReal @@ -2283,15 +2280,15 @@ function math_eigenvalues33(M) end function math_eigenvalues33 -!********************************************************************** -!**** HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M - +!-------------------------------------------------------------------------------------------------- +!> @brief HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M +!-------------------------------------------------------------------------------------------------- pure subroutine math_hi(M,HI1M,HI2M,HI3M) - + implicit none - real(pReal), intent(in) :: M(3,3) - real(pReal), intent(out) :: HI1M, HI2M, HI3M + 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.0_pReal/2.0_pReal- (M(1,1)**2.0_pReal+M(2,2)**2.0_pReal+M(3,3)**2.0_pReal)& @@ -2360,8 +2357,8 @@ subroutine get_seed(seed) end subroutine get_seed -!******************************************************************************* -! HALTON computes the next element in the Halton sequence. +!-------------------------------------------------------------------------------------------------- +!> @brief HALTON computes the next element in the Halton sequence. ! ! Parameters: ! Input, integer NDIM, the dimension of the element. @@ -2372,7 +2369,7 @@ end subroutine get_seed ! ! Modified: 29 April 2005 ! Author: Franz Roters -! +!-------------------------------------------------------------------------------------------------- subroutine halton(ndim, r) implicit none @@ -2395,8 +2392,8 @@ subroutine halton(ndim, r) end subroutine halton -!******************************************************************************* -! HALTON_MEMORY sets or returns quantities associated with the Halton sequence. +!-------------------------------------------------------------------------------------------------- +!> @brief HALTON_MEMORY sets or returns quantities associated with the Halton sequence. ! ! Parameters: ! Input, character (len = *) action_halton, the desired action. @@ -2507,8 +2504,8 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton) end subroutine halton_memory -!******************************************************************************* -! HALTON_NDIM_SET sets the dimension for a Halton sequence. +!-------------------------------------------------------------------------------------------------- +!> @brief HALTON_NDIM_SET sets the dimension for a Halton sequence. ! ! Parameters: ! Input, integer NDIM, the dimension of the Halton vectors. @@ -2518,7 +2515,7 @@ end subroutine halton_memory ! ! Modified: 29 April 2005 ! Author: Franz Roters -! +!-------------------------------------------------------------------------------------------------- subroutine halton_ndim_set (ndim) implicit none @@ -2531,19 +2528,19 @@ subroutine halton_ndim_set (ndim) end subroutine halton_ndim_set -!******************************************************************************* -! HALTON_SEED_SET sets the "seed" for the Halton sequence. + +!> HALTON_SEED_SET sets the "seed" for the Halton sequence. ! ! 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. ! ! Parameters: @@ -2554,7 +2551,7 @@ end subroutine halton_ndim_set ! ! Modified: 29 April 2005 ! Author: Franz Roters -! +!-------------------------------------------------------------------------------------------------- subroutine halton_seed_set (seed) implicit none @@ -2568,24 +2565,24 @@ subroutine halton_seed_set (seed) end subroutine halton_seed_set -!******************************************************************************* -! I_TO_HALTON computes an element of a Halton sequence. +!-------------------------------------------------------------------------------------------------- +!> @brief 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. -! +! ! 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. ! @@ -2626,35 +2623,35 @@ subroutine i_to_halton (seed, base, ndim, r) end subroutine i_to_halton -!******************************************************************************* -! PRIME returns any of the first PRIME_MAX prime numbers. +!-------------------------------------------------------------------------------------------------- +!> @brief PRIME returns any of the first PRIME_MAX prime numbers. ! ! Note: ! PRIME_MAX is 1500, and the largest prime stored is 12553. ! 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: 21 June 2002 ! Author: John Burkardt -! +! ! Modified: 29 April 2005 ! Author: Franz Roters -! +!-------------------------------------------------------------------------------------------------- function prime(n) - + use IO, only: IO_error implicit none @@ -2666,8 +2663,8 @@ function prime(n) if (icall == 0_pInt) then icall = 1_pInt - - npvec(1:100) = (/& + + npvec = [& 2_pInt, 3_pInt, 5_pInt, 7_pInt, 11_pInt, 13_pInt, 17_pInt, 19_pInt, 23_pInt, 29_pInt, & 31_pInt, 37_pInt, 41_pInt, 43_pInt, 47_pInt, 53_pInt, 59_pInt, 61_pInt, 67_pInt, 71_pInt, & 73_pInt, 79_pInt, 83_pInt, 89_pInt, 97_pInt, 101_pInt, 103_pInt, 107_pInt, 109_pInt, 113_pInt, & @@ -2677,9 +2674,8 @@ function prime(n) 283_pInt, 293_pInt, 307_pInt, 311_pInt, 313_pInt, 317_pInt, 331_pInt, 337_pInt, 347_pInt, 349_pInt, & 353_pInt, 359_pInt, 367_pInt, 373_pInt, 379_pInt, 383_pInt, 389_pInt, 397_pInt, 401_pInt, 409_pInt, & 419_pInt, 421_pInt, 431_pInt, 433_pInt, 439_pInt, 443_pInt, 449_pInt, 457_pInt, 461_pInt, 463_pInt, & - 467_pInt, 479_pInt, 487_pInt, 491_pInt, 499_pInt, 503_pInt, 509_pInt, 521_pInt, 523_pInt, 541_pInt/) - - npvec(101:200) = (/ & + 467_pInt, 479_pInt, 487_pInt, 491_pInt, 499_pInt, 503_pInt, 509_pInt, 521_pInt, 523_pInt, 541_pInt, & + ! 101:200 547_pInt, 557_pInt, 563_pInt, 569_pInt, 571_pInt, 577_pInt, 587_pInt, 593_pInt, 599_pInt, 601_pInt, & 607_pInt, 613_pInt, 617_pInt, 619_pInt, 631_pInt, 641_pInt, 643_pInt, 647_pInt, 653_pInt, 659_pInt, & 661_pInt, 673_pInt, 677_pInt, 683_pInt, 691_pInt, 701_pInt, 709_pInt, 719_pInt, 727_pInt, 733_pInt, & @@ -2689,9 +2685,8 @@ function prime(n) 947_pInt, 953_pInt, 967_pInt, 971_pInt, 977_pInt, 983_pInt, 991_pInt, 997_pInt, 1009_pInt, 1013_pInt, & 1019_pInt, 1021_pInt, 1031_pInt, 1033_pInt, 1039_pInt, 1049_pInt, 1051_pInt, 1061_pInt, 1063_pInt, 1069_pInt, & 1087_pInt, 1091_pInt, 1093_pInt, 1097_pInt, 1103_pInt, 1109_pInt, 1117_pInt, 1123_pInt, 1129_pInt, 1151_pInt, & - 1153_pInt, 1163_pInt, 1171_pInt, 1181_pInt, 1187_pInt, 1193_pInt, 1201_pInt, 1213_pInt, 1217_pInt, 1223_pInt/) - - npvec(201:300) = (/ & + 1153_pInt, 1163_pInt, 1171_pInt, 1181_pInt, 1187_pInt, 1193_pInt, 1201_pInt, 1213_pInt, 1217_pInt, 1223_pInt, & + ! 201:300 1229_pInt, 1231_pInt, 1237_pInt, 1249_pInt, 1259_pInt, 1277_pInt, 1279_pInt, 1283_pInt, 1289_pInt, 1291_pInt, & 1297_pInt, 1301_pInt, 1303_pInt, 1307_pInt, 1319_pInt, 1321_pInt, 1327_pInt, 1361_pInt, 1367_pInt, 1373_pInt, & 1381_pInt, 1399_pInt, 1409_pInt, 1423_pInt, 1427_pInt, 1429_pInt, 1433_pInt, 1439_pInt, 1447_pInt, 1451_pInt, & @@ -2701,9 +2696,8 @@ function prime(n) 1663_pInt, 1667_pInt, 1669_pInt, 1693_pInt, 1697_pInt, 1699_pInt, 1709_pInt, 1721_pInt, 1723_pInt, 1733_pInt, & 1741_pInt, 1747_pInt, 1753_pInt, 1759_pInt, 1777_pInt, 1783_pInt, 1787_pInt, 1789_pInt, 1801_pInt, 1811_pInt, & 1823_pInt, 1831_pInt, 1847_pInt, 1861_pInt, 1867_pInt, 1871_pInt, 1873_pInt, 1877_pInt, 1879_pInt, 1889_pInt, & - 1901_pInt, 1907_pInt, 1913_pInt, 1931_pInt, 1933_pInt, 1949_pInt, 1951_pInt, 1973_pInt, 1979_pInt, 1987_pInt/) - - npvec(301:400) = (/ & + 1901_pInt, 1907_pInt, 1913_pInt, 1931_pInt, 1933_pInt, 1949_pInt, 1951_pInt, 1973_pInt, 1979_pInt, 1987_pInt, & + ! 301:400 1993_pInt, 1997_pInt, 1999_pInt, 2003_pInt, 2011_pInt, 2017_pInt, 2027_pInt, 2029_pInt, 2039_pInt, 2053_pInt, & 2063_pInt, 2069_pInt, 2081_pInt, 2083_pInt, 2087_pInt, 2089_pInt, 2099_pInt, 2111_pInt, 2113_pInt, 2129_pInt, & 2131_pInt, 2137_pInt, 2141_pInt, 2143_pInt, 2153_pInt, 2161_pInt, 2179_pInt, 2203_pInt, 2207_pInt, 2213_pInt, & @@ -2713,9 +2707,8 @@ function prime(n) 2437_pInt, 2441_pInt, 2447_pInt, 2459_pInt, 2467_pInt, 2473_pInt, 2477_pInt, 2503_pInt, 2521_pInt, 2531_pInt, & 2539_pInt, 2543_pInt, 2549_pInt, 2551_pInt, 2557_pInt, 2579_pInt, 2591_pInt, 2593_pInt, 2609_pInt, 2617_pInt, & 2621_pInt, 2633_pInt, 2647_pInt, 2657_pInt, 2659_pInt, 2663_pInt, 2671_pInt, 2677_pInt, 2683_pInt, 2687_pInt, & - 2689_pInt, 2693_pInt, 2699_pInt, 2707_pInt, 2711_pInt, 2713_pInt, 2719_pInt, 2729_pInt, 2731_pInt, 2741_pInt/) - - npvec(401:500) = (/ & + 2689_pInt, 2693_pInt, 2699_pInt, 2707_pInt, 2711_pInt, 2713_pInt, 2719_pInt, 2729_pInt, 2731_pInt, 2741_pInt, & + ! 401:500 2749_pInt, 2753_pInt, 2767_pInt, 2777_pInt, 2789_pInt, 2791_pInt, 2797_pInt, 2801_pInt, 2803_pInt, 2819_pInt, & 2833_pInt, 2837_pInt, 2843_pInt, 2851_pInt, 2857_pInt, 2861_pInt, 2879_pInt, 2887_pInt, 2897_pInt, 2903_pInt, & 2909_pInt, 2917_pInt, 2927_pInt, 2939_pInt, 2953_pInt, 2957_pInt, 2963_pInt, 2969_pInt, 2971_pInt, 2999_pInt, & @@ -2725,9 +2718,8 @@ function prime(n) 3259_pInt, 3271_pInt, 3299_pInt, 3301_pInt, 3307_pInt, 3313_pInt, 3319_pInt, 3323_pInt, 3329_pInt, 3331_pInt, & 3343_pInt, 3347_pInt, 3359_pInt, 3361_pInt, 3371_pInt, 3373_pInt, 3389_pInt, 3391_pInt, 3407_pInt, 3413_pInt, & 3433_pInt, 3449_pInt, 3457_pInt, 3461_pInt, 3463_pInt, 3467_pInt, 3469_pInt, 3491_pInt, 3499_pInt, 3511_pInt, & - 3517_pInt, 3527_pInt, 3529_pInt, 3533_pInt, 3539_pInt, 3541_pInt, 3547_pInt, 3557_pInt, 3559_pInt, 3571_pInt/) - - npvec(501:600) = (/ & + 3517_pInt, 3527_pInt, 3529_pInt, 3533_pInt, 3539_pInt, 3541_pInt, 3547_pInt, 3557_pInt, 3559_pInt, 3571_pInt, & + ! 501:600 3581_pInt, 3583_pInt, 3593_pInt, 3607_pInt, 3613_pInt, 3617_pInt, 3623_pInt, 3631_pInt, 3637_pInt, 3643_pInt, & 3659_pInt, 3671_pInt, 3673_pInt, 3677_pInt, 3691_pInt, 3697_pInt, 3701_pInt, 3709_pInt, 3719_pInt, 3727_pInt, & 3733_pInt, 3739_pInt, 3761_pInt, 3767_pInt, 3769_pInt, 3779_pInt, 3793_pInt, 3797_pInt, 3803_pInt, 3821_pInt, & @@ -2737,9 +2729,8 @@ function prime(n) 4073_pInt, 4079_pInt, 4091_pInt, 4093_pInt, 4099_pInt, 4111_pInt, 4127_pInt, 4129_pInt, 4133_pInt, 4139_pInt, & 4153_pInt, 4157_pInt, 4159_pInt, 4177_pInt, 4201_pInt, 4211_pInt, 4217_pInt, 4219_pInt, 4229_pInt, 4231_pInt, & 4241_pInt, 4243_pInt, 4253_pInt, 4259_pInt, 4261_pInt, 4271_pInt, 4273_pInt, 4283_pInt, 4289_pInt, 4297_pInt, & - 4327_pInt, 4337_pInt, 4339_pInt, 4349_pInt, 4357_pInt, 4363_pInt, 4373_pInt, 4391_pInt, 4397_pInt, 4409_pInt/) - - npvec(601:700) = (/ & + 4327_pInt, 4337_pInt, 4339_pInt, 4349_pInt, 4357_pInt, 4363_pInt, 4373_pInt, 4391_pInt, 4397_pInt, 4409_pInt, & + ! 601:700 4421_pInt, 4423_pInt, 4441_pInt, 4447_pInt, 4451_pInt, 4457_pInt, 4463_pInt, 4481_pInt, 4483_pInt, 4493_pInt, & 4507_pInt, 4513_pInt, 4517_pInt, 4519_pInt, 4523_pInt, 4547_pInt, 4549_pInt, 4561_pInt, 4567_pInt, 4583_pInt, & 4591_pInt, 4597_pInt, 4603_pInt, 4621_pInt, 4637_pInt, 4639_pInt, 4643_pInt, 4649_pInt, 4651_pInt, 4657_pInt, & @@ -2749,9 +2740,8 @@ function prime(n) 4943_pInt, 4951_pInt, 4957_pInt, 4967_pInt, 4969_pInt, 4973_pInt, 4987_pInt, 4993_pInt, 4999_pInt, 5003_pInt, & 5009_pInt, 5011_pInt, 5021_pInt, 5023_pInt, 5039_pInt, 5051_pInt, 5059_pInt, 5077_pInt, 5081_pInt, 5087_pInt, & 5099_pInt, 5101_pInt, 5107_pInt, 5113_pInt, 5119_pInt, 5147_pInt, 5153_pInt, 5167_pInt, 5171_pInt, 5179_pInt, & - 5189_pInt, 5197_pInt, 5209_pInt, 5227_pInt, 5231_pInt, 5233_pInt, 5237_pInt, 5261_pInt, 5273_pInt, 5279_pInt/) - - npvec(701:800) = (/ & + 5189_pInt, 5197_pInt, 5209_pInt, 5227_pInt, 5231_pInt, 5233_pInt, 5237_pInt, 5261_pInt, 5273_pInt, 5279_pInt, & + ! 701:800 5281_pInt, 5297_pInt, 5303_pInt, 5309_pInt, 5323_pInt, 5333_pInt, 5347_pInt, 5351_pInt, 5381_pInt, 5387_pInt, & 5393_pInt, 5399_pInt, 5407_pInt, 5413_pInt, 5417_pInt, 5419_pInt, 5431_pInt, 5437_pInt, 5441_pInt, 5443_pInt, & 5449_pInt, 5471_pInt, 5477_pInt, 5479_pInt, 5483_pInt, 5501_pInt, 5503_pInt, 5507_pInt, 5519_pInt, 5521_pInt, & @@ -2761,9 +2751,8 @@ function prime(n) 5801_pInt, 5807_pInt, 5813_pInt, 5821_pInt, 5827_pInt, 5839_pInt, 5843_pInt, 5849_pInt, 5851_pInt, 5857_pInt, & 5861_pInt, 5867_pInt, 5869_pInt, 5879_pInt, 5881_pInt, 5897_pInt, 5903_pInt, 5923_pInt, 5927_pInt, 5939_pInt, & 5953_pInt, 5981_pInt, 5987_pInt, 6007_pInt, 6011_pInt, 6029_pInt, 6037_pInt, 6043_pInt, 6047_pInt, 6053_pInt, & - 6067_pInt, 6073_pInt, 6079_pInt, 6089_pInt, 6091_pInt, 6101_pInt, 6113_pInt, 6121_pInt, 6131_pInt, 6133_pInt/) - - npvec(801:900) = (/ & + 6067_pInt, 6073_pInt, 6079_pInt, 6089_pInt, 6091_pInt, 6101_pInt, 6113_pInt, 6121_pInt, 6131_pInt, 6133_pInt, & + ! 801:900 6143_pInt, 6151_pInt, 6163_pInt, 6173_pInt, 6197_pInt, 6199_pInt, 6203_pInt, 6211_pInt, 6217_pInt, 6221_pInt, & 6229_pInt, 6247_pInt, 6257_pInt, 6263_pInt, 6269_pInt, 6271_pInt, 6277_pInt, 6287_pInt, 6299_pInt, 6301_pInt, & 6311_pInt, 6317_pInt, 6323_pInt, 6329_pInt, 6337_pInt, 6343_pInt, 6353_pInt, 6359_pInt, 6361_pInt, 6367_pInt, & @@ -2773,9 +2762,8 @@ function prime(n) 6679_pInt, 6689_pInt, 6691_pInt, 6701_pInt, 6703_pInt, 6709_pInt, 6719_pInt, 6733_pInt, 6737_pInt, 6761_pInt, & 6763_pInt, 6779_pInt, 6781_pInt, 6791_pInt, 6793_pInt, 6803_pInt, 6823_pInt, 6827_pInt, 6829_pInt, 6833_pInt, & 6841_pInt, 6857_pInt, 6863_pInt, 6869_pInt, 6871_pInt, 6883_pInt, 6899_pInt, 6907_pInt, 6911_pInt, 6917_pInt, & - 6947_pInt, 6949_pInt, 6959_pInt, 6961_pInt, 6967_pInt, 6971_pInt, 6977_pInt, 6983_pInt, 6991_pInt, 6997_pInt/) - - npvec(901:1000) = (/ & + 6947_pInt, 6949_pInt, 6959_pInt, 6961_pInt, 6967_pInt, 6971_pInt, 6977_pInt, 6983_pInt, 6991_pInt, 6997_pInt, & + ! 901:1000 7001_pInt, 7013_pInt, 7019_pInt, 7027_pInt, 7039_pInt, 7043_pInt, 7057_pInt, 7069_pInt, 7079_pInt, 7103_pInt, & 7109_pInt, 7121_pInt, 7127_pInt, 7129_pInt, 7151_pInt, 7159_pInt, 7177_pInt, 7187_pInt, 7193_pInt, 7207_pInt, & 7211_pInt, 7213_pInt, 7219_pInt, 7229_pInt, 7237_pInt, 7243_pInt, 7247_pInt, 7253_pInt, 7283_pInt, 7297_pInt, & @@ -2785,9 +2773,8 @@ function prime(n) 7573_pInt, 7577_pInt, 7583_pInt, 7589_pInt, 7591_pInt, 7603_pInt, 7607_pInt, 7621_pInt, 7639_pInt, 7643_pInt, & 7649_pInt, 7669_pInt, 7673_pInt, 7681_pInt, 7687_pInt, 7691_pInt, 7699_pInt, 7703_pInt, 7717_pInt, 7723_pInt, & 7727_pInt, 7741_pInt, 7753_pInt, 7757_pInt, 7759_pInt, 7789_pInt, 7793_pInt, 7817_pInt, 7823_pInt, 7829_pInt, & - 7841_pInt, 7853_pInt, 7867_pInt, 7873_pInt, 7877_pInt, 7879_pInt, 7883_pInt, 7901_pInt, 7907_pInt, 7919_pInt/) - - npvec(1001:1100) = (/ & + 7841_pInt, 7853_pInt, 7867_pInt, 7873_pInt, 7877_pInt, 7879_pInt, 7883_pInt, 7901_pInt, 7907_pInt, 7919_pInt, & + ! 1001:1100 7927_pInt, 7933_pInt, 7937_pInt, 7949_pInt, 7951_pInt, 7963_pInt, 7993_pInt, 8009_pInt, 8011_pInt, 8017_pInt, & 8039_pInt, 8053_pInt, 8059_pInt, 8069_pInt, 8081_pInt, 8087_pInt, 8089_pInt, 8093_pInt, 8101_pInt, 8111_pInt, & 8117_pInt, 8123_pInt, 8147_pInt, 8161_pInt, 8167_pInt, 8171_pInt, 8179_pInt, 8191_pInt, 8209_pInt, 8219_pInt, & @@ -2797,9 +2784,8 @@ function prime(n) 8513_pInt, 8521_pInt, 8527_pInt, 8537_pInt, 8539_pInt, 8543_pInt, 8563_pInt, 8573_pInt, 8581_pInt, 8597_pInt, & 8599_pInt, 8609_pInt, 8623_pInt, 8627_pInt, 8629_pInt, 8641_pInt, 8647_pInt, 8663_pInt, 8669_pInt, 8677_pInt, & 8681_pInt, 8689_pInt, 8693_pInt, 8699_pInt, 8707_pInt, 8713_pInt, 8719_pInt, 8731_pInt, 8737_pInt, 8741_pInt, & - 8747_pInt, 8753_pInt, 8761_pInt, 8779_pInt, 8783_pInt, 8803_pInt, 8807_pInt, 8819_pInt, 8821_pInt, 8831_pInt/) - - npvec(1101:1200) = (/ & + 8747_pInt, 8753_pInt, 8761_pInt, 8779_pInt, 8783_pInt, 8803_pInt, 8807_pInt, 8819_pInt, 8821_pInt, 8831_pInt, & + ! 1101:1200 8837_pInt, 8839_pInt, 8849_pInt, 8861_pInt, 8863_pInt, 8867_pInt, 8887_pInt, 8893_pInt, 8923_pInt, 8929_pInt, & 8933_pInt, 8941_pInt, 8951_pInt, 8963_pInt, 8969_pInt, 8971_pInt, 8999_pInt, 9001_pInt, 9007_pInt, 9011_pInt, & 9013_pInt, 9029_pInt, 9041_pInt, 9043_pInt, 9049_pInt, 9059_pInt, 9067_pInt, 9091_pInt, 9103_pInt, 9109_pInt, & @@ -2809,9 +2795,8 @@ function prime(n) 9391_pInt, 9397_pInt, 9403_pInt, 9413_pInt, 9419_pInt, 9421_pInt, 9431_pInt, 9433_pInt, 9437_pInt, 9439_pInt, & 9461_pInt, 9463_pInt, 9467_pInt, 9473_pInt, 9479_pInt, 9491_pInt, 9497_pInt, 9511_pInt, 9521_pInt, 9533_pInt, & 9539_pInt, 9547_pInt, 9551_pInt, 9587_pInt, 9601_pInt, 9613_pInt, 9619_pInt, 9623_pInt, 9629_pInt, 9631_pInt, & - 9643_pInt, 9649_pInt, 9661_pInt, 9677_pInt, 9679_pInt, 9689_pInt, 9697_pInt, 9719_pInt, 9721_pInt, 9733_pInt/) - - npvec(1201:1300) = (/ & + 9643_pInt, 9649_pInt, 9661_pInt, 9677_pInt, 9679_pInt, 9689_pInt, 9697_pInt, 9719_pInt, 9721_pInt, 9733_pInt, & + ! 1201:1300 9739_pInt, 9743_pInt, 9749_pInt, 9767_pInt, 9769_pInt, 9781_pInt, 9787_pInt, 9791_pInt, 9803_pInt, 9811_pInt, & 9817_pInt, 9829_pInt, 9833_pInt, 9839_pInt, 9851_pInt, 9857_pInt, 9859_pInt, 9871_pInt, 9883_pInt, 9887_pInt, & 9901_pInt, 9907_pInt, 9923_pInt, 9929_pInt, 9931_pInt, 9941_pInt, 9949_pInt, 9967_pInt, 9973_pInt,10007_pInt, & @@ -2821,9 +2806,8 @@ function prime(n) 10273_pInt,10289_pInt,10301_pInt,10303_pInt,10313_pInt,10321_pInt,10331_pInt,10333_pInt,10337_pInt,10343_pInt, & 10357_pInt,10369_pInt,10391_pInt,10399_pInt,10427_pInt,10429_pInt,10433_pInt,10453_pInt,10457_pInt,10459_pInt, & 10463_pInt,10477_pInt,10487_pInt,10499_pInt,10501_pInt,10513_pInt,10529_pInt,10531_pInt,10559_pInt,10567_pInt, & - 10589_pInt,10597_pInt,10601_pInt,10607_pInt,10613_pInt,10627_pInt,10631_pInt,10639_pInt,10651_pInt,10657_pInt/) - - npvec(1301:1400) = (/ & + 10589_pInt,10597_pInt,10601_pInt,10607_pInt,10613_pInt,10627_pInt,10631_pInt,10639_pInt,10651_pInt,10657_pInt, & + ! 1301:1400 10663_pInt,10667_pInt,10687_pInt,10691_pInt,10709_pInt,10711_pInt,10723_pInt,10729_pInt,10733_pInt,10739_pInt, & 10753_pInt,10771_pInt,10781_pInt,10789_pInt,10799_pInt,10831_pInt,10837_pInt,10847_pInt,10853_pInt,10859_pInt, & 10861_pInt,10867_pInt,10883_pInt,10889_pInt,10891_pInt,10903_pInt,10909_pInt,19037_pInt,10939_pInt,10949_pInt, & @@ -2833,9 +2817,8 @@ function prime(n) 11257_pInt,11261_pInt,11273_pInt,11279_pInt,11287_pInt,11299_pInt,11311_pInt,11317_pInt,11321_pInt,11329_pInt, & 11351_pInt,11353_pInt,11369_pInt,11383_pInt,11393_pInt,11399_pInt,11411_pInt,11423_pInt,11437_pInt,11443_pInt, & 11447_pInt,11467_pInt,11471_pInt,11483_pInt,11489_pInt,11491_pInt,11497_pInt,11503_pInt,11519_pInt,11527_pInt, & - 11549_pInt,11551_pInt,11579_pInt,11587_pInt,11593_pInt,11597_pInt,11617_pInt,11621_pInt,11633_pInt,11657_pInt/) - - npvec(1401:1500) = (/ & + 11549_pInt,11551_pInt,11579_pInt,11587_pInt,11593_pInt,11597_pInt,11617_pInt,11621_pInt,11633_pInt,11657_pInt, & + ! 1401:1500 11677_pInt,11681_pInt,11689_pInt,11699_pInt,11701_pInt,11717_pInt,11719_pInt,11731_pInt,11743_pInt,11777_pInt, & 11779_pInt,11783_pInt,11789_pInt,11801_pInt,11807_pInt,11813_pInt,11821_pInt,11827_pInt,11831_pInt,11833_pInt, & 11839_pInt,11863_pInt,11867_pInt,11887_pInt,11897_pInt,11903_pInt,11909_pInt,11923_pInt,11927_pInt,11933_pInt, & @@ -2845,8 +2828,7 @@ function prime(n) 12227_pInt,12239_pInt,12241_pInt,12251_pInt,12253_pInt,12263_pInt,12269_pInt,12277_pInt,12281_pInt,12289_pInt, & 12301_pInt,12323_pInt,12329_pInt,12343_pInt,12347_pInt,12373_pInt,12377_pInt,12379_pInt,12391_pInt,12401_pInt, & 12409_pInt,12413_pInt,12421_pInt,12433_pInt,12437_pInt,12451_pInt,12457_pInt,12473_pInt,12479_pInt,12487_pInt, & - 12491_pInt,12497_pInt,12503_pInt,12511_pInt,12517_pInt,12527_pInt,12539_pInt,12541_pInt,12547_pInt,12553_pInt/) - + 12491_pInt,12497_pInt,12503_pInt,12511_pInt,12517_pInt,12527_pInt,12539_pInt,12541_pInt,12547_pInt,12553_pInt] endif if(n == -1_pInt) then @@ -2861,10 +2843,10 @@ function prime(n) end function prime -!************************************************************************** -! volume of tetrahedron given by four vertices -!************************************************************************** -pure function math_volTetrahedron(v1,v2,v3,v4) +!-------------------------------------------------------------------------------------------------- +!> @brief volume of tetrahedron given by four vertices +!-------------------------------------------------------------------------------------------------- +pure function math_volTetrahedron(v1,v2,v3,v4) implicit none @@ -2876,47 +2858,46 @@ pure function math_volTetrahedron(v1,v2,v3,v4) m(1:3,2) = v2-v3 m(1:3,3) = v3-v4 - math_volTetrahedron = math_det33(m)/6.0_pReal + math_volTetrahedron = math_det33(m)/6.0_pReal end function math_volTetrahedron -!************************************************************************** -! rotate 33 tensor forward -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief rotate 33 tensor forward +!-------------------------------------------------------------------------------------------------- pure function math_rotate_forward33(tensor,rot_tensor) implicit none real(pReal), dimension(3,3) :: math_rotate_forward33 real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - + math_rotate_forward33 = math_mul33x33(rot_tensor,& math_mul33x33(tensor,math_transpose33(rot_tensor))) - + end function math_rotate_forward33 -!************************************************************************** -! rotate 33 tensor backward -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief rotate 33 tensor backward +!-------------------------------------------------------------------------------------------------- pure function math_rotate_backward33(tensor,rot_tensor) implicit none real(pReal), dimension(3,3) :: math_rotate_backward33 real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor - + math_rotate_backward33 = math_mul33x33(math_transpose33(rot_tensor),& math_mul33x33(tensor,rot_tensor)) - + end function math_rotate_backward33 -!************************************************************************** -! rotate 3333 tensor -! C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop -!************************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop +!-------------------------------------------------------------------------------------------------- pure function math_rotate_forward3333(tensor,rot_tensor) implicit none @@ -2925,7 +2906,7 @@ pure function math_rotate_forward3333(tensor,rot_tensor) real(pReal), dimension(3,3), intent(in) :: rot_tensor real(pReal), dimension(3,3,3,3), intent(in) :: tensor integer(pInt) :: i,j,k,l,m,n,o,p - + math_rotate_forward3333= 0.0_pReal do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt @@ -2933,7 +2914,7 @@ pure function math_rotate_forward3333(tensor,rot_tensor) math_rotate_forward3333(i,j,k,l) = tensor(i,j,k,l)+rot_tensor(m,i)*rot_tensor(n,j)*& rot_tensor(o,k)*rot_tensor(p,l)*tensor(m,n,o,p) enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo - + end function math_rotate_forward3333 @@ -3473,7 +3454,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) field_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order 1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions field_fourier,(/res(3),res(2) ,res1_red/),& - 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) + 1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag) fftw_back = fftw_plan_many_dft_c2r(3_pInt,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& curl_fourier,(/res(3),res(2) ,res1_red/),& @@ -3487,7 +3468,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) enddo; enddo; enddo call fftw_execute_dft_r2c(fftw_forth, field_real, field_fourier) - + !remove highest frequency in each direction if(res(1)>1_pInt) & field_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,& @@ -3498,18 +3479,18 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl) if(res(3)>1_pInt) & field_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,& 1:vec_tens,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) - + do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) k_s(3) = k - 1_pInt if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) do j = 1_pInt, res(2) k_s(2) = j - 1_pInt - if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) + if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) do i = 1_pInt, res1_red k_s(1) = i - 1_pInt xi(i,j,k,1:3) = real(k_s, pReal)/geomdim enddo; enddo; enddo - + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do l = 1_pInt, vec_tens curl_fourier(i,j,k,l,1) = ( field_fourier(i,j,k,l,3)*xi(i,j,k,2)& @@ -3601,23 +3582,23 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) divergence_fourier,(/res(3),res(2) ,res1_red/),& 1_pInt, res(3)*res(2)* res1_red,& divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),& - 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding + 1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) field_real(i,j,k,1:vec_tens,1:3) = field(i,j,k,1:vec_tens,1:3) ! ensure that data is aligned properly (fftw_alloc) enddo; enddo; enddo - + call fftw_execute_dft_r2c(fftw_forth, field_real, field_fourier) do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) k_s(3) = k - 1_pInt if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) do j = 1_pInt, res(2) k_s(2) = j - 1_pInt - if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) + if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) do i = 1_pInt, res1_red k_s(1) = i - 1_pInt xi(i,j,k,1:3) = real(k_s, pReal)/geomdim enddo; enddo; enddo - + !remove highest frequency in each direction if(res(1)>1_pInt) & field_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,& @@ -3628,7 +3609,7 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) if(res(3)>1_pInt) & field_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,& 1:vec_tens,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) - + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do l = 1_pInt, vec_tens divergence_fourier(i,j,k,l)=sum(field_fourier(i,j,k,l,1:3)*cmplx(xi(i,j,k,1:3),0.0_pReal,pReal))&