From ee7a8ad52a2480151d67baa1612c06e57c225c2c Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 24 Jan 2011 16:21:17 +0000 Subject: [PATCH] general speed up spectral format records physical coordinates tessellation now based on phys coords instead of discretization. --- processing/pre/mentat_spectralBox | 131 +++++++++++++++++++---- processing/pre/voronoi_randomSeeding.f90 | 5 +- processing/pre/voronoi_tessellation.f90 | 69 ++++++------ 3 files changed, 149 insertions(+), 56 deletions(-) diff --git a/processing/pre/mentat_spectralBox b/processing/pre/mentat_spectralBox index 895a85dfd..a23fbee87 100755 --- a/processing/pre/mentat_spectralBox +++ b/processing/pre/mentat_spectralBox @@ -45,7 +45,9 @@ def output(cmds,locals,dest): +#-------------------- def init(): +#-------------------- return ["*new_model yes", "*reset", "*select_clear", @@ -59,7 +61,9 @@ def init(): ] -def mesh(N,d): +#-------------------- +def mesh(r,d): +#-------------------- return [ "*add_nodes", "%f %f %f"%(0.0,0.0,0.0), @@ -73,11 +77,11 @@ def mesh(N,d): "*add_elements", range(1,9), "*sub_divisions", - "%i %i %i"%(N[2],N[1],N[0]), + "%i %i %i"%(r[2],r[1],r[0]), "*subdivide_elements", "all_existing", "*set_sweep_tolerance", - "%f"%(float(min(d))/max(N)/2.0), + "%f"%(float(min(d))/max(r)/2.0), "*sweep_all", "*renumber_all", "*set_move_scale_factor x -1", @@ -89,7 +93,9 @@ def mesh(N,d): ] +#-------------------- def material(): +#-------------------- cmds = [\ "*new_mater standard", "*mater_option general:state:solid", @@ -107,7 +113,9 @@ def material(): return cmds +#-------------------- def geometry(): +#-------------------- cmds = [\ "*geometry_type mech_three_solid", "*geometry_option red_integ_capacity:on", @@ -120,13 +128,13 @@ def geometry(): return cmds -def initial_conditions(N,data): +#-------------------- +def initial_conditions(homogenization,grains): +#-------------------- elements = [] element = 0 - for line in data: + for id in grains: element += 1 - phi1,phi,phi2,x,y,z,id,phase = line.split() - id = int(id) if len(elements) < id: for i in range(id-len(elements)): elements.append([]) @@ -134,17 +142,17 @@ def initial_conditions(N,data): cmds = [\ "*new_icond", - "*icond_name temperature", + "*icond_name _temperature", "*icond_type state_variable", "*icond_param_value state_var_id 1", "*icond_dof_value var 300", "*add_icond_elements", "all_existing", "*new_icond", - "*icond_name homogenization", + "*icond_name _homogenization", "*icond_type state_variable", "*icond_param_value state_var_id 2", - "*icond_dof_value var 1", + "*icond_dof_value var %i"%homogenization, "*add_icond_elements", "all_existing", ] @@ -163,10 +171,67 @@ def initial_conditions(N,data): return cmds +#-------------------- +def parse_geomFile(content,homog): +#-------------------- + (skip,key) = content[0].split()[:2] + if key[:4].lower() == 'head': + skip = int(skip)+1 + else: + skip = 0 + + res = [0,0,0] + dim = [0.0,0.0,0.0] + homog = 0 + + for line in content[:skip]: + data = line.split() + if data[0].lower() == 'resolution': + res = map(int,data[2:8:2]) + if data[0].lower() == 'dimension': + dim = map(float,data[2:8:2]) + if data[0].lower() == 'homogenization': + homog = int(data[1]) + + grains = [] + for line in content[skip:]: + grains.append(int(line.split()[0])) + + return (res,dim,homog,grains) + +#-------------------- +def parse_spectralFile(content,homog): +#-------------------- + + coords = [{},{},{}] + maxBox = [-1.0e20,-1.0e20,-1.0e20] + minBox = [ 1.0e20, 1.0e20, 1.0e20] + dim = [0.0,0.0,0.0] + res = [0,0,0] + grains = [] + + for line in content: + data = line.split()[3:7] + grains.append(int(data[3])) + for i in range(3): + maxBox[i] = max(maxBox[i],float(data[i])) + minBox[i] = min(minBox[i],float(data[i])) + coords[i][data[i]] = True + + for i in range(3): + res[i] = len(coords[i]) + dim[i] = (maxBox[i]-minBox[i])*res[i]/(res[i]-1.0) + + return (res,dim,homog,grains) + # ----------------------- MAIN ------------------------------- parser = OptionParser(usage='%prog [options] spectral.datafile', description = """ -Generate FE hexahedral mesh from spectral description file of format: phi1,Phi,phi2,x,y,z,id,phase. +Generate FE hexahedral mesh from spectral description file. + +Acceptable formats are +geom: header plus list of grain numbers or +spectral: phi1,Phi,phi2,x,y,z,id,phase. $Id$ """) @@ -176,12 +241,20 @@ parser.add_option("-p", "--port", type="int",\ parser.add_option("-d", "--dimension", type="int", nargs=3,\ dest="d",\ help="physical dimension") -parser.add_option("-N", "--subdivisions", type="int", nargs=3,\ - dest="N",\ - help="number of subdivisions") +parser.add_option("--homogenization", type="int",\ + dest="homogenization",\ + help="index of homogenization to be chosen from material.config file") +parser.add_option("-s", "--spectral", action="store_const", const="spectral",\ + dest="filetype",\ + help="file has 'spectral' format") +parser.add_option("-g", "--geom", action="store_const", const="geom",\ + dest="filetype",\ + help="file has 'geom' format") + +parser.set_defaults(homogenization = 1) + +(options, args) = parser.parse_args() -parser.set_defaults(d = (16,16,16)) -parser.set_defaults(N = (16,16,16)) try: file = open('%s/../MSCpath'%os.path.dirname(os.path.realpath(sys.argv[0]))) @@ -203,10 +276,8 @@ for release,subdirs in sorted(releases.items(),reverse=True): try: from py_mentat import * except: - print('error: no valid Mentat release found in %s'%MSCpath) - sys.exit(-1) - -(options, args) = parser.parse_args() + print('no valid Mentat release found in %s'%MSCpath) + if options.port != None: sys.exit(-1) if not os.path.isfile(args[0]): parser.error("cannot open %s"%args[0]) @@ -215,12 +286,28 @@ file = open(args[0]) content = file.readlines() file.close() +if options.filetype not in ['spectral','geom']: + options.filetype = os.path.splitext(args[0])[1][1:] + +print '\nparsing %s...'%options.filetype, +sys.stdout.flush() + +(res,dim,homog,grains) = {\ + 'geom': parse_geomFile, + 'spectral': parse_spectralFile, + }[options.filetype](content,options.homogenization) + +print '%i grains in %s with resolution %s and homogenization %i\n'%(len(list(set(grains))),str(dim),str(res),homog) + + cmds = [\ init(), - mesh(options.N,options.d), + mesh(res,dim), material(), geometry(), - initial_conditions(options.N,content), + initial_conditions(homog,grains), + '*identify_sets', + '*redraw', ] outputLocals = {} diff --git a/processing/pre/voronoi_randomSeeding.f90 b/processing/pre/voronoi_randomSeeding.f90 index 4a0e94a29..015edae0c 100644 --- a/processing/pre/voronoi_randomSeeding.f90 +++ b/processing/pre/voronoi_randomSeeding.f90 @@ -27,6 +27,7 @@ program voronoi print*, '******************************************************************************' print*, ' Voronoi description file' print*, '******************************************************************************' + print*, '$Id$' print*, '' print*, 'generates:' print*, ' * description file "_OUTPUT_.seeds":' @@ -35,11 +36,11 @@ program voronoi read(*, *), filename write(*, '(A)', advance = 'NO') 'Please enter value random seed: ' read(*, *), randomSeed; randomSeed = max(0_pInt,randomSeed) - write(*, '(A)', advance = 'NO') 'Please enter value for first resolution: ' + write(*, '(A)', advance = 'NO') 'Please enter value for first resolution: ' read(*, *), a write(*, '(A)', advance = 'NO') 'Please enter value for second resolution: ' read(*, *), b - write(*, '(A)', advance = 'NO') 'Please enter value for third resolution: ' + write(*, '(A)', advance = 'NO') 'Please enter value for third resolution: ' read(*, *), c write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: ' read(*, *), N_Seeds diff --git a/processing/pre/voronoi_tessellation.f90 b/processing/pre/voronoi_tessellation.f90 index 68662b0d7..3b3329ca2 100644 --- a/processing/pre/voronoi_tessellation.f90 +++ b/processing/pre/voronoi_tessellation.f90 @@ -172,16 +172,17 @@ program voronoi logical gotN_Seeds, gotResolution logical, dimension(:), allocatable :: grainCheck character(len=1024) input_name, output_name, format1, format2, N_Digits, line, key - integer(pInt) a, b, c, N_Seeds, seedPoint, theGrain, minDistance, myDistance, i, j, k, l, m - integer(pInt), dimension(3) :: coordinates + integer(pInt) a, b, c, N_Seeds, theGrain, i, j, k, l, m integer(pInt), dimension (1+2*7) :: posGeom real(pReal), dimension(:,:), allocatable :: grainEuler, seeds - real(pReal), dimension(3) :: dim + real(pReal), dimension(3) :: step,dim,delta + real(pReal) minDist, theDist real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal print*, '******************************************************************************' print*, ' Spectral Method Problem Set-up' print*, '******************************************************************************' + print*, '$Id$' print*, '' print*, 'generates:' print*, ' * geom file "_OUTPUT_.geom": Geometrical information for solver' @@ -220,23 +221,27 @@ program voronoi enddo 100 allocate(grainEuler(N_Seeds,3)) - allocate(Seeds(N_Seeds,3)) + allocate(seeds(N_Seeds,3)) allocate(grainCheck(N_Seeds)) grainCheck = .false. print*, 'resolution: ' ,a,b,c - write(*, '(A)', advance = 'NO') 'New first resolution: ' + write(*, '(A)', advance = 'NO') 'New first resolution: ' read(*, *), a write(*, '(A)', advance = 'NO') 'New second resolution: ' read(*, *), b - write(*, '(A)', advance = 'NO') 'New third resolution: ' + write(*, '(A)', advance = 'NO') 'New third resolution: ' read(*, *), c - write(*, '(A)', advance = 'NO') 'First dimension: ' + step(1) = 1.0_pReal/real(a,pReal) + step(2) = 1.0_pReal/real(b,pReal) + step(3) = 1.0_pReal/real(c,pReal) + + write(*, '(A)', advance = 'NO') 'First dimension: ' read(*, *), dim(1) write(*, '(A)', advance = 'NO') 'Second dimension: ' read(*, *), dim(2) - write(*, '(A)', advance = 'NO') 'Third dimension: ' + write(*, '(A)', advance = 'NO') 'Third dimension: ' read(*, *), dim(3) rewind(20) @@ -251,19 +256,14 @@ program voronoi do i=1, N_seeds read(20,'(a1024)') line if (IO_isBlank(line)) cycle ! skip empty lines - posGeom = IO_stringPos(line,6) - Seeds(i,1)=IO_floatValue(line,posGeom,1) - Seeds(i,2)=IO_floatValue(line,posGeom,2) - Seeds(i,3)=IO_floatValue(line,posGeom,3) - grainEuler(i,1)=IO_floatValue(line,posGeom,4) - grainEuler(i,2)=IO_floatValue(line,posGeom,5) - grainEuler(i,3)=IO_floatValue(line,posGeom,6) + posGeom = IO_stringPos(line,6) ! split line + do j=1,3 + seeds(i,j) = IO_floatValue(line,posGeom,j) + grainEuler(i,j) = IO_floatValue(line,posGeom,j+3) + enddo enddo close(20) - seeds(:,1) = seeds(:,1)*real(a, pReal) - seeds(:,2) = seeds(:,2)*real(b, pReal) - seeds(:,3) = seeds(:,3)*real(c, pReal) ! calculate No. of digits needed for name of the grains i = 1 + int( log10(real( N_Seeds ))) @@ -302,33 +302,38 @@ program voronoi write(20, '(A)'), 'homogenization 1' format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' ! geom format - format2 = '(3(tr2, f6.2), 3(I10), I10, a)' ! spectral (Lebensohn) format + format2 = '(3(tr2, f6.2), 3(tr2,g10.5), I10, a)' ! spectral (Lebensohn) format ! perform voronoi tessellation and write result to files do i = 0, a*b*c-1 - minDistance = a*a+b*b+c*c + minDist = dim(1)*dim(1)+dim(2)*dim(2)+dim(3)*dim(3) ! diagonal of rve do j = 1, N_Seeds + delta(1) = step(1)*(mod(i , a)+0.5_pReal) - seeds(j,1) + delta(2) = step(2)*(mod(i/a , b)+0.5_pReal) - seeds(j,2) + delta(3) = step(3)*(mod(i/a/b, c)+0.5_pReal) - seeds(j,3) do k = -1, 1 ! left, me, right image do l = -1, 1 ! front, me, back image do m = -1, 1 ! lower, me, upper image - myDistance = (( mod(i, a)+1 -seeds(j,1)-m*a)**2 + & - (mod((i/a), b)+1 -seeds(j,2)-l*b)**2 + & - (mod((i/(a*b)), c)+1 -seeds(j,3)-k*c)**2) - if (myDistance < minDistance) then - minDistance = myDistance + theDist = ( dim(1) * ( delta(1)-real(k,pReal) ) )**2 + & + ( dim(2) * ( delta(2)-real(l,pReal) ) )**2 + & + ( dim(3) * ( delta(3)-real(m,pReal) ) )**2 + if (theDist < minDist) then + minDist = theDist theGrain = j - end if - end do - end do - end do - end do + endif + enddo + enddo + enddo + enddo grainCheck(theGrain) = .true. write(20, trim(format1)), theGrain write(21, trim(format2)), grainEuler(theGrain,1), grainEuler(theGrain,2), grainEuler(theGrain,3), & - mod(i, a)+1, mod((i/a), b)+1, mod((i/(a*b)), c)+1, & + dim(1)*step(1)*(mod(i , a)+0.5_pReal), & + dim(2)*step(2)*(mod(i/a , b)+0.5_pReal), & + dim(3)*step(3)*(mod(i/a/b, c)+0.5_pReal), & theGrain, ' 1' - end do + enddo close(20) close(21) print*, 'voronoi tesselation done.'