223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch '29-rename-mesh_element-array' into development

This commit is contained in:
Martin Diehl 2018-10-10 15:55:23 +02:00
parent 5038578188 f0b9c0caf7
commit 7217cdac1b
9 changed files with 247 additions and 252 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/CPFEM.f90
View File

@ -98,7 +98,7 @@ subroutine CPFEM_initAll(el,ip)
call config_init call config_init
call math_init call math_init
call FE_init call FE_init
call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip call mesh_init(ip, el)
call lattice_init call lattice_init
call material_init call material_init
call constitutive_init call constitutive_init
@ -314,7 +314,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
thermal_type, & thermal_type, &
THERMAL_conduction_ID, & THERMAL_conduction_ID, &
phase_Nsources, & phase_Nsources, &
material_homog material_homogenizationAt
use config, only: & use config, only: &
material_Nhomogenization material_Nhomogenization
use crystallite, only: & use crystallite, only: &
@ -503,7 +503,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (.not. parallelExecution) then if (.not. parallelExecution) then
chosenThermal1: select case (thermal_type(mesh_element(3,elCP))) chosenThermal1: select case (thermal_type(mesh_element(3,elCP)))
case (THERMAL_conduction_ID) chosenThermal1 case (THERMAL_conduction_ID) chosenThermal1
temperature(material_homog(ip,elCP))%p(thermalMapping(material_homog(ip,elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature_inp temperature_inp
end select chosenThermal1 end select chosenThermal1
materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn
@ -516,7 +516,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6)
chosenThermal2: select case (thermal_type(mesh_element(3,elCP))) chosenThermal2: select case (thermal_type(mesh_element(3,elCP)))
case (THERMAL_conduction_ID) chosenThermal2 case (THERMAL_conduction_ID) chosenThermal2
temperature(material_homog(ip,elCP))%p(thermalMapping(material_homog(ip,elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature_inp temperature_inp
end select chosenThermal2 end select chosenThermal2
materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn

6
src/CPFEM2.f90
View File

@ -18,7 +18,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief call (thread safe) all module initializations !> @brief call (thread safe) all module initializations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll(el,ip) subroutine CPFEM_initAll()
use prec, only: & use prec, only: &
pInt pInt
use prec, only: & use prec, only: &
@ -55,8 +55,6 @@ subroutine CPFEM_initAll(el,ip)
#endif #endif
implicit none implicit none
integer(pInt), intent(in) :: el, & !< FE el number
ip !< FE integration point number
call DAMASK_interface_init ! Spectral and FEM interface to commandline call DAMASK_interface_init ! Spectral and FEM interface to commandline
call prec_init call prec_init
@ -69,7 +67,7 @@ subroutine CPFEM_initAll(el,ip)
call config_init call config_init
call math_init call math_init
call FE_init call FE_init
call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip call mesh_init
call lattice_init call lattice_init
call material_init call material_init
call constitutive_init call constitutive_init

2
src/DAMASK_FEM.f90
View File

@ -116,7 +116,7 @@ program DAMASK_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"

2
src/DAMASK_spectral.f90
View File

@ -151,7 +151,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'
write(6,'(/,a,/)') ' Roters et al., Computational Materials Science, 2018' write(6,'(/,a,/)') ' Roters et al., Computational Materials Science, 2018'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()

2
src/IO.f90
View File

@ -1488,6 +1488,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'no microstructure specified via State Variable 3' msg = 'no microstructure specified via State Variable 3'
case (190_pInt) case (190_pInt)
msg = 'unknown element type:' msg = 'unknown element type:'
case (191_pInt)
msg = 'mesh consists of more than one element type'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! plasticity error messages ! plasticity error messages

21
src/constitutive.f90
View File

@ -386,7 +386,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase, & material_phase, &
material_homog, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
@ -413,7 +413,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
real(pReal), intent(in), dimension(:,:,:,:) :: & real(pReal), intent(in), dimension(:,:,:,:) :: &
orientations !< crystal orientations as quaternions orientations !< crystal orientations as quaternions
ho = material_homog(ip,el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -444,7 +444,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
phase_plasticity, & phase_plasticity, &
phase_plasticityInstance, & phase_plasticityInstance, &
material_phase, & material_phase, &
material_homog, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
@ -494,7 +494,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
integer(pInt) :: & integer(pInt) :: &
i, j, instance, of i, j, instance, of
ho = material_homog(ip,el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
S = math_Mandel6to33(S6) S = math_Mandel6to33(S6)
@ -760,7 +760,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
math_I3 math_I3
use material, only: & use material, only: &
material_phase, & material_phase, &
material_homog, & material_homogenizationAt, &
phase_NstiffnessDegradations, & phase_NstiffnessDegradations, &
phase_stiffnessDegradation, & phase_stiffnessDegradation, &
damage, & damage, &
@ -791,8 +791,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
integer(pInt) :: & integer(pInt) :: &
i, j i, j
ho = material_homog(ip,el) ho = material_homogenizationAt(el)
C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el)) C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el)) DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
@ -843,7 +842,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
phase_source, & phase_source, &
phase_Nsources, & phase_Nsources, &
material_phase, & material_phase, &
material_homog, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
@ -903,7 +902,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
s, & !< counter in source loop s, & !< counter in source loop
instance, of instance, of
ho = material_homog( ip,el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
@ -1062,7 +1061,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
phase_source, & phase_source, &
phase_Nsources, & phase_Nsources, &
material_phase, & material_phase, &
material_homog, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
@ -1125,7 +1124,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
ho = material_homog( ip,el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
startPos = 1_pInt startPos = 1_pInt

70
src/material.f90
View File

@ -169,6 +169,7 @@ module material
homogenization_maxNgrains !< max number of grains in any USED homogenization homogenization_maxNgrains !< max number of grains in any USED homogenization
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
material_homogenizationAt, & !< homogenization ID of each element (copy of mesh_homogenizationAt)
phase_Nsources, & !< number of source mechanisms active in each phase phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
@ -199,8 +200,10 @@ module material
integer(pInt), dimension(:,:,:), allocatable, public :: & integer(pInt), dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element material_phase !< phase (index) of each grain,IP,element
! BEGIN DEPRECATED: use material_homogenizationAt
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public :: &
material_homog !< homogenization (index) of each IP,element material_homog !< homogenization (index) of each IP,element
! END DEPRECATED
type(tPlasticState), allocatable, dimension(:), public :: & type(tPlasticState), allocatable, dimension(:), public :: &
plasticState plasticState
type(tSourceState), allocatable, dimension(:), public :: & type(tSourceState), allocatable, dimension(:), public :: &
@ -362,10 +365,10 @@ subroutine material_init()
phase_name, & phase_name, &
texture_name texture_name
use mesh, only: & use mesh, only: &
mesh_homogenizationAt, &
mesh_NipsPerElem, &
mesh_maxNips, & mesh_maxNips, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_element, &
FE_Nips, &
FE_geomtype FE_geomtype
implicit none implicit none
@ -480,11 +483,11 @@ subroutine material_init()
allocate(CrystallitePosition (size(config_phase)), source=0_pInt) allocate(CrystallitePosition (size(config_phase)), source=0_pInt)
ElemLoop:do e = 1_pInt,mesh_NcpElems ElemLoop:do e = 1_pInt,mesh_NcpElems
myHomog = mesh_element(3,e) myHomog = mesh_homogenizationAt(e)
IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) IPloop:do i = 1_pInt, mesh_NipsPerElem
HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt
mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog] mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog]
GrainLoop:do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) GrainLoop:do g = 1_pInt,homogenization_Ngrains(myHomog)
phase = material_phase(g,i,e) phase = material_phase(g,i,e)
ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase
phaseAt(g,i,e) = phase phaseAt(g,i,e) = phase
@ -519,10 +522,10 @@ end subroutine material_init
subroutine material_parseHomogenization subroutine material_parseHomogenization
use config, only : & use config, only : &
config_homogenization config_homogenization
use mesh, only: &
mesh_homogenizationAt
use IO, only: & use IO, only: &
IO_error IO_error
use mesh, only: &
mesh_element
implicit none implicit none
integer(pInt) :: h integer(pInt) :: h
@ -549,7 +552,8 @@ subroutine material_parseHomogenization
allocate(porosity_initialPhi(size(config_homogenization)), source=1.0_pReal) allocate(porosity_initialPhi(size(config_homogenization)), source=1.0_pReal)
allocate(hydrogenflux_initialCh(size(config_homogenization)), source=0.0_pReal) allocate(hydrogenflux_initialCh(size(config_homogenization)), source=0.0_pReal)
forall (h = 1_pInt:size(config_homogenization)) homogenization_active(h) = any(mesh_element(3,:) == h) forall (h = 1_pInt:size(config_homogenization)) &
homogenization_active(h) = any(mesh_homogenizationAt == h)
do h=1_pInt, size(config_homogenization) do h=1_pInt, size(config_homogenization)
@ -685,7 +689,7 @@ subroutine material_parseMicrostructure
config_microstructure, & config_microstructure, &
microstructure_name microstructure_name
use mesh, only: & use mesh, only: &
mesh_element, & mesh_microstructureAt, &
mesh_NcpElems mesh_NcpElems
implicit none implicit none
@ -701,10 +705,11 @@ subroutine material_parseMicrostructure
allocate(microstructure_active(size(config_microstructure)), source=.false.) allocate(microstructure_active(size(config_microstructure)), source=.false.)
allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.) allocate(microstructure_elemhomo(size(config_microstructure)), source=.false.)
if(any(mesh_element(4,1:mesh_NcpElems) > size(config_microstructure))) & if(any(mesh_microstructureAt > size(config_microstructure))) &
call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config') call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config')
forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements forall (e = 1_pInt:mesh_NcpElems) &
microstructure_active(mesh_microstructureAt(e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements
do m=1_pInt, size(config_microstructure) do m=1_pInt, size(config_microstructure)
microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)') microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)')
@ -1082,11 +1087,13 @@ subroutine material_populateGrains
math_sampleFiberOri, & math_sampleFiberOri, &
math_symmetricEulers math_symmetricEulers
use mesh, only: & use mesh, only: &
mesh_element, & mesh_NipsPerElem, &
mesh_elemType, &
mesh_homogenizationAt, &
mesh_microstructureAt, &
mesh_maxNips, & mesh_maxNips, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_ipVolume, & mesh_ipVolume, &
FE_Nips, &
FE_geomtype FE_geomtype
use config, only: & use config, only: &
config_homogenization, & config_homogenization, &
@ -1127,6 +1134,7 @@ subroutine material_populateGrains
allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0.0_pReal)
allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_homog(mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_homog(mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_homogenizationAt,source=mesh_homogenizationAt)
allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
@ -1136,18 +1144,14 @@ subroutine material_populateGrains
! populating homogenization schemes in each ! populating homogenization schemes in each
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
do e = 1_pInt, mesh_NcpElems do e = 1_pInt, mesh_NcpElems
material_homog(1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))),e) = mesh_element(3,e) material_homog(1_pInt:mesh_NipsPerElem,e) = mesh_homogenizationAt(e)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! precounting of elements for each homog/micro pair ! precounting of elements for each homog/micro pair
do e = 1_pInt, mesh_NcpElems do e = 1_pInt, mesh_NcpElems
homog = mesh_element(3,e) homog = mesh_homogenizationAt(e)
micro = mesh_element(4,e) micro = mesh_microstructureAt(e)
if (homog < 1_pInt .or. homog > size(Nelems,1)) &
call IO_error(154_pInt,e,0_pInt,0_pInt)
if (micro < 1_pInt .or. micro > size(Nelems,2)) &
call IO_error(155_pInt,e,0_pInt,0_pInt)
Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
enddo enddo
allocate(elemsOfHomogMicro(size(config_homogenization),size(config_microstructure))) allocate(elemsOfHomogMicro(size(config_homogenization),size(config_microstructure)))
@ -1164,13 +1168,17 @@ subroutine material_populateGrains
! identify maximum grain count per IP (from element) and find grains per homog/micro pair ! identify maximum grain count per IP (from element) and find grains per homog/micro pair
Nelems = 0_pInt ! reuse as counter Nelems = 0_pInt ! reuse as counter
elementLooping: do e = 1_pInt,mesh_NcpElems elementLooping: do e = 1_pInt,mesh_NcpElems
t = FE_geomtype(mesh_element(2,e)) t = mesh_elemType
homog = mesh_element(3,e) homog = mesh_homogenizationAt(e)
micro = mesh_element(4,e) micro = mesh_microstructureAt(e)
if (homog < 1_pInt .or. homog > size(config_homogenization)) & ! out of bounds
call IO_error(154_pInt,e,0_pInt,0_pInt)
if (micro < 1_pInt .or. micro > size(config_microstructure)) & ! out of bounds
call IO_error(155_pInt,e,0_pInt,0_pInt)
if (microstructure_elemhomo(micro)) then ! how many grains are needed at this element? if (microstructure_elemhomo(micro)) then ! how many grains are needed at this element?
dGrains = homogenization_Ngrains(homog) ! only one set of Ngrains (other IPs are plain copies) dGrains = homogenization_Ngrains(homog) ! only one set of Ngrains (other IPs are plain copies)
else else
dGrains = homogenization_Ngrains(homog) * FE_Nips(t) ! each IP has Ngrains dGrains = homogenization_Ngrains(homog) * mesh_NipsPerElem ! each IP has Ngrains
endif endif
Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains ! total grain count Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains ! total grain count
Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt ! total element count Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt ! total element count
@ -1204,16 +1212,16 @@ subroutine material_populateGrains
do hme = 1_pInt, Nelems(homog,micro) do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(homog,micro)%p(hme) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex e = elemsOfHomogMicro(homog,micro)%p(hme) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
t = FE_geomtype(mesh_element(2,e)) t = mesh_elemType
if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs
volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(t),e))/& volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:mesh_NipsPerElem,e))/&
real(dGrains,pReal) ! each grain combines size of all IPs in that element real(dGrains,pReal) ! each grain combines size of all IPs in that element
grain = grain + dGrains ! wind forward by Ngrains@IP grain = grain + dGrains ! wind forward by Ngrains@IP
else else
forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs forall (i = 1_pInt:mesh_NipsPerElem) & ! loop over IPs
volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = & volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
mesh_ipVolume(i,e)/real(dGrains,pReal) ! assign IPvolume/Ngrains@IP to all grains of IP mesh_ipVolume(i,e)/real(dGrains,pReal) ! assign IPvolume/Ngrains@IP to all grains of IP
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*Ngrains@IP grain = grain + mesh_NipsPerElem * dGrains ! wind forward by Nips*Ngrains@IP
endif endif
enddo enddo
@ -1367,11 +1375,11 @@ subroutine material_populateGrains
do hme = 1_pInt, Nelems(homog,micro) do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(homog,micro)%p(hme) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex e = elemsOfHomogMicro(homog,micro)%p(hme) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
t = FE_geomtype(mesh_element(2,e)) t = mesh_elemType
if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs
m = 1_pInt ! process only first IP m = 1_pInt ! process only first IP
else else
m = FE_Nips(t) ! process all IPs m = mesh_NipsPerElem
endif endif
do i = 1_pInt, m ! loop over necessary IPs do i = 1_pInt, m ! loop over necessary IPs
@ -1409,7 +1417,7 @@ subroutine material_populateGrains
enddo enddo
do i = i, FE_Nips(t) ! loop over IPs to (possibly) distribute copies from first IP do i = i, mesh_NipsPerElem ! loop over IPs to (possibly) distribute copies from first IP
material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e) material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e)
material_phase (1_pInt:dGrains,i,e) = material_phase (1_pInt:dGrains,1,e) material_phase (1_pInt:dGrains,i,e) = material_phase (1_pInt:dGrains,1,e)
material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e) material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e)

