write out and calculate cell and element connectivity

src/CPFEM.f90

@@ -5,9 +5,21 @@
 !--------------------------------------------------------------------------------------------------
 module CPFEM
  use prec
+ use numerics
+ use debug
+ use FEsolving
+ use math
+ use mesh
+ use material
+ use config
+ use crystallite
+ use homogenization
+ use IO
+ use DAMASK_interface
 
  implicit none
  private
+ 
  real(pReal),                      parameter,   private :: &
    CPFEM_odd_stress    = 1e15_pReal, &                                                              !< return value for stress in case of ping pong dummy cycle
    CPFEM_odd_jacobian  = 1e50_pReal                                                                 !< return value for jacobian in case of ping pong dummy cycle
@@ -238,86 +250,6 @@ end subroutine CPFEM_init
 !> @brief perform initialization at first call, update variables and call the actual material model
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
- use numerics, only: &
-   defgradTolerance, &
-   iJacoStiffness
- use debug, only: &
-   debug_level, &
-   debug_CPFEM, &
-   debug_levelBasic, &
-   debug_levelExtensive, &
-   debug_levelSelective, &
-   debug_stressMaxLocation, &
-   debug_stressMinLocation, &
-   debug_jacobianMaxLocation, &
-   debug_jacobianMinLocation, &
-   debug_stressMax, &
-   debug_stressMin, &
-   debug_jacobianMax, &
-   debug_jacobianMin, &
-   debug_e, &
-   debug_i
- use FEsolving, only: &
-   terminallyIll, &
-   FEsolving_execElem, &
-   FEsolving_execIP, &
-   restartWrite
- use math, only: &
-   math_identity2nd, &
-   math_det33, &
-   math_delta, &
-   math_sym3333to66, &
-   math_66toSym3333, &
-   math_sym33to6, &
-   math_6toSym33
- use mesh, only: &
-   mesh_FEasCP, &
-   theMesh, &
-   mesh_element
- use material, only: &
-   microstructure_elemhomo, &
-   plasticState, &
-   sourceState, &
-   homogState, &
-   thermalState, &
-   damageState, &
-   phaseAt, phasememberAt, &
-   material_phase, &
-   phase_plasticity, &
-   temperature, &
-   thermalMapping, &
-   thermal_type, &
-   THERMAL_conduction_ID, &
-   phase_Nsources, &
-   material_homogenizationAt
- use config, only: &
-   material_Nhomogenization
- use crystallite, only: &
-   crystallite_partionedF,&
-   crystallite_F0, &
-   crystallite_Fp0, &
-   crystallite_Fp, &
-   crystallite_Fi0, &
-   crystallite_Fi, &
-   crystallite_Lp0, &
-   crystallite_Lp, &
-   crystallite_Li0, &
-   crystallite_Li, &
-   crystallite_dPdF, &
-   crystallite_S0, &
-   crystallite_S
- use homogenization, only: &  
-   materialpoint_F, &
-   materialpoint_F0, &
-   materialpoint_P, &
-   materialpoint_dPdF, &
-   materialpoint_results, &
-   materialpoint_sizeResults, &
-   materialpoint_stressAndItsTangent, &
-   materialpoint_postResults
- use IO, only: &
-   IO_warning
- use DAMASK_interface
 
  integer(pInt), intent(in) ::                        elFE, &                                        !< FE element number
                                                      ip                                             !< integration point number
@@ -520,15 +452,12 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
      if (.not. parallelExecution) then
        FEsolving_execElem(1)     = elCP
        FEsolving_execElem(2)     = elCP
-       if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. &                               ! calculate unless homogeneous
+       FEsolving_execIP(1,elCP) = ip
-                (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then            ! and then only first ip
+       FEsolving_execIP(2,elCP) = ip
-         FEsolving_execIP(1,elCP) = ip
+       if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /=  0_pInt) &
-         FEsolving_execIP(2,elCP) = ip
+         write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
-         if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /=  0_pInt) &
+       call materialpoint_stressAndItsTangent(updateJaco, dt)                                     ! calculate stress and its tangent
-           write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
+       call materialpoint_postResults()
-         call materialpoint_stressAndItsTangent(updateJaco, dt)                                     ! calculate stress and its tangent
-         call materialpoint_postResults()
-       endif
 
      !* parallel computation and calulation not yet done
 
