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

Merge branch 'grid-mesh-cleanup' into 'development' 

Grid mesh cleanup

See merge request damask/DAMASK!80
This commit is contained in:
Franz Roters 2019-05-31 12:34:46 +02:00
parent 45efa9cc91 30826c5c86
commit 7042815542
7 changed files with 140 additions and 519 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/SpectralMethod/Polycrystal/material.config
View File

@ -11,7 +11,6 @@ mech none
[almostAll]
(output) phase
(output) texture
(output) volume
(output) orientation # quaternion
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4)
(output) f # deformation gradient tensor; synonyms: "defgrad"

1
src/commercialFEM_fileList.f90
View File

@ -14,6 +14,7 @@
#include "Lambert.f90"
#include "rotations.f90"
#include "FEsolving.f90"
#include "geometry_plastic_nonlocal.f90"
#include "element.f90"
#include "mesh_base.f90"
#ifdef Abaqus

10
src/crystallite.f90
View File

@ -93,7 +93,6 @@ module crystallite
enumerator :: undefined_ID, &
phase_ID, &
texture_ID, &
volume_ID, &
orientation_ID, &
grainrotation_ID, &
defgrad_ID, &
@ -286,8 +285,6 @@ subroutine crystallite_init
crystallite_outputID(o,c) = phase_ID
case ('texture') outputName
crystallite_outputID(o,c) = texture_ID
case ('volume') outputName
crystallite_outputID(o,c) = volume_ID
case ('orientation') outputName
crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName
@ -336,7 +333,7 @@ subroutine crystallite_init
do r = 1,size(config_crystallite)
do o = 1,crystallite_Noutput(r)
select case(crystallite_outputID(o,r))
case(phase_ID,texture_ID,volume_ID)
case(phase_ID,texture_ID)
mySize = 1
case(orientation_ID,grainrotation_ID)
mySize = 4
@ -914,11 +911,6 @@ function crystallite_postResults(ipc, ip, el)
case (texture_ID)
mySize = 1
crystallite_postResults(c+1) = real(material_texture(ipc,ip,el),pReal) ! textureID of grain
case (volume_ID)
mySize = 1
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
/ real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
case (orientation_ID)
mySize = 4
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion()

52
src/geometry_plastic_nonlocal.f90 Normal file
View File

@ -0,0 +1,52 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Geometric information about the IP cells needed for the nonlocal
! plasticity model
!--------------------------------------------------------------------------------------------------
module geometry_plastic_nonlocal
use prec
implicit none
private
logical, dimension(3), public, parameter :: &
geometry_plastic_nonlocal_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) NEEDED?
integer, dimension(:,:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
real(pReal), dimension(:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPvolume !< volume associated with IP (initially!)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IParea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable, public, protected :: &
geometry_plastic_nonlocal_IPareaNormal !< area normal of interface to neighboring IP (initially!)
public :: &
geometry_plastic_nonlocal_set_IPneighborhood, &
geometry_plastic_nonlocal_set_IPvolume
contains
subroutine geometry_plastic_nonlocal_set_IPneighborhood(IPneighborhood)
integer, dimension(:,:,:,:), intent(in) :: IPneighborhood
geometry_plastic_nonlocal_IPneighborhood = IPneighborhood
end subroutine geometry_plastic_nonlocal_set_IPneighborhood
subroutine geometry_plastic_nonlocal_set_IPvolume(IPvolume)
real(pReal), dimension(:,:), intent(in) :: IPvolume
geometry_plastic_nonlocal_IPvolume = IPvolume
end subroutine geometry_plastic_nonlocal_set_IPvolume
end module geometry_plastic_nonlocal

22
src/grid/spectral_utilities.f90
View File

@ -196,7 +196,6 @@ subroutine utilities_init
grid3Offset, &
geomSize
implicit none
PetscErrorCode :: ierr
integer :: i, j, k, &
FFTW_planner_flag
@ -425,7 +424,6 @@ subroutine utilities_updateGamma(C,saveReference)
math_det33, &
math_invert2
implicit none
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
logical , intent(in) :: saveReference !< save reference stiffness to file for restart
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
@ -473,7 +471,6 @@ end subroutine utilities_updateGamma
!> @details Does an unweighted filtered FFT transform from real to complex
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward
implicit none
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
@ -485,7 +482,6 @@ end subroutine utilities_FFTtensorForward
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
tensorField_real = tensorField_real * wgt ! normalize the result by number of elements
@ -497,7 +493,6 @@ end subroutine utilities_FFTtensorBackward
!> @details Does an unweighted filtered FFT transform from real to complex
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward
implicit none
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
@ -509,7 +504,6 @@ end subroutine utilities_FFTscalarForward
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarBackward
implicit none
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
scalarField_real = scalarField_real * wgt ! normalize the result by number of elements
@ -522,7 +516,6 @@ end subroutine utilities_FFTscalarBackward
!> @details Does an unweighted filtered FFT transform from real to complex.
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward
implicit none
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
@ -534,7 +527,6 @@ end subroutine utilities_FFTvectorForward
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
vectorField_real = vectorField_real * wgt ! normalize the result by number of elements
@ -554,7 +546,6 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
grid, &
grid3Offset
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
@ -615,7 +606,6 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
grid, &
grid3
implicit none
real(pReal), dimension(3,3), intent(in) :: D_ref
real(pReal), intent(in) :: mobility_ref, deltaT
complex(pReal) :: GreenOp_hat
@ -644,7 +634,6 @@ real(pReal) function utilities_divergenceRMS()
grid, &
grid3
implicit none
integer :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom
@ -694,7 +683,6 @@ real(pReal) function utilities_curlRMS()
grid, &
grid3
implicit none
integer :: i, j, k, l, ierr
complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom
@ -766,7 +754,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
math_rotate_forward33, &
math_invert2
implicit none
real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance
real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness
real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame
@ -861,7 +848,6 @@ subroutine utilities_fourierScalarGradient()
grid3, &
grid
implicit none
integer :: i, j, k
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -879,7 +865,6 @@ subroutine utilities_fourierVectorDivergence()
grid3, &
grid
implicit none
integer :: i, j, k
scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -898,7 +883,6 @@ subroutine utilities_fourierVectorGradient()
grid3, &
grid
implicit none
integer :: i, j, k, m, n
tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -919,7 +903,6 @@ subroutine utilities_fourierTensorDivergence()
grid3, &
grid
implicit none
integer :: i, j, k, m, n
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -954,7 +937,6 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
materialpoint_dPdF, &
materialpoint_stressAndItsTangent
implicit none
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
@ -1032,7 +1014,6 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
grid3, &
grid
implicit none
real(pReal), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon
real(pReal), intent(in) :: &
@ -1063,7 +1044,6 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
grid3, &
grid
implicit none
real(pReal), intent(in) :: &
timeinc !< timeinc of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
@ -1100,7 +1080,6 @@ pure function utilities_getFreqDerivative(k_s)
geomSize, &
grid
implicit none
integer, intent(in), dimension(3) :: k_s !< indices of frequency
complex(pReal), dimension(3) :: utilities_getFreqDerivative
@ -1158,7 +1137,6 @@ subroutine utilities_updateIPcoords(F)
grid3Offset, &
geomSize, &
mesh_ipCoordinates
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
integer :: i, j, k, m, ierr

562
src/mesh_grid.f90
View File

@ -7,18 +7,14 @@
!--------------------------------------------------------------------------------------------------
module mesh
use, intrinsic :: iso_c_binding
use prec, only: pReal, pInt
use prec
use geometry_plastic_nonlocal
use mesh_base
implicit none
private
integer(pInt), public, protected :: &
mesh_Nnodes, & !< total number of nodes in mesh
mesh_Ncellnodes, & !< total number of cell nodes in mesh (including duplicates)
mesh_Ncells, & !< total number of cells in mesh
mesh_maxNipNeighbors, & !< max number of IP neighbors in any CP element
mesh_maxNsharedElems !< max number of CP elements sharing a node
mesh_Nnodes
integer(pInt), dimension(:), allocatable, private :: &
microGlobal
@ -34,9 +30,9 @@ module mesh
real(pReal), public, protected :: &
mesh_unitlength !< physical length of one unit in mesh
real(pReal), dimension(:,:), allocatable, public :: &
mesh_node, & !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
mesh_cellnode !< cell node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, private :: &
mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
real(pReal), dimension(:,:), allocatable, public, protected :: &
mesh_ipVolume, & !< volume associated with IP (initially!)
@ -53,56 +49,8 @@ module mesh
logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes)
integer(pInt), dimension(:,:), allocatable, private :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
integer(pInt),dimension(:,:,:), allocatable, private :: &
mesh_cell !< cell connectivity for each element,ip/cell
integer(pInt), dimension(:,:,:), allocatable, private :: &
FE_cellface !< list of intra-cell cell node IDs that constitute the cell faces of a specific type of cell
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_"
integer(pInt), parameter, private :: &
FE_Ngeomtypes = 10_pInt, &
FE_Ncelltypes = 4_pInt, &
FE_maxNmatchingNodesPerFace = 4_pInt, &
FE_maxNfaces = 6_pInt, &
FE_maxNcellnodesPerCell = 8_pInt, &
FE_maxNcellfaces = 6_pInt, &
FE_maxNcellnodesPerCellface = 4_pInt
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCell = & !< number of cell nodes in a specific cell type
int([ &
3, & ! (2D 3node)
4, & ! (2D 4node)
4, & ! (3D 4node)
8 & ! (3D 8node)
],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type
int([&
2, & ! (2D 3node)
2, & ! (2D 4node)
3, & ! (3D 4node)
4 & ! (3D 8node)
],pInt)
integer(pInt), dimension(FE_Ncelltypes), parameter, private :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([&
3, & ! (2D 3node)
4, & ! (2D 4node)
4, & ! (3D 4node)
6 & ! (3D 8node)
],pInt)
! grid specific
integer(pInt), dimension(3), public, protected :: &
grid !< (global) grid
integer(pInt), public, protected :: &
@ -116,18 +64,14 @@ integer(pInt), dimension(:,:), allocatable, private :: &
size3offset !< (local) size offset in 3rd direction
public :: &
mesh_init, &
mesh_cellCenterCoordinates
mesh_init
private :: &
mesh_build_cellconnectivity, &
mesh_build_ipAreas, &
mesh_build_FEdata, &
mesh_build_ipNormals, &
mesh_spectral_build_nodes, &
mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, &
mesh_build_cellnodes, &
mesh_build_ipVolumes, &
mesh_build_ipCoordinates
type, public, extends(tMesh) :: tMesh_grid
@ -190,9 +134,8 @@ subroutine mesh_init(ip,el)
implicit none
include 'fftw3-mpi.f03'
integer(C_INTPTR_T) :: devNull, local_K, local_K_offset
integer :: ierr, worldsize
integer :: ierr, worldsize, j
integer(pInt), intent(in), optional :: el, ip
integer(pInt) :: j
logical :: myDebug
write(6,'(/,a)') ' <<<+- mesh init -+>>>'
@ -225,31 +168,31 @@ subroutine mesh_init(ip,el)
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
call mesh_spectral_build_nodes()
mesh_node0 = mesh_spectral_build_nodes()
mesh_node = mesh_node0
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call theMesh%init(mesh_node)
call theMesh%setNelems(product(grid(1:2))*grid3)
mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3) ! reallocate/shrink in case of MPI
mesh_maxNipNeighbors = theMesh%elem%nIPneighbors
call mesh_spectral_build_elements()
mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3Offset+1: &
product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI
if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
call mesh_build_FEdata ! get properties of the different types of elements
call mesh_build_cellconnectivity
if (myDebug) write(6,'(a)') ' Built cell connectivity'; flush(6)
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes)
if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6)
call mesh_build_ipCoordinates
mesh_ipCoordinates = mesh_build_ipCoordinates()
if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6)
call mesh_build_ipVolumes
allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3])))
if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6)
call mesh_build_ipAreas
mesh_ipArea = mesh_build_ipAreas()
mesh_ipAreaNormal = mesh_build_ipNormals()
if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6)
call mesh_spectral_build_ipNeighborhood
call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood)
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
@ -264,13 +207,10 @@ subroutine mesh_init(ip,el)
!!!! COMPATIBILITY HACK !!!!
! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
! hence, xxPerElem instead of maxXX
! better name
theMesh%homogenizationAt = mesh_element(3,:)
theMesh%microstructureAt = mesh_element(4,:)
!!!!!!!!!!!!!!!!!!!!!!!!
deallocate(mesh_cell)
end subroutine mesh_init
@ -394,7 +334,7 @@ subroutine mesh_spectral_read_grid()
allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant)
!--------------------------------------------------------------------------------------------------
! read and interprete content
! read and interpret content
e = 1_pInt
do while (startPos < len(rawData))
endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt
@ -429,44 +369,53 @@ subroutine mesh_spectral_read_grid()
end subroutine mesh_spectral_read_grid
!--------------------------------------------------------------------------------------------------
!> @brief Store x,y,z coordinates of all nodes in mesh.
!! Allocates global arrays 'mesh_node0' and 'mesh_node'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_nodes()
!---------------------------------------------------------------------------------------------------
!> @brief Calculates position of nodes (pretend to be an element)
!---------------------------------------------------------------------------------------------------
pure function mesh_spectral_build_nodes()
implicit none
integer(pInt) :: n
real(pReal), dimension(3,mesh_Nnodes) :: mesh_spectral_build_nodes
integer :: n,a,b,c
allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal)
n = 0
do c = 0, grid3
do b = 0, grid(2)
do a = 0, grid(1)
n = n + 1
mesh_spectral_build_nodes(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
enddo
enddo
enddo
forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = mesh_unitlength * &
geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) &
/ real(grid(1),pReal)
mesh_node0(2,n+1_pInt) = mesh_unitlength * &
geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) &
/ real(grid(2),pReal)
mesh_node0(3,n+1_pInt) = mesh_unitlength * &
size3*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid3+1_pInt)),pReal) &
/ real(grid3,pReal) + &
size3offset
end forall
end function mesh_spectral_build_nodes
mesh_node = mesh_node0
end subroutine mesh_spectral_build_nodes
!---------------------------------------------------------------------------------------------------
!> @brief Calculates position of IPs/cell centres (pretend to be an element)
!---------------------------------------------------------------------------------------------------
function mesh_build_ipCoordinates()
real(pReal), dimension(3,1,theMesh%nElems) :: mesh_build_ipCoordinates
integer :: n,a,b,c
n = 0
do c = 1, grid3
do b = 1, grid(2)
do a = 1, grid(1)
n = n + 1
mesh_build_ipCoordinates(1:3,1,n) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
enddo
enddo
enddo
end function mesh_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Store FEid, type, material, texture, and node list per element.
!! Allocates global array 'mesh_element'
!> @todo does the IO_error makes sense?
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_elements()
use IO, only: &
IO_error
implicit none
integer(pInt) :: &
e, &
elemOffset
@ -475,11 +424,9 @@ subroutine mesh_spectral_build_elements()
allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt)
elemOffset = product(grid(1:2))*grid3Offset
e = 0_pInt
do while (e < theMesh%nElems) ! fill expected number of elements, stop at end of data
e = e+1_pInt ! valid element entry
do e=1, theMesh%nElems
mesh_element( 1,e) = -1_pInt ! DEPRECATED
mesh_element( 2,e) = 10_pInt
mesh_element( 2,e) = -1_pInt ! DEPRECATED
mesh_element( 3,e) = mesh_homogenizationAt(e)
mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &
@ -493,8 +440,6 @@ subroutine mesh_spectral_build_elements()
mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt
enddo
if (e /= theMesh%nElems) call IO_error(880_pInt,e)
end subroutine mesh_spectral_build_elements
@ -508,7 +453,7 @@ subroutine mesh_spectral_build_ipNeighborhood
integer(pInt) :: &
x,y,z, &
e
allocate(mesh_ipNeighborhood(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems),source=0_pInt)
allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt)
e = 0_pInt
do z = 0_pInt,grid3-1_pInt
@ -562,7 +507,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
debug_level, &
debug_levelBasic
implicit none
real(pReal), intent(in), dimension(:,:,:,:) :: &
centres
real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: &
@ -641,385 +585,35 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
end function mesh_nodesAroundCentres
!#################################################################################################################
!#################################################################################################################
!#################################################################################################################
! The following routines are not solver specific and should be included in mesh_base (most likely in modified form)
!#################################################################################################################
!#################################################################################################################
!#################################################################################################################
!--------------------------------------------------------------------------------------------------
!> @brief Split CP elements into cells.
!> @details Build a mapping between cells and the corresponding cell nodes ('mesh_cell').
!> Cell nodes that are also matching nodes are unique in the list of cell nodes,
!> all others (currently) might be stored more than once.
!> Also allocates the 'mesh_node' array.
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_cellconnectivity
pure function mesh_build_ipAreas()
implicit none
integer(pInt), dimension(:), allocatable :: &
matchingNode2cellnode
integer(pInt), dimension(:,:), allocatable :: &
cellnodeParent
integer(pInt), dimension(theMesh%elem%Ncellnodes) :: &
localCellnode2globalCellnode
integer(pInt) :: &
e,n,i, &
matchingNodeID, &
localCellnodeID
integer(pInt), dimension(FE_Ngeomtypes), parameter :: FE_NmatchingNodes = & !< number of nodes that are needed for face matching in a specific type of element geometry
int([ &
3, & ! element 6 (2D 3node 1ip)
3, & ! element 125 (2D 6node 3ip)
4, & ! element 11 (2D 4node 4ip)
4, & ! element 27 (2D 8node 9ip)
4, & ! element 134 (3D 4node 1ip)
4, & ! element 127 (3D 10node 4ip)
6, & ! element 136 (3D 6node 6ip)
8, & ! element 117 (3D 8node 1ip)
8, & ! element 7 (3D 8node 8ip)
8 & ! element 21 (3D 20node 27ip)
],pInt)
real(pReal), dimension(6,1,theMesh%nElems) :: mesh_build_ipAreas
allocate(mesh_cell(FE_maxNcellnodesPerCell,theMesh%elem%nIPs,theMesh%nElems), source=0_pInt)
allocate(matchingNode2cellnode(theMesh%nNodes), source=0_pInt)
allocate(cellnodeParent(2_pInt,theMesh%elem%Ncellnodes*theMesh%nElems), source=0_pInt)
mesh_Ncells = theMesh%nElems*theMesh%elem%nIPs
!--------------------------------------------------------------------------------------------------
! Count cell nodes (including duplicates) and generate cell connectivity list
mesh_Ncellnodes = 0_pInt
do e = 1_pInt,theMesh%nElems
localCellnode2globalCellnode = 0_pInt
do i = 1_pInt,theMesh%elem%nIPs
do n = 1_pInt,theMesh%elem%NcellnodesPerCell
localCellnodeID = theMesh%elem%cell(n,i)
if (localCellnodeID <= FE_NmatchingNodes(theMesh%elem%geomType)) then ! this cell node is a matching node
matchingNodeID = mesh_element(4_pInt+localCellnodeID,e)
if (matchingNode2cellnode(matchingNodeID) == 0_pInt) then ! if this matching node does not yet exist in the glbal cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
matchingNode2cellnode(matchingNodeID) = mesh_Ncellnodes ! ... and remember its global ID
cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and where it belongs to
cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID
endif
mesh_cell(n,i,e) = matchingNode2cellnode(matchingNodeID)
else ! this cell node is no matching node
if (localCellnode2globalCellnode(localCellnodeID) == 0_pInt) then ! if this local cell node does not yet exist in the global cell node list ...
mesh_Ncellnodes = mesh_Ncellnodes + 1_pInt ! ... count it as cell node ...
localCellnode2globalCellnode(localCellnodeID) = mesh_Ncellnodes ! ... and remember its global ID ...
cellnodeParent(1_pInt,mesh_Ncellnodes) = e ! ... and it belongs to
cellnodeParent(2_pInt,mesh_Ncellnodes) = localCellnodeID
endif
mesh_cell(n,i,e) = localCellnode2globalCellnode(localCellnodeID)
endif
enddo
enddo
enddo
allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes))
allocate(mesh_cellnode(3_pInt,mesh_Ncellnodes))
forall(n = 1_pInt:mesh_Ncellnodes)
mesh_cellnodeParent(1,n) = cellnodeParent(1,n)
mesh_cellnodeParent(2,n) = cellnodeParent(2,n)
endforall
end subroutine mesh_build_cellconnectivity
!--------------------------------------------------------------------------------------------------
!> @brief Calculate position of cellnodes from the given position of nodes
!> Build list of cellnodes' coordinates.
!> Cellnode coordinates are calculated from a weighted sum of node coordinates.
!--------------------------------------------------------------------------------------------------
function mesh_build_cellnodes(nodes,Ncellnodes)
implicit none
integer(pInt), intent(in) :: Ncellnodes !< requested number of cellnodes
real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes
real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes
integer(pInt) :: &
e,n,m, &
localCellnodeID
real(pReal), dimension(3) :: &
myCoords
mesh_build_cellnodes = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(e,localCellnodeID,myCoords)
do n = 1_pInt,Ncellnodes ! loop over cell nodes
e = mesh_cellnodeParent(1,n)
localCellnodeID = mesh_cellnodeParent(2,n)
myCoords = 0.0_pReal
do m = 1_pInt,theMesh%elem%nNodes
myCoords = myCoords + nodes(1:3,mesh_element(4_pInt+m,e)) &
* theMesh%elem%cellNodeParentNodeWeights(m,localCellnodeID)
enddo
mesh_build_cellnodes(1:3,n) = myCoords / sum(theMesh%elem%cellNodeParentNodeWeights(:,localCellnodeID))
enddo
!$OMP END PARALLEL DO
end function mesh_build_cellnodes
!--------------------------------------------------------------------------------------------------
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of one in order to calculate the volume.
!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
!> shape with a cell face as basis and the central ip at the tip. This subvolume is
!> calculated as an average of four tetrahedals with three corners on the cell face
!> and one corner at the central ip.
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipVolumes
use math, only: &
math_volTetrahedron, &
math_areaTriangle
implicit none
integer(pInt) :: e,t,g,c,i,m,f,n
real(pReal), dimension(FE_maxNcellnodesPerCellface,FE_maxNcellfaces) :: subvolume
allocate(mesh_ipVolume(theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal)
!$OMP PARALLEL DO PRIVATE(t,g,c,m,subvolume)
do e = 1_pInt,theMesh%nElems ! loop over cpElems
select case (theMesh%elem%cellType)
case (1_pInt) ! 2D 3node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), &
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e)))
case (2_pInt) ! 2D 4node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_areaTriangle(mesh_cellnode(1:3,mesh_cell(1,i,e)), & ! here we assume a planar shape, so division in two triangles suffices
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e))) &
+ math_areaTriangle(mesh_cellnode(1:3,mesh_cell(3,i,e)), &
mesh_cellnode(1:3,mesh_cell(4,i,e)), &
mesh_cellnode(1:3,mesh_cell(1,i,e)))
case (3_pInt) ! 3D 4node
forall (i = 1_pInt:theMesh%elem%nIPs) & ! loop over ips=cells in this element
mesh_ipVolume(i,e) = math_volTetrahedron(mesh_cellnode(1:3,mesh_cell(1,i,e)), &
mesh_cellnode(1:3,mesh_cell(2,i,e)), &
mesh_cellnode(1:3,mesh_cell(3,i,e)), &
mesh_cellnode(1:3,mesh_cell(4,i,e)))
case (4_pInt)
c = theMesh%elem%cellType ! 3D 8node
m = FE_NcellnodesPerCellface(c)
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
subvolume = 0.0_pReal
forall(f = 1_pInt:FE_NipNeighbors(c), n = 1_pInt:FE_NcellnodesPerCellface(c)) &
subvolume(n,f) = math_volTetrahedron(&
mesh_cellnode(1:3,mesh_cell(FE_cellface( n ,f,c),i,e)), &
mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n ,m),f,c),i,e)), &
mesh_cellnode(1:3,mesh_cell(FE_cellface(1+mod(n+1,m),f,c),i,e)), &
mesh_ipCoordinates(1:3,i,e))
mesh_ipVolume(i,e) = 0.5_pReal * sum(subvolume) ! each subvolume is based on four tetrahedrons, altough the face consists of only two triangles -> averaging factor two
enddo
end select
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipVolumes
!--------------------------------------------------------------------------------------------------
!> @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 CENTER 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
implicit none
integer(pInt) :: e,c,i,n
real(pReal), dimension(3) :: myCoords
if (.not. allocated(mesh_ipCoordinates)) &
allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal)
!$OMP PARALLEL DO PRIVATE(c,myCoords)
do e = 1_pInt,theMesh%nElems ! loop over cpElems
c = theMesh%elem%cellType
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
myCoords = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
enddo
mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal)
enddo
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
!> @brief Calculates cell center coordinates.
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
integer(pInt) :: c,n
c = theMesh%elem%cellType
mesh_cellCenterCoordinates = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal)
end function mesh_cellCenterCoordinates
mesh_build_ipAreas(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3))
mesh_build_ipAreas(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1))
mesh_build_ipAreas(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2))
end function mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipAreas
use math, only: &
math_cross
pure function mesh_build_ipNormals()
implicit none
integer(pInt) :: e,t,g,c,i,f,n,m
real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals
real(pReal), dimension(3) :: normal
real, dimension(3,6,1,theMesh%nElems) :: mesh_build_ipNormals
allocate(mesh_ipArea(theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
allocate(mesh_ipAreaNormal(3_pInt,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal)
!$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals)
do e = 1_pInt,theMesh%nElems ! loop over cpElems
c = theMesh%elem%cellType
select case (c)
case (1_pInt,2_pInt) ! 2D 3 or 4 node
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
normal(1) = nodePos(2,2) - nodePos(2,1) ! x_normal = y_connectingVector
normal(2) = -(nodePos(1,2) - nodePos(1,1)) ! y_normal = -x_connectingVector
normal(3) = 0.0_pReal
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
enddo
enddo
case (3_pInt) ! 3D 4node
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
nodePos(1:3,3) - nodePos(1:3,1))
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal) ! ensure unit length of area normal
enddo
enddo
case (4_pInt) ! 3D 8node
! for this cell type we get the normal of the quadrilateral face as an average of
! four normals of triangular subfaces; since the face consists only of two triangles,
! the sum has to be divided by two; this whole prcedure tries to compensate for
! probable non-planar cell surfaces
m = FE_NcellnodesPerCellface(c)
do i = 1_pInt,theMesh%elem%nIPs ! loop over ips=cells in this element
do f = 1_pInt,FE_NipNeighbors(c) ! loop over cell faces
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
normals(1:3,n) = 0.5_pReal &
* math_cross(nodePos(1:3,1+mod(n ,m)) - nodePos(1:3,n), &
nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n))
normal = 0.5_pReal * sum(normals,2)
mesh_ipArea(f,i,e) = norm2(normal)
mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal)
enddo
enddo
end select
enddo
!$OMP END PARALLEL DO
end subroutine mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
!> @brief get properties of different types of finite elements
!> @details assign globals: FE_nodesAtIP, FE_ipNeighbor, FE_subNodeOnIPFace
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_FEdata
implicit none
integer(pInt) :: me
allocate(FE_cellface(FE_maxNcellnodesPerCellface,FE_maxNcellfaces,FE_Ncelltypes), source=0_pInt)
! *** FE_cellface ***
me = 0_pInt
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 3node, VTK_TRIANGLE (5)
reshape(int([&
2,3, &
3,1, &
1,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 2D 4node, VTK_QUAD (9)
reshape(int([&
2,3, &
4,1, &
3,4, &
1,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 4node, VTK_TETRA (10)
reshape(int([&
1,3,2, &
1,2,4, &
2,3,4, &
1,4,3 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
me = me + 1_pInt
FE_cellface(1:FE_NcellnodesPerCellface(me),1:FE_NipNeighbors(me),me) = & ! 3D 8node, VTK_HEXAHEDRON (12)
reshape(int([&
2,3,7,6, &
4,1,5,8, &
3,4,8,7, &
1,2,6,5, &
5,6,7,8, &
1,4,3,2 &
],pInt),[FE_NcellnodesPerCellface(me),FE_NipNeighbors(me)])
mesh_build_ipNormals(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,theMesh%nElems)
mesh_build_ipNormals(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,theMesh%nElems)
end subroutine mesh_build_FEdata
end function mesh_build_ipNormals
end module mesh

11
src/plastic_nonlocal.f90
View File

@ -5,10 +5,15 @@
!> @brief material subroutine for plasticity including dislocation flux
!--------------------------------------------------------------------------------------------------
module plastic_nonlocal
use prec, only: &
pReal
use prec
use future
use geometry_plastic_nonlocal, only: &
periodicSurface => geometry_plastic_nonlocal_periodicSurface, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
IPvolume => geometry_plastic_nonlocal_IPvolume, &
IParea => geometry_plastic_nonlocal_IParea, &
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal
implicit none
private
real(pReal), parameter, private :: &