numpy array to data object

Asked by Tom

Hi ,

I want to specify an initial temperature field. I will generate a numpy array with a shape matching the Solution Space.

I am stuck getting each of theses data points from the numpy array to the corresponding positions in the Solution Space.

heres a snippet of code.
# create simple domain
domWidth = 1.*U.m

domHeight = 1*U.m
numDivsX = 9

mydat1 = Data(value=np.ones((nx,ny)), what=Function(dom),expanded=True)
numDivsY = 9

dom = Rectangle(numDivsX, numDivsY, l0=domWidth, l1=domHeight)
# The number of data points in this Domain will correspond to
nx= numDivsX+1
ny=numDivsY+1

# create a Data object from a numpy array
mydat1 = Data(value=np.ones((nx,ny)), what=Function(dom),expanded=True)

I get an error that python argument types dont match the C++ signature.

Im trying to get the array elements to the corresponding positions in the domain.

Cheers

Tom

Question information

Language:
English Edit question
Status:
Solved
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Solved by:
Tom
Solved:
Last query:
Last reply:

This question was reopened

  • by Tom
Revision history for this message
Joel Fenwick (j-fenwick1) said :
#1

escript does not support directly loading an array into a Data object. Unless you want all points to have the same value.

You will need to use table interpolation for this.
That is, interpolate from your prepared array onto the domain.
The user guide has details.

Revision history for this message
Tom (chienenrage3030) said :
#2

thanks Joel

Revision history for this message
Tom (chienenrage3030) said :
#3

