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

switched element library to geomType based. 

saves to copy same geometry description for different elements that are essentially similar regarding the IP number but differ in total node count.

introduced quadratic tetrahedron (Marc element 127 -- element 157 might also work, but did not perform well in fully elastic calc so far)
This commit is contained in:
Philip Eisenlohr 2012-11-15 22:45:20 +00:00
parent 36b890fba3
commit d9a98417ca
9 changed files with 729 additions and 660 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
code/CPFEM.f90
View File

@ -277,7 +277,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
mesh_maxNips, &
mesh_element, &
FE_Nips, &
FE_Nnodes
FE_Nnodes, &
FE_geomtype
use material, only: homogenization_maxNgrains, &
microstructure_elemhomo, &
material_phase
@ -532,7 +533,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
forall (i = 2:FE_Nips(FE_geomtype(mesh_element(2,e)))) ! copy results of first IP to all others
materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e)
materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e)
materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e)

5
code/DAMASK_marc.f90
View File

@ -242,7 +242,8 @@ subroutine hypela2(&
mesh_node, &
mesh_build_subNodeCoords, &
mesh_build_ipCoordinates, &
FE_Nnodes
FE_Nnodes, &
FE_geomtype
use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done
!$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS
@ -361,7 +362,7 @@ subroutine hypela2(&
else
computationMode = 3 ! plain collect
endif
do node = 1,FE_Nnodes(mesh_element(2,cp_en))
do node = 1,FE_Nnodes(FE_geomtype(mesh_element(2,cp_en)))
FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en))
mesh_node(1:3,FEnodeID) = mesh_node0(1:3,FEnodeID) + numerics_unitlength * dispt(1:3,node)
enddo

11
code/constitutive.f90
View File

@ -93,9 +93,10 @@ subroutine constitutive_init
use mesh, only: mesh_maxNips, &
mesh_NcpElems, &
mesh_element, &
mesh_ipNeighborhood, &
FE_Nips, &
FE_NipNeighbors, &
mesh_ipNeighborhood
FE_geomtype
use material, only: material_phase, &
material_Nphase, &
material_localFileExt, &
@ -221,7 +222,7 @@ endif
do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs
do g = 1_pInt,myNgrains ! loop over grains
select case(phase_elasticity(material_phase(g,i,e)))
@ -375,8 +376,8 @@ do e = 1_pInt,mesh_NcpElems ! loop over element
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(myInstance)
case (constitutive_nonlocal_label)
select case(mesh_element(2,e))
case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt)
select case(FE_geomtype(mesh_element(2,e)))
case (7_pInt,8_pInt,9_pInt,10_pInt)
! all fine
case default
call IO_error(253_pInt,e,i,g)
@ -417,7 +418,7 @@ enddo
call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax))
do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall(i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:myNgrains)
forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains)
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init
endforall

47
code/constitutive_nonlocal.f90
View File

