general speed up

spectral format records physical coordinates
tessellation now based on phys coords instead of discretization.

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,\
+parser.add_option("--homogenization", type="int",\
-        dest="N",\
+        dest="homogenization",\
-        help="number of subdivisions")
+        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)
+  print('no valid Mentat release found in %s'%MSCpath)
-  sys.exit(-1)
+  if options.port != None: sys.exit(-1)
-
-(options, args) = parser.parse_args()
 
 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 = {}

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

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) a, b, c, N_Seeds, theGrain, i, j, k, l, m
- integer(pInt), dimension(3) :: coordinates
  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)             
+  posGeom = IO_stringPos(line,6)                         ! split line
-  Seeds(i,1)=IO_floatValue(line,posGeom,1)
+  do j=1,3
-  Seeds(i,2)=IO_floatValue(line,posGeom,2)
+    seeds(i,j) =      IO_floatValue(line,posGeom,j)
-  Seeds(i,3)=IO_floatValue(line,posGeom,3)
+    grainEuler(i,j) = IO_floatValue(line,posGeom,j+3)
-  grainEuler(i,1)=IO_floatValue(line,posGeom,4)
+  enddo
-  grainEuler(i,2)=IO_floatValue(line,posGeom,5)
-  grainEuler(i,3)=IO_floatValue(line,posGeom,6)
  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 + &
+            theDist = ( dim(1) * ( delta(1)-real(k,pReal) ) )**2 + &
-                               (mod((i/a), b)+1 -seeds(j,2)-l*b)**2 + &
+                      ( dim(2) * ( delta(2)-real(l,pReal) ) )**2 + &
-                           (mod((i/(a*b)), c)+1 -seeds(j,3)-k*c)**2)
+                      ( dim(3) * ( delta(3)-real(m,pReal) ) )**2
-            if (myDistance < minDistance) then
+            if (theDist < minDist) then
-              minDistance = myDistance
+              minDist = theDist
               theGrain = j
-            end if
+            endif
-          end do
+          enddo
-        end do
+        enddo
-      end do
+      enddo
-    end do
+    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.'