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

introduced distinction between ip coordinates and cell center coordinates; the former can be globally calculated by subroutine mesh_build_ipCoordinates, the latter locally by the function mesh_cellCenterCoordinates; renamed mesh_ipCenterOfGravity to mesh_ipCordinates

This commit is contained in:
Christoph Kords 2012-11-06 14:37:13 +00:00
parent bb033c5fe7
commit 5b6baa7c0d
3 changed files with 79 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/CPFEM.f90
View File

@ -279,7 +279,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
mesh_element, &
mesh_node0, &
mesh_node, &
mesh_ipCenterOfGravity, &
mesh_ipCoordinates, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates, &
@ -612,7 +612,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_calc_done = .false.
#ifndef Marc
mesh_ipCenterOfGravity(1:3,IP,cp_en) = numerics_unitlength * coords(1:3,1)
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))

24
code/constitutive_nonlocal.f90
View File

@ -1180,7 +1180,7 @@ use mesh, only: mesh_NcpElems, &
FE_NipNeighbors, &
FE_maxNipNeighbors, &
mesh_ipNeighborhood, &
mesh_ipCenterOfGravity, &
mesh_ipCoordinates, &
mesh_ipVolume, &
mesh_ipAreaNormal
use material, only: homogenization_maxNgrains, &
@ -1338,7 +1338,7 @@ tauBack = 0.0_pReal
if (.not. phase_localPlasticity(phase) .and. constitutive_nonlocal_shortRangeStressCorrection(instance)) then
call math_invert33(Fe, invFe, detFe, inversionError)
call math_invert33(Fp, invFp, detFp, inversionError)
ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)
ipCoords = mesh_ipCoordinates(1:3,ip,el)
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4)
FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal)
@ -1355,7 +1355,7 @@ if (.not. phase_localPlasticity(phase) .and. constitutive_nonlocal_shortRangeStr
neighboring_instance = phase_plasticityInstance(neighboring_phase)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el)
neighboring_ipCoords = mesh_ipCoordinates(1:3,neighboring_ip,neighboring_el)
if (.not. phase_localPlasticity(neighboring_phase) &
.and. neighboring_latticeStruct == latticeStruct &
.and. neighboring_instance == instance) then
@ -2767,7 +2767,7 @@ use mesh, only: mesh_NcpElems, &
mesh_element, &
mesh_node0, &
FE_Nips, &
mesh_ipCenterOfGravity, &
mesh_cellCenterCoordinates, &
mesh_ipVolume, &
mesh_periodicSurface
use material, only: homogenization_maxNgrains, &
@ -2827,8 +2827,8 @@ real(pReal), dimension(3) :: connection, & ! connection vecto
connection_neighboringSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor
maxCoord, minCoord, &
meshSize, &
ipCoords, &
neighboring_ipCoords
coords, & ! x,y,z coordinates of cell center of ip volume
neighboring_coords ! x,y,z coordinates of cell center of neighboring ip volume
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in neighboring material point's slip system frame
Tdislo_neighboringLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighboring material point
invFe, & ! inverse of my elastic deformation gradient
@ -2876,12 +2876,12 @@ if (.not. phase_localPlasticity(phase)) then
minCoord(dir) = minval(mesh_node0(dir,:))
enddo
meshSize = maxCoord - minCoord
ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)
coords = mesh_cellCenterCoordinates(ip,el)
periodicImages = 0_pInt
do dir = 1_pInt,3_pInt
if (mesh_periodicSurface(dir)) then
periodicImages(1,dir) = floor((ipCoords(dir) - constitutive_nonlocal_R(instance) - minCoord(dir)) / meshSize(dir), pInt)
periodicImages(2,dir) = ceiling((ipCoords(dir) + constitutive_nonlocal_R(instance) - maxCoord(dir)) / meshSize(dir), pInt)
periodicImages(1,dir) = floor((coords(dir) - constitutive_nonlocal_R(instance) - minCoord(dir)) / meshSize(dir), pInt)
periodicImages(2,dir) = ceiling((coords(dir) + constitutive_nonlocal_R(instance) - maxCoord(dir)) / meshSize(dir), pInt)
endif
enddo
@ -2921,9 +2921,9 @@ ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el))
if (neighboring_el /= el .or. neighboring_ip /= ip &
.or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then
neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) &
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
connection = neighboring_ipCoords - ipCoords
neighboring_coords = mesh_cellCenterCoordinates(neighboring_ip,neighboring_el) &
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
connection = neighboring_coords - coords
distance = sqrt(sum(connection * connection))
if (distance > constitutive_nonlocal_R(instance)) then
cycle

73
code/mesh.f90
View File

