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

removed coordinates from call to CPFEM_general and introduced direct storage of coordinates to mesh_ipCoordinates

This commit is contained in:
Martin Diehl 2012-11-14 14:38:10 +00:00
parent 1896c07c70
commit 7e683ca7c3
10 changed files with 295 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
code/CPFEM.f90
View File

@ -230,15 +230,14 @@ end subroutine CPFEM_init
!*** perform initialization at first call, update variables and *** !*** perform initialization at first call, update variables and ***
!*** call the actual material model *** !*** call the actual material model ***
!*********************************************************************** !***********************************************************************
subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,&
& jacobian, pstress, dPdF) & jacobian, pstress, dPdF)
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt use prec, only: pInt
use numerics, only: defgradTolerance, & use numerics, only: defgradTolerance, &
iJacoStiffness, & iJacoStiffness
numerics_unitlength
use debug, only: debug_level, & use debug, only: debug_level, &
debug_CPFEM, & debug_CPFEM, &
debug_levelBasic, & debug_levelBasic, &
@ -277,12 +276,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_element, & mesh_element, &
mesh_node0, &
mesh_node, &
mesh_ipCoordinates, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
FE_Nips, & FE_Nips, &
FE_Nnodes FE_Nnodes
use material, only: homogenization_maxNgrains, & use material, only: homogenization_maxNgrains, &
@ -324,9 +317,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
IP ! FE integration point number IP ! FE integration point number
real(pReal), intent(inout) :: Temperature ! temperature real(pReal), intent(inout) :: Temperature ! temperature
real(pReal), intent(in) :: dt ! time increment real(pReal), intent(in) :: dt ! time increment
real(pReal), dimension (3,*), intent(in) :: coords ! MARC: displacements for each node of the current element
! ABAQUS: coordinates of the current material point (IP)
! SPECTRAL: coordinates of the current material point (IP)
real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0
ffn1 ! deformation gradient for t=t1 ffn1 ! deformation gradient for t=t1
integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results
@ -356,9 +346,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
jacobian3333 ! jacobian in Matrix notation jacobian3333 ! jacobian in Matrix notation
integer(pInt) cp_en, & ! crystal plasticity element number integer(pInt) cp_en, & ! crystal plasticity element number
i, j, k, l, m, n, e i, j, k, l, m, n, e
#ifdef Marc
integer(pInt):: node, FEnodeID
#endif
logical updateJaco ! flag indicating if JAcobian has to be updated logical updateJaco ! flag indicating if JAcobian has to be updated
@ -535,11 +522,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
#ifdef Marc
! marc returns nodal coordinates, whereas Abaqus and spectral solver return ip coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates.
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
#endif
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) write(6,'(a,i8,a,i8)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
@ -611,14 +593,6 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_calc_done = .false. CPFEM_calc_done = .false.
#ifndef Marc
mesh_ipCoordinates(1:3,IP,cp_en) = numerics_unitlength * coords(1:3,1)
#else
do node = 1,FE_Nnodes(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 * coords(1:3,node)
enddo
#endif
! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+--

12
code/DAMASK_abaqus_exp.f
View File

@ -170,6 +170,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use numerics, only: numerics_unitlength
use FEsolving, only: cycleCounter, & use FEsolving, only: cycleCounter, &
theTime, & theTime, &
outdatedByNewInc, & outdatedByNewInc, &
@ -182,7 +183,8 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
debug_levelBasic, & debug_levelBasic, &
debug_level, & debug_level, &
debug_abaqus debug_abaqus
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP, &
mesh_ipCoordinates
use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
@ -211,7 +213,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
real(pReal), dimension(6) :: stress real(pReal), dimension(6) :: stress
real(pReal), dimension(6,6) :: ddsdde real(pReal), dimension(6,6) :: ddsdde
real(pReal) temp, timeInc real(pReal) temp, timeInc
integer(pInt) computationMode, n, i integer(pInt) computationMode, n, i, cp_en
logical :: cutBack logical :: cutBack
do n = 1,nblock ! loop over vector of IPs do n = 1,nblock ! loop over vector of IPs
@ -289,8 +291,10 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
defgrd1(3,2) = defgradNew(n,8) defgrd1(3,2) = defgradNew(n,8)
endif endif
coordinates = coordMp(n,1:3) cp_en = mesh_FEasCP('elem',nElement(n))
call CPFEM_general(computationMode,coordinates,defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF) mesh_ipCoordinates(1:3,n,cp_en) = numerics_unitlength * coordMp(n,1:3)
call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,cp_en,nMatPoint(n),stress,ddsdde, pstress, dPdF)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! straight: 11, 22, 33, 12, 23, 13 ! straight: 11, 22, 33, 12, 23, 13

7
code/DAMASK_abaqus_std.f
View File

@ -130,6 +130,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use numerics, only: numerics_unitlength
use FEsolving, only: cycleCounter, & use FEsolving, only: cycleCounter, &
theInc, & theInc, &
calcMode, & calcMode, &
@ -147,7 +148,8 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
debug_levelBasic, & debug_levelBasic, &
debug_level, & debug_level, &
debug_abaqus debug_abaqus
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP, &
mesh_ipCoordinates
use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
@ -253,6 +255,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
else else
computationMode = 3 ! plain collect computationMode = 3 ! plain collect
endif endif
mesh_ipCoordinates(1:3,npt,cp_en) = numerics_unitlength * COORDS
endif endif
theTime = time(2) ! record current starting time theTime = time(2) ! record current starting time
@ -266,7 +269,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call CPFEM_general(computationMode,COORDS,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! straight: 11, 22, 33, 12, 23, 13 ! straight: 11, 22, 33, 12, 23, 13

22
code/DAMASK_marc.f90
View File

@ -221,6 +221,7 @@ subroutine hypela2(&
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use numerics, only: numerics_unitlength
use FEsolving, only: cycleCounter, & use FEsolving, only: cycleCounter, &
theInc, & theInc, &
calcMode, & calcMode, &
@ -235,7 +236,13 @@ subroutine hypela2(&
use math, only: invnrmMandel use math, only: invnrmMandel
use debug, only: debug_info, & use debug, only: debug_info, &
debug_reset debug_reset
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP, &
mesh_element, &
mesh_node0, &
mesh_node, &
mesh_build_subNodeCoords, &
mesh_build_ipCoordinates, &
FE_Nnodes
use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done
!$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS !$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS
@ -265,7 +272,8 @@ subroutine hypela2(&
real(pReal), dimension (3,3) :: pstress ! dummy argument for call of cpfem_general (used by mpie_spectral) real(pReal), dimension (3,3) :: pstress ! dummy argument for call of cpfem_general (used by mpie_spectral)
real(pReal), dimension (3,3,3,3) :: dPdF ! dummy argument for call of cpfem_general (used by mpie_spectral) real(pReal), dimension (3,3,3,3) :: dPdF ! dummy argument for call of cpfem_general (used by mpie_spectral)
integer(pInt) computationMode, i, cp_en integer(pInt) computationMode, i, cp_en
integer(pInt) node, FEnodeID
! OpenMP variable ! OpenMP variable
!$ integer(pInt) defaultNumThreadsInt ! default value set by Marc !$ integer(pInt) defaultNumThreadsInt ! default value set by Marc
@ -351,6 +359,10 @@ subroutine hypela2(&
else else
computationMode = 3 ! plain collect computationMode = 3 ! plain collect
endif endif
do node = 1,FE_Nnodes(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
endif endif
theTime = cptim ! record current starting time theTime = cptim ! record current starting time
@ -359,7 +371,11 @@ subroutine hypela2(&
lastMode = calcMode(nn,cp_en) ! record calculationMode lastMode = calcMode(nn,cp_en) ! record calculationMode
endif endif
call CPFEM_general(computationMode,dispt,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) ! marc returns nodal coordinates. So for marc we have to calculate the ip coordinates from the nodal coordinates.
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipCoordinates() ! update ip coordinates
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13

22
code/DAMASK_spectral.f90
View File

@ -84,7 +84,8 @@ program DAMASK_spectral
wgt, & wgt, &
geomdim, & geomdim, &
virt_dim, & virt_dim, &
deformed_FFT deformed_FFT, &
mesh_ipCoordinates
use CPFEM, only: & use CPFEM, only: &
CPFEM_general, & CPFEM_general, &
@ -519,15 +520,20 @@ program DAMASK_spectral
if (debugRestart) write(6,'(a)') 'Data read in' if (debugRestart) write(6,'(a)') 'Data read in'
endif endif
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem +1_pInt
mesh_ipCoordinates(1:3,1,ielem) = coordinates(i,j,k,1:3)
enddo; enddo; enddo
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& call CPFEM_general(3_pInt,F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF)
enddo; enddo; enddo enddo; enddo; enddo
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),& call CPFEM_general(2_pInt,F(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3),temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF)
C = C + dPdF C = C + dPdF
enddo; enddo; enddo enddo; enddo; enddo
@ -686,7 +692,11 @@ program DAMASK_spectral
enddo; enddo; enddo enddo; enddo; enddo
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),& ! calculate current coordinates call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,bc(loadcase)%rotation),& ! calculate current coordinates
1.0_pReal,F_lastInc,coordinates) 1.0_pReal,F_lastInc,coordinates)
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem +1_pInt
mesh_ipCoordinates(1:3,1,ielem) = coordinates(i,j,k,1:3)
enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate reduced compliance ! calculate reduced compliance
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
@ -753,7 +763,7 @@ program DAMASK_spectral
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle call CPFEM_general(3_pInt,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3), & F_lastInc(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3), &
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,& temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,&
P_real(i,j,k,1:3,1:3),dPdF) P_real(i,j,k,1:3,1:3),dPdF)
enddo; enddo; enddo enddo; enddo; enddo
@ -765,7 +775,7 @@ program DAMASK_spectral
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), F(i,j,k,1:3,1:3), & ! others get 2 (saves winding forward effort) F_lastInc(i,j,k,1:3,1:3), F(i,j,k,1:3,1:3), & ! others get 2 (saves winding forward effort)
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, & temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, &
P_real(i,j,k,1:3,1:3),dPdF) P_real(i,j,k,1:3,1:3),dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt

