Merge remote-tracking branch 'origin/development' into grid-mesh-cleanup

.gitlab-ci.yml

@@ -506,7 +506,7 @@ Processing:
     - rm abq_addUserOutput.py marc_addUserOutput.py
     - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
     - cd $DAMASKROOT/processing/post
-    - rm marc_to_vtk.py vtk2ang.py DAD*.py
+    - rm vtk2ang.py DAD*.py
     - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
   except:
     - master

VERSION

@@ -1 +1 @@
-v2.0.3-252-g0e77be25
+v2.0.3-261-g99878952

processing/post/marc_to_vtk.py

@@ -1,199 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os,sys,re
-import argparse
-import damask
-import vtk, numpy as np
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName, damask.version])
-
-parser = argparse.ArgumentParser(description='Convert from Marc input file format (.dat) to VTK format (.vtu)', version = scriptID)
-parser.add_argument('filename', type=str, help='file to convert')
-parser.add_argument('-t', '--table',  type=str, help='ASCIItable file containing nodal data to subdivide and interpolate')
-
-args = parser.parse_args()
-
-with open(args.filename, 'r') as marcfile:
-  marctext = marcfile.read();
-  
-# Load table (if any)
-if args.table is not None:
-  try:
-    table = damask.ASCIItable(
-      name=args.table,
-      outname='subdivided_{}'.format(args.table),
-      buffered=True
-    )
-    
-    table.head_read()
-    table.data_readArray()
-    
-    # Python list is faster for appending
-    nodal_data = list(table.data)
-  except: args.table = None
-
-# Extract connectivity chunk from file...
-connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
-connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL))
-connectivity_header = connectivity_lines[0]
-connectivity_lines = connectivity_lines[1:]
-
-# Construct element map
-elements = dict(map(lambda line: 
-  (
-    int(line[0:10]), # index
-    { 
-      'type':  int(line[10:20]), 
-      'verts': list(map(int, re.split(r' +', line[20:].strip())))
-    }
-  ), connectivity_lines))
-
-# Extract coordinate chunk from file
-coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
-coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
-coordinates_header = coordinates_lines[0]
-coordinates_lines = coordinates_lines[1:]
-
-# marc input file does not use "e" in scientific notation, this adds it and converts
-fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string))
-# Construct coordinate map
-coordinates = dict(map(lambda line:
-  (
-    int(line[0:10]),
-    np.array([
-      fl_format(line[10:30]), 
-      fl_format(line[30:50]), 
-      fl_format(line[50:70])
-    ])
-  ), coordinates_lines))
-
-# Subdivide volumes
-grid = vtk.vtkUnstructuredGrid()
-vertex_count = len(coordinates)
-edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here
-ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional
-
-def subdivide_edge(vert1, vert2):
-  edge = ordered_pair(vert1, vert2)
-  
-  if edge in edge_to_vert:
-    return edge_to_vert[edge]
-  
-  # Vertex does not exist, create it
-  newvert = len(coordinates) + 1
-  coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average
-  edge_to_vert[edge] = newvert;
-  
-  # Interpolate nodal data
-  if args.table is not None:
-    nodal_data.append(0.5 * (nodal_data[vert1 - 1] + nodal_data[vert2 - 1]))
-  return newvert;
-
-for el_id in range(1, len(elements) + 1): # Marc starts counting at 1
-  el = elements[el_id]
-  if el['type'] == 7:
-    # Hexahedron, subdivided
-    
-    # There may be a better way to iterate over these, but this is consistent
-    # with the ordering scheme provided at https://damask.mpie.de/pub/Documentation/ElementType
-    
-    subverts = np.zeros((3,3,3), dtype=int)
-    # Get corners
-    subverts[0, 0, 0] = el['verts'][0]
-    subverts[2, 0, 0] = el['verts'][1]
-    subverts[2, 2, 0] = el['verts'][2]
-    subverts[0, 2, 0] = el['verts'][3]
-    subverts[0, 0, 2] = el['verts'][4]
-    subverts[2, 0, 2] = el['verts'][5]
-    subverts[2, 2, 2] = el['verts'][6]
-    subverts[0, 2, 2] = el['verts'][7]
-    
-    # lower edges
-    subverts[1, 0, 0] = subdivide_edge(subverts[0, 0, 0], subverts[2, 0, 0])
-    subverts[2, 1, 0] = subdivide_edge(subverts[2, 0, 0], subverts[2, 2, 0])
-    subverts[1, 2, 0] = subdivide_edge(subverts[2, 2, 0], subverts[0, 2, 0])
-    subverts[0, 1, 0] = subdivide_edge(subverts[0, 2, 0], subverts[0, 0, 0])
-    
-    # middle edges
-    subverts[0, 0, 1] = subdivide_edge(subverts[0, 0, 0], subverts[0, 0, 2])
-    subverts[2, 0, 1] = subdivide_edge(subverts[2, 0, 0], subverts[2, 0, 2])
-    subverts[2, 2, 1] = subdivide_edge(subverts[2, 2, 0], subverts[2, 2, 2])
-    subverts[0, 2, 1] = subdivide_edge(subverts[0, 2, 0], subverts[0, 2, 2])
-    
-    # top edges
-    subverts[1, 0, 2] = subdivide_edge(subverts[0, 0, 2], subverts[2, 0, 2])
-    subverts[2, 1, 2] = subdivide_edge(subverts[2, 0, 2], subverts[2, 2, 2])
-    subverts[1, 2, 2] = subdivide_edge(subverts[2, 2, 2], subverts[0, 2, 2])
-    subverts[0, 1, 2] = subdivide_edge(subverts[0, 2, 2], subverts[0, 0, 2])
-    
-    # then faces...  The edge_to_vert addition is due to there being two ways
-    # to calculate a face vertex, depending on which opposite vertices are used to subdivide.
-    # This way, we avoid creating duplicate vertices.
-    subverts[1, 1, 0] = subdivide_edge(subverts[1, 0, 0], subverts[1, 2, 0])
-    edge_to_vert[ordered_pair(subverts[0, 1, 0], subverts[2, 1, 0])] = subverts[1, 1, 0]
-    
-    subverts[1, 0, 1] = subdivide_edge(subverts[1, 0, 0], subverts[1, 0, 2])
-    edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[2, 0, 1])] = subverts[1, 0, 1]
-    
-    subverts[2, 1, 1] = subdivide_edge(subverts[2, 1, 0], subverts[2, 1, 2])
-    edge_to_vert[ordered_pair(subverts[2, 0, 1], subverts[2, 2, 1])] = subverts[2, 1, 1]
-    
-    subverts[1, 2, 1] = subdivide_edge(subverts[1, 2, 0], subverts[1, 2, 2])
-    edge_to_vert[ordered_pair(subverts[0, 2, 1], subverts[2, 2, 1])] = subverts[1, 2, 1]
-    
-    subverts[0, 1, 1] = subdivide_edge(subverts[0, 1, 0], subverts[0, 1, 2])
-    edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[0, 2, 1])] = subverts[0, 1, 1]
-    
-    subverts[1, 1, 2] = subdivide_edge(subverts[1, 0, 2], subverts[1, 2, 2])
-    edge_to_vert[ordered_pair(subverts[0, 1, 2], subverts[2, 1, 2])] = subverts[1, 1, 2]
-    
-    # and finally the center.  There are three ways to calculate, but elements should
-    # not intersect, so the edge_to_vert part isn't needed here.
-    subverts[1, 1, 1] = subdivide_edge(subverts[1, 1, 0], subverts[1, 1, 2])
-    
-    
-    # Now make the hexahedron subelements
-    # order in which vtk expects vertices for a hexahedron
-    order = np.array([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)])
-    for z in range(2):
-      for y in range(2):
-        for x in range(2):
-          hex_ = vtk.vtkHexahedron()
-          for vert_id in range(8):
-            coord = order[vert_id] + (x, y, z)
-            # minus one, since vtk starts at zero but marc starts at one
-            hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) 
-          grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
-    
-  else:
-    damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type']))
-
-# Load all points
-points = vtk.vtkPoints()
-for point in range(1, len(coordinates) + 1): # marc indices start at 1
-  points.InsertNextPoint(coordinates[point].tolist())
-
-grid.SetPoints(points)
-
-# grid now contains the elements from the given marc file
-writer = vtk.vtkXMLUnstructuredGridWriter()
-writer.SetFileName(re.sub(r'\..+', ".vtu", args.filename))  # *.vtk extension does not work in paraview
-
-if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
-else:                          writer.SetInputData(grid)
-writer.Write()
-
-if args.table is not None:  
-  table.info_append([
-    scriptID + ' ' + ' '.join(sys.argv[1:]),
-  ])
-  table.head_write()
-  table.output_flush()
-  
-  table.data = np.array(nodal_data)
-  
-  table.data_writeArray()
-
-table.close()

src/CPFEM2.f90

@@ -299,6 +299,8 @@ subroutine CPFEM_results(inc,time)
    pInt
  use results
  use HDF5_utilities
+  use homogenization, only: &
+   homogenization_results
  use constitutive, only: &
    constitutive_results
  use crystallite, only: &
@@ -312,6 +314,7 @@ subroutine CPFEM_results(inc,time)
  call results_addIncrement(inc,time)
  call constitutive_results
  call crystallite_results
+ call homogenization_results
  call results_removeLink('current') ! ToDo: put this into closeJobFile
  call results_closeJobFile
 

src/DAMASK_interface.f90

@@ -15,6 +15,9 @@
 #define PETSC_MINOR_MIN 10
 #define PETSC_MINOR_MAX 11
 module DAMASK_interface
+  use prec
+  use system_routines
+  
   implicit none
   private
   logical, public, protected :: &
@@ -39,6 +42,7 @@ module DAMASK_interface
     getLoadCaseFile, &
     rectifyPath, &
     makeRelativePath
+
 contains
 
 !--------------------------------------------------------------------------------------------------
@@ -46,18 +50,8 @@ contains
 !! information on computation to screen
 !--------------------------------------------------------------------------------------------------
 subroutine DAMASK_interface_init
-  use, intrinsic :: &
+  use, intrinsic :: iso_fortran_env
-    iso_fortran_env
-  use, intrinsic :: &
-    iso_c_binding
   use PETScSys
-  use prec, only: &
-    pReal
-  use system_routines, only: &
-    signalusr1_C, &
-    signalusr2_C, &
-    getHostName, &
-    getCWD
 
 #include <petsc/finclude/petscsys.h>
 #if defined(__GFORTRAN__) &&  __GNUC__<GCC_MIN
@@ -96,7 +90,6 @@ subroutine DAMASK_interface_init
 ===================================================================================================
 #endif
 
-  implicit none
   character(len=1024) :: &
     commandLine, &                                                                                  !< command line call as string
     arg, &                                                                                          !< individual argument
