Campos vectoriales en Blender

En este post se enseña cómo representar en Blender archivos de datos que contengan las componentes de los vectores que forman un campo vectorial. El resultado podría ser el siguiente:

En este ejemplo se enseña cómo cargar archivos generados por OOMMF, pero vale cualquier archivo de texto que tenga las componentes de los vectores en tres columnas, ordenadas en orden FORTRAN. Es decir:
x11 y11 z11
...  ...  ...
x1N y1N z1N
x2N y2N z2N
...  ...  ...
xMN yMN zMN

Se puede extender a 3D, pero en este ejemplo se usa 2D, una matriz NxM. Además, en este caso la matriz es cuadrada (M=N). En python se usa la función reshape del módulo numpy para reordenar estas listas de ongitud NxM en matrices de tamaño (N,M).
Importante: Este script necesita numpy, que no viene instalado en las versiones de Blenders más antiguas. En la versión 2.7 parece que sí viene, desde luego en la actual (2.7.2 al momento de escribir esto), sí viene.

Lo primero es en Blender crear una flecha, o el objeto que se quiera utilizar para representar los vectores. (En el script se puede descomentar una línea que ahora hay comentada para usar cilindros). Pero si hemos creado la flecha, seleccionándola con el ratón (si no está seleccionada hay un error al ejecutar el script), se ejecuta en el "Text editor" de Blender  el script que hay más abajo. A la flecha original hay que llamarla "Arrow" (sin comillas) para que funcione el script.

Flecha seleccionada antes de ejecutar el script

Este script abre el archivo, lee los contenidos, separa las componentes x, y, z, las reordena en matrices y para cada entrada de la matriz coloca en una posición dada por el índice "i, j" una flecha, cilindo, o cualquier objeto que se quiera, orientándolo según las componentes de los vectores. Es importante usar math.atan2 para el calcular los ángulos de rotación, y no math.atan, ya que ocurren cosas raras y algunas flechas apuntan en el sentido opuesto al correcto. En esta otra entrada se explica la idea de cómo convertir coordenadas angulares a rotaciones en Blender, que no es trivial.

 Por lo que se ve, cuantas más flechas hay ya colocadas, más le cuesta crear nuevas flechas, como si Blender perdiera bastante tiempo contando todos los objetos que ya hay creados , por lo que este script tiene un cuello de botella tipo O(N).

Tras crear las flechas, se crea un plano debajo, que en el ejemplo tiene un material con textura de tablero de ajedrez. Para ello, hay que seleccionar con botón derecho el plano, y en el material, siendo difusivo, por ejemplo, hacer click en el botón que hay junto al campo de color. En el desplegable, seleccionar checker texture. Luego se puede cambiar el color de cada uno de los dos tipos de cuadros, o poner tantos cuadros como flechas (en el campo "scale"). Debajo de este plano hay colocado (manualmente, no lo crea este script) un plano blanco para que los cuadros blancos del primer plano sean del mismo color del fondo y no se distingan de este.

El script hay que cargarlo en el text editor, y ejecutarlo con el botón "run script".



Cosas que hay que cambiar en el script para cada caso:

1. El nombre del archivo del que se leen los datos. Ahora se llama "a.omf". Para que blender lo encuentre, hay que situarse en la carpeta donde está el archivo. Una forma es ejecutar blender desde una terminal, habiéndose situado antes esta en el directorio donde está guardado el archivo.
2. La longitud de la cabecera, para leer solo números e ignorar las líneas que haya al principio del archivo de datos. Si no hay, se pone 0. En este ejemplo, hay 34 (skiprows=34). 
3. El delimitador que separa las columnas x, y, z del archivo. Ahora hay espacios (delimiter= ' '). Para poner tabuladores hay que usar delimiter='\t', y para usar comas, delimiter=','.
4. El número de filas y columnas en que hay que reordenar estas matrices. Son las variables "rows" y "cols" (sin las comillas).
5. En la gráfica final se interpola el número de flechas, para que no sea un "caos" donde no se vea nada, reduciendo el número de columnas. Esto está dado por las variables "rows_red" y "cols_red". Si se pone un número muy grande será la mayor causa de lentitud de ejecución del script (50x50 empieza a ser un número grande en términos de esperar a que se ejecute el script).

