Merge remote-tracking branch 'origin/development' into 32_NewStyleNonlocal-3

.gitlab-ci.yml

@@ -7,9 +7,9 @@ stages:
   - compilePETScGNU
   - prepareSpectral
   - spectral
-  - compileMarc2018_1
+  - compileMarc
   - marc
-  - compileAbaqus2017
+  - compileAbaqus
   - example
   - performance
   - createPackage
@@ -51,39 +51,37 @@ variables:
   # Names of module files to load
   # ===============================================================================================
   # ++++++++++++ Compiler ++++++++++++++++++++++++++++++++++++++++++++++
-  IntelCompiler16_0:               "Compiler/Intel/16.0 Libraries/IMKL/2016"
+  IntelCompiler16_4:               "Compiler/Intel/16.4 Libraries/IMKL/2016"
-  IntelCompiler16_4:               "Compiler/Intel/16.4 Libraries/IMKL/2016-4"
+  IntelCompiler17_8:               "Compiler/Intel/17.8 Libraries/IMKL/2017"
-  IntelCompiler17_0:               "Compiler/Intel/17.0 Libraries/IMKL/2017"
+  IntelCompiler18_4:               "Compiler/Intel/18.4 Libraries/IMKL/2018"
-  IntelCompiler18_1:               "Compiler/Intel/18.1 Libraries/IMKL/2018"
+  GNUCompiler8_2:                  "Compiler/GNU/8.2"
-  GNUCompiler7_3:                  "Compiler/GNU/7.3"
   # ------------ Defaults ----------------------------------------------
-  IntelCompiler:                   "$IntelCompiler18_1"
+  IntelCompiler:                   "$IntelCompiler18_4"
-  GNUCompiler:                     "$GNUCompiler7_3"
+  GNUCompiler:                     "$GNUCompiler8_2"
   # ++++++++++++ MPI +++++++++++++++++++++++++++++++++++++++++++++++++++
-  MPICH3_2Intel18_1:               "MPI/Intel/18.1/MPICH/3.2.1"
+  IMPI2018Intel18_4:               "MPI/Intel/18.4/IntelMPI/2018"
-  MPICH3_2GNU7_3:                  "MPI/GNU/7.3/MPICH/3.2.1"
+  MPICH3_3GNU8_2:                  "MPI/GNU/8.2/MPICH/3.3"
   # ------------ Defaults ----------------------------------------------
-  MPICH_Intel:                     "$MPICH3_2Intel18_1"
+  MPICH_Intel:                     "$IMPI2018Intel18_4"
-  MPICH_GNU:                       "$MPICH3_2GNU7_3"
+  MPICH_GNU:                       "$MPICH3_3GNU8_2"
   # ++++++++++++ PETSc +++++++++++++++++++++++++++++++++++++++++++++++++
-  PETSc3_10_0MPICH3_2Intel18_1:    "Libraries/PETSc/3.10.0/Intel-18.1-MPICH-3.2.1"
+  PETSc3_10_3IMPI2018Intel18_4:    "Libraries/PETSc/3.10.3/Intel-18.4-IntelMPI-2018"
-  PETSc3_10_0MPICH3_2GNU7_3:       "Libraries/PETSc/3.10.0/GNU-7.3-MPICH-3.2.1"
+  PETSc3_10_3MPICH3_3GNU8_2:       "Libraries/PETSc/3.10.3/GNU-8.2-MPICH-3.3"
   # ------------ Defaults ----------------------------------------------
-  PETSc_MPICH_Intel:               "$PETSc3_10_0MPICH3_2Intel18_1"
+  PETSc_MPICH_Intel:               "$PETSc3_10_3IMPI2018Intel18_4"
-  PETSc_MPICH_GNU:                 "$PETSc3_10_0MPICH3_2GNU7_3"
+  PETSc_MPICH_GNU:                 "$PETSc3_10_3MPICH3_3GNU8_2"
   # ++++++++++++ FEM +++++++++++++++++++++++++++++++++++++++++++++++++++
-  Abaqus2017:                      "FEM/Abaqus/2017"
+  Abaqus2019:                      "FEM/Abaqus/2019"
   MSC2018_1:                       "FEM/MSC/2018.1"
-  MSC2017:                         "FEM/MSC/2017"
   # ------------ Defaults ----------------------------------------------
-  Abaqus:                          "$Abaqus2017"
+  Abaqus:                          "$Abaqus2019"
   MSC:                             "$MSC2018_1"
-  IntelMarc:                       "$IntelCompiler17_0"
+  IntelMarc:                       "$IntelCompiler17_8"
   IntelAbaqus:                     "$IntelCompiler16_4"
   # ++++++++++++ Documentation +++++++++++++++++++++++++++++++++++++++++
-  Doxygen1_8_13:                   "Documentation/Doxygen/1.8.13"
+  Doxygen1_8_15:                   "Documentation/Doxygen/1.8.15"
   # ------------ Defaults ----------------------------------------------
-  Doxygen:                         "$Doxygen1_8_13"
+  Doxygen:                         "$Doxygen1_8_15"
 
 
 ###################################################################################################
@@ -158,6 +156,13 @@ Post_AverageDown:
     - master
     - release
 
+Post_ASCIItable:
+  stage: postprocessing
+  script: ASCIItable/test.py
+  except:
+    - master
+    - release
+
 Post_General:
   stage: postprocessing
   script: PostProcessing/test.py
@@ -383,9 +388,9 @@ TextureComponents:
 
 ###################################################################################################
 Marc_compileIfort2018_1:
-  stage: compileMarc2018_1
+  stage: compileMarc
   script:
-    - module load $IntelCompiler17_0 $MSC2018_1
+    - module load $IntelMarc $MSC
     - Marc_compileIfort/test.py -m 2018.1
   except:
     - master
@@ -430,11 +435,11 @@ J2_plasticBehavior:
     - release
 
 ###################################################################################################
-Abaqus_compile2017:
+Abaqus_compile:
-  stage: compileAbaqus2017
+  stage: compileAbaqus
   script:
-    - module load $IntelCompiler16_4 $Abaqus2017
+    - module load $IntelAbaqus $Abaqus
-    - Abaqus_compileIfort/test.py -a 2017
+    - Abaqus_compileIfort/test.py
   except:
     - master
     - release

CONFIG

@@ -8,6 +8,6 @@ set DAMASK_NUM_THREADS    = 4
 set MSC_ROOT              = /opt/msc
 set MARC_VERSION          = 2018.1
 
-set ABAQUS_VERSION        = 2017
+set ABAQUS_VERSION        = 2019
 
 set DAMASK_HDF5           = OFF

PRIVATE

@@ -1 +1 @@
-Subproject commit beb9682fff7d4d6c65aba12ffd04c7441dc6ba6b
+Subproject commit 18ba1ba6a5e9ba446dc9311acf2acf2781614db1

README

@@ -10,3 +10,4 @@ Germany
 
 Email: DAMASK@mpie.de
 https://damask.mpie.de
+https://magit1.mpie.de

VERSION

@@ -1 +1 @@
-v2.0.2-1689-g1a471bcd
+v2.0.2-1789-g524bfb8c

processing/post/addCumulative.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys

processing/post/addDerivative.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys

processing/post/addLinked.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys
@@ -21,7 +21,7 @@ Add data of selected column(s) from (first) row of linked ASCIItable that shares
 parser.add_option('--link',
                   dest = 'link', nargs = 2,
                   type = 'string', metavar = 'string string',
-                  help = 'column labels containing linked values')
+                  help = 'column labels of table and linked table containing linking values')
 parser.add_option('-l','--label',
                   dest = 'label',
                   action = 'extend', metavar = '<string LIST>',
@@ -105,7 +105,8 @@ for name in filenames:
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
     try:
-      table.data_append(data[np.argwhere(np.all((map(float,table.data[myLink:myLink+myLinkDim]) - index)==0,axis=1))[0]])  # add data of first matching line
+      table.data_append(data[np.argwhere(np.all((list(map(float,table.data[myLink:myLink+myLinkDim])) - index)==0,
+                                                 axis=1))[0]])                                      # add data of first matching line
     except IndexError:
       table.data_append(np.nan*np.ones_like(data[0]))                                               # or add NaNs
     outputAlive = table.data_write()                                                                # output processed line

processing/post/addTable.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys

processing/post/blowUp.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys

processing/post/vtk_addGridData.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,vtk

processing/post/vtk_addPointcloudData.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,vtk

processing/post/vtk_addRectilinearGridData.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,vtk

processing/post/vtk_pointcloud.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys,vtk

processing/post/vtk_rectilinearGrid.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys,vtk

processing/pre/geom_grainGrowth.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys,math
@@ -49,7 +49,7 @@ parser.set_defaults(d = 1,
 
 (options, filenames) = parser.parse_args()
 
-options.immutable = map(int,options.immutable)
+options.immutable = list(map(int,options.immutable))
 
 getInterfaceEnergy = lambda A,B: np.float32((A*B != 0)*(A != B)*1.0)                               # 1.0 if A & B are distinct & nonzero, 0.0 otherwise
 struc = ndimage.generate_binary_structure(3,1)                                                     # 3D von Neumann neighborhood
@@ -70,9 +70,9 @@ for name in filenames:
   table.head_read()
   info,extra_header = table.head_getGeom()
   
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(map(str,info['grid']))),
+  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(list(map(str,info['grid'])))),
-                     'size     x y z:  {}'.format(' x '.join(map(str,info['size']))),
+                     'size     x y z:  {}'.format(' x '.join(list(map(str,info['size'])))),
-                     'origin   x y z:  {}'.format(' : '.join(map(str,info['origin']))),
+                     'origin   x y z:  {}'.format(' : '.join(list(map(str,info['origin'])))),
                      'homogenization:  {}'.format(info['homogenization']),
                      'microstructures: {}'.format(info['microstructures']),
                     ])
@@ -102,9 +102,9 @@ for name in filenames:
     gauss = np.exp(-(X*X + Y*Y + Z*Z)/(2.0*options.d*options.d),dtype=np.float32) \
             /np.power(2.0*np.pi*options.d*options.d,(3.0 - np.count_nonzero(info['grid'] == 1))/2.,dtype=np.float32)
           
-    gauss[:,:,:grid[2]/2:-1] = gauss[:,:,1:(grid[2]+1)/2]     # trying to cope with uneven (odd) grid size
+    gauss[:,:,:grid[2]//2:-1] = gauss[:,:,1:(grid[2]+1)//2]     # trying to cope with uneven (odd) grid size
-    gauss[:,:grid[1]/2:-1,:] = gauss[:,1:(grid[1]+1)/2,:]
+    gauss[:,:grid[1]//2:-1,:] = gauss[:,1:(grid[1]+1)//2,:]
-    gauss[:grid[0]/2:-1,:,:] = gauss[1:(grid[0]+1)/2,:,:]
+    gauss[:grid[0]//2:-1,:,:] = gauss[1:(grid[0]+1)//2,:,:]
     gauss = np.fft.rfftn(gauss).astype(np.complex64)
 
   for smoothIter in range(options.N):
@@ -119,9 +119,9 @@ for name in filenames:
                                                           microstructure,i,axis=0), j,axis=1), k,axis=2)))
 
     # periodically extend interfacial energy array by half a grid size in positive and negative directions
-    periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]/2:-grid[0]/2,
+    periodic_interfaceEnergy = np.tile(interfaceEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
-                                                                grid[1]/2:-grid[1]/2,
+                                                                grid[1]//2:-grid[1]//2,
-                                                                grid[2]/2:-grid[2]/2]
+                                                                grid[2]//2:-grid[2]//2]
 
     # transform bulk volume (i.e. where interfacial energy remained zero), store index of closest boundary voxel
     index = ndimage.morphology.distance_transform_edt(periodic_interfaceEnergy == 0.,                                            
@@ -148,15 +148,15 @@ for name in filenames:
                          ndimage.morphology.binary_dilation(interfaceEnergy > 0.,
                                                             structure = struc,
                                                             iterations = int(round(options.d*2.))-1),# fat boundary
-                         periodic_bulkEnergy[grid[0]/2:-grid[0]/2,                                 # retain filled energy on fat boundary...
+                         periodic_bulkEnergy[grid[0]//2:-grid[0]//2,                                 # retain filled energy on fat boundary...
-                                             grid[1]/2:-grid[1]/2,
+                                             grid[1]//2:-grid[1]//2,
-                                             grid[2]/2:-grid[2]/2],                                # ...and zero everywhere else
+                                             grid[2]//2:-grid[2]//2],                                # ...and zero everywhere else
                          0.)).astype(np.complex64) *
                          gauss).astype(np.float32)
 
-      periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]/2:-grid[0]/2,
+      periodic_diffusedEnergy = np.tile(diffusedEnergy,(3,3,3))[grid[0]//2:-grid[0]//2,
-                                                                grid[1]/2:-grid[1]/2,
+                                                                grid[1]//2:-grid[1]//2,
-                                                                grid[2]/2:-grid[2]/2]              # periodically extend the smoothed bulk energy
+                                                                grid[2]//2:-grid[2]//2]              # periodically extend the smoothed bulk energy
 
 
     # transform voxels close to interface region
@@ -164,15 +164,15 @@ for name in filenames:
                                                       return_distances = False,
                                                       return_indices = True)                       # want index of closest bulk grain
 
-    periodic_microstructure = np.tile(microstructure,(3,3,3))[grid[0]/2:-grid[0]/2,
+    periodic_microstructure = np.tile(microstructure,(3,3,3))[grid[0]//2:-grid[0]//2,
-                                                              grid[1]/2:-grid[1]/2,
+                                                              grid[1]//2:-grid[1]//2,
-                                                              grid[2]/2:-grid[2]/2]                # periodically extend the microstructure
+                                                              grid[2]//2:-grid[2]//2]                # periodically extend the microstructure
 
     microstructure = periodic_microstructure[index[0],
                                              index[1],
-                                             index[2]].reshape(2*grid)[grid[0]/2:-grid[0]/2,
+                                             index[2]].reshape(2*grid)[grid[0]//2:-grid[0]//2,
-                                                                       grid[1]/2:-grid[1]/2,
+                                                                       grid[1]//2:-grid[1]//2,
-                                                                       grid[2]/2:-grid[2]/2]       # extent grains into interface region
+                                                                       grid[2]//2:-grid[2]//2]       # extent grains into interface region
 
     # replace immutable microstructures with closest mutable ones
     index = ndimage.morphology.distance_transform_edt(np.in1d(microstructure,options.immutable).reshape(grid),
@@ -236,3 +236,4 @@ for name in filenames:
 # --- output finalization --------------------------------------------------------------------------
 
   table.close()     
+  

processing/pre/geom_toTable.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import os,sys
@@ -48,11 +48,11 @@ for name in filenames:
   table.head_read()
   info,extra_header = table.head_getGeom()
   
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
+  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(list(map(str,info['grid'])))),
-                     'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
+                     'size     x y z:  {}'.format(' x '.join(list(map(str,info['size'])))),
-                     'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
+                     'origin   x y z:  {}'.format(' : '.join(list(map(str,info['origin'])))),
-                     'homogenization:  %i'%info['homogenization'],
+                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: %i'%info['microstructures'],
+                     'microstructures: {}'.format(info['microstructures']),
                     ])
 
   errors = []

