Merge branch 'HDF5-out-homog-2' into 'development'

Hdf5 out homog 2

See merge request damask/DAMASK!78

.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

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/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,12 +20,6 @@ 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 :: &

src/material.f90

@@ -53,55 +53,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