@@ -551,13 +480,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
 
      else terminalIllness
 
-       if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then                    ! me homogenous? --> copy from first ip
-         materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP)
-         materialpoint_F(1:3,1:3,ip,elCP) = materialpoint_F(1:3,1:3,1,elCP)
-         materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,elCP) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,elCP)
-         materialpoint_results(1:materialpoint_sizeResults,ip,elCP) = &
-         materialpoint_results(1:materialpoint_sizeResults,1,elCP)
-       endif
 
        ! translate from P to CS
        Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))

src/geometry_plastic_nonlocal.f90

@@ -10,8 +10,6 @@ module geometry_plastic_nonlocal
 
   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]

src/material.f90

@@ -152,7 +152,6 @@ module material
 
  logical, dimension(:), allocatable, public, protected :: &
    microstructure_active, &
-   microstructure_elemhomo, &                                                                       !< flag to indicate homogeneous microstructure distribution over element's IPs
    phase_localPlasticity                                                                            !< flags phases with local constitutive law
 
  integer, private :: &
@@ -323,12 +322,11 @@ subroutine material_init
    do h = 1,size(config_homogenization)
      write(6,'(1x,a32,1x,a16,1x,i6)') homogenization_name(h),homogenization_type(h),homogenization_Ngrains(h)
    enddo
-   write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents','homogeneous'
+   write(6,'(/,a14,18x,1x,a11,1x,a12,1x,a13)') 'microstructure','crystallite','constituents'
    do m = 1,size(config_microstructure)
-     write(6,'(1x,a32,1x,i11,1x,i12,1x,l13)') microstructure_name(m), &
+     write(6,'(1x,a32,1x,i11,1x,i12)') microstructure_name(m), &
-                                        microstructure_crystallite(m), &
+                                       microstructure_crystallite(m), &
-                                        microstructure_Nconstituents(m), &
+                                       microstructure_Nconstituents(m)
-                                        microstructure_elemhomo(m)
      if (microstructure_Nconstituents(m) > 0) then
        do c = 1,microstructure_Nconstituents(m)
          write(6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(c,m)),&
@@ -538,7 +536,6 @@ subroutine material_parseMicrostructure
  allocate(microstructure_crystallite(size(config_microstructure)),          source=0)
  allocate(microstructure_Nconstituents(size(config_microstructure)),        source=0)
  allocate(microstructure_active(size(config_microstructure)),               source=.false.)
- allocate(microstructure_elemhomo(size(config_microstructure)),             source=.false.)
 
  if(any(theMesh%microstructureAt > size(config_microstructure))) &
   call IO_error(155,ext_msg='More microstructures in geometry than sections in material.config')
@@ -549,7 +546,6 @@ subroutine material_parseMicrostructure
  do m=1, size(config_microstructure)
    microstructure_Nconstituents(m) =  config_microstructure(m)%countKeys('(constituent)')
    microstructure_crystallite(m)   =  config_microstructure(m)%getInt('crystallite')
-   microstructure_elemhomo(m)      =  config_microstructure(m)%keyExists('/elementhomogeneous/')
  enddo
 
  microstructure_maxNconstituents = maxval(microstructure_Nconstituents)

src/mesh_marc.f90

@@ -229,7 +229,6 @@ integer, dimension(:,:), allocatable, private :: &
 
 
  private :: &
-   mesh_get_damaskOptions, &
    mesh_build_cellconnectivity, &
    mesh_build_ipAreas, &
    FE_mapElemtype, &
@@ -341,8 +340,6 @@ subroutine mesh_init(ip,el)
  
  call mesh_marc_build_elements(FILEUNIT)
  if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
- call mesh_get_damaskOptions(mesh_periodicSurface,FILEUNIT)
- if (myDebug) write(6,'(a)') ' Got DAMASK options'; flush(6)
  close (FILEUNIT)
   
  
@@ -808,8 +805,8 @@ subroutine mesh_marc_build_elements(fileUnit)
      exit
    endif
  enddo
-
 620 rewind(fileUnit)                                                                                ! just in case "initial state" appears before "connectivity"