23
code/DAMASK_spectral_solverAL.f90
View File

@ -46,7 +46,6 @@ module DAMASK_spectral_solverAL
F_lambda_lastInc, & !< field of previous incompatible deformation gradient F_lambda_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient
F_lambdaDot !< field of assumed rate of incopatible deformation gradient F_lambdaDot !< field of assumed rate of incopatible deformation gradient
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private :: temperature !< temperature, no spatial quantity at the moment real(pReal), private :: temperature !< temperature, no spatial quantity at the moment
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -94,7 +93,8 @@ subroutine AL_init()
use mesh, only: & use mesh, only: &
res, & res, &
geomdim, & geomdim, &
mesh_NcpElems mesh_NcpElems, &
mesh_ipCoordinates
use math, only: & use math, only: &
math_invSym3333 math_invSym3333
@ -119,7 +119,6 @@ subroutine AL_init()
! allocate (Fdot,source = F_lastInc) somethin like that should be possible ! allocate (Fdot,source = F_lastInc) somethin like that should be possible
allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Init ! PETSc Init
@ -151,10 +150,6 @@ subroutine AL_init()
F_lambda_lastInc = F_lastInc F_lambda_lastInc = F_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F_lambda = F F_lambda = F
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
@ -180,9 +175,10 @@ subroutine AL_init()
call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc))
read (777,rec=1) F_aim_lastInc read (777,rec=1) F_aim_lastInc
close (777) close (777)
coordinates = 0.0 ! change it later!!!
endif endif
call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -219,7 +215,7 @@ type(tSolutionState) function &
use mesh, only: & use mesh, only: &
res,& res,&
geomdim,& geomdim,&
deformed_fft mesh_ipCoordinates
use IO, only: & use IO, only: &
IO_write_JobBinaryFile IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
@ -292,7 +288,8 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
@ -300,8 +297,6 @@ type(tSolutionState) function &
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)])
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
endif endif
F_aim = F_aim + f_aimDot * timeinc F_aim = F_aim + f_aimDot * timeinc
@ -406,7 +401,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & call Utilities_constitutiveResponse(F_lastInc,F,temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC) residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False. ForwardData = .False.