312
src/mesh.f90
View File

@ -3,7 +3,6 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH !> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module mesh module mesh
@ -14,35 +13,27 @@ module mesh
private private
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
mesh_NcpElems, & !< total number of CP elements in local mesh mesh_NcpElems, & !< total number of CP elements in local mesh
mesh_NelemSets, & mesh_elemType, & !< Element type of the mesh (only support homogeneous meshes)
mesh_maxNelemInSet, &
mesh_Nmaterials, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates) mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates)
mesh_Ncells, & !< total number of cells in mesh mesh_Ncells, & !< total number of cells in mesh
mesh_maxNnodes, & !< max number of nodes in any CP element mesh_NipsPerElem, & !< number of IPs in per element
mesh_maxNips, & !< max number of IPs in any CP element mesh_NcellnodesPerElem, & !< number of cell nodes per element
mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element
mesh_maxNsharedElems, & !< max number of CP elements sharing a node mesh_maxNsharedElems !< max number of CP elements sharing a node
mesh_maxNcellnodes, & !< max number of cell nodes in any CP element !!!! BEGIN DEPRECATED !!!!!
mesh_Nelems !< total number of elements in mesh
#ifdef Spectral
integer(pInt), dimension(3), public, protected :: &
grid !< (global) grid
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh mesh_maxNips, & !< max number of IPs in any CP element
grid3, & !< (local) grid in 3rd direction mesh_maxNcellnodes !< max number of cell nodes in any CP element
grid3Offset !< (local) grid offset in 3rd direction !!!! BEGIN DEPRECATED !!!!!
real(pReal), dimension(3), public, protected :: &
geomSize integer(pInt), dimension(:), allocatable, public, protected :: &
real(pReal), public, protected :: & mesh_homogenizationAt, & !< homogenization ID of each element
size3, & !< (local) size in 3rd direction mesh_microstructureAt !< microstructure ID of each element
size3offset !< (local) size offset in 3rd direction
#endif
integer(pInt), dimension(:,:), allocatable, public, protected :: & integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element, & !< FEid, type(internal representation), material, texture, node indices as CP IDs mesh_CPnodeID, & !< nodes forming an element
mesh_element, & !DEPRECATED
mesh_sharedElem, & !< entryCount and list of elements containing node mesh_sharedElem, & !< entryCount and list of elements containing node
mesh_nodeTwins !< node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions) mesh_nodeTwins !< node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions)
@ -71,36 +62,18 @@ module mesh
logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes) logical, dimension(3), public, protected :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
#ifdef Marc4DAMASK #if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), private :: & integer(pInt), private :: &
MarcVersion, & !< Version of input file format (Marc only) mesh_maxNelemInSet, &
hypoelasticTableStyle, & !< Table style (Marc only) mesh_Nmaterials
initialcondTableStyle !< Table style (Marc only)
integer(pInt), dimension(:), allocatable, private :: &
Marc_matNumber !< array of material numbers for hypoelastic material (Marc only)
#endif #endif
integer(pInt), dimension(2), private :: & integer(pInt), dimension(2), private :: &
mesh_maxValStateVar = 0_pInt mesh_maxValStateVar = 0_pInt
#ifndef Spectral integer(pInt), dimension(:,:), allocatable, private :: &
character(len=64), dimension(:), allocatable, private :: &
mesh_nameElemSet, & !< names of elementSet
mesh_nameMaterial, & !< names of material in solid section
mesh_mapMaterial !< name of elementSet for material
integer(pInt), dimension(:,:), allocatable, private :: &
mesh_mapElemSet !< list of elements in elementSet
#endif
integer(pInt), dimension(:,:), allocatable, private :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
#if defined(Marc4DAMASK) || defined(Abaqus)
integer(pInt), dimension(:,:), allocatable, target, private :: &
mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
#endif
integer(pInt),dimension(:,:,:), allocatable, private :: & integer(pInt),dimension(:,:,:), allocatable, private :: &
mesh_cell !< cell connectivity for each element,ip/cell mesh_cell !< cell connectivity for each element,ip/cell
@ -116,10 +89,6 @@ module mesh
integer(pInt), dimension(:,:,:,:), allocatable, private :: & integer(pInt), dimension(:,:,:,:), allocatable, private :: &
FE_subNodeOnIPFace FE_subNodeOnIPFace
#ifdef Abaqus
logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
#endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) ! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_" ! Hence, I suggest to prefix with "FE_"
@ -375,60 +344,81 @@ module mesh
4 & ! element 21 (3D 20node 27ip) 4 & ! element 21 (3D 20node 27ip)
],pInt) ],pInt)
#if defined(Spectral)
! integer(pInt), dimension(FE_Nelemtypes), parameter, private :: MESH_VTKELEMTYPE = & integer(pInt), dimension(3), public, protected :: &
! int([ & grid !< (global) grid
! 5, & ! element 6 (2D 3node 1ip) integer(pInt), public, protected :: &
! 22, & ! element 125 (2D 6node 3ip) mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh
! 9, & ! element 11 (2D 4node 4ip) grid3, & !< (local) grid in 3rd direction
! 23, & ! element 27 (2D 8node 9ip) grid3Offset !< (local) grid offset in 3rd direction
! 23, & ! element 54 (2D 8node 4ip) real(pReal), dimension(3), public, protected :: &
! 10, & ! element 134 (3D 4node 1ip) geomSize
! 10, & ! element 157 (3D 5node 4ip) real(pReal), public, protected :: &
! 24, & ! element 127 (3D 10node 4ip) size3, & !< (local) size in 3rd direction
! 13, & ! element 136 (3D 6node 6ip) size3offset !< (local) size offset in 3rd direction
! 12, & ! element 117 (3D 8node 1ip) #elif defined(Marc4DAMASK) || defined(Abaqus)
! 12, & ! element 7 (3D 8node 8ip) integer(pInt), private :: &
! 25, & ! element 57 (3D 20node 8ip) mesh_Nelems, & !< total number of elements in mesh (including non-DAMASK elements)
! 25 & ! element 21 (3D 20node 27ip) mesh_maxNnodes, & !< max number of nodes in any CP element
! ],pInt) mesh_NelemSets
! character(len=64), dimension(:), allocatable, private :: &
! integer(pInt), dimension(FE_Ncelltypes), parameter, private :: MESH_VTKCELLTYPE = & mesh_nameElemSet, & !< names of elementSet
! int([ & mesh_nameMaterial, & !< names of material in solid section
! 5, & ! (2D 3node) mesh_mapMaterial !< name of elementSet for material
! 9, & ! (2D 4node) integer(pInt), dimension(:,:), allocatable, private :: &
! 10, & ! (3D 4node) mesh_mapElemSet !< list of elements in elementSet
! 12 & ! (3D 8node) integer(pInt), dimension(:,:), allocatable, target, private :: &
! ],pInt) mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid]
mesh_mapFEtoCPnode !< [sorted FEid, corresponding CPid]
#endif
#if defined(Marc4DAMASK)
integer(pInt), private :: &
MarcVersion, & !< Version of input file format (Marc only)
hypoelasticTableStyle, & !< Table style (Marc only)
initialcondTableStyle !< Table style (Marc only)
integer(pInt), dimension(:), allocatable, private :: &
Marc_matNumber !< array of material numbers for hypoelastic material (Marc only)
#elif defined(Abaqus)
logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
#endif
public :: & public :: &
mesh_init, & mesh_init, &
#if defined(Marc4DAMASK) || defined(Abaqus)
mesh_FEasCP, &
#endif
mesh_build_cellnodes, & mesh_build_cellnodes, &
mesh_build_ipVolumes, & mesh_build_ipVolumes, &
mesh_build_ipCoordinates, & mesh_build_ipCoordinates, &
mesh_cellCenterCoordinates, & mesh_cellCenterCoordinates, &
mesh_get_Ncellnodes, & mesh_get_Ncellnodes, &
mesh_get_unitlength, & mesh_get_unitlength, &
mesh_get_nodeAtIP mesh_get_nodeAtIP, &
#ifdef Spectral #if defined(Spectral)
public :: &
mesh_spectral_getGrid, & mesh_spectral_getGrid, &
mesh_spectral_getSize mesh_spectral_getSize
#elif defined(Marc4DAMASK) || defined(Abaqus)
mesh_FEasCP
#endif #endif
private :: & private :: &
#ifdef Spectral mesh_get_damaskOptions, &
mesh_build_cellconnectivity, &
mesh_build_ipAreas, &
mesh_tell_statistics, &
FE_mapElemtype, &
mesh_faceMatch, &
mesh_build_FEdata, &
#if defined(Spectral)
mesh_spectral_getHomogenization, & mesh_spectral_getHomogenization, &
mesh_spectral_count, & mesh_spectral_count, &
mesh_spectral_count_cpSizes, & mesh_spectral_count_cpSizes, &
mesh_spectral_build_nodes, & mesh_spectral_build_nodes, &
mesh_spectral_build_elements, & mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, & mesh_spectral_build_ipNeighborhood
#elif defined Marc4DAMASK #elif defined(Marc4DAMASK) || defined(Abaqus)
mesh_build_nodeTwins, &
mesh_build_sharedElems, &
mesh_build_ipNeighborhood, &
#endif
#if defined(Marc4DAMASK)
mesh_marc_get_fileFormat, & mesh_marc_get_fileFormat, &
mesh_marc_get_tableStyles, & mesh_marc_get_tableStyles, &
mesh_marc_get_matNumber, & mesh_marc_get_matNumber, &
@ -440,8 +430,8 @@ module mesh
mesh_marc_map_nodes, & mesh_marc_map_nodes, &
mesh_marc_build_nodes, & mesh_marc_build_nodes, &
mesh_marc_count_cpSizes, & mesh_marc_count_cpSizes, &
mesh_marc_build_elements, & mesh_marc_build_elements
#elif defined Abaqus #elif defined(Abaqus)
mesh_abaqus_count_nodesAndElements, & mesh_abaqus_count_nodesAndElements, &
mesh_abaqus_count_elementSets, & mesh_abaqus_count_elementSets, &
mesh_abaqus_count_materials, & mesh_abaqus_count_materials, &
@ -452,20 +442,8 @@ module mesh
mesh_abaqus_map_nodes, & mesh_abaqus_map_nodes, &
mesh_abaqus_build_nodes, & mesh_abaqus_build_nodes, &
mesh_abaqus_count_cpSizes, & mesh_abaqus_count_cpSizes, &
mesh_abaqus_build_elements, & mesh_abaqus_build_elements
#endif #endif
#ifndef Spectral
mesh_build_nodeTwins, &
mesh_build_sharedElems, &
mesh_build_ipNeighborhood, &
#endif
mesh_get_damaskOptions, &
mesh_build_cellconnectivity, &
mesh_build_ipAreas, &
mesh_tell_statistics, &
FE_mapElemtype, &
mesh_faceMatch, &
mesh_build_FEdata
contains contains
@ -509,12 +487,12 @@ subroutine mesh_init(ip,el)
numerics_unitlength, & numerics_unitlength, &
worldrank worldrank
use FEsolving, only: & use FEsolving, only: &
FEsolving_execElem, &
#ifndef Spectral #ifndef Spectral
modelName, & modelName, &
calcMode, &
#endif #endif
FEsolving_execIP, & FEsolving_execElem, &
calcMode FEsolving_execIP
implicit none implicit none
#ifdef Spectral #ifdef Spectral
@ -523,7 +501,7 @@ subroutine mesh_init(ip,el)
integer :: ierr, worldsize integer :: ierr, worldsize
#endif #endif
integer(pInt), parameter :: FILEUNIT = 222_pInt integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip integer(pInt), intent(in), optional :: el, ip
integer(pInt) :: j integer(pInt) :: j
logical :: myDebug logical :: myDebug
@ -546,8 +524,12 @@ subroutine mesh_init(ip,el)
if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)') if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)')
geomSize = mesh_spectral_getSize(fileUnit) geomSize = mesh_spectral_getSize(fileUnit)
devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T),int(grid(2),C_INTPTR_T),& devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
int(grid(1),C_INTPTR_T)/2+1,PETSC_COMM_WORLD,local_K,local_K_offset) int(grid(2),C_INTPTR_T), &
int(grid(1),C_INTPTR_T)/2+1, &
PETSC_COMM_WORLD, &
local_K, & ! domain grid size along z
local_K_offset) ! domain grid offset along z
grid3 = int(local_K,pInt) grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt) grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
@ -647,25 +629,36 @@ subroutine mesh_init(ip,el)
call mesh_tell_statistics call mesh_tell_statistics
endif endif
#if defined(Marc4DAMASK) || defined(Abaqus)
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) & if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
#endif
if (debug_e < 1 .or. debug_e > mesh_NcpElems) & if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist call IO_error(602_pInt,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) & if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP
FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements FEsolving_execElem = [ 1_pInt,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP... allocate(FEsolving_execIP(2_pInt,mesh_NcpElems), source=1_pInt) ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
#if defined(Marc4DAMASK) || defined(Abaqus)
allocate(calcMode(mesh_maxNips,mesh_NcpElems)) allocate(calcMode(mesh_maxNips,mesh_NcpElems))
calcMode = .false. ! pretend to have collected what first call is asking (F = I) calcMode = .false. ! pretend to have collected what first call is asking (F = I)
#if defined(Marc4DAMASK) || defined(Abaqus)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
#else
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc"
#endif #endif
!!!! COMPATIBILITY HACK !!!!
! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
! hence, xxPerElem instead of maxXX
mesh_NipsPerElem = mesh_maxNips
mesh_NcellnodesPerElem = mesh_maxNcellnodes
! better name
mesh_homogenizationAt = mesh_element(3,:)
mesh_microstructureAt = mesh_element(4,:)
mesh_CPnodeID = mesh_element(5:4+mesh_NipsPerElem,:)
!!!!!!!!!!!!!!!!!!!!!!!!
end subroutine mesh_init end subroutine mesh_init
@ -1184,8 +1177,7 @@ subroutine mesh_spectral_count()
implicit none implicit none
mesh_Nelems = product(grid(1:2))*grid3 mesh_NcpElems= product(grid(1:2))*grid3
mesh_NcpElems= mesh_Nelems
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
mesh_NcpElemsGlobal = product(grid) mesh_NcpElemsGlobal = product(grid)
@ -1195,7 +1187,7 @@ end subroutine mesh_spectral_count
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements. !> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
!! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors', !! Sets global values 'mesh_maxNips', 'mesh_maxNipNeighbors',
!! and 'mesh_maxNcellnodes' !! and 'mesh_maxNcellnodes'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_count_cpSizes subroutine mesh_spectral_count_cpSizes
@ -1207,7 +1199,6 @@ subroutine mesh_spectral_count_cpSizes
g = FE_geomtype(t) g = FE_geomtype(t)
c = FE_celltype(g) c = FE_celltype(g)
mesh_maxNnodes = FE_Nnodes(t)
mesh_maxNips = FE_Nips(g) mesh_maxNips = FE_Nips(g)
mesh_maxNipNeighbors = FE_NipNeighbors(c) mesh_maxNipNeighbors = FE_NipNeighbors(c)
mesh_maxNcellnodes = FE_Ncellnodes(g) mesh_maxNcellnodes = FE_Ncellnodes(g)
@ -1268,13 +1259,13 @@ subroutine mesh_spectral_build_elements(fileUnit)
integer(pInt) :: & integer(pInt) :: &
e, i, & e, i, &
headerLength = 0_pInt, & headerLength = 0_pInt, &
maxIntCount, & maxDataPerLine, &
homog, & homog, &
elemType, & elemType, &
elemOffset elemOffset
integer(pInt), dimension(:), allocatable :: & integer(pInt), dimension(:), allocatable :: &
microstructures, & microstructures, &
mesh_microGlobal microGlobal
integer(pInt), dimension(1,1) :: & integer(pInt), dimension(1,1) :: &
dummySet = 0_pInt dummySet = 0_pInt
character(len=65536) :: & character(len=65536) :: &
@ -1304,16 +1295,16 @@ subroutine mesh_spectral_build_elements(fileUnit)
read(fileUnit,'(a65536)') line read(fileUnit,'(a65536)') line
enddo enddo
maxIntCount = 0_pInt maxDataPerLine = 0_pInt
i = 1_pInt i = 1_pInt
do while (i > 0_pInt) do while (i > 0_pInt)
i = IO_countContinuousIntValues(fileUnit) i = IO_countContinuousIntValues(fileUnit)
maxIntCount = max(maxIntCount, i) maxDataPerLine = max(maxDataPerLine, i) ! found a longer line?
enddo enddo
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems), source = 0_pInt) allocate(mesh_element (4_pInt+8_pInt,mesh_NcpElems), source = 0_pInt)
allocate (microstructures (1_pInt+maxIntCount), source = 1_pInt) allocate(microstructures (1_pInt+maxDataPerLine), source = 1_pInt) ! prepare to receive counter and max data size
allocate (mesh_microGlobal(mesh_NcpElemsGlobal), source = 1_pInt) allocate(microGlobal (mesh_NcpElemsGlobal), source = 1_pInt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read in microstructures ! read in microstructures
@ -1324,10 +1315,10 @@ subroutine mesh_spectral_build_elements(fileUnit)
e = 0_pInt e = 0_pInt
do while (e < mesh_NcpElemsGlobal .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) do while (e < mesh_NcpElemsGlobal .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 microstructures = IO_continuousIntValues(fileUnit,maxDataPerLine,dummyName,dummySet,0_pInt) ! get affected elements
do i = 1_pInt,microstructures(1_pInt) do i = 1_pInt,microstructures(1_pInt)
e = e+1_pInt ! valid element entry e = e+1_pInt ! valid element entry
mesh_microGlobal(e) = microstructures(1_pInt+i) microGlobal(e) = microstructures(1_pInt+i)
enddo enddo
enddo enddo
@ -1336,10 +1327,10 @@ subroutine mesh_spectral_build_elements(fileUnit)
e = 0_pInt e = 0_pInt
do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!) do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!)
e = e+1_pInt ! valid element entry e = e+1_pInt ! valid element entry
mesh_element( 1,e) = e ! FE id mesh_element( 1,e) = -1_pInt ! DEPRECATED
mesh_element( 2,e) = elemType ! elem type mesh_element( 2,e) = elemType ! elem type
mesh_element( 3,e) = homog ! homogenization mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = mesh_microGlobal(e+elemOffset) ! microstructure mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &
((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node ((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt mesh_element( 6,e) = mesh_element(5,e) + 1_pInt
@ -1715,8 +1706,8 @@ subroutine mesh_marc_map_elementSets(fileUnit)
character(len=300) :: line character(len=300) :: line
integer(pInt) :: elemSet = 0_pInt integer(pInt) :: elemSet = 0_pInt
allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' allocate (mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = ''
allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets), source=0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -1813,7 +1804,7 @@ subroutine mesh_marc_map_elements(fileUnit)
integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts
integer(pInt) :: i,cpElem = 0_pInt integer(pInt) :: i,cpElem = 0_pInt
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems), source = 0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -1883,7 +1874,7 @@ subroutine mesh_marc_map_nodes(fileUnit)
integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt) :: i integer(pInt) :: i
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes),source=0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -1930,8 +1921,8 @@ subroutine mesh_marc_build_nodes(fileUnit)
character(len=300) :: line character(len=300) :: line
integer(pInt) :: i,j,m integer(pInt) :: i,j,m
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal allocate ( mesh_node0 (3,mesh_Nnodes), source=0.0_pReal)
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal allocate ( mesh_node (3,mesh_Nnodes), source=0.0_pReal)
610 FORMAT(A300) 610 FORMAT(A300)
@ -2023,7 +2014,8 @@ subroutine mesh_marc_build_elements(fileUnit)
IO_skipChunks, & IO_skipChunks, &
IO_stringPos, & IO_stringPos, &
IO_intValue, & IO_intValue, &
IO_continuousIntValues IO_continuousIntValues, &
IO_error
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -2034,7 +2026,8 @@ subroutine mesh_marc_build_elements(fileUnit)
integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts integer(pInt), dimension(1_pInt+mesh_NcpElems) :: contInts
integer(pInt) :: i,j,t,sv,myVal,e,nNodesAlreadyRead integer(pInt) :: i,j,t,sv,myVal,e,nNodesAlreadyRead
allocate (mesh_element(4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt allocate(mesh_element(4_pInt+mesh_maxNnodes,mesh_NcpElems), source=0_pInt)
mesh_elemType = -1_pInt
610 FORMAT(A300) 610 FORMAT(A300)
@ -2049,8 +2042,11 @@ subroutine mesh_marc_build_elements(fileUnit)
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1_pInt)) e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1_pInt))
if (e /= 0_pInt) then ! disregard non CP elems if (e /= 0_pInt) then ! disregard non CP elems
mesh_element(1,e) = IO_IntValue (line,chunkPos,1_pInt) ! FE id mesh_element(1,e) = -1_pInt ! DEPRECATED
t = FE_mapElemtype(IO_StringValue(line,chunkPos,2_pInt)) ! elem type t = FE_mapElemtype(IO_StringValue(line,chunkPos,2_pInt)) ! elem type
if (mesh_elemType /= t .and. mesh_elemType /= -1_pInt) &
call IO_error(191,el=t,ip=mesh_elemType)
mesh_elemType = t
mesh_element(2,e) = t mesh_element(2,e) = t
nNodesAlreadyRead = 0_pInt nNodesAlreadyRead = 0_pInt
do j = 1_pInt,chunkPos(1)-2_pInt do j = 1_pInt,chunkPos(1)-2_pInt
@ -2280,8 +2276,8 @@ subroutine mesh_abaqus_map_elementSets(fileUnit)
integer(pInt) :: elemSet = 0_pInt,i integer(pInt) :: elemSet = 0_pInt,i
logical :: inPart = .false. logical :: inPart = .false.
allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = '' allocate (mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = ''
allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets)) ; mesh_mapElemSet = 0_pInt allocate (mesh_mapElemSet(1_pInt+mesh_maxNelemInSet,mesh_NelemSets),source=0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -2332,8 +2328,8 @@ subroutine mesh_abaqus_map_materials(fileUnit)
logical :: inPart = .false. logical :: inPart = .false.
character(len=64) :: elemSetName,materialName character(len=64) :: elemSetName,materialName
allocate (mesh_nameMaterial(mesh_Nmaterials)) ; mesh_nameMaterial = '' allocate (mesh_nameMaterial(mesh_Nmaterials)); mesh_nameMaterial = ''
allocate (mesh_mapMaterial(mesh_Nmaterials)) ; mesh_mapMaterial = '' allocate (mesh_mapMaterial(mesh_Nmaterials)); mesh_mapMaterial = ''
610 FORMAT(A300) 610 FORMAT(A300)
@ -2450,7 +2446,7 @@ subroutine mesh_abaqus_map_elements(fileUnit)
logical :: materialFound = .false. logical :: materialFound = .false.
character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS? character (len=64) materialName,elemSetName ! why limited to 64? ABAQUS?
allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt allocate (mesh_mapFEtoCPelem(2,mesh_NcpElems), source = 0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -2513,7 +2509,7 @@ subroutine mesh_abaqus_map_nodes(fileUnit)
integer(pInt) :: i,c,cpNode = 0_pInt integer(pInt) :: i,c,cpNode = 0_pInt
logical :: inPart = .false. logical :: inPart = .false.
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source=0_pInt)
610 FORMAT(A300) 610 FORMAT(A300)
@ -2575,8 +2571,8 @@ subroutine mesh_abaqus_build_nodes(fileUnit)
integer(pInt) :: i,j,m,c integer(pInt) :: i,j,m,c
logical :: inPart logical :: inPart
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal allocate ( mesh_node0 (3,mesh_Nnodes), source=0.0_pReal)
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal allocate ( mesh_node (3,mesh_Nnodes), source=0.0_pReal)
610 FORMAT(A300) 610 FORMAT(A300)
@ -2688,8 +2684,8 @@ subroutine mesh_abaqus_build_elements(fileUnit)
IO_intValue, & IO_intValue, &
IO_extractValue, & IO_extractValue, &
IO_floatValue, & IO_floatValue, &
IO_error, & IO_countDataLines, &
IO_countDataLines IO_error
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -2701,7 +2697,8 @@ subroutine mesh_abaqus_build_elements(fileUnit)
character (len=64) :: materialName,elemSetName character (len=64) :: materialName,elemSetName
character(len=300) :: line character(len=300) :: line
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt allocate(mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems), source=0_pInt)
mesh_elemType = -1_pInt
610 FORMAT(A300) 610 FORMAT(A300)
@ -2730,7 +2727,10 @@ subroutine mesh_abaqus_build_elements(fileUnit)
chunkPos = IO_stringPos(line) ! limit to 64 nodes max chunkPos = IO_stringPos(line) ! limit to 64 nodes max
e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1_pInt)) e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1_pInt))
if (e /= 0_pInt) then ! disregard non CP elems if (e /= 0_pInt) then ! disregard non CP elems
mesh_element(1,e) = IO_intValue(line,chunkPos,1_pInt) ! FE id mesh_element(1,e) = -1_pInt ! DEPRECATED
if (mesh_elemType /= t .and. mesh_elemType /= -1_pInt) &
call IO_error(191,el=t,ip=mesh_elemType)
mesh_elemType = t
mesh_element(2,e) = t ! elem type mesh_element(2,e) = t ! elem type
nNodesAlreadyRead = 0_pInt nNodesAlreadyRead = 0_pInt
do j = 1_pInt,chunkPos(1)-1_pInt do j = 1_pInt,chunkPos(1)-1_pInt
@ -3010,7 +3010,7 @@ subroutine mesh_build_sharedElems
myDim, & ! dimension index myDim, & ! dimension index
nodeTwin ! node twin in the specified dimension nodeTwin ! node twin in the specified dimension
integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt), dimension (:), allocatable :: node_seen integer(pInt), dimension(:), allocatable :: node_seen
allocate(node_seen(maxval(FE_NmatchingNodes))) allocate(node_seen(maxval(FE_NmatchingNodes)))
@ -3035,8 +3035,7 @@ subroutine mesh_build_sharedElems
mesh_maxNsharedElems = int(maxval(node_count),pInt) ! most shared node mesh_maxNsharedElems = int(maxval(node_count),pInt) ! most shared node
allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes)) allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes),source=0_pInt)
mesh_sharedElem = 0_pInt
do e = 1_pInt,mesh_NcpElems do e = 1_pInt,mesh_NcpElems
g = FE_geomtype(mesh_element(2,e)) ! get elemGeomType g = FE_geomtype(mesh_element(2,e)) ! get elemGeomType
@ -3258,7 +3257,7 @@ subroutine mesh_tell_statistics
if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified
if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified
allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2)),source = 0_pInt)
do e = 1_pInt,mesh_NcpElems do e = 1_pInt,mesh_NcpElems
if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,el=e) ! no homogenization specified if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,el=e) ! no homogenization specified
if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,el=e) ! no microstructure specified if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,el=e) ! no microstructure specified
@ -3268,13 +3267,8 @@ subroutine mesh_tell_statistics
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
write(6,'(/,a,/)') ' Input Parser: STATISTICS' write(6,'(/,a,/)') ' Input Parser: STATISTICS'
write(6,*) mesh_Nelems, ' : total number of elements in mesh'
write(6,*) mesh_NcpElems, ' : total number of CP elements in mesh' write(6,*) mesh_NcpElems, ' : total number of CP elements in mesh'
write(6,*) mesh_Nnodes, ' : total number of nodes in mesh' write(6,*) mesh_Nnodes, ' : total number of nodes in mesh'
write(6,*) mesh_maxNnodes, ' : max number of nodes in any CP element'
write(6,*) mesh_maxNips, ' : max number of IPs in any CP element'
write(6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element'
write(6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node'
write(6,'(/,a,/)') ' Input Parser: HOMOGENIZATION/MICROSTRUCTURE' write(6,'(/,a,/)') ' Input Parser: HOMOGENIZATION/MICROSTRUCTURE'
write(6,*) mesh_maxValStateVar(1), ' : maximum homogenization index' write(6,*) mesh_maxValStateVar(1), ' : maximum homogenization index'
write(6,*) mesh_maxValStateVar(2), ' : maximum microstructure index' write(6,*) mesh_maxValStateVar(2), ' : maximum microstructure index'
@ -3527,11 +3521,11 @@ subroutine mesh_build_FEdata
implicit none implicit none
integer(pInt) :: me integer(pInt) :: me
allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes)); FE_nodesAtIP = 0_pInt allocate(FE_nodesAtIP(FE_maxmaxNnodesAtIP,FE_maxNips,FE_Ngeomtypes), source=0_pInt)
allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes)); FE_ipNeighbor = 0_pInt allocate(FE_ipNeighbor(FE_maxNipNeighbors,FE_maxNips,FE_Ngeomtypes), source=0_pInt)
allocate(FE_cell(FE_maxNcellnodesPerCell,FE_maxNips,FE_Ngeomtypes)); FE_cell = 0_pInt allocate(FE_cell(FE_maxNcellnodesPerCell,FE_maxNips,FE_Ngeomtypes), source=0_pInt)
allocate(FE_cellnodeParentnodeWeights(FE_maxNnodes,FE_maxNcellnodes,FE_Nelemtypes)); FE_cellnodeParentnodeWeights = 0.0_pReal allocate(FE_cellnodeParentnodeWeights(FE_maxNnodes,FE_maxNcellnodes,FE_Nelemtypes), source=0.0_pReal)
allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes)); FE_cellface = 0_pInt allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0_pInt)
!*** fill FE_nodesAtIP with data *** !*** fill FE_nodesAtIP with data ***

