Dear Paraviewers, I have very recently begun to use PV (been a VisItor up to this point) so there are a lot of features with which I am not entirely familiar. I am trying to run a python script which does some math then stream tracing on two vectors variables "e" and "b". My code dumps out components of each vector field , which I read into PV where they appear as point data.. I t hen combine this data into vectors "e" and "b" with PV's calculator feature. After this step I run my python script from Tools->Python Shell. The script (attached it as fieldl_line.py -- looks terrible when opened with Vim) does not get very far and complains about various things. I attached the log file entitled PVlog. I cannot decipher some of these messages but I have a strong suspicion the inputs variable does not get populated correctly or the script cannot find "e" and "b". Any ideas where I have gone wrong? Is there an intermediate step after defining new variables on the calculator to make the script aware of the existence of "e" and "b". I plotted these new vector variables on PV without any issue.
Thank you very much for your time and interest. -- Cihan -- Cihan -- Cihan
PVlog
Description: Binary data
# create a new 'Programmable Filter'
programmableFilter1 = ProgrammableFilter()
programmableFilter1.OutputDataSetType = 'vtkMultiblockDataSet'
programmableFilter1.Script = 'import numpy\n\nDIMS = (30, 1, 30)\norigins = numpy.array([\n (-1, 0, 0), (1, 0, 0),\n (0, 1, 0),\n (0, 0, -1), (0, 0, 1),\n (-1, 1, 0), (1, 1, 0),\n (0, 1, -1), (0, 1, 1)\n])\nnxplanes = 10\n\nfrom vtk.vtkFiltersFlowPaths import vtkStreamTracer\nfrom vtk.vtkFiltersParallelFlowPaths import vtkPStreamTracer\nfrom vtk.vtkFiltersSources import vtkPlaneSource\nfrom vtk.vtkCommonExecutionModel import vtkStreamingDemandDrivenPipeline\nfrom vtk.vtkCommonSystem import vtkTimerLog\nfrom vtk.numpy_interface import dataset_adapter as dsa\nfrom vtk.numpy_interface import algorithms as algs\nfrom mpi4py import MPI\ncomm = MPI.COMM_WORLD\n\nipt = inputs[0].VTKObject\nif ipt.IsA("vtkMultiBlockDataSet"):\n cp0 = ipt.GetBlock(0)\n ssample = 4\nelse:\n cp0 = ipt\n ssample = 1\n\nmb = vtk.vtkMultiBlockDataSet()\nmb.SetBlock(0, cp0)\n\nspacing = cp0.GetSpacing()\nextent = self.GetInputInformation(0, 0).Get(vtkStreamingDemandDrivenPipeline.WHOLE_EXTENT())\nextent = [e*ssample for e in extent]\nbounds = [s * e for (s,e) in zip(spacing, extent[1:6:2])]\n\n# Change where the y plane is.\nyplane = 0\n\nidx = 1\nfor origin in origins:\n cp = vtk.vtkImageData()\n cp.ShallowCopy(cp0)\n cp.SetOrigin(origin * bounds)\n mb.SetBlock(idx, cp)\n idx += 1\n\nfrom vtk import vtkParallelCore as pc\ngc = pc.vtkMultiProcessController.GetGlobalController()\n\ndef calculate_integral(xbeg, xend, nx):\n st = vtkPStreamTracer()\n st.SetInputData(mb)\n st.SetInputArrayToProcess(0, 0, 0, 0, "b")\n st.SetIntegratorTypeToRungeKutta4()\n st.SetInitialIntegrationStep(0.4)\n st.SetIntegrationStepUnit(st.CELL_LENGTH_UNIT)\n st.SetMaximumPropagation(250)\n st.SetIntegrationDirectionToForward()\n\n\n img2 = vtk.vtkImageData()\n img2.SetDimensions(nx, DIMS[1], DIMS[2])\n img2.SetOrigin(xbeg, yplane, 1)\n img2.SetSpacing((xend-xbeg) / (nx-1), 1, (bounds[2]-2)/(DIMS[2]-1))\n\n pl = vtkPlaneSource()\n pl.SetOrigin(xbeg, yplane, 1)\n pl.SetPoint1(xend , yplane, 1)\n pl.SetPoint2(xbeg, yplane, bounds[2] - 1)\n pl.SetResolution(nx - 1, DIMS[2]-1)\n\n st.SetSourceConnection(pl.GetOutputPort())\n vtkTimerLog.MarkStartEvent("Streamline integration")\n st.Update()\n vtkTimerLog.MarkEndEvent("Streamline integration")\n opt = dsa.WrapDataObject(st.GetOutput())\n\n sline = st.GetOutput()\n\n # Change the equation here\n # If ncomps changes, also change SetTuple3() and GetTuple3()\n # below\n e = opt.PointData["e"]\n b = opt.PointData["b"]\n print algs.dot(e,b)\n ia = algs.dot(e, algs.norm(b))\n ncomps = 1\n\n integral = vtk.vtkFloatArray()\n integral.SetName("integral")\n integral.SetNumberOfComponents(ncomps)\n integral.SetNumberOfTuples(img2.GetNumberOfPoints())\n for i in range(ncomps):\n integral.FillComponent(i, 0)\n\n lengtha = vtk.vtkFloatArray()\n lengtha.SetName("length")\n lengtha.SetNumberOfTuples(img2.GetNumberOfPoints())\n lengtha.FillComponent(0, 0)\n\n img2.GetPointData().AddArray(integral)\n img2.GetPointData().AddArray(lengtha)\n\n pts = opt.Points\n sids = opt.CellData["SeedIds"]\n lines = opt.GetLines()\n lines = dsa.vtkDataArrayToVTKArray(lines.GetData())\n length = len(lines)\n index = 0\n cellid = 0\n\n vtkTimerLog.MarkStartEvent("Integration")\n while index < length:\n npts = lines[index]\n index += 1\n if npts > 0:\n cell = lines[index:index+npts]\n index += npts\n delta = pts[cell[1:]] - pts[cell[0:-1]]\n itg = numpy.sum(algs.mag(delta) * ((ia[cell[1:]] + ia[cell[0:-1]]) / 2), axis=0)\n itg += integral.GetTuple1(sids[cellid])\n integral.SetTuple1(sids[cellid], itg)\n\n ln = numpy.sum(algs.mag(delta), axis=0)\n ln += lengtha.GetTuple1(sids[cellid])\n lengtha.SetTuple1(sids[cellid], ln)\n cellid += 1\n vtkTimerLog.MarkEndEvent("Integration")\n\n for array in (integral, lengtha):\n itg = dsa.vtkDataArrayToVTKArray(array)\n red = numpy.empty(itg.shape, dtype=numpy.float32)\n comm.Reduce([itg, MPI.FLOAT], [red, MPI.FLOAT], op=MPI.SUM)\n\n itg[:] = red\n\n return sline, img2\n\ntmpOutput = vtk.vtkMultiBlockDataSet()\n\nxbeg = 1.\nnseeds = DIMS[0] / nxplanes\ndx = (bounds[0] - 1 - xbeg) / (DIMS[0] - 1)\nidx = 0\ntmpOutput.SetNumberOfBlocks(nxplanes)\nfor plane in range(nxplanes):\n xend = xbeg + (nseeds - 1) * dx\n sline, img = calculate_integral(xbeg, xend, nseeds)\n# tmpOutput.SetBlock(idx, sline)\n# idx += 1\n if gc.GetLocalProcessId() == 0:\n tmpOutput.SetBlock(idx, img)\n else:\n tmpOutput.SetBlock(idx, None)\n idx += 1\n xbeg += nseeds * dx\n\nif gc.GetLocalProcessId() != 0:\n output.SetNumberOfBlocks(1)\n output.SetBlock(0, None)\n return\n\ntmpOutput = dsa.WrapDataObject(tmpOutput)\n\nbds = [1,0,1,0,1,0]\nn = 0\nfor ds in tmpOutput:\n if ds.IsA("vtkImageData"):\n n += 1\n spacing = ds.GetSpacing()\n dbds = ds.GetBounds()\n for i in range(3):\n bds[2*i+1] = max(bds[2*i+1], dbds[2*i+1])\n bds[2*i] = min(bds[2*i], dbds[2*i])\nbds = numpy.array(bds)\ndims = numpy.round((bds[1:6:2] - bds[0:6:2])/ spacing).astype(numpy.int32)\n\noarray = numpy.empty((dims[0]+1, 1, dims[2]+1), order=\'F\', dtype=numpy.float64)\nolength = numpy.empty((dims[0]+1, 1, dims[2]+1), order=\'F\', dtype=numpy.float64)\n\n\nfor ds in tmpOutput:\n if ds.IsA("vtkImageData"):\n origin = numpy.array(ds.GetOrigin())\n origin = numpy.round((origin - bds[0:6:2]) / spacing).astype(numpy.int32)\n dims_ = ds.GetDimensions()\n integral = ds.PointData[\'integral\']\n integral = integral.ravel().reshape(dims_, order=\'F\')\n oarray[origin[0]:origin[0]+dims_[0], :, origin[2]:origin[2]+dims_[2]] = integral\n length = ds.PointData[\'length\']\n length = length.ravel().reshape(dims_, order=\'F\')\n olength[origin[0]:origin[0]+dims_[0], :, origin[2]:origin[2]+dims_[2]] = length\n\nimage = vtk.vtkImageData()\nimage.SetOrigin(bds[0:6:2])\nimage.SetSpacing(spacing)\nimage.SetDimensions(dims[0]+1, 1, dims[2]+1)\n\nwimage = dsa.WrapDataObject(image)\nwimage.PointData.append(oarray.ravel(\'F\'), \'integral\')\nwimage.PointData.append(olength.ravel(\'F\'), \'length\')\n\noutput.SetNumberOfBlocks(1)\noutput.SetBlock(0, image)'
programmableFilter1.RequestInformationScript = ''
programmableFilter1.RequestUpdateExtentScript = 'io = self.GetInputInformation(0, 0)\nfrom vtk import vtkCommonExecutionModel as em\net = em.vtkExtentTranslator\nio.Set(et.UPDATE_SPLIT_MODE(), et.Z_SLAB_MODE)\n'
programmableFilter1.PythonPath = ''
_______________________________________________ Powered by www.kitware.com Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html Please keep messages on-topic and check the ParaView Wiki at: http://paraview.org/Wiki/ParaView Search the list archives at: http://markmail.org/search/?q=ParaView Follow this link to subscribe/unsubscribe: http://public.kitware.com/mailman/listinfo/paraview