@@ -304,11 +297,7 @@ end subroutine DAMASK_interface_init
 !!        possibly converting relative arguments to absolut path
 !--------------------------------------------------------------------------------------------------
 subroutine setWorkingDirectory(workingDirectoryArg)
-  use system_routines, only: &
-    getCWD, &
-    setCWD
 
-  implicit none
   character(len=*),  intent(in) :: workingDirectoryArg                                              !< working directory argument
   character(len=1024)           :: workingDirectory                                                 !< working directory argument
   logical                       :: error
@@ -336,22 +325,17 @@ end subroutine setWorkingDirectory
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getSolverJobName()
 
-  implicit none
   integer :: posExt,posSep
-  character(len=1024) :: tempString
 
+  posExt = scan(geometryFile,'.',back=.true.)
+  posSep = scan(geometryFile,'/',back=.true.)
 
-  tempString = geometryFile
+  getSolverJobName = geometryFile(posSep+1:posExt-1)
-  posExt = scan(tempString,'.',back=.true.)
-  posSep = scan(tempString,'/',back=.true.)
 
-  getSolverJobName = tempString(posSep+1:posExt-1)
+  posExt = scan(loadCaseFile,'.',back=.true.)
+  posSep = scan(loadCaseFile,'/',back=.true.)
 
-  tempString = loadCaseFile
+  getSolverJobName = trim(getSolverJobName)//'_'//loadCaseFile(posSep+1:posExt-1)
-  posExt = scan(tempString,'.',back=.true.)
-  posSep = scan(tempString,'/',back=.true.)
-
-  getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1)
 
 end function getSolverJobName
 
@@ -360,10 +344,7 @@ end function getSolverJobName
 !> @brief basename of geometry file with extension from command line arguments
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getGeometryFile(geometryParameter)
-  use system_routines, only: &
-    getCWD
 
-  implicit none
   character(len=1024), intent(in) :: geometryParameter
   logical                         :: file_exists
   external                        :: quit
@@ -385,10 +366,7 @@ end function getGeometryFile
 !> @brief relative path of loadcase from command line arguments
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getLoadCaseFile(loadCaseParameter)
-  use system_routines, only: &
-    getCWD
 
-  implicit none
   character(len=1024), intent(in) :: loadCaseParameter
   logical                         :: file_exists
   external                        :: quit
@@ -412,7 +390,6 @@ end function getLoadCaseFile
 !--------------------------------------------------------------------------------------------------
 function rectifyPath(path)
 
-  implicit none
   character(len=*)    :: path
   character(len=1024) :: rectifyPath
   integer :: i,j,k,l
@@ -457,7 +434,6 @@ end function rectifyPath
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function makeRelativePath(a,b)
 
-  implicit none
   character (len=*), intent(in) :: a,b
   character (len=1024)          :: a_cleaned,b_cleaned
   integer :: i,posLastCommonSlash,remainingSlashes
@@ -484,9 +460,7 @@ end function makeRelativePath
 !> @brief sets global variable SIGTERM to .true.
 !--------------------------------------------------------------------------------------------------
 subroutine catchSIGTERM(signal) bind(C)
-  use :: iso_c_binding
 
-  implicit none
   integer(C_INT), value :: signal
   SIGTERM = .true.
 
@@ -500,7 +474,6 @@ end subroutine catchSIGTERM
 !--------------------------------------------------------------------------------------------------
 subroutine setSIGTERM(state)
 
-  implicit none
   logical, intent(in) :: state
   SIGTERM = state
 
@@ -511,9 +484,7 @@ end subroutine setSIGTERM
 !> @brief sets global variable SIGUSR1 to .true.
 !--------------------------------------------------------------------------------------------------
 subroutine catchSIGUSR1(signal) bind(C)
-  use :: iso_c_binding
 
-  implicit none
   integer(C_INT), value :: signal
   SIGUSR1 = .true.
 
@@ -527,7 +498,6 @@ end subroutine catchSIGUSR1
 !--------------------------------------------------------------------------------------------------
 subroutine setSIGUSR1(state)
 
-  implicit none
   logical, intent(in) :: state
   SIGUSR1 = state
 
@@ -538,9 +508,7 @@ end subroutine setSIGUSR1
 !> @brief sets global variable SIGUSR2 to .true. if program receives SIGUSR2
 !--------------------------------------------------------------------------------------------------
 subroutine catchSIGUSR2(signal) bind(C)
-  use :: iso_c_binding
 
-  implicit none
   integer(C_INT), value :: signal
   SIGUSR2 = .true.
 
@@ -554,7 +522,6 @@ end subroutine catchSIGUSR2
 !--------------------------------------------------------------------------------------------------
 subroutine setSIGUSR2(state)
 
-  implicit none
   logical, intent(in) :: state
   SIGUSR2 = state
 

src/Lambert.f90

@@ -38,11 +38,7 @@
 !> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
 !--------------------------------------------------------------------------
 module Lambert
-  use prec, only: &
+  use math
-    pReal
-  use math, only: &
-    PI
-  use future
 
   implicit none
   private
@@ -73,11 +69,7 @@ contains
 !> @brief map from 3D cubic grid to 3D ball 
 !--------------------------------------------------------------------------
 function LambertCubeToBall(cube) result(ball)
-  use, intrinsic :: IEEE_ARITHMETIC
-  use prec, only: &
-    dEq0
  
-  implicit none
   real(pReal), intent(in), dimension(3) :: cube
   real(pReal),             dimension(3) :: ball, LamXYZ, XYZ
   real(pReal),             dimension(2) :: T
@@ -133,15 +125,7 @@ end function LambertCubeToBall
 !> @brief map from 3D ball to 3D cubic grid  
 !--------------------------------------------------------------------------
 pure function LambertBallToCube(xyz) result(cube)
-  use, intrinsic :: IEEE_ARITHMETIC, only:&
-    IEEE_positive_inf, &
-    IEEE_value
-  use prec, only: &
-    dEq0
-   use math, only: &
-    math_clip
  
-  implicit none
   real(pReal), intent(in), dimension(3) :: xyz
   real(pReal),             dimension(3) :: cube, xyz1, xyz3
   real(pReal),             dimension(2) :: Tinv, xyz2
@@ -196,7 +180,6 @@ end function LambertBallToCube
 !--------------------------------------------------------------------------
 pure function GetPyramidOrder(xyz)
 
-  implicit none
   real(pReal),intent(in),dimension(3) :: xyz
   integer,               dimension(3) :: GetPyramidOrder
  

src/constitutive.f90

@@ -4,8 +4,7 @@
 !> @brief elasticity, plasticity, internal microstructure state
 !--------------------------------------------------------------------------------------------------
 module constitutive
- use prec, only: &
+ use math
-   pInt
 
  implicit none
  private
@@ -38,8 +37,6 @@ contains
 !> @brief allocates arrays pointing to array of the various constitutive modules
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_init
- use prec, only: &
-   pReal
  use debug, only: &
    debug_constitutive, &
    debug_levelBasic
@@ -109,7 +106,6 @@ subroutine constitutive_init
  use kinematics_slipplane_opening
  use kinematics_thermal_expansion
 
- implicit none
  integer, parameter :: FILEUNIT = 204
  integer :: &
    o, &                                                                                             !< counter in output loop
@@ -285,8 +281,6 @@ end subroutine constitutive_init
 !> ToDo: homogenizedC66 would be more consistent
 !--------------------------------------------------------------------------------------------------
 function constitutive_homogenizedC(ipc,ip,el)
- use prec, only: &
-   pReal
  use material, only: &
    phase_plasticity, &
    material_phase, &
@@ -297,7 +291,6 @@ function constitutive_homogenizedC(ipc,ip,el)
  use lattice, only: &
    lattice_C66
 
- implicit none
  real(pReal), dimension(6,6) :: constitutive_homogenizedC
  integer, intent(in) :: &
    ipc, &                                                                                            !< component-ID of integration point
@@ -317,8 +310,6 @@ end function constitutive_homogenizedC
 !> @brief calls microstructure function of the different constitutive models
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
- use prec, only: &
-   pReal
  use material, only: &
    phasememberAt, &
    phase_plasticity, &
@@ -337,7 +328,6 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
  use plastic_disloUCLA, only: &
    plastic_disloUCLA_dependentState
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -376,8 +366,6 @@ end subroutine constitutive_microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
                                          S, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
  use material, only: &
    phasememberAt, &
    phase_plasticity, &
@@ -408,7 +396,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
  use plastic_nonlocal, only: &
    plastic_nonlocal_LpAndItsTangent
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -488,12 +475,6 @@ end subroutine constitutive_LpAndItsTangents
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
                                          S, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
- use math, only: &
-   math_I3, &
-   math_inv33, &
-   math_det33
  use material, only: &
    phasememberAt, &
    phase_plasticity, &
@@ -515,7 +496,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
  use kinematics_thermal_expansion, only: &
    kinematics_thermal_expansion_LiAndItsTangent
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -593,10 +573,6 @@ end subroutine constitutive_LiAndItsTangents
 !> @brief  collects initial intermediate deformation gradient
 !--------------------------------------------------------------------------------------------------
 pure function constitutive_initialFi(ipc, ip, el)
- use prec, only: &
-   pReal
- use math, only: &
-   math_I3
  use material, only: &
    material_phase, &
    material_homogenizationAt, &
@@ -608,7 +584,6 @@ pure function constitutive_initialFi(ipc, ip, el)
  use kinematics_thermal_expansion, only: &
    kinematics_thermal_expansion_initialStrain
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -643,10 +618,7 @@ end function constitutive_initialFi
 !! (so far no case switch because only Hooke is implemented)
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -672,12 +644,6 @@ end subroutine constitutive_SandItsTangents
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
                                               Fe, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
- use math, only : &
-   math_mul3333xx33, &
-   math_66toSym3333, &
-   math_I3
  use material, only: &
    material_phase, &
    material_homogenizationAt, &
@@ -687,7 +653,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
    damageMapping, &
    STIFFNESS_DEGRADATION_damage_ID
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -735,8 +700,6 @@ end subroutine constitutive_hooke_SandItsTangents
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el)
- use prec, only: &
-   pReal
  use debug, only: &
    debug_level, &
    debug_constitutive, &
@@ -786,7 +749,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
  use source_thermal_externalheat, only: &
    source_thermal_externalheat_dotState
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -873,8 +835,6 @@ end subroutine constitutive_collectDotState
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
  use debug, only: &
    debug_level, &
    debug_constitutive, &
@@ -896,7 +856,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
  use source_damage_isoBrittle, only: &
    source_damage_isoBrittle_deltaState
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -944,8 +903,6 @@ end subroutine constitutive_collectDeltaState
 !> @brief returns array of constitutive results
 !--------------------------------------------------------------------------------------------------
 function constitutive_postResults(S, Fi, ipc, ip, el)
