Merge branch 'nonlocal-again-with-Marc' into 'development'

Nonlocal again with marc

See merge request damask/DAMASK!99

src/IO.f90

@@ -44,7 +44,6 @@ module IO
     IO_extractValue, &
     IO_countDataLines
 #elif defined(Marc4DAMASK)
-    IO_skipChunks, &
     IO_fixedNoEFloatValue, &
     IO_fixedIntValue, &
     IO_countNumericalDataLines
@@ -403,7 +402,7 @@ pure function IO_stringPos(string)
    left  = right + verify(string(right+1:),SEP)
    right = left + scan(string(left:),SEP) - 2
    if ( string(left:left) == '#' ) exit
-   IO_stringPos = [IO_stringPos,left, right]
+   IO_stringPos = [IO_stringPos,left,right]
    IO_stringPos(1) = IO_stringPos(1)+1
    endOfString: if (right < left) then
      IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string)
@@ -1023,27 +1022,6 @@ integer function IO_countNumericalDataLines(fileUnit)
   backspace(fileUnit)
 
 end function IO_countNumericalDataLines
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief reads file to skip (at least) N chunks (may be over multiple lines)
-!--------------------------------------------------------------------------------------------------
-subroutine IO_skipChunks(fileUnit,N)
-
-  integer, intent(in)  :: fileUnit, &                                                               !< file handle
-                          N                                                                         !< minimum number of chunks to skip
- 
-  integer              :: remainingChunks
-  character(len=65536) :: line
- 
-  line = ''
-  remainingChunks = N
- 
-  do while (trim(line) /= IO_EOF .and. remainingChunks > 0)
-    line = IO_read(fileUnit)
-    remainingChunks = remainingChunks - (size(IO_stringPos(line))-1)/2
-  enddo
-end subroutine IO_skipChunks
 #endif
 
 

src/commercialFEM_fileList.f90

@@ -14,11 +14,11 @@
 #include "Lambert.f90"
 #include "rotations.f90"
 #include "FEsolving.f90"
-#include "geometry_plastic_nonlocal.f90"
 #include "element.f90"
 #include "mesh_base.f90"
 #include "HDF5_utilities.f90"
 #include "results.f90"
+#include "geometry_plastic_nonlocal.f90"
 #include "discretization.f90"
 #ifdef Abaqus
 #include "mesh_abaqus.f90"

src/element.f90

@@ -146,11 +146,11 @@ module element
        8  & ! 3D 8node
     ]                                                                                               !< number of cell nodes in a specific cell type
 
-  ! *** FE_ipNeighbor ***
+  ! *** IPneighbor ***
-  ! is a list of the neighborhood of each IP.
+  ! list of the neighborhood of each IP.
   ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction.
-  ! Positive integers denote an intra-FE IP identifier.
+  ! Positive integers denote an intra-element IP identifier.
-  ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located.
+  ! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located.
  
   integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = &
     reshape([&
@@ -256,10 +256,6 @@ module element
       -3,26,-4,24,-6,18  &
      ],[nIPneighbor(cellType(10)),nIP(10)])
 
-! MD: probably not needed END
-! --------------------------------------------------------------------------------------------------
-
-
  
   integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
     reshape([&
@@ -757,7 +753,7 @@ subroutine tElement_init(self,elemType)
        self%cell        = CELL10
    end select
 
-   self%NcellNodesPerCell = NCELLNODEPERCELL(self%cellType)
+   self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType)
   
    select case(self%cellType)
      case(1)

src/geometry_plastic_nonlocal.f90

@@ -7,6 +7,7 @@
 !--------------------------------------------------------------------------------------------------
 module geometry_plastic_nonlocal
   use prec
+  use results
 
   implicit none
   private
@@ -32,6 +33,7 @@ module geometry_plastic_nonlocal
     geometry_plastic_nonlocal_setIPvolume, &
     geometry_plastic_nonlocal_setIParea, &
     geometry_plastic_nonlocal_setIPareaNormal, &
+    geometry_plastic_nonlocal_results, &
     geometry_plastic_nonlocal_disable
     
 contains
@@ -112,4 +114,45 @@ subroutine geometry_plastic_nonlocal_disable
   
 end subroutine geometry_plastic_nonlocal_disable
 
+
+!---------------------------------------------------------------------------------------------------
+!> @brief Writes geometry data to results file
+!---------------------------------------------------------------------------------------------------
+subroutine geometry_plastic_nonlocal_results
+  
+  integer,     dimension(:),   allocatable :: shp
+
+#if defined(DAMASK_HDF5)
+  call results_openJobFile
+
+  writeVolume: block
+    real(pReal), dimension(:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IPvolume0)
+    temp = reshape(geometry_plastic_nonlocal_IPvolume0,[shp(1)*shp(2)])
+    call results_writeDataset('geometry',temp,'v_0',&
+                              'initial cell volume','m³')
+  end block writeVolume
+
+  writeAreas: block
+    real(pReal), dimension(:,:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IParea0)
+    temp = reshape(geometry_plastic_nonlocal_IParea0,[shp(1),shp(2)*shp(3)])
+    call results_writeDataset('geometry',temp,'a_0',&
+                              'initial cell face area','m²')
+  end block writeAreas
+
+  writeNormals: block
+    real(pReal), dimension(:,:,:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IPareaNormal0)
+    temp = reshape(geometry_plastic_nonlocal_IPareaNormal0,[shp(1),shp(2),shp(3)*shp(4)])
+    call results_writeDataset('geometry',temp,'n_0',&
+                              'initial cell face normals','-',transposed=.false.)
+  end block writeNormals
+
+
+  call results_closeJobFile
+#endif
+  
+end subroutine geometry_plastic_nonlocal_results
+
 end module geometry_plastic_nonlocal

src/math.f90

@@ -1307,10 +1307,10 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4)
   real(pReal), dimension (3,3) :: m
 
   m(1:3,1) = v1-v2
