From 94a24e45eeaca801c7c1e7dd0404d165865a646c Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Sat, 2 Feb 2019 11:46:38 +0100
Subject: [PATCH] bugfixes:

theMesh%Nelems need to be set (using an intermediate function until a
routine does that)
spectral.geom file can have "N+n to N" (backwards counting)
---
 src/mesh_FEM.f90    |  1 +
 src/mesh_abaqus.f90 |  1 +
 src/mesh_base.f90   | 12 ++++++++++++
 src/mesh_grid.f90   | 20 +++++++++++++-------
 src/mesh_marc.f90   |  2 +-
 5 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90
index b9d171fc3..e2b08db4c 100644
--- a/src/mesh_FEM.f90
+++ b/src/mesh_FEM.f90
@@ -276,6 +276,7 @@ subroutine mesh_init()
 !!!!!!!!!!!!!!!!!!!!!!!!
  allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
  call theMesh%init(dimplex,integrationOrder,mesh_node0)
+ call theMesh%setNelems(mesh_NcpElems)
  
 end subroutine mesh_init
 
diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90
index 5c761ad7e..05e1d7c7d 100644
--- a/src/mesh_abaqus.f90
+++ b/src/mesh_abaqus.f90
@@ -421,6 +421,7 @@ subroutine tMesh_abaqus_init(self,elemType,nodes)
  integer(pInt), intent(in) :: elemType
  
  call self%tMesh%init('mesh',elemType,nodes)
+  call theMesh%setNelems(mesh_NcpElems)
  
 end subroutine tMesh_abaqus_init
 
diff --git a/src/mesh_base.f90 b/src/mesh_base.f90
index e0ca78c03..c0f012256 100644
--- a/src/mesh_base.f90
+++ b/src/mesh_base.f90
@@ -45,6 +45,7 @@ module mesh_base
      connectivity
    contains
    procedure, pass(self) :: tMesh_base_init
+   procedure :: setNelems =>  tMesh_base_setNelems                                                  ! not needed once we compute the cells from the connectivity
    generic, public :: init => tMesh_base_init
  end type tMesh
 
@@ -68,4 +69,15 @@ subroutine tMesh_base_init(self,meshType,elemType,nodes)
 
 end subroutine tMesh_base_init
 
+
+subroutine tMesh_base_setNelems(self,Nelems)
+ 
+  implicit none
+  class(tMesh) :: self
+  integer(pInt), intent(in) :: Nelems
+
+  self%Nelems = Nelems
+
+end subroutine tMesh_base_setNelems
+
 end module mesh_base
diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90
index 942611b1f..88484a693 100644
--- a/src/mesh_grid.f90
+++ b/src/mesh_grid.f90
@@ -442,7 +442,7 @@ subroutine mesh_init(ip,el)
  mesh_microstructureAt  = mesh_element(4,:)
  mesh_CPnodeID          = mesh_element(5:4+mesh_NipsPerElem,:)
 !!!!!!!!!!!!!!!!!!!!!!!!
- 
+ call theMesh%setNelems(mesh_NcpElems)
 end subroutine mesh_init
 
 !--------------------------------------------------------------------------------------------------
@@ -686,6 +686,8 @@ end function mesh_cellCenterCoordinates
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Parses geometry file
+!> @details important variables have an implicit "save" attribute. Therefore, this function is 
+! supposed to be called only once!
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_spectral_read_grid()
  use IO, only: &
@@ -706,12 +708,16 @@ subroutine mesh_spectral_read_grid()
   real(pReal), dimension(3) :: s = -1_pInt
   integer(pInt) :: h =- 1_pInt
   integer(pInt) ::  &
-    headerLength = -1_pInt, &
-    fileLength, &
+    headerLength = -1_pInt, &                                                                       !< length of header (in lines)
+    fileLength, &                                                                                   !< lenght of the geom file (in characters)
     fileUnit, &
     startPos, endPos, &
     myStat, &
-    l, i, j, e, c
+    l, &                                                                                            !< line counter
+    c, &                                                                                            !< counter for # microstructures in line
+    o, &                                                                                            !< order of "to" packing
+    e, &                                                                                            !< "element", i.e. spectral collocation point 
+    i, j
   logical :: &
     gotGrid = .false., &
     gotSize = .false., &
@@ -807,8 +813,9 @@ subroutine mesh_spectral_read_grid()
         c = IO_intValue(line,chunkPos,1)
         microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,3),i = 1_pInt,IO_intValue(line,chunkPos,1))]
       else if (IO_lc(IO_stringValue(line,chunkPos,2))  == 'to') then
-        c = IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1) + 1_pInt
-        microGlobal(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3))]
+        c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1_pInt
+        o = merge(+1_pInt, -1_pInt, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1))
+        microGlobal(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)]
       else
         c = chunkPos(1)
         do i = 0_pInt, c - 1_pInt
@@ -822,7 +829,6 @@ subroutine mesh_spectral_read_grid()
       enddo
 
     endif
-
     e = e+c
   end do
 
diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90
index 5607791fb..dd4098879 100644
--- a/src/mesh_marc.f90
+++ b/src/mesh_marc.f90
@@ -560,7 +560,7 @@ subroutine mesh_init(ip,el)
  mesh_CPnodeID          = mesh_element(5:4+mesh_NipsPerElem,:)
 !!!!!!!!!!!!!!!!!!!!!!!!
  call theMesh%init(mesh_element(2,1),mesh_node0)
- 
+ call theMesh%setNelems(mesh_NcpElems)
 end subroutine mesh_init