¿Cómo calculo la diferencia numérica entre dos campos almacenados en dos archivos VTK diferentes con la misma estructura?

Supongamos que tengo dos archivos VTK, ambos en formato de cuadrícula estructurada. Las cuadrículas estructuradas son las mismas (tienen la misma lista de puntos, en el mismo orden), y hay un campo, llamado "Phi", en cada archivo VTK. Quiero crear un tercer archivo VTK, nuevamente con la misma cuadrícula estructurada, y trazar un campo que sea la diferencia entre Phi en el primer archivo VTK y Phi en el segundo archivo VTK.

Sé cómo hacer esto manualmente; Puedo analizar el texto sin formato en los dos archivos VTK, copiar los datos en matrices, restar una matriz de la otra y luego volcar los datos en el formato correcto en un nuevo archivo. ¿Hay una mejor manera de calcular esta diferencia y exportarla a VTK? Una solución en Python, o en software de visualización como VisIt o Paraview sería preferible a usar un lenguaje compilado como C ++.

El propósito de calcular esta diferencia es comparar diferentes métodos numéricos para calcular la solución de un PDE; Como estoy usando el mismo software para generar las soluciones, puedo garantizar que todos los datos, excepto el campo Phi, serán los mismos en cada archivo que genere.

La forma más fácil que podría encontrar para restar dos campos de diferentes archivos VTK con la misma cuadrícula estructurada es usar un filtro programable en Paraview, que le permite manipular datos utilizando scripts de Python.

En el cuadro de diálogo de filtro programable, puede restar las dos matrices y escribir en la salida con el código:

   phi_0 = inputs[0].CellData['Phi']
   phi_1 = inputs[1].CellData['Phi']
   output.CellData.append(phi_1 - phi_0, 'difference')

En este caso, el campo Phi pasa a ser datos de celda. Si su campo son datos de puntos, reemplace CellDataen todas partes en el script con PointData. Consulte http://public.kitware.com/pipermail/paraview/2010-April/016667.html para obtener más detalles.

En ParaView existe el filtro Atributos de Anexos que puede usarse para esto. Requiere que haya el mismo número de puntos en el conjunto de datos para agregar datos de puntos correctamente y el mismo número de celdas en el conjunto de datos para agregar datos de celdas correctamente. Sin embargo, tendrá problemas con las matrices del mismo nombre (es decir, Phi en su ejemplo). Sin embargo, puede copiar esa matriz fácilmente con el filtro Calculadora antes de usar el filtro Atributos Anexar. Luego puede usar otro filtro de Calculadora para hacer la resta. Sin embargo, esto es probablemente menos eficiente que usar el filtro programable Python de ParaView. Además de eso, podría usar el ejecutable vtkpython para hacerlo manualmente, ya que tendría acceso directo a las cuadrículas y sus atributos.

No tengo un enfoque particularmente bueno, pero copiaría el campo 'phi' de un archivo VTK al otro y lo llamaría 'phiprime' o algo así. Tanto en Paraview como en Visit, tiene la opción de definir nuevos campos mediante una fórmula que utiliza los valores de otros campos. Luego puede definir un "error" de campo como "error = phip-phiprime" en el editor de campo, y trazar este "error" de campo como una superficie, un diagrama de contorno o lo que le interese.

El paso de copiar el bloque de datos de un archivo a otro es claramente incómodo, pero es lo mejor que se me ocurre.

Me di cuenta de que esto es un poco más antiguo, pero pensé que podría estar interesado en la solución VisIt:

Puede hacer esto en VisIt con algo llamado Expresión de campo de malla cruzada basada en conectividad. Es un bocado, pero es básicamente una maquinaria para mapear campos entre bases de datos (en su caso, archivos VTK).

El "Basado en conectividad" (conn_cmfe) se usa cuando la topología es la misma entre los archivos, como en su caso.

También hay un "Basado en posición" (pos_cmfe) que muestrea entre mallas con diferentes topologías.

Para su caso, abra el primer archivo y use la ventana Expresiones para definir una expresión (MyPhi_Diff):

Phi - conn_cmfe(<file2.vtk:Phi>, mesh)

Luego puede trazar "MyPhi_Diff" con un diagrama de Pseudocolor.

También hay un asistente que puede usar para ayudar a definir la expresión

(Menú de opciones -> "Comparaciones de nivel de datos")

Aquí hay más información:

http://visitusers.org/index.php?title=Cmfe