I have, using Fenics solved a pde on a spherical mesh, and saved that solution in a pvd file. I can plot it on a sphere (it just becomes a sphere when I open it in paraview). The problem is that a plot on a sphere is not very useful, especially if the interactiveness is removed. I am completely new to plotting in paraview, but what I would like to do is either to move to polar coordinates and plot my solution in the phi-theta plane, or use a cartographic projection of some sort to obtain a plot on the plane. Can this be achieved somehow? I have tried to use both the calculator and python calculator filter, but nothing works. I have also tried to create a square at one of the poles of the sphere and then run the “resample from data” filter, but nothing seems to work! Any suggestions?
So, I found this, https://blog.kitware.com/paraviews-python-programmable-filters-in-geophysics/, and I try to use the first algorithm. I paste the code into a programmable filter, and paraview just crashes.
See code below.
import numpy as np
from vtk.numpy_interface import algorithms as algs
from vtk.numpy_interface
import dataset_adapter as dsa
executive = self.GetExecutive()
outInfo = executive.GetOutputInformation(0)
exts = [executive.UPDATE_EXTENT().Get(outInfo, i) for i in xrange(6)]
whole = [executive.WHOLE_EXTENT().Get(outInfo, i) for i in xrange(6)]
global_dims = ([whole[1]-whole[0]+1,
whole[3]-whole[2]+1, whole[5]-whole[4]+1])
sg0 = vtk.vtkStructuredGrid()
sg0.SetExtent(exts)
Radius = 0.8
ThetaAxis = np.linspace(-np.pi,np.pi,
global_dims[0])
PhiAxis = np.linspace(-np.pi/2,np.pi/2,
global_dims[1])
xc, zc = np.meshgrid(ThetaAxis,PhiAxis,
indexing="xy")
Xc = Radius * np.sin(xc)
Yc = Radius * np.cos(xc)
coordinates = algs.make_vector(Xc.ravel(),
Yc.ravel(), Radius * zc.ravel())
pts = vtk.vtkPoints()
pts.SetData(dsa.numpyTovtkDataArray(coordinates,"Points"))
sg0.SetPoints(pts)
sg1 = vtk.vtkStructuredGrid()
sg1.SetExtent(exts)
zc_other = np.zeros(xc.size).reshape(xc.shape)
coordinates = algs.make_vector(xc.ravel(),
zc.ravel(), zc_other.ravel())
pts = vtk.vtkPoints()
pts.SetData(dsa.numpyTovtkDataArray(coordinates,"Points"))
sg1.SetPoints(pts)
output.SetBlock(0, sg0)
output.SetBlock(1, sg1)
I have also tried using the code in a python script, but it won’t even run.