非定常レーザ計測シミュレーション

Programmable Sourceで時刻データにアクセスする方法についてはProducing Data with Timesteps (Source) | ParaView Wikiを参照のこと。

動作確認バージョン

ParaView 4.2.0
ParaView 5.0.0 RC1

Output Data Set Type

vtkPolyData

Script

from paraview import vtk
import numpy as np

pps = 700000 # points per second
frq = 15 # frequency
dsr = 32 # DSRs
s_ang = 90.0
s_len = 120.0
s_gnd = -2.7
s_hil = -2.5

def GetUpdateTimestep(algorithm):
   """Returns the requested time value, or None if not present"""
   executive = algorithm.GetExecutive()
   outInfo = executive.GetOutputInformation(0)
   if not outInfo.Has(executive.UPDATE_TIME_STEP()):
       return None
   return outInfo.Get(executive.UPDATE_TIME_STEP())

s_len -= GetUpdateTimestep(self) * 60.0 / 3.6

dsr_ang_v = np.array([(8.0 - x) * 4.0 / 3.0 / 180.0 * np.pi for x in range(9, 32)]).reshape((-1, 1))
#dsr_ang_h = np.array([2.0 * np.pi / pps * frq * dsr * x for x in range(0, 1458)]).reshape((1, -1))
dsr_ang_h = np.array([2.0 * np.pi / pps * frq * dsr * x for x in range(-364, 365)]).reshape((1, -1))
dsr_ang_h[np.where(dsr_ang_h == 0.0)] = 1.0e-100

ray = np.array((
np.sin(dsr_ang_v) * np.sin(dsr_ang_h),
np.sin(dsr_ang_v) * np.cos(dsr_ang_h),
np.cos(dsr_ang_v) * np.sin(dsr_ang_h) * np.cos(dsr_ang_h),
))

def getCloud(ang, length, height):
    t = np.tan(-ang / 180.0 * np.pi)
    plane = np.array((t, 0.0, 1.0, -length * t - height))
    x = -plane[3] / (plane[0] + ray[0] / ray[1] * plane[1] + ray[0] / ray[2] * plane[2])
    y = ray[0] / ray[1] * x
    z = ray[0] / ray[2] * x
    a = np.ones(np.shape(x)) * ang
    id = np.indices(np.shape(x)) + (32 - np.shape(x)[0])
    return np.vstack((x.reshape(-1), y.reshape(-1), z.reshape(-1), a.reshape(-1), id[0].reshape(-1)))

cloud_g = getCloud(0.0, s_len, s_gnd)
cloud_s = getCloud(s_ang, s_len, s_gnd)
cloud_h = getCloud(0.0, s_len, s_hil)

for i in range(np.shape(cloud_g)[1]):
    if cloud_g[0, i] <= s_len:
        continue
    elif cloud_s[2, i] <= s_hil:
        cloud_g[:, i] = cloud_s[:, i]
    else:
        cloud_g[:, i] = cloud_h[:, i]

cloud_g[np.where((-1.0e-20 < cloud_g) & (cloud_g < 1.0e-20))] = 0.0

pdo = self.GetPolyDataOutput() # pdo = self.GetOutput()
pts = vtk.vtkPoints()
for i in range(np.shape(cloud_g)[1]):
    pts.InsertPoint(i, cloud_g[0, i], cloud_g[1, i], cloud_g[2, i] - s_gnd)
pdo.SetPoints(pts)

pdo.Allocate(np.shape(cloud_g)[1], 1)
ang = vtk.vtkDoubleArray()
ang.SetName('angle')
ang.SetNumberOfComponents(1)
ang.SetNumberOfTuples(np.shape(cloud_g)[1])
id = vtk.vtkDoubleArray()
id.SetName('id')
id.SetNumberOfComponents(1)
id.SetNumberOfTuples(np.shape(cloud_g)[1])
for i in range(np.shape(cloud_g)[1]):
    vrt = vtk.vtkVertex()
    vrt.GetPointIds().SetNumberOfIds(1)
    vrt.GetPointIds().SetId(0, i)
    pdo.InsertNextCell(vrt.GetCellType(), vrt.GetPointIds())
    ang.SetValue(i, cloud_g[3, i])
    id.SetValue(i, cloud_g[4, i])
pdo.GetPointData().AddArray(ang)
pdo.GetPointData().AddArray(id)

Script (RequestInformation)

import numpy as np

def SetOutputTimesteps(algorithm, timesteps):
    executive = algorithm.GetExecutive()
    outInfo = executive.GetOutputInformation(0)
    outInfo.Remove(executive.TIME_STEPS())
    for timestep in timesteps:
        outInfo.Append(executive.TIME_STEPS(), timestep)
    outInfo.Remove(executive.TIME_RANGE())
    outInfo.Append(executive.TIME_RANGE(), timesteps[0])
    outInfo.Append(executive.TIME_RANGE(), timesteps[-1])
SetOutputTimesteps(self, np.arange(0, 8, 0.1))