Un par de detalles en cuanto a las flechas. Queda mejor si se marcan todas las flechas y con Ctrl-Alt-Shift-c se mueve el centro de cada flecha al "center of mass". Y para que tengan el mismo material que la original hay que marcar la original (botón derecho) y luego con la tecla b arrastrar un rectángulo para seleccionar todas las demás, pero con la original seleccionada (el borde se ve resaltado) y dar a Ctrl-l (letra L minúscula). En el submenú que sale, darle a Materials, para linkar los materiales de todas las flechas al de la original. Por último, la flecha original yo la pongo por encima del plano Z=0, para que no salga al renderizar, pero siempre está ahi.

Este es el script que hay que ejecutar con la flecha seleccionada. Copiar respetando las indentaciones


import numpy as np
import math
import bpy
import time

for ob in bpy.data.objects:
    if ob.name == "Arrow":
        model_obj = bpy.data.objects["Arrow"]
        print('xxx  ok')
        #break
       
    if "Arrow." in ob.name:
        bpy.data.objects[ob.name].select = True
        bpy.ops.object.delete()

#if ob.name == "Arrow":
model_obj = bpy.data.objects["Arrow"]
print('xxx  ok')

#bpy.data.objects["Arrow"].select = True


def place_obj(x, y, z, theta, phi, model_obj):
    #place cylinder
    #bpy.ops.mesh.primitive_cylinder_add(radius=0.1,depth=0.7,location=(i,j,0),rotation=(0,theta,math.pi/2-phi)) #rotation va en rads
    #place custom object
    pi = math.pi
    bpy.ops.object.duplicate()
    bpy.context.object.location = (i, j, 0)
    bpy.context.object.rotation_euler = (0, theta, pi/2-phi)

def interpolate_array_2d(rows, rows_red, cols, cols_red, A):
        A_red_1 = np.zeros((rows_red, cols))
        A_red_2 = np.zeros((rows_red, cols_red))
       
        for i in range(cols):
                A_red_1[:,i] = np.interp(np.linspace(0,rows,rows_red),range(rows),A[:,i])
       
        for i in range(rows_red):
                A_red_2[i,:] = np.interp(np.linspace(0,cols,cols_red),range(cols),A_red_1[i,:])
        return A_red_2


rows = 240
cols = 240

a = np.loadtxt('a.omf',skiprows=34,delimiter=' ')
mx = a[:,0]
my = a[:,1]
mz = a[:,2]

Mx = np.reshape(mx,(rows,cols))
My = np.reshape(my,(rows,cols))
Mz = np.reshape(mz,(rows,cols))

rows_red = 30
cols_red = 30

Mx_red = interpolate_array_2d(rows, rows_red, cols, cols_red, Mx)
My_red = interpolate_array_2d(rows, rows_red, cols, cols_red, My)
Mz_red = interpolate_array_2d(rows, rows_red, cols, cols_red, Mz)

for i in range(rows_red):
    print(i)
    for j in range(cols_red):
        if Mx_red[i,j] != 0 and My_red[i,j] != 0 and Mz_red[i,j] != 0:
            theta = math.atan2(math.sqrt(Mx_red[i,j]**2 + My_red[i,j]**2),Mz_red[i,j])
            phi = math.atan2(My_red[i,j],Mx_red[i,j])
            place_obj(i, j, 0, theta, phi, model_obj)
           

bpy.ops.mesh.primitive_plane_add(location=(rows_red/2-0.5,cols_red/2-0.5,0))
ob = bpy.context.active_object
ob.scale=(rows_red/2,cols_red/2,0)
# Get material
mat = bpy.data.materials.get("Material")
if mat is None:
    # create material
    mat = bpy.data.materials.new(name="Material")

# Assign it to object
if ob.data.materials:
    # assign to 1st material slot
    ob.data.materials[0] = mat
else:
    # no slots
    ob.data.materials.append(mat)