-  m(1:3,2) = v2-v3
+  m(1:3,2) = v1-v3
-  m(1:3,3) = v3-v4
+  m(1:3,3) = v1-v4
 
-  math_volTetrahedron = math_det33(m)/6.0_pReal
+  math_volTetrahedron = abs(math_det33(m))/6.0_pReal
 
 end function math_volTetrahedron
 
@@ -1404,6 +1404,7 @@ subroutine unitTest
   integer,     dimension(5) :: range_out_ = [1,2,3,4,5]
 
   real(pReal)                 :: det
+  real(pReal), dimension(3)   :: v3_1,v3_2,v3_3,v3_4
   real(pReal), dimension(6)   :: v6
   real(pReal), dimension(9)   :: v9
   real(pReal), dimension(3,3) :: t33,t33_2
@@ -1452,6 +1453,15 @@ subroutine unitTest
   if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     call IO_error(401,ext_msg='math_symmetric33')
  
+  call random_number(v3_1)
+  call random_number(v3_2)
+  call random_number(v3_3)
+  call random_number(v3_4)
+  
+  if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
+          math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
+  call IO_error(401,ext_msg='math_volTetrahedron')
+    
   call random_number(t33)
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     call IO_error(401,ext_msg='math_det33/math_detSym33')

src/mesh_FEM.f90

@@ -28,8 +28,7 @@ module mesh
  integer, public, protected :: &
    mesh_Nboundaries, &
    mesh_NcpElems, &                                                                                 !< total number of CP elements in mesh
-   mesh_NcpElemsGlobal, &
+   mesh_NcpElemsGlobal
-   mesh_Nnodes                                                                                      !< total number of nodes in mesh
 
 !!!! BEGIN DEPRECATED !!!!!
  integer, public, protected :: &
@@ -69,8 +68,7 @@ module mesh
  public :: &
    mesh_init, &
    mesh_FEM_build_ipVolumes, &
-   mesh_FEM_build_ipCoordinates, &
+   mesh_FEM_build_ipCoordinates
-   mesh_cellCenterCoordinates
 
 contains
 
@@ -107,7 +105,8 @@ subroutine mesh_init
  integer, parameter :: FILEUNIT = 222
  integer :: j
  integer, allocatable, dimension(:) :: chunkPos
- integer :: dimPlex
+ integer :: dimPlex, &
+   mesh_Nnodes                                                                                      !< total number of nodes in mesh
  integer, parameter :: &
    mesh_ElemType=1                                                                             !< Element type of the mesh (only support homogeneous meshes)
  character(len=512) :: &
@@ -228,26 +227,8 @@ subroutine mesh_init
 end subroutine mesh_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Calculates cell center coordinates.
-!--------------------------------------------------------------------------------------------------
-pure function mesh_cellCenterCoordinates(ip,el)
- 
- integer, 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
-
-end function mesh_cellCenterCoordinates
-
-
 !--------------------------------------------------------------------------------------------------
 !> @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_FEM_build_ipVolumes(dimPlex)
  

src/mesh_marc.f90

@@ -44,7 +44,6 @@ module mesh
     mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
 !--------------------------------------------------------------------------------------------------
 
-
  public :: &
    mesh_init, &
    mesh_FEasCP
@@ -60,30 +59,17 @@ subroutine mesh_init(ip,el)
 
   integer, intent(in) :: el, ip
    
-  integer, parameter  :: FILEUNIT = 222
-  character(len=pStringLen), dimension(:), allocatable :: inputFile                                 !< file content, separated per lines
- integer :: &
-   mesh_NelemSets
   real(pReal), dimension(:,:), allocatable :: &
    node0_elem, &                                                                                       !< node x,y,z coordinates (initially!)
    node0_cell
+  type(tElement) :: elem
 
-  integer :: j, fileFormatVersion, elemType, &
+  integer :: nElems
-    mesh_maxNelemInSet, &
-    mesh_nElems, &
-    hypoelasticTableStyle, &
-    initialcondTableStyle
-  integer, dimension(:), allocatable :: &
-   marc_matNumber                                                                                   !< array of material numbers for hypoelastic material (Marc only)
   integer,     dimension(:),   allocatable :: &
     microstructureAt, &
     homogenizationAt
- character(len=64), dimension(:), allocatable :: &
-   mesh_nameElemSet
- integer, dimension(:,:), allocatable :: &
-   mesh_mapElemSet                                                                                  !< list of elements in elementSet
   integer:: &
