def computeFractionalAnisotropy(
        farray_e1,
        farray_e2,
        farray_e3,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeFractionalAnisotropy ***")

    n_tuples = farray_e1.GetNumberOfTuples()

    farray_FA   = myVTK.createFloatArray("FA"   , 1, n_tuples)
    farray_FA12 = myVTK.createFloatArray("FA_12", 1, n_tuples)
    farray_FA23 = myVTK.createFloatArray("FA_23", 1, n_tuples)

    for k_tuple in xrange(n_tuples):
        e1 = farray_e1.GetTuple1(k_tuple)
        e2 = farray_e2.GetTuple1(k_tuple)
        e3 = farray_e3.GetTuple1(k_tuple)
        FA   = ((e1-e2)**2+(e1-e3)**2+(e2-e3)**2)**(0.5) / (2*(e1**2+e2**2+e3**2))**(0.5)
        FA12 = ((e1-e2)**2)**(0.5) / (e1**2+e2**2)**(0.5)
        FA23 = ((e2-e3)**2)**(0.5) / (e2**2+e3**2)**(0.5)

        farray_FA.SetTuple1(k_tuple, FA)
        farray_FA12.SetTuple1(k_tuple, FA12)
        farray_FA23.SetTuple1(k_tuple, FA23)

    return (farray_FA,
            farray_FA12,
            farray_FA23)

def readAbaqusFibersFromINP(
        filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readAbaqusFibersFromINP: " + filename + " ***")

    eF_array = myVTK.createFloatArray('eF', 3)
    eS_array = myVTK.createFloatArray('eS', 3)
    eN_array = myVTK.createFloatArray('eN', 3)

    file = open(filename, 'r')
    file.readline()

    for line in file:
        line = line.split(', ')
        #print line

        eF = [float(item) for item in line[1:4]]
        eS = [float(item) for item in line[4:7]]
        eN = numpy.cross(eF,eS)
        #print "eF =", eF
        #print "eS =", eS
        #print "eN =", eN

        eF_array.InsertNextTuple(eF)
        eS_array.InsertNextTuple(eS)
        eN_array.InsertNextTuple(eN)

    file.close()

    myVTK.myPrint(verbose, "n_tuples = " + str(eF_array.GetNumberOfTuples()))

    return (eF_array,
            eS_array,
            eN_array)

def computeHelixTransverseSheetAngles(
        farray_eRR,
        farray_eCC,
        farray_eLL,
        farray_eF,
        farray_eS,
        farray_eN,
        use_new_definition=False,
        verbose=0):

    myVTK.myPrint(verbose, "*** computeHelixTransverseSheetAngles ***")

    n_tuples = farray_eRR.GetNumberOfTuples()
    assert (farray_eCC.GetNumberOfTuples() == n_tuples)
    assert (farray_eLL.GetNumberOfTuples() == n_tuples)
    assert (farray_eF.GetNumberOfTuples() == n_tuples)
    assert (farray_eS.GetNumberOfTuples() == n_tuples)
    assert (farray_eN.GetNumberOfTuples() == n_tuples)

    farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_tuples)
    farray_angle_trans = myVTK.createFloatArray("angle_trans", 1, n_tuples)
    farray_angle_sheet = myVTK.createFloatArray("angle_sheet", 1, n_tuples)

    eRR = numpy.empty(3)
    eCC = numpy.empty(3)
    eLL = numpy.empty(3)
    eF = numpy.empty(3)
    for k_tuple in xrange(n_tuples):
        farray_eRR.GetTuple(k_tuple, eRR)
        farray_eCC.GetTuple(k_tuple, eCC)
        farray_eLL.GetTuple(k_tuple, eLL)

        farray_eF.GetTuple(k_tuple, eF)
        eF -= numpy.dot(eF, eRR) * eRR
        eF /= numpy.linalg.norm(eF)
        angle_helix = math.copysign(1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eLL)))) * (180./math.pi)
        farray_angle_helix.SetTuple1(k_tuple, angle_helix)

        eF  = numpy.array(farray_eF.GetTuple(k_tuple))
        eF -= numpy.dot(eF, eLL) * eLL
        eF /= numpy.linalg.norm(eF)
        angle_trans = math.copysign(-1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eRR)))) * (180./math.pi)
        farray_angle_trans.SetTuple1(k_tuple, angle_trans)

        #if (use_new_definition):
            #assert 0, "TODO"
        #else:
            #assert 0, "TODO"

    return (farray_angle_helix,
            farray_angle_trans,
            farray_angle_sheet)