processing/pre/seeds_fromDistribution.py

@@ -1,10 +1,11 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 # -*- coding: UTF-8 no BOM -*-
 
 import threading,time,os,sys,random
 import numpy as np
 from optparse import OptionParser
-from cStringIO import StringIO
+from io import StringIO
+import binascii
 import damask
 
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
@@ -96,7 +97,7 @@ class myThread (threading.Thread):
       perturbedGeomVFile = StringIO()
       perturbedSeedsVFile.reset()
       perturbedGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
-                     ' -g '+' '.join(map(str, options.grid)),streamIn=perturbedSeedsVFile)[0])
+                     ' -g '+' '.join(list(map(str, options.grid))),streamIn=perturbedSeedsVFile)[0])
       perturbedGeomVFile.reset()
 
 #--- evaluate current seeds file ----------------------------------------------------------------------
@@ -214,7 +215,7 @@ options = parser.parse_args()[0]
 damask.util.report(scriptName,options.seedFile)
 
 if options.randomSeed is None:
-  options.randomSeed = int(os.urandom(4).encode('hex'), 16)
+  options.randomSeed = int(binascii.hexlify(os.urandom(4)),16)
 damask.util.croak(options.randomSeed)
 delta = (options.scale/options.grid[0],options.scale/options.grid[1],options.scale/options.grid[2])
 baseFile=os.path.splitext(os.path.basename(options.seedFile))[0]
@@ -240,17 +241,17 @@ if os.path.isfile(os.path.splitext(options.seedFile)[0]+'.seeds'):
     for line in initialSeedFile: bestSeedsVFile.write(line)
 else:
   bestSeedsVFile.write(damask.util.execute('seeds_fromRandom'+\
-                                ' -g '+' '.join(map(str, options.grid))+\
+                                ' -g '+' '.join(list(map(str, options.grid)))+\
                                 ' -r {:d}'.format(options.randomSeed)+\
                                 ' -N '+str(nMicrostructures))[0])
 bestSeedsUpdate = time.time()
 
 # ----------- tessellate initial seed file to get and evaluate geom file
-bestSeedsVFile.reset()
+bestSeedsVFile.seek(0)
 initialGeomVFile = StringIO()
 initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
-                               ' -g '+' '.join(map(str, options.grid)),bestSeedsVFile)[0])
+                               ' -g '+' '.join(list(map(str, options.grid))),bestSeedsVFile)[0])
-initialGeomVFile.reset()
+initialGeomVFile.seek(0)
 initialGeomTable = damask.ASCIItable(initialGeomVFile,None,labeled=False,readonly=True)
 initialGeomTable.head_read()
 info,devNull =  initialGeomTable.head_getGeom()

python/damask/solver/abaqus.py

@@ -2,7 +2,7 @@
 
 from .solver import Solver
 import damask
-import subprocess,re
+import subprocess
 
 class Abaqus(Solver):
 
@@ -15,14 +15,13 @@ class Abaqus(Solver):
 
   def return_run_command(self,model):
     env=damask.Environment()
-    shortVersion = re.sub('[\.,-]', '',self.version)
     try:
-      cmd='abq'+shortVersion
+      cmd='abq'+self.version
-      subprocess.check_output(['abq'+shortVersion,'information=release'])
+      subprocess.check_output([cmd,'information=release'])
     except OSError:                                                                                 # link to abqXXX not existing
       cmd='abaqus'
       process = subprocess.Popen(['abaqus','information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE)
-      detectedVersion = process.stdout.readlines()[1].split()[1]
+      detectedVersion = process.stdout.readlines()[1].split()[1].decode('utf-8')
       if self.version != detectedVersion:
-        raise Exception('found Abaqus version %s, but requested %s'%(detectedVersion,self.version))
+        raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version))
-    return '%s -job %s -user %s/src/DAMASK_abaqus interactive'%(cmd,model,env.rootDir())
+    return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,env.rootDir())

src/CMakeLists.txt

@@ -7,6 +7,7 @@ endif()
 
 # The dependency detection in CMake is not functioning for Fortran,
 # hence we declare the dependencies from top to bottom in the following
+
 add_library(C_ROUTINES OBJECT "C_routines.c")
 set(OBJECTFILES $<TARGET_OBJECTS:C_ROUTINES>)
 
@@ -17,6 +18,10 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
 add_library(PREC OBJECT "prec.f90")
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
 
+add_library(ELEMENT OBJECT "element.f90")
+add_dependencies(ELEMENT PREC)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:ELEMENT>)
+
 add_library(QUIT OBJECT "quit.f90")
 add_dependencies(QUIT PREC)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUIT>)
@@ -34,7 +39,7 @@ add_dependencies(NUMERICS IO)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:NUMERICS>)
 
 add_library(DEBUG OBJECT "debug.f90")
-add_dependencies(DEBUG NUMERICS)
+add_dependencies(DEBUG IO)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:DEBUG>)
 
 add_library(DAMASK_CONFIG OBJECT "config.f90")
@@ -42,7 +47,7 @@ add_dependencies(DAMASK_CONFIG DEBUG)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CONFIG>)
 
 add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90")
-add_dependencies(HDF5_UTILITIES DAMASK_CONFIG)
+add_dependencies(HDF5_UTILITIES DAMASK_CONFIG NUMERICS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:HDF5_UTILITIES>)
 
 add_library(RESULTS OBJECT "results.f90")
@@ -50,24 +55,28 @@ add_dependencies(RESULTS HDF5_UTILITIES)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:RESULTS>)
 
 add_library(FEsolving OBJECT "FEsolving.f90")
-add_dependencies(FEsolving RESULTS)
+add_dependencies(FEsolving DEBUG)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEsolving>)
 
-add_library(DAMASK_MATH OBJECT "math.f90")
+add_library(MATH OBJECT "math.f90")
-add_dependencies(DAMASK_MATH FEsolving)
+add_dependencies(MATH NUMERICS)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_MATH>)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>)
+
+add_library(MESH_BASE OBJECT "mesh_base.f90")
+add_dependencies(MESH_BASE ELEMENT)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>)
 
 # SPECTRAL solver and FEM solver use different mesh files
 if (PROJECT_NAME STREQUAL "DAMASK_spectral")
-      add_library(MESH OBJECT "mesh.f90")
+      add_library(MESH OBJECT "mesh_grid.f90")
-      add_dependencies(MESH DAMASK_MATH)
+      add_dependencies(MESH MESH_BASE MATH FEsolving)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
 elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
       add_library(FEZoo OBJECT "FEM_zoo.f90")
-      add_dependencies(FEZoo DAMASK_MATH)
+      add_dependencies(FEZoo IO)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
-      add_library(MESH  OBJECT "meshFEM.f90")
+      add_library(MESH  OBJECT "mesh_FEM.f90")
-      add_dependencies(MESH FEZoo)
+      add_dependencies(MESH FEZoo MESH_BASE MATH FEsolving)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
 endif()
 
@@ -75,9 +84,9 @@ add_library(MATERIAL OBJECT "material.f90")
 add_dependencies(MATERIAL MESH DAMASK_CONFIG)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>)
 