-   mesh_Nnodes                                                                                      !< total number of nodes in mesh
+    Nnodes                                                                                      !< total number of nodes in mesh
    
   real(pReal), dimension(:,:), allocatable :: &
     ip_reshaped
@@ -91,89 +77,55 @@ subroutine mesh_init(ip,el)
     connectivity_cell                                                                                !< cell connectivity for each element,ip/cell
   integer, dimension(:,:), allocatable :: &
     connectivity_elem
-   
+  real(pReal), dimension(:,:,:,:),allocatable :: &
-
+    unscaledNormals
- type(tMesh) :: theMesh
- 
 
   write(6,'(/,a)')   ' <<<+-  mesh init  -+>>>'
  
   mesh_unitlength = numerics_unitlength                                                             ! set physical extent of a length unit in mesh
-  inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension))
 
-  ! parsing Marc input file
+  call inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
-  fileFormatVersion = inputRead_fileFormat(inputFile)
+  nElems = size(connectivity_elem,2)
-  call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle,inputFile)
-  if (fileFormatVersion > 12) &
-    marc_matNumber = inputRead_matNumber(hypoelasticTableStyle,inputFile)
-  call inputRead_NnodesAndElements(mesh_nNodes, mesh_nElems, inputFile)
 
-  call IO_open_inputFile(FILEUNIT,modelName)
+  if (debug_e < 1 .or. debug_e > nElems)    call IO_error(602,ext_msg='element')
-  call inputRead_NelemSets(mesh_NelemSets,mesh_maxNelemInSet,FILEUNIT) 
+  if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
-  allocate(mesh_nameElemSet(mesh_NelemSets)); mesh_nameElemSet = 'n/a'
-  allocate(mesh_mapElemSet(1+mesh_maxNelemInSet,mesh_NelemSets),source=0)
-  call inputRead_mapElemSets(mesh_nameElemSet,mesh_mapElemSet,FILEUNIT)
-  allocate (mesh_mapFEtoCPelem(2,mesh_nElems), source = 0)
-  call inputRead_mapElems(hypoelasticTableStyle,mesh_nameElemSet,mesh_mapElemSet,&
-                              mesh_nElems,fileFormatVersion,marc_matNumber,FILEUNIT)  
-  allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes),source=0)
-  call inputRead_mapNodes(mesh_Nnodes,inputFile)                                                                  !ToDo: don't work on global variables
  
-  node0_elem = inputRead_elemNodes(mesh_Nnodes,inputFile)
+  FEsolving_execElem = [ 1,nElems ]                                                                 ! parallel loop bounds set to comprise all DAMASK elements
+  allocate(FEsolving_execIP(2,nElems), source=1)                                                    ! parallel loop bounds set to comprise from first IP...
+  FEsolving_execIP(2,:) = elem%nIPs
+ 
+  allocate(calcMode(elem%nIPs,nElems),source=.false.)                                               ! pretend to have collected what first call is asking (F = I)
+  calcMode(ip,mesh_FEasCP('elem',el)) = .true.                                                      ! first ip,el needs to be already pingponged to "calc"
  
 
-  elemType = inputRead_elemType(mesh_nElems,FILEUNIT)
+  allocate(mesh_ipCoordinates(3,elem%nIPs,nElems),source=0.0_pReal)                                 ! deprecated
 
-  call theMesh%init('mesh',elemType,node0_elem)
+  allocate(cellNodeDefinition(elem%nNodes-1))
-  call theMesh%setNelems(mesh_nElems)
+  allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
-  
-  allocate(microstructureAt(theMesh%nElems), source=0)
-  allocate(homogenizationAt(theMesh%nElems), source=0)
- 
-  connectivity_elem = inputRead_connectivityElem(mesh_nElems,theMesh%elem%nNodes,FILEUNIT)
-  call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
-                               mesh_nElems,theMesh%elem%nNodes,mesh_nameElemSet,mesh_mapElemSet,&
-                               initialcondTableStyle,FILEUNIT)
-  close (FILEUNIT)
-
-
-  allocate(mesh_ipCoordinates(3,theMesh%elem%nIPs,theMesh%nElems),source=0.0_pReal) 
-
-  allocate(cellNodeDefinition(theMesh%elem%nNodes-1))
-  allocate(connectivity_cell(theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs,theMesh%nElems))
   call buildCells(connectivity_cell,cellNodeDefinition,&
-                  theMesh%elem,connectivity_elem)
+                  elem,connectivity_elem)
   allocate(node0_cell(3,maxval(connectivity_cell)))
   call buildCellNodes(node0_cell,&
                       cellNodeDefinition,node0_elem)
