expressions/tensor_expr.py
# ----------------------------------------------------------------------------
#  CLASSES: nightly
#
#  Test Case:  tensor_expr.py
#
#  Tests:      tensor expressions using simple, constant valued tensors for
#              which answers are known. The known answers were obtained by
#              scouring the internet for example problems that were worked.
#
#  Mark C. Miller, Mon Nov 11 14:52:02 PST 2019
#
# ----------------------------------------------------------------------------
import math, re

#
# Scans a string for possible python iterables, builds a list of them
# and returns them with their values rounded to specific numbers of
# digits.
#
def ExtractIterablesFromString(s,pair='()',rnd=4):
    retval = []
    for q in s.split(pair[0]):
        for r in q.split(pair[1]):
            if re.match('^[ 0-9eE.,+-]+$', r):
                retval.append([round(float(x),rnd) for x in r.split(',')])
    return retval

def EqualEigVecs(a,b):
    r = []
    for i in range(len(b)):
        if b[i] != 0:
            r.append(float(a[i])/float(b[i]))
    if min(r) == 0:
        return max(r) < pow(10,-Prec+1)
    else:
        return abs((max(r)/min(r))-1) < pow(10,-Prec+1)

#
# Creates a group of related expressions; 9 constant scalar expressions and
# from them a constant tensor expression with specified centering all in a
# sub-menu with of the given name. For 2D meshes, you would still create a
# 9 component (3x3) tensor but the z-dimension values would all be zeros.
#
def CreateConstantTensorExpr(name, meshName, constType, vals):
    comps = ("s11","s12","s13","s21","s22","s23","s31","s32","s33")
    for i in range(len(vals)):
        DefineScalarExpression("%s/%s"%(name,comps[i]),\
            "%s_constant(<%s>, %g)"%(constType,meshName,vals[i]))
    DefineTensorExpression("%s/tensor"%name,
        "{{<%s/s11>,<%s/s12>,<%s/s13>},\
          {<%s/s21>,<%s/s22>,<%s/s23>},\
          {<%s/s31>,<%s/s32>,<%s/s33>}}"%(name, name, name, name,\
              name, name, name, name, name))

#
# Set precision for rounding operations
#
Prec = 5

#  Since we use the expression system to construct mesh-wide constant values,
#  all we need as far as a database is a simple, small mesh. Maybe for both
#  two and three dimensions.
OpenDatabase(silo_data_path("arbpoly.silo"))
meshName = "clipped_hex"

TestSection("2D Tensor Maximum Shear")
CreateConstantTensorExpr("max_shear_2d", meshName, "nodal",\
   (50,  30,   0,\
    30, -20,   0,\
     0,   0,   0))
DefineScalarExpression("max_shear_2d/result", "tensor_maximum_shear(<max_shear_2d/tensor>)")
AddPlot("Pseudocolor", "max_shear_2d/result")
DrawPlots()
Query("MinMax")
q = GetQueryOutputObject()
TestValueEQ("Maximum Shear 2D", q['min'], 46.09772, Prec)
TestValueEQ("Maximum Shear 2D", q['max'], 46.09772, Prec)
DeleteAllPlots()

TestSection("3D Tensor Maximum Shear")
CreateConstantTensorExpr("max_shear", meshName, "nodal",\
   (5,   0,   0,\
    0,  -6, -12,\
    0, -12,   1))
DefineScalarExpression("max_shear/result", "tensor_maximum_shear(<max_shear/tensor>)")
AddPlot("Pseudocolor", "max_shear/result")
DrawPlots()
Query("MinMax")
q = GetQueryOutputObject()
TestValueEQ("Maximum Shear 3D", q['min'], 12.5, Prec)
TestValueEQ("Maximum Shear 3D", q['max'], 12.5, Prec)
DeleteAllPlots()

TestSection("2D Effective Tensor")
CreateConstantTensorExpr("eff_tensor_2d", meshName, "nodal",\
   (50,  30,   0,\
    30, -20,   0,\
     0,   0,   0))
DefineScalarExpression("eff_tensor_2d/result", "effective_tensor(<eff_tensor_2d/tensor>)")
AddPlot("Pseudocolor", "eff_tensor_2d/result")
DrawPlots()
Query("MinMax")
q = GetQueryOutputObject()
TestValueEQ("2D Effective Tensor", q['min'], 81.24039, Prec)
TestValueEQ("2D Effective Tensor", q['max'], 81.24039, Prec)
DeleteAllPlots()