- use prec, only: &
-   pReal
  use material, only: &
    phasememberAt, &
    phase_plasticityInstance, &
@@ -990,7 +947,6 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
  use source_damage_anisoDuctile, only: &
    source_damage_anisoDuctile_postResults
 
- implicit none
  integer, intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point

src/crystallite.f90

@@ -937,8 +937,6 @@ end function crystallite_push33ToRef
 !--------------------------------------------------------------------------------------------------
 function crystallite_postResults(ipc, ip, el)
  use math, only: &
-   math_qToEuler, &
-   math_qToEulerAxisAngle, &
    math_det33, &
    math_I3, &
    inDeg

src/homogenization.f90

@@ -7,6 +7,7 @@
 module homogenization
  use prec, only: &
    pReal
+ use material
 
 !--------------------------------------------------------------------------------------------------
 ! General variables for the homogenization at a  material point
@@ -66,7 +67,8 @@ module homogenization
  public ::  &
    homogenization_init, &
    materialpoint_stressAndItsTangent, &
-   materialpoint_postResults
+   materialpoint_postResults, &
+   homogenization_results
  private :: &
    partitionDeformation, &
    updateState, &
@@ -100,7 +102,6 @@ subroutine homogenization_init
   config_deallocate, &
   config_homogenization, &
   homogenization_name
- use material
  use homogenization_mech_RGC
  use thermal_isothermal
  use thermal_adiabatic
@@ -291,7 +292,6 @@ subroutine homogenization_init
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_requested:        ', shape(materialpoint_requested)
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_converged:        ', shape(materialpoint_converged)
    write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy:     ', shape(materialpoint_doneAndHappy)
-   write(6,'(a32,1x,7(i8,1x))')   'maxSizePostResults: ', homogenization_maxSizePostResults
  endif
  flush(6)
 
@@ -316,17 +316,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
    terminallyIll
  use mesh, only: &
    mesh_element
- use material, only: &
-   plasticState, &
-   sourceState, &
-   homogState, &
-   thermalState, &
-   damageState, &
-   phase_Nsources, &
-   material_homogenizationAt, &
-   mappingHomogenization, &
-   phaseAt, phasememberAt, &
-   homogenization_Ngrains
  use crystallite, only: &
    crystallite_F0, &
    crystallite_Fp0, &
@@ -658,16 +647,6 @@ subroutine materialpoint_postResults
    FEsolving_execIP
  use mesh, only: &
    mesh_element
- use material, only: &
-   material_homogenizationAt, &
-   homogState, &
-   thermalState, &
-   damageState, &
-   plasticState, &
-   sourceState, &
-   material_phase, &
-   homogenization_Ngrains, &
-   microstructure_crystallite
  use crystallite, only: &
    crystallite_sizePostResults, &
    crystallite_postResults
@@ -723,12 +702,6 @@ end subroutine materialpoint_postResults
 subroutine partitionDeformation(ip,el)
  use mesh, only: &
    mesh_element
- use material, only: &
-   homogenization_type, &
-   homogenization_Ngrains, &
-   HOMOGENIZATION_NONE_ID, &
-   HOMOGENIZATION_ISOSTRAIN_ID, &
-   HOMOGENIZATION_RGC_ID
  use crystallite, only: &
    crystallite_partionedF
  use homogenization_mech_RGC, only: &
@@ -767,14 +740,6 @@ end subroutine partitionDeformation
 function updateState(ip,el)
  use mesh, only: &
    mesh_element
- use material, only: &
-   homogenization_type, &
-   thermal_type, &
-   damage_type, &
-   homogenization_Ngrains, &
-   HOMOGENIZATION_RGC_ID, &
-   THERMAL_adiabatic_ID, &
-   DAMAGE_local_ID
  use crystallite, only: &
    crystallite_P, &
    crystallite_dPdF, &
@@ -835,13 +800,6 @@ end function updateState
 subroutine averageStressAndItsTangent(ip,el)
  use mesh, only: &
    mesh_element
- use material, only: &
-   homogenization_type, &
-   homogenization_typeInstance, &
-   homogenization_Ngrains, &
-   HOMOGENIZATION_NONE_ID, &
-   HOMOGENIZATION_ISOSTRAIN_ID, &
-   HOMOGENIZATION_RGC_ID
  use crystallite, only: &
    crystallite_P,crystallite_dPdF
  use homogenization_mech_RGC, only: &
@@ -884,27 +842,6 @@ end subroutine averageStressAndItsTangent
 function postResults(ip,el)
  use mesh, only: &
    mesh_element
- use material, only: &
-   thermalMapping, &
-   thermal_typeInstance, &
-   material_homogenizationAt, &
-   homogenization_typeInstance,&
-   mappingHomogenization, &
-   homogState, &
-   thermalState, &
-   damageState, &
-   homogenization_type, &
-   thermal_type, &
-   damage_type, &
-   HOMOGENIZATION_NONE_ID, &
-   HOMOGENIZATION_ISOSTRAIN_ID, &
-   HOMOGENIZATION_RGC_ID, &
-   THERMAL_isothermal_ID, &
-   THERMAL_adiabatic_ID, &
-   THERMAL_conduction_ID, &
-   DAMAGE_none_ID, &
-   DAMAGE_local_ID, &
-   DAMAGE_nonlocal_ID
  use homogenization_mech_RGC, only: &
    homogenization_RGC_postResults
  use thermal_adiabatic, only: &
@@ -969,4 +906,41 @@ function postResults(ip,el)
 
 end function postResults
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief writes homogenization results to HDF5 output file
+!--------------------------------------------------------------------------------------------------
+subroutine homogenization_results
+#if defined(PETSc) || defined(DAMASK_HDF5)
+  use results
+  use homogenization_mech_RGC
+  use HDF5_utilities
+  use config, only: &
+    config_name_homogenization => homogenization_name                                               ! anticipate logical name
+   
+  use material, only: &
+    homogenization_typeInstance, &
+    material_homogenization_type => homogenization_type
+    
+  integer :: p
+  character(len=256) :: group
+                                             
+  do p=1,size(config_name_homogenization)
+    group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))
+    call HDF5_closeGroup(results_addGroup(group))
+    
+    group = trim(group)//'/mech'
+    
+    call HDF5_closeGroup(results_addGroup(group))  
+    select case(material_homogenization_type(p))
+      case(HOMOGENIZATION_rgc_ID)
+        call mech_RGC_results(homogenization_typeInstance(p),group)
+    end select
+  
+ enddo   
+#endif
+
+
+end subroutine homogenization_results
+
 end module homogenization

src/homogenization_mech_RGC.f90

@@ -9,6 +9,7 @@
 module homogenization_mech_RGC
  use prec, only: &
    pReal
+ use material
 
  implicit none
  private
@@ -75,7 +76,8 @@ module homogenization_mech_RGC
    homogenization_RGC_partitionDeformation, &
    homogenization_RGC_averageStressAndItsTangent, &
    homogenization_RGC_updateState, &
-   homogenization_RGC_postResults
+   homogenization_RGC_postResults, &
+   mech_RGC_results ! name suited for planned submodule situation
  private :: &
    relaxationVector, &
    interfaceNormal, &
@@ -104,18 +106,6 @@ subroutine homogenization_RGC_init()
    INRAD
  use IO, only: &
    IO_error
- use material, only: &
-#ifdef DEBUG
-   mappingHomogenization, &
-#endif
-   homogenization_type, &
-   material_homogenizationAt, &
-   homogState, &
-   HOMOGENIZATION_RGC_ID, &
-   HOMOGENIZATION_RGC_LABEL, &
-   homogenization_typeInstance, &
-   homogenization_Noutput, &
-   homogenization_Ngrains
  use config, only: &
    config_homogenization
 
@@ -322,10 +312,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
 #endif
  use math, only: &
    math_invert2
- use material, only: &
-   material_homogenizationAt, &
-   homogenization_typeInstance, &
-   mappingHomogenization
  use numerics, only: &
    absTol_RGC, &
    relTol_RGC, &
@@ -1109,6 +1095,54 @@ pure function homogenization_RGC_postResults(instance,of) result(postResults)
 end function homogenization_RGC_postResults
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief writes results to HDF5 output file
+! ToDo: check wheter units are correct
+!--------------------------------------------------------------------------------------------------
+subroutine mech_RGC_results(instance,group)
+#if defined(PETSc) || defined(DAMASK_HDF5)
+  use results, only: &
+    results_writeDataset
+
+  integer, intent(in) :: instance
+  character(len=*) :: group
+  integer :: o
+  
+  associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
+
+  outputsLoop: do o = 1,size(prm%outputID)
+    select case(prm%outputID(o))
+    
+      case (constitutivework_ID)
+        call results_writeDataset(group,stt%work,'W',&
+                                  'work density','J/m³')
+      case (magnitudemismatch_ID)
+        call results_writeDataset(group,dst%mismatch,'N',&
+                                  'average mismatch tensor','1')
+      case (penaltyenergy_ID)
+        call results_writeDataset(group,stt%penaltyEnergy,'R',&
+                                  'mismatch penalty density','J/m³')
+      case (volumediscrepancy_ID)
+        call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',&
+                                  'volume discrepancy','m³')
+      case (maximumrelaxrate_ID)
+        call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',&
+                                  'maximum relaxation rate','m/s')
+      case (averagerelaxrate_ID)
+        call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',&
+                                  'average relaxation rate','m/s')
+    end select
+  enddo outputsLoop
+  end associate
+  
+#else
+  integer, intent(in) :: instance
+  character(len=*) :: group
+#endif
+
+end subroutine mech_RGC_results
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief collect relaxation vectors of an interface
 !--------------------------------------------------------------------------------------------------

src/homogenization_mech_isostrain.f90

@@ -36,13 +36,6 @@ module subroutine mech_isostrain_init
     debug_levelBasic
   use IO, only: &
     IO_error
-  use material, only: &
-    homogenization_type, &
-    material_homogenizationAt, &
-    homogState, &
-    HOMOGENIZATION_ISOSTRAIN_ID, &
-    HOMOGENIZATION_ISOSTRAIN_LABEL, &
-    homogenization_typeInstance
   use config, only: &
     config_homogenization
   

src/homogenization_mech_none.f90

@@ -20,13 +20,7 @@ module subroutine mech_none_init
     debug_levelBasic
   use config, only: &
     config_homogenization
-  use material, only: &
+
-    homogenization_type, &
-    material_homogenizationAt, &
-    homogState, &
-    HOMOGENIZATION_NONE_LABEL, &
-    HOMOGENIZATION_NONE_ID
- 
   implicit none
   integer :: &
     Ninstance, &

src/list.f90

@@ -70,7 +70,6 @@ subroutine add(this,string)
     IO_lc, &
     IO_stringPos
 
-  implicit none
   class(tPartitionedStringList),  target, intent(in) :: this
   character(len=*),                       intent(in) :: string
   type(tPartitionedStringList),   pointer            :: new, temp
