From 715e7fd91870f1bffa0b6f9a7949b80705be11b8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 15 May 2014 08:52:16 +0000 Subject: [PATCH] some small iprovements: reading in geometry in mesh, readability of math, output formatting in fesolving and debug, hdf5 in constitutive --- code/FEsolving.f90 | 2 +- code/constitutive_j2.f90 | 2 +- code/debug.f90 | 2 +- code/math.f90 | 22 +++++++++++----------- code/mesh.f90 | 10 ++++++---- code/prec.f90 | 4 ++-- 6 files changed, 22 insertions(+), 20 deletions(-) diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index dff449057..6334d78f7 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -87,7 +87,7 @@ subroutine FE_init write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' write(6,'(a)') ' $Id$' - write(6,'(a16,a)') ' Current time : ',IO_timeStamp() + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" modelName = getSolverJobName() diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 225ead1f7..6f772c719 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -327,7 +327,7 @@ subroutine constitutive_j2_init(fileUnit) if (any(numerics_integrator == 4_pInt)) & allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal) if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (sizeDotState,NofMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) #endif endif enddo initializeInstances diff --git a/code/debug.f90 b/code/debug.f90 index 96a8866f5..b1b74b5c2 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -131,7 +131,7 @@ subroutine debug_init write(6,'(/,a)') ' <<<+- debug init -+>>>' write(6,'(a)') ' $Id$' - write(6,'(a16,a)') ' Current time : ',IO_timeStamp() + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" if (allocated(debug_StressLoopDistribution)) & diff --git a/code/math.f90 b/code/math.f90 index 155d3b043..5d2fc185d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -1308,9 +1308,9 @@ pure function math_qRot(Q,v) 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 = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), & + v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), & + 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 @@ -1337,8 +1337,8 @@ pure function math_RtoEuler(R) ! 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 + myVal = min(myVal, 1.0_pReal) + myVal = max(myVal,-1.0_pReal) math_RtoEuler(2) = acos(myVal) @@ -1347,23 +1347,23 @@ pure function math_RtoEuler(R) math_RtoEuler(3) = 0.0_pReal ! calculate phi1 myVal=R(1,1)/squvw - if(myVal > 1.0_pReal) myVal = 1.0_pReal - if(myVal < -1.0_pReal) myVal = -1.0_pReal + myVal = min(myVal, 1.0_pReal) + myVal = max(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 ! calculate phi2 myVal=R(2,3)/sqhk - if(myVal > 1.0_pReal) myVal = 1.0_pReal - if(myVal < -1.0_pReal) myVal = -1.0_pReal + myVal = min(myVal, 1.0_pReal) + myVal = max(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 + myVal = min(myVal, 1.0_pReal) + myVal = max(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) diff --git a/code/mesh.f90 b/code/mesh.f90 index b59ffba75..62448081c 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -1251,7 +1251,8 @@ subroutine mesh_spectral_build_elements(fileUnit) e, i, & headerLength = 0_pInt, & maxIntCount, & - homog + homog, & + elemType integer(pInt), dimension(:), allocatable :: & microstructures integer(pInt), dimension(3) :: & @@ -1295,16 +1296,17 @@ subroutine mesh_spectral_build_elements(fileUnit) read(fileUnit,'(a65536)') line enddo - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt - allocate (microstructures (1_pInt+maxIntCount)); microstructures = 2_pInt + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems), source = 0_pInt) + allocate (microstructures (1_pInt+maxIntCount), source = 1_pInt) + elemType = FE_mapElemtype('C3D8R') e = 0_pInt do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) microstructures = IO_continuousIntValues(fileUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements do i = 1_pInt,microstructures(1_pInt) e = e+1_pInt ! valid element entry mesh_element( 1,e) = e ! FE id - mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! elem type + mesh_element( 2,e) = elemType ! elem type mesh_element( 3,e) = homog ! homogenization mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & diff --git a/code/prec.f90 b/code/prec.f90 index e68dea7de..04d53396b 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -71,8 +71,8 @@ module prec previousDotState2, & dotState_backup, & RK4dotState, & - aTolState, & - RKCK45dotState + aTolState + real(pReal), pointer, dimension(:,:,:) :: RKCK45dotState end type #endif