TestSection("3D Effective Tensor")
CreateConstantTensorExpr("eff_tensor", meshName, "nodal",\
   (2,  -3,   4,\
   -3,  -5,   1,\
    4,   1,   6))
DefineScalarExpression("eff_tensor/result", "effective_tensor(<eff_tensor/tensor>)")
AddPlot("Pseudocolor", "eff_tensor/result")
DrawPlots()
Query("MinMax")
q = GetQueryOutputObject()
TestValueEQ("Effective Tensor", q['min'], 13.0767, Prec)
TestValueEQ("Effective Tensor", q['max'], 13.0767, Prec)
DeleteAllPlots()

TestSection("3D, Symmetric Eigenvalues and Eigenvectors")
CreateConstantTensorExpr("eigvals_symm2", meshName, "nodal",\
    (2, -1,  0,\
    -1,  2, -1,\
     0, -1,  2))
DefineVectorExpression("eigvals_symm2/result", "eigenvalue(<eigvals_symm2/tensor>)")
AddPlot("Vector", "eigvals_symm2/result")
DrawPlots()
p = PickByNode(0,('eigvals_symm2/result',))
eigvals = list(p['eigvals_symm2/result'])
TestValueIN("First Eigenvalue of 2", eigvals, 2, Prec)
TestValueIN("Second Eigenvalue of 2+sqrt(2)", eigvals, 2+math.sqrt(2), Prec)
TestValueIN("Third Eigenvalue of 2-sqrt(2)",  eigvals, 2-math.sqrt(2), Prec)
DeleteAllPlots()
DefineTensorExpression("eigvals_symm2/result2", "transpose(eigenvector(<eigvals_symm2/tensor>))")
AddPlot("Tensor", "eigvals_symm2/result2")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
vecs = ExtractIterablesFromString(s, '()', Prec)
TestValueIN("First Eigenvector of (1,0,-1)", vecs, (1,0,-1), Prec, EqualEigVecs)
TestValueIN("Second Eigenvector of (1,-sqrt(2),1)", vecs, (1,-math.sqrt(2),1), Prec, EqualEigVecs)
TestValueIN("Third Eigenvector of (1,sqrt(2),1)", vecs, (1,math.sqrt(2),1), Prec, EqualEigVecs)

TestSection("3D, Symmetric Eigenvalues and Eigenvectors with Repeated values")
CreateConstantTensorExpr("eigvals_symm", meshName, "nodal",\
    (3,2,4,\
     2,0,2,\
     4,2,3))
DefineVectorExpression("eigvals_symm/result", "eigenvalue(<eigvals_symm/tensor>)")
AddPlot("Vector", "eigvals_symm/result")
DrawPlots()
p = PickByNode(0)
eigvals = list(p['eigvals_symm/result'])
TestValueIN("First Eigenvalue of -1", eigvals, -1, Prec)
eigvals.remove(-1)
TestValueIN("Second Eigenvalue of -1", eigvals, -1, Prec)
TestValueIN("Third Eigenvalue of 8", eigvals, 8, Prec)
DeleteAllPlots()
DefineTensorExpression("eigvals_symm/result2", "transpose(eigenvector(<eigvals_symm/tensor>))")
AddPlot("Tensor", "eigvals_symm/result2")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
vecs = ExtractIterablesFromString(s, '()', Prec+1)
TestValueIN("First Eigenvector of (1,-2,0)", vecs, (1,-2,0), Prec, EqualEigVecs)
TestValueIN("Second Eigenvector of (4,2,-5)", vecs, (4,2,-5), Prec, EqualEigVecs)
TestValueIN("Third Eigenvector of (2,1,2)", vecs, (2,1,2), Prec, EqualEigVecs)
DeleteAllPlots()

# Confirm principal_tensor function gives same result as above
TestSection("Cross Principal Stresses and Eigenvalues")
DefineVectorExpression("pcomps_symm/result", "principal_tensor(<eigvals_symm/tensor>)")
AddPlot("Vector", "pcomps_symm/result")
DrawPlots()
p = PickByNode(0)
pcomps = list(p['pcomps_symm/result'])
TestValueIN("First principal component is first eigenvalue", pcomps, -1, Prec)
pcomps.remove(-1) # elim the first of the expected two -1 eigvals
TestValueIN("Second principal component is second eigenvalue", pcomps, -1, Prec)
TestValueIN("Third principal component is third eigenvalue", pcomps, 8, Prec)

