dbeurle / neon

A finite element code
Other
12 stars 8 forks source link

Python script for testing output #73

Closed dbeurle closed 6 years ago

dbeurle commented 6 years ago

Complete integration tests would be useful in ensuring the results are actually correct. This could be done using a python script that checks against analytical solutions.

dbeurle commented 6 years ago

Something similar to


# Script taken from 
# https://stackoverflow.com/questions/21630987/probing-sampling-interpolating-vtk-data-using-python-tvtk-or-mayavi

import sys

import numpy as np
import matplotlib.pyplot as plt
import vtk
from vtk.util.numpy_support import vtk_to_numpy

def read_vtk_file(filename):
    # Read the vtk file with an unstructured grid
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(filename)
    reader.Update()
    return reader

def create_line(p1, p2, sample_points):
    # Create the line along which you want to sample
    line = vtk.vtkLineSource()
    line.SetResolution(sample_points)
    line.SetPoint1(p1)
    line.SetPoint2(p2)
    line.Update()
    return line

def interpolate_over_line(line, reader):

    # Interpolate the data from the VTK-file on the created line.

    # vtkProbeFilter, the probe line is the input, and the underlying dataset
    # is the source.
    probe = vtk.vtkProbeFilter()
    probe.SetInputConnection(line.GetOutputPort())
    probe.SetSourceData(reader.GetOutput())
    probe.Update()

    # Get the data from the VTK-object (probe) to an numpy array
    q = vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('Displacements'))

    samples_on_line = probe.GetOutput().GetNumberOfPoints() 

    # Initialise the points on the line    
    x = np.zeros(samples_on_line)
    y = np.zeros(samples_on_line)
    z = np.zeros(samples_on_line)
    points = np.zeros((samples_on_line , 3))

    # Get the coordinates of the points on the line
    for i in range(samples_on_line):
        x[i], y[i], z[i] = probe.GetOutput().GetPoint(i)
        points[i, 0] = x[i]
        points[i, 1] = y[i]
        points[i, 2] = z[i]
    return points,q

#
#
#

if len(sys.argv) < 3:
    print("A file name and a list of fields must be provided")
    print('Number of arguments:', len(sys.argv), 'arguments.')
    print('Argument List:', sys.argv[1::])
    quit()

file_name = sys.argv[1]
field_requests = sys.argv[2::]

print(field_requests)

start_point = [0.0, 0.0, 10.0]
end_point = [100.0, 0.0, 10.0]

interpolation_points = 10

reader = read_vtk_file(file_name)
line = create_line(start_point, end_point, interpolation_points)

if "Displacements" in field_requests:

    points, results = interpolate_over_line(line, reader)

    results[results == 0] = np.nan 

    # Plot the data over the coordinates and the numerical results
    plt.figure()
    plt.plot(points[:,0], results[:,:])
    plt.title('Displacement over line')
    plt.xlabel('')
    plt.ylabel('')
    plt.legend(['x', 'y', 'z'])
    plt.show()
dbeurle commented 6 years ago

We can also use https://crascit.com/2016/10/18/test-fixtures-with-cmake-ctest/ to add in the order of running using the freshly build executable, then the python script can perform the checks.

dbeurle commented 6 years ago

Implemented for one test case. The rest will follow over time.