def createArray(
        array_name,
        n_components=1,
        n_tuples=0,
        array_type="float",
        init_to_zero=0,
        verbose=1):

    assert (type(array_type) in [type, str]), "array_type must be a type or a str. Aborting."

    if (type(array_type) is type):
        assert (array_type in [float, int]), "if a type, array_type must be equal to float or int. Aborting."

        if (array_type == float):
            return myVTK.createFloatArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == int):
            return myVTK.createIntArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
    elif (type(array_type) is str):
        assert (array_type in ["double", "float", "int", "short"]), "if a str, array_type must be equal to double, float, int or short. Aborting."

        if (array_type == "float") or (array_type == "double"):
            return myVTK.createFloatArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == "int"):
            return myVTK.createIntArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == "short"):
            return myVTK.createShortArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)

def readAbaqusStressFromDAT(
        data_filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readAbaqusStressFromDAT: " + data_filename + " ***")

    s_array = myVTK.createFloatArray("", 6)

    data_file = open(data_filename, 'r')
    context = ""
    k_cell = 0
    for line in data_file:
        if (context == "reading stresses"):
            #print line
            if ("MAXIMUM" in line):
                context = ""
                continue
            if ("OR" in line):
                splitted_line = line.split()
                assert (int(splitted_line[0]) == k_cell+1), "Wrong element number. Aborting."
                s_list = [float(splitted_line[k]) for k in xrange(3,9)]
                s_array.InsertNextTuple(s_list)
                k_cell += 1

        if (line == "    ELEMENT  PT FOOT-       S11         S22         S33         S12         S13         S23     \n"):
            context = "reading stresses"

    data_file.close()

    myVTK.myPrint(verbose, "n_tuples = " + str(s_array.GetNumberOfTuples()))

    return s_array

def readDynaDeformationGradients(
        mesh,
        hystory_files_basename,
        array_name,
        verbose=1):

    myVTK.myPrint(verbose, "*** readDynaDeformationGradients ***")

    n_cells = mesh.GetNumberOfCells()

    history_files_names = [hystory_files_basename + '.history#' + str(num) for num in xrange(11,20)]

    F_list = [[0. for k_component in xrange(9)] for k_cell in xrange(n_cells)]

    for k_component in xrange(9):
        history_file = open(history_files_names[k_component], 'r')
        for line in history_file:
            if line.startswith('*') or line.startswith('$'): continue
            line = line.split()
            F_list[int(line[0])-1][k_component] = float(line[1])
        history_file.close()

    F_array = myVTK.createFloatArray(array_name, 9, n_cells)

    for k_cell in xrange(n_cells):
        F_array.InsertTuple(k_cell, F_list[k_cell])

    myVTK.myPrint(verbose, "n_tuples = " + str(F_array.GetNumberOfTuples()))

    mesh.GetCellData().AddArray(F_array)

def readAbaqusDeformationGradientsFromDAT(
        data_filename,
        verbose=0):

    myVTK.myPrint(verbose, "*** readAbaqusDeformationGradientsFromDAT: "+data_filename+" ***")

    farray_F = myVTK.createFloatArray("F", 9)

    data_file = open(data_filename, 'r')
    context = ""
    k_cell = 0
    for line in data_file:
        if (context == "reading deformation gradients"):
            #print line
            if ("MAXIMUM" in line):
                context = ""
                continue
            if ("OR" in line):
                splitted_line = line.split()
                assert (int(splitted_line[0]) == k_cell+1), "Wrong element number. Aborting."
                F_list = [float(splitted_line[ 3]), float(splitted_line[ 6]), float(splitted_line[7]),
                          float(splitted_line[ 9]), float(splitted_line[ 4]), float(splitted_line[8]),
                          float(splitted_line[10]), float(splitted_line[11]), float(splitted_line[5])]
                farray_F.InsertNextTuple(F_list)
                k_cell += 1

        if (line == "    ELEMENT  PT FOOT-       DG11        DG22        DG33        DG12        DG13        DG23        DG21        DG31        DG32    \n"):
            context = "reading deformation gradients"

    data_file.close()

    myVTK.myPrint(verbose-1, "n_tuples = "+str(farray_F.GetNumberOfTuples()))

    return farray_F

def computeHelixAngles(
        farray_eRR,
        farray_eCC,
        farray_eLL,
        farray_eF,
        verbose=0):

    myVTK.myPrint(verbose, "*** computeHelixAngles ***")

    n_tuples = farray_eRR.GetNumberOfTuples()
    assert (farray_eCC.GetNumberOfTuples() == n_tuples)
    assert (farray_eLL.GetNumberOfTuples() == n_tuples)
    assert (farray_eF.GetNumberOfTuples() == n_tuples)

    farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_tuples)

    eRR = numpy.empty(3)
    eCC = numpy.empty(3)
    eLL = numpy.empty(3)
    eF = numpy.empty(3)
    for k_tuple in xrange(n_tuples):
        farray_eRR.GetTuple(k_tuple, eRR)
        farray_eCC.GetTuple(k_tuple, eCC)
        farray_eLL.GetTuple(k_tuple, eLL)
        farray_eF.GetTuple(k_tuple, eF)
        eF -= numpy.dot(eF, eRR) * eRR
        eF /= numpy.linalg.norm(eF)
        helix_angle = math.copysign(1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eLL)))) * (180./math.pi)
        farray_angle_helix.SetTuple1(k_tuple, helix_angle)

    return farray_angle_helix