@@ -95,7 +94,6 @@ end subroutine add
 !--------------------------------------------------------------------------------------------------
 subroutine show(this)
 
-  implicit none
   class(tPartitionedStringList), target, intent(in) :: this
   type(tPartitionedStringList),  pointer            :: item
 
@@ -114,7 +112,6 @@ end subroutine show
 !--------------------------------------------------------------------------------------------------
 subroutine free(this)
 
-  implicit none
   class(tPartitionedStringList),  intent(inout) :: this
 
   if(associated(this%next)) deallocate(this%next)
@@ -128,7 +125,6 @@ end subroutine free
 !--------------------------------------------------------------------------------------------------
 recursive subroutine finalize(this)
 
-  implicit none
   type(tPartitionedStringList),  intent(inout) :: this
 
   if(associated(this%next)) deallocate(this%next)
@@ -142,7 +138,6 @@ end subroutine finalize
 !--------------------------------------------------------------------------------------------------
 subroutine finalizeArray(this)
 
-  implicit none
   integer :: i
   type(tPartitionedStringList),  intent(inout), dimension(:) :: this
   type(tPartitionedStringList),  pointer :: temp ! bug in Gfortran?
@@ -165,7 +160,6 @@ logical function keyExists(this,key)
   use IO, only: &
     IO_stringValue
 
-  implicit none
   class(tPartitionedStringList), target, intent(in) :: this
   character(len=*),                      intent(in) :: key
   type(tPartitionedStringList),  pointer            :: item
@@ -189,8 +183,6 @@ integer function countKeys(this,key)
   use IO, only: &
     IO_stringValue
 
-  implicit none
-
   class(tPartitionedStringList), target, intent(in) :: this
   character(len=*),                      intent(in) :: key
   type(tPartitionedStringList),  pointer            :: item
@@ -218,7 +210,6 @@ real(pReal) function getFloat(this,key,defaultVal)
     IO_stringValue, &
     IO_FloatValue
 
-  implicit none
   class(tPartitionedStringList), target, intent(in)           :: this
   character(len=*),                      intent(in)           :: key
   real(pReal),                           intent(in), optional :: defaultVal
@@ -255,7 +246,6 @@ integer function getInt(this,key,defaultVal)
     IO_stringValue, &
     IO_IntValue
 
-  implicit none
   class(tPartitionedStringList), target, intent(in)           :: this
   character(len=*),                      intent(in)           :: key
   integer,                               intent(in), optional :: defaultVal
@@ -292,7 +282,6 @@ character(len=65536) function getString(this,key,defaultVal,raw)
     IO_error, &
     IO_stringValue
  
-  implicit none
   class(tPartitionedStringList), target, intent(in)           :: this
   character(len=*),                      intent(in)           :: key
   character(len=*),                      intent(in), optional :: defaultVal
@@ -343,7 +332,6 @@ function getFloats(this,key,defaultVal,requiredSize)
     IO_stringValue, &
     IO_FloatValue
  
-  implicit none
   real(pReal),     dimension(:), allocatable          :: getFloats
   class(tPartitionedStringList), target, intent(in)   :: this
   character(len=*),              intent(in)           :: key
@@ -393,7 +381,6 @@ function getInts(this,key,defaultVal,requiredSize)
     IO_stringValue, &
     IO_IntValue
  
-  implicit none
   integer, dimension(:), allocatable                          :: getInts
   class(tPartitionedStringList), target, intent(in)           :: this
   character(len=*),                      intent(in)           :: key
@@ -443,7 +430,6 @@ function getStrings(this,key,defaultVal,raw)
     IO_error, &
     IO_StringValue
  
-  implicit none
   character(len=65536),dimension(:), allocatable             :: getStrings
   class(tPartitionedStringList),target, intent(in)           :: this
   character(len=*),                     intent(in)           :: key

src/material.f90

@@ -8,15 +8,9 @@
 !! 'phase', 'texture', and 'microstucture'
 !--------------------------------------------------------------------------------------------------
 module material
- use prec, only: &
+ use prec
-   pReal, &
+ use math
-   pInt, &
+ use config
-   tState, &
-   tPlasticState, &
-   tSourceState, &
-   tHomogMapping, &
-   group_float, &
-   group_int
 
  implicit none
  private
@@ -53,55 +47,35 @@ module material
 
  enum, bind(c)
    enumerator :: ELASTICITY_undefined_ID, &
-                 ELASTICITY_hooke_ID
+                 ELASTICITY_hooke_ID, &
- end enum
+                 PLASTICITY_undefined_ID, &
- enum, bind(c)
-   enumerator :: PLASTICITY_undefined_ID, &
                  PLASTICITY_none_ID, &
                  PLASTICITY_isotropic_ID, &
                  PLASTICITY_phenopowerlaw_ID, &
                  PLASTICITY_kinehardening_ID, &
                  PLASTICITY_dislotwin_ID, &
                  PLASTICITY_disloucla_ID, &
-                 PLASTICITY_nonlocal_ID
+                 PLASTICITY_nonlocal_ID, &
- end enum
+                 SOURCE_undefined_ID, &
-
- enum, bind(c)
-   enumerator :: SOURCE_undefined_ID, &
                  SOURCE_thermal_dissipation_ID, &
                  SOURCE_thermal_externalheat_ID, &
                  SOURCE_damage_isoBrittle_ID, &
                  SOURCE_damage_isoDuctile_ID, &
                  SOURCE_damage_anisoBrittle_ID, &
-                 SOURCE_damage_anisoDuctile_ID
+                 SOURCE_damage_anisoDuctile_ID, &
- end enum
+                 KINEMATICS_undefined_ID, &
-
- enum, bind(c)
-   enumerator :: KINEMATICS_undefined_ID, &
                  KINEMATICS_cleavage_opening_ID, &
                  KINEMATICS_slipplane_opening_ID, &
-                 KINEMATICS_thermal_expansion_ID
+                 KINEMATICS_thermal_expansion_ID, &
- end enum
+                 STIFFNESS_DEGRADATION_undefined_ID, &
-
+                 STIFFNESS_DEGRADATION_damage_ID, &
- enum, bind(c)
+                 THERMAL_isothermal_ID, &
-   enumerator :: STIFFNESS_DEGRADATION_undefined_ID, &
-                 STIFFNESS_DEGRADATION_damage_ID
- end enum
-
- enum, bind(c)
-   enumerator :: THERMAL_isothermal_ID, &
                  THERMAL_adiabatic_ID, &
-                 THERMAL_conduction_ID
+                 THERMAL_conduction_ID, &
- end enum
+                 DAMAGE_none_ID, &
-
- enum, bind(c)
-   enumerator :: DAMAGE_none_ID, &
                  DAMAGE_local_ID, &
-                 DAMAGE_nonlocal_ID
+                 DAMAGE_nonlocal_ID, &
- end enum
+                 HOMOGENIZATION_undefined_ID, &
-
- enum, bind(c)
-   enumerator :: HOMOGENIZATION_undefined_ID, &
                  HOMOGENIZATION_none_ID, &
                  HOMOGENIZATION_isostrain_ID, &
                  HOMOGENIZATION_rgc_ID
@@ -279,20 +253,9 @@ subroutine material_init
    debug_material, &
    debug_levelBasic, &
    debug_levelExtensive
- use config, only: &
-   config_crystallite, &
-   config_homogenization, &
-   config_microstructure, &
-   config_phase, &
-   config_texture, &
-   homogenization_name, &
-   microstructure_name, &
-   phase_name, &
-   texture_name
  use mesh, only: &
    theMesh
 
- implicit none
  integer(pInt), parameter :: FILEUNIT = 210_pInt
  integer(pInt)            :: m,c,h, myDebug, myPhase, myHomog
  integer(pInt) :: &
@@ -461,14 +424,11 @@ end subroutine material_init
 !> @brief parses the homogenization part from the material configuration
 !--------------------------------------------------------------------------------------------------
 subroutine material_parseHomogenization
- use config, only : &
-   config_homogenization
  use mesh, only: &
    theMesh
  use IO, only: &
    IO_error
 
- implicit none
  integer(pInt)        :: h
  character(len=65536) :: tag
 
@@ -559,21 +519,15 @@ end subroutine material_parseHomogenization
 !> @brief parses the microstructure part in the material configuration file
 !--------------------------------------------------------------------------------------------------
 subroutine material_parseMicrostructure
- use prec, only: &
-  dNeq
  use IO, only: &
    IO_floatValue, &
    IO_intValue, &
    IO_stringValue, &
    IO_stringPos, &
    IO_error
- use config, only: &
-   config_microstructure, &
-   microstructure_name
  use mesh, only: &
    theMesh
 
- implicit none
  character(len=65536), dimension(:), allocatable :: &
    strings
  integer(pInt), allocatable, dimension(:) :: chunkPos
@@ -637,10 +591,7 @@ end subroutine material_parseMicrostructure
 !> @brief parses the crystallite part in the material configuration file
 !--------------------------------------------------------------------------------------------------
 subroutine material_parseCrystallite
- use config, only: &
-   config_crystallite
 
- implicit none
  integer(pInt)        :: c
 
  allocate(crystallite_Noutput(size(config_crystallite)),source=0_pInt)
@@ -659,10 +610,7 @@ subroutine material_parsePhase
    IO_error, &
    IO_getTag, &
    IO_stringValue
- use config, only: &
-   config_phase
 
- implicit none
  integer(pInt) :: sourceCtr, kinematicsCtr, stiffDegradationCtr, p
  character(len=65536), dimension(:), allocatable ::  str 
 
@@ -785,23 +733,12 @@ end subroutine material_parsePhase
 !> @brief parses the texture part in the material configuration file
 !--------------------------------------------------------------------------------------------------
 subroutine material_parseTexture
- use prec, only: &
-   dNeq
  use IO, only: &
    IO_error, &
    IO_stringPos, &
    IO_floatValue, &
    IO_stringValue
- use config, only: &
-   config_deallocate, &
-   config_texture
- use math, only: &
-   inRad, &
-   math_sampleRandomOri, &
-   math_I3, &
-   math_det33
 
- implicit none
  integer(pInt) :: section, gauss, j, t, i
  character(len=65536), dimension(:), allocatable ::  strings                                     ! Values for given key in material config 
  integer(pInt), dimension(:), allocatable :: chunkPos
@@ -880,7 +817,6 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
  use numerics, only: &
    numerics_integrator
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    NofMyPhase, &
@@ -928,7 +864,6 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
  use numerics, only: &
    numerics_integrator
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    of, &
@@ -967,47 +902,15 @@ end subroutine material_allocateSourceState
 !! calculates the volume of the grains and deals with texture components
 !--------------------------------------------------------------------------------------------------
 subroutine material_populateGrains