66
src/meshFEM.f90
View File

@ -19,21 +19,27 @@ use PETScis
implicit none implicit none
private private
integer(pInt), public, parameter :: &
mesh_ElemType=1_pInt !< Element type of the mesh (only support homogeneous meshes)
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
mesh_Nboundaries, & mesh_Nboundaries, &
mesh_NcpElems, & !< total number of CP elements in mesh mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElemsGlobal, & mesh_NcpElemsGlobal, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
mesh_maxNnodes, & !< max number of nodes in any CP element mesh_NipsPerElem, & !< number of IPs in per element
mesh_maxNips, & !< max number of IPs in any CP element mesh_maxNipNeighbors
mesh_maxNipNeighbors, & !!!! BEGIN DEPRECATED !!!!!
mesh_Nelems !< total number of elements in mesh integer(pInt), public, protected :: &
mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!!
real(pReal), public, protected :: charLength integer(pInt), dimension(:), allocatable, public, protected :: &
mesh_homogenizationAt, & !< homogenization ID of each element
mesh_microstructureAt !< microstructure ID of each element
integer(pInt), dimension(:,:), allocatable, public, protected :: & integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element !< FEid, type(internal representation), material, texture, node indices as CP IDs mesh_element !DEPRECATED
real(pReal), dimension(:,:), allocatable, public :: & real(pReal), dimension(:,:), allocatable, public :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
@ -62,26 +68,16 @@ use PETScis
mesh_boundaries mesh_boundaries
integer(pInt), parameter, public :: & integer(pInt), dimension(1_pInt), parameter, public :: FE_geomtype = & !< geometry type of particular element type
FE_Nelemtypes = 1_pInt, &
FE_Ngeomtypes = 1_pInt, &
FE_Ncelltypes = 1_pInt, &
FE_maxNnodes = 1_pInt, &
FE_maxNips = 14_pInt
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = & !< geometry type of particular element type
int([1],pInt) int([1],pInt)
integer(pInt), dimension(FE_Ngeomtypes), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type integer(pInt), dimension(1_pInt), parameter, public :: FE_celltype = & !< cell type that is used by each geometry type
int([1],pInt) int([1],pInt)
integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = & !< number of nodes that constitute a specific type of element integer(pInt), dimension(1_pInt), public :: FE_Nips = & !< number of IPs in a specific type of element
int([0],pInt) int([0],pInt)
integer(pInt), dimension(FE_Ngeomtypes), public :: FE_Nips = & !< number of IPs in a specific type of element integer(pInt), dimension(1_pInt), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([0],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([6],pInt) int([6],pInt)
@ -98,7 +94,7 @@ contains
!> @brief initializes the mesh by calling all necessary private routines the mesh module !> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver !! Order and routines strongly depend on type of solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el) subroutine mesh_init()
use DAMASK_interface use DAMASK_interface
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use IO, only: & use IO, only: &
@ -120,15 +116,13 @@ subroutine mesh_init(ip,el)
worldsize worldsize
use FEsolving, only: & use FEsolving, only: &
FEsolving_execElem, & FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP
calcMode
use FEM_Zoo, only: & use FEM_Zoo, only: &
FEM_Zoo_nQuadrature, & FEM_Zoo_nQuadrature, &
FEM_Zoo_QuadraturePoints FEM_Zoo_QuadraturePoints
implicit none implicit none
integer(pInt), parameter :: FILEUNIT = 222_pInt integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip
integer(pInt) :: j integer(pInt) :: j
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer :: dimPlex integer :: dimPlex
@ -212,29 +206,25 @@ subroutine mesh_init(ip,el)
endif endif
call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_Nelems,ierr) call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
mesh_NcpElems = mesh_Nelems
FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNnodes = FE_Nnodes(1_pInt)
mesh_maxNips = FE_Nips(1_pInt) mesh_maxNips = FE_Nips(1_pInt)
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex) call mesh_FEM_build_ipVolumes(dimPlex)
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt allocate (mesh_element (4_pInt,mesh_NcpElems)); mesh_element = 0_pInt
do j = 1, mesh_NcpElems do j = 1, mesh_NcpElems
mesh_element( 1,j) = j mesh_element( 1,j) = -1_pInt ! DEPRECATED
mesh_element( 2,j) = 1_pInt ! elem type mesh_element( 2,j) = mesh_elemType ! elem type
mesh_element( 3,j) = 1_pInt ! homogenization mesh_element( 3,j) = 1_pInt ! homogenization
call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
end do end do
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
if (debug_e < 1 .or. debug_e > mesh_NcpElems) & if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
call IO_error(602_pInt,ext_msg='element') ! selected element does not exist call IO_error(602_pInt,ext_msg='element') ! selected element does not exist
if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) & if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
@ -245,10 +235,14 @@ subroutine mesh_init(ip,el)
allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP... allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt ! parallel loop bounds set to comprise from first IP...
forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j))) ! ...up to own IP count for each element
if (allocated(calcMode)) deallocate(calcMode) !!!! COMPATIBILITY HACK !!!!
allocate(calcMode(mesh_maxNips,mesh_NcpElems)) ! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
calcMode = .false. ! pretend to have collected what first call is asking (F = I) ! hence, xxPerElem instead of maxXX
calcMode(ip,el) = .true. ! first ip,el needs to be already pingponged to "calc" mesh_NipsPerElem = mesh_maxNips
! better name
mesh_homogenizationAt = mesh_element(3,:)
mesh_microstructureAt = mesh_element(4,:)
!!!!!!!!!!!!!!!!!!!!!!!!
end subroutine mesh_init end subroutine mesh_init