def computeSyntheticHelixAngles(
        farray_rr,
        helix_angle_end,
        helix_angle_epi,
        farray_angle_helix=None,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeSyntheticHelixAngles ***")

    n_cells = farray_rr.GetNumberOfTuples()

    if (farray_angle_helix is None):
        farray_angle_helix = myVTK.createFloatArray(
            name="angle_helix",
            n_components=1,
            n_tuples=n_cells)

    for k_cell in xrange(n_cells):
        rr = farray_rr.GetTuple1(k_cell)

        helix_angle_in_degrees = (1.-rr) * helix_angle_end \
                               +     rr  * helix_angle_epi
        farray_angle_helix.SetTuple1(
            k_cell,
            helix_angle_in_degrees)

    return farray_angle_helix

def computeSystolicStrainsFromEndDiastolicAndEndSystolicStates(
        farray_F_dia,
        farray_F_sys,
        verbose=0):

    myVTK.myPrint(verbose, "*** computeSystolicStrainsFromEndDiastolicAndEndSystolicStates ***")

    n_tuples = farray_F_dia.GetNumberOfTuples()
    assert (farray_F_sys.GetNumberOfTuples() == n_tuples)

    farray_E_dia = myVTK.createFloatArray('E_dia', 6, n_tuples)
    farray_E_sys = myVTK.createFloatArray('E_sys', 6, n_tuples)
    farray_F_num = myVTK.createFloatArray('F_num', 6, n_tuples)
    farray_E_num = myVTK.createFloatArray('E_num', 6, n_tuples)

    I = numpy.eye(3)
    E_vec = numpy.empty(6)
    for k_tuple in xrange(n_tuples):
        F_dia = numpy.reshape(farray_F_dia.GetTuple(k_tuple), (3,3), order='C')
        F_sys = numpy.reshape(farray_F_sys.GetTuple(k_tuple), (3,3), order='C')
        #print 'F_dia =', F_dia
        #print 'F_sys =', F_sys

        C = numpy.dot(numpy.transpose(F_dia), F_dia)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_dia.SetTuple(k_tuple, E_vec)

        C = numpy.dot(numpy.transpose(F_sys), F_sys)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_sys.SetTuple(k_tuple, E_vec)

        F = numpy.dot(F_sys, numpy.linalg.inv(F_dia))
        farray_F_num.SetTuple(k_tuple, numpy.reshape(F, 9, order='C'))
        #print 'F =', F

        C = numpy.dot(numpy.transpose(F), F)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_num.SetTuple(k_tuple, E_vec)

    return (farray_E_dia,
            farray_E_sys,
            farray_F_num,
            farray_E_num)

def generateWarpedMesh(
        mesh_folder,
        mesh_basename,
        images,
        structure,
        deformation,
        evolution,
        verbose=0):

    myVTK.myPrint(verbose, "*** generateWarpedMesh ***")

    mesh = myVTK.readUGrid(
        filename=mesh_folder+"/"+mesh_basename+".vtk",
        verbose=verbose-1)
    n_points = mesh.GetNumberOfPoints()
    n_cells = mesh.GetNumberOfCells()

    if os.path.exists(mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk"):
        ref_mesh = myVTK.readUGrid(
            filename=mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk",
            verbose=verbose-1)
    else:
        ref_mesh = None

    farray_disp = myVTK.createFloatArray(
        name="displacement",
        n_components=3,
        n_tuples=n_points,
        verbose=verbose-1)
    mesh.GetPointData().AddArray(farray_disp)

    mapping = Mapping(images, structure, deformation, evolution)

    X = numpy.empty(3)
    x = numpy.empty(3)
    U = numpy.empty(3)
    if ("zfill" not in images.keys()):
        images["zfill"] = len(str(images["n_frames"]))
    for k_frame in xrange(images["n_frames"]):
        t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
        mapping.init_t(t)

        for k_point in xrange(n_points):
            mesh.GetPoint(k_point, X)
            mapping.x(X, x)
            U = x - X
            farray_disp.SetTuple(k_point, U)

        myVTK.computeStrainsFromDisplacements(
            mesh=mesh,
            disp_array_name="displacement",
            ref_mesh=ref_mesh,
            verbose=verbose-1)

        myVTK.writeUGrid(
            ugrid=mesh,
            filename=mesh_folder+"/"+mesh_basename+"_"+str(k_frame).zfill(images["zfill"])+".vtk",
            verbose=verbose-1)

def computeCartesianCoordinates(
        points,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeCartesianCoordinates ***")

    n_points = points.GetNumberOfPoints()

    [xmin, xmax, ymin, ymax, zmin, zmax] = points.GetBounds()
    dx = xmax-xmin
    dy = ymax-ymin
    dz = zmax-zmin
    if (verbose >= 2): print "xmin = "+str(xmin)
    if (verbose >= 2): print "xmax = "+str(xmax)
    if (verbose >= 2): print "dx = "+str(dx)
    if (verbose >= 2): print "ymin = "+str(ymin)
    if (verbose >= 2): print "ymax = "+str(ymax)
    if (verbose >= 2): print "dy = "+str(dy)
    if (verbose >= 2): print "zmin = "+str(zmin)
    if (verbose >= 2): print "zmax = "+str(zmax)
    if (verbose >= 2): print "dz = "+str(dz)

    farray_xx = myVTK.createFloatArray("xx", 1, n_points)
    farray_yy = myVTK.createFloatArray("yy", 1, n_points)
    farray_zz = myVTK.createFloatArray("zz", 1, n_points)

    point = numpy.empty(3)
    for k_point in xrange(n_points):
        if (verbose >= 2): print "k_point = "+str(k_point)

        points.GetPoint(k_point, point)
        if (verbose >= 2): print "point = "+str(point)

        xx = (point[0] - xmin) / dx
        yy = (point[1] - ymin) / dy
        zz = (point[2] - zmin) / dz

        farray_xx.SetTuple1(k_point, xx)
        farray_yy.SetTuple1(k_point, yy)
        farray_zz.SetTuple1(k_point, zz)

    return (farray_xx,
            farray_yy,
            farray_zz)

def computeCartesianCoordinates(
        points,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeCartesianCoordinates ***")

    n_points = points.GetNumberOfPoints()

    [xmin, xmax, ymin, ymax, zmin, zmax] = points.GetBounds()
    dx = xmax-xmin
    dy = ymax-ymin
    dz = zmax-zmin
    if (verbose >= 2): print "xmin = " + str(xmin)
    if (verbose >= 2): print "xmax = " + str(xmax)
    if (verbose >= 2): print "dx = " + str(dx)
    if (verbose >= 2): print "ymin = " + str(ymin)
    if (verbose >= 2): print "ymax = " + str(ymax)
    if (verbose >= 2): print "dy = " + str(dy)
    if (verbose >= 2): print "zmin = " + str(zmin)
    if (verbose >= 2): print "zmax = " + str(zmax)
    if (verbose >= 2): print "dz = " + str(dz)

    farray_norm_x = myVTK.createFloatArray("norm_x", 1, n_points)
    farray_norm_y = myVTK.createFloatArray("norm_y", 1, n_points)
    farray_norm_z = myVTK.createFloatArray("norm_z", 1, n_points)

    for k_point in xrange(n_points):
        if (verbose >= 2): print "k_point = " + str(k_point)

        point = points.GetPoint(k_point)
        if (verbose >= 2): print "point = " + str(point)

        norm_x = (point[0] - xmin) / dx
        norm_y = (point[1] - ymin) / dy
        norm_z = (point[2] - zmin) / dz

        farray_norm_x.InsertTuple(k_point, [norm_x])
        farray_norm_y.InsertTuple(k_point, [norm_y])
        farray_norm_z.InsertTuple(k_point, [norm_z])

    return (farray_norm_x,
            farray_norm_y,
            farray_norm_z)

def computeSystolicStrains(
        farray_F_dia,
        farray_F_sys,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeSystolicStrains ***")

    n_tuples = farray_F_dia.GetNumberOfTuples()

    farray_E_dia = myVTK.createFloatArray('E_dia', 6, n_tuples)
    farray_E_sys = myVTK.createFloatArray('E_sys', 6, n_tuples)
    farray_F_num = myVTK.createFloatArray('F_num', 6, n_tuples)
    farray_E_num = myVTK.createFloatArray('E_num', 6, n_tuples)

    for k_tuple in xrange(n_tuples):
        F_dia = numpy.reshape(farray_F_dia.GetTuple(k_tuple), (3,3), order='C')
        F_sys = numpy.reshape(farray_F_sys.GetTuple(k_tuple), (3,3), order='C')
        #print 'F_dia =', F_dia
        #print 'F_sys =', F_sys

        C = numpy.dot(numpy.transpose(F_dia), F_dia)
        E = (C - numpy.eye(3))/2
        farray_E_dia.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

        C = numpy.dot(numpy.transpose(F_sys), F_sys)
        E = (C - numpy.eye(3))/2
        farray_E_sys.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

        F = numpy.dot(F_sys, numpy.linalg.inv(F_dia))
        farray_F_num.InsertTuple(k_tuple, numpy.reshape(F, 9, order='C'))
        #print 'F =', F

        C = numpy.dot(numpy.transpose(F), F)
        E = (C - numpy.eye(3))/2
        farray_E_num.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

    return (farray_E_dia,
            farray_E_sys,
            farray_F_num,
            farray_E_num)

def rotateSymmetricMatrix(
        old_array,
        old_array_storage="vec",
        in_vecs=None,
        out_vecs=None,
        verbose=1):

    myVTK.myPrint(verbose, "*** rotateSymmetricMatrix ***")

    n_tuples = old_array.GetNumberOfTuples()
    new_array = myVTK.createFloatArray("", 6, n_tuples)

    for k_tuple in xrange(n_tuples):
        old_matrix = old_array.GetTuple(k_tuple)
        if (old_array_storage == "vec"):
            old_matrix = vec_col_to_mat_sym(old_matrix)
        elif (old_array_storage == "Cmat"):
            old_matrix = numpy.reshape(old_matrix, (3,3), order='C')
        elif (old_array_storage == "Fmat"):
            old_matrix = numpy.reshape(old_matrix, (3,3), order='F')

        if (in_vecs == None):
            in_R = numpy.eye(3)
        else:
            in_R = numpy.transpose(numpy.array([in_vecs[0].GetTuple(k_tuple),
                                                in_vecs[1].GetTuple(k_tuple),
                                                in_vecs[2].GetTuple(k_tuple)]))

        if (out_vecs == None):
            out_R = numpy.eye(3)
        else:
            out_R = numpy.transpose(numpy.array([out_vecs[0].GetTuple(k_tuple),
                                                 out_vecs[1].GetTuple(k_tuple),
                                                 out_vecs[2].GetTuple(k_tuple)]))

        R = numpy.dot(numpy.transpose(in_R), out_R)

        new_matrix = numpy.dot(numpy.dot(numpy.transpose(R), old_matrix), R)

        if (old_array_storage == "vec"):
            new_matrix = mat_sym_to_vec_col(new_matrix)
        elif (old_array_storage == "Cmat"):
            new_matrix = numpy.reshape(new_matrix, 9, order='C')
        elif (old_array_storage == "Fmat"):
            new_matrix = numpy.reshape(new_matrix, 9, order='F')

        new_array.InsertTuple(k_tuple, new_matrix)

    return new_array

def maskImageWithMesh(
    image,
    mesh,
    filter_with_field=None,
    verbose=0):

    myVTK.myPrint(verbose, "*** maskImageWithMesh ***")

    n_points = image.GetNumberOfPoints()
    farray_scalars_image = image.GetPointData().GetArray("scalars") # note that the field is defined at the points, not the cells

    (cell_locator,
     closest_point,
     generic_cell,
     k_cell,
     subId,
     dist) = myVTK.getCellLocator(
        mesh=mesh,
        verbose=verbose-1)

    if (filter_with_field != None):
        field = mesh.GetCellData().GetArray(filter_with_field[0])
        field_values = filter_with_field[1]

    points = vtk.vtkPoints()

    farray_scalars = myVTK.createFloatArray(
        name="scalars",
        n_components=1)

    point = numpy.empty(3)
    for k_point in xrange(n_points):
        image.GetPoint(k_point, point)

        k_cell = cell_locator.FindCell(point)
        if (k_cell == -1): continue
        if (filter_with_field != None) and (field.GetTuple1(k_cell) not in field_values): continue

        points.InsertNextPoint(point)
        farray_scalars.InsertNextTuple(farray_scalars_image.GetTuple(k_point))

    ugrid = vtk.vtkUnstructuredGrid()
    ugrid.SetPoints(points)
    ugrid.GetPointData().AddArray(farray_scalars)
    myVTK.addVertices(
        ugrid)

    return ugrid

def computeSyntheticHelixAngles(farray_rr, helix_angle_end, helix_angle_epi, farray_angle_helix=None, verbose=1):

    myVTK.myPrint(verbose, "*** computeSyntheticHelixAngles ***")

    n_cells = farray_rr.GetNumberOfTuples()

    if farray_angle_helix == None:
        farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_cells)

    for k_cell in xrange(n_cells):
        rr = farray_rr.GetTuple(k_cell)[0]

        helix_angle_in_degrees = (1.0 - rr) * helix_angle_end + rr * helix_angle_epi
        farray_angle_helix.InsertTuple(k_cell, [helix_angle_in_degrees])

    return farray_angle_helix

def addDeformationGradients(
        mesh,
        disp_array_name="Displacement",
        defo_grad_array_name="DeformationGradient",
        verbose=0):

    mypy.my_print(verbose, "*** addDeformationGradients ***")

    mypy.my_print(min(verbose,1), "*** Warning: at some point the ordering of vector gradient components has changed, and uses C ordering instead of F. ***")
    if   (vtk.vtkVersion.GetVTKMajorVersion() >= 8):
        ordering = "C"
    elif (vtk.vtkVersion.GetVTKMajorVersion() == 7) and ((vtk.vtkVersion.GetVTKMinorVersion() > 0) or (vtk.vtkVersion.GetVTKBuildVersion() > 0)):
        ordering = "C"
    else:
        ordering = "F"

    n_cells = mesh.GetNumberOfCells()

    assert (mesh.GetPointData().HasArray(disp_array_name))
    mesh.GetPointData().SetActiveVectors(disp_array_name)
    cell_derivatives = vtk.vtkCellDerivatives()
    cell_derivatives.SetVectorModeToPassVectors()
    cell_derivatives.SetTensorModeToComputeGradient()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        cell_derivatives.SetInputData(mesh)
    else:
        cell_derivatives.SetInput(mesh)
    cell_derivatives.Update()
    farray_gu = cell_derivatives.GetOutput().GetCellData().GetArray("VectorGradient")

    farray_defo_grad = myvtk.createFloatArray(
        name=defo_grad_array_name,
        n_components=9,
        n_tuples=n_cells)
    mesh.GetCellData().AddArray(farray_defo_grad)
    I = numpy.eye(3)
    for k_cell in range(n_cells):
        GU = numpy.reshape(farray_gu.GetTuple(k_cell), (3,3), order=ordering)
        F = I + GU
        farray_defo_grad.SetTuple(k_cell, numpy.reshape(F, 9, order='C'))

def rotateTensors(
        old_tensors,
        R,
        storage="vec",
        verbose=1):

    myVTK.myPrint(verbose, "*** rotateTensors ***")

    assert (storage in ["vec", "Cmat", "Fmat"]), "\"storage\" must be \"vec\", \"Cmat\" or \"Fmat\". Aborting."

    n_tuples = old_tensors.GetNumberOfTuples()
    n_components = old_tensors.GetNumberOfComponents()

    new_tensors = myVTK.createFloatArray("tensors", n_components, n_tuples)

    for k_tuple in xrange(n_tuples):
        old_tensor = old_tensors.GetTuple(k_tuple)
        if (storage == "vec"):
            old_tensor = vec_col_to_mat_sym(old_tensor)
        elif (storage == "Cmat"):
            old_tensor = numpy.reshape(old_tensor, (3,3), order='C')
        elif (storage == "Fmat"):
            old_tensor = numpy.reshape(old_tensor, (3,3), order='F')

        new_tensor = numpy.dot(R, numpy.dot(old_tensor, numpy.transpose(R)))

        if (storage == "vec"):
            new_tensor = mat_sym_to_vec_col(new_tensor)
        elif (storage == "Cmat"):
            new_tensor = numpy.reshape(new_tensor, 9, order='C')
        elif (storage == "Fmat"):
            new_tensor = numpy.reshape(new_tensor, 9, order='F')

        new_tensors.InsertTuple(k_tuple, new_tensor)

    return new_tensors

def computeHelixTransverseSheetAngles2(
        farray_eRR,
        farray_eCC,
        farray_eLL,
        farray_eF,
        farray_eS,
        farray_eN,
        use_new_definition=False,
        ref_vectors_are_material_basis=False,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeHelixTransverseSheetAngles2 ***")

    n_tuples = farray_eRR.GetNumberOfTuples()
    assert (farray_eCC.GetNumberOfTuples() == n_tuples)
    assert (farray_eLL.GetNumberOfTuples() == n_tuples)
    assert (farray_eF.GetNumberOfTuples() == n_tuples)
    assert (farray_eS.GetNumberOfTuples() == n_tuples)
    assert (farray_eN.GetNumberOfTuples() == n_tuples)

    farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_tuples)
    farray_angle_trans = myVTK.createFloatArray("angle_trans", 1, n_tuples)
    farray_angle_sheet = myVTK.createFloatArray("angle_sheet", 1, n_tuples)

    for k_tuple in xrange(n_tuples):
        #print "k_tuple = " + str(k_tuple)

        eRR = numpy.array(farray_eRR.GetTuple(k_tuple))
        eCC = numpy.array(farray_eCC.GetTuple(k_tuple))
        eLL = numpy.array(farray_eLL.GetTuple(k_tuple))
        eF  = numpy.array(farray_eF.GetTuple(k_tuple))
        eS  = numpy.array(farray_eS.GetTuple(k_tuple))
        eN  = numpy.array(farray_eN.GetTuple(k_tuple))

        #print "eRR = " + str(eRR)
        #print "eCC = " + str(eCC)
        #print "eLL = " + str(eLL)
        #print "eF = " + str(eF)
        #print "eS = " + str(eS)
        #print "eN = " + str(eN)

        #print "|eRR| = " + str(numpy.linalg.norm(eRR))
        #print "|eCC| = " + str(numpy.linalg.norm(eCC))
        #print "|eLL| = " + str(numpy.linalg.norm(eLL))

        #print "eRR.eCC = " + str(numpy.dot(eRR, eCC))
        #print "eCC.eLL = " + str(numpy.dot(eCC, eLL))
        #print "eLL.eRR = " + str(numpy.dot(eLL, eRR))

        #print "|eF| = " + str(numpy.linalg.norm(eF))
        #print "|eS| = " + str(numpy.linalg.norm(eS))
        #print "|eN| = " + str(numpy.linalg.norm(eN))

        #print "eF.eS = " + str(numpy.dot(eF, eS))
        #print "eS.eN = " + str(numpy.dot(eS, eN))
        #print "eN.eF = " + str(numpy.dot(eN, eF))

        assert (round(numpy.linalg.norm(eRR),1) == 1.0),\
            "|eRR| = " + str(numpy.linalg.norm(eRR)) + " ≠ 1. Aborting."
        assert (round(numpy.linalg.norm(eCC),1) == 1.0),\
            "|eCC| = " + str(numpy.linalg.norm(eCC)) + " ≠ 1. Aborting."
        assert (round(numpy.linalg.norm(eLL),1) == 1.0),\
            "|eLL| = " + str(numpy.linalg.norm(eLL)) + " ≠ 1. Aborting."

        assert (round(numpy.dot(eRR,eCC),1) == 0.0),\
            "eRR.eCC = " + str(numpy.dot(eRR,eCC)) + " ≠ 0. Aborting."
        assert (round(numpy.dot(eCC,eLL),1) == 0.0),\
            "eCC.eLL = " + str(numpy.dot(eCC,eLL)) + " ≠ 0. Aborting."
        assert (round(numpy.dot(eLL,eRR),1) == 0.0),\
            "eLL.eRR = " + str(numpy.dot(eLL,eRR)) + " ≠ 0. Aborting."

        assert (round(numpy.linalg.det([eRR, eCC, eLL]),1) == 1.0),\
            "det([eRR, eCC, eLL]) = " + str(numpy.linalg.det([eRR, eCC, eLL])) + " ≠ 1. Aborting."

        assert (round(numpy.linalg.norm(eF),1) == 1.0),\
            "|eF| = " + str(numpy.linalg.norm(eF)) + " ≠ 1. Aborting."
        assert (round(numpy.linalg.norm(eS),1) == 1.0),\
            "|eS| = " + str(numpy.linalg.norm(eS)) + " ≠ 1. Aborting."
        assert (round(numpy.linalg.norm(eN),1) == 1.0),\
            "|eN| = " + str(numpy.linalg.norm(eN)) + " ≠ 1. Aborting."

        assert (round(numpy.dot(eF,eS),1) == 0.0),\
            "eF.eS = " + str(numpy.dot(eF,eS)) + " ≠ 0. Aborting."
        assert (round(numpy.dot(eS,eN),1) == 0.0),\
            "eS.eN = " + str(numpy.dot(eS,eN)) + " ≠ 0. Aborting."
        assert (round(numpy.dot(eN,eF),1) == 0.0),\
            "eN.eF = " + str(numpy.dot(eN,eF)) + " ≠ 0. Aborting."

        assert (round(numpy.linalg.det([eF, eS, eN]),1) == 1.0),\
            "det([eF, eS, eN]) = " + str(numpy.linalg.det([eF, eS, eN])) + " ≠ 1. Aborting."

        # reference basis
        if (ref_vectors_are_material_basis):
            ref = numpy.array([+eRR, +eCC, +eLL])
        else:
            if (use_new_definition):
                ref = numpy.array([+eCC, +eLL, +eRR])
            else:
                ref = numpy.array([+eCC, +eRR, -eLL])

#.........这里部分代码省略.........