- use math, only: &
-   math_EulertoR, &
-   math_RtoEuler, &
-   math_mul33x33, &
-   math_range
  use mesh, only: &
    theMesh
- use config, only: &
-   config_homogenization, &
-   config_microstructure, &
-   config_deallocate
- use IO, only: &
-   IO_error
-
- implicit none
- integer(pInt), dimension (:,:), allocatable :: Ngrains
- integer(pInt), dimension (microstructure_maxNconstituents)  :: &
-   NgrainsOfConstituent, &
-   currentGrainOfConstituent, &
-   randomOrder
- real(pReal), dimension (microstructure_maxNconstituents)  :: &
-   rndArray
- real(pReal), dimension (:,:),   allocatable :: orientationOfGrain
- real(pReal), dimension (3)                  :: orientation
- integer(pInt), dimension (:),   allocatable :: phaseOfGrain, textureOfGrain
- integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, &
-                  phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, &
-                  grain,constituentGrain,ipGrain,ip
- real(pReal) :: deviation,extreme,rnd
- integer(pInt),         dimension (:,:), allocatable :: Nelems                                      ! counts number of elements in homog, micro array
- type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro                                 ! lists element number in homog, micro array
 
+ integer(pInt) :: e,i,c,homog,micro
 
  allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),        source=0_pInt)
  allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),      source=0_pInt)
  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)
-
-
   do e = 1, theMesh%Nelems
     do i = 1, theMesh%elem%nIPs
       homog = theMesh%homogenizationAt(e)
@@ -1017,7 +920,7 @@ subroutine material_populateGrains
         material_texture(c,i,e) = microstructure_texture(c,micro)
         material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e))
         material_EulerAngles(1:3,c,i,e) = math_RtoEuler( &                                       ! translate back to Euler angles
-                                             math_mul33x33( &                                       ! pre-multiply
+                                             matmul( &                                       ! pre-multiply
                                                math_EulertoR(material_EulerAngles(1:3,c,i,e)), &    ! face-value orientation
                                                texture_transformation(1:3,1:3,material_texture(c,i,e)) &          ! and transformation matrix
                                              ) &
@@ -1026,209 +929,6 @@ subroutine material_populateGrains
     enddo
   enddo
 
- return
-!--------------------------------------------------------------------------------------------------
-! precounting of elements for each homog/micro pair
- do e = 1_pInt, theMesh%Nelems
-   homog = theMesh%homogenizationAt(e)
-   micro = theMesh%microstructureAt(e)
-   Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
- enddo
- allocate(elemsOfHomogMicro(size(config_homogenization),size(config_microstructure)))
- do homog = 1,size(config_homogenization)
-   do micro = 1,size(config_microstructure)
-     if (Nelems(homog,micro) > 0_pInt) then
-       allocate(elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)))
-       elemsOfHomogMicro(homog,micro)%p = 0_pInt
-    endif
-   enddo
- enddo
-
-!--------------------------------------------------------------------------------------------------
-! 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,theMesh%Nelems
-   homog = theMesh%homogenizationAt(e)
-   micro = theMesh%microstructureAt(e)
-   if (homog < 1_pInt .or. homog > size(config_homogenization)) &                                      ! out of bounds
-     call IO_error(154_pInt,e,0_pInt,0_pInt)
-   if (micro < 1_pInt .or. micro > size(config_microstructure)) &                                      ! out of bounds
-     call IO_error(155_pInt,e,0_pInt,0_pInt)
-   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) * 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
-   elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e                                        ! remember elements active in this homog/micro pair
- enddo elementLooping
-
- allocate(phaseOfGrain(maxval(Ngrains)),        source=0_pInt)                                      ! reserve memory for maximum case
- allocate(textureOfGrain(maxval(Ngrains)),      source=0_pInt)                                      ! reserve memory for maximum case
- allocate(orientationOfGrain(3,maxval(Ngrains)),source=0.0_pReal)                                   ! reserve memory for maximum case
-
-
- homogenizationLoop: do homog = 1_pInt,size(config_homogenization)
-   dGrains = homogenization_Ngrains(homog)                                                          ! grain number per material point
-   microstructureLoop: do micro = 1_pInt,size(config_microstructure)                                ! all pairs of homog and micro
-     activePair: if (Ngrains(homog,micro) > 0_pInt) then
-       myNgrains = Ngrains(homog,micro)                                                             ! assign short name for total number of grains to populate
-       myNconstituents = microstructure_Nconstituents(micro)                                        ! assign short name for number of constituents
-
-!--------------------------------------------------------------------------------------------------
-! divide myNgrains as best over constituents
-!
-! example: three constituents with fractions of 0.25, 0.25, and 0.5 distributed over 20 (microstructure) grains
-!
-!                       ***** ***** **********
-! NgrainsOfConstituent: 5,    5,    10
-! counters:
-!                      |-----> grain (if constituent == 2)
-!                            |--> constituentGrain (of constituent 2)
-!
-
-       NgrainsOfConstituent = 0_pInt                                                                ! reset counter of grains per constituent
-       forall (i = 1_pInt:myNconstituents) &
-         NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro)*real(myNgrains,pReal),pInt)! do rounding integer conversion
-       do while (sum(NgrainsOfConstituent) /= myNgrains)                                            ! total grain count over constituents wrong?
-         sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent))                                  ! direction of required change
-         extreme = 0.0_pReal
-         t = 0_pInt
-         do i = 1_pInt,myNconstituents                                                              ! find largest deviator
-           deviation = real(sgn,pReal)*log( microstructure_fraction(i,micro) / &
-                                           !-------------------------------- &
-                                           (real(NgrainsOfConstituent(i),pReal)/real(myNgrains,pReal) ) )
-           if (deviation > extreme) then
-             extreme = deviation
-             t = i
-           endif
-         enddo
-         NgrainsOfConstituent(t) = NgrainsOfConstituent(t) + sgn                                    ! change that by one
-       enddo
-
-!--------------------------------------------------------------------------------------------------
-! assign phase and texture info
-
-       phaseOfGrain = 0_pInt
-       textureOfGrain = 0_pInt
-       orientationOfGrain = 0.0_pReal
-
-       texture: do i = 1_pInt,myNconstituents                                                       ! loop over constituents
-         grain            = sum(NgrainsOfConstituent(1_pInt:i-1_pInt))                              ! set microstructure grain index of current constituent
-                                                                                                    ! "grain" points to start of this constituent's grain population
-         constituentGrain = 0_pInt                                                                  ! constituent grain index
-
-         phaseID   = microstructure_phase(i,micro)
-         textureID = microstructure_texture(i,micro)
-         phaseOfGrain  (grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID                       ! assign resp. phase
-         textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID                     ! assign resp. texture
-
-         myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/1.0,pInt)                    ! max number of unique orientations (excl. symmetry)
-
-!--------------------------------------------------------------------------------------------------
-! has texture components
-         gauss: do t = 1_pInt,texture_Ngauss(textureID)                                             ! loop over Gauss components
-           do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt)         ! loop over required grain count
-             orientationOfGrain(:,grain+constituentGrain+g) = texture_Gauss(1:3,t,textureID)
-           enddo
-           constituentGrain = &
-           constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID))         ! advance counter for grains of current constituent
-         enddo gauss
-
-
-!--------------------------------------------------------------------------------------------------
-! ...texture transformation
-
-         do j = 1_pInt,myNorientations                                                              ! loop over each "real" orientation
-           orientationOfGrain(1:3,grain+j) = math_RtoEuler( &                                       ! translate back to Euler angles
-                                             math_mul33x33( &                                       ! pre-multiply
-                                               math_EulertoR(orientationOfGrain(1:3,grain+j)), &    ! face-value orientation
-                                               texture_transformation(1:3,1:3,textureID) &          ! and transformation matrix
-                                             ) &
-                                             )
-         enddo
-
-!--------------------------------------------------------------------------------------------------
-! shuffle grains within current constituent
-
-         do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt                                               ! walk thru grains of current constituent
-           call random_number(rnd)
-           t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt)         ! select a grain in remaining list
-           m                               = phaseOfGrain(grain+t)                                  ! exchange current with random
-           phaseOfGrain(grain+t)           = phaseOfGrain(grain+j)
-           phaseOfGrain(grain+j)           = m
-           m                               = textureOfGrain(grain+t)                                ! exchange current with random
-           textureOfGrain(grain+t)         = textureOfGrain(grain+j)
-           textureOfGrain(grain+j)         = m
-           orientation                     = orientationOfGrain(1:3,grain+t)                        ! exchange current with random
-           orientationOfGrain(1:3,grain+t) = orientationOfGrain(1:3,grain+j)
-           orientationOfGrain(1:3,grain+j) = orientation
-         enddo
-
-       enddo texture
-!< @todo calc fraction after weighing with volumePerGrain, exchange in MC steps to improve result (humbug at the moment)
-
-
-
-!--------------------------------------------------------------------------------------------------
-! distribute grains of all constituents as accurately as possible to given constituent fractions
-
-       ip = 0_pInt
-       currentGrainOfConstituent = 0_pInt
-
-       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
-         if (microstructure_elemhomo(micro)) then                                                   ! homogeneous distribution of grains over each element's IPs
-           m = 1_pInt                                                                               ! process only first IP
-         else
-           m = theMesh%elem%nIPs
-         endif
-
-         do i = 1_pInt, m                                                                           ! loop over necessary IPs
-           ip = ip + 1_pInt                                                                         ! keep track of total ip count
-           ipGrain = 0_pInt                                                                         ! count number of grains assigned at this IP
-           randomOrder = math_range(microstructure_maxNconstituents)                                ! start out with ordered sequence of constituents
-           call random_number(rndArray)                                                             ! as many rnd numbers as (max) constituents
-           do j = 1_pInt, myNconstituents - 1_pInt                                                  ! loop over constituents ...
-             r = nint(rndArray(j)*real(myNconstituents-j,pReal)+real(j,pReal)+0.5_pReal,pInt)       ! ... select one in remaining list
-             c = randomOrder(r)                                                                     ! ... call it "c"
-             randomOrder(r) = randomOrder(j)                                                        ! ... and exchange with present position in constituent list
-             grain = sum(NgrainsOfConstituent(1:c-1_pInt))                                          ! figure out actual starting index in overall/consecutive grain population
-             do g = 1_pInt, min(dGrains-ipGrain, &                                                  ! leftover number of grains at this IP
-                                max(0_pInt, &                                                       ! no negative values
-                                    nint(real(ip * dGrains * NgrainsOfConstituent(c)) / &           ! fraction of grains scaled to this constituent...
-                                         real(myNgrains),pInt) - &                                  ! ...minus those already distributed
-                                         currentGrainOfConstituent(c)))
-               ipGrain = ipGrain + 1_pInt                                                           ! advance IP grain counter
-               currentGrainOfConstituent(c)  = currentGrainOfConstituent(c) + 1_pInt                ! advance index of grain population for constituent c
-               material_phase(ipGrain,i,e)   = phaseOfGrain(grain+currentGrainOfConstituent(c))
-               material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c))
-               material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c))
-           enddo; enddo
-
-           c = randomOrder(microstructure_Nconstituents(micro))                                     ! look up constituent remaining after random shuffling
-           grain = sum(NgrainsOfConstituent(1:c-1_pInt))                                            ! figure out actual starting index in overall/consecutive grain population
-           do ipGrain = ipGrain + 1_pInt, dGrains                                                   ! ensure last constituent fills up to dGrains
-             currentGrainOfConstituent(c)  = currentGrainOfConstituent(c) + 1_pInt
-             material_phase(ipGrain,i,e)   = phaseOfGrain(grain+currentGrainOfConstituent(c))
-             material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c))
-             material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c))
-           enddo
-
-         enddo
-
-         do i = i, theMesh%elem%nIPs                                                                 ! loop over IPs to (possibly) distribute copies from first IP
-           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)
-           material_EulerAngles(1:3,1_pInt:dGrains,i,e) = material_EulerAngles(1:3,1_pInt:dGrains,1,e)
-         enddo
-
-       enddo
-     endif activePair
-   enddo microstructureLoop
- enddo homogenizationLoop
-
  deallocate(texture_transformation)
 
  call config_deallocate('material.config/microstructure')

