From 650b7ef4acda9aa2cc9c7913b6114a87a1789d75 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Tue, 23 Apr 2013 18:30:56 +0000 Subject: [PATCH] added python script that produces vtk files with the (node based and ip based) deformed mesh from marc output file --- code/DAMASK_marc.f90 | 3 +- code/damask.core.pyf | 12 +- code/mesh.f90 | 29 +++-- processing/post/marc_deformedGeometry.py | 152 +++++++++++++++++++++++ processing/setup/symLink_Processing.py | 1 + 5 files changed, 185 insertions(+), 12 deletions(-) create mode 100755 processing/post/marc_deformedGeometry.py diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index 92f7494a7..51b79c2eb 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -203,6 +203,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & mesh_element, & mesh_node0, & mesh_node, & + mesh_Ncellnodes, & mesh_cellnode, & mesh_build_cellnodes, & mesh_build_ipCoordinates, & @@ -352,7 +353,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & call debug_reset() ! resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1_pInt - mesh_cellnode = mesh_build_cellnodes(mesh_node) ! update cell node coordinates + mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) ! update cell node coordinates call mesh_build_ipCoordinates() ! update ip coordinates endif if ( outdatedByNewInc ) then diff --git a/code/damask.core.pyf b/code/damask.core.pyf index fa007414d..32b98b6dd 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -193,12 +193,16 @@ python module core ! in character(len=*), intent(in) :: filepath end function mesh_init_postprocessing - function mesh_build_cellnodes(nodes) ! in :mesh:mesh.f90 - integer - real*8, dimension(3,:), intent(in) :: nodes - real*8, dimension(3,:) :: mesh_build_cellnodes + function mesh_build_cellnodes(nodes,Ncellnodes) ! in :mesh:mesh.f90 + integer, intent(in) :: Ncellnodes + real*8, dimension(3,:), intent(in) :: nodes + real*8, dimension(3,Ncellnodes), depend(Ncellnodes) :: mesh_build_cellnodes end function mesh_build_cellnodes + function mesh_get_Ncellnodes() ! in :mesh:mesh.f90 + integer :: mesh_get_Ncellnodes + end function mesh_get_Ncellnodes + end module mesh end interface end python module core diff --git a/code/mesh.f90 b/code/mesh.f90 index d572826b8..097fcfc70 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -419,7 +419,8 @@ module mesh mesh_build_ipVolumes, & mesh_build_ipCoordinates, & mesh_cellCenterCoordinates, & - mesh_init_postprocessing + mesh_init_postprocessing, & + mesh_get_Ncellnodes #ifdef Spectral public :: & mesh_regrid, & @@ -584,7 +585,7 @@ subroutine mesh_init(ip,el) call mesh_get_damaskOptions(fileUnit) close (fileUnit) call mesh_build_cellconnectivity - mesh_cellnode = mesh_build_cellnodes(mesh_node) + mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates call mesh_build_ipVolumes call mesh_build_ipAreas @@ -605,7 +606,7 @@ subroutine mesh_init(ip,el) call mesh_get_damaskOptions(fileUnit) close (fileUnit) call mesh_build_cellconnectivity - mesh_cellnode = mesh_build_cellnodes(mesh_node) + mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates call mesh_build_ipVolumes call mesh_build_ipAreas @@ -630,7 +631,7 @@ subroutine mesh_init(ip,el) call mesh_get_damaskOptions(fileUnit) close (fileUnit) call mesh_build_cellconnectivity - mesh_cellnode = mesh_build_cellnodes(mesh_node) + mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes) call mesh_build_ipCoordinates call mesh_build_ipVolumes call mesh_build_ipAreas @@ -790,12 +791,13 @@ end subroutine mesh_build_cellconnectivity !> Build list of cellnodes' coordinates. !> Cellnode coordinates are calculated from a weighted sum of node coordinates. !-------------------------------------------------------------------------------------------------- -function mesh_build_cellnodes(nodes) +function mesh_build_cellnodes(nodes,Ncellnodes) implicit none + integer(pInt), intent(in) :: Ncellnodes !< requested number of cellnodes real(pReal), dimension(3,mesh_Nnodes), intent(in) :: nodes - real(pReal), dimension(3,mesh_Ncellnodes) :: mesh_build_cellnodes + real(pReal), dimension(3,Ncellnodes) :: mesh_build_cellnodes integer(pInt) & e,t,n,m, & @@ -805,7 +807,7 @@ function mesh_build_cellnodes(nodes) mesh_build_cellnodes = 0.0_pReal !$OMP PARALLEL DO PRIVATE(e,localCellnodeID,t,myCoords) - do n = 1_pInt,mesh_Ncellnodes ! loop over cell nodes + do n = 1_pInt,Ncellnodes ! loop over cell nodes e = mesh_cellnodeParent(1,n) localCellnodeID = mesh_cellnodeParent(2,n) t = mesh_element(2,e) ! get element type @@ -5300,5 +5302,18 @@ integer function mesh_init_postprocessing(filepath) end function mesh_init_postprocessing + +!-------------------------------------------------------------------------------------------------- +!> @brief just returns global variable mesh_Ncellnodes +!-------------------------------------------------------------------------------------------------- +integer(pInt) function mesh_get_Ncellnodes() + + implicit none + + mesh_get_Ncellnodes = mesh_Ncellnodes + +end function mesh_get_Ncellnodes + + end module mesh diff --git a/processing/post/marc_deformedGeometry.py b/processing/post/marc_deformedGeometry.py new file mode 100755 index 000000000..1528c8dfb --- /dev/null +++ b/processing/post/marc_deformedGeometry.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python + +import os, sys, math, string, numpy, shutil +import damask +from optparse import OptionParser, Option + + + +# ----------------------------- +class MyOption(Option): +# ----------------------------- +# used for definition of new option parser action 'extend', which enables to take multiple option arguments +# taken from online tutorial http://docs.python.org/library/optparse.html + + ACTIONS = Option.ACTIONS + ("extend",) + STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) + TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) + ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) + + def take_action(self, action, dest, opt, value, values, parser): + if action == "extend": + lvalue = value.split(",") + values.ensure_value(dest, []).extend(lvalue) + else: + Option.take_action(self, action, dest, opt, value, values, parser) + + +# ----------------------------- +# MAIN FUNCTION STARTS HERE +# ----------------------------- + +# --- input parsing + +parser = OptionParser(option_class=MyOption, usage='%prog [options] resultfile', description = """ +Extract data from a .t16 (MSC.Marc) results file. +""" + string.replace('$Id: postResults.py 2107 2013-01-28 16:25:43Z MPIE\p.eisenlohr $','\n','\\n') +) + +parser.add_option('-d','--dir', dest='dir', \ + help='name of subdirectory to hold output [%default]') +parser.add_option('-r','--range', dest='range', type='int', nargs=3, \ + help='range of positions (or increments) to output (start, end, step) [all]') +parser.add_option('--increments', action='store_true', dest='getIncrements', \ + help='switch to increment range [%default]') + + +parser.set_defaults(dir = 'vtk') +parser.set_defaults(getIncrements= False) + +(options, files) = parser.parse_args() + +# --- basic sanity checks + +if files == []: + parser.print_help() + parser.error('no file specified...') + +filename = os.path.splitext(files[0])[0] +if not os.path.exists(filename+'.t16'): + parser.print_help() + parser.error('invalid file "%s" specified...'%filename+'.t16') + + +# --- more sanity checks + +sys.path.append(damask.solver.Marc().libraryPath('../../')) +try: + from py_post import * +except: + print('error: no valid Mentat release found') + sys.exit(-1) + + +# --------------------------- open results file and initialize mesh ---------- + +p = post_open(filename+'.t16') +p.moveto(0) +Nnodes = p.nodes() +Nincrements = p.increments() - 1 # t16 contains one "virtual" increment (at 0) +if damask.core.mesh.mesh_init_postprocessing(filename+'.mesh') > 0: + print('error: init not successful') + sys.exit(-1) +Ncellnodes = damask.core.mesh.mesh_get_Ncellnodes() + + +# --------------------------- create output dir -------------------------------- + +dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) +if not os.path.isdir(dirname): + os.mkdir(dirname,0755) + + +# --------------------------- get positions -------------------------------- + +incAtPosition = {} +positionOfInc = {} + +for position in range(Nincrements): + p.moveto(position+1) + incAtPosition[position] = p.increment # remember "real" increment at this position + positionOfInc[p.increment] = position # remember position of "real" increment + +if not options.range: + options.getIncrements = False + locations = range(Nincrements) # process all positions +else: + options.range = list(options.range) # convert to list + if options.getIncrements: + locations = [positionOfInc[x] for x in range(options.range[0],options.range[1]+1,options.range[2]) + if x in positionOfInc] + else: + locations = range( max(0,options.range[0]), + min(Nincrements,options.range[1]+1), + options.range[2] ) + +increments = [incAtPosition[x] for x in locations] # build list of increments to process + + + +# --------------------------- loop over positions -------------------------------- + +for incCount,position in enumerate(locations): # walk through locations + + p.moveto(position+1) # wind to correct position + + # --- get displacements + + node_displacement = [[0,0,0] for i in range(Nnodes)] + for n in range(Nnodes): + if p.node_displacements(): + node_displacement[n] = list(p.node_displacement(n)) + c = damask.core.mesh.mesh_build_cellnodes(numpy.array(node_displacement).T,Ncellnodes) + cellnode_displacement = [[c[i][n] for i in range(3)] for n in range(Ncellnodes)] + + + # --- append displacements to corresponding files + + for geomtype in ['nodebased', 'ipbased']: + outFilename = eval('"'+eval("'%%s_%%s_inc%%0%ii.vtk'%(math.log10(max(increments+[1]))+1)")+'"%(dirname + os.sep + os.path.split(filename)[1],geomtype,increments[incCount])') + print outFilename + shutil.copyfile('%s_%s.vtk'%(filename,geomtype),outFilename) + + with open(outFilename,'a') as myfile: + myfile.write("POINT_DATA %i\n"%{'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype]) + myfile.write("VECTORS displacement float\n") + coordinates = {'nodebased':node_displacement,'ipbased':cellnode_displacement}[geomtype] + for n in range({'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype]): + myfile.write("%.8e %.8e %.8e\n"%(coordinates[n][0],coordinates[n][1],coordinates[n][2])) + + + +# --------------------------- DONE -------------------------------- diff --git a/processing/setup/symLink_Processing.py b/processing/setup/symLink_Processing.py index 6503f3937..daaa72f6e 100755 --- a/processing/setup/symLink_Processing.py +++ b/processing/setup/symLink_Processing.py @@ -55,6 +55,7 @@ bin_link = { \ 'deleteColumn.py', 'deleteInfo.py', 'filterTable.py', + 'marc_deformedGeometry.py', 'mentat_colorMap.py', 'nodesFromCentroids.py', 'perceptualUniformColorMap.py',