+ call calcCells(theMesh,mesh_element(5:,:))
  read (fileUnit,'(A300)',END=630) line
  do
    chunkPos = IO_stringPos(line)
@@ -847,45 +844,6 @@ subroutine mesh_marc_build_elements(fileUnit)
 630 end subroutine mesh_marc_build_elements
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief get any additional damask options from input file, sets mesh_periodicSurface
-!--------------------------------------------------------------------------------------------------
-subroutine mesh_get_damaskOptions(periodic_surface,fileUnit)
-
- integer, intent(in) :: fileUnit
-
- integer, allocatable, dimension(:) :: chunkPos
- character(len=300) :: line
- integer :: myStat
- integer :: chunk, Nchunks
- character(len=300) ::  v
- logical, dimension(3) :: periodic_surface
- 
-
- periodic_surface = .false.
- myStat = 0
- rewind(fileUnit)
- do while(myStat == 0)
-   read (fileUnit,'(a300)',iostat=myStat) line
-   chunkPos = IO_stringPos(line)
-   Nchunks = chunkPos(1)
-   if (IO_lc(IO_stringValue(line,chunkPos,1)) == '$damask' .and. Nchunks > 1) then        ! found keyword for damask option and there is at least one more chunk to read
-     select case(IO_lc(IO_stringValue(line,chunkPos,2)))
-       case('periodic')                                                                             ! damask Option that allows to specify periodic fluxes
-         do chunk = 3,Nchunks                                                                  ! loop through chunks (skipping the keyword)
-            v = IO_lc(IO_stringValue(line,chunkPos,chunk))                                          ! chunk matches keyvalues x,y, or z?
-            mesh_periodicSurface(1) = mesh_periodicSurface(1) .or. v == 'x'
-            mesh_periodicSurface(2) = mesh_periodicSurface(2) .or. v == 'y'
-            mesh_periodicSurface(3) = mesh_periodicSurface(3) .or. v == 'z'
-         enddo
-     endselect
-   endif
- enddo
-
-end subroutine mesh_get_damaskOptions
-
-
-
 subroutine calcCells(thisMesh,connectivity_elem)
 
   class(tMesh) :: thisMesh
@@ -893,11 +851,17 @@ subroutine calcCells(thisMesh,connectivity_elem)
   integer(pInt),dimension(:,:), allocatable :: con_elem,temp,con,parentsAndWeights,candidates_global
   integer(pInt),dimension(:), allocatable :: l, nodes, candidates_local
   integer(pInt),dimension(:,:,:), allocatable :: con_cell,connectivity_cell
-  integer(pInt),dimension(:,:), allocatable :: sorted,test
+  integer(pInt),dimension(:,:), allocatable :: sorted,test,connectivity_cell_reshape
   real(pReal), dimension(:,:), allocatable :: coordinates,nodes5
   integer(pInt) :: e, n, c, p, s,u,i,m,j,nParentNodes,nCellNode,ierr
 
-  !nodes5 = thisMesh%node_0
+#if defined(DAMASK_HDF5)
+ call results_openJobFile
+ call HDF5_closeGroup(results_addGroup('geometry'))
+ call results_writeDataset('geometry',connectivity_elem,'connectivity_element',&
+                           'connectivity of the elements','-')
+#endif
+
 !---------------------------------------------------------------------------------------------------
 ! initialize global connectivity to negative local connectivity
   allocate(connectivity_cell(thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs,thisMesh%Nelems))
@@ -1020,8 +984,13 @@ subroutine calcCells(thisMesh,connectivity_elem)
      if (p/=0) nodes5 = reshape([nodes5,coordinates],[3,nCellNode])
   enddo
   thisMesh%node_0 = nodes5
-
+ connectivity_cell_reshape = reshape(connectivity_cell,[thisMesh%elem%NcellNodesPerCell,thisMesh%elem%nIPs*thisMesh%Nelems])
-  
+ 
+#if defined(DAMASK_HDF5) 
+ call results_writeDataset('geometry',connectivity_cell_reshape,'connectivity_cell',&
+                           'connectivity of the cells','-')
+ call results_closeJobFile
+#endif  
 end subroutine calcCells
 
 !--------------------------------------------------------------------------------------------------

src/plastic_nonlocal.f90

@@ -8,7 +8,6 @@ module plastic_nonlocal
   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, &