32
code/DAMASK_spectral_solverBasic.f90
View File

@ -27,7 +27,6 @@ module DAMASK_spectral_SolverBasic
F, & !< deformation gradient field F, & !< deformation gradient field
F_lastInc, & !< deformation gradient field last increment F_lastInc, & !< deformation gradient field last increment
Fdot !< assumed rate for F n to F n+1 Fdot !< assumed rate for F n to F n+1
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private :: temperature !< not pointwise real(pReal), private :: temperature !< not pointwise
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -63,7 +62,10 @@ subroutine basic_init()
debugRestart debugRestart
use mesh, only: & use mesh, only: &
res, & res, &
geomdim geomdim, &
mesh_ipCoordinates, &
mesh_NcpElems, &
mesh_deformedCoordsFFT
implicit none implicit none
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
@ -86,18 +88,12 @@ subroutine basic_init()
allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields and average quantities ! init fields and average quantities
if (restartInc == 1_pInt) then ! no deformation (no restart) if (restartInc == 1_pInt) then ! no deformation (no restart)
F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lastInc = F F_lastInc = F
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') & if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'Reading values of increment', restartInc - 1_pInt, 'from file' 'Reading values of increment', restartInc - 1_pInt, 'from file'
@ -127,11 +123,10 @@ subroutine basic_init()
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
read (777,rec=1) temp3333_Real read (777,rec=1) temp3333_Real
close (777) close (777)
coordinates = 0.0 ! change it later!!!
endif endif
call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C temp3333_Real = C
endif endif
@ -158,8 +153,10 @@ type(tSolutionState) function &
use mesh, only: & use mesh, only: &
res,& res,&
geomdim, & geomdim, &
deformed_fft, & wgt, &
wgt mesh_ipCoordinates,&
mesh_NcpElems, &
mesh_deformedCoordsFFT
use IO, only: & use IO, only: &
IO_write_JobBinaryFile, & IO_write_JobBinaryFile, &
IO_intOut IO_intOut
@ -252,6 +249,8 @@ type(tSolutionState) function &
C = C_lastInc C = C_lastInc
else else
C_lastInc = C C_lastInc = C
mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
@ -259,9 +258,6 @@ type(tSolutionState) function &
endif endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & ! update coordinates and rate and forward last inc
1.0_pReal,F_lastInc,coordinates)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,F) timeinc,timeinc_old,guess,F_lastInc,F)
F_lastInc = F F_lastInc = F
@ -291,7 +287,7 @@ type(tSolutionState) function &
! evaluate constitutive response ! evaluate constitutive response
F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC)
basic_solution%termIll = .false. basic_solution%termIll = .false.
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& call Utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
P,C,P_av,ForwardData,rotation_BC) P,C,P_av,ForwardData,rotation_BC)
basic_solution%termIll = terminallyIll basic_solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.