If I am reading the function description correctly,
(at http://esys.esscc.uq.edu.au/esys13/3.2/epydoc/esys.escript.escriptcpp.Data-class.html#interpolateTable)cify

This method will only allow me to specify 2D data at each point on the domain. I want to specfiy a velocity field in 3D, so I need to interpolate a 3D vector field to a Data object. Is this possible ?

Thanks

Tom

Revision history for this message
Joel Fenwick (j-fenwick1) said :
#4

The 2D version was added first.

The combined interface is the function (not method) interpolateTable().

This will deal with both 2D and 3D depending on what dimension of table you pass in.

>>> help(interpolateTable)
>>> Help on function interpolateTable in module esys.escript.util:

interpolateTable(tab, dat, start, step, undef=1.0000000000000001e+50, check_boundaries=False)

Revision history for this message
Tom (chienenrage3030) said :
#5

Hi Joel,

thanks for your help on this.

I dont think I was clear, I need to specify a vector field, so at each point in the domain, i need a vector, so over a 2D rectangular domain with vx and vy, i need to specify the vx and vy at each point.

I want to load in a vector field, but all the functions I have experimented with so far seem to be limited to specifying just rank 0 or 1dimensional data at each point in the domain.

if interpolate table cant help me, what about setting slices? below is my attempt to get this working.

thanks

Tom
#######################################################

# import

import os, sys

from esys.escript import *

from esys.escript.linearPDEs import LinearPDE

from esys.escript.datamanager import DataManager

import esys.escript.unitsSI as U

from esys.finley import Rectangle

from esys.escript.linearPDEs import TransportPDE

from esys.escript.models import DarcyFlow

from esys.weipa import saveVTK

import numpy as np

from esys.escript import saveDataCSV
##########################################
# set up the domain

domWidth = 1.*U.m

domHeight = 1*U.m

numDivsX = 1

numDivsY = 1

dom = Rectangle(numDivsX, numDivsY, l0=domWidth, l1=domHeight)
##########################################################################

def toXYarray(coords):

 coords = np.array(coords.toListOfTuples()) # convert to tuple

 coordX = coords[:,0] # X components

 coordY = coords[:,1] # Y components

 return coordX,coordY

############################################################################
# Describe the Domain
locations = Solution(dom).getX()#
nx= numDivsX+1

ny=numDivsY+1

#mydat1 = Data(value=np.ones((nx,ny)), what=Function(dom),expanded=True)

#mydat2 = Data(value=[[1,1],[1,1],[1,1]],what=ContinuousFunction(dom))
vt3 = Data(0,what=Function(dom),shape=(1,2))
x = dom.getX()

# get Slice
print vt3.getShape()
print vt3.getRank()
coordX,coordY = toXYarray(locations)
print coordX,coordY
vt3.getFunctionSpace().getX()
print 'slice'
print vt3[0,]
print 'slice2'
print vt3[:,:]
test = np.array([1,0])
vt3[:,:] = [test]
locations = dom.getX()
saveDataCSV('vt3.csv',v=vt3)
saveDataCSV('locations.csv',xs = locations)
vec = Vector(0,ContinuousFunction(dom))
print 'vec shape'
print vec.getShape()
print 'vec rank'
print vec.getRank()
print vec[0,:]
vec[:] = np.array([8,78])

print
#vec[0:]= np.array([0,1,2,3,4,5,6,7])
#vec[0:]= np.array([[0],[1],[2],[3],[4],[5],[6],[7]])

saveDataCSV('vec.csv',v=vec)

t = Data(1.,(4,4,6,6), Function(dom))
print 't shape is'
print t.getShape()
print 't rank is'
print t.getRank()#######################################################

# import

import os, sys

from esys.escript import *

from esys.escript.linearPDEs import LinearPDE

from esys.escript.datamanager import DataManager

import esys.escript.unitsSI as U

from esys.finley import Rectangle

from esys.escript.linearPDEs import TransportPDE

from esys.escript.models import DarcyFlow

from esys.weipa import saveVTK

import numpy as np

from esys.escript import saveDataCSV
##########################################
# set up the domain

domWidth = 1.*U.m

domHeight = 1*U.m

numDivsX = 1

numDivsY = 1

dom = Rectangle(numDivsX, numDivsY, l0=domWidth, l1=domHeight)
##########################################################################

def toXYarray(coords):

 coords = np.array(coords.toListOfTuples()) # convert to tuple

 coordX = coords[:,0] # X components

 coordY = coords[:,1] # Y components

 return coordX,coordY

############################################################################
# Describe the Domain
locations = Solution(dom).getX()#
nx= numDivsX+1

ny=numDivsY+1

#mydat1 = Data(value=np.ones((nx,ny)), what=Function(dom),expanded=True)

#mydat2 = Data(value=[[1,1],[1,1],[1,1]],what=ContinuousFunction(dom))
vt3 = Data(0,what=Function(dom),shape=(1,2))
x = dom.getX()

# get Slice
print vt3.getShape()
print vt3.getRank()
coordX,coordY = toXYarray(locations)
print coordX,coordY
vt3.getFunctionSpace().getX()
print 'slice'
print vt3[0,]
print 'slice2'
print vt3[:,:]
test = np.array([1,0])
vt3[:,:] = [test]
locations = dom.getX()
saveDataCSV('vt3.csv',v=vt3)
saveDataCSV('locations.csv',xs = locations)
vec = Vector(0,ContinuousFunction(dom))
print 'vec shape'
print vec.getShape()
print 'vec rank'
print vec.getRank()
print vec[0,:]
vec[:] = np.array([8,78])

print
#vec[0:]= np.array([0,1,2,3,4,5,6,7])
#vec[0:]= np.array([[0],[1],[2],[3],[4],[5],[6],[7]])

saveDataCSV('vec.csv',v=vec)

t = Data(1.,(4,4,6,6), Function(dom))
print 't shape is'
print t.getShape()
print 't rank is'
print t.getRank()
t[1,1,1,0] =9
s = t[:2,:,2:6,5]
print s.getRank()
print s
t[:2,:,2:6,5] = np.zeros((2,4,4))
saveDataCSV('t.csv',r=t)
t[1,1,1,0] =9
s = t[:2,:,2:6,5]
print s.getRank()
print s
t[:2,:,2:6,5] = np.zeros((2,4,4))
saveDataCSV('t.csv',r=t)

Revision history for this message
Joel Fenwick (j-fenwick1) said :
#6

Hmm,
   I'm not really a fan of using slice if we don't have to.
How about:

using separate interpolation steps to produce vx and vy then

V=fx*[1,0]+fy*[0,1]

Revision history for this message
Tom (chienenrage3030) said :
#7

Hi Joel,

that will work. But,,

interpolateTable is limited to interpolating a function which is constant across the perpendicular dimension.

What I mean is, the function f(x) must be constant in the y-direction. Even for 2D interpolations, it will only work if the functions are uniform in the x and y direction.

The field I am working with varies in both x and y. it is composed of sin(x)cos(y) terms and similar.

The thing about slice i was misunderstanding was that it would give me access to an array of values, it doesn't as far as I can figure. such as all the x0 components of a table of [x0's,x1's]

my script is below.

1.Could you show me if it is possible to take a list of values and assign it to the corresponding positions of a 2D domain, by adjusting the domain spacing and value spacing accordingly?

2. can i append data objects, and make the field by lists one row at a time...
from esys.escript import saveDataCSV, sup, interpolateTable
import numpy as np
from esys.finley import Rectangle
from math import *
from esys.escript import whereZero
from esys.escript import *

# create a rectangular domain
n = 8
r = Rectangle(n,n)
numDivsX = 2
numDivsY = 2
# will need lists of the same length
lengthX = 1.
lengthY = 1.
dom = Rectangle(numDivsX, numDivsY, l0=lengthX, l1=lengthY)
xdom = dom.getX()
toobig = 100
xlin = np.linspace(0,1,numDivsX)
ylin = np.linspace(0,1,numDivsY)
# create a table of values to interpolate from

table = xlin
###############################################
# create the function for vx
x = np.linspace(0,lengthX,numDivsX+1) # will need normalised x,y if the domain is greater than 1 by 1
y = np.linspace(0,lengthY,numDivsY+1)
A = 1
B = 1

u1 = -3*sqrt(2)*pi**2*A*np.sin(x*pi/sqrt(2))*np.cos(pi*y) - 3*sqrt(2)*pi**2*B*np.sin(pi*sqrt(2)*x)*np.cos(pi*y)

v1 = -3*pi**2*A*np.cos(x*pi/sqrt(2))*np.sin(pi*y) - 6*pi**2*B*np.cos(pi*sqrt(2)*x)*np.sin(pi*y)

T1 = [] # I want T1 to be a 2D table
for j in range(len(y)):
 tmp = []
 for i in range(len(x)):
  tmp.append(2*A*np.cos(x[i]*pi/sqrt(2))*np.sin(pi*y[j]) + 2*B*np.cos(pi*sqrt(2)*x[i])*np.sin(pi*y[j]))
 tmp = np.array(tmp)
 T1 = tmp # need an x array and a y array
 T2 = [tmp,tmp]

testarray = np.linspace(0,8,9)
print 'test array'
print [testarray,testarray]
###############################################
#create the function for vy
step = lengthX/(numDivsX)
numslices=len(T1)-1
minval= 0
maxval= 1
#Tstep = (maxval - minval)/numslices #The width of the gap between entries in the table
xstep = step
ystep = step

# interpolate these initial conditions to a table over the domain
Tmin = 0
Ttoobig = 1

T1table = interpolateTable(T1, xdom[0], minval, step, toobig)
testTable = interpolateTable(testarray, xdom[0], 0, 1, toobig)
T1table2 = interpolateTable(T2,xdom,(minval,minval),(xstep,ystep),toobig)
testTable2 = interpolateTable([testarray,testarray],xdom,(0,0),(1,1),toobig)

###############################################

result=interpolateTable(table, xdom[0], minval, xstep, toobig)
vx = result*[1,0]
vy = result*[0,1]
v = vx + vy
saveDataCSV('x.csv',locs = xdom)
saveDataCSV('result.csv', output = v) # save something
saveDataCSV('Temp1.csv', output = T1table) # save T1 initial condition
saveDataCSV('Temp2.csv', output = T1table2) # save T1 initial condition
saveDataCSV('test2.csv', output = testTable2) # save T1 initial condition

Revision history for this message
Tom (chienenrage3030) said :
#8

Problem solved.