TestSection("2D, Symmetric Eigenvalues and Eigenvectors")
CreateConstantTensorExpr("eigvals_symm_2d", meshName, "nodal",\
    (2,  3,  0,\
     3,  2,  0,\
     0,  0,  0))
DefineVectorExpression("eigvals_symm_2d/result", "eigenvalue(<eigvals_symm_2d/tensor>)")
AddPlot("Vector", "eigvals_symm_2d/result")
DrawPlots()
p = PickByNode(0,('eigvals_symm_2d/result',))
eigvals = list(p['eigvals_symm_2d/result'])
TestValueIN("First Eigenvalue of -1", eigvals, -1, Prec)
TestValueIN("Second Eigenvalue of 5", eigvals, 5, Prec)
DeleteAllPlots()
DefineTensorExpression("eigvals_symm_2d/result2", "transpose(eigenvector(<eigvals_symm_2d/tensor>))")
AddPlot("Tensor", "eigvals_symm_2d/result2")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
vecs = ExtractIterablesFromString(s, '()', Prec)
vecs.remove([0,0,1]) # we have to take out the Z eigenvector
TestValueIN("First Eigenvector of (1,-1)", vecs, (1,-1), Prec, EqualEigVecs)
TestValueIN("Second Eigenvector of (1,1)", vecs, (1,1), Prec, EqualEigVecs)

#
# Test a case where eigenvalues are complex (e.g. imaginary)
# The real eigenvalues are 2, (4+3i)/5, (4-3i)/5 but what you
# get from VisIt is 2, 7/5, 1/5 (as though i==1 in the above).
#
TestSection("3D, Complex Eigenvalues and Eigenvectors")
CreateConstantTensorExpr("eigvals_complex", meshName, "nodal",\
    (4.0/5.0, -3.0/5.0,     0,\
     3.0/5.0,  4.0/5.0,     0,\
        1,        2,        2))
DefineVectorExpression("eigvals_complex/result", "eigenvalue(<eigvals_complex/tensor>)")
AddPlot("Vector", "eigvals_complex/result")
DrawPlots()
p = PickByNode(0)
eigvals = list(p['eigvals_complex/result'])
TestValueIN("First Eigenvalue of 2", eigvals, 2, Prec)
TestValueIN("Second Eigenvalue of (4+3i)/5", eigvals, float(4+3)/5.0, Prec)
TestValueIN("Third Eigenvalue of (4-3i)/5", eigvals, float(4-3)/5.0, Prec)
DeleteAllPlots()
DefineTensorExpression("eigvals_complex/result2", "transpose(eigenvector(<eigvals_complex/tensor>))")
AddPlot("Tensor", "eigvals_complex/result2")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
vecs = ExtractIterablesFromString(s, '()', Prec)
TestValueIN("First Eigenvector of (0,0,1)", vecs, (0,0,1), Prec, EqualEigVecs)
DeleteAllPlots()

# Re-use eigvals_symm here
TestSection("Cross Check Deviatoric and Principal Stresses")
DefineVectorExpression("eigvals_symm/dev", "principal_deviatoric_tensor(<eigvals_symm/tensor>)")
AddPlot("Vector", "eigvals_symm/dev")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
dev_vec = ExtractIterablesFromString(s, '()', Prec)
DeleteAllPlots()
DefineScalarExpression("eigvals_symm/tr3", "trace(<eigvals_symm/tensor>)/3.0")
DefineTensorExpression("eigvals_symm/tensor3",\
"""
    {{<eigvals_symm/tensor>[0][0]-<eigvals_symm/tr3>, <eigvals_symm/tensor>[0][1], <eigvals_symm/tensor>[0][2]},
     {<eigvals_symm/tensor>[1][0], <eigvals_symm/tensor>[1][1]-<eigvals_symm/tr3>, <eigvals_symm/tensor>[1][2]},
     {<eigvals_symm/tensor>[2][0], <eigvals_symm/tensor>[2][1], <eigvals_symm/tensor>[2][2]-<eigvals_symm/tr3>}}
"""
)
DefineVectorExpression("eigvals_symm/dev2", "principal_tensor(<eigvals_symm/tensor3>)")
AddPlot("Vector", "eigvals_symm/dev2")
DrawPlots()
PickByNode(0)
s = GetPickOutput()
dev2_vec = ExtractIterablesFromString(s, '()', Prec)
DeleteAllPlots()
TestValueEQ("Principal deviatoric and principal-tr()/3 agree", dev_vec[0], dev2_vec[0])

Exit()