src/mesh/FEM_mech.f90

@@ -55,8 +55,7 @@ module FEM_mech
  public :: &
    FEM_mech_init, &
    FEM_mech_solution ,&
-   FEM_mech_forward, &
+   FEM_mech_forward
-   FEM_mech_destroy
 contains
 
 !--------------------------------------------------------------------------------------------------
@@ -583,6 +582,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
 
 end subroutine FEM_mech_formJacobian
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief forwarding routine
 !--------------------------------------------------------------------------------------------------
@@ -655,7 +655,6 @@ end subroutine FEM_mech_forward
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
  use numerics, only: &
-   worldrank, &
    err_struct_tolAbs, &
    err_struct_tolRel
  use IO, only: &
@@ -677,30 +676,13 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
  call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
  CHKERRQ(ierr)
  if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
- if (worldrank == 0) then 
+ write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), &
-   write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), &
+                 ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
-                   ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
+                                                 int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
-                                                   int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
+ write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
-   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
                                                        transpose(P_av)*1.e-6_pReal
-   flush(6) 
+ flush(6)
- endif
  
 end subroutine FEM_mech_converged
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief destroy routine
-!--------------------------------------------------------------------------------------------------
-subroutine FEM_mech_destroy()
-
- implicit none
- PetscErrorCode :: ierr
-
- call VecDestroy(solution,ierr); CHKERRQ(ierr)
- call VecDestroy(solution_rate,ierr); CHKERRQ(ierr)
- call SNESDestroy(mech_snes,ierr); CHKERRQ(ierr)
-
-end subroutine FEM_mech_destroy
-
 end module FEM_mech

src/mesh/FEM_utilities.f90

@@ -23,7 +23,6 @@ use PETScis
 !--------------------------------------------------------------------------------------------------
 ! grid related information information
  real(pReal),   public                :: wgt                                                        !< weighting factor 1/Nelems
- real(pReal),   public                :: wgtDof                                                     !< weighting factor 1/Nelems
   
 !--------------------------------------------------------------------------------------------------
 ! output data
@@ -35,33 +34,18 @@ use PETScis
 
  enum, bind(c)
    enumerator :: FIELD_UNDEFINED_ID, &
-                 FIELD_MECH_ID, &
+                 FIELD_MECH_ID
-                 FIELD_THERMAL_ID, &
-                 FIELD_DAMAGE_ID, &
-                 FIELD_SOLUTE_ID, &
-                 FIELD_MGTWIN_ID
  end enum
  enum, bind(c)
    enumerator :: COMPONENT_UNDEFINED_ID, &
                  COMPONENT_MECH_X_ID, &
                  COMPONENT_MECH_Y_ID, &
-                 COMPONENT_MECH_Z_ID, &
+                 COMPONENT_MECH_Z_ID
-                 COMPONENT_THERMAL_T_ID, &
-                 COMPONENT_DAMAGE_PHI_ID, &
-                 COMPONENT_SOLUTE_CV_ID, &
-                 COMPONENT_SOLUTE_CVPOT_ID, &
-                 COMPONENT_SOLUTE_CH_ID, &
-                 COMPONENT_SOLUTE_CHPOT_ID, &
-                 COMPONENT_SOLUTE_CVaH_ID, &
-                 COMPONENT_SOLUTE_CVaHPOT_ID, &
-                 COMPONENT_MGTWIN_PHI_ID
  end enum
  
 !--------------------------------------------------------------------------------------------------
 ! variables controlling debugging
  logical, private :: &
-   debugGeneral, &                                                                                  !< general debugging of FEM solver
-   debugRotation, &                                                                                 !< also printing out results in lab frame
    debugPETSc                                                                                       !< use some in debug defined options for more verbose PETSc solution
 
 !--------------------------------------------------------------------------------------------------
@@ -111,36 +95,26 @@ use PETScis
  public :: &
    utilities_init, &
    utilities_constitutiveResponse, &
-   utilities_indexBoundaryDofs, &
    utilities_projectBCValues, &
-   utilities_indexActiveSet, &
-   utilities_destroy, &
    FIELD_MECH_ID, &
    COMPONENT_MECH_X_ID, &
    COMPONENT_MECH_Y_ID, &
-   COMPONENT_MECH_Z_ID, &
+   COMPONENT_MECH_Z_ID
-   COMPONENT_THERMAL_T_ID
 
 contains 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all neccessary fields, sets debug flags
 !--------------------------------------------------------------------------------------------------
-subroutine utilities_init()
+subroutine utilities_init
  use numerics, only: &                        
    structOrder, &
-   integrationOrder, &
-   worldsize, & 
-   worldrank, &
    petsc_defaultOptions, &
    petsc_options
  use debug, only: &
    debug_level, &
    debug_SPECTRAL, &
-   debug_LEVELBASIC, &
+   debug_SPECTRALPETSC,&
-   debug_SPECTRALPETSC, &
-   debug_SPECTRALROTATION
- use debug, only: &
    PETSCDEBUG
  use math                                                                                           ! must use the whole module for use of FFTW
  use mesh, only: &
@@ -151,16 +125,12 @@ subroutine utilities_init()
  implicit none
 
  character(len=1024)                :: petsc_optionsPhysics
- integer(pInt)                      :: dimPlex
- PetscInt                           :: dim
  PetscErrorCode                     :: ierr
 
  write(6,'(/,a)')   ' <<<+-  DAMASK_FEM_utilities init  -+>>>'
  
 !--------------------------------------------------------------------------------------------------
 ! set debugging parameters
- debugGeneral    = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC)         /= 0
- debugRotation   = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION)   /= 0
  debugPETSc      = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC)      /= 0
  if(debugPETSc) write(6,'(3(/,a),/)') &
                 ' Initializing PETSc with debug options: ', &
@@ -180,7 +150,6 @@ subroutine utilities_init()
  
  wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
 
- call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr) 
 
 end subroutine utilities_init
 
@@ -188,21 +157,13 @@ end subroutine utilities_init
 !> @brief calculates constitutive response
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
- use debug, only: &
-   debug_reset, &
-   debug_info
  use math, only: &
-   math_transpose33, &
-   math_rotate_forward33, &
    math_det33
  use FEsolving, only: &
    restartWrite
  use homogenization, only: &
-   materialpoint_F, &
    materialpoint_P, &
    materialpoint_stressAndItsTangent
- use mesh, only: &
-   mesh_NcpElems  
  
  implicit none
  real(pReal), intent(in)                                         :: timeinc                         !< loading time
@@ -213,9 +174,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
  logical :: &
    age
 
- integer(pInt) :: &
-   j
- real(pReal)   :: defgradDetMin, defgradDetMax, defgradDet
  PetscErrorCode :: ierr
 
  write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
@@ -227,27 +185,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
  if (cutBack) then                                                                                  ! restore saved variables
    age = .False.
  endif
- call debug_reset()
-
-!--------------------------------------------------------------------------------------------------
-! calculate bounds of det(F) and report
- if(debugGeneral) then
-   defgradDetMax = -huge(1.0_pReal)
-   defgradDetMin = +huge(1.0_pReal)
-   do j = 1_pInt, mesh_NcpElems
-     defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
-     defgradDetMax = max(defgradDetMax,defgradDet)
-     defgradDetMin = min(defgradDetMin,defgradDet) 
-   end do
-   write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
-   write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
-   flush(6)
- endif
   
  call materialpoint_stressAndItsTangent(.true.,timeinc)                                             ! calculate P field
 
- call debug_info()
- 
  restartWrite = .false.                                                                             ! reset restartWrite status
  cutBack = .false.                                                                                  ! reset cutBack status
  
@@ -257,97 +197,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 end subroutine utilities_constitutiveResponse
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Create index sets of boundary dofs (in local and global numbering)
-!--------------------------------------------------------------------------------------------------
-subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS)
-
- implicit none
- 
- DM                                    :: dm_local
- ISLocalToGlobalMapping                :: local2global
- PetscSection                          :: section
- PetscInt                              :: nFaceSets, numFields, nDof
- IS,    dimension(nFaceSets,numFields) :: localIS, globalIS
- PetscInt                              :: field, faceSet, point, dof, offset
- PetscInt                              :: localSize, storageSize, ISSize
- PetscInt, dimension(:)  , allocatable :: localIndices
- IS                                    :: faceSetIS, BC_IS, dummyIS
- PetscInt, dimension(:)  ,     pointer :: pFaceSets, pBCvertex, pBCvertexlc
- DMLabel                               :: BCLabel
- PetscErrorCode                        :: ierr
-
- call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
- call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr)
- call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr)
- call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr)
- call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
- allocate(localIndices (storageSize))
- do faceSet = 1, nFaceSets
-   call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr)
-   CHKERRQ(ierr)
-   call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr)
-   CHKERRQ(ierr)
-   if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr)
-   do field = 1, numFields
-     localSize  = 0
-     do point = 1, ISSize
-       call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr)
-       CHKERRQ(ierr)
-       call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr)
-       CHKERRQ(ierr)
-       do dof = 1, nDof
-         localSize = localSize + 1
-         localIndices(localSize) = offset + dof - 1
-       enddo  
-     enddo
-     call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
-                          localIS(faceSet,field),ierr)
-     CHKERRQ(ierr)
-     call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), &
-                                        globalIS(faceSet,field),ierr)
-     CHKERRQ(ierr)
-   enddo
-   if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr)
-   call ISDestroy(BC_IS,ierr); CHKERRQ(ierr)
- enddo  
- call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
- call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr)
- 
- do faceSet = 1, nFaceSets; do field = 1, numFields
-   call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr)
-   if (ISSize > 0) then
-     call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
-     call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
-   endif  
-   localSize = 0
-   do point = 1, ISSize
-     if (pBCvertex(point) >= 0) then
-       localSize = localSize + 1
-       localIndices(localSize) = pBCvertexlc(point)
-     endif  
-   enddo
-   if (ISSize > 0) then
-     call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
-     call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
-   endif  
-   call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr) 
-   call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
-                        globalIS(faceSet,field),ierr) 
-   CHKERRQ(ierr)
-   if (ISSize > 0) then
-     call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr)
-     call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr)
-     call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr)
-     CHKERRQ(ierr)
-     call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
-   endif
- enddo; enddo                       
- deallocate(localIndices)
- 
-end subroutine utilities_indexBoundaryDofs
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Project BC values to local vector
 !--------------------------------------------------------------------------------------------------