-add_library(DAMASK_HELPERS OBJECT "lattice.f90")
+add_library(LATTICE OBJECT "lattice.f90")
-add_dependencies(DAMASK_HELPERS MATERIAL)
+add_dependencies(LATTICE MATERIAL)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_HELPERS>)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:LATTICE>)
 
 # For each modular section
 add_library (PLASTIC OBJECT
@@ -88,14 +97,14 @@ add_library (PLASTIC OBJECT
   "plastic_kinematichardening.f90"
   "plastic_nonlocal.f90"
   "plastic_none.f90")
-add_dependencies(PLASTIC DAMASK_HELPERS)
+add_dependencies(PLASTIC LATTICE RESULTS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:PLASTIC>)
 
 add_library (KINEMATICS OBJECT
   "kinematics_cleavage_opening.f90"
   "kinematics_slipplane_opening.f90"
   "kinematics_thermal_expansion.f90")
-add_dependencies(KINEMATICS DAMASK_HELPERS)
+add_dependencies(KINEMATICS LATTICE RESULTS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:KINEMATICS>)
 
 add_library (SOURCE OBJECT
@@ -105,7 +114,7 @@ add_library (SOURCE OBJECT
   "source_damage_isoDuctile.f90"
   "source_damage_anisoBrittle.f90"
   "source_damage_anisoDuctile.f90")
-add_dependencies(SOURCE DAMASK_HELPERS)
+add_dependencies(SOURCE LATTICE RESULTS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:SOURCE>)
 
 add_library(CONSTITUTIVE OBJECT "constitutive.f90")

src/CPFEM.f90

@@ -140,8 +140,7 @@ subroutine CPFEM_init
    restartRead, &
    modelName
  use mesh, only: &
-   mesh_NcpElems, &
+   theMesh
-   mesh_maxNips
  use material, only: &
    material_phase, &
    homogState, &
@@ -168,10 +167,9 @@ subroutine CPFEM_init
    flush(6)
  endif mainProcess
 
- ! initialize stress and jacobian to zero
+ allocate(CPFEM_cs(               6,theMesh%elem%nIPs,theMesh%Nelems), source= 0.0_pReal)
- allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ;                CPFEM_cs              = 0.0_pReal
+ allocate(CPFEM_dcsdE(          6,6,theMesh%elem%nIPs,theMesh%Nelems), source= 0.0_pReal)
- allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ;           CPFEM_dcsdE           = 0.0_pReal
+ allocate(CPFEM_dcsdE_knownGood(6,6,theMesh%elem%nIPs,theMesh%Nelems), source= 0.0_pReal)
- allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE_knownGood = 0.0_pReal
 
  ! *** restore the last converged values of each essential variable from the binary file
  if (restartRead) then
@@ -289,8 +287,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
    math_6toSym33
  use mesh, only: &
    mesh_FEasCP, &
-   mesh_NcpElems, &
+   theMesh, &
-   mesh_maxNips, &
    mesh_element
  use material, only: &
    microstructure_elemhomo, &
@@ -401,7 +398,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
    enddo; enddo
    if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
      write(6,'(a)') '<< CPFEM >> aging states'
-     if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then
+     if (debug_e <= theMesh%Nelems .and. debug_i <= theMesh%elem%nIPs) then
        write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') &
              '<< CPFEM >> aged state of elFE ip grain',debug_e, debug_i, 1, &
               plasticState(phaseAt(1,debug_i,debug_e))%state(:,phasememberAt(1,debug_i,debug_e))

src/CPFEM2.f90

@@ -95,8 +95,6 @@ subroutine CPFEM_init
  use prec, only: &
    pInt, pReal, pLongInt
  use IO, only: &
-   IO_read_realFile,&
-   IO_read_intFile, &
    IO_timeStamp, &
    IO_error
  use numerics, only: &

src/DAMASK_abaqus.f

@@ -30,6 +30,11 @@ contains
 !> @brief reports and sets working directory
 !--------------------------------------------------------------------------------------------------
 subroutine DAMASK_interface_init
+#if __INTEL_COMPILER >= 1800
+ use, intrinsic :: iso_fortran_env, only: &
+   compiler_version, &
+   compiler_options
+#endif
  use ifport, only: &
    CHDIR
  
@@ -40,16 +45,25 @@ subroutine DAMASK_interface_init
  character(len=256) :: wd
 
  call date_and_time(values = dateAndTime)
- write(6,'(/,a)') ' <<<+-  DAMASK_abaqus_std  -+>>>'
+ write(6,'(/,a)') ' <<<+-  DAMASK_abaqus  -+>>>'
- write(6,'(/,a)') ' Roters et al., Computational Materials Science, 2018'
+ write(6,'(/,a)') ' Roters et al., Computational Materials Science 158, 2018, 420-478'
- write(6,'(/,a)')              ' Version: '//DAMASKVERSION
+ write(6,'(a,/)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
- write(6,'(a,2(i2.2,a),i4.4)') ' Date:    ',dateAndTime(3),'/',&
+
-                                            dateAndTime(2),'/',&
+ write(6,'(a,/)')              ' Version: '//DAMASKVERSION
-                                            dateAndTime(1) 
+
- write(6,'(a,2(i2.2,a),i2.2)') ' Time:    ',dateAndTime(5),':',&
+! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
-                                            dateAndTime(6),':',&
+#if __INTEL_COMPILER >= 1800
-                                            dateAndTime(7)  
+ write(6,*) 'Compiled with: ', compiler_version()
- write(6,'(/,a)') ' <<<+-  DAMASK_interface init  -+>>>'
+ write(6,*) 'Compiler options: ', compiler_options()
+#else
+ write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
+                                                    ', build date :', __INTEL_COMPILER_BUILD_DATE
+#endif
+
+ write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__
+
+ write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
+ write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
 
  call getoutdir(wd, lenOutDir)
  ierr = CHDIR(wd)

src/DAMASK_interface.f90

@@ -143,16 +143,27 @@ subroutine DAMASK_interface_init()
 
  call date_and_time(values = dateAndTime)
  write(6,'(/,a)') ' <<<+-  DAMASK_interface init  -+>>>'
- write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018'
+ write(6,'(/,a)') ' Roters et al., Computational Materials Science 158, 2018, 420-478'
- write(6,'(/,a)')              ' Version: '//DAMASKVERSION
+ write(6,'(a,/)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
- write(6,'(a,2(i2.2,a),i4.4)') ' Date:    ',dateAndTime(3),'/',&
+
-                                            dateAndTime(2),'/',&
+ write(6,'(a,/)')              ' Version: '//DAMASKVERSION
-                                            dateAndTime(1) 
+
- write(6,'(a,2(i2.2,a),i2.2)') ' Time:    ',dateAndTime(5),':',&
+! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
-                                            dateAndTime(6),':',&
+#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
-                                            dateAndTime(7)  
+ write(6,*) 'Compiled with: ', compiler_version()
- write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
+ write(6,*) 'Compiler options: ', compiler_options()
-#include "compilation_info.f90"
+#elif defined(__INTEL_COMPILER)
+ write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
+                                                    ', build date :', __INTEL_COMPILER_BUILD_DATE
+#elif defined(__PGI)
+ write(6,'(a,i4.4,a,i8.8)') ' Compiled with PGI fortran version :', __PGIC__,&
+                                                               '.', __PGIC_MINOR__
+#endif
+
+ write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__
+
+ write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
+ write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
  
  call get_command(commandLine)
  chunkPos = IIO_stringPos(commandLine)
@@ -219,9 +230,11 @@ subroutine DAMASK_interface_init()
 
  call get_environment_variable('USER',userName)
  ! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
+ write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
  write(6,'(a,a)')      ' Host name: ', trim(getHostName())
  write(6,'(a,a)')      ' User name: ', trim(userName)
- write(6,'(a,a)')      ' Command line call:      ', trim(commandLine)
+
+ write(6,'(/a,a)')     ' Command line call:      ', trim(commandLine)
  if (len(trim(workingDirArg)) > 0) &
    write(6,'(a,a)')    ' Working dir argument:   ', trim(workingDirArg)
  write(6,'(a,a)')      ' Geometry argument:      ', trim(geometryArg)

src/DAMASK_marc.f90

@@ -43,6 +43,11 @@ contains
 !> @brief reports and sets working directory
 !--------------------------------------------------------------------------------------------------
 subroutine DAMASK_interface_init
+#if __INTEL_COMPILER >= 1800
+ use, intrinsic :: iso_fortran_env, only: &
+   compiler_version, &
+   compiler_options
+#endif
  use ifport, only: &
    CHDIR
 
@@ -53,17 +58,26 @@ subroutine DAMASK_interface_init
  character(len=1024) :: wd
 
  call date_and_time(values = dateAndTime)
- write(6,'(/,a)') ' <<<+-  DAMASK_Marc  -+>>>'
+ write(6,'(/,a)') ' <<<+-  DAMASK_abaqus  -+>>>'
- write(6,'(/,a)') ' Roters et al., Computational Materials Science, 2018'
+ write(6,'(/,a)') ' Roters et al., Computational Materials Science 158, 2018, 420-478'
- write(6,'(/,a)')              ' Version: '//DAMASKVERSION
+ write(6,'(a,/)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
- write(6,'(a,2(i2.2,a),i4.4)') ' Date:    ',dateAndTime(3),'/',&
+
-                                            dateAndTime(2),'/',&
+ write(6,'(a,/)')              ' Version: '//DAMASKVERSION
-                                            dateAndTime(1)
+
- write(6,'(a,2(i2.2,a),i2.2)') ' Time:    ',dateAndTime(5),':',&
+! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
-                                            dateAndTime(6),':',&
+#if __INTEL_COMPILER >= 1800
-                                            dateAndTime(7)
+ write(6,*) 'Compiled with: ', compiler_version()
- write(6,'(/,a)') ' <<<+-  DAMASK_interface init  -+>>>'
+ write(6,*) 'Compiler options: ', compiler_options()
-#include "compilation_info.f90"
+#else
+ write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
+                                                    ', build date :', __INTEL_COMPILER_BUILD_DATE
+#endif
+
+ write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__
+
+ write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
+ write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
+
  inquire(5, name=wd)                                                                                ! determine inputputfile
  wd = wd(1:scan(wd,'/',back=.true.))
  ierr = CHDIR(wd)
@@ -134,6 +148,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
    debug_info, &
    debug_reset
  use mesh, only: &
+   theMesh, &
    mesh_FEasCP, &
    mesh_element, &
    mesh_node0, &
@@ -141,8 +156,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
    mesh_Ncellnodes, &
    mesh_cellnode, &
    mesh_build_cellnodes, &
-   mesh_build_ipCoordinates, &
+   mesh_build_ipCoordinates
-   FE_Nnodes
  use CPFEM, only: &
    CPFEM_general, &
    CPFEM_init_done, &
@@ -314,7 +328,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
          computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN)                                ! collect and backup Jacobian after convergence
          lastIncConverged = .false.                                                                 ! reset flag
        endif
-       do node = 1,FE_Nnodes(mesh_element(2,cp_en))
+       do node = 1,theMesh%elem%nNodes
          CPnodeID = mesh_element(4_pInt+node,cp_en)
          mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node)
        enddo

src/FEM_utilities.f90

@@ -495,7 +495,6 @@ subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,global
    CHKERRQ(ierr)
    call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
  endif                      
- deallocate(localIndices)
  
 end subroutine utilities_indexActiveSet
 

src/FEM_zoo.f90

@@ -9,11 +9,11 @@ module FEM_Zoo
  private
  integer(pInt), parameter, public:: &
    maxOrder = 5                                                                            !< current max interpolation set at cubic (intended to be arbitrary)
- real(pReal),   dimension(2,3),       private, protected :: &
+ real(pReal),   dimension(2,3),       private, parameter :: &
    triangle    = reshape([-1.0_pReal, -1.0_pReal, &
                            1.0_pReal, -1.0_pReal, &
                           -1.0_pReal,  1.0_pReal], shape=[2,3])
- real(pReal),   dimension(3,4),       private, protected :: &
+ real(pReal),   dimension(3,4),       private, parameter :: &
    tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
                            1.0_pReal, -1.0_pReal, -1.0_pReal, &
                           -1.0_pReal,  1.0_pReal, -1.0_pReal, &

src/commercialFEM_fileList.f90

@@ -12,7 +12,14 @@
 #endif
 #include "math.f90"
 #include "FEsolving.f90"
-#include "mesh.f90"
+#include "element.f90"
+#include "mesh_base.f90"
+#ifdef Abaqus
+#include "mesh_abaqus.f90"
+#endif
+#ifdef Marc4DAMASK
+#include "mesh_marc.f90"
+#endif
 #include "material.f90"
 #include "lattice.f90"
 #include "source_thermal_dissipation.f90"

src/compilation_info.f90

@@ -1,14 +0,0 @@
-! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
-#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
-  write(6,*) 'Compiled with ', compiler_version()
-  write(6,*) 'With options  ', compiler_options()
-#elif defined(__INTEL_COMPILER)
-  write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version ', __INTEL_COMPILER,&
-                                                     ', build date ', __INTEL_COMPILER_BUILD_DATE
-#elif defined(__PGI)
-  write(6,'(a,i4.4,a,i8.8)') ' Compiled with PGI fortran version ', __PGIC__,&
-                                                               '.', __PGIC_MINOR__
-#endif
-write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__
-write(6,*)
-flush(6)

src/constitutive.f90

@@ -56,12 +56,9 @@ subroutine constitutive_init()
    IO_checkAndRewind, &
    IO_open_jobFile_stat, &
    IO_write_jobFile, &
-   IO_write_jobIntFile, &
    IO_timeStamp
  use config, only: &
    config_phase
- use mesh, only: &
-   FE_geomtype
  use config, only: &
    material_Nphase, &
    material_localFileExt, &
@@ -790,8 +787,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
    math_sym33to6, &
    math_mul33x33
  use mesh, only: &
-   mesh_NcpElems, &
+   theMesh
-   mesh_maxNips
  use material, only: &
    phasememberAt, &
    phase_plasticityInstance, &
@@ -842,9 +838,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
    el                                                                                               !< element
  real(pReal),  intent(in) :: &
    subdt                                                                                            !< timestep
- real(pReal),  intent(in), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ real(pReal),  intent(in), dimension(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    subfracArray                                                                                     !< subfraction of timestep
- real(pReal),  intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ real(pReal),  intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    FeArray, &                                                                                       !< elastic deformation gradient
    FpArray                                                                                          !< plastic deformation gradient
  real(pReal),  intent(in), dimension(3,3) :: &
@@ -1004,8 +1000,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
   math_6toSym33, &
   math_mul33x33
  use mesh, only: &
-   mesh_NcpElems, &
+   theMesh
-   mesh_maxNips
  use material, only: &
    phasememberAt, &
    phase_plasticityInstance, &
@@ -1061,7 +1056,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
    constitutive_postResults
  real(pReal),  intent(in), dimension(3,3) :: &
    Fi                                                                                               !< intermediate deformation gradient
- real(pReal),  intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ real(pReal),  intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    FeArray                                                                                          !< elastic deformation gradient
  real(pReal),  intent(in), dimension(6) :: &
    S6                                                                                               !< 2nd Piola Kirchhoff stress (vector notation)

src/crystallite.f90

@@ -12,8 +12,6 @@ module crystallite
  use FEsolving, only:  &
    FEsolving_execElem, &
    FEsolving_execIP
- use mesh, only: &
-   mesh_element
  use material, only: &
    homogenization_Ngrains
  use prec, only: &
@@ -155,10 +153,8 @@ subroutine crystallite_init
    math_inv33, &
    math_mul33x33
  use mesh, only: &
-   mesh_element, &
+   theMesh, &
-   mesh_NcpElems, &
+   mesh_element
-   mesh_maxNips, &
-   mesh_maxNipNeighbors
  use IO, only: &
    IO_timeStamp, &
    IO_stringValue, &
@@ -196,8 +192,8 @@ subroutine crystallite_init
 #include "compilation_info.f90"
 
  cMax = homogenization_maxNgrains
- iMax = mesh_maxNips
+ iMax = theMesh%elem%nIPs
- eMax = mesh_NcpElems
+ eMax = theMesh%nElems
 
 ! ---------------------------------------------------------------------------
 ! ToDo (when working on homogenization): should be 3x3 tensor called S
@@ -333,7 +329,7 @@ subroutine crystallite_init
        case(elasmatrix_ID)
          mySize = 36_pInt
        case(neighboringip_ID,neighboringelement_ID)
-         mySize = mesh_maxNipNeighbors
+         mySize = theMesh%elem%nIPneighbors
        case default
          mySize = 0_pInt
      end select
@@ -415,7 +411,7 @@ subroutine crystallite_init
    write(6,'(a42,1x,i10)') '    # of elements:                       ', eMax
    write(6,'(a42,1x,i10)') 'max # of integration points/element:     ', iMax
    write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
-   write(6,'(a42,1x,i10)') 'max # of neigbours/integration point:    ', mesh_maxNipNeighbors
+   write(6,'(a42,1x,i10)') 'max # of neigbours/integration point:    ', theMesh%elem%nIPneighbors
    write(6,'(a42,1x,i10)') '    # of nonlocal constituents:          ',count(.not. crystallite_localPlasticity)
    flush(6)
  endif
@@ -430,7 +426,7 @@ end subroutine crystallite_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-function crystallite_stress()
+function crystallite_stress(a)
  use prec, only: &
    tol_math_check, &
    dNeq0
@@ -458,10 +454,8 @@ function crystallite_stress()
    math_6toSym33, &
    math_sym33to6
  use mesh, only: &
-   mesh_NcpElems, &
+   theMesh, &
-   mesh_element, &
+   mesh_element
-   mesh_maxNips, &
-   FE_geomtype
  use material, only: &
    homogenization_Ngrains, &
    plasticState, &
@@ -474,7 +468,8 @@ function crystallite_stress()
    constitutive_LiAndItsTangents
 
  implicit none
- logical, dimension(mesh_maxNips,mesh_NcpElems) :: crystallite_stress
+ logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
+ real(pReal), intent(in), optional :: a                                                             !ToDo: for some reason this prevents an internal compiler error in GNU. Very strange
  real(pReal) :: &
    formerSubStep
  integer(pInt) :: &
@@ -541,7 +536,7 @@ function crystallite_stress()
    endIP   = startIP
  else singleRun
    startIP = 1_pInt
-   endIP = mesh_maxNips
+   endIP = theMesh%elem%nIPs
  endif singleRun
 
  NiterationCrystallite = 0_pInt
@@ -727,8 +722,7 @@ subroutine crystallite_stressTangent()
    math_invert2, &
    math_det33
  use mesh, only: &
-   mesh_element, &
+   mesh_element
-   FE_geomtype
  use material, only: &
    homogenization_Ngrains
  use constitutive, only:  &
@@ -929,7 +923,7 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
   math_inv33, &
   math_EulerToR
  use material, only: &
-  material_EulerAngles
+  material_EulerAngles                                                                              ! ToDo: Why stored? We also have crystallite_orientation0
 
  implicit none
  real(pReal), dimension(3,3) :: crystallite_push33ToRef
@@ -960,13 +954,10 @@ function crystallite_postResults(ipc, ip, el)
    inDeg, &
    math_6toSym33
  use mesh, only: &
+   theMesh, &
    mesh_element, &
    mesh_ipVolume, &
-   mesh_maxNipNeighbors, &
+   mesh_ipNeighborhood
-   mesh_ipNeighborhood, &
-   FE_NipNeighbors, &
-   FE_geomtype, &
-   FE_celltype
  use material, only: &
    plasticState, &
    sourceState, &
@@ -1070,14 +1061,14 @@ function crystallite_postResults(ipc, ip, el)
        mySize = 36_pInt
        crystallite_postResults(c+1:c+mySize) = reshape(constitutive_homogenizedC(ipc,ip,el),[mySize])
      case(neighboringelement_ID)
-       mySize = mesh_maxNipNeighbors
+       mySize = theMesh%elem%nIPneighbors
        crystallite_postResults(c+1:c+mySize) = 0.0_pReal
-       forall (n = 1_pInt:FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) &
+       forall (n = 1_pInt:mySize) &
          crystallite_postResults(c+n) = real(mesh_ipNeighborhood(1,n,ip,el),pReal)
      case(neighboringip_ID)
-       mySize = mesh_maxNipNeighbors
+       mySize = theMesh%elem%nIPneighbors
        crystallite_postResults(c+1:c+mySize) = 0.0_pReal
-       forall (n = 1_pInt:FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) &
+       forall (n = 1_pInt:mySize) &
          crystallite_postResults(c+n) = real(mesh_ipNeighborhood(2,n,ip,el),pReal)
    end select
    c = c + mySize
@@ -1754,9 +1745,8 @@ end subroutine integrateStateEuler
 !--------------------------------------------------------------------------------------------------
 subroutine integrateStateAdaptiveEuler()
  use mesh, only: &
-   mesh_element, &
+   theMesh, &
-   mesh_NcpElems, &
+   mesh_element
-   mesh_maxNips
  use material, only: &
    homogenization_Ngrains, &
    plasticState, &
@@ -1780,11 +1770,11 @@ subroutine integrateStateAdaptiveEuler()
    
   ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler
  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
-                        homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+                        homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    residuum_plastic
  real(pReal), dimension(constitutive_source_maxSizeDotState,&
                         maxval(phase_Nsources), &
-                        homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+                        homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    residuum_source
 
 !--------------------------------------------------------------------------------------------------
@@ -1931,8 +1921,7 @@ end subroutine integrateStateRK4
 subroutine integrateStateRKCK45()
  use mesh, only: &
    mesh_element, &
-   mesh_NcpElems, &
+   theMesh
-   mesh_maxNips
  use material, only: &
    homogenization_Ngrains, &
    plasticState, &
@@ -1979,11 +1968,11 @@ subroutine integrateStateRKCK45()
    ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45
 
  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
-                        homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+                        homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    residuum_plastic                                                                         ! relative residuum from evolution in microstructure
  real(pReal), dimension(constitutive_source_maxSizeDotState, &
                         maxval(phase_Nsources), &
-                        homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+                        homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
    residuum_source                                                                   ! relative residuum from evolution in microstructure
 
 
@@ -2128,7 +2117,8 @@ end subroutine nonlocalConvergenceCheck
 !> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
 !--------------------------------------------------------------------------------------------------
 subroutine setConvergenceFlag()
-
+ use mesh, only: &
+   mesh_element
  implicit none
  integer(pInt) :: &
    e, &                                                                                             !< element index in element loop
@@ -2168,7 +2158,8 @@ end subroutine setConvergenceFlag
 !> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
 !--------------------------------------------------------------------------------------------------
 subroutine update_stress(timeFraction)
-
+ use mesh, only: &
+   mesh_element
  implicit none
  real(pReal), intent(in) :: &
    timeFraction
@@ -2200,6 +2191,8 @@ end subroutine update_stress
 !> @brief tbd
 !--------------------------------------------------------------------------------------------------
 subroutine update_dependentState()
+ use mesh, only: &
+   mesh_element
  use constitutive, only: &
    constitutive_dependentState => constitutive_microstructure
 
@@ -2232,6 +2225,8 @@ subroutine update_state(timeFraction)
    sourceState, &
    phase_Nsources, &
    phaseAt, phasememberAt
+ use mesh, only: &
+   mesh_element
 
  implicit none
  real(pReal), intent(in) :: &
@@ -2281,6 +2276,8 @@ subroutine update_dotState(timeFraction)
    sourceState, &
    phaseAt, phasememberAt, &
    phase_Nsources
+ use mesh, only: &
+   mesh_element
  use constitutive, only: &
    constitutive_collectDotState
 
@@ -2334,6 +2331,8 @@ subroutine update_deltaState
    IEEE_arithmetic
  use prec, only: &
    dNeq0
+ use mesh, only: &
+   mesh_element
  use material, only: &
    plasticState, &
    sourceState, &
@@ -2429,6 +2428,8 @@ logical function stateJump(ipc,ip,el)
    sourceState, &
    phase_Nsources, &
    phaseAt, phasememberAt
+ use mesh, only: &
+   mesh_element
  use constitutive, only: &
    constitutive_collectDeltaState
  use math, only: &

src/element.f90

@@ -0,0 +1,921 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Christoph Koords, Max-Planck-Institut für Eisenforschung GmbH
+!--------------------------------------------------------------------------------------------------
+module element
+ use prec, only: &
+   pInt, &
+   pReal
+
+ implicit none
+ private 
+
+!---------------------------------------------------------------------------------------------------
+!> Properties of a single element (the element used in the mesh)
+!---------------------------------------------------------------------------------------------------
+ type, public :: tElement
+   integer(pInt) :: &
+     elemType, &
+     geomType, &                                                                                    ! geometry type (same for same dimension and same number of integration points)
+     cellType, &
+     Nnodes, &
+     Ncellnodes, &
+     NcellnodesPerCell, &
+     nIPs, &
+     nIPneighbors, &                                                                                   ! ToDo: MD: Do all IPs in one element type have the same number of neighbors?
+     maxNnodeAtIP
+   integer(pInt), dimension(:,:), allocatable :: &
+     Cell, &                                                                                        ! intra-element (cell) nodes that constitute a cell
+     NnodeAtIP, &
+     IPneighbor, &
+     cellFace
+   real(pReal), dimension(:,:), allocatable :: &
+ ! center of gravity of the weighted nodes gives the position of the cell node.
+ ! example: face-centered cell node with face nodes 1,2,5,6 to be used in,
+ !          e.g., an 8 node element, would be encoded:
+ !          1, 1, 0, 0, 1, 1, 0, 0
+     cellNodeParentNodeWeights
+   contains
+     procedure :: init => tElement_init
+ end type
+
+ integer(pInt), parameter, private :: &
+   NELEMTYPE = 13_pInt
+
+ integer(pInt), dimension(NelemType), parameter, private :: NNODE = &
+ int([ &
+      3, & ! 2D 3node 1ip
+      6, & ! 2D 6node 3ip
+      4, & ! 2D 4node 4ip
+      8, & ! 2D 8node 9ip
+      8, & ! 2D 8node 4ip
+      !--------------------
+      4, & ! 3D 4node 1ip
+      5, & ! 3D 5node 4ip
+     10, & ! 3D 10node 4ip
+      6, & ! 3D 6node 6ip
+      8, & ! 3D 8node 1ip
+      8, & ! 3D 8node 8ip
+     20, & ! 3D 20node 8ip
+     20  & ! 3D 20node 27ip
+  ],pInt)                                                                                           !< number of nodes that constitute a specific type of element
+
+ integer(pInt), dimension(NelemType), parameter, public :: GEOMTYPE = &
+ int([ &
+      1, & ! 2D 3node 1ip
+      2, & ! 2D 6node 3ip
+      3, & ! 2D 4node 4ip
+      4, & ! 2D 8node 9ip
+      3, & ! 2D 8node 4ip
+      !--------------------
+      5, & ! 3D 4node 1ip
+      6, & ! 3D 5node 4ip
+      6, & ! 3D 10node 4ip
+      7, & ! 3D 6node 6ip
+      8, & ! 3D 8node 1ip
+      9, & ! 3D 8node 8ip
+      9, & ! 3D 20node 8ip
+     10  & ! 3D 20node 27ip
+  ],pInt)                                                                                           !< geometry type of particular element type
+
+ !integer(pInt), dimension(maxval(geomType)), parameter, private :: NCELLNODE = &                   ! Intel 16.0 complains
+ integer(pInt), dimension(10), parameter, private :: NCELLNODE = &
+ int([ &
+      3, &
+      7, &
+      9, &
+     16, &
+      4, &
+     15, &
+     21, &
+      8, &
+     27, &
+     64  &
+  ],pInt)                                                                                           !< number of cell nodes in a specific geometry type
+
+ !integer(pInt), dimension(maxval(geomType)), parameter, private :: NIP = &                         ! Intel 16.0 complains
+ integer(pInt), dimension(10), parameter, private :: NIP = &
+ int([ &
+      1, & 
+      3, & 
+      4, & 
+      9, & 
+      1, & 
+      4, & 
+      6, & 
+      1, & 
+      8, & 
+     27  & 
+  ],pInt)                                                                                           !< number of IPs in a specific geometry type
+
+ !integer(pInt), dimension(maxval(geomType)), parameter, private  :: CELLTYPE = &                   ! Intel 16.0 complains
+ integer(pInt), dimension(10), parameter, private  :: CELLTYPE = &                                  !< cell type that is used by each geometry type
+ int([ &
+      1, & ! 2D 3node
+      2, & ! 2D 4node
+      2, & ! 2D 4node
+      2, & ! 2D 4node
+      3, & ! 3D 4node
+      4, & ! 3D 8node
+      4, & ! 3D 8node
+      4, & ! 3D 8node
+      4, & ! 3D 8node
+      4  & ! 3D 8node
+  ],pInt)
+ 
+ !integer(pInt), dimension(maxval(cellType)), parameter, private :: nIPNeighbor = &                 ! causes problem with Intel 16.0
+ integer(pInt), dimension(4), parameter, private :: NIPNEIGHBOR = &                                 !< 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)
+
+ !integer(pInt), dimension(maxval(cellType)), parameter, private :: NCELLNODESPERCELLFACE = &            
+ integer(pInt), dimension(4), parameter, private :: NCELLNODEPERCELLFACE = &                        !< number of cell nodes in a specific cell type
+ int([ &
+      2, & ! 2D 3node
+      2, & ! 2D 4node
+      3, & ! 3D 4node
+      4  & ! 3D 8node
+  ],pInt)
+
+ !integer(pInt), dimension(maxval(geomType)), parameter, private :: maxNodeAtIP = &                 ! causes problem with Intel 16.0
+ integer(pInt), dimension(10), parameter, private :: maxNnodeAtIP = &                               !< maximum number of parent nodes that belong to an IP for a specific type of element
+ int([ &
+      3, &
+      1, &
+      1, &
+      2, &
+      4, &
+      1, &
+      1, &
+      8, &
+      1, &
+      4  &
+  ],pInt)
+
+
+ !integer(pInt), dimension(maxval(CELLTYPE)), parameter, private  :: NCELLNODEPERCELL = &           ! Intel 16.0 complains
+ integer(pInt), dimension(4), parameter, private :: NCELLNODEPERCELL = &                            !< 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(maxNnodeAtIP(1),nIP(1)), parameter, private :: NnodeAtIP1 = &
+    reshape(int([&
+    1,2,3   &
+    ],pInt),[maxNnodeAtIP(1),nIP(1)])
+
+ integer(pInt), dimension(maxNnodeAtIP(2),nIP(2)), parameter, private :: NnodeAtIP2 = &
+    reshape(int([&
+    1,  &
+    2,  &
+    3   &
+    ],pInt),[maxNnodeAtIP(2),nIP(2)])
+
+ integer(pInt), dimension(maxNnodeAtIP(3),nIP(3)), parameter, private :: NnodeAtIP3 = &
+    reshape(int([&
+    1,  &
+    2,  &
+    4,  &
+    3   &
+    ],pInt),[maxNnodeAtIP(3),nIP(3)])
+
+ integer(pInt), dimension(maxNnodeAtIP(4),nIP(4)), parameter, private :: NnodeAtIP4 = &
+    reshape(int([&
+    1,0,  &
+    1,2,  &
+    2,0,  &
+    1,4,  &
+    0,0,  &
+    2,3,  &
+    4,0,  &
+    3,4,  &
+    3,0   &
+    ],pInt),[maxNnodeAtIP(4),nIP(4)])
+
+ integer(pInt), dimension(maxNnodeAtIP(5),nIP(5)), parameter, private :: NnodeAtIP5 = &
+    reshape(int([&
+    1,2,3,4   &
+    ],pInt),[maxNnodeAtIP(5),nIP(5)])
+
+ integer(pInt), dimension(maxNnodeAtIP(6),nIP(6)), parameter, private :: NnodeAtIP6 = &
+    reshape(int([&
+    1,  &
+    2,  &
+    3,  &
+    4   &
+    ],pInt),[maxNnodeAtIP(6),nIP(6)])
+
+ integer(pInt), dimension(maxNnodeAtIP(7),nIP(7)), parameter, private :: NnodeAtIP7 = &
+    reshape(int([&
+    1,  &
+    2,  &
+    3,  &
+    4,  &
+    5,  &
+    6   &
+    ],pInt),[maxNnodeAtIP(7),nIP(7)])
+
+ integer(pInt), dimension(maxNnodeAtIP(8),nIP(8)), parameter, private :: NnodeAtIP8 = &
+    reshape(int([&
+    1,2,3,4,5,6,7,8   &
+    ],pInt),[maxNnodeAtIP(8),nIP(8)])
+
+ integer(pInt), dimension(maxNnodeAtIP(9),nIP(9)), parameter, private :: NnodeAtIP9 = &
+    reshape(int([&
+    1,  &
+    2,  &
+    4,  &
+    3,  &
+    5,  &
+    6,  &
+    8,  &
+    7   &
+    ],pInt),[maxNnodeAtIP(9),nIP(9)])
+
+ integer(pInt), dimension(maxNnodeAtIP(10),nIP(10)), parameter, private :: NnodeAtIP10 = &
+    reshape(int([&
+    1,0, 0,0,  &
+    1,2, 0,0,  &
+    2,0, 0,0,  &
+    1,4, 0,0,  &
+    1,3, 2,4,  &
+    2,3, 0,0,  &
+    4,0, 0,0,  &
+    3,4, 0,0,  &
+    3,0, 0,0,  &
+    1,5, 0,0,  &
+    1,6, 2,5,  &
+    2,6, 0,0,  &
+    1,8, 4,5,  &
+    0,0, 0,0,  &
+    2,7, 3,6,  &
+    4,8, 0,0,  &
+    3,8, 4,7,  &
+    3,7, 0,0,  &
+    5,0, 0,0,  &
+    5,6, 0,0,  &
+    6,0, 0,0,  &
+    5,8, 0,0,  &
+    5,7, 6,8,  &
+    6,7, 0,0,  &
+    8,0, 0,0,  &
+    7,8, 0,0,  &
+    7,0, 0,0   &
+    ],pInt),[maxNnodeAtIP(10),nIP(10)])
+
+ ! *** FE_ipNeighbor ***
+ ! is a 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.
+ ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located.
+
+
+ integer(pInt), dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = &
+    reshape(int([&
+    -2,-3,-1   &
+    ],pInt),[nIPneighbor(cellType(1)),nIP(1)])
+
+ integer(pInt), dimension(nIPneighbor(cellType(2)),nIP(2)), parameter, private :: IPneighbor2 = &
+     reshape(int([&
+      2,-3, 3,-1,  &
+     -2, 1, 3,-1,  &
+      2,-3,-2, 1   &
+    ],pInt),[nIPneighbor(cellType(2)),nIP(2)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(3)),nIP(3)), parameter, private :: IPneighbor3 = &
+     reshape(int([&
+      2,-4, 3,-1,  &
+     -2, 1, 4,-1,  &
+      4,-4,-3, 1,  &
+     -2, 3,-3, 2   &
+    ],pInt),[nIPneighbor(cellType(3)),nIP(3)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(4)),nIP(4)), parameter, private :: IPneighbor4 = &
+     reshape(int([&
+      2,-4, 4,-1,  &
+      3, 1, 5,-1,  &
+     -2, 2, 6,-1,  &
+      5,-4, 7, 1,  &
+      6, 4, 8, 2,  &
+     -2, 5, 9, 3,  &
+      8,-4,-3, 4,  &
+      9, 7,-3, 5,  &
+     -2, 8,-3, 6   &
+    ],pInt),[nIPneighbor(cellType(4)),nIP(4)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(5)),nIP(5)), parameter, private :: IPneighbor5 = &
+     reshape(int([&
+     -1,-2,-3,-4   &
+    ],pInt),[nIPneighbor(cellType(5)),nIP(5)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(6)),nIP(6)), parameter, private :: IPneighbor6 = &
+     reshape(int([&
+      2,-4, 3,-2, 4,-1,  &
+     -2, 1, 3,-2, 4,-1,  &
+      2,-4,-3, 1, 4,-1,  &
+      2,-4, 3,-2,-3, 1   &
+    ],pInt),[nIPneighbor(cellType(6)),nIP(6)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(7)),nIP(7)), parameter, private :: IPneighbor7 = &
+     reshape(int([&
+      2,-4, 3,-2, 4,-1,  &
+     -3, 1, 3,-2, 5,-1,  &
+      2,-4,-3, 1, 6,-1,  &
+      5,-4, 6,-2,-5, 1,  &
+     -3, 4, 6,-2,-5, 2,  &
+      5,-4,-3, 4,-5, 3   &
+    ],pInt),[nIPneighbor(cellType(7)),nIP(7)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(8)),nIP(8)), parameter, private :: IPneighbor8 = &
+     reshape(int([&
+     -3,-5,-4,-2,-6,-1   &
+    ],pInt),[nIPneighbor(cellType(8)),nIP(8)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(9)),nIP(9)), parameter, private :: IPneighbor9 = &
+     reshape(int([&
+      2,-5, 3,-2, 5,-1,  &
+     -3, 1, 4,-2, 6,-1,  &
+      4,-5,-4, 1, 7,-1,  &
+     -3, 3,-4, 2, 8,-1,  &
+      6,-5, 7,-2,-6, 1,  &
+     -3, 5, 8,-2,-6, 2,  &
+      8,-5,-4, 5,-6, 3,  &
+     -3, 7,-4, 6,-6, 4   &
+    ],pInt),[nIPneighbor(cellType(9)),nIP(9)])
+ 
+ integer(pInt), dimension(nIPneighbor(cellType(10)),nIP(10)), parameter, private :: IPneighbor10 = &
+     reshape(int([&
+      2,-5, 4,-2,10,-1,  &
+      3, 1, 5,-2,11,-1,  &
+     -3, 2, 6,-2,12,-1,  &
+      5,-5, 7, 1,13,-1,  &
+      6, 4, 8, 2,14,-1,  &
+     -3, 5, 9, 3,15,-1,  &
+      8,-5,-4, 4,16,-1,  &
+      9, 7,-4, 5,17,-1,  &
+     -3, 8,-4, 6,18,-1,  &
+     11,-5,13,-2,19, 1,  &
+     12,10,14,-2,20, 2,  &
+     -3,11,15,-2,21, 3,  &
+     14,-5,16,10,22, 4,  &
+     15,13,17,11,23, 5,  &
+     -3,14,18,12,24, 6,  &
+     17,-5,-4,13,25, 7,  &
+     18,16,-4,14,26, 8,  &
+     -3,17,-4,15,27, 9,  &
+     20,-5,22,-2,-6,10,  &
+     21,19,23,-2,-6,11,  &
+     -3,20,24,-2,-6,12,  &
+     23,-5,25,19,-6,13,  &
+     24,22,26,20,-6,14,  &
+     -3,23,27,21,-6,15,  &
+     26,-5,-4,22,-6,16,  &
+     27,25,-4,23,-6,17,  &
+     -3,26,-4,24,-6,18   &
+    ],pInt),[nIPneighbor(cellType(10)),nIP(10)])
+
+
+ real(pReal), dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
+    reshape(real([&
+    1, 0, 0,  &
+    0, 1, 0,  &
+    0, 0, 1   &
+    ],pReal),[nNode(1),NcellNode(geomType(1))]) ! 2D 3node 1ip
+
+ real(pReal), dimension(nNode(2),NcellNode(geomType(2))), parameter :: cellNodeParentNodeWeights2 = &
+     reshape(real([&
+     1, 0, 0, 0, 0, 0,  &
+     0, 1, 0, 0, 0, 0,  &
+     0, 0, 1, 0, 0, 0,  &
+     0, 0, 0, 1, 0, 0,  &
+     0, 0, 0, 0, 1, 0,  &
+     0, 0, 0, 0, 0, 1,  &
+     1, 1, 1, 2, 2, 2   &
+    ],pReal),[nNode(2),NcellNode(geomType(2))]) ! 2D 6node 3ip
+
+ real(pReal), dimension(nNode(3),NcellNode(geomType(3))), parameter :: cellNodeParentNodeWeights3 = &
+    reshape(real([&
+    1, 0, 0, 0,  &
+    0, 1, 0, 0,  &
+    0, 0, 1, 0,  &
+    0, 0, 0, 1,  &
+    1, 1, 0, 0,  &
+    0, 1, 1, 0,  &
+    0, 0, 1, 1,  &
+    1, 0, 0, 1,  &
+    1, 1, 1, 1   &
+    ],pReal),[nNode(3),NcellNode(geomType(3))]) ! 2D 6node 3ip
+
+ real(pReal), dimension(nNode(4),NcellNode(geomType(4))), parameter :: cellNodeParentNodeWeights4 = &
+    reshape(real([&
+    1, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 1, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 1, 0, 0, 0, 0,  &
+    1, 0, 0, 0, 2, 0, 0, 0,  &
+    0, 1, 0, 0, 2, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 2, 0, 0,  &
+    0, 0, 1, 0, 0, 2, 0, 0,  &
+    0, 0, 1, 0, 0, 0, 2, 0,  &
+    0, 0, 0, 1, 0, 0, 2, 0,  &
+    0, 0, 0, 1, 0, 0, 0, 2,  &
+    1, 0, 0, 0, 0, 0, 0, 2,  &
+    4, 1, 1, 1, 8, 2, 2, 8,  &
+    1, 4, 1, 1, 8, 8, 2, 2,  &
+    1, 1, 4, 1, 2, 8, 8, 2,  &
+    1, 1, 1, 4, 2, 2, 8, 8   &
+    ],pReal),[nNode(4),NcellNode(geomType(4))])  ! 2D 8node 9ip
+
+ real(pReal), dimension(nNode(5),NcellNode(geomType(5))), parameter :: cellNodeParentNodeWeights5 = &
+    reshape(real([&
+    1, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 1, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 1, 0, 0, 0, 0,  &
+    0, 0, 0, 0, 1, 0, 0, 0,  &
+    0, 0, 0, 0, 0, 1, 0, 0,  &
+    0, 0, 0, 0, 0, 0, 1, 0,  &
+    0, 0, 0, 0, 0, 0, 0, 1,  &
+    1, 1, 1, 1, 2, 2, 2, 2   &
+    ],pReal),[nNode(5),NcellNode(geomType(5))])  ! 2D 8node 4ip
+
+ real(pReal), dimension(nNode(6),NcellNode(geomType(6))), parameter :: cellNodeParentNodeWeights6 = &
+    reshape(real([&
+    1, 0, 0, 0,  &
+    0, 1, 0, 0,  &
+    0, 0, 1, 0,  &
+    0, 0, 0, 1   &
+    ],pReal),[nNode(6),NcellNode(geomType(6))]) ! 3D 4node 1ip
+
+ real(pReal), dimension(nNode(7),NcellNode(geomType(7))), parameter :: cellNodeParentNodeWeights7 = &
+    reshape(real([&
+    1, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0,  &
+    0, 0, 1, 0, 0,  &
+    0, 0, 0, 1, 0,  &
+    1, 1, 0, 0, 0,  &
+    0, 1, 1, 0, 0,  &
+    1, 0, 1, 0, 0,  &
+    1, 0, 0, 1, 0,  &
+    0, 1, 0, 1, 0,  &
+    0, 0, 1, 1, 0,  &
+    1, 1, 1, 0, 0,  &
+    1, 1, 0, 1, 0,  &
+    0, 1, 1, 1, 0,  &
+    1, 0, 1, 1, 0,  &
+    0, 0, 0, 0, 1   &
+    ],pReal),[nNode(7),NcellNode(geomType(7))])  ! 3D 5node 4ip
+
+ real(pReal), dimension(nNode(8),NcellNode(geomType(8))), parameter :: cellNodeParentNodeWeights8 = &
+    reshape(real([&
+    1, 0, 0, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 1, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 1, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 0, 1, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 0, 0, 1, 0, 0, 0, 0,  &
+    0, 0, 0, 0, 0, 0, 1, 0, 0, 0,  &
+    0, 0, 0, 0, 0, 0, 0, 1, 0, 0,  &
+    0, 0, 0, 0, 0, 0, 0, 0, 1, 0,  &
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 1,  &
+    1, 1, 1, 0, 2, 2, 2, 0, 0, 0,  &
+    1, 1, 0, 1, 2, 0, 0, 2, 2, 0,  &
+    0, 1, 1, 1, 0, 2, 0, 0, 2, 2,  &
+    1, 0, 1, 1, 0, 0, 2, 2, 0, 2,  &
+    3, 3, 3, 3, 4, 4, 4, 4, 4, 4   &
+    ],pReal),[nNode(8),NcellNode(geomType(8))])  ! 3D 10node 4ip
+
+ real(pReal), dimension(nNode(9),NcellNode(geomType(9))), parameter :: cellNodeParentNodeWeights9 = &
+    reshape(real([&
+    1, 0, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 0,  &
+    0, 0, 1, 0, 0, 0,  &
+    0, 0, 0, 1, 0, 0,  &
+    0, 0, 0, 0, 1, 0,  &
+    0, 0, 0, 0, 0, 1,  &
+    1, 1, 0, 0, 0, 0,  &
+    0, 1, 1, 0, 0, 0,  &
+    1, 0, 1, 0, 0, 0,  &
+    1, 0, 0, 1, 0, 0,  &
+    0, 1, 0, 0, 1, 0,  &
+    0, 0, 1, 0, 0, 1,  &
+    0, 0, 0, 1, 1, 0,  &
+    0, 0, 0, 0, 1, 1,  &
+    0, 0, 0, 1, 0, 1,  &
+    1, 1, 1, 0, 0, 0,  &
+    1, 1, 0, 1, 1, 0,  &
+    0, 1, 1, 0, 1, 1,  &
+    1, 0, 1, 1, 0, 1,  &
+    0, 0, 0, 1, 1, 1,  &
+    1, 1, 1, 1, 1, 1   &
+    ],pReal),[nNode(9),NcellNode(geomType(9))])  ! 3D 6node 6ip
+
+ real(pReal), dimension(nNode(10),NcellNode(geomType(10))), parameter :: cellNodeParentNodeWeights10 = &
+    reshape(real([&
+    1, 0, 0, 0, 0, 0, 0, 0,  &
+    0, 1, 0, 0, 0, 0, 0, 0,  &
+    0, 0, 1, 0, 0, 0, 0, 0,  &
+    0, 0, 0, 1, 0, 0, 0, 0,  &
+    0, 0, 0, 0, 1, 0, 0, 0,  &
+    0, 0, 0, 0, 0, 1, 0, 0,  &
+    0, 0, 0, 0, 0, 0, 1, 0,  &
+    0, 0, 0, 0, 0, 0, 0, 1   &
+    ],pReal),[nNode(10),NcellNode(geomType(10))])  ! 3D 8node 1ip
+ 
+ real(pReal), dimension(nNode(11),NcellNode(geomType(11))), parameter :: cellNodeParentNodeWeights11 = &
+    reshape(real([&
+    1, 0, 0, 0,  0, 0, 0, 0,  &   !
+    0, 1, 0, 0,  0, 0, 0, 0,  &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  &   !
+    0, 0, 0, 0,  1, 0, 0, 0,  &   !  5
+    0, 0, 0, 0,  0, 1, 0, 0,  &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  &   !
+    0, 0, 0, 0,  0, 0, 0, 1,  &   !
+    1, 1, 0, 0,  0, 0, 0, 0,  &   !
+    0, 1, 1, 0,  0, 0, 0, 0,  &   ! 10
+    0, 0, 1, 1,  0, 0, 0, 0,  &   !
+    1, 0, 0, 1,  0, 0, 0, 0,  &   !
+    1, 0, 0, 0,  1, 0, 0, 0,  &   !
+    0, 1, 0, 0,  0, 1, 0, 0,  &   !
+    0, 0, 1, 0,  0, 0, 1, 0,  &   ! 15
+    0, 0, 0, 1,  0, 0, 0, 1,  &   !
+    0, 0, 0, 0,  1, 1, 0, 0,  &   !
+    0, 0, 0, 0,  0, 1, 1, 0,  &   !
+    0, 0, 0, 0,  0, 0, 1, 1,  &   !
+    0, 0, 0, 0,  1, 0, 0, 1,  &   ! 20
+    1, 1, 1, 1,  0, 0, 0, 0,  &   !
+    1, 1, 0, 0,  1, 1, 0, 0,  &   !
+    0, 1, 1, 0,  0, 1, 1, 0,  &   !
+    0, 0, 1, 1,  0, 0, 1, 1,  &   !
+    1, 0, 0, 1,  1, 0, 0, 1,  &   ! 25
+    0, 0, 0, 0,  1, 1, 1, 1,  &   !
+    1, 1, 1, 1,  1, 1, 1, 1   &   !
+    ],pReal),[nNode(11),NcellNode(geomType(11))])  ! 3D 8node 8ip
+ 
+ real(pReal), dimension(nNode(12),NcellNode(geomType(12))), parameter :: cellNodeParentNodeWeights12 = &
+    reshape(real([&
+    1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !  5
+    0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   ! 10
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  1, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 1, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 1, 0, &   ! 15
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 1, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0, &   ! 20
+    1, 1, 1, 1,  0, 0, 0, 0,  2, 2, 2, 2,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    1, 1, 0, 0,  1, 1, 0, 0,  2, 0, 0, 0,  2, 0, 0, 0,  2, 2, 0, 0, &   !
+    0, 1, 1, 0,  0, 1, 1, 0,  0, 2, 0, 0,  0, 2, 0, 0,  0, 2, 2, 0, &   !
+    0, 0, 1, 1,  0, 0, 1, 1,  0, 0, 2, 0,  0, 0, 2, 0,  0, 0, 2, 2, &   !
+    1, 0, 0, 1,  1, 0, 0, 1,  0, 0, 0, 2,  0, 0, 0, 2,  2, 0, 0, 2, &   ! 25
+    0, 0, 0, 0,  1, 1, 1, 1,  0, 0, 0, 0,  2, 2, 2, 2,  0, 0, 0, 0, &   !
+    3, 3, 3, 3,  3, 3, 3, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4  &   !
+    ],pReal),[nNode(12),NcellNode(geomType(12))])  ! 3D 20node 8ip
+ 
+ real(pReal), dimension(nNode(13),NcellNode(geomType(13))), parameter :: cellNodeParentNodeWeights13 = &
+    reshape(real([&
+    1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !  5
+    0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    1, 0, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 1, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   ! 10
+    0, 1, 0, 0,  0, 0, 0, 0,  0, 2, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  0, 2, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 2, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 2, 0,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 2,  0, 0, 0, 0,  0, 0, 0, 0, &   ! 15
+    1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 2,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0, &   !
+    0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 2, 0, 0, &   !
+    0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 2, 0, &   !
+    0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 2, &   ! 20
+    0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 2, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 2, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 2, &   !
+    0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0,  0, 0, 0, 0, &   ! 25
+    0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  2, 0, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 1, 0, 0,  0, 0, 0, 0,  0, 2, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 2, 0, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 1, 0,  0, 0, 0, 0,  0, 0, 2, 0,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 2, 0,  0, 0, 0, 0, &   ! 30
+    0, 0, 0, 0,  0, 0, 0, 1,  0, 0, 0, 0,  0, 0, 0, 2,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  1, 0, 0, 0,  0, 0, 0, 0,  0, 0, 0, 2,  0, 0, 0, 0, &   !
+    4, 1, 1, 1,  0, 0, 0, 0,  8, 2, 2, 8,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    1, 4, 1, 1,  0, 0, 0, 0,  8, 8, 2, 2,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    1, 1, 4, 1,  0, 0, 0, 0,  2, 8, 8, 2,  0, 0, 0, 0,  0, 0, 0, 0, &   ! 35
+    1, 1, 1, 4,  0, 0, 0, 0,  2, 2, 8, 8,  0, 0, 0, 0,  0, 0, 0, 0, &   !
+    4, 1, 0, 0,  1, 1, 0, 0,  8, 0, 0, 0,  2, 0, 0, 0,  8, 2, 0, 0, &   !
+    1, 4, 0, 0,  1, 1, 0, 0,  8, 0, 0, 0,  2, 0, 0, 0,  2, 8, 0, 0, &   !
+    0, 4, 1, 0,  0, 1, 1, 0,  0, 8, 0, 0,  0, 2, 0, 0,  0, 8, 2, 0, &   !
+    0, 1, 4, 0,  0, 1, 1, 0,  0, 8, 0, 0,  0, 2, 0, 0,  0, 2, 8, 0, &   ! 40
+    0, 0, 4, 1,  0, 0, 1, 1,  0, 0, 8, 0,  0, 0, 2, 0,  0, 0, 8, 2, &   !
+    0, 0, 1, 4,  0, 0, 1, 1,  0, 0, 8, 0,  0, 0, 2, 0,  0, 0, 2, 8, &   !
+    1, 0, 0, 4,  1, 0, 0, 1,  0, 0, 0, 8,  0, 0, 0, 2,  2, 0, 0, 8, &   !
+    4, 0, 0, 1,  1, 0, 0, 1,  0, 0, 0, 8,  0, 0, 0, 2,  8, 0, 0, 2, &   !
+    1, 1, 0, 0,  4, 1, 0, 0,  2, 0, 0, 0,  8, 0, 0, 0,  8, 2, 0, 0, &   ! 45
+    1, 1, 0, 0,  1, 4, 0, 0,  2, 0, 0, 0,  8, 0, 0, 0,  2, 8, 0, 0, &   !
+    0, 1, 1, 0,  0, 4, 1, 0,  0, 2, 0, 0,  0, 8, 0, 0,  0, 8, 2, 0, &   !
+    0, 1, 1, 0,  0, 1, 4, 0,  0, 2, 0, 0,  0, 8, 0, 0,  0, 2, 8, 0, &   !
+    0, 0, 1, 1,  0, 0, 4, 1,  0, 0, 2, 0,  0, 0, 8, 0,  0, 0, 8, 2, &   !
+    0, 0, 1, 1,  0, 0, 1, 4,  0, 0, 2, 0,  0, 0, 8, 0,  0, 0, 2, 8, &   ! 50
+    1, 0, 0, 1,  1, 0, 0, 4,  0, 0, 0, 2,  0, 0, 0, 8,  2, 0, 0, 8, &   !
+    1, 0, 0, 1,  4, 0, 0, 1,  0, 0, 0, 2,  0, 0, 0, 8,  8, 0, 0, 2, &   !
+    0, 0, 0, 0,  4, 1, 1, 1,  0, 0, 0, 0,  8, 2, 2, 8,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  1, 4, 1, 1,  0, 0, 0, 0,  8, 8, 2, 2,  0, 0, 0, 0, &   !
+    0, 0, 0, 0,  1, 1, 4, 1,  0, 0, 0, 0,  2, 8, 8, 2,  0, 0, 0, 0, &   ! 55
+    0, 0, 0, 0,  1, 1, 1, 4,  0, 0, 0, 0,  2, 2, 8, 8,  0, 0, 0, 0, &   !
+   24, 8, 4, 8,  8, 4, 3, 4, 32,12,12,32, 12, 4, 4,12, 32,12, 4,12, &   !
+    8,24, 8, 4,  4, 8, 4, 3, 32,32,12,12, 12,12, 4, 4, 12,32,12, 4, &   !
+    4, 8,24, 8,  3, 4, 8, 4, 12,32,32,12,  4,12,12, 4,  4,12,32,12, &   !
+    8, 4, 8,24,  4, 3, 4, 8, 12,12,32,32,  4, 4,12,12, 12, 4,12,32, &   ! 60
+    8, 4, 3, 4, 24, 8, 4, 8, 12, 4, 4,12, 32,12,12,32, 32,12, 4,12, &   !
+    4, 8, 4, 3,  8,24, 8, 4, 12,12, 4, 4, 32,32,12,12, 12,32,12, 4, &   !
+    3, 4, 8, 4,  4, 8,24, 8,  4,12,12, 4, 12,32,32,12,  4,12,32,12, &   !
+    4, 3, 4, 8,  8, 4, 8,24,  4, 4,12,12, 12,12,32,32, 12, 4,12,32  &   !
+    ],pReal),[nNode(13),NcellNode(geomType(13))])  ! 3D 20node 27ip
+
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)), parameter :: CELL1 = &
+    reshape(int([&
+    1,2,3   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(1)),NIP(1)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)), parameter :: CELL2 = &
+    reshape(int([&
+    1, 4, 7, 6,   &
+    2, 5, 7, 4,   &
+    3, 6, 7, 5    &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(2)),NIP(2)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)), parameter :: CELL3 = &
+    reshape(int([&
+    1, 5, 9, 8,   &
+    5, 2, 6, 9,   &
+    8, 9, 7, 4,   &
+    9, 6, 3, 7    &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(3)),NIP(3)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)), parameter :: CELL4 = &
+    reshape(int([&
+    1, 5,13,12,   &
+    5, 6,14,13,   &
+    6, 2, 7,14,   &
+   12,13,16,11,   &
+   13,14,15,16,   &
+   14, 7, 8,15,   &
+   11,16,10, 4,   &
+   16,15, 9,10,   &
+   15, 8, 3, 9    &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(4)),NIP(4)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)), parameter :: CELL5 = &
+    reshape(int([&
+    1, 2, 3, 4   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(5)),NIP(5)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)), parameter :: CELL6 = &
+    reshape(int([&
+    1, 5,11, 7, 8,12,15,14,  &
+    5, 2, 6,11,12, 9,13,15,  &
+    7,11, 6, 3,14,15,13,10,  &
+    8,12,15, 4, 4, 9,13,10   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(6)),NIP(6)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)), parameter :: CELL7 = &
+    reshape(int([&
+    1, 7,16, 9,10,17,21,19,  &
+    7, 2, 8,16,17,11,18,21,  &
+    9,16, 8, 3,19,21,18,12,  &
+   10,17,21,19, 4,13,20,15,  &
+   17,11,18,21,13, 5,14,20,  &
+   19,21,18,12,15,20,14, 6   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(7)),NIP(7)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)), parameter :: CELL8 = &
+    reshape(int([&
+    1, 2, 3, 4, 5, 6, 7, 8   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(8)),NIP(8)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)), parameter :: CELL9 = &
+    reshape(int([&
+    1, 9,21,12,13,22,27,25,  &
+    9, 2,10,21,22,14,23,27,  &
+   12,21,11, 4,25,27,24,16,  &
+   21,10, 3,11,27,23,15,24,  &
+   13,22,27,25, 5,17,26,20,  &
+   22,14,23,27,17, 6,18,26,  &
+   25,27,24,16,20,26,19, 8,  &
+   27,23,15,24,26,18, 7,19   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(9)),NIP(9)])
+
+ integer(pInt), dimension(NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)), parameter :: CELL10 = &
+    reshape(int([&
+    1, 9,33,16,17,37,57,44,  &
+    9,10,34,33,37,38,58,57,  &
+   10, 2,11,34,38,18,39,58,  &
+   16,33,36,15,44,57,60,43,  &
+   33,34,35,36,57,58,59,60,  &
+   34,11,12,35,58,39,40,59,  &
+   15,36,14, 4,43,60,42,20,  &
+   36,35,13,14,60,59,41,42,  &
+   35,12, 3,13,59,40,19,41,  &
+   17,37,57,44,21,45,61,52,  &
+   37,38,58,57,45,46,62,61,  &
+   38,18,39,58,46,22,47,62,  &
+   44,57,60,43,52,61,64,51,  &
+   57,58,59,60,61,62,63,64,  &
+   58,39,40,59,62,47,48,63,  &
+   43,60,42,20,51,64,50,24,  &
+   60,59,41,42,64,63,49,50,  &
+   59,40,19,41,63,48,23,49,  &
+   21,45,61,52, 5,25,53,32,  &
+   45,46,62,61,25,26,54,53,  &
+   46,22,47,62,26, 6,27,54,  &
+   52,61,64,51,32,53,56,31,  &
+   61,62,63,64,53,54,55,56,  &
+   62,47,48,63,54,27,28,55,  &
+   51,64,50,24,31,56,30, 8,  &
+   64,63,49,50,56,55,29,30,  &
+   63,48,23,49,55,28, 7,29   &
+    ],pInt),[NCELLNODEPERCELL(CELLTYPE(10)),NIP(10)])
+
+
+  integer(pInt), dimension(NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)), parameter :: CELLFACE1 = &
+    reshape(int([&
+    2,3,  &
+    3,1,  &
+    1,2   &
+    ],pInt),[NCELLNODEPERCELLFACE(1),NIPNEIGHBOR(1)])                   ! 2D 3node, VTK_TRIANGLE (5)
+
+ integer(pInt), dimension(NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)), parameter :: CELLFACE2 = &
+    reshape(int([&
+    2,3,  &
+    4,1,  &
+    3,4,  &
+    1,2   &
+    ],pInt),[NCELLNODEPERCELLFACE(2),NIPNEIGHBOR(2)])              ! 2D 4node, VTK_QUAD (9)
+
+ integer(pInt), dimension(NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)), parameter :: CELLFACE3 = &
+    reshape(int([&
+    1,3,2,  &
+    1,2,4,  &
+    2,3,4,  &
+    1,4,3   &
+    ],pInt),[NCELLNODEPERCELLFACE(3),NIPNEIGHBOR(3)])                     ! 3D 4node, VTK_TETRA (10)
+
+ integer(pInt), dimension(NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)), parameter :: CELLFACE4 = &
+    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),[NCELLNODEPERCELLFACE(4),NIPNEIGHBOR(4)])              ! 3D 8node, VTK_HEXAHEDRON (12)
+
+
+contains
+
+ subroutine tElement_init(self,elemType)
+   implicit none
+   class(tElement) :: self
+   integer(pInt), intent(in) :: elemType
+   self%elemType = elemType
+
+   self%Nnodes     = Nnode    (self%elemType)
+   self%geomType   = geomType (self%elemType)
+   select case (self%elemType)
+     case(1_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights1
+     case(2_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights2
+     case(3_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights3
+     case(4_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights4
+     case(5_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights5
+     case(6_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights6
+     case(7_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights7
+     case(8_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights8
+     case(9_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights9
+     case(10_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights10
+     case(11_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights11
+     case(12_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights12
+     case(13_pInt)
+       self%cellNodeParentNodeWeights  = cellNodeParentNodeWeights13
+     case default
+       print*, 'Mist'
+  end select
+  
+
+   self%NcellNodes   = NcellNode    (self%geomType)
+   self%maxNnodeAtIP = maxNnodeAtIP (self%geomType) 
+   self%nIPs         = nIP          (self%geomType)
+   self%cellType     = cellType     (self%geomType)
+
+
+   select case (self%geomType)
+     case(1_pInt)
+       self%NnodeAtIP   = NnodeAtIP1
+       self%IPneighbor  = IPneighbor1
+       self%cell        = CELL1
+     case(2_pInt)
+       self%NnodeAtIP   = NnodeAtIP2
+       self%IPneighbor  = IPneighbor2
+       self%cell        = CELL2
+     case(3_pInt)
+       self%NnodeAtIP   = NnodeAtIP3
+       self%IPneighbor  = IPneighbor3
+       self%cell        = CELL3
+     case(4_pInt)
+       self%NnodeAtIP   = NnodeAtIP4
+       self%IPneighbor  = IPneighbor4
+       self%cell        = CELL4
+     case(5_pInt)
+       self%NnodeAtIP   = NnodeAtIP5
+       self%IPneighbor  = IPneighbor5
+       self%cell        = CELL5
+     case(6_pInt)
+       self%NnodeAtIP   = NnodeAtIP6
+       self%IPneighbor  = IPneighbor6
+       self%cell        = CELL6
+     case(7_pInt)
+       self%NnodeAtIP   = NnodeAtIP7
+       self%IPneighbor  = IPneighbor7
+       self%cell        = CELL7
+     case(8_pInt)
+       self%NnodeAtIP   = NnodeAtIP8
+       self%IPneighbor  = IPneighbor8
+       self%cell        = CELL8
+     case(9_pInt)
+       self%NnodeAtIP   = NnodeAtIP9
+       self%IPneighbor  = IPneighbor9
+       self%cell        = CELL9
+     case(10_pInt)
+       self%NnodeAtIP   = NnodeAtIP10
+       self%IPneighbor  = IPneighbor10
+       self%cell        = CELL10
+  end select
+   self%NcellNodesPerCell = NCELLNODEPERCELL(self%cellType)
+ 
+  select case(self%cellType)
+    case(1_pInt)
+      self%cellFace = CELLFACE1
+    case(2_pInt)
+      self%cellFace = CELLFACE2
+    case(3_pInt)
+      self%cellFace = CELLFACE3
+    case(4_pInt)
+      self%cellFace = CELLFACE4
+   end select
+   
+   self%nIPneighbors = size(self%IPneighbor,1)
+   
+   write(6,'(/,a)')   ' <<<+-  element_init  -+>>>'
+   
+   write(6,*)' element type        ',self%elemType
+   write(6,*)'   geom type         ',self%geomType
+   write(6,*)'   cell type         ',self%cellType
+   write(6,*)'   # node            ',self%Nnodes
+   write(6,*)'   # IP              ',self%nIPs
+   write(6,*)'   # cellnode        ',self%Ncellnodes
+   write(6,*)'   # cellnode/cell   ',self%NcellnodesPerCell
+   write(6,*)'   # IP neighbor     ',self%nIPneighbors
+   write(6,*)'   max # node at IP  ',self%maxNnodeAtIP
+   
+ end subroutine tElement_init
+
+end module element

src/homogenization.f90

@@ -71,11 +71,8 @@ subroutine homogenization_init
    debug_e, &
    debug_g
  use mesh, only: &
-   mesh_maxNips, &
+   theMesh, &
-   mesh_NcpElems, &
+   mesh_element
-   mesh_element, &
-   FE_Nips, &
-   FE_geomtype
  use constitutive, only: &
    constitutive_plasticity_maxSizePostResults, &
    constitutive_source_maxSizePostResults
@@ -244,20 +241,20 @@ subroutine homogenization_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate and initialize global variables
- allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems),       source=0.0_pReal)
+ allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems),       source=0.0_pReal)
- allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems),             source=0.0_pReal)
+ allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems),             source=0.0_pReal)
- materialpoint_F0 = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems)                          ! initialize to identity
+ materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems)                          ! initialize to identity
- allocate(materialpoint_F(3,3,mesh_maxNips,mesh_NcpElems),              source=0.0_pReal)
+ allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems),              source=0.0_pReal)
  materialpoint_F = materialpoint_F0                                                                 ! initialize to identity
- allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems),          source=0.0_pReal)
+ allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems),          source=0.0_pReal)
- allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems),           source=0.0_pReal)
+ allocate(materialpoint_subF(3,3,theMesh%elem%nIPs,theMesh%nElems),           source=0.0_pReal)
- allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems),              source=0.0_pReal)
+ allocate(materialpoint_P(3,3,theMesh%elem%nIPs,theMesh%nElems),              source=0.0_pReal)
- allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems),            source=0.0_pReal)
+ allocate(materialpoint_subFrac(theMesh%elem%nIPs,theMesh%nElems),            source=0.0_pReal)
- allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems),            source=0.0_pReal)
+ allocate(materialpoint_subStep(theMesh%elem%nIPs,theMesh%nElems),            source=0.0_pReal)
- allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems),              source=0.0_pReal)
+ allocate(materialpoint_subdt(theMesh%elem%nIPs,theMesh%nElems),              source=0.0_pReal)
- allocate(materialpoint_requested(mesh_maxNips,mesh_NcpElems),          source=.false.)
+ allocate(materialpoint_requested(theMesh%elem%nIPs,theMesh%nElems),          source=.false.)
- allocate(materialpoint_converged(mesh_maxNips,mesh_NcpElems),          source=.true.)
+ allocate(materialpoint_converged(theMesh%elem%nIPs,theMesh%nElems),          source=.true.)
- allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems),     source=.true.)
+ allocate(materialpoint_doneAndHappy(2,theMesh%elem%nIPs,theMesh%nElems),     source=.true.)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate and initialize global state and postresutls variables
@@ -277,7 +274,7 @@ subroutine homogenization_init
                            + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults &      ! crystallite size & crystallite results
                                                         + 1 + constitutive_plasticity_maxSizePostResults &     ! constitutive size & constitutive results
                                                             + constitutive_source_maxSizePostResults)
- allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
+ allocate(materialpoint_results(materialpoint_sizeResults,theMesh%elem%nIPs,theMesh%nElems))
 
  write(6,'(/,a)')   ' <<<+-  homogenization init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@@ -346,7 +343,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
    crystallite_Lp, &
    crystallite_Li0, &
    crystallite_Li, &
-   crystallite_dPdF, &
    crystallite_Tstar0_v, &
    crystallite_Tstar_v, &
    crystallite_partionedF0, &
@@ -614,6 +610,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
 ! crystallite integration
 ! based on crystallite_partionedF0,.._partionedF
 ! incrementing by crystallite_dt
+    
     materialpoint_converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
 
 !--------------------------------------------------------------------------------------------------

src/kinematics_thermal_expansion.f90

@@ -24,10 +24,10 @@ module kinematics_thermal_expansion
  integer(pInt),                       dimension(:),           allocatable, target, public  :: &
    kinematics_thermal_expansion_Noutput                                                                           !< number of outputs per instance of this damage 
 
-! enum, bind(c)                                                                                                   ! ToDo kinematics need state machinery to deal with sizePostResult
+ enum, bind(c)                                                                                                   ! ToDo kinematics need state machinery to deal with sizePostResult
-!   enumerator :: undefined_ID, &                                                                                 ! possible remedy is to decouple having state vars from having output
+   enumerator :: undefined_ID, &                                                                                 ! possible remedy is to decouple having state vars from having output
-!                 thermalexpansionrate_ID                                                                         ! which means to separate user-defined types tState + tOutput...
+                 thermalexpansionrate_ID                                                                         ! which means to separate user-defined types tState + tOutput...
-! end enum
+ end enum
  public :: &
    kinematics_thermal_expansion_init, &
    kinematics_thermal_expansion_initialStrain, &