@ -60,7 +60,7 @@ module mesh
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:,:), allocatable, public :: &
mesh_ipCenterOfGravity, & !< center of gravity of IP (after deformation!)
mesh_ipCoordinates, & !< IP x,y,z coordinates (after deformation!)
mesh_ipArea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public :: &
@ -279,7 +279,8 @@ module mesh
mesh_FEasCP, &
mesh_build_subNodeCoords, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates
mesh_build_ipCoordinates, &
mesh_cellCenterCoordinates
#ifdef Spectral
public :: mesh_regrid, &
mesh_regular_grid, &
@ -384,7 +385,7 @@ subroutine mesh_init(ip,element)
if (allocated(mesh_node)) deallocate(mesh_node)
if (allocated(mesh_element)) deallocate(mesh_element)
if (allocated(mesh_subNodeCoord)) deallocate(mesh_subNodeCoord)
if (allocated(mesh_ipCenterOfGravity)) deallocate(mesh_ipCenterOfGravity)
if (allocated(mesh_ipCoordinates)) deallocate(mesh_ipCoordinates)
if (allocated(mesh_ipArea)) deallocate(mesh_ipArea)
if (allocated(mesh_ipAreaNormal)) deallocate(mesh_ipAreaNormal)
if (allocated(mesh_sharedElem)) deallocate(mesh_sharedElem)
@ -589,7 +590,7 @@ subroutine mesh_build_ipVolumes
volume(j,n) = math_volTetrahedron(nPos(:,n), & ! calc volume of respective tetrahedron to CoG
nPos(:,1_pInt+mod(n-1_pInt +j ,FE_NipFaceNodes)),& ! start at offset j
nPos(:,1_pInt+mod(n-1_pInt +j+1_pInt,FE_NipFaceNodes)),& ! and take j's neighbor
mesh_ipCenterOfGravity(:,i,e))
mesh_cellCenterCoordinates(i,e))
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
@ -600,7 +601,16 @@ end subroutine mesh_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCenterOfGravity'
!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
! Called by all solvers in mesh_init in order to initialize the ip coordinates.
! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
! so in this case the ip coordinates are always calculated on the basis of this subroutine.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CELLENTER AND NOT THE IP COORDINATES,
! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME.
! HAS TO BE CHANGED IN A LATER VERSION.
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipCoordinates
@ -612,7 +622,7 @@ subroutine mesh_build_ipCoordinates
real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav
real(pReal), dimension(3) :: centerOfGravity
if (.not. allocated(mesh_ipCenterOfGravity)) allocate(mesh_ipCenterOfGravity(3,mesh_maxNips,mesh_NcpElems))
if (.not. allocated(mesh_ipCoordinates)) allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems))
do e = 1_pInt,mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get elemType
@ -638,13 +648,60 @@ subroutine mesh_build_ipCoordinates
endif
enddo
centerOfGravity = sum(gravityNodePos,2)/real(count(gravityNode),pReal)
mesh_ipCenterOfGravity(:,i,e) = centerOfGravity
mesh_ipCoordinates(:,i,e) = centerOfGravity
enddo
enddo
end subroutine mesh_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(i,e)
use prec, only: tol_gravityNodePos
implicit none
!*** input variables
integer(pInt), intent(in) :: e, & ! element number
i ! integration point number
!*** output variables
real(pReal), dimension(3) :: mesh_cellCenterCoordinates ! x,y,z coordinates of the cell center of the requested IP cell
!*** local variables
integer(pInt) :: f,t,j,k,n
logical, dimension(mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNode ! flagList to find subnodes determining center of grav
real(pReal), dimension(3,mesh_maxNnodes+mesh_maxNsubNodes) :: gravityNodePos ! coordinates of subnodes determining center of grav
t = mesh_element(2,e) ! get elemType
gravityNode = .false. ! reset flagList
gravityNodePos = 0.0_pReal ! reset coordinates
do f = 1_pInt,FE_NipNeighbors(t) ! loop over interfaces of IP
do n = 1_pInt,FE_NipFaceNodes ! loop over nodes on interface
gravityNode(FE_subNodeOnIPFace(n,f,i,t)) = .true.
gravityNodePos(:,FE_subNodeOnIPFace(n,f,i,t)) = mesh_subNodeCoord(:,FE_subNodeOnIPFace(n,f,i,t),e)
enddo
enddo
do j = 1_pInt,mesh_maxNnodes+mesh_maxNsubNodes-1_pInt ! walk through entire flagList except last
if (gravityNode(j)) then ! valid node index
do k = j+1_pInt,mesh_maxNnodes+mesh_maxNsubNodes ! walk through remainder of list
if (gravityNode(k) .and. all(abs(gravityNodePos(:,j) - gravityNodePos(:,k)) < tol_gravityNodePos)) then ! found duplicate
gravityNode(j) = .false. ! delete first instance
gravityNodePos(:,j) = 0.0_pReal
exit ! continue with next suspect
endif
enddo
endif
enddo
mesh_cellCenterCoordinates = sum(gravityNodePos,2)/real(count(gravityNode),pReal)
endfunction mesh_cellCenterCoordinates
#ifdef Spectral
!--------------------------------------------------------------------------------------------------
!> @brief Reads resolution information from geometry file. If fileUnit is given,
@ -3709,7 +3766,7 @@ enddo
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_e /= e) cycle
do i = 1_pInt,FE_Nips(mesh_element(2,e))
if (iand(myDebug,debug_levelSelective) /= 0_pInt .and. debug_i /= i) cycle
write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCoordinates(:,i,e)
enddo
enddo
write (6,*)