@ -952,8 +952,9 @@ use math, only: math_sampleGaussVar
use mesh, only: mesh_ipVolume, &
mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
FE_Nips, &
mesh_element
FE_geomtype
use material, only: material_phase, &
phase_plasticityInstance, &
phase_plasticity
@ -1012,7 +1013,7 @@ do myInstance = 1_pInt,maxNinstance
minimumIpVolume = 1e99_pReal
do el = 1_pInt,mesh_NcpElems
do ip = 1_pInt,FE_Nips(mesh_element(2,el))
do ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,el)))
if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) &
.and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then
totalVolume = totalVolume + mesh_ipVolume(ip,el)
@ -1029,7 +1030,7 @@ do myInstance = 1_pInt,maxNinstance
do while(meanDensity < constitutive_nonlocal_rhoSglRandom(myInstance))
call random_number(rnd)
el = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt)
ip = nint(rnd(2)*real(FE_Nips(mesh_element(2,el)),pReal)+0.5_pReal,pInt)
ip = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,el))),pReal)+0.5_pReal,pInt)
if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) &
.and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then
s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt)
@ -1042,7 +1043,7 @@ do myInstance = 1_pInt,maxNinstance
! homogeneous distribution of density with some noise
else
do el = 1_pInt,mesh_NcpElems
do ip = 1_pInt,FE_Nips(mesh_element(2,el))
do ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,el)))
if (constitutive_nonlocal_label == phase_plasticity(material_phase(1,ip,el)) &
.and. myInstance == phase_plasticityInstance(material_phase(1,ip,el))) then
do f = 1_pInt,lattice_maxNslipFamily
@ -1177,12 +1178,13 @@ use debug, only: debug_level, &
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
FE_NipNeighbors, &
FE_maxNipNeighbors, &
mesh_ipNeighborhood, &
mesh_ipCoordinates, &
mesh_ipVolume, &
mesh_ipAreaNormal
mesh_ipAreaNormal, &
FE_NipNeighbors, &
FE_maxNipNeighbors, &
FE_geomtype
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_localPlasticity, &
@ -1347,7 +1349,7 @@ if (.not. phase_localPlasticity(phase) .and. constitutive_nonlocal_shortRangeStr
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
do n = 1_pInt,FE_NipNeighbors(mesh_element(2,el))
do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,el)))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if (neighboring_el > 0 .and. neighboring_ip > 0) then
@ -2031,11 +2033,12 @@ use math, only: math_norm3, &
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
FE_NipNeighbors, &
mesh_ipNeighborhood, &
mesh_ipVolume, &
mesh_ipArea, &
mesh_ipAreaNormal
mesh_ipAreaNormal, &
FE_NipNeighbors, &
FE_geomtype
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance, &
@ -2300,7 +2303,7 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then
my_Fe = Fe(1:3,1:3,g,ip,el)
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el))
do n = 1_pInt,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,el))) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
neighboring_n = mesh_ipNeighborhood(3,n,ip,el)
@ -2561,9 +2564,10 @@ use material, only: material_phase, &
homogenization_maxNgrains
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors, &
mesh_maxNips, &
mesh_NcpElems
mesh_NcpElems, &
FE_NipNeighbors, &
FE_geomtype
use lattice, only: lattice_sn, &
lattice_sd
@ -2594,7 +2598,7 @@ integer(pInt) Nneighbors, & !
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
FE_NipNeighbors(mesh_element(2,e))) :: &
FE_NipNeighbors(FE_geomtype(mesh_element(2,e)))) :: &
compatibility ! compatibility for current element and ip
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &
slipNormal, &
@ -2606,10 +2610,10 @@ logical, dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(mat
belowThreshold
Nneighbors = FE_NipNeighbors(mesh_element(2,e))
my_phase = material_phase(1,i,e)
my_texture = material_texture(1,i,e)
my_instance = phase_plasticityInstance(my_phase)
Nneighbors = FE_NipNeighbors(FE_geomtype(mesh_element(2,e)))
my_phase = material_phase(1,i,e)
my_texture = material_texture(1,i,e)
my_instance = phase_plasticityInstance(my_phase)
my_structure = constitutive_nonlocal_structure(my_instance)
ns = constitutive_nonlocal_totalNslip(my_instance)
slipNormal(1:3,1:ns) = lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
@ -2766,10 +2770,11 @@ use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
mesh_node0, &
FE_Nips, &
mesh_cellCenterCoordinates, &
mesh_ipVolume, &
mesh_periodicSurface
mesh_periodicSurface, &
FE_Nips, &
FE_geomtype
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_localPlasticity, &
@ -2890,7 +2895,7 @@ if (.not. phase_localPlasticity(phase)) then
!* but only consider nonlocal neighbors within a certain cutoff radius R
do neighboring_el = 1_pInt,mesh_NcpElems
ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el))
ipLoop: do neighboring_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighboring_el)))
neighboring_phase = material_phase(g,neighboring_ip,neighboring_el)
if (phase_localPlasticity(neighboring_phase)) then
cycle

9
code/crystallite.f90
View File

@ -3038,7 +3038,8 @@ use material, only: material_phase, &
phase_plasticityInstance
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors
FE_NipNeighbors, &
FE_geomtype
use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility
@ -3104,19 +3105,19 @@ logical error
! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,e))) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity
neighboringInstance = phase_plasticityInstance(neighboringPhase)
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = &
math_QuaternionDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
crystallite_symmetryID(1,i,e)) ! calculate disorientation
crystallite_symmetryID(1,i,e)) ! calculate disorientation
else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis
endif

4
code/homogenization.f90
View File

