From fa39a7423b9d55fc69f51fb363e7b593db063628 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 Jan 2020 17:15:02 +0100 Subject: [PATCH 01/12] preparing for actual use --- python/damask/grid_filters.py | 18 +++++++++++++++--- python/tests/test_grid_filters.py | 6 ++++++ 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 36ce3b8e2..3850bf70b 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -338,7 +338,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): Parameters ---------- coord0 : numpy.ndarray - array of undeformed nodal coordinates + array of undeformed nodal coordinates. ordered : bool, optional expect coord0 data to be ordered (x fast, z slow). @@ -365,7 +365,19 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): def regrid(size,F,new_grid): - """tbd.""" + """ + Return mapping from coordinates in deformed configuration to a regular grid + + Parameters + ---------- + size : numpy.ndarray + physical size + F : numpy.ndarray + deformation gradient field + new_grid : numpy.ndarray + new grid for undeformed coordinates + + """ c = cell_coord0(F.shape[:3][::-1],size) \ + cell_displacement_avg(size,F) \ + cell_displacement_fluct(size,F) @@ -376,4 +388,4 @@ def regrid(size,F,new_grid): c[np.where(c[:,:,:,d]>outer[d])] -= outer[d] tree = spatial.cKDTree(c.reshape((-1,3)),boxsize=outer) - return tree.query(cell_coord0(new_grid,outer))[1] + return tree.query(cell_coord0(new_grid,outer))[1].flatten() diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index a5455e1ae..ed1336842 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -77,3 +77,9 @@ class TestGridFilters: grid = np.random.randint(8,32,(3)) F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) # noqa assert np.allclose(eval('grid_filters.{}_displacement_fluct(size,F)'.format(mode)),0.0) + + def test_regrid(self): + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + F = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3)) + assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod())) From b938f1a98df60515de73800d1cebd2a2f9f1369d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 14:01:14 +0100 Subject: [PATCH 02/12] polishing * constants in CAPITALS * more tests * 'forall' is deprecated in Fortran 2018 --- src/math.f90 | 181 +++++++++++++++++++++++++++++---------------------- 1 file changed, 104 insertions(+), 77 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 4be2d4fc9..c63f30881 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -8,14 +8,16 @@ module math use prec use IO - use debug use numerics + implicit none public #if __INTEL_COMPILER >= 1900 ! do not make use associated entities available to other modules private :: & - prec + prec, & + IO, & + numerics #endif real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter @@ -24,24 +26,24 @@ module math complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) real(pReal), dimension(3,3), parameter :: & - MATH_I3 = reshape([& + math_I3 = reshape([& 1.0_pReal,0.0_pReal,0.0_pReal, & 0.0_pReal,1.0_pReal,0.0_pReal, & 0.0_pReal,0.0_pReal,1.0_pReal & ],[3,3]) !< 3x3 Identity real(pReal), dimension(6), parameter, private :: & - nrmMandel = [& + NRMMANDEL = [& 1.0_pReal, 1.0_pReal, 1.0_pReal, & sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward) real(pReal), dimension(6), parameter , private :: & - invnrmMandel = [& + INVNRMMANDEL = [& 1.0_pReal, 1.0_pReal, 1.0_pReal, & 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward) integer, dimension (2,6), parameter, private :: & - mapNye = reshape([& + MAPNYE = reshape([& 1,1, & 2,2, & 3,3, & @@ -51,7 +53,7 @@ module math ],[2,6]) !< arrangement in Nye notation. integer, dimension (2,6), parameter, private :: & - mapVoigt = reshape([& + MAPVOIGT = reshape([& 1,1, & 2,2, & 3,3, & @@ -61,7 +63,7 @@ module math ],[2,6]) !< arrangement in Voigt notation integer, dimension (2,9), parameter, private :: & - mapPlain = reshape([& + MAPPLAIN = reshape([& 1,1, & 1,2, & 1,3, & @@ -129,13 +131,13 @@ recursive subroutine math_sort(a, istart, iend, sortDim) integer, dimension(:,:), intent(inout) :: a integer, intent(in),optional :: istart,iend, sortDim integer :: ipivot,s,e,d - + if(present(istart)) then s = istart else s = lbound(a,2) endif - + if(present(iend)) then e = iend else @@ -147,7 +149,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim) else d = 1 endif - + if (s < e) then ipivot = qsort_partition(a,s, e, d) call math_sort(a, s, ipivot-1, d) @@ -232,14 +234,14 @@ end function math_range !-------------------------------------------------------------------------------------------------- !> @brief second rank identity tensor of specified dimension !-------------------------------------------------------------------------------------------------- -pure function math_identity2nd(dimen) +pure function math_identity2nd(d) - integer, intent(in) :: dimen !< tensor dimension + integer, intent(in) :: d !< tensor dimension integer :: i - real(pReal), dimension(dimen,dimen) :: math_identity2nd + real(pReal), dimension(d,d) :: math_identity2nd math_identity2nd = 0.0_pReal - do i=1, dimen + do i=1,d math_identity2nd(i,i) = 1.0_pReal enddo @@ -250,16 +252,17 @@ end function math_identity2nd !> @brief symmetric fourth rank identity tensor of specified dimension ! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself !-------------------------------------------------------------------------------------------------- -pure function math_identity4th(dimen) +pure function math_identity4th(d) - integer, intent(in) :: dimen !< tensor dimension + integer, intent(in) :: d !< tensor dimension integer :: i,j,k,l - real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th - real(pReal), dimension(dimen,dimen) :: identity2nd + real(pReal), dimension(d,d,d,d) :: math_identity4th + real(pReal), dimension(d,d) :: identity2nd - identity2nd = math_identity2nd(dimen) - forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) & + identity2nd = math_identity2nd(d) + do i=1,d; do j=1,d; do k=1,d; do l=1,d math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) + enddo; enddo; enddo; enddo end function math_identity4th @@ -324,7 +327,9 @@ pure function math_outer(A,B) real(pReal), dimension(size(A,1),size(B,1)) :: math_outer integer :: i,j - forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) + do i=1,size(A,1); do j=1,size(B,1) + math_outer(i,j) = A(i)*B(j) + enddo; enddo end function math_outer @@ -347,11 +352,13 @@ end function math_inner !-------------------------------------------------------------------------------------------------- real(pReal) pure function math_mul33xx33(A,B) - real(pReal), dimension(3,3), intent(in) :: A,B + real(pReal), dimension(3,3), intent(in) :: A,B integer :: i,j - real(pReal), dimension(3,3) :: C + real(pReal), dimension(3,3) :: C - forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j) + do i=1,3; do j=1,3 + C(i,j) = A(i,j) * B(i,j) + enddo; enddo math_mul33xx33 = sum(C) end function math_mul33xx33 @@ -362,12 +369,14 @@ end function math_mul33xx33 !-------------------------------------------------------------------------------------------------- pure function math_mul3333xx33(A,B) + real(pReal), dimension(3,3,3,3), intent(in) :: A + real(pReal), dimension(3,3), intent(in) :: B real(pReal), dimension(3,3) :: math_mul3333xx33 - real(pReal), dimension(3,3,3,3), intent(in) :: A - real(pReal), dimension(3,3), intent(in) :: B integer :: i,j - forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) + do i=1,3; do j=1,3 + math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) + enddo; enddo end function math_mul3333xx33 @@ -378,12 +387,13 @@ end function math_mul3333xx33 pure function math_mul3333xx3333(A,B) integer :: i,j,k,l - real(pReal), dimension(3,3,3,3), intent(in) :: A - real(pReal), dimension(3,3,3,3), intent(in) :: B + real(pReal), dimension(3,3,3,3), intent(in) :: A + real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 - forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) & + do i=1,3; do j=1,3; do k=1,3; do l=1,3 math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) + enddo; enddo; enddo; enddo end function math_mul3333xx3333 @@ -393,12 +403,11 @@ end function math_mul3333xx3333 !-------------------------------------------------------------------------------------------------- pure function math_exp33(A,n) - integer :: i - integer, intent(in), optional :: n real(pReal), dimension(3,3), intent(in) :: A + integer, intent(in), optional :: n real(pReal), dimension(3,3) :: B, math_exp33 real(pReal) :: invFac - integer :: n_ + integer :: n_,i if (present(n)) then n_ = n @@ -646,7 +655,7 @@ pure function math_33to9(m33) integer :: i do i = 1, 9 - math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) + math_33to9(i) = m33(MAPPLAIN(1,i),MAPPLAIN(2,i)) enddo end function math_33to9 @@ -663,7 +672,7 @@ pure function math_9to33(v9) integer :: i do i = 1, 9 - math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) + math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i) enddo end function math_9to33 @@ -685,13 +694,13 @@ pure function math_sym33to6(m33,weighted) integer :: i if(present(weighted)) then - w = merge(nrmMandel,1.0_pReal,weighted) + w = merge(NRMMANDEL,1.0_pReal,weighted) else - w = nrmMandel + w = NRMMANDEL endif do i = 1, 6 - math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) + math_sym33to6(i) = w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)) enddo end function math_sym33to6 @@ -713,14 +722,14 @@ pure function math_6toSym33(v6,weighted) integer :: i if(present(weighted)) then - w = merge(invnrmMandel,1.0_pReal,weighted) + w = merge(INVNRMMANDEL,1.0_pReal,weighted) else - w = invnrmMandel + w = INVNRMMANDEL endif do i=1,6 - math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i) - math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i) + math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i) + math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i) enddo end function math_6toSym33 @@ -737,7 +746,7 @@ pure function math_3333to99(m3333) integer :: i,j do i=1,9; do j=1,9 - math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) + math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) enddo; enddo end function math_3333to99 @@ -754,7 +763,7 @@ pure function math_99to3333(m99) integer :: i,j do i=1,9; do j=1,9 - math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j) + math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) enddo; enddo end function math_99to3333 @@ -776,13 +785,13 @@ pure function math_sym3333to66(m3333,weighted) integer :: i,j if(present(weighted)) then - w = merge(nrmMandel,1.0_pReal,weighted) + w = merge(NRMMANDEL,1.0_pReal,weighted) else - w = nrmMandel + w = NRMMANDEL endif do i=1,6; do j=1,6 - math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) + math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) enddo; enddo end function math_sym3333to66 @@ -804,16 +813,16 @@ pure function math_66toSym3333(m66,weighted) integer :: i,j if(present(weighted)) then - w = merge(invnrmMandel,1.0_pReal,weighted) + w = merge(INVNRMMANDEL,1.0_pReal,weighted) else - w = invnrmMandel + w = INVNRMMANDEL endif do i=1,6; do j=1,6 - math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) - math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j) + math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j) enddo; enddo end function math_66toSym3333 @@ -829,10 +838,10 @@ pure function math_Voigt66to3333(m66) integer :: i,j do i=1,6; do j=1, 6 - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) - math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) + math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) + math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) + math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) + math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) enddo; enddo end function math_Voigt66to3333 @@ -1201,7 +1210,7 @@ end function math_invariantsSym33 integer pure function math_factorial(n) integer, intent(in) :: n - + math_factorial = product(math_range(n)) end function math_factorial @@ -1233,7 +1242,7 @@ integer pure function math_multinomial(alpha) integer, intent(in), dimension(:) :: alpha integer :: i - + math_multinomial = 1 do i = 1, size(alpha) math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) @@ -1293,12 +1302,12 @@ end function math_clip !> @brief check correctness of (some) math functions !-------------------------------------------------------------------------------------------------- subroutine unitTest - + integer, dimension(2,4) :: & sort_in_ = reshape([+1,+5, +5,+6, -1,-1, +3,-2],[2,4]) integer, dimension(2,4), parameter :: & sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4]) - + integer, dimension(5) :: range_out_ = [1,2,3,4,5] real(pReal) :: det @@ -1308,8 +1317,12 @@ subroutine unitTest real(pReal), dimension(3,3) :: t33,t33_2 real(pReal), dimension(6,6) :: t66 real(pReal), dimension(9,9) :: t99,t99_2 + real(pReal), dimension(:,:), & + allocatable :: txx,txx_2 + real(pReal) :: r + integer :: d logical :: e - + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) & call IO_error(0,ext_msg='math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]') @@ -1322,15 +1335,19 @@ subroutine unitTest math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) & call IO_error(0,ext_msg='math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]') - call math_sort(sort_in_,1,3,2) if(any(sort_in_ /= sort_out_)) & call IO_error(0,ext_msg='math_sort') - + if(any(math_range(5) /= range_out_)) & call IO_error(0,ext_msg='math_range') - - + + if(any(dNeq(math_exp33(math_I3,0),math_I3))) & + call IO_error(0,ext_msg='math_exp33(math_I3,1)') + + if(any(dNeq(math_exp33(math_I3,256),exp(1.0_pReal)*math_I3))) & + call IO_error(0,ext_msg='math_exp33(math_I3,256)') + call random_number(v9) if(any(dNeq(math_33to9(math_9to33(v9)),v9))) & call IO_error(0,ext_msg='math_33to9/math_9to33') @@ -1338,56 +1355,66 @@ subroutine unitTest call random_number(t99) if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & call IO_error(0,ext_msg='math_3333to99/math_99to3333') - + call random_number(v6) if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & call IO_error(0,ext_msg='math_sym33to6/math_6toSym33') - + call random_number(t66) if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) & call IO_error(0,ext_msg='math_sym3333to66/math_66toSym3333') - + call random_number(v6) if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & call IO_error(0,ext_msg='math_symmetric33') - + call random_number(v3_1) call random_number(v3_2) call random_number(v3_3) call random_number(v3_4) - + if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & call IO_error(0,ext_msg='math_volTetrahedron') - + call random_number(t33) if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & call IO_error(0,ext_msg='math_det33/math_detSym33') + if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) & + call IO_error(0,ext_msg='math_inv33(math_I3)') + do while(abs(math_det33(t33))<1.0e-9_pReal) call random_number(t33) enddo if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_identity2nd(3),tol=1.0e-9_pReal))) & call IO_error(0,ext_msg='math_inv33') - + call math_invert33(t33_2,det,e,t33) if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) & call IO_error(0,ext_msg='math_invert33: T:T^-1 != I') if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & call IO_error(0,ext_msg='math_invert33 (determinant)') - + call math_invert(t33_2,e,t33) if(any(dNeq0(matmul(t33,t33_2) - math_identity2nd(3),tol=1.0e-9_pReal)) .or. e) & call IO_error(0,ext_msg='math_invert t33') - + t33_2 = transpose(math_rotationalPart33(t33)) if(any(dNeq0(math_rotationalPart33(matmul(t33_2,t33)) - MATH_I3,tol=5.0e-4_pReal))) & call IO_error(0,ext_msg='math_rotationalPart33') + call random_number(r) + d = int(r*5.0_pReal) + 1 + txx = math_identity2nd(d) + allocate(txx_2(d,d)) + call math_invert(txx_2,e,txx) + if(any(dNeq0(txx_2,txx) .or. e)) & + call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd') call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) & - call IO_error(0,ext_msg='math_invert t99') + call IO_error(0,ext_msg='math_invert(t99)') if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & call IO_error(0,ext_msg='math_clip') From 64be6a277d5f30b05d77da2fa093334c7794a462 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 14:48:15 +0100 Subject: [PATCH 03/12] it's the surface, not the edge --- src/mesh_grid.f90 | 84 +++++++++++++++++++++++------------------------ 1 file changed, 42 insertions(+), 42 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index b1ba6b590..6fe326fd0 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -18,10 +18,10 @@ module mesh_grid use discretization use geometry_plastic_nonlocal use FEsolving - + implicit none private - + integer, dimension(3), public, protected :: & grid !< (global) grid integer, public, protected :: & @@ -32,10 +32,10 @@ module mesh_grid real(pReal), public, protected :: & size3, & !< (local) size in 3rd direction size3offset !< (local) size offset in 3rd direction - + public :: & mesh_init - + contains @@ -45,7 +45,7 @@ contains subroutine mesh_init(ip,el) integer, intent(in), optional :: el, ip ! for compatibility reasons - + include 'fftw3-mpi.f03' real(pReal), dimension(3) :: & mySize, & !< domain size of this process @@ -93,7 +93,7 @@ subroutine mesh_init(ip,el) call discretization_init(homogenizationAt,microstructureAt, & IPcoordinates0(myGrid,mySize,grid3Offset), & Nodes0(myGrid,mySize,grid3Offset),& - merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer + merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer (grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1) worldrank<1)) @@ -113,8 +113,8 @@ subroutine mesh_init(ip,el) ! geometry information required by the nonlocal CP model call geometry_plastic_nonlocal_setIPvolume(reshape([(product(mySize/real(myGrid,pReal)),j=1,product(myGrid))], & [1,product(myGrid)])) - call geometry_plastic_nonlocal_setIParea (cellEdgeArea(mySize,myGrid)) - call geometry_plastic_nonlocal_setIPareaNormal (cellEdgeNormal(product(myGrid))) + call geometry_plastic_nonlocal_setIParea (cellSurfaceArea(mySize,myGrid)) + call geometry_plastic_nonlocal_setIPareaNormal (cellSurfaceNormal(product(myGrid))) call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(myGrid)) !-------------------------------------------------------------------------------------------------- @@ -127,7 +127,7 @@ end subroutine mesh_init !-------------------------------------------------------------------------------------------------- !> @brief Parses geometry file -!> @details important variables have an implicit "save" attribute. Therefore, this function is +!> @details important variables have an implicit "save" attribute. Therefore, this function is ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) @@ -140,7 +140,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) integer, dimension(:), intent(out), allocatable :: & microstructure, & homogenization - + character(len=:), allocatable :: rawData character(len=65536) :: line integer, allocatable, dimension(:) :: chunkPos @@ -154,9 +154,9 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) l, & !< line counter c, & !< counter for # microstructures in line o, & !< order of "to" packing - e, & !< "element", i.e. spectral collocation point + e, & !< "element", i.e. spectral collocation point i, j - + grid = -1 geomSize = -1.0_pReal @@ -169,7 +169,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) allocate(character(len=fileLength)::rawData) read(fileUnit) rawData close(fileUnit) - + !-------------------------------------------------------------------------------------------------- ! get header length endPos = index(rawData,IO_EOL) @@ -196,7 +196,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) chunkPos = IO_stringPos(trim(line)) if (chunkPos(1) < 2) cycle ! need at least one keyword value pair - + select case (IO_lc(IO_StringValue(trim(line),chunkPos,1)) ) case ('grid') if (chunkPos(1) > 6) then @@ -211,7 +211,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) end select enddo endif - + case ('size') if (chunkPos(1) > 6) then do j = 2,6,2 @@ -225,7 +225,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) end select enddo endif - + case ('origin') if (chunkPos(1) > 6) then do j = 2,6,2 @@ -258,7 +258,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) allocate(microstructure(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant) allocate(homogenization(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) - + !-------------------------------------------------------------------------------------------------- ! read and interpret content e = 1 @@ -269,7 +269,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization) startPos = endPos + 1 l = l + 1 chunkPos = IO_stringPos(trim(line)) - + noCompression: if (chunkPos(1) /= 3) then c = chunkPos(1) microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] @@ -309,7 +309,7 @@ function IPcoordinates0(grid,geomSize,grid3Offset) integer :: & a,b,c, & i - + i = 0 do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) i = i + 1 @@ -346,37 +346,37 @@ end function nodes0 !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas !-------------------------------------------------------------------------------------------------- -pure function cellEdgeArea(geomSize,grid) - +pure function cellSurfaceArea(geomSize,grid) + real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) integer, dimension(3), intent(in) :: grid ! grid (for this process!) - - real(pReal), dimension(6,1,product(grid)) :: cellEdgeArea - cellEdgeArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) - cellEdgeArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) - cellEdgeArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) - -end function cellEdgeArea + real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea + + cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) + cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) + cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + +end function cellSurfaceArea !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas normals !-------------------------------------------------------------------------------------------------- -pure function cellEdgeNormal(nElems) +pure function cellSurfaceNormal(nElems) integer, intent(in) :: nElems - - real(pReal), dimension(3,6,1,nElems) :: cellEdgeNormal - cellEdgeNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) - cellEdgeNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) - cellEdgeNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems) - cellEdgeNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems) - cellEdgeNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems) - cellEdgeNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems) - -end function cellEdgeNormal + real(pReal), dimension(3,6,1,nElems) :: cellSurfaceNormal + + cellSurfaceNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) + cellSurfaceNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) + cellSurfaceNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems) + cellSurfaceNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems) + cellSurfaceNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems) + cellSurfaceNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems) + +end function cellSurfaceNormal !-------------------------------------------------------------------------------------------------- @@ -385,9 +385,9 @@ end function cellEdgeNormal pure function IPneighborhood(grid) integer, dimension(3), intent(in) :: grid ! grid (for this process!) - + integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - + integer :: & x,y,z, & e @@ -432,7 +432,7 @@ pure function IPneighborhood(grid) IPneighborhood(3,6,1,e) = 5 enddo; enddo; enddo - + end function IPneighborhood From 9c138c87f00fe2590e2460a54a44d58cb0cfcc9f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 14:52:09 +0100 Subject: [PATCH 04/12] not needed --- src/FEsolving.f90 | 4 ++-- src/config.f90 | 2 +- src/element.f90 | 1 - 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index f81110414..fc04499fa 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -4,10 +4,10 @@ !> @brief global variables for flow control !-------------------------------------------------------------------------------------------------- module FEsolving - use prec implicit none - + public + logical :: & terminallyIll = .false. !< at least one material point is terminally ill diff --git a/src/config.f90 b/src/config.f90 index bcf9b4b3d..58b0d69a1 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -27,7 +27,7 @@ module config config_numerics, & config_debug - character(len=pStringLen), dimension(:), allocatable, public, protected :: & + character(len=pStringLen), public, protected, allocatable, dimension(:) :: & config_name_phase, & !< name of each phase config_name_homogenization, & !< name of each homogenization config_name_crystallite, & !< name of each crystallite setting diff --git a/src/element.f90 b/src/element.f90 index a38036f93..de4000700 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -3,7 +3,6 @@ !> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH !-------------------------------------------------------------------------------------------------- module element - use prec use IO implicit none From 4d07ab75198e529982b513aafec44362377968a4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 16:31:05 +0100 Subject: [PATCH 05/12] prospector complained --- python/damask/grid_filters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 3850bf70b..db9f243c9 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -366,7 +366,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): def regrid(size,F,new_grid): """ - Return mapping from coordinates in deformed configuration to a regular grid + Return mapping from coordinates in deformed configuration to a regular grid. Parameters ---------- From 9fbdb6b757a2ee7019665201b546f9bbdcf81dae Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 17:49:35 +0100 Subject: [PATCH 06/12] [skip ci] need to follow DAMASK paper --- examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config index 87b4e2312..5af1eee11 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config @@ -18,4 +18,4 @@ tau0_slip 95.e6 97.e6 # per family, optimization long tausat_slip 222.e6 412.7e6 # per family, optimization long simplex 109 h0_slipslip 1000.0e6 interaction_slipslip 1 1 1.4 1.4 1.4 1.4 -w0_slip 2.0 +a_slip 2.0 From f5bd544b36a1b2bbcd039758be44f59a5207c336 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 18:15:49 +0100 Subject: [PATCH 07/12] [skip ci] was too general --- src/mesh_grid.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 6fe326fd0..bce392d78 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -386,7 +386,7 @@ pure function IPneighborhood(grid) integer, dimension(3), intent(in) :: grid ! grid (for this process!) - integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID] integer :: & x,y,z, & From 9690f170e1bf97c13b45abad6d98f4e5fc27a845 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 22:44:42 +0100 Subject: [PATCH 08/12] 4 space indentation --- CMakeLists.txt | 16 ++++++++-------- src/CMakeLists.txt | 34 +++++++++++++++++----------------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7b4203ae2..0f488ec37 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,7 @@ cmake_minimum_required (VERSION 3.10.0 FATAL_ERROR) # Find PETSc from system environment set(PETSC_DIR $ENV{PETSC_DIR}) if (PETSC_DIR STREQUAL "") - message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined") + message (FATAL_ERROR "PETSc location (PETSC_DIR) is not defined") endif () set (petsc_conf_variables "${PETSC_DIR}/lib/petsc/conf/variables") @@ -118,7 +118,7 @@ endif () list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) if (CMAKE_BUILD_TYPE STREQUAL "") - set (CMAKE_BUILD_TYPE "RELEASE") + set (CMAKE_BUILD_TYPE "RELEASE") endif () # Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE @@ -165,13 +165,13 @@ add_definitions (-DDAMASKVERSION="${DAMASK_V}") add_definitions (-DPETSc) if (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel") - include(Compiler-Intel) + include (Compiler-Intel) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") - include(Compiler-GNU) + include (Compiler-GNU) elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "PGI") - include(Compiler-PGI) + include (Compiler-PGI) else () - message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized") + message (FATAL_ERROR "Compiler type (CMAKE_Fortran_COMPILER_ID) not recognized") endif () @@ -179,8 +179,8 @@ set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${ set (CMAKE_Fortran_LINK_EXECUTABLE "${BUILDCMD_PRE} ${PETSC_LINKER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}") if (CMAKE_BUILD_TYPE STREQUAL "DEBUG") - set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") - set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}") + set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${DEBUG_FLAGS}") + set (CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}") endif () set (CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b2eee5561..2e4336698 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,37 +4,37 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") endif() -file(GLOB damask-sources *.f90 *.c) +file (GLOB damask-sources *.f90 *.c) # probably we should have subfolders for abaqus and MSC.Marc -list(FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90") -list(FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90") -list(FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90") -list(FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90") -list(FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90") +list (FILTER damask-sources EXCLUDE REGEX ".*CPFEM.f90") +list (FILTER damask-sources EXCLUDE REGEX ".*DAMASK_marc.*.f90") +list (FILTER damask-sources EXCLUDE REGEX ".*mesh_marc.*.f90") +list (FILTER damask-sources EXCLUDE REGEX ".*mesh_abaqus.*.f90") +list (FILTER damask-sources EXCLUDE REGEX ".*commercialFEM_fileList.*.f90") if (PROJECT_NAME STREQUAL "damask-grid") - list(FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90") - file(GLOB grid-sources grid/*.f90) + list (FILTER damask-sources EXCLUDE REGEX ".*mesh_FEM.*.f90") + file (GLOB grid-sources grid/*.f90) - if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") - add_executable(DAMASK_spectral ${damask-sources} ${grid-sources}) + if (NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") + add_executable (DAMASK_spectral ${damask-sources} ${grid-sources}) install (TARGETS DAMASK_spectral RUNTIME DESTINATION bin) - else() - add_library(DAMASK_spectral OBJECT ${damask-sources} ${grid-sources}) + else() + add_library (DAMASK_spectral OBJECT ${damask-sources} ${grid-sources}) exec_program (mktemp OUTPUT_VARIABLE nothing) exec_program (mktemp ARGS -d OUTPUT_VARIABLE black_hole) install (PROGRAMS ${nothing} DESTINATION ${black_hole}) - endif() + endif() elseif (PROJECT_NAME STREQUAL "damask-mesh") - list(FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90") - file(GLOB mesh-sources mesh/*.f90) + list (FILTER damask-sources EXCLUDE REGEX ".*mesh_grid.*.f90") + file (GLOB mesh-sources mesh/*.f90) - add_executable(DAMASK_FEM ${damask-sources} ${mesh-sources}) - install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin) + add_executable (DAMASK_FEM ${damask-sources} ${mesh-sources}) + install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin) endif() From 59fe9d06b0f8a915521deabcad296d1c832346a1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 22:53:19 +0100 Subject: [PATCH 09/12] shortening --- src/IO.f90 | 21 ++++++++++----------- src/rotations.f90 | 4 ++-- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 9406c349b..5c83876c5 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -400,15 +400,14 @@ pure function IO_lc(string) character(len=*), intent(in) :: string !< string to convert character(len=len(string)) :: IO_lc - character(26), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' - character(26), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' + character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' - integer :: i,n + integer :: i,n do i=1,len(string) - IO_lc(i:i) = string(i:i) n = index(UPPER,IO_lc(i:i)) - if (n/=0) IO_lc(i:i) = LOWER(n:n) + IO_lc(i:i) = merge(LOWER(n:n),string(i:i),n/=0) enddo end function IO_lc @@ -549,6 +548,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) ! math errors case (400) msg = 'matrix inversion error' + case (401) + msg = 'error in Eigenvalue calculation' case (402) msg = 'invalid orientation specified' @@ -876,11 +877,10 @@ integer function verifyIntValue(string) character(len=*), intent(in) :: string !< string for conversion to int value - integer :: readStatus, invalidWhere + integer :: readStatus character(len=*), parameter :: VALIDCHARS = '0123456789+- ' - invalidWhere = verify(string,VALIDCHARS) - valid: if (invalidWhere == 0) then + valid: if (verify(string,VALIDCHARS) == 0) then read(string,*,iostat=readStatus) verifyIntValue if (readStatus /= 0) call IO_error(111,ext_msg=string) else valid @@ -898,11 +898,10 @@ real(pReal) function verifyFloatValue(string) character(len=*), intent(in) :: string !< string for conversion to float value - integer :: readStatus, invalidWhere + integer :: readStatus character(len=*), parameter :: VALIDCHARS = '0123456789eE.+- ' - invalidWhere = verify(string,VALIDCHARS) - valid: if (invalidWhere == 0) then + valid: if (verify(string,VALIDCHARS) == 0) then read(string,*,iostat=readStatus) verifyFloatValue if (readStatus /= 0) call IO_error(112,ext_msg=string) else valid diff --git a/src/rotations.f90 b/src/rotations.f90 index 16e63c180..c1fb1f2cf 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -608,13 +608,13 @@ function om2ax(om) result(ax) ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ] else call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr) - if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax: DGEEV return not zero') + if (ierr /= 0) call IO_error(401,ext_msg='Error in om2ax: DGEEV return not zero') #if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI) i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) #else i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0) #endif - if (i == 0) call IO_error(0,ext_msg='Error in om2ax Real: eigenvalue not found') + if (i == 0) call IO_error(401,ext_msg='Error in om2ax Real: eigenvalue not found') ax(1:3) = VR(1:3,i) where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) From c2cdcb17f7ffb73d1bb839045bc4c985a54e1dfd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 23:04:15 +0100 Subject: [PATCH 10/12] wrong string --- src/IO.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IO.f90 b/src/IO.f90 index 5c83876c5..9f23fd6d1 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -406,7 +406,7 @@ pure function IO_lc(string) integer :: i,n do i=1,len(string) - n = index(UPPER,IO_lc(i:i)) + n = index(UPPER,string(i:i)) IO_lc(i:i) = merge(LOWER(n:n),string(i:i),n/=0) enddo From d54b8714e1a3724f52aaddadc6718158e4e3f07d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Jan 2020 23:42:50 +0100 Subject: [PATCH 11/12] avoid invalid error access --- src/IO.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/IO.f90 b/src/IO.f90 index 9f23fd6d1..da3dcbf3c 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -407,7 +407,11 @@ pure function IO_lc(string) do i=1,len(string) n = index(UPPER,string(i:i)) - IO_lc(i:i) = merge(LOWER(n:n),string(i:i),n/=0) + if(n/=0) then + IO_lc(i:i) = LOWER(n:n) + else + IO_lc(i:i) = string(i:i) + endif enddo end function IO_lc From 0b8ff64884f2c0509edefc520192deb2f5c64ecf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Jan 2020 23:33:33 +0100 Subject: [PATCH 12/12] store mapping MARC/FEM2DAMASK mapping do not calculate the mapping for elements and nodes per call on the fly, rather store it. Not memory efficient in the case that numbers are not consequtive (order does not matter, but missing nodes/elements would waste some 2 integers per missing number). However, this seem to cause problems anyway when range indicators like '1 to 10' are used. --- src/CPFEM.f90 | 2 +- src/DAMASK_marc.f90 | 4 +- src/mesh_marc.f90 | 126 +++++++++++++------------------------------- 3 files changed, 41 insertions(+), 91 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 1ef5f1389..fd406219f 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -152,7 +152,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt i, j, k, l, m, n, ph, homog, mySource logical updateJaco ! flag indicating if JAcobian has to be updated - elCP = mesh_FEasCP('elem',elFE) + elCP = mesh_FEM2DAMASK_elem(elFE) if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt & .and. elCP == debug_e .and. ip == debug_i) then diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ba7751372..f7efb5279 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -269,7 +269,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback computationMode = CPFEM_RESTOREJACOBIAN elseif (lovl == 6) then ! stress requested by marc - cp_en = mesh_FEasCP('elem',m(1)) + cp_en = mesh_FEM2DAMASK_elem(m(1)) if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" terminallyIll = .false. cycleCounter = -1 ! first calc step increments this to cycle = 0 @@ -391,7 +391,7 @@ subroutine flux(f,ts,n,time) real(pReal), dimension(2), intent(out) :: & f - call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1))) + call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 05bb77b8a..14e001eaf 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -32,17 +32,12 @@ module mesh real(pReal), public, protected :: & mesh_unitlength !< physical length of one unit in mesh - integer, dimension(:,:), allocatable, target :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEnode, corresponding CPnode] - - integer, dimension(:), allocatable :: & - mapMarc2DAMASK_elem, & !< DAMASK element ID for Marc element ID - mapMarc2DAMASK_node !< DAMASK node ID for Marc node ID + integer, dimension(:), allocatable, public :: & + mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID + mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID public :: & - mesh_init, & - mesh_FEasCP + mesh_init contains @@ -89,7 +84,7 @@ subroutine mesh_init(ip,el) FEsolving_execIP = [1,elem%nIPs] allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" + calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc" allocate(cellNodeDefinition(elem%nNodes-1)) @@ -215,13 +210,11 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_elemType(elem, & nElems,inputFile) - allocate (mesh_mapFEtoCPelem(2,nElems), source=0) - call inputRead_mapElems(mapMarc2DAMASK_elem,& - elem%nNodes,inputFile) + call inputRead_mapElems(mesh_FEM2DAMASK_elem,& + nElems,elem%nNodes,inputFile) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) - call inputRead_mapNodes(mapMarc2DAMASK_node,& - inputFile) + call inputRead_mapNodes(mesh_FEM2DAMASK_node,& + nNodes,inputFile) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) @@ -432,14 +425,16 @@ end subroutine inputRead_mapElemSets !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapElems(map, & - nNodes,fileContent) +subroutine inputRead_mapElems(FEM2DAMASK, & + nElems,nNodesPerElem,fileContent) - integer, allocatable, dimension(:), intent(out) :: map + integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK - integer, intent(in) :: nNodes !< number of nodes per element + integer, intent(in) :: nElems, & !< number of elements + nNodesPerElem !< number of nodes per element character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, dimension(2,nElems) :: map_unsorted integer, allocatable, dimension(:) :: chunkPos integer :: i,j,l,nNodesAlreadyRead @@ -448,11 +443,11 @@ subroutine inputRead_mapElems(map, & if(chunkPos(1) < 1) cycle if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity') then j = 0 - do i = 1,size(mesh_mapFEtoCPelem,2) + do i = 1,nElems chunkPos = IO_stringPos(fileContent(l+1+i+j)) - mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i] + map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i] nNodesAlreadyRead = chunkPos(1) - 2 - do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line + do while(nNodesAlreadyRead < nNodesPerElem) ! read on if not all nodes in one line j = j + 1 chunkPos = IO_stringPos(fileContent(l+1+i+j)) nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) @@ -462,10 +457,10 @@ subroutine inputRead_mapElems(map, & endif enddo - call math_sort(mesh_mapFEtoCPelem) - allocate(map(minval(mesh_mapFEtoCPelem(1,:)):maxval(mesh_mapFEtoCPelem(1,:))),source=-1) - do i = 1,size(mesh_mapFEtoCPelem,2) - map(mesh_mapFEtoCPelem(1,i)) = mesh_mapFEtoCPelem(2,i) + call math_sort(map_unsorted) + allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1) + do i = 1, nElems + FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i) enddo end subroutine inputRead_mapElems @@ -474,13 +469,15 @@ end subroutine inputRead_mapElems !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapNodes(map, & - fileContent) +subroutine inputRead_mapNodes(FEM2DAMASK, & + nNodes,fileContent) - integer, allocatable, dimension(:), intent(out) :: map + integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK + integer, intent(in) :: nNodes !< number of nodes character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines + integer, dimension(2,nNodes) :: map_unsorted integer, allocatable, dimension(:) :: chunkPos integer :: i, l @@ -488,18 +485,18 @@ subroutine inputRead_mapNodes(map, & chunkPos = IO_stringPos(fileContent(l)) if(chunkPos(1) < 1) cycle if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then - do i = 1,size(mesh_mapFEtoCPnode,2) + do i = 1,nNodes chunkPos = IO_stringPos(fileContent(l+1+i)) - mesh_mapFEtoCPnode(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i] + map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i] enddo exit endif enddo - call math_sort(mesh_mapFEtoCPnode) - allocate(map(minval(mesh_mapFEtoCPnode(1,:)):maxval(mesh_mapFEtoCPnode(1,:))),source=-1) - do i = 1,size(mesh_mapFEtoCPnode,2) - map(mesh_mapFEtoCPnode(1,i)) = mesh_mapFEtoCPnode(2,i) + call math_sort(map_unsorted) + allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1) + do i = 1, nNodes + FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i) enddo end subroutine inputRead_mapNodes @@ -526,7 +523,7 @@ subroutine inputRead_elemNodes(nodes, & if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then do i=1,nNode chunkPos = IO_stringPos(fileContent(l+1+i)) - m = mesh_FEasCP('node',IO_intValue(fileContent(l+1+i),chunkPos,1)) + m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1)) do j = 1,3 nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) enddo @@ -650,11 +647,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent) j = 0 do i = 1,nElem chunkPos = IO_stringPos(fileContent(l+1+i+j)) - e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1)) + e = mesh_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1)) if (e /= 0) then ! disregard non CP elems do k = 1,chunkPos(1)-2 inputRead_connectivityElem(k,e) = & - mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) + mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) enddo nNodesAlreadyRead = chunkPos(1) - 2 do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line @@ -662,7 +659,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent) chunkPos = IO_stringPos(fileContent(l+1+i+j)) do k = 1,chunkPos(1) inputRead_connectivityElem(nNodesAlreadyRead+k,e) = & - mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) + mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) enddo nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -717,7 +714,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza if (initialcondTableStyle == 2) m = m + 2 contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements do i = 1,contInts(1) - e = mesh_FEasCP('elem',contInts(1+i)) + e = mesh_FEM2DAMASK_elem(contInts(1+i)) if (sv == 2) microstructureAt(e) = myVal if (sv == 3) homogenizationAt(e) = myVal enddo @@ -1050,53 +1047,6 @@ function IPareaNormal(elem,nElem,connectivity,node) end function IPareaNormal -!-------------------------------------------------------------------------------------------------- -!> @brief Gives the FE to CP ID mapping by binary search through lookup array -!! valid questions (what) are 'elem', 'node' -!-------------------------------------------------------------------------------------------------- -integer function mesh_FEasCP(what,myID) - - character(len=*), intent(in) :: what - integer, intent(in) :: myID - - integer, dimension(:,:), pointer :: lookupMap - integer :: lower,upper,center - - mesh_FEasCP = 0 - select case(IO_lc(what(1:4))) - case('elem') - lookupMap => mesh_mapFEtoCPelem - case('node') - lookupMap => mesh_mapFEtoCPnode - case default - return - endselect - - lower = 1 - upper = int(size(lookupMap,2),pInt) - - if (lookupMap(1,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds? - mesh_FEasCP = lookupMap(2,lower) - return - elseif (lookupMap(1,upper) == myID) then - mesh_FEasCP = lookupMap(2,upper) - return - endif - binarySearch: do while (upper-lower > 1) - center = (lower+upper)/2 - if (lookupMap(1,center) < myID) then - lower = center - elseif (lookupMap(1,center) > myID) then - upper = center - else - mesh_FEasCP = lookupMap(2,center) - exit - endif - enddo binarySearch - -end function mesh_FEasCP - - !-------------------------------------------------------------------------------------------------- !> @brief return integer list corresponding to items in consecutive lines. !! First integer in array is counter