-  allocate(ip_reshaped(3,theMesh%elem%nIPs*theMesh%nElems),source=0.0_pReal)
+  allocate(ip_reshaped(3,elem%nIPs*nElems),source=0.0_pReal)
-  call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,&
+  call buildIPcoordinates(ip_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,&
-                                    theMesh%elem%nIPs*theMesh%nElems]),node0_cell)
+                          elem%nIPs*nElems]),node0_cell)
- 
- 
-  if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) &
-    call IO_error(600)                                                                          ! ping-pong must be disabled when having non-DAMASK elements
-  if (debug_e < 1 .or. debug_e > theMesh%nElems) &
-    call IO_error(602,ext_msg='element')                                                        ! selected element does not exist
-  if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) &
-    call IO_error(602,ext_msg='IP')                                                             ! selected element does not have requested IP
- 
-  FEsolving_execElem = [ 1,theMesh%nElems ]                                                      ! parallel loop bounds set to comprise all DAMASK elements
-  allocate(FEsolving_execIP(2,theMesh%nElems), source=1)                                    ! parallel loop bounds set to comprise from first IP...
-  FEsolving_execIP(2,:) = theMesh%elem%nIPs
- 
-  allocate(calcMode(theMesh%elem%nIPs,theMesh%nElems))
-  calcMode = .false.                                                                                 ! pretend to have collected what first call is asking (F = I)
-  calcMode(ip,mesh_FEasCP('elem',el)) = .true.                                                       ! first ip,el needs to be already pingponged to "calc"
 
   call discretization_init(microstructureAt,homogenizationAt,&
                            ip_reshaped,&
                            node0_elem)
                            
   call writeGeometry(0,connectivity_elem,&
-                           reshape(connectivity_cell,[theMesh%elem%NcellNodesPerCell,theMesh%elem%nIPs*theMesh%nElems]),&
+                     reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*nElems]),&
                      node0_cell,ip_reshaped)
 
+  call geometry_plastic_nonlocal_setIPvolume(IPvolume(elem,node0_cell,connectivity_cell))
+  unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell)
+  call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1))
+  call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3))
+  call geometry_plastic_nonlocal_results
+  
+
 end subroutine mesh_init
 
 !--------------------------------------------------------------------------------------------------
@@ -216,17 +168,85 @@ subroutine writeGeometry(elemType, &
   call results_writeDataset('geometry',coordinates_temp,'x_p', &
                             'coordinates of the material points','m')
                       
-  call results_closeJobFile()
+  call results_closeJobFile
 #endif  
 
 end subroutine writeGeometry
 
 
+subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
+
+  type(tElement), intent(out) :: elem
+  real(pReal), dimension(:,:), allocatable, intent(out) :: &
+   node0_elem                                                                                       !< node x,y,z coordinates (initially!)
+  integer, dimension(:,:), allocatable, intent(out) :: &
+    connectivity_elem
+
+  integer,     dimension(:),   allocatable, intent(out) :: &
+    microstructureAt, &
+    homogenizationAt
+   
+  integer :: &
+    fileFormatVersion, &
+    hypoelasticTableStyle, &
+    initialcondTableStyle, &
+    nNodes, &
+    nElems
+  integer, parameter :: &
+    FILEUNIT = 222
+  integer, dimension(:), allocatable :: &
+    matNumber                                                                                       !< material numbers for hypoelastic material
+  character(len=pStringLen), dimension(:), allocatable :: inputFile                                 !< file content, separated per lines
+  
+  character(len=64), dimension(:), allocatable :: &
+    nameElemSet
+  integer, dimension(:,:), allocatable :: &
+    mapElemSet                                                                                  !< list of elements in elementSet
+
+  inputFile = IO_read_ASCII(trim(modelName)//trim(InputFileExtension))
+  call inputRead_fileFormat(fileFormatVersion, &
+                            inputFile)
+  call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &
+                             inputFile)
+  if (fileFormatVersion > 12) &
+    call inputRead_matNumber(matNumber, &
+                             hypoelasticTableStyle,inputFile)
+  call inputRead_NnodesAndElements(nNodes,nElems,&
+                                   inputFile)
+  
+  call IO_open_inputFile(FILEUNIT,modelName)                                                        ! ToDo: It would be better to use fileContent
+  
+  call inputRead_mapElemSets(nameElemSet,mapElemSet,&
+                             FILEUNIT)
+
+  allocate (mesh_mapFEtoCPelem(2,nElems), source=0)
+  call inputRead_mapElems(hypoelasticTableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,FILEUNIT)
+
+  allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0)
+  call inputRead_mapNodes(inputFile)
+
+  call inputRead_elemType(elem, &
+                          nElems,inputFile)
+  call inputRead_elemNodes(node0_elem, &
+                           Nnodes,inputFile)
+
+  connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile)
+
+  call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
+                                                 nElems,elem%nNodes,nameElemSet,mapElemSet,&
+                                                 initialcondTableStyle,FILEUNIT)
+  close(FILEUNIT)
+   
+end subroutine inputRead
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Figures out version of Marc input file format
 !--------------------------------------------------------------------------------------------------
-integer function inputRead_fileFormat(fileContent)
+subroutine inputRead_fileFormat(fileFormat,fileContent)
  
+  integer,                                 intent(out) :: fileFormat
   character(len=pStringLen), dimension(:), intent(in)  :: fileContent                                !< file content, separated per lines
   
   integer, allocatable, dimension(:) :: chunkPos
@@ -235,12 +255,12 @@ integer function inputRead_fileFormat(fileContent)
   do l = 1, size(fileContent)
     chunkPos = IO_stringPos(fileContent(l))
     if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then
