diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 5ca5f187c..20668d740 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -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)) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 57ebaceed..18dc55800 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -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 diff --git a/code/mesh.f90 b/code/mesh.f90 index 6a0c0750a..871c04b05 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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,*)