src/material.f90

@@ -305,9 +305,7 @@ subroutine material_init()
    texture_name
  use mesh, only: &
    mesh_homogenizationAt, &
-   mesh_NipsPerElem, &
+   theMesh
-   mesh_NcpElems, &
-   FE_geomtype
 
  implicit none
  integer(pInt), parameter :: FILEUNIT = 210_pInt
@@ -399,10 +397,10 @@ subroutine material_init()
  call material_populateGrains
 
 ! BEGIN DEPRECATED
- allocate(phaseAt                   (  homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
+ allocate(phaseAt                   (  homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
- allocate(phasememberAt             (  homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
+ allocate(phasememberAt             (  homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
- allocate(mappingHomogenization     (2,                          mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
+ allocate(mappingHomogenization     (2,                          theMesh%elem%nIPs,theMesh%Nelems),source=0_pInt)
- allocate(mappingHomogenizationConst(                            mesh_nIPsPerElem,mesh_NcpElems),source=1_pInt)
+ allocate(mappingHomogenizationConst(                            theMesh%elem%nIPs,theMesh%Nelems),source=1_pInt)
 ! END DEPRECATED
 
  allocate(material_homogenizationAt,source=mesh_homogenizationAt)
@@ -410,9 +408,9 @@ subroutine material_init()
  allocate(CounterHomogenization(size(config_homogenization)),source=0_pInt)
 
 ! BEGIN DEPRECATED
- do e = 1_pInt,mesh_NcpElems
+ do e = 1_pInt,theMesh%Nelems
  myHomog = mesh_homogenizationAt(e)
-   do i = 1_pInt, mesh_NipsPerElem
+   do i = 1_pInt, theMesh%elem%nIPs
      CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1_pInt
      mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),myHomog]
      do g = 1_pInt,homogenization_Ngrains(myHomog)
@@ -553,7 +551,7 @@ subroutine material_parseMicrostructure
    microstructure_name
  use mesh, only: &
    mesh_microstructureAt, &
-   mesh_NcpElems
+   theMesh
 
  implicit none
  character(len=65536), dimension(:), allocatable :: &
@@ -571,7 +569,7 @@ subroutine material_parseMicrostructure
  if(any(mesh_microstructureAt > size(config_microstructure))) &
   call IO_error(155_pInt,ext_msg='More microstructures in geometry than sections in material.config')
 
- forall (e = 1_pInt:mesh_NcpElems) &
+ forall (e = 1_pInt:theMesh%Nelems) &
    microstructure_active(mesh_microstructureAt(e)) = .true.                                         ! current microstructure used in model? Elementwise view, maximum N operations for N elements
 
  do m=1_pInt, size(config_microstructure)
@@ -984,13 +982,10 @@ subroutine material_populateGrains
    math_sampleFiberOri, &
    math_symmetricEulers
  use mesh, only: &
-   mesh_NipsPerElem, &
-   mesh_elemType, &
    mesh_homogenizationAt, &
    mesh_microstructureAt, &
-   mesh_NcpElems, &
+   theMesh, &
-   mesh_ipVolume, &
+   mesh_ipVolume
-   FE_geomtype
  use config, only: &
    config_homogenization, &
    config_microstructure, &
@@ -1026,24 +1021,24 @@ subroutine material_populateGrains
 
  myDebug = debug_level(debug_material)
 
- allocate(material_volume(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),       source=0.0_pReal)
+ allocate(material_volume(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),       source=0.0_pReal)
- allocate(material_phase(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),        source=0_pInt)
+ allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),        source=0_pInt)
- allocate(material_homog(mesh_nIPsPerElem,mesh_NcpElems),                                  source=0_pInt)
+ allocate(material_homog(theMesh%elem%nIPs,theMesh%Nelems),                                  source=0_pInt)
- allocate(material_texture(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),      source=0_pInt)
+ allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),      source=0_pInt)
- allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0.0_pReal)
+ allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal)
 
  allocate(Ngrains(size(config_homogenization),size(config_microstructure)),            source=0_pInt)
  allocate(Nelems (size(config_homogenization),size(config_microstructure)),            source=0_pInt)
 
 ! populating homogenization schemes in each
 !--------------------------------------------------------------------------------------------------
- do e = 1_pInt, mesh_NcpElems
+ do e = 1_pInt, theMesh%Nelems
-   material_homog(1_pInt:mesh_NipsPerElem,e) = mesh_homogenizationAt(e)
+   material_homog(1_pInt:theMesh%elem%nIPs,e) = mesh_homogenizationAt(e)
  enddo
 
 !--------------------------------------------------------------------------------------------------
 ! precounting of elements for each homog/micro pair
- do e = 1_pInt, mesh_NcpElems
+ do e = 1_pInt, theMesh%Nelems
    homog = mesh_homogenizationAt(e)
    micro = mesh_microstructureAt(e)
    Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
@@ -1061,8 +1056,7 @@ subroutine material_populateGrains
 !--------------------------------------------------------------------------------------------------
 ! identify maximum grain count per IP (from element) and find grains per homog/micro pair
  Nelems = 0_pInt                                                                                    ! reuse as counter
- elementLooping: do e = 1_pInt,mesh_NcpElems
+ elementLooping: do e = 1_pInt,theMesh%Nelems
-   t     = mesh_elemType
    homog = mesh_homogenizationAt(e)
    micro = mesh_microstructureAt(e)
    if (homog < 1_pInt .or. homog > size(config_homogenization)) &                                      ! out of bounds
@@ -1072,7 +1066,7 @@ subroutine material_populateGrains
    if (microstructure_elemhomo(micro)) then                                                         ! how many grains are needed at this element?
      dGrains = homogenization_Ngrains(homog)                                                        ! only one set of Ngrains (other IPs are plain copies)
    else
-     dGrains = homogenization_Ngrains(homog) * mesh_NipsPerElem                                     ! each IP has Ngrains
+     dGrains = homogenization_Ngrains(homog) * theMesh%elem%nIPs                                     ! each IP has Ngrains
    endif
    Ngrains(homog,micro) = Ngrains(homog,micro) + dGrains                                            ! total grain count
    Nelems(homog,micro)  = Nelems(homog,micro) + 1_pInt                                              ! total element count
@@ -1106,16 +1100,15 @@ subroutine material_populateGrains
 
        do hme = 1_pInt, Nelems(homog,micro)
          e = elemsOfHomogMicro(homog,micro)%p(hme)                                                  ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
-         t = mesh_elemType
          if (microstructure_elemhomo(micro)) then                                                   ! homogeneous distribution of grains over each element's IPs
-           volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:mesh_NipsPerElem,e))/&
+           volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:theMesh%elem%nIPs,e))/&
                                                                          real(dGrains,pReal)        ! each grain combines size of all IPs in that element
            grain = grain + dGrains                                                                  ! wind forward by Ngrains@IP
          else
