test_cellOnVertex_winding.py

#!/usr/bin/python
"""
    Module to test cell properties

    Phillip Wolfram
    LANL
    03/26/2014
"""

# IMPORT PACKAGES ###############
import sys
import shutil

# Need to edit this to include the path to the netcdf4 site packages for
# your python version
sys.path.append('PATH_TO_NETCDF4')

import netCDF4
import numpy as np


# DEFINE FUNCTIONS ###############

def append_cell_winding(netcdfDB):
    """
    Compute and append cell winding information
    to the NETCDF database netcdfDB
    Phillip Wolfram 03/26/2014
    """
    # build new variable

    cellOnVertexWinding = initialize_new_variable(netcdfDB, 'cellOnVertexWinding','i',('nVertices'))

    compute_winding(cellOnVertexWinding, netcdfDB)

def compute_winding(cellOnVertexWinding, netcdfDB):
    """
    compute the winding angle and store sign in cellOnVertexWinding
    Phillip Wolfram 03/26/2014
    """

    # get the datastructures
    nVertices = len(netcdfDB.dimensions['nVertices'])
    cellsOnVertex = netcdfDB.variables['cellsOnVertex']
    xCell = netcdfDB.variables['xCell']
    yCell = netcdfDB.variables['yCell']
    zCell = netcdfDB.variables['zCell']

    for avertex in np.arange(0,nVertices):
        triCells = cellsOnVertex[avertex,:]-1
        if(any(triCells < 0)):
            cellOnVertexWinding[avertex] = 0
        else:
            x = xCell[triCells]
            y = yCell[triCells]
            z = zCell[triCells]
            points = np.vstack((x,y,z))

            AB = points[:,1] - points[:,0]
            BC = points[:,2] - points[:,1]

            pos = np.mean(points,axis=1)

            if(np.dot(np.cross(AB,BC),pos) > 0.0):
                cellOnVertexWinding[avertex] = 1
            else:
                cellOnVertexWinding[avertex] = -1


def initialize_new_variable(netcdfDB, name, datatype, dimsize):
    """
    create a new variable in grid file, name is name of variable,
    datatype (e.g., 'f8'), dimsize, (e.g., 'nVertex')
    Phillip Wolfram 03/26/2014
    """
    try:
        # make new variable
        variable = netcdfDB.createVariable(name, datatype, dimsize)
    except:
        # dimension already exisits
        variable = netcdfDB.variables[name]
    return variable


def save_additions(f_inname, f_outname, append_fields):
    """
    append_fields is a function pointer which appends the
    proper fields to the output database (f_outname) after operating on
    a copied (f_inname) database
    Phillip Wolfram 03/26/2014
    """
    # Read fields from netcdf file
    f_in = netCDF4.Dataset(f_inname, 'r')

    # build the new fields to populate with the returned pointers
    # copy file and open it to append
    shutil.copyfile(options.filename, f_outname)
    f_out = netCDF4.Dataset(f_outname, 'a')

    # add the data onto the new fields
    append_fields(f_out)

    # particle data appended to netCDF file
    f_out.close()
    print 'Finished appending data to ', f_outname, ' file'


# running module as a script...
# TESTING SCRIPT ###############
if __name__ == "__main__":
    from optparse import OptionParser

    # Get command line parameters
    parser = OptionParser()
    parser.add_option("-f", "--file", dest="filename",
                      help="file to open for appending \
                      particle data 'particle_' extension",
                      metavar="FILE")
    options, args = parser.parse_args()

    if not options.filename:
        parser.error("Filename is a required input.")

    save_additions(options.filename, 'winding_' + options.filename, append_cell_winding)