49
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -10,10 +10,8 @@ module DAMASK_spectral_SolverBasicPETSc
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal
use math, only: & use math, only: &
math_I3 math_I3
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
tSolutionState tSolutionState
@ -44,7 +42,6 @@ module DAMASK_spectral_SolverBasicPETSc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal) :: temperature real(pReal) :: temperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -72,33 +69,27 @@ module DAMASK_spectral_SolverBasicPETSc
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info !> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine basicPETSc_init() subroutine basicPETSc_init()
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: &
IO_read_JobBinaryFile, & IO_read_JobBinaryFile, &
IO_write_JobBinaryFile IO_write_JobBinaryFile
use FEsolving, only: & use FEsolving, only: &
restartInc restartInc
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
Utilities_init, & Utilities_init, &
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
debugRestart debugRestart
use numerics, only: & use numerics, only: &
petsc_options petsc_options
use mesh, only: & use mesh, only: &
res, & res, &
geomdim, & geomdim, &
mesh_NcpElems mesh_NcpElems, &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
use math, only: & use math, only: &
math_invSym3333 math_invSym3333
@ -106,7 +97,7 @@ subroutine basicPETSc_init()
#include <finclude/petscdmda.h90> #include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90> #include <finclude/petscsnes.h90>
integer(pInt) :: i,j,k integer(pInt) :: i,j,k
real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, dimension(:,:,:,:), pointer :: F
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscObject :: dummy PetscObject :: dummy
@ -119,9 +110,8 @@ subroutine basicPETSc_init()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc(3,3,res(1),res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot (3,3,res(1),res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
@ -150,10 +140,6 @@ subroutine basicPETSc_init()
if (restartInc == 1_pInt) then ! no deformation (no restart) if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
@ -172,11 +158,10 @@ subroutine basicPETSc_init()
call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc))
read (777,rec=1) F_aim_lastInc read (777,rec=1) F_aim_lastInc
close (777) close (777)
coordinates = 0.0 ! change it later!!!
endif endif
mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
call Utilities_constitutiveResponse(coordinates,& !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
temperature,0.0_pReal,P,C,P_av,.false.,math_I3) temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
@ -211,7 +196,9 @@ type(tSolutionState) function &
use mesh, only: & use mesh, only: &
res,& res,&
geomdim,& geomdim,&
deformed_fft mesh_ipCoordinates,&
mesh_NcpElems, &
mesh_deformedCoordsFFT
use IO, only: & use IO, only: &
IO_write_JobBinaryFile IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
@ -260,14 +247,16 @@ type(tSolutionState) function &
close(777) close(777)
endif endif
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(F,[3,3,res(1),res(2),res(3)])),&
[3,1,mesh_NcpElems])
if ( cutBack) then if ( cutBack) then
F_aim = F_aim_lastInc F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc C = C_lastInc
else else
C_lastInc = C C_lastInc = C
mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
@ -280,8 +269,6 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
@ -292,7 +279,6 @@ type(tSolutionState) function &
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
@ -322,7 +308,6 @@ end function BasicPETSc_solution
!> @brief forms the AL residual vector !> @brief forms the AL residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr) subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin itmin
@ -366,7 +351,7 @@ subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, & call Utilities_constitutiveResponse(F_lastInc,x_scal,temperature,params%timeinc, &
f_scal,C,P_av,ForwardData,params%rotation_BC) f_scal,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .false. ForwardData = .false.