-           forall (i = 1_pInt:mesh_NipsPerElem) &                                                   ! loop over IPs
+           forall (i = 1_pInt:theMesh%elem%nIPs) &                                                   ! loop over IPs
              volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
                mesh_ipVolume(i,e)/real(dGrains,pReal)                                               ! assign IPvolume/Ngrains@IP to all grains of IP
-           grain = grain + mesh_NipsPerElem * dGrains                                               ! wind forward by Nips*Ngrains@IP
+           grain = grain + theMesh%elem%nIPs * dGrains                                               ! wind forward by Nips*Ngrains@IP
          endif
        enddo
 
@@ -1261,11 +1254,10 @@ subroutine material_populateGrains
 
        do hme = 1_pInt, Nelems(homog,micro)
          e = elemsOfHomogMicro(homog,micro)%p(hme)                                                  ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
-         t = mesh_elemType
          if (microstructure_elemhomo(micro)) then                                                   ! homogeneous distribution of grains over each element's IPs
            m = 1_pInt                                                                               ! process only first IP
          else
-           m = mesh_NipsPerElem
+           m = theMesh%elem%nIPs
          endif
 
          do i = 1_pInt, m                                                                           ! loop over necessary IPs
@@ -1303,7 +1295,7 @@ subroutine material_populateGrains
 
          enddo
 
-         do i = i, mesh_NipsPerElem                                                                 ! loop over IPs to (possibly) distribute copies from first IP
+         do i = i, theMesh%elem%nIPs                                                                 ! loop over IPs to (possibly) distribute copies from first IP
            material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e)
            material_phase  (1_pInt:dGrains,i,e) = material_phase  (1_pInt:dGrains,1,e)
            material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e)

src/mesh_FEM.f90

@@ -12,7 +12,7 @@ module mesh
 #include <petsc/finclude/petscis.h>
 #include <petsc/finclude/petscdmda.h>
  use prec, only: pReal, pInt
-
+ use mesh_base
 use PETScdmplex
 use PETScdmda
 use PETScis
@@ -27,7 +27,6 @@ use PETScis
    mesh_NcpElems, &                                                                                 !< total number of CP elements in mesh
    mesh_NcpElemsGlobal, &
    mesh_Nnodes, &                                                                                   !< total number of nodes in mesh
-   mesh_NipsPerElem, &                                                                              !< number of IPs in per element
    mesh_maxNipNeighbors
 !!!! BEGIN DEPRECATED !!!!!
  integer(pInt), public, protected :: &
@@ -80,6 +79,16 @@ use PETScis
  integer(pInt), dimension(1_pInt), parameter, public :: FE_NipNeighbors = &                  !< number of ip neighbors / cell faces in a specific cell type
  int([6],pInt)
  
+  type, public, extends(tMesh) :: tMesh_FEM
+
+   
+   contains
+   procedure, pass(self) :: tMesh_FEM_init
+   generic, public :: init => tMesh_FEM_init
+ end type tMesh_FEM
+ 
+ type(tMesh_FEM), public, protected :: theMesh
+ 
 
  public :: &
    mesh_init, &
@@ -89,6 +98,25 @@ use PETScis
 
 contains
 
+subroutine tMesh_FEM_init(self,dimen,order,nodes)
+ 
+ implicit none
+  integer, intent(in) :: dimen
+ integer(pInt), intent(in) :: order
+  real(pReal), intent(in), dimension(:,:) :: nodes
+ class(tMesh_FEM) :: self
+ 
+ if (dimen == 2_pInt) then
+   if (order == 1_pInt) call  self%tMesh%init('mesh',1_pInt,nodes)
+   if (order == 2_pInt) call  self%tMesh%init('mesh',2_pInt,nodes)
+ elseif(dimen == 3_pInt) then
+   if (order == 1_pInt) call self%tMesh%init('mesh',6_pInt,nodes)
+   if (order == 2_pInt) call self%tMesh%init('mesh',8_pInt,nodes)
+ endif
+
+ end subroutine tMesh_FEM_init
+
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief initializes the mesh by calling all necessary private routines the mesh module
@@ -213,6 +241,8 @@ subroutine mesh_init()
 
  FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
  mesh_maxNips = FE_Nips(1_pInt)