-      inputRead_fileFormat = IO_intValue(fileContent(l),chunkPos,2)
+      fileFormat = IO_intValue(fileContent(l),chunkPos,2)
       exit
     endif
   enddo
 
-end function inputRead_fileFormat
+end subroutine inputRead_fileFormat
 
 
 !--------------------------------------------------------------------------------------------------
@@ -272,12 +292,13 @@ end subroutine inputRead_tableStyles
 !--------------------------------------------------------------------------------------------------
 !> @brief Figures out material number of hypoelastic material
 !--------------------------------------------------------------------------------------------------
-function inputRead_matNumber(tableStyle,fileContent)
+subroutine inputRead_matNumber(matNumber, &
+                               tableStyle,fileContent)
  
+  integer, allocatable,      dimension(:), intent(out) :: matNumber
   integer,                                 intent(in)  :: tableStyle
   character(len=pStringLen), dimension(:), intent(in)  :: fileContent                               !< file content, separated per lines
 
-  integer, dimension(:), allocatable :: inputRead_matNumber
 
   integer, allocatable, dimension(:) :: chunkPos
   integer :: i, j, data_blocks, l
@@ -291,17 +312,17 @@ function inputRead_matNumber(tableStyle,fileContent)
       else
         data_blocks = 1
       endif
-      allocate(inputRead_matNumber(data_blocks), source = 0)
+      allocate(matNumber(data_blocks), source = 0)
       do i = 0, data_blocks - 1
         j = i*(2+tableStyle) + 1
         chunkPos = IO_stringPos(fileContent(l+1+j))
-        inputRead_matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1)
+        matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1)
       enddo
       exit
     endif
   enddo
 
-end function inputRead_matNumber
+end subroutine inputRead_matNumber
 
 
 !--------------------------------------------------------------------------------------------------
@@ -367,14 +388,18 @@ subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
 !--------------------------------------------------------------------------------------------------
 subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit)
  
-  character(len=64), dimension(:),   intent(out) :: nameElemSet
+  character(len=64), dimension(:),   allocatable, intent(out) :: nameElemSet
-  integer,           dimension(:,:), intent(out) :: mapElemSet
+  integer,           dimension(:,:), allocatable, intent(out) :: mapElemSet
   integer,                                        intent(in)  :: fileUnit
 
   integer, allocatable, dimension(:) :: chunkPos
   character(len=300) :: line
-  integer :: elemSet
+  integer :: elemSet, NelemSets, maxNelemInSet
 
+
+  call inputRead_NelemSets(NelemSets,maxNelemInSet,fileUnit) 
+  allocate(nameElemSet(NelemSets)); nameElemSet = 'n/a'
+  allocate(mapElemSet(1+maxNelemInSet,NelemSets),source=0)
   elemSet = 0
 
   rewind(fileUnit)
@@ -396,9 +421,9 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,fileUnit)
 !--------------------------------------------------------------------------------------------------
 !> @brief Maps elements from FE ID to internal (consecutive) representation.
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileFormatVersion,matNumber,fileUnit)
+subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,fileFormatVersion,matNumber,fileUnit)
  
- integer, intent(in) :: fileUnit,tableStyle,nElems,fileFormatVersion
+  integer, intent(in) :: fileUnit,tableStyle,fileFormatVersion
   integer, dimension(:), intent(in) :: matNumber
   character(len=64), intent(in), dimension(:) :: nameElemSet
   integer, dimension(:,:), intent(in) :: &
@@ -408,25 +433,27 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileForma
   character(len=300) :: line, &
                         tmp
  
- integer, dimension (1+nElems) :: contInts
+  integer, dimension(:), allocatable :: contInts
   integer :: i,cpElem 
  
+  allocate(contInts(size(mesh_mapFEtoCPelem,2)+1))
+ 
   cpElem = 0
   contInts = 0
   rewind(fileUnit)
   do
     read (fileUnit,'(A300)',END=620) line
     chunkPos = IO_stringPos(line)
-   if (fileFormatVersion < 13) then                                                                       ! Marc 2016 or earlier
+    Marc2016andEarlier: if (fileFormatVersion < 13) then
       if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'hypoelastic' ) then
-       do i=1,3+TableStyle                                                                          ! skip three (or four if new table style!) lines
+        skipLines: do i=1,3+TableStyle
           read (fileUnit,'(A300)') line
-       enddo
+        enddo skipLines
-       contInts = IO_continuousIntValues(fileUnit,nElems,nameElemSet,&
+        contInts = IO_continuousIntValues(fileUnit,size(mesh_mapFEtoCPelem,2),nameElemSet,&
                                           mapElemSet,size(nameElemSet))
         exit
       endif  
-   else                                                                                             ! Marc2017 and later
+    else Marc2016andEarlier
       if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity') then
         read (fileUnit,'(A300)',END=620) line
         chunkPos = IO_stringPos(line)
@@ -444,7 +471,7 @@ subroutine inputRead_mapElems(tableStyle,nameElemSet,mapElemSet,nElems,fileForma
           enddo
         endif  
       endif
-   endif    
+    endif Marc2016andEarlier
   enddo    
 620 do i = 1,contInts(1)
     cpElem = cpElem+1
@@ -460,9 +487,8 @@ end subroutine inputRead_mapElems
 !--------------------------------------------------------------------------------------------------
 !> @brief Maps node from FE ID to internal (consecutive) representation.
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead_mapNodes(nNodes,fileContent)
+subroutine inputRead_mapNodes(fileContent)
 
-  integer, intent(in)                                 :: nNodes
   character(len=pStringLen), dimension(:), intent(in) :: fileContent                                !< file content, separated per lines
 
   integer, allocatable, dimension(:) :: chunkPos
@@ -471,7 +497,7 @@ subroutine inputRead_mapNodes(nNodes,fileContent)
   do l = 1, size(fileContent)
     chunkPos = IO_stringPos(fileContent(l))
     if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
-      do i = 1,nNodes
+      do i = 1,size(mesh_mapFEtoCPnode,2)
         mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (fileContent(l+1+i),[0,10],1)
         mesh_mapFEtoCPnode(2,i) = i
       enddo
@@ -487,16 +513,19 @@ end subroutine inputRead_mapNodes
 !--------------------------------------------------------------------------------------------------
 !> @brief store x,y,z coordinates of all nodes in mesh.
 !--------------------------------------------------------------------------------------------------
-function inputRead_elemNodes(nNode,fileContent) result(nodes)
+subroutine inputRead_elemNodes(nodes, &
+                               nNode,fileContent)
 
+  real(pReal), allocatable,  dimension(:,:), intent(out) :: nodes 
   integer,                                   intent(in)  :: nNode
   character(len=pStringLen), dimension(:),   intent(in)  :: fileContent                             !< file content, separated per lines
   
-  real(pReal), dimension(3,nNode)    :: nodes 
   integer, dimension(5), parameter   :: node_ends = [0,10,30,50,70]
   integer, allocatable, dimension(:) :: chunkPos
   integer :: i,j,m,l
 
+  allocate(nodes(3,nNode))
+
   do l = 1, size(fileContent)
     chunkPos = IO_stringPos(fileContent(l))
     if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates' ) then
@@ -510,42 +539,41 @@ function inputRead_elemNodes(nNode,fileContent) result(nodes)
     endif
   enddo
 
-end function inputRead_elemNodes
+end subroutine inputRead_elemNodes
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Gets element type (and checks if the whole mesh comprises of only one type)
 !--------------------------------------------------------------------------------------------------
-integer function inputRead_elemType(nElem,fileUnit)
+subroutine inputRead_elemType(elem, &
+                              nElem,fileContent)
 
-  integer, intent(in) :: &
+  type(tElement),                          intent(out) :: elem
-    nElem, &
+  integer,                                 intent(in)  :: nElem
-    fileUnit
+  character(len=pStringLen), dimension(:), intent(in)  :: fileContent                               !< file content, separated per lines
 
-  type(tElement) :: tempEl
   integer, allocatable, dimension(:) :: chunkPos
-  character(len=300) :: line
+  integer :: i,j,t,l,remainingChunks
-  integer :: i,t
 
   t = -1
-
+  do l = 1, size(fileContent)
-  rewind(fileUnit)
+    chunkPos = IO_stringPos(fileContent(l))
-  do
+    if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then
-    read (fileUnit,'(A300)',END=620) line
+      j = 0
-    chunkPos = IO_stringPos(line)
-    if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then
-      read (fileUnit,'(A300)') line                                                                 ! Garbage line
       do i=1,nElem                                                                                  ! read all elements
-        read (fileUnit,'(A300)') line
+        chunkPos = IO_stringPos(fileContent(l+1+i+j))
-        chunkPos = IO_stringPos(line)
         if (t == -1) then
-          t = mapElemtype(IO_stringValue(line,chunkPos,2))
+          t = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))
-          call tempEl%init(t)
+          call elem%init(t)
-          inputRead_elemType = t
         else
-          if (t /= mapElemtype(IO_stringValue(line,chunkPos,2))) call IO_error(191,el=t,ip=i)
+          if (t /= mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))) call IO_error(191,el=t,ip=i)
         endif
-        call IO_skipChunks(fileUnit,tempEl%nNodes-(chunkPos(1)-2))
+        remainingChunks = elem%nNodes - (chunkPos(1) - 2)
+        do while(remainingChunks > 0)
+          j = j + 1
+          chunkPos = IO_stringPos(fileContent(l+1+i+j))
+          remainingChunks = remainingChunks - chunkPos(1)
+        enddo
       enddo
       exit
     endif
@@ -594,54 +622,49 @@ integer function inputRead_elemType(nElem,fileUnit)
         call IO_error(error_ID=190,ext_msg=IO_lc(what))
    end select
 
-end function mapElemtype
+  end function mapElemtype
 
 
-620 end function inputRead_elemType
+end subroutine inputRead_elemType
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Stores node IDs
 !--------------------------------------------------------------------------------------------------
-function inputRead_connectivityElem(nElem,nNodes,fileUnit)
+function inputRead_connectivityElem(nElem,nNodes,fileContent)
  
   integer, intent(in) :: &
     nElem, &
