Merge branch 'MiscImprovements' of magit1.mpie.de:/damask/DAMASK into MiscImprovements

This commit is contained in:
Martin Diehl 2020-01-31 00:13:25 +01:00
commit 6ce4ce523a
15 changed files with 255 additions and 258 deletions

View File

@ -165,11 +165,11 @@ 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")
endif ()

View File

@ -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

View File

@ -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()

View File

@ -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()))

View File

@ -4,26 +4,26 @@ 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})
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})
@ -31,10 +31,10 @@ if (PROJECT_NAME STREQUAL "damask-grid")
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})
add_executable (DAMASK_FEM ${damask-sources} ${mesh-sources})
install (TARGETS DAMASK_FEM RUNTIME DESTINATION bin)
endif()

View File

@ -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

View File

@ -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

View File

@ -4,9 +4,9 @@
!> @brief global variables for flow control
!--------------------------------------------------------------------------------------------------
module FEsolving
use prec
implicit none
public
logical :: &
terminallyIll = .false. !< at least one material point is terminally ill

View File

@ -400,15 +400,18 @@ 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
do i=1,len(string)
n = index(UPPER,string(i:i))
if(n/=0) then
IO_lc(i:i) = LOWER(n:n)
else
IO_lc(i:i) = string(i:i)
n = index(UPPER,IO_lc(i:i))
if (n/=0) IO_lc(i:i) = LOWER(n:n)
endif
enddo
end function IO_lc
@ -549,6 +552,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 +881,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 +902,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

View File

@ -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

View File

@ -3,7 +3,6 @@
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!--------------------------------------------------------------------------------------------------
module element
use prec
use IO
implicit none

View File

@ -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, &
@ -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
@ -351,7 +356,9 @@ real(pReal) pure function math_mul33xx33(A,B)
integer :: i,j
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) :: math_mul3333xx33
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
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
@ -382,8 +391,9 @@ pure function math_mul3333xx3333(A,B)
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
@ -1308,6 +1317,10 @@ 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] - &
@ -1322,7 +1335,6 @@ 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')
@ -1330,6 +1342,11 @@ subroutine unitTest
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))) &
@ -1364,6 +1381,9 @@ subroutine unitTest
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
@ -1384,10 +1404,17 @@ subroutine unitTest
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')

View File

@ -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))
!--------------------------------------------------------------------------------------------------
@ -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
real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea
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))
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 cellEdgeArea
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
real(pReal), dimension(3,6,1,nElems) :: cellSurfaceNormal
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)
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 cellEdgeNormal
end function cellSurfaceNormal
!--------------------------------------------------------------------------------------------------
@ -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, &

View File

@ -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

View File

@ -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)])