+ 
+ write(6,*) 'mesh_maxNips',mesh_maxNips
  call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
  call mesh_FEM_build_ipVolumes(dimPlex)
  
@@ -238,11 +268,13 @@ subroutine mesh_init()
 !!!! COMPATIBILITY HACK !!!!
 ! for a homogeneous mesh, all elements have the same number of IPs and and cell nodes.
 ! hence, xxPerElem instead of maxXX
- mesh_NipsPerElem      = mesh_maxNips
 ! better name
  mesh_homogenizationAt = mesh_element(3,:)
  mesh_microstructureAt = mesh_element(4,:)
 !!!!!!!!!!!!!!!!!!!!!!!!
+ 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
 

src/mesh_base.f90

@@ -0,0 +1,85 @@
+
+!--------------------------------------------------------------------------------------------------
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @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 Sets up the mesh for the solvers MSC.Marc,FEM, Abaqus and the spectral solver
+!--------------------------------------------------------------------------------------------------
+module mesh_base
+
+ use, intrinsic :: iso_c_binding
+ use prec, only: &
+   pStringLen, &
+   pReal, &
+   pInt
+ use element, only: &
+   tElement
+
+  implicit none
+
+!---------------------------------------------------------------------------------------------------
+!> Properties of a the whole mesh (consisting of one type of elements)
+!---------------------------------------------------------------------------------------------------
+ type, public :: tMesh
+   type(tElement) :: &
+     elem
+   real(pReal), dimension(:,:), allocatable, public :: &
+     ipVolume, &                                                                                 !< volume associated with each IP (initially!)
+     node0, &                                                                                    !< node x,y,z coordinates (initially)
+     node                                                                                        !< node x,y,z coordinates (deformed)
+   integer(pInt), dimension(:,:), allocatable, public :: &    
+     cellnodeParent                                                                               !< cellnode's parent element ID, cellnode's intra-element ID 
+   character(pStringLen) :: type = "n/a"
+   integer(pInt)         :: &
+     Nnodes, &                                                                                   !< total number of nodes in mesh
+     Nelems = -1_pInt, &
+     elemType, &
+     Ncells, &
+     nIPneighbors, &
+     NcellNodes, &
+     maxElemsPerNode
+   integer(pInt), dimension(:), allocatable, public :: &
+     homogenizationAt, &
+     microstructureAt
+   integer(pInt), dimension(:,:), allocatable, public :: &
+     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
+
+contains
+subroutine tMesh_base_init(self,meshType,elemType,nodes)
+ 
+ implicit none
+ class(tMesh) :: self
+ character(len=*), intent(in) :: meshType
+ integer(pInt), intent(in) :: elemType
+ real(pReal), dimension(:,:), intent(in) :: nodes
+ 
+ write(6,'(/,a)')   ' <<<+-  mesh_base_init  -+>>>'
+ 
+ write(6,*)' mesh type ',meshType
+ write(6,*)' # node    ',size(nodes,2)
+
+ self%type = meshType
+ call self%elem%init(elemType)
+ self%node0 = nodes
+ self%nNodes = size(nodes,2)
+
+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

src/numerics.f90

@@ -177,13 +177,8 @@ subroutine numerics_init
 #include <petsc/finclude/petscsys.h>
    use petscsys
 #endif
-#if !defined(Marc4DAMASK)
+!$ use OMP_LIB, only: omp_set_num_threads
-!$ use OMP_LIB, only: omp_set_num_threads                                                           ! Standard conforming module
  implicit none
-#else  
- implicit none
-!$ include "omp_lib.h"                                                                              ! MSC.Marc includes this file on !its own, avoid conflict with the OMP_LIB module
-#endif
  integer(pInt), parameter ::                 FILEUNIT = 300_pInt
 !$ integer ::                                gotDAMASK_NUM_THREADS = 1
  integer :: i, ierr                                                                                 ! no pInt

src/plastic_nonlocal.f90

@@ -309,9 +309,7 @@ use IO,       only: IO_read, &
 use debug,    only: debug_level, &
                     debug_constitutive, &
                     debug_levelBasic
-use mesh,     only: mesh_NcpElems, &
+use mesh,     only: theMesh
-                    mesh_maxNips, &
-                    mesh_maxNipNeighbors
 use material, only: phase_plasticity, &
                     homogenization_maxNgrains, &
                     phase_plasticityInstance, &
@@ -430,6 +428,7 @@ allocate(minDipoleHeightPerSlipFamily(lattice_maxNslipFamily,2,maxNinstances), s
 allocate(peierlsStressPerSlipFamily(lattice_maxNslipFamily,2,maxNinstances),   source=0.0_pReal)
 allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances),                  source=0.0_pReal)
 
+
  rewind(fileUnit)
  phase = 0_pInt
  do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
@@ -697,23 +696,23 @@ allocate(forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstances),
 allocate(forestProjectionScrew(maxTotalNslip,maxTotalNslip,maxNinstances),            source=0.0_pReal)
 allocate(interactionMatrixSlipSlip(maxTotalNslip,maxTotalNslip,maxNinstances),        source=0.0_pReal)
 allocate(lattice2slip(1:3, 1:3, maxTotalNslip,maxNinstances),                         source=0.0_pReal)
-allocate(sourceProbability(maxTotalNslip,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(sourceProbability(maxTotalNslip,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=2.0_pReal)
 
-allocate(rhoDotFluxOutput(maxTotalNslip,8,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotFluxOutput(maxTotalNslip,8,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
-allocate(rhoDotMultiplicationOutput(maxTotalNslip,2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotMultiplicationOutput(maxTotalNslip,2,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
-allocate(rhoDotSingle2DipoleGlideOutput(maxTotalNslip,2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotSingle2DipoleGlideOutput(maxTotalNslip,2,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
-allocate(rhoDotAthermalAnnihilationOutput(maxTotalNslip,2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotAthermalAnnihilationOutput(maxTotalNslip,2,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
-allocate(rhoDotThermalAnnihilationOutput(maxTotalNslip,2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotThermalAnnihilationOutput(maxTotalNslip,2,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
-allocate(rhoDotEdgeJogsOutput(maxTotalNslip,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), &
+allocate(rhoDotEdgeJogsOutput(maxTotalNslip,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
 
-allocate(compatibility(2,maxTotalNslip,maxTotalNslip,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems), &
+allocate(compatibility(2,maxTotalNslip,maxTotalNslip,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=0.0_pReal)
 allocate(peierlsStress(maxTotalNslip,2,maxNinstances),                                source=0.0_pReal)
 allocate(colinearSystem(maxTotalNslip,maxNinstances),                                 source=0_pInt)
@@ -798,7 +797,6 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
      
                 
 
-
      plasticState(phase)%nonlocal = .true.
      call material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState,sizeDeltaState, &
                                         totalNslip(instance),0_pInt,0_pInt)
@@ -999,10 +997,8 @@ use IO,       only: IO_error
 use lattice,  only: lattice_maxNslipFamily
 use math,     only: math_sampleGaussVar
 use mesh,     only: mesh_ipVolume, &
-                    mesh_NcpElems, &
+                    theMesh, &
-                    mesh_element, &
+                    mesh_element
-                    FE_Nips, &
-                    FE_geomtype
 use material, only: material_phase, &
                     phase_plasticityInstance, &
                     plasticState, &
@@ -1041,8 +1037,8 @@ do instance = 1_pInt,maxNinstances
 
     minimumIpVolume = huge(1.0_pReal)
     totalVolume = 0.0_pReal
-    do e = 1_pInt,mesh_NcpElems
+    do e = 1_pInt,theMesh%nElems
-      do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
+      do i = 1_pInt,theMesh%elem%nIPs
         if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) &
             .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then
           totalVolume = totalVolume + mesh_ipVolume(i,e)
@@ -1057,8 +1053,8 @@ do instance = 1_pInt,maxNinstances
     meanDensity = 0.0_pReal
     do while(meanDensity < rhoSglRandom(instance))
       call random_number(rnd)
-      e = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt)
+      e = nint(rnd(1)*real(theMesh%nElems,pReal)+0.5_pReal,pInt)
-      i = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,e))),pReal)+0.5_pReal,pInt)
+      i = nint(rnd(2)*real(theMesh%elem%nIPs,pReal)+0.5_pReal,pInt)
       if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) &
           .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then
         s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt)
@@ -1071,8 +1067,8 @@ do instance = 1_pInt,maxNinstances
     enddo
   ! homogeneous distribution of density with some noise
   else
-    do e = 1_pInt,mesh_NcpElems
+    do e = 1_pInt,theMesh%nElems
-      do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
+      do i = 1_pInt,theMesh%elem%nIPs
         if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) &
             .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then
           do f = 1_pInt,lattice_maxNslipFamily
@@ -1154,16 +1150,13 @@ use debug, only: &
   debug_i, &
   debug_e
 use mesh, only: &
+  theMesh, &
   mesh_element, &
   mesh_ipNeighborhood, &
   mesh_ipCoordinates, &
   mesh_ipVolume, &
   mesh_ipAreaNormal, &
-  mesh_ipArea, &
+  mesh_ipArea
-  FE_NipNeighbors, &
-  mesh_maxNipNeighbors, &
-  FE_geomtype, &
-  FE_celltype
 use material, only: &
   material_phase, &
   phase_localPlasticity, &
@@ -1223,7 +1216,7 @@ real(pReal), dimension(3,3) ::  invFe, &                      ! inverse of elast
                                 invFp, &                      ! inverse of plastic deformation gradient
                                 connections, &
                                 invConnections
-real(pReal), dimension(3,mesh_maxNipNeighbors) :: &
+real(pReal), dimension(3,theMesh%elem%nIPneighbors) :: &
                                 connection_latticeConf
 real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: &
                                 rhoExcess
@@ -1234,7 +1227,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt
 real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))), &
                        totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el)))) :: &
                                 myInteractionMatrix           ! corrected slip interaction matrix
-real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: &
+real(pReal), dimension(2,maxval(totalNslip),theMesh%elem%nIPneighbors) :: &
                                 neighbor_rhoExcess, &      ! excess density at neighboring material point
                                 neighbor_rhoTotal          ! total density at neighboring material point
 real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: &
@@ -1309,7 +1302,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
   
   nRealNeighbors = 0_pInt
   neighbor_rhoTotal = 0.0_pReal
-  do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))
+  do n = 1_pInt,theMesh%elem%nIPneighbors
     neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
     neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
     np = phaseAt(1,neighbor_ip,neighbor_el)
@@ -1995,16 +1988,12 @@ use math,     only: math_mul6x6, &
                     math_det33, &
                     math_transpose33, &  
                     pi
-use mesh,     only: mesh_NcpElems, &
+use mesh,     only: theMesh, &
-                    mesh_maxNips, &
                     mesh_element, &
                     mesh_ipNeighborhood, &
                     mesh_ipVolume, &
                     mesh_ipArea, &
-                    mesh_ipAreaNormal, &
+                    mesh_ipAreaNormal
-                    FE_NipNeighbors, &
-                    FE_geomtype, &
-                    FE_celltype
 use material, only: homogenization_maxNgrains, &
                     material_phase, &
                     phase_plasticityInstance, &
@@ -2030,9 +2019,9 @@ integer(pInt), intent(in) ::                ip, &
 real(pReal), intent(in) ::                  Temperature, &                                          !< temperature
                                             timestep                                                !< substepped crystallite time increment
 real(pReal), dimension(6), intent(in) ::    Tstar_v                                                 !< current 2nd Piola-Kirchhoff stress in Mandel notation
-real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
+real(pReal), dimension(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
                                             subfrac                                                 !< fraction of timestep at the beginning of the substepped crystallite time increment 
-real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
+real(pReal), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
                                             Fe, &                                                   !< elastic deformation gradient
                                             Fp                                                      !< plastic deformation gradient
 
@@ -2311,8 +2300,8 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
   my_Fe = Fe(1:3,1:3,1_pInt,ip,el)
   my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1_pInt,ip,el))
   
-  do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))                       ! loop through my neighbors
+  do n = 1_pInt,theMesh%elem%nIPneighbors
-!       write(6,*) 'c'
+
     neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
     neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
     neighbor_n  = mesh_ipNeighborhood(3,n,ip,el)
@@ -2611,11 +2600,7 @@ use material, only:   material_phase, &
                       homogenization_maxNgrains
 use mesh, only:       mesh_element, &
                       mesh_ipNeighborhood, &
-                      mesh_maxNips, &
+                      theMesh
-                      mesh_NcpElems, &
-                      FE_NipNeighbors, &
-                      FE_geomtype, &
-                      FE_celltype
 use lattice, only:    lattice_sn, &
                       lattice_sd, &
                       lattice_qDisorientation
@@ -2625,7 +2610,7 @@ implicit none
 !* input variables
 integer(pInt), intent(in) ::                    i, &                                                ! ip index
                                                 e                                                   ! element index
-real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
+real(pReal), dimension(4,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
                                                 orientation                                         ! crystal orientation in quaternions
                                             
 !* local variables
@@ -2644,7 +2629,7 @@ integer(pInt)                                   Nneighbors, &
 real(pReal), dimension(4) ::                    absoluteMisorientation                              ! absolute misorientation (without symmetry) between me and my neighbor
 real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
                          totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
-                         FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))) :: &  
+                         theMesh%elem%nIPneighbors) :: &  
                                                 my_compatibility                                    ! my_compatibility for current element and ip
 real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &  
                                                 slipNormal, &
@@ -2656,7 +2641,7 @@ logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e))))
                                                 belowThreshold
 
 
-Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
+Nneighbors = theMesh%elem%nIPneighbors
 ph = material_phase(1,i,e)
 textureID = material_texture(1,i,e)
 instance = phase_plasticityInstance(ph)
@@ -2771,8 +2756,7 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
    math_mul33x33, &
    pi
  use mesh, only: &
-   mesh_NcpElems, &
+   theMesh
-   mesh_maxNips
  use material, only: &
    homogenization_maxNgrains, &
    material_phase, &
@@ -2790,7 +2774,7 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
  implicit none
  real(pReal),   dimension(6),  intent(in) :: &
    Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal),   dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
+ real(pReal),   dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
    Fe                                                                                               !< elastic deformation gradient
  integer(pInt),                intent(in) :: &
    ip, &                                                                                            !< integration point
@@ -2900,8 +2884,8 @@ forall (s = 1_pInt:ns) &
                                       lattice_sn(1:3,slipSystemLattice(s,instance),ph))
 
 
- outputsLoop: do o = 1_pInt,size(param(instance)%outputID)
+outputsLoop: do o = 1_pInt,plastic_nonlocal_Noutput(instance)
-   select case(param(instance)%outputID(o))
+  select case(plastic_nonlocal_outputID(o,instance))
       
     case (rho_sgl_edge_pos_mobile_ID)
       plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1)