@@ -384,104 +233,4 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
  
 end subroutine utilities_projectBCValues
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Create index sets of boundary dofs (in local and global numbering)
-!--------------------------------------------------------------------------------------------------
-subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS)
- use mesh, only: &
-   geomMesh
-
- implicit none
- 
- ISLocalToGlobalMapping                   :: local2global
- PetscSection                             :: section
- Vec                                      :: x_local, f_local
- PetscInt                                 :: field
- IS                                       :: localIS, globalIS, dummyIS
- PetscScalar, dimension(:)  ,     pointer :: x_scal, f_scal
- PetscInt                                 :: ISSize
- PetscInt                                 :: chart, chartStart, chartEnd, nDof, dof, offset
- PetscInt                                 :: localSize
- PetscInt,    dimension(:)  , allocatable :: localIndices
- PetscInt,    dimension(:)  ,     pointer :: pBCvertex, pBCvertexlc
- PetscErrorCode                           :: ierr
-
- call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr)
- CHKERRQ(ierr)
- call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr)
- CHKERRQ(ierr)
- call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
- call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
- localSize = 0
- do chart = chartStart, chartEnd-1
-   call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
-   call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
-   do dof = offset+1, offset+nDof
-     if (((x_scal(dof) <       1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
-         ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1
-   enddo
- enddo
- allocate(localIndices(localSize))        
- localSize = 0
- do chart = chartStart, chartEnd-1
-   call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
-   call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
-   do dof = offset+1, offset+nDof
-     if (((x_scal(dof) <       1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
-         ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then
-       localSize = localSize + 1
-       localIndices(localSize) = dof-1
-     endif  
-   enddo
- enddo
- call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
- call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
- call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr)
- CHKERRQ(ierr)
- call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr)
- CHKERRQ(ierr)
- call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr)
- if (ISSize > 0) then
-   call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
-   call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
- endif  
- localSize = 0
- do chart = 1, ISSize
-   if (pBCvertex(chart) >= 0) then
-     localSize = localSize + 1
-     localIndices(localSize) = pBCvertexlc(chart)
-   endif  
- enddo
- if (ISSize > 0) then
-   call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
-   call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
- endif  
- call ISDestroy(globalIS,ierr); CHKERRQ(ierr) 
- call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr) 
- CHKERRQ(ierr)
- if (ISSize > 0) then
-   call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr)
-   call ISDestroy(localIS,ierr); CHKERRQ(ierr)
-   call ISDifference(dummyIS,globalIS,localIS,ierr)
-   CHKERRQ(ierr)
-   call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
- endif                      
- 
-end subroutine utilities_indexActiveSet
-
-!--------------------------------------------------------------------------------------------------
-!> @brief cleans up
-!--------------------------------------------------------------------------------------------------
-subroutine utilities_destroy()
-
- !implicit none
- !PetscInt       :: homog, cryst, grain, phase 
- !PetscErrorCode :: ierr
-
- !call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
- !call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
-
-end subroutine utilities_destroy
-
-
 end module FEM_utilities

src/mesh_abaqus.f90

@@ -884,7 +884,7 @@ end subroutine mesh_abaqus_count_cpElements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_abaqus_map_elements(fileUnit)
 
- use math, only: math_qsort
+ use math, only: math_sort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -935,7 +935,7 @@ subroutine mesh_abaqus_map_elements(fileUnit)
    endselect
  enddo
 
- call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))               ! should be mesh_NcpElems
+ call math_sort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))               ! should be mesh_NcpElems
 
  if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt)
 
@@ -948,7 +948,7 @@ end subroutine mesh_abaqus_map_elements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_abaqus_map_nodes(fileUnit)
 
- use math, only: math_qsort
+ use math, only: math_sort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -999,7 +999,7 @@ subroutine mesh_abaqus_map_nodes(fileUnit)
    endif
  enddo
 
- call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
+ call math_sort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
 
  if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt)
 
@@ -1541,7 +1541,7 @@ pure function mesh_cellCenterCoordinates(ip,el)
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_build_ipAreas
  use math, only: &
-   math_crossproduct
+   math_cross
 
  implicit none
  integer(pInt) :: e,t,g,c,i,f,n,m
@@ -1576,7 +1576,7 @@ subroutine mesh_build_ipAreas
            do f = 1_pInt,FE_NipNeighbors(c)                                                         ! loop over cell faces
              forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
                nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
-             normal = math_crossproduct(nodePos(1:3,2) - nodePos(1:3,1), &
+             normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
                                          nodePos(1:3,3) - nodePos(1:3,1))
              mesh_ipArea(f,i,e) = norm2(normal)
              mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal)                             ! ensure unit length of area normal
@@ -1595,7 +1595,7 @@ subroutine mesh_build_ipAreas
                nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
              forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
                normals(1:3,n) = 0.5_pReal &
-                              * math_crossproduct(nodePos(1:3,1+mod(n  ,m)) - nodePos(1:3,n), &
+                              * math_cross(nodePos(1:3,1+mod(n  ,m)) - nodePos(1:3,n), &
                                                    nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n))
              normal = 0.5_pReal * sum(normals,2)
              mesh_ipArea(f,i,e) = norm2(normal)

src/mesh_marc.f90

@@ -640,7 +640,7 @@ subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit)
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileUnit)
 
- use math, only: math_qsort
+ use math, only: math_sort
  use IO,   only: IO_lc, &
                  IO_intValue, &
                  IO_stringValue, &
@@ -701,7 +701,7 @@ subroutine mesh_marc_map_elements(tableStyle,nameElemSet,mapElemSet,nElems,fileU
       mesh_mapFEtoCPelem(2,cpElem) = cpElem
     enddo
  
-call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                ! should be mesh_NcpElems
+call math_sort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                ! should be mesh_NcpElems
 
 end subroutine mesh_marc_map_elements
 
@@ -711,7 +711,7 @@ end subroutine mesh_marc_map_elements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_marc_map_nodes(nNodes,fileUnit)
 
- use math, only: math_qsort
+ use math, only: math_sort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -743,7 +743,7 @@ subroutine mesh_marc_map_nodes(nNodes,fileUnit)
    endif
  enddo
 
-650 call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
+650 call math_sort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
 
 end subroutine mesh_marc_map_nodes
 
@@ -1230,7 +1230,7 @@ pure function mesh_cellCenterCoordinates(ip,el)
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_build_ipAreas
  use math, only: &
-   math_crossproduct
+   math_cross
 
  implicit none
  integer(pInt) :: e,t,g,c,i,f,n,m
@@ -1265,7 +1265,7 @@ subroutine mesh_build_ipAreas
            do f = 1_pInt,FE_NipNeighbors(c)                                                         ! loop over cell faces
              forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
                nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
-             normal = math_crossproduct(nodePos(1:3,2) - nodePos(1:3,1), &
+             normal = math_cross(nodePos(1:3,2) - nodePos(1:3,1), &
                                          nodePos(1:3,3) - nodePos(1:3,1))
              mesh_ipArea(f,i,e) = norm2(normal)
              mesh_ipAreaNormal(1:3,f,i,e) = normal / norm2(normal)                             ! ensure unit length of area normal
@@ -1284,7 +1284,7 @@ subroutine mesh_build_ipAreas
                nodePos(1:3,n) = mesh_cellnode(1:3,mesh_cell(FE_cellface(n,f,c),i,e))
              forall(n = 1_pInt:FE_NcellnodesPerCellface(c)) &
                normals(1:3,n) = 0.5_pReal &
-                              * math_crossproduct(nodePos(1:3,1+mod(n  ,m)) - nodePos(1:3,n), &
+                              * math_cross(nodePos(1:3,1+mod(n  ,m)) - nodePos(1:3,n), &
                                                    nodePos(1:3,1+mod(n+1,m)) - nodePos(1:3,n))
              normal = 0.5_pReal * sum(normals,2)
              mesh_ipArea(f,i,e) = norm2(normal)

src/prec.f90

@@ -7,11 +7,9 @@
 !> @brief    setting precision for real and int type
 !--------------------------------------------------------------------------------------------------
 module prec
-  use, intrinsic :: IEEE_arithmetic, only:&
+  use, intrinsic :: IEEE_arithmetic
-    IEEE_selected_real_kind
 
   implicit none
-  private
   ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
 #ifdef Abaqus
   integer,     parameter, public :: pReal      = selected_real_kind(15,307)                         !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
@@ -102,7 +100,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
 
-  implicit none
   integer, allocatable, dimension(:) :: realloc_lhs_test
 
   external :: &
@@ -131,7 +128,6 @@ end subroutine prec_init
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq(a,b,tol)
 
-  implicit none
   real(pReal), intent(in)           :: a,b
   real(pReal), intent(in), optional :: tol
   real(pReal)                       :: eps
@@ -155,7 +151,6 @@ end function dEq
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dNeq(a,b,tol)
 
-  implicit none
   real(pReal), intent(in)           :: a,b
   real(pReal), intent(in), optional :: tol
   real(pReal)                       :: eps
@@ -179,7 +174,6 @@ end function dNeq
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq0(a,tol)
 
-  implicit none
   real(pReal), intent(in)           :: a
   real(pReal), intent(in), optional :: tol
   real(pReal)                       :: eps
@@ -203,7 +197,6 @@ end function dEq0
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dNeq0(a,tol)
 
-  implicit none
   real(pReal), intent(in)           :: a
   real(pReal), intent(in), optional :: tol
   real(pReal)                       :: eps
@@ -228,7 +221,6 @@ end function dNeq0
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function cEq(a,b,tol)
 
-  implicit none
   complex(pReal), intent(in)           :: a,b
   real(pReal),    intent(in), optional :: tol
   real(pReal)                          :: eps
@@ -253,7 +245,6 @@ end function cEq
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function cNeq(a,b,tol)
 
-  implicit none
   complex(pReal), intent(in)           :: a,b
   real(pReal),    intent(in), optional :: tol
   real(pReal)                          :: eps

src/source_damage_anisoBrittle.f90

@@ -91,8 +91,6 @@ subroutine source_damage_anisoBrittle_init
    lattice_SchmidMatrix_cleavage, &
    lattice_maxNcleavageFamily
 
- implicit none
-
  integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
  integer(pInt) :: NofMyPhase,p   ,i
  integer(pInt),          dimension(0), parameter :: emptyIntArray    = [integer(pInt)::]
@@ -219,7 +217,6 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
    lattice_maxNcleavageFamily, &
    lattice_NcleavageSystem
 
- implicit none
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -279,7 +276,6 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent
@@ -307,7 +303,6 @@ function source_damage_anisoBrittle_postResults(phase, constituent)
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent

src/source_damage_anisoDuctile.f90

@@ -82,7 +82,6 @@ subroutine source_damage_anisoDuctile_init
    config_phase
 
    
- implicit none
  integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
  integer(pInt) :: NofMyPhase,p ,i
 
@@ -194,7 +193,6 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
    damage, &
    damageMapping
 
- implicit none
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -231,7 +229,6 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent
@@ -259,7 +256,6 @@ function source_damage_anisoDuctile_postResults(phase, constituent)
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent

src/source_damage_isoBrittle.f90

@@ -74,7 +74,6 @@ subroutine source_damage_isoBrittle_init
    config_phase, &
    material_Nphase
 
- implicit none
 
  integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
  integer(pInt) :: NofMyPhase,p,i   
@@ -176,7 +175,6 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
    math_sym33to6, &
    math_I3
 
- implicit none
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -221,7 +219,6 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent
@@ -251,7 +248,6 @@ function source_damage_isoBrittle_postResults(phase, constituent)
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent

src/source_damage_isoDuctile.f90

@@ -74,7 +74,6 @@ subroutine source_damage_isoDuctile_init
    config_phase, &
    material_Nphase
 
- implicit none
 
  integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
  integer(pInt) :: NofMyPhase,p,i
@@ -177,7 +176,6 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
    damage, &
    damageMapping
 
- implicit none
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
@@ -206,7 +204,6 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent
@@ -234,7 +231,6 @@ function source_damage_isoDuctile_postResults(phase, constituent)
  use material, only: &
    sourceState
 
- implicit none
  integer(pInt), intent(in) :: &
    phase, &
    constituent

src/source_thermal_dissipation.f90

@@ -56,7 +56,6 @@ subroutine source_thermal_dissipation_init
     config_phase, &
     material_Nphase
  
-  implicit none
   integer :: Ninstance,instance,source,sourceOffset
   integer :: NofMyPhase,p   
  
@@ -103,7 +102,6 @@ end subroutine source_thermal_dissipation_init
 !--------------------------------------------------------------------------------------------------
 subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase)
 
-  implicit none
   integer, intent(in) :: &
     phase
   real(pReal),  intent(in), dimension(3,3) :: &

src/source_thermal_externalheat.f90

@@ -63,7 +63,6 @@ subroutine source_thermal_externalheat_init
     config_phase, &
     material_Nphase
  
-  implicit none
  
   integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p   
  
@@ -120,7 +119,6 @@ subroutine source_thermal_externalheat_dotState(phase, of)
   use material, only: &
     sourceState
  
-  implicit none
   integer, intent(in) :: &
     phase, &
     of
@@ -140,7 +138,6 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas
   use material, only: &
     sourceState
  
-  implicit none
   integer, intent(in) :: &
     phase, &
     of

src/system_routines.f90

@@ -3,13 +3,9 @@
 !> @brief    provides wrappers to C routines
 !--------------------------------------------------------------------------------------------------
 module system_routines
-  use, intrinsic :: ISO_C_Binding, only: &
+  use, intrinsic :: ISO_C_Binding
-    C_INT, &
-    C_CHAR, &
-    C_NULL_CHAR
  
   implicit none
-  private
   
   public :: &
     signalterm_C, &
@@ -81,16 +77,15 @@ contains
 !--------------------------------------------------------------------------------------------------
 logical function isDirectory(path)
 
- implicit none
+  character(len=*), intent(in) :: path
- character(len=*), intent(in) :: path
+  character(kind=C_CHAR), dimension(1024) :: strFixedLength                                         ! C string as array
- character(kind=C_CHAR), dimension(1024) :: strFixedLength                                         ! C string as array
+  integer :: i 
- integer :: i 
+ 
-
+  strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
- strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
+  do i=1,len(path)                                                                                  ! copy array components
- do i=1,len(path)                                                                                   ! copy array components
+    strFixedLength(i)=path(i:i)
-   strFixedLength(i)=path(i:i)
+  enddo
- enddo
+  isDirectory=merge(.True.,.False.,isDirectory_C(strFixedLength) /= 0_C_INT)
- isDirectory=merge(.True.,.False.,isDirectory_C(strFixedLength) /= 0_C_INT)
 
 end function isDirectory
 
@@ -100,24 +95,23 @@ end function isDirectory
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getCWD()
 
- implicit none
+  character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
- character(kind=C_CHAR), dimension(1024) :: charArray                                               ! C string is an array
+  integer(C_INT) :: stat
- integer(C_INT) :: stat
+  integer :: i
- integer :: i
+ 
-
+  call getCurrentWorkDir_C(charArray,stat)
- call getCurrentWorkDir_C(charArray,stat)
+  if (stat /= 0_C_INT) then
- if (stat /= 0_C_INT) then
+    getCWD = 'Error occured when getting currend working directory'
-   getCWD = 'Error occured when getting currend working directory'
+  else
- else
+    getCWD = repeat('',len(getCWD))
-   getCWD = repeat('',len(getCWD))
+    arrayToString: do i=1,len(getCWD)
-   arrayToString: do i=1,len(getCWD)
+      if (charArray(i) /= C_NULL_CHAR) then
-     if (charArray(i) /= C_NULL_CHAR) then
+        getCWD(i:i)=charArray(i)
-       getCWD(i:i)=charArray(i)
+      else
-     else
+        exit
-       exit
+      endif
-     endif
+    enddo arrayToString
-   enddo arrayToString
+  endif
- endif
 
 end function getCWD
 
@@ -126,24 +120,24 @@ end function getCWD
 !> @brief gets the current host name
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getHostName()
- implicit none
- character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
- integer(C_INT) :: stat
- integer :: i
 
- call getHostName_C(charArray,stat)
+  character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
- if (stat /= 0_C_INT) then
+  integer(C_INT) :: stat
-   getHostName = 'Error occured when getting host name'
+  integer :: i
- else
+ 
-   getHostName = repeat('',len(getHostName))
+  call getHostName_C(charArray,stat)
-   arrayToString: do i=1,len(getHostName)
+  if (stat /= 0_C_INT) then
-     if (charArray(i) /= C_NULL_CHAR) then
+    getHostName = 'Error occured when getting host name'
-       getHostName(i:i)=charArray(i)
+  else
-     else
+    getHostName = repeat('',len(getHostName))
-       exit
+    arrayToString: do i=1,len(getHostName)
-     endif
+      if (charArray(i) /= C_NULL_CHAR) then
-   enddo arrayToString
+        getHostName(i:i)=charArray(i)
- endif
+      else
+        exit
+      endif
+    enddo arrayToString
+  endif
 
 end function getHostName
 
@@ -152,16 +146,16 @@ end function getHostName
 !> @brief changes the current working directory
 !--------------------------------------------------------------------------------------------------
 logical function setCWD(path)
- implicit none
- character(len=*), intent(in) :: path
- character(kind=C_CHAR), dimension(1024) :: strFixedLength                                          ! C string is an array
- integer :: i
 
- strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
+  character(len=*), intent(in) :: path
- do i=1,len(path)                                                                                   ! copy array components
+  character(kind=C_CHAR), dimension(1024) :: strFixedLength                                         ! C string is an array
-   strFixedLength(i)=path(i:i)
+  integer :: i
- enddo
+ 
- setCWD=merge(.True.,.False.,chdir_C(strFixedLength) /= 0_C_INT)
+  strFixedLength = repeat(C_NULL_CHAR,len(strFixedLength))
+  do i=1,len(path)                                                                                  ! copy array components
+    strFixedLength(i)=path(i:i)
+  enddo
+  setCWD=merge(.True.,.False.,chdir_C(strFixedLength) /= 0_C_INT)
 
 end function setCWD
 

src/thermal_adiabatic.f90

@@ -57,7 +57,6 @@ subroutine thermal_adiabatic_init
   use config, only: &
     config_homogenization
  
-  implicit none
   integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog   
   character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
   character(len=65536), dimension(:), allocatable :: outputs
@@ -124,7 +123,6 @@ function thermal_adiabatic_updateState(subdt, ip, el)
     temperatureRate, &
     thermalMapping
   
-  implicit none
   integer,     intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -181,7 +179,6 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
     crystallite_S, &
     crystallite_Lp  
  
-  implicit none
   integer,     intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -246,7 +243,6 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
   use mesh, only: &
     mesh_element
  
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -282,7 +278,6 @@ function thermal_adiabatic_getMassDensity(ip,el)
   use mesh, only: &
     mesh_element
     
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -312,7 +307,6 @@ function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
   use material, only: &
     temperature
  
-  implicit none
   integer, intent(in) :: &
     homog, &
     instance, &

src/thermal_conduction.f90

@@ -58,7 +58,6 @@ subroutine thermal_conduction_init
   use config, only: &
     config_homogenization
  
-  implicit none
   integer :: maxNinstance,section,instance,i
   integer :: sizeState
   integer :: NofMyHomog   
@@ -135,7 +134,6 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
     crystallite_S, &
     crystallite_Lp  
  
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -205,7 +203,6 @@ function thermal_conduction_getConductivity33(ip,el)
   use crystallite, only: &
     crystallite_push33ToRef
  
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -239,7 +236,6 @@ function thermal_conduction_getSpecificHeat(ip,el)
   use mesh, only: &
     mesh_element
  
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -273,7 +269,6 @@ function thermal_conduction_getMassDensity(ip,el)
   use mesh, only: &
     mesh_element
  
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -306,7 +301,6 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el)
     temperatureRate, &
     thermalMapping
 
-  implicit none
   integer, intent(in) :: &
     ip, &                                                                                           !< integration point number
     el                                                                                              !< element number
@@ -332,7 +326,6 @@ function thermal_conduction_postResults(homog,instance,of) result(postResults)
   use material, only: &
     temperature
  
-  implicit none
   integer,              intent(in) :: &
     homog, &
     instance, &

src/thermal_isothermal.f90

@@ -22,7 +22,6 @@ subroutine thermal_isothermal_init()
    material_Nhomogenization
  use material
  
- implicit none
  integer :: &
    homog, &
    NofMyHomog