@ -79,7 +79,7 @@ subroutine homogenization_init(Temperature)
use debug, only: debug_level, debug_homogenization, debug_levelBasic
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile, &
IO_write_jobBinaryIntFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips,FE_geomtype
use material
use constitutive, only: constitutive_maxSizePostResults
use crystallite, only: crystallite_maxSizePostResults
@ -178,7 +178,7 @@ subroutine homogenization_init(Temperature)
!$OMP PARALLEL DO PRIVATE(myInstance)
do e = 1,mesh_NcpElems ! loop over elements
myInstance = homogenization_typeInstance(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs
select case(homogenization_type(mesh_element(3,e)))
case (homogenization_isostrain_label)
if (homogenization_isostrain_sizeState(myInstance) > 0_pInt) then

6
code/homogenization_RGC.f90
View File

@ -77,7 +77,7 @@ subroutine homogenization_RGC_init(&
math_sampleRandomOri,&
math_EulerToR,&
INRAD
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips,FE_geomtype
use IO
use material
@ -164,7 +164,7 @@ subroutine homogenization_RGC_init(&
myInstance = homogenization_typeInstance(mesh_element(3,e))
if (all (homogenization_RGC_angles(:,myInstance) >= 399.9_pReal)) then
homogenization_RGC_orientation(:,:,1,e) = math_EulerToR(math_sampleRandomOri())
do i = 1_pInt,FE_Nips(mesh_element(2,e))
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
if (microstructure_elemhomo(mesh_element(4,e))) then
homogenization_RGC_orientation(:,:,i,e) = homogenization_RGC_orientation(:,:,1,e)
else
@ -172,7 +172,7 @@ subroutine homogenization_RGC_init(&
endif
enddo
else
do i = 1_pInt,FE_Nips(mesh_element(2,e))
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(homogenization_RGC_angles(:,myInstance)*inRad)
enddo
endif

20
code/material.f90
View File

@ -672,7 +672,8 @@ subroutine material_populateGrains
mesh_maxNips, &
mesh_NcpElems, &
mesh_ipVolume, &
FE_Nips
FE_Nips, &
FE_geomtype
use IO, only: IO_error, &
IO_hybridIA
use FEsolving, only: FEsolving_execIP
@ -720,6 +721,7 @@ subroutine material_populateGrains
! identify maximum grain count per IP (from element) and find grains per homog/micro pair
do e = 1_pInt,mesh_NcpElems
t = FE_geomtype(mesh_element(2,e))
homog = mesh_element(3,e)
micro = mesh_element(4,e)
if (homog < 1_pInt .or. homog > material_Nhomogenization) & ! out of bounds
@ -729,7 +731,7 @@ subroutine material_populateGrains
if (microstructure_elemhomo(micro)) then
dGrains = homogenization_Ngrains(homog)
else
dGrains = homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
dGrains = homogenization_Ngrains(homog) * FE_Nips(t)
endif
Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains
Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
@ -767,15 +769,16 @@ subroutine material_populateGrains
grain = 0_pInt
do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(hme,homog,micro) ! 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))
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(mesh_element(2,e)),e))/&
volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(t),e))/&
real(dGrains,pReal)
grain = grain + dGrains ! wind forward by NgrainsPerIP
else
forall (i = 1_pInt:FE_Nips(mesh_element(2,e))) & ! loop over IPs
forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs
volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*NgrainsPerIP
endif
enddo
@ -886,8 +889,9 @@ subroutine material_populateGrains
grain = 0_pInt
do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(hme,homog,micro) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
t = FE_geomtype(mesh_element(2,e))
if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs
forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains
forall (i = 1_pInt:FE_Nips(t), g = 1_pInt:dGrains) ! loop over IPs and grains
material_volume(g,i,e) = volumeOfGrain(grain+g)
material_phase(g,i,e) = phaseOfGrain(grain+g)
material_texture(g,i,e) = textureOfGrain(grain+g)
@ -896,13 +900,13 @@ subroutine material_populateGrains
FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this
grain = grain + dGrains ! wind forward by NgrainsPerIP
else
forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains
forall (i = 1_pInt:FE_Nips(t), g = 1_pInt:dGrains) ! loop over IPs and grains
material_volume(g,i,e) = volumeOfGrain(grain+(i-1_pInt)*dGrains+g)
material_phase(g,i,e) = phaseOfGrain(grain+(i-1_pInt)*dGrains+g)
material_texture(g,i,e) = textureOfGrain(grain+(i-1_pInt)*dGrains+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1_pInt)*dGrains+g)
end forall
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*NgrainsPerIP
endif
enddo
endif ! active homog,micro pair

1282
code/mesh.f90
View File

File diff suppressed because it is too large Load Diff