-    nNodes, &                                                                                       !< number of nodes per element
+    nNodes                                                                                          !< number of nodes per element
-    fileUnit
+  character(len=pStringLen), dimension(:), intent(in)  :: fileContent                               !< file content, separated per lines
     
   integer, dimension(nNodes,nElem) :: &
     inputRead_connectivityElem
  
   integer, allocatable, dimension(:) :: chunkPos
-  character(len=300) line
  
   integer, dimension(1+nElem) :: contInts
-  integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead
+  integer :: i,k,j,t,e,l,nNodesAlreadyRead
  
-  rewind(fileUnit)
+  do l = 1, size(fileContent)
-  do
+    chunkPos = IO_stringPos(fileContent(l))
-    read (fileUnit,'(A300)',END=620) line
+    if( IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity' ) then
-    chunkPos = IO_stringPos(line)
+      j = 0
-    if( IO_lc(IO_stringValue(line,chunkPos,1)) == 'connectivity' ) then
-      read (fileUnit,'(A300)',END=620) line                                                         ! garbage line
       do i = 1,nElem
-        read (fileUnit,'(A300)',END=620) line
+        chunkPos = IO_stringPos(fileContent(l+1+i+j))
-        chunkPos = IO_stringPos(line)
+        e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1))
-        e = mesh_FEasCP('elem',IO_intValue(line,chunkPos,1))
         if (e /= 0) then                                                                            ! disregard non CP elems
-          nNodesAlreadyRead = 0
+          do k = 1,chunkPos(1)-2
-          do j = 1,chunkPos(1)-2
+            inputRead_connectivityElem(k,e) = &
-            inputRead_connectivityElem(j,e) = &
+              mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
-            mesh_FEasCP('node',IO_IntValue(line,chunkPos,j+2))
           enddo
           nNodesAlreadyRead = chunkPos(1) - 2
           do while(nNodesAlreadyRead < nNodes)                                                      ! read on if not all nodes in one line
-            read (fileUnit,'(A300)',END=620) line
+            j = j + 1
-            chunkPos = IO_stringPos(line)
+            chunkPos = IO_stringPos(fileContent(l+1+i+j))
-            do j = 1,chunkPos(1)
+            do k = 1,chunkPos(1)
-              inputRead_connectivityElem(nNodesAlreadyRead+j,e) = &
+              inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
-              mesh_FEasCP('node',IO_IntValue(line,chunkPos,j))
+                mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
             enddo
             nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
           enddo
@@ -651,7 +674,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileUnit)
     endif
   enddo
 
-620 end function inputRead_connectivityElem
+end function inputRead_connectivityElem
 
 
 !--------------------------------------------------------------------------------------------------
@@ -660,7 +683,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileUnit)
 subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
                                     nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileUnit)
 
-  integer, dimension(:), intent(out) :: &
+  integer, dimension(:), allocatable, intent(out) :: &
     microstructureAt, &
     homogenizationAt
   integer, intent(in) :: &
@@ -679,6 +702,10 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
   integer, dimension(1+nElem) :: contInts
   integer :: i,j,t,sv,myVal,e,nNodesAlreadyRead
 
+
+  allocate(microstructureAt(nElem),source=0)
+  allocate(homogenizationAt(nElem),source=0)
+
   rewind(fileUnit)
   read (fileUnit,'(A300)',END=630) line
   do
@@ -718,7 +745,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
 630 end subroutine inputRead_microstructureAndHomogenization
 
 
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculates cell node coordinates from element node coordinates
 !--------------------------------------------------------------------------------------------------
@@ -930,6 +956,115 @@ subroutine buildIPcoordinates(IPcoordinates, &
 end subroutine buildIPcoordinates
 
 
+!---------------------------------------------------------------------------------------------------
+!> @brief Calculates IP volume.
+!> @details The IP volume is calculated differently depending on the cell type.
+!> 2D cells assume an element depth of 1.0
+!---------------------------------------------------------------------------------------------------
+function IPvolume(elem,node,connectivity)
+ 
+  type(tElement),                intent(in) :: elem
+  real(pReal), dimension(:,:),   intent(in) :: node
+  integer,     dimension(:,:,:), intent(in) :: connectivity
+  
+  real(pReal), dimension(elem%nIPs,size(connectivity,3)) :: IPvolume
+  real(pReal), dimension(3) :: x0,x1,x2,x3,x4,x5,x6,x7
+
+  integer :: e,i
+
+  do e = 1,size(connectivity,3)
+    do i = 1,elem%nIPs
+
+      select case (elem%cellType)
+        case (1)                                                                                    ! 2D 3node
+          IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), &
+                                            node(1:3,connectivity(2,i,e)), &
+                                            node(1:3,connectivity(3,i,e)))
+
+        case (2)                                                                                    ! 2D 4node
+          IPvolume(i,e) = math_areaTriangle(node(1:3,connectivity(1,i,e)), &                        ! assume planar shape, division in two triangles suffices
+                                            node(1:3,connectivity(2,i,e)), &
+                                            node(1:3,connectivity(3,i,e))) &
+                        + math_areaTriangle(node(1:3,connectivity(3,i,e)), &
+                                            node(1:3,connectivity(4,i,e)), &
+                                            node(1:3,connectivity(1,i,e)))
+        case (3)                                                                                    ! 3D 4node
+          IPvolume(i,e) = math_volTetrahedron(node(1:3,connectivity(1,i,e)), &
+                                              node(1:3,connectivity(2,i,e)), &
+                                              node(1:3,connectivity(3,i,e)), &
+                                              node(1:3,connectivity(4,i,e)))
+        case (4)                                                                                    ! 3D 8node
+          ! J. Grandy, Efficient Calculation of Volume of Hexahedral Cells
+          ! Lawrence Livermore National Laboratory
+          ! https://www.osti.gov/servlets/purl/632793
+          x0 = node(1:3,connectivity(1,i,e))
+          x1 = node(1:3,connectivity(2,i,e))
+          x2 = node(1:3,connectivity(4,i,e))
+          x3 = node(1:3,connectivity(3,i,e))
+          x4 = node(1:3,connectivity(5,i,e))
+          x5 = node(1:3,connectivity(6,i,e))
+          x6 = node(1:3,connectivity(8,i,e))
+          x7 = node(1:3,connectivity(7,i,e))
+          IPvolume(i,e) = dot_product((x7-x1)+(x6-x0),math_cross((x7-x2),        (x3-x0))) &
+                        + dot_product((x6-x0),        math_cross((x7-x2)+(x5-x0),(x7-x4))) &
+                        + dot_product((x7-x1),        math_cross((x5-x0),        (x7-x4)+(x3-x0)))
+          IPvolume(i,e) = IPvolume(i,e)/12.0_pReal
+      end select
+    enddo
+  enddo
+
+end function IPvolume
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculation of IP interface areas
+!--------------------------------------------------------------------------------------------------
+function IPareaNormal(elem,nElem,connectivity,node)
+
+  type(tElement),                intent(in) :: elem
+  integer,                       intent(in) :: nElem
+  integer,     dimension(:,:,:), intent(in) :: connectivity
+  real(pReal), dimension(:,:),   intent(in) :: node
+  
+  real(pReal), dimension(3,elem%nIPneighbors,elem%nIPs,nElem) :: ipAreaNormal
+
+  real(pReal), dimension (3,size(elem%cellFace,1)) :: nodePos
+  integer :: e,i,f,n,m
+ 
+  m = size(elem%cellFace,1)
+
+  do e = 1,nElem
+    do i = 1,elem%nIPs
+      do f = 1,elem%nIPneighbors
+        nodePos = node(1:3,connectivity(elem%cellface(1:m,f),i,e))
+        
+        select case (elem%cellType)
+          case (1,2)                                                                                ! 2D 3 or 4 node
+            IPareaNormal(1,f,i,e) =   nodePos(2,2) - nodePos(2,1)                                   ! x_normal =  y_connectingVector
+            IPareaNormal(2,f,i,e) = -(nodePos(1,2) - nodePos(1,1))                                  ! y_normal = -x_connectingVector
+            IPareaNormal(3,f,i,e) = 0.0_pReal
+          case (3)                                                                                  ! 3D 4node
+            IPareaNormal(1:3,f,i,e) = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
+                                                 nodePos(1:3,3) - nodePos(1:3,1))
+          case (4)                                                                                  ! 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
+            IPareaNormal(1:3,f,i,e) = 0.0_pReal
+            do n = 1, m 
+              IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) &
+                                      + math_cross(nodePos(1:3,mod(n+0,m)+1) - nodePos(1:3,n), &
+                                                   nodePos(1:3,mod(n+1,m)+1) - nodePos(1:3,n)) * 0.5_pReal
+            enddo
+        end select
+      enddo
+    enddo
+  enddo
+
+end function IPareaNormal
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Gives the FE to CP ID mapping by binary search through lookup array
 !! valid questions (what) are 'elem', 'node'

src/results.f90

@@ -296,21 +296,34 @@ end subroutine results_writeVectorDataset_real
 !--------------------------------------------------------------------------------------------------
 !> @brief stores a tensor dataset in a group
 !--------------------------------------------------------------------------------------------------
-subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit)
+subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit,transposed)
 
   character(len=*), intent(in)                   :: label,group,description
   character(len=*), intent(in), optional         :: SIunit
+  logical,          intent(in), optional         :: transposed
   real(pReal),      intent(in), dimension(:,:,:) :: dataset
   
   integer :: i 
+  logical :: T
   integer(HID_T) :: groupHandle
   real(pReal), dimension(:,:,:), allocatable :: dataset_transposed
 
   
+  if(present(transposed)) then
+    T = transposed
+  else
+    T = .true.
+  endif
+
+  if(T) then
+    if(size(dataset,1) /= size(dataset,2)) call IO_error(0,ext_msg='transpose non-symmetric tensor')
     allocate(dataset_transposed,mold=dataset)
-  do i=1,size(dataset,3)
+    do i=1,size(dataset_transposed,3)
-    dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i))
+      dataset_transposed(:,:,i) = transpose(dataset(:,:,i))
     enddo
+  else
+    allocate(dataset_transposed,source=dataset)
+  endif
 
   groupHandle = results_openGroup(group)