7
code/DAMASK_spectral_utilities.f90
View File

@ -630,7 +630,7 @@ end function utilities_maskedCompliance
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response !> @brief calculates constitutive response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
P,C,P_av,forwardData,rotation_BC) P,C,P_av,forwardData,rotation_BC)
use debug, only: & use debug, only: &
debug_reset, & debug_reset, &
@ -648,7 +648,6 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
implicit none implicit none
real(pReal), intent(inout) :: temperature !< temperature (no field) real(pReal), intent(inout) :: temperature !< temperature (no field)
real(pReal), intent(in), dimension(res(1),res(2),res(3),3) :: coordinates !< coordinates field
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: &
F_lastInc, & !< target deformation gradient F_lastInc, & !< target deformation gradient
F !< previous deformation gradient F !< previous deformation gradient
@ -702,7 +701,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(collectMode,& ! collect cycle call CPFEM_general(collectMode,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), &
temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
collectMode = 3_pInt collectMode = 3_pInt
enddo; enddo; enddo enddo; enddo; enddo
@ -714,7 +713,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort)
temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
calcMode = 2_pInt calcMode = 2_pInt
C = C + dPdF C = C + dPdF

253
code/mesh.f90
View File

@ -1,7 +1,7 @@
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH ! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
! !
! This file is part of DAMASK, ! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit. ! the Düsseldorf Advanced Material Simulation Kit.
! !
! DAMASK is free software: you can redistribute it and/or modify ! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by ! it under the terms of the GNU General Public License as published by
@ -19,11 +19,11 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!* $Id$ !* $Id$
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !! Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH !! Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !! Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!! Krishna Komerla, Max-Planck-Institut für Eisenforschung GmbH !! 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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -286,6 +286,7 @@ module mesh
mesh_regular_grid, & mesh_regular_grid, &
deformed_linear, & deformed_linear, &
deformed_fft, & deformed_fft, &
mesh_deformedCoordsFFT, &
volume_compare, & volume_compare, &
shape_compare shape_compare
#endif #endif
@ -1184,7 +1185,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
integer(pInt), dimension(:,:,:), allocatable :: & integer(pInt), dimension(:,:,:), allocatable :: &
material_phase, material_phaseNew, & material_phase, material_phaseNew, &
sizeStateConst sizeStateConst
write(6,*) 'Regridding geometry' write(6,*) 'Regridding geometry'
if (adaptive) then if (adaptive) then
write(6,*) 'adaptive resolution determination' write(6,*) 'adaptive resolution determination'
@ -1355,7 +1356,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
allocate(F_Linear(3,3,mesh_NcpElems)) allocate(F_Linear(3,3,mesh_NcpElems))
allocate(F_Linear_New(3,3,NpointsNew)) allocate(F_Linear_New(3,3,NpointsNew))
allocate(FNew(resNew(1),resNew(2),resNew(3),3,3)) allocate(FNew(resNew(1),resNew(2),resNew(3),3,3))
ielem = 0_pInt ielem = 0_pInt
do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1) do k=1_pInt,res(3); do j=1_pInt, res(2); do i=1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
@ -1365,7 +1366,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
do i=1_pInt, NpointsNew do i=1_pInt, NpointsNew
F_Linear_New(1:3,1:3,i) = F_Linear(1:3,1:3,indices(i)) ! -- mapping old to new ...based on indices F_Linear_New(1:3,1:3,i) = F_Linear(1:3,1:3,indices(i)) ! -- mapping old to new ...based on indices
enddo enddo
ielem = 0_pInt ielem = 0_pInt
do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
@ -1377,11 +1378,11 @@ function mesh_regrid(adaptive,resNewInput,minRes)
enddo; enddo enddo; enddo
deltaF = Favg - FavgNew deltaF = Favg - FavgNew
do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
FNew(i,j,k,1:3,1:3) = FNew(i,j,k,1:3,1:3) + deltaF FNew(i,j,k,1:3,1:3) = FNew(i,j,k,1:3,1:3) + deltaF
enddo; enddo; enddo enddo; enddo; enddo
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(FNew)) call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(FNew))
write (777,rec=1) FNew write (777,rec=1) FNew
close (777) close (777)
@ -1394,12 +1395,12 @@ function mesh_regrid(adaptive,resNewInput,minRes)
allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3)) allocate(F_lastIncNew(resNew(1),resNew(2),resNew(3),3,3))
allocate(F_Linear(3,3,mesh_NcpElems)) allocate(F_Linear(3,3,mesh_NcpElems))
allocate(F_Linear_New(3,3,NpointsNew)) allocate(F_Linear_New(3,3,NpointsNew))
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc', & call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc', &
trim(getSolverJobName()),size(F_lastInc)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc
close (777) close (777)
call IO_read_jobBinaryFile(777,'F_aim_lastInc', & call IO_read_jobBinaryFile(777,'F_aim_lastInc', &
trim(getSolverJobName()),size(Favg_LastInc)) trim(getSolverJobName()),size(Favg_LastInc))
read (777,rec=1) Favg_LastInc read (777,rec=1) Favg_LastInc
@ -1433,7 +1434,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1) do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
F_LastIncNew(i,j,k,1:3,1:3) = F_LastIncNew(i,j,k,1:3,1:3) + deltaF_lastInc F_LastIncNew(i,j,k,1:3,1:3) = F_LastIncNew(i,j,k,1:3,1:3) + deltaF_lastInc
enddo; enddo; enddo enddo; enddo; enddo
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_LastIncNew)) call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_LastIncNew))
write (777,rec=1) F_LastIncNew write (777,rec=1) F_LastIncNew
close (777) close (777)
@ -1441,7 +1442,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
deallocate(F_Linear_New) deallocate(F_Linear_New)
deallocate(F_lastInc) deallocate(F_lastInc)
deallocate(F_lastIncNew) deallocate(F_lastIncNew)
! relocating data of material subroutine --------------------------------------------------------- ! relocating data of material subroutine ---------------------------------------------------------
allocate(material_phase (1,1, mesh_NcpElems)) allocate(material_phase (1,1, mesh_NcpElems))
allocate(material_phaseNew (1,1, NpointsNew)) allocate(material_phaseNew (1,1, NpointsNew))
@ -1924,6 +1925,156 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
end subroutine deformed_fft end subroutine deformed_fft
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function mesh_deformedCoordsFFT(geomdim,F,scalingIn,FavgIn)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Routine to calculate coordinates in current configuration for given defgrad
! using integration in Fourier space (more accurate than deformed(...))
!
use IO, only: &
IO_error
use numerics, only: &
fftw_timelimit, &
fftw_planner_flag
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
PI
implicit none
real(pReal), intent(in), dimension(3) :: geomdim
real(pReal), intent(in), dimension(:,:,:,:,:) :: F
real(pReal), intent(in), dimension(3,3), optional :: FavgIn
real(pReal), intent(in), optional :: scalingIn
! function
real(pReal), dimension(3,size(F,3),size(F,4),size(F,5)) :: mesh_deformedCoordsFFT
! allocatable arrays for fftw c routines
type(C_PTR) :: fftw_forth, fftw_back
type(C_PTR) :: coords_fftw, defgrad_fftw
real(pReal), dimension(:,:,:,:,:), pointer :: F_real
complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier
real(pReal), dimension(:,:,:,:), pointer :: coords_real
complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier
! other variables
integer(pInt) :: i, j, k, m, res1_red
integer(pInt), dimension(3) :: k_s, res
real(pReal), dimension(3) :: step, offset_coords, integrator
real(pReal), dimension(3,3) :: Favg
real(pReal) :: scaling
if (present(scalingIn)) then
if (scalingIn < 0.0_pReal) then !the f2py way to tell it is not present
scaling = 1.0_pReal
else
scaling = scalingIn
endif
else
scaling = 1.0_pReal
endif
res = [size(F,3),size(F,4),size(F,5)]
integrator = geomdim / 2.0_pReal / pi ! see notes where it is used
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
print*, 'Restore geometry using FFT-based integration'
print '(a,3(e12.5))', ' Dimension: ', geomdim
print '(a,3(i5))', ' Resolution:', res
endif
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = geomdim/real(res, pReal)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit)
defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, F_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt])
call c_f_pointer(defgrad_fftw, F_fourier,[res1_red ,res(2),res(3),3_pInt,3_pInt])
coords_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(coords_fftw, coords_real, [res(1)+2_pInt,res(2),res(3),3_pInt])
call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,res(2),res(3),3_pInt])
fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),9_pInt,& ! dimensions , length in each dimension in reversed order
F_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order
1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
F_fourier,(/res(3),res(2) ,res1_red/),&
1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag)
fftw_back = fftw_plan_many_dft_c2r(3_pInt,(/res(3),res(2) ,res(1)/),3_pInt,&
coords_fourier,(/res(3),res(2) ,res1_red/),&
1_pInt, res(3)*res(2)* res1_red,&
coords_real,(/res(3),res(2) ,res(1)+2_pInt/),&
1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
F_real(i,j,k,1:3,1:3) = F(1:3,1:3,i,j,k) ! ensure that data is aligned properly (fftw_alloc)
enddo; enddo; enddo
call fftw_execute_dft_r2c(fftw_forth, F_real, F_fourier)
if (present(FavgIn)) then
if (all(FavgIn < 0.0_pReal)) then
Favg = real(F_fourier(1,1,1,1:3,1:3)*real((res(1)*res(2)*res(3)),pReal),pReal) !the f2py way to tell it is not present
else
Favg = FavgIn
endif
else
Favg = real(F_fourier(1,1,1,1:3,1:3)*real((res(1)*res(2)*res(3)),pReal),pReal)
endif
!remove highest frequency in each direction
if(res(1)>1_pInt) &
F_fourier( res(1)/2_pInt+1_pInt,1:res(2) ,1:res(3) ,&
1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(2)>1_pInt) &
F_fourier(1:res1_red ,res(2)/2_pInt+1_pInt,1:res(3) ,&
1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)>1_pInt) &
F_fourier(1:res1_red ,1:res(2) ,res(3)/2_pInt+1_pInt,&
1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, res(3)
k_s(3) = k-1_pInt
if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3)
do j = 1_pInt, res(2)
k_s(2) = j-1_pInt
if(j > res(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-res(2)
do i = 1_pInt, res1_red
k_s(1) = i-1_pInt
do m = 1_pInt,3_pInt
coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal))
enddo
if (k_s(3) /= 0_pInt .or. k_s(2) /= 0_pInt .or. k_s(1) /= 0_pInt) &
coords_fourier(i,j,k,1:3) = coords_fourier(i,j,k,1:3) / real(-sum(k_s*k_s),pReal)
enddo; enddo; enddo
call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real)
coords_real = coords_real/real(res(1)*res(2)*res(3),pReal)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
mesh_deformedCoordsFFT(1:3,i,j,k) = coords_real(i,j,k,1:3) ! ensure that data is aligned properly (fftw_alloc)
enddo; enddo; enddo
offset_coords = matmul(F(1:3,1:3,1,1,1),step/2.0_pReal) - scaling*mesh_deformedCoordsFFT(1:3,1,1,1)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
mesh_deformedCoordsFFT(1:3,i,j,k) = scaling*mesh_deformedCoordsFFT(1:3,i,j,k) &
+ offset_coords + matmul(Favg,&
(/step(1)*real(i-1_pInt,pReal),&
step(2)*real(j-1_pInt,pReal),&
step(3)*real(k-1_pInt,pReal)/))
enddo; enddo; enddo
call fftw_destroy_plan(fftw_forth)
call fftw_destroy_plan(fftw_back)
call fftw_free(defgrad_fftw)
call fftw_free(coords_fftw)
end function mesh_deformedCoordsFFT
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -3711,28 +3862,28 @@ do e = 1_pInt,mesh_NcpElems
enddo enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) then if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
write (6,*) write(6,*)
write (6,*) 'Input Parser: STATISTICS' write(6,*) 'Input Parser: STATISTICS'
write (6,*) write(6,*)
write (6,*) mesh_Nelems, ' : total number of elements in mesh' 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_maxNnodes, ' : max number of nodes in any CP element'
write (6,*) mesh_maxNips, ' : max number of IPs 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_maxNipNeighbors, ' : max number of IP neighbors in any CP element'
write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element' write(6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element'
write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node' write(6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node'
write (6,*) write(6,*)
write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE' write(6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE'
write (6,*) write(6,*)
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'
write (6,*) write(6,*)
write (myFmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))' write (myFmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
write (6,myFmt) '+-',math_range(mesh_maxValStateVar(2)) write(6,myFmt) '+-',math_range(mesh_maxValStateVar(2))
write (myFmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))' write (myFmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations do i=1_pInt,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations
write (6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures write(6,myFmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures
enddo enddo
write(6,*) write(6,*)
write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS' write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS'
@ -3743,9 +3894,9 @@ enddo
endif endif
if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then
write (6,*) write(6,*)
write (6,*) 'Input Parser: SUBNODE COORDINATES' write(6,*) 'Input Parser: SUBNODE COORDINATES'
write (6,*) write(6,*)
write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')& write(6,'(a8,1x,a5,1x,2(a15,1x),a20,3(1x,a12))')&
'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z'
do e = 1_pInt,mesh_NcpElems ! loop over cpElems do e = 1_pInt,mesh_NcpElems ! loop over cpElems
@ -3770,28 +3921,28 @@ enddo
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle
do i = 1_pInt,FE_Nips(mesh_element(2,e)) do i = 1_pInt,FE_Nips(mesh_element(2,e))
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle
write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e) write(6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e)
enddo enddo
enddo enddo
write (6,*) write(6,*)
write (6,*) 'Input Parser: ELEMENT VOLUME' write(6,*) 'Input Parser: ELEMENT VOLUME'
write (6,*) write(6,*)
write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume) write(6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume)
write (6,*) write(6,*)
write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --' write(6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --'
do e = 1_pInt,mesh_NcpElems do e = 1_pInt,mesh_NcpElems
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle
do i = 1_pInt,FE_Nips(mesh_element(2,e)) do i = 1_pInt,FE_Nips(mesh_element(2,e))
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle
write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e) write(6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e)
do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) do f = 1_pInt,FE_NipNeighbors(mesh_element(2,e))
write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) write(6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
enddo enddo
enddo enddo
enddo enddo
write (6,*) write(6,*)
write (6,*) 'Input Parser: NODE TWINS' write(6,*) 'Input Parser: NODE TWINS'
write (6,*) write(6,*)
write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z' write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z'
do n = 1_pInt,mesh_Nnodes ! loop over cpNodes do n = 1_pInt,mesh_Nnodes ! loop over cpNodes
if (debug_e <= mesh_NcpElems) then if (debug_e <= mesh_NcpElems) then
@ -3810,7 +3961,7 @@ enddo
do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem do i = 1_pInt,FE_Nips(t) ! loop over IPs of elem
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle
do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP do n = 1_pInt,FE_NipNeighbors(t) ! loop over neighbors of IP
write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) write(6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
enddo enddo
enddo enddo
enddo enddo