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

ip coordinates are now updated every cycle; this is needed for the nonlocal internal stress fields 

* Marc: node displacements are added to initial node coordinates (mesh_node0) to get current node positions (mesh_node), then ip coordinates are deduced
* Abaqus: ip coordinates are directly updated, no update of node coordinates!
* Spectral: for the moment no update of either ip or node coordinates! passing only dummy values with initial ip coordinates
This commit is contained in:
Christoph Kords 2011-05-24 15:57:59 +00:00
parent 73b1dbc5f7
commit cac45cff96
8 changed files with 174 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
code/CPFEM.f90
View File

@ -224,7 +224,7 @@ endsubroutine
!*** perform initialization at first call, update variables and ***
!*** call the actual material model ***
!***********************************************************************
subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,&
subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,&
& jacobian, pstress, dPdF)
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE
@ -270,7 +270,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
FE_Nips
mesh_node0, &
mesh_node, &
mesh_ipCenterOfGravity, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
FE_Nips, &
FE_Nnodes
use material, only: homogenization_maxNgrains, &
microstructure_elemhomo, &
material_phase
@ -284,7 +290,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v
crystallite_Tstar_v, &
crystallite_localConstitution
use homogenization, only: homogenization_sizeState, &
homogenization_state, &
homogenization_state0, &
@ -303,11 +310,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
implicit none
!*** input variables ***!
integer(pInt), intent(in) :: element, & ! FE element number
IP ! FE integration point number
real(pReal), intent(inout) :: Temperature ! temperature
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
ffn1 ! deformation gradient for t=t1
integer(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results
@ -317,13 +329,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! 5: restore tangent from former converged inc
! 6: recycling of former results (MARC speciality)
!*** output variables ***!
real(pReal), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation
real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation
real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation
real(pReal), dimension (3,3,3,3), intent(out) :: dPdF !
!*** local variables ***!
real(pReal) J_inverse, & ! inverse of Jacobian
rnd
real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress in Matrix notation
@ -338,17 +354,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
l, &
m, &
n, &
e
e, &
node, &
FEnodeID
logical updateJaco ! flag indicating if JAcobian has to be updated
!*** global variables ***!
! CPFEM_cs, &
! CPFEM_dcsdE, &
! CPFEM_dcsdE_knownGood, &
! CPFEM_init_done, &
! CPFEM_calc_done, &
! CPFEM_odd_stress, &
! CPFEM_odd_jacobian
cp_en = mesh_FEasCP('elem',element)
@ -367,12 +377,18 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
!$OMP END CRITICAL (write2out)
endif
! according to our "mode" we decide what to do
!*** according to our "mode" we decide what to do
select case (mode)
! --+>> REGULAR COMPUTATION (WITH AGING OF RESULTS IF MODE == 1) <<+--
case (1,2,8,9)
! age results if mode == 1
!*** age results
if (mode == 1 .or. mode == 8) then
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
@ -402,8 +418,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
enddo
!$OMP END PARALLEL DO
! *** dump the last converged values of each essential variable to a binary file
! * dump the last converged values of each essential variable to a binary file
if (restartWrite) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
@ -459,7 +475,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
close (777)
endif
endif
! *** end of dumping
! * end of dumping
endif
if (mode == 8 .or. mode == 9) then ! Abaqus explicit skips collect
@ -468,7 +484,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
materialpoint_F(1:3,1:3,IP,cp_en) = ffn1
endif
! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then
if (.not. terminallyIll .and. .not. outdatedFFN1) then
if (debug_verbosity > 0) then
@ -485,14 +503,15 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
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)
! deformation gradient is not outdated
!*** deformation gradient is not outdated
else
! set flag for Jacobian update
updateJaco = mod(cycleCounter,iJacoStiffness) == 0
! no parallel computation
!* no parallel computation, so we use just one single element and IP for computation
if (.not. parallelExecution) then
! we just take one single element and IP
FEsolving_execElem(1) = cp_en
FEsolving_execElem(2) = cp_en
FEsolving_execIP(1,cp_en) = IP
@ -505,13 +524,21 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
call materialpoint_postResults(dt) ! post results
! parallel computation and calulation not yet done
!* parallel computation and calulation not yet done
elseif (.not. CPFEM_calc_done) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i8)') '<< CPFEM >> Calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
endif
if (any(.not. crystallite_localConstitution) .and. FEsolver == 'Marc') then
call mesh_build_subNodeCoords() ! update subnodal coordinates
call mesh_build_ipVolumes() ! update ip center of gravity
endif
!$OMP CRITICAL (write2out)
write(6,'(a,i5,a,i5)') '<< CPFEM >> Start stress and tangent ',FEsolving_execElem(1),' to ',FEsolving_execElem(2)
!$OMP END CRITICAL (write2out)
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults(dt) ! post results
!$OMP PARALLEL DO
@ -529,6 +556,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_calc_done = .true.
endif
!*** map stress and stiffness (or return odd values if terminally ill)
if ( terminallyIll ) then
call random_number(rnd)
rnd = 2.0_pReal * rnd - 1.0_pReal
@ -558,7 +588,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
endif
! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+--
case (3,4,5)
if (mode == 4) then
CPFEM_dcsde_knownGood = CPFEM_dcsde ! --+>> BACKUP JACOBIAN FROM FORMER CONVERGED INC
@ -573,29 +605,41 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
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_calc_done = .false.
select case (FEsolver)
case ('Abaqus','Spectral')
mesh_ipCenterOfGravity(1:3,IP,cp_en) = coords(1:3,1)
case ('Marc')
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) + coords(1:3,node)
enddo
end select
! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+--
case (6)
! do nothing
! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC
case (7)
CPFEM_dcsde = CPFEM_dcsde_knownGood
end select
! return the local stress and the jacobian from storage
!*** fill output variables with computed values
cauchyStress = CPFEM_cs(1:6,IP,cp_en)
jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en)
! copy P and dPdF to the output variables
pstress = materialpoint_P(1:3,1:3,IP,cp_en)
dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en)
! warning for zero stiffness
if (all(abs(jacobian) < 1e-10_pReal)) then
call IO_warning(601,cp_en,IP)
if (theTime > 0.0_pReal) then
Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
endif
if (mode < 6 .and. debug_verbosity > 0 .and. ((debug_e == cp_en .and. debug_i == IP) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,x,i2,/,12(x),6(f10.3,x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1e6
@ -604,7 +648,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
!$OMP END CRITICAL (write2out)
endif
! remember extreme values of stress and jacobian
!*** warn if stiffness close to zero
if (all(abs(jacobian) < 1e-10_pReal)) then
call IO_warning(601,cp_en,IP)
endif
!*** remember extreme values of stress and jacobian
if (mode < 3) then
cauchyStress33 = math_Mandel6to33(cauchyStress)
if (maxval(cauchyStress33) > debug_stressMax) then
@ -626,9 +679,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
endif
! return temperature
if (theTime > 0.0_pReal) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
end subroutine
END MODULE CPFEM

4
code/DAMASK_abaqus_exp.f
View File

@ -193,7 +193,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
stateOld(nblock,nstatev), enerInternOld(nblock), &
enerInelasOld(nblock), tempNew(nblock), tempOld(nblock), &
stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), &
stressNew(nblock,ndir+nshr)
stressNew(nblock,ndir+nshr), coordMp(nblock,3)
dimension enerInelasNew(nblock),stateNew(nblock,nstatev),enerInternNew(nblock)
dimension nElement(nblock),nMatPoint(nblock)
@ -283,7 +283,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
defgrd1(3,2) = defgradNew(n,8)
endif
call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF)
call CPFEM_general(computationMode,coordMp(n,1:3),defgrd0,defgrd1,temp,timeInc,nElement(n),nMatPoint(n),stress,ddsdde, pstress, dPdF)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! straight: 11, 22, 33, 12, 23, 13

2
code/DAMASK_abaqus_std.f
View File

@ -261,7 +261,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
!$OMP END CRITICAL (write2out)
endif
call CPFEM_general(computationMode,dfgrd0,dfgrd1,temp,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF)
call CPFEM_general(computationMode,COORDS,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
! straight: 11, 22, 33, 12, 23, 13

2
code/DAMASK_marc.f90
View File

@ -353,7 +353,7 @@ subroutine hypela2(&
lastMode = calcMode(nn,cp_en) ! record calculationMode
endif
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF)
call CPFEM_general(computationMode,dispt,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
! Marc: 11, 22, 33, 12, 23, 13

13
code/DAMASK_spectral.f90
View File

@ -91,6 +91,7 @@ program DAMASK_spectral
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation
real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors
real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold
real(pReal), dimension(:,:,:,:), allocatable :: coordinates
! variables storing information for spectral method
complex(pReal) :: img
@ -261,6 +262,7 @@ program DAMASK_spectral
allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal
wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal)
defgradAim = math_I3
@ -275,7 +277,8 @@ program DAMASK_spectral
defgradold(i,j,k,:,:) = math_I3 ! no deformation at the beginning
defgrad(i,j,k,:,:) = math_I3
ielem = ielem +1
call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates (SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!!)
call CPFEM_general(2,coordinates(1:3,i,j,k),math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c066 = c066 + dsde
enddo; enddo; enddo
c066 = c066 * wgt
@ -410,7 +413,7 @@ program DAMASK_spectral
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1
call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
call CPFEM_general(3, coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
enddo; enddo; enddo
@ -419,6 +422,7 @@ program DAMASK_spectral
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(1:3,i,j,k),
defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& ! others get 2 (saves winding forward effort)
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
@ -469,7 +473,7 @@ program DAMASK_spectral
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
call CPFEM_general(3, coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
enddo; enddo; enddo
@ -477,6 +481,7 @@ program DAMASK_spectral
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(1:3,i,j,k),
defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
@ -594,4 +599,4 @@ subroutine quit(id)
integer(pInt) id
stop
end subroutine
end subroutine

13
code/DAMASK_spectral_single.f90
View File

@ -91,6 +91,7 @@ program DAMASK_spectral
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation
real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors
real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold
real(pReal), dimension(:,:,:,:), allocatable :: coordinates
! variables storing information for spectral method
complex(pReal) :: img
@ -261,6 +262,7 @@ program DAMASK_spectral
allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal
wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal)
defgradAim = math_I3
@ -275,7 +277,8 @@ program DAMASK_spectral
defgradold(i,j,k,:,:) = math_I3 !no deformation at the beginning
defgrad(i,j,k,:,:) = math_I3
ielem = ielem +1
call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates (SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!!)
call CPFEM_general(2,coordinates(1:3,i,j,k),math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c066 = c066 + dsde
enddo; enddo; enddo
c066 = c066 * wgt
@ -411,7 +414,7 @@ program DAMASK_spectral
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1
call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
call CPFEM_general(3, coordinates(1:3,i,j,k),defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
enddo; enddo; enddo
@ -420,6 +423,7 @@ program DAMASK_spectral
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(1:3,i,j,k),
defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& ! others get 2 (saves winding forward effort)
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
@ -470,7 +474,7 @@ program DAMASK_spectral
ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
call CPFEM_general(3, coordinates(1:3,i,j,k),defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
enddo; enddo; enddo
@ -478,6 +482,7 @@ program DAMASK_spectral
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(1:3,i,j,k),
defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),&
temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF)
@ -595,4 +600,4 @@ subroutine quit(id)
integer(pInt) id
stop
end subroutine
end subroutine

6
code/constitutive_nonlocal.f90
View File

@ -837,7 +837,7 @@ use debug, only: debug_verbosity, &
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
mesh_node, &
mesh_node0, &
FE_Nips, &
mesh_ipCenterOfGravity, &
mesh_ipVolume, &
@ -968,8 +968,8 @@ if (.not. phase_localConstitution(phase)) then
!* in case of periodic surfaces we have to find out how many periodic images in each direction we need
do dir = 1,3
maxCoord(dir) = maxval(mesh_node(dir,:))
minCoord(dir) = minval(mesh_node(dir,:))
maxCoord(dir) = maxval(mesh_node0(dir,:))
minCoord(dir) = minval(mesh_node0(dir,:))
enddo
meshSize = maxCoord - minCoord
ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)

109
code/mesh.f90
View File

@ -35,7 +35,8 @@
! _maxNsharedElems : max number of CP elements sharing a node
!
! _element : FEid, type(internal representation), material, texture, node indices
! _node : x,y,z coordinates (initially!)
! _node0 : x,y,z coordinates (initially!)
! _node : x,y,z coordinates (after deformation!)
! _sharedElem : entryCount and list of elements containing node
!
! _mapFEtoCPelem : [sorted FEid, corresponding CPid]
@ -56,7 +57,6 @@
! _nodeOnFace : list of node indices on each face of a specific type of element
! _maxNnodesAtIP : max number of (equivalent) nodes attached to an IP
! _nodesAtIP : map IP index to two node indices in a specific type of element
! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index]
! _NsubNodes : # subnodes required to fully define all IP volumes
! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype
@ -70,19 +70,20 @@
mesh_mapMaterial ! name of elementSet for material
integer(pInt), dimension(:,:), allocatable :: mesh_mapElemSet ! list of elements in elementSet
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem, mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable :: mesh_element, &
mesh_sharedElem, &
integer(pInt), dimension(:,:), allocatable :: mesh_element, & ! FEid, type(internal representation), material, texture, node indices
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)
integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood
integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood ! 6 or less neighboring IPs as [element_num, IP_index]
real(pReal), dimension(:,:,:), allocatable :: mesh_subNodeCoord ! coordinates of subnodes per element
real(pReal), dimension(:,:), allocatable :: mesh_node, &
mesh_ipVolume ! volume associated with IP
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP
mesh_ipCenterOfGravity ! center of gravity of IP
real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP
real(pReal), dimension(:,:), allocatable :: mesh_node0, & ! node coordinates (initially!)
mesh_node, & ! node coordinates (after deformation! ONLY FOR MARC!!!)
mesh_ipVolume ! volume associated with IP (initially!)
real(pReal), dimension(:,:,:), allocatable :: mesh_ipArea, & ! area of interface to neighboring IP (initially!)
mesh_ipCenterOfGravity ! center of gravity of IP (after deformation!)
real(pReal), dimension(:,:,:,:), allocatable :: mesh_ipAreaNormal ! area normal of interface to neighboring IP (initially!)
integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP
integer(pInt), dimension(:,:,:), allocatable :: FE_nodesAtIP ! map IP index to two node indices in a specific type of element
integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor
integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent
integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace
@ -776,7 +777,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
/),(/FE_NipNeighbors(10),FE_Nips(10)/))
! *** FE_subNodeParent ***
! lists the group of IPs for which the center of gravity
! lists the group of nodes for which the center of gravity
! corresponds to the location of a each subnode.
! example: face-centered subnode with faceNodes 1,2,3,4 to be used in,
! e.g., a 8 IP grid, would be encoded:
@ -2417,6 +2418,7 @@ subroutine mesh_marc_count_cpSizes (unit)
character(len=64) tag
character(len=1024) line
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt
a = 1_pInt
@ -2473,10 +2475,12 @@ subroutine mesh_marc_count_cpSizes (unit)
if (x <= 0.0_pReal .or. y <= 0.0_pReal .or. z <= 0.0_pReal) call IO_error(44)
forall (n = 0:mesh_Nnodes-1)
mesh_node(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal))
mesh_node(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal))
mesh_node(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal))
mesh_node0(1,n+1) = x * dble(mod(n,a) / (a-1.0_pReal))
mesh_node0(2,n+1) = y * dble(mod(n/a,b) / (b-1.0_pReal))
mesh_node0(3,n+1) = z * dble(mod(n/a/b,c) / (c-1.0_pReal))
end forall
mesh_node = mesh_node0
100 return
@ -2502,6 +2506,7 @@ subroutine mesh_marc_count_cpSizes (unit)
integer(pInt) unit,i,j,m
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt
610 FORMAT(A300)
@ -2515,13 +2520,14 @@ subroutine mesh_marc_count_cpSizes (unit)
do i=1,mesh_Nnodes
read (unit,610,END=670) line
m = mesh_FEasCP('node',IO_fixedIntValue(line,node_ends,1))
forall (j = 1:3) mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1)
forall (j = 1:3) mesh_node0(j,m) = IO_fixedNoEFloatValue(line,node_ends,j+1)
enddo
exit
endif
enddo
670 return
670 mesh_node = mesh_node0
return
endsubroutine
@ -2545,6 +2551,7 @@ subroutine mesh_marc_count_cpSizes (unit)
integer(pInt) unit,i,j,m,count
logical inPart
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0_pInt
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0_pInt
610 FORMAT(A300)
@ -2573,12 +2580,13 @@ subroutine mesh_marc_count_cpSizes (unit)
read (unit,610,END=670) line
pos = IO_stringPos(line,maxNchunks)
m = mesh_FEasCP('node',IO_intValue(line,pos,1))
forall (j=1:3) mesh_node(j,m) = IO_floatValue(line,pos,j+1)
forall (j=1:3) mesh_node0(j,m) = IO_floatValue(line,pos,j+1)
enddo
endif
enddo
670 if (size(mesh_node,2) /= mesh_Nnodes) call IO_error(909)
670 if (size(mesh_node0,2) /= mesh_Nnodes) call IO_error(909)
mesh_node = mesh_node0
return
endsubroutine
@ -3098,25 +3106,26 @@ endsubroutine
integer(pInt) e,t,n,p
allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems)) ; mesh_subNodeCoord = 0.0_pReal
if (.not. allocated(mesh_subNodeCoord)) then
allocate(mesh_subNodeCoord(3,mesh_maxNnodes+mesh_maxNsubNodes,mesh_NcpElems))
endif
mesh_subNodeCoord = 0.0_pReal
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
do n = 1,FE_Nnodes(t)
mesh_subNodeCoord(:,n,e) = mesh_node(:,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type
mesh_subNodeCoord(1:3,n,e) = mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+n,e))) ! loop over nodes of this element type
enddo
do n = 1,FE_NsubNodes(t) ! now for the true subnodes
do p = 1,FE_Nips(t) ! loop through possible parent nodes
if (FE_subNodeParent(p,n,t) > 0) & ! valid parent node
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = &
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) + &
mesh_node(:,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents
mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) = mesh_subNodeCoord(1:3,FE_Nnodes(t)+n,e) &
+ mesh_node(1:3,mesh_FEasCP('node',mesh_element(4+FE_subNodeParent(p,n,t),e))) ! add up parents
enddo
mesh_subNodeCoord(:,n+FE_Nnodes(t),e) = mesh_subNodeCoord(:,n+FE_Nnodes(t),e) / count(FE_subNodeParent(:,n,t) > 0)
mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) = mesh_subNodeCoord(1:3,n+FE_Nnodes(t),e) &
/ count(FE_subNodeParent(:,n,t) > 0)
enddo
enddo
return
endsubroutine
@ -3140,9 +3149,15 @@ endsubroutine
real(pReal), dimension(3,FE_NipFaceNodes) :: nPos ! coordinates of nodes on IP face
real(pReal), dimension(Ntriangles,FE_NipFaceNodes) :: volume ! volumes of possible tetrahedra
real(pReal), dimension(3) :: centerOfGravity
logical :: calcIPvolume = .false.
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems)) ; mesh_ipVolume = 0.0_pReal
allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems)) ; mesh_ipCenterOfGravity = 0.0_pReal
if (.not. allocated(mesh_ipVolume)) then
allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems))
mesh_ipVolume = 0.0_pReal
mesh_ipCenterOfGravity = 0.0_pReal
calcIPvolume = .true.
endif
do e = 1,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
@ -3168,21 +3183,23 @@ endsubroutine
endif
enddo
centerOfGravity = sum(gravityNodePos,2)/count(gravityNode)
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j
nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor
centerOfGravity)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
mesh_ipCenterOfGravity(:,i,e) = centerOfGravity
if (calcIPvolume) then
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP and add tetrahedra which connect to CoG
forall (n = 1:FE_NipFaceNodes) nPos(:,n) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
forall (n = 1:FE_NipFaceNodes, j = 1:Ntriangles) & ! start at each interface node and build valid triangles to cover interface
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1+mod(n-1 +j ,FE_NipFaceNodes)), & ! start at offset j
nPos(:,1+mod(n-1 +j+1,FE_NipFaceNodes)), & ! and take j's neighbor
centerOfGravity)
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) + sum(volume) ! add contribution from this interface
enddo
mesh_ipVolume(i,e) = mesh_ipVolume(i,e) / FE_NipFaceNodes ! renormalize with interfaceNodeNum due to loop over them
endif
enddo
enddo
return
endsubroutine
@ -3268,13 +3285,13 @@ do dir = 1,3 ! check periodicity in direction
minimumNodes = 0_pInt
maximumNodes = 0_pInt
minCoord = minval(mesh_node(dir,:))
maxCoord = maxval(mesh_node(dir,:))
minCoord = minval(mesh_node0(dir,:))
maxCoord = maxval(mesh_node0(dir,:))
do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes
if (abs(mesh_node(dir,node) - minCoord) <= tolerance) then
if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then
minimumNodes(1) = minimumNodes(1) + 1_pInt
minimumNodes(minimumNodes(1)+1) = node
elseif (abs(mesh_node(dir,node) - maxCoord) <= tolerance) then
elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then
maximumNodes(1) = maximumNodes(1) + 1_pInt
maximumNodes(maximumNodes(1)+1) = node
endif
@ -3290,7 +3307,7 @@ do dir = 1,3 ! check periodicity in direction
cycle
do n2 = 1,maximumNodes(1)
maximumNode = maximumNodes(n2+1)
distance = abs(mesh_node(:,minimumNode) - mesh_node(:,maximumNode))
distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode))
if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance)
mesh_nodeTwins(dir,minimumNode) = maximumNode
mesh_nodeTwins(dir,maximumNode) = minimumNode