diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7b33574de..7d8f5ad33 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -498,7 +498,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 diff --git a/processing/post/marc_to_vtk.py b/processing/post/marc_to_vtk.py deleted file mode 100755 index 05f7a6908..000000000 --- a/processing/post/marc_to_vtk.py +++ /dev/null @@ -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() diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 45e423fe3..3a0954e29 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -298,6 +298,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: & @@ -311,6 +313,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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 06da6ab2e..d4d30d5da 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -66,7 +66,8 @@ module homogenization public :: & homogenization_init, & materialpoint_stressAndItsTangent, & - materialpoint_postResults + materialpoint_postResults, & + homogenization_results private :: & partitionDeformation, & updateState, & @@ -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) @@ -969,4 +969,39 @@ 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 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)) + + end select + + enddo +#endif + + +end subroutine homogenization_results + end module homogenization diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 012ad6086..57b573506 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -1109,6 +1109,36 @@ pure function homogenization_RGC_postResults(instance,of) result(postResults) end function homogenization_RGC_postResults +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +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)) + + + 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 !--------------------------------------------------------------------------------------------------