Cómo hacer un diagrama vectorial VTK en Nutils

El siguiente código genera un vecs campo de vector 2D aleatorio y luego produce un gráfico de carcaj Matplotlib y un archivo de datos VTK del mismo campo.

desde nutils import mesh, plot
importar numpy

dominio, geom = mesh.rectilinear ([rango (4), rango (4)])
base = dominio.basis (‘estándar’, grado = 1) .vector (2)
vecs = base.dot (numpy.random.normal (tamaño = len (base)))

# matplotlib quiver plot
geom_, vecs_ = dominio.elem_eval ([geom, vecs], ischeme = ‘uniform2’, separe = Falso)
con plot.PyPlot (‘vecs’) como plt:
plt.vectors (geom_, vecs_)

# vtk archivo de datos
geom_, pvecs_ = dominio.elem_eval ([geom, vecs], ischeme = ‘vtk’, separe = Verdadero)
cvecs_ = dominio.elem_eval (vecs, ischeme = ‘uniform1’, separe = Falso)
con plot.VTKFile (‘vecs’) como vtk:
vtk.unstructuredgrid (geom_)
vtk.pointdataarray (‘pvecs’, pvecs_)
vtk.celldataarray (‘cvecs’, cvecs_)

Algunas observaciones:

  1. El archivo VTK resultante contiene el campo vectorial como datos de punto ( pvecs ) y datos de celda ( cvecs ), este último simplemente contiene el valor muestreado en el punto medio del elemento ( 'uniform1' ). Cuál de los dos es más apropiado depende de la aplicación.
  2. Los vectores en pvecs se muestrean para cada elemento y, por lo tanto, se repiten en caso de que los nodos sean compartidos por elementos. Esto no es visible si el campo vectorial es continuo, como con vecs en el ejemplo, pero para campos vectoriales discontinuos dará como resultado varios vectores que se originan desde el mismo punto.
  3. Debido a que el formato de datos VTK es inherentemente tridimensional, al cargar el archivo vecs.vtk en paraview se muestran pvecs y cvecs con los componentes x, y y z. Este último componente es cero.
  4. Al modificar las líneas 4 y 5, el código puede hacerse en 3D. No será necesario realizar cambios en la generación VTK, que en este caso generará un archivo de datos con componentes z distintos de cero.