corrected mixup of Fortran (1,2,3,...) indexing to Python (0,1,2,...) style.
fixed erroneous array indexing in neighborhood calculation.
This commit is contained in:
Philip Eisenlohr
parent
5814262f55
commit
f598e3b6ed
|
|
@ -115,7 +115,7 @@ def vtk_writeASCII_mesh(mesh,data,res):
|
|||
|
||||
cmds = [\
|
||||
'# vtk DataFile Version 3.1',
|
||||
'powered by 3Dvisualize',
|
||||
'powered by $Id$',
|
||||
'ASCII',
|
||||
'DATASET UNSTRUCTURED_GRID',
|
||||
'POINTS %i float'%N1,
|
||||
|
|
@ -124,19 +124,19 @@ def vtk_writeASCII_mesh(mesh,data,res):
|
|||
]
|
||||
|
||||
# cells
|
||||
for i in range (res[2]):
|
||||
for j in range (res[1]):
|
||||
for k in range (res[0]):
|
||||
base = i*(res[1]+1)*(res[2]+1)+j*(res[1]+1)+k
|
||||
for z in range (res[2]):
|
||||
for y in range (res[1]):
|
||||
for x in range (res[0]):
|
||||
base = z*(res[1]+1)*(res[0]+1)+y*(res[0]+1)+x
|
||||
cmds.append('8 '+'\t'.join(map(str,[ \
|
||||
base,
|
||||
base+1,
|
||||
base+res[1]+2,
|
||||
base+res[1]+1,
|
||||
base+(res[1]+1)*(res[2]+1),
|
||||
base+(res[1]+1)*(res[2]+1)+1,
|
||||
base+(res[1]+1)*(res[2]+1)+res[1]+2,
|
||||
base+(res[1]+1)*(res[2]+1)+res[1]+1,
|
||||
base+res[0]+2,
|
||||
base+res[0]+1,
|
||||
base+(res[1]+1)*(res[0]+1),
|
||||
base+(res[1]+1)*(res[0]+1)+1,
|
||||
base+(res[1]+1)*(res[0]+1)+res[0]+2,
|
||||
base+(res[1]+1)*(res[0]+1)+res[0]+1,
|
||||
])))
|
||||
cmds += [\
|
||||
'CELL_TYPES %i'%N,
|
||||
|
|
@ -177,24 +177,26 @@ def gmsh_writeASCII_mesh(mesh,data,res):
|
|||
'%i'%N,
|
||||
]
|
||||
|
||||
# cells
|
||||
n_elem = 0
|
||||
for i in range (res[2]):
|
||||
for j in range (res[1]):
|
||||
for k in range (res[0]):
|
||||
base = i*(res[1]+1)*(res[2]+1)+j*(res[1]+1)+k
|
||||
for z in range (res[2]):
|
||||
for y in range (res[1]):
|
||||
for x in range (res[0]):
|
||||
base = z*(res[1]+1)*(res[0]+1)+y*(res[0]+1)+x
|
||||
n_elem +=1
|
||||
cmds.append('\t'.join(map(str,[ \
|
||||
n_elem,
|
||||
'5',
|
||||
base,
|
||||
base+1,
|
||||
base+res[1]+2,
|
||||
base+res[1]+1,
|
||||
base+(res[1]+1)*(res[2]+1),
|
||||
base+(res[1]+1)*(res[2]+1)+1,
|
||||
base+(res[1]+1)*(res[2]+1)+res[1]+2,
|
||||
base+(res[1]+1)*(res[2]+1)+res[1]+1,
|
||||
base+res[0]+2,
|
||||
base+res[0]+1,
|
||||
base+(res[1]+1)*(res[0]+1),
|
||||
base+(res[1]+1)*(res[0]+1)+1,
|
||||
base+(res[1]+1)*(res[0]+1)+res[0]+2,
|
||||
base+(res[1]+1)*(res[0]+1)+res[0]+1,
|
||||
])))
|
||||
|
||||
cmds += [\
|
||||
'ElementData',
|
||||
'1',
|
||||
|
|
@ -228,7 +230,7 @@ def vtk_writeASCII_points(coordinates,data,res):
|
|||
|
||||
cmds = [\
|
||||
'# vtk DataFile Version 3.1',
|
||||
'powered by 3Dvisualize',
|
||||
'powered by $Id$',
|
||||
'ASCII',
|
||||
'DATASET UNSTRUCTURED_GRID',
|
||||
'POINTS %i float'%N,
|
||||
|
|
@ -267,7 +269,7 @@ def vtk_writeASCII_box(diag,defgrad):
|
|||
|
||||
cmds = [\
|
||||
'# vtk DataFile Version 3.1',
|
||||
'powered by 3Dvisualize',
|
||||
'powered by $Id$',
|
||||
'ASCII',
|
||||
'DATASET UNSTRUCTURED_GRID',
|
||||
'POINTS 8 float',
|
||||
|
|
@ -285,7 +287,7 @@ def vtk_writeASCII_box(diag,defgrad):
|
|||
# ----------------------- MAIN -------------------------------
|
||||
|
||||
parser = OptionParser(option_class=extendedOption, usage='%prog [options] datafile', description = """
|
||||
Produce VTK file from data field.
|
||||
Produce VTK file from data field. Coordinates are taken from (consecutive) ip.x, ip.y, and ip.z columns.
|
||||
|
||||
$Id$
|
||||
""")
|
||||
|
|
@ -361,33 +363,37 @@ for filename in args:
|
|||
maxcol = max(maxcol,col+1)
|
||||
break
|
||||
|
||||
values = numpy.array(sorted([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],
|
||||
key=lambda x:(x[ipcol+0],x[ipcol+1],x[ipcol+2])), # sort with z as fastest and x as slowest index
|
||||
'd')
|
||||
|
||||
values = numpy.array([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],'d')
|
||||
N = len(values)
|
||||
grid = [{},{},{}]
|
||||
for i in range(N):
|
||||
grid[0][str(values[i,ipcol+0])] = True
|
||||
grid[1][str(values[i,ipcol+1])] = True
|
||||
grid[2][str(values[i,ipcol+2])] = True
|
||||
for j in range(3):
|
||||
for i in range(N):
|
||||
grid[j][str(values[i,ipcol+j])] = True
|
||||
|
||||
res = numpy.array([len(grid[0]),\
|
||||
len(grid[1]),\
|
||||
len(grid[2]),],'i')
|
||||
dim = numpy.array([max(map(float,grid[0].keys()))-min(map(float,grid[0].keys())),\
|
||||
max(map(float,grid[1].keys()))-min(map(float,grid[1].keys())),\
|
||||
max(map(float,grid[1].keys()))-min(map(float,grid[1].keys())),]*res/(res-numpy.ones(3)), 'd')
|
||||
max(map(float,grid[2].keys()))-min(map(float,grid[2].keys())),]*res/(res-numpy.ones(3)), 'd')
|
||||
|
||||
print 'resolution',res
|
||||
print 'dimension',dim
|
||||
defgrad_av = postprocessingMath.tensor_avg(res[0],res[1],res[2],\
|
||||
numpy.reshape(values[:,column['tensor'][options.defgrad]:
|
||||
column['tensor'][options.defgrad]+9],
|
||||
(res[0],res[1],res[2],3,3)))
|
||||
numpy.reshape(values[:,column['tensor'][options.defgrad]:
|
||||
column['tensor'][options.defgrad]+9],
|
||||
(res[0],res[1],res[2],3,3)))
|
||||
print 'avg F\n',defgrad_av
|
||||
centroids = postprocessingMath.deformed_fft(res[0],res[1],res[2],dim,\
|
||||
numpy.reshape(values[:,column['tensor'][options.defgrad]:
|
||||
column['tensor'][options.defgrad]+9],
|
||||
(res[0],res[1],res[2],3,3)),defgrad_av,1.0)
|
||||
numpy.reshape(values[:,column['tensor'][options.defgrad]:
|
||||
column['tensor'][options.defgrad]+9],
|
||||
(res[0],res[1],res[2],3,3)),defgrad_av,1.0)
|
||||
print 'centroids\n',centroids
|
||||
ms = postprocessingMath.mesh(res[0],res[1],res[2],dim,defgrad_av,centroids)
|
||||
print 'mesh\n',ms
|
||||
|
||||
fields = {\
|
||||
'tensors': {},\
|
||||
|
|
|
|||
|
|
@ -274,6 +274,7 @@ end function math_mul33x33
|
|||
subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes)
|
||||
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
||||
implicit none
|
||||
|
||||
real*8 geomdim(3)
|
||||
integer res_x, res_y, res_z
|
||||
real*8 wrappedCentroids(res_x+2,res_y+2,res_z+2,3)
|
||||
|
|
@ -281,14 +282,14 @@ subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes)
|
|||
real*8 centroids(res_x ,res_y ,res_z ,3)
|
||||
|
||||
integer, dimension(3,8) :: neighbor = reshape((/ &
|
||||
0, 0, 0,&
|
||||
1, 0, 0,&
|
||||
1, 1, 0,&
|
||||
0, 1, 0,&
|
||||
0, 0, 1,&
|
||||
1, 0, 1,&
|
||||
1, 1, 1,&
|
||||
0, 1, 1 &
|
||||
0, 0, 0, &
|
||||
1, 0, 0, &
|
||||
1, 1, 0, &
|
||||
0, 1, 0, &
|
||||
0, 0, 1, &
|
||||
1, 0, 1, &
|
||||
1, 1, 1, &
|
||||
0, 1, 1 &
|
||||
/), &
|
||||
(/3,8/))
|
||||
|
||||
|
|
@ -296,33 +297,36 @@ subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes)
|
|||
real*8, dimension(3,3) :: defgrad_av
|
||||
integer, dimension(3) :: diag, shift, lookup, me, res
|
||||
|
||||
nodes = 0.0
|
||||
diag = 1
|
||||
shift = 0
|
||||
lookup = 0
|
||||
|
||||
res = (/res_x,res_y,res_z/)
|
||||
|
||||
wrappedCentroids=0.0
|
||||
wrappedCentroids = 0.0
|
||||
wrappedCentroids(2:res_x+1,2:res_y+1,2:res_z+1,:) = centroids
|
||||
|
||||
do k=0, res_z+1
|
||||
do j=0, res_y+1
|
||||
do i=0, res_x+1
|
||||
if (k==0 .or. k==res_z+1 .or. &
|
||||
j==0 .or. j==res_y+1 .or. &
|
||||
i==0 .or. i==res_x+1 ) then
|
||||
me = (/i,j,k/)
|
||||
do k = 0,res_z+1
|
||||
do j = 0,res_y+1
|
||||
do i = 0,res_x+1
|
||||
if (k==0 .or. k==res_z+1 .or. & ! z skin
|
||||
j==0 .or. j==res_y+1 .or. & ! y skin
|
||||
i==0 .or. i==res_x+1 ) then ! x skin
|
||||
me = (/i,j,k/) ! me on skin
|
||||
shift = sign(abs(res+diag-2*me)/(res+diag),res+diag-2*me)
|
||||
lookup = me-diag+shift*res
|
||||
wrappedCentroids(i+1,j+1,k+1,:) = centroids(lookup(1)+1,lookup(2)+1,lookup(3)+1,:)- &
|
||||
matmul(defgrad_av, shift*geomdim)
|
||||
wrappedCentroids(i+1,j+1,k+1,:) = centroids(lookup(1)+1,lookup(2)+1,lookup(3)+1,:) - &
|
||||
matmul(defgrad_av, shift*geomdim)
|
||||
endif
|
||||
enddo; enddo; enddo
|
||||
do k=0, res_z
|
||||
do j=0, res_y
|
||||
do i=0, res_x
|
||||
do n=1,8
|
||||
nodes(i+1,j+1,k+1,:) = nodes(i+1,j+1,k+1,:) + wrappedCentroids(i+1+neighbor(n,1),j+1+neighbor(n,2),k+1+neighbor(n,3),:)
|
||||
do k = 0,res_z
|
||||
do j = 0,res_y
|
||||
do i = 0,res_x
|
||||
do n = 1,8
|
||||
nodes(i+1,j+1,k+1,:) = nodes(i+1,j+1,k+1,:) + wrappedCentroids(i+1+neighbor(1,n), &
|
||||
j+1+neighbor(2,n), &
|
||||
k+1+neighbor(3,n), :)
|
||||
enddo; enddo; enddo; enddo
|
||||
nodes = nodes/8.0
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue