File: GreenLagrangeStrain.py

package info (click to toggle)
vtk7 7.1.1%2Bdfsg2-8
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 127,396 kB
  • sloc: cpp: 1,539,584; ansic: 124,382; python: 78,038; tcl: 47,013; xml: 8,142; yacc: 5,040; java: 4,439; perl: 3,132; lex: 1,926; sh: 1,500; makefile: 126; objc: 83
file content (134 lines) | stat: -rwxr-xr-x 4,651 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#!/usr/bin/env python
#coding=utf8

from __future__ import print_function

import sys
import math
try:
    import numpy
except ImportError:
    print("WARNING: This test requires NumPy: http://http://www.numpy.org/")
    sys.exit(0)
import vtk


# Create a simple hexahedron
points = vtk.vtkPoints()
points.SetNumberOfPoints(8)
points.InsertPoint(0, 0, 0, 0)
points.InsertPoint(1, 1, 0, 0)
points.InsertPoint(2, 1, 1, 0)
points.InsertPoint(3, 0, 1, 0)
points.InsertPoint(4, 0, 0, 1)
points.InsertPoint(5, 1, 0, 1)
points.InsertPoint(6, 1, 1, 1)
points.InsertPoint(7, 0, 1, 1)
#for k_point in xrange(8): print points.GetPoint(k_point)

hexahedron = vtk.vtkHexahedron()
for k_point in xrange(8):
    hexahedron.GetPointIds().SetId(k_point, k_point)

cell_array = vtk.vtkCellArray()
cell_array.InsertNextCell(hexahedron)

farray_disp = vtk.vtkFloatArray()
farray_disp.SetName("displacement")
farray_disp.SetNumberOfComponents(3)
farray_disp.SetNumberOfTuples(8)

ugrid_hex = vtk.vtkUnstructuredGrid()
ugrid_hex.SetPoints(points)
ugrid_hex.SetCells(hexahedron.GetCellType(), cell_array)
ugrid_hex.GetPointData().AddArray(farray_disp)
ugrid_hex.GetPointData().SetActiveVectors("displacement")

# Test a few deformation gradients
F_lst  = []

# Simple extentions
F_lst += [numpy.array([[2,0,0],
                       [0,1,0],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [0,2,0],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [0,1,0],
                       [0,0,2]])]

# Simple shears
F_lst += [numpy.array([[1,1,0],
                       [0,1,0],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,1],
                       [0,1,0],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [0,1,1],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [1,1,0],
                       [0,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [0,1,0],
                       [1,0,1]])]
F_lst += [numpy.array([[1,0,0],
                       [0,1,0],
                       [0,1,1]])]

# Rotation
alpha = math.pi/2
F_lst += [numpy.array([[+math.cos(alpha),-math.sin(alpha),0],
                       [+math.sin(alpha),+math.cos(alpha),0],
                       [               0,               0,1]])]

for F in F_lst:
    print("F = " + str(F))

    # Assuming an homogeneous deformation gradient, compute the displacement and displacement gradient: u = x(X)-X = F.X-X = (F-1).X = GU.X with GU = F-1
    GU = F - numpy.eye(3)
    print("GU = " + str(GU))

    for k_point in xrange(8):
        farray_disp.SetTuple(k_point, numpy.dot(GU, points.GetPoint(k_point)))
        #print farray_disp.GetTuple(k_point)

    # Compute the small (linearized) strain tensor: e = (GU)_sym = (GU+GU^T)/2
    e = (GU + numpy.transpose(GU))/2
    print("e = " + str(e))

    # Compute the Green-Lagrange strain tensor: E = (F.F^T-1)/2 = (GU+GU^T+GU^T.GU)/2 = e + GU^T.GU/2
    E = e + numpy.dot(numpy.transpose(GU), GU)/2
    print("E = " + str(E))

    # Compute displacement gradient using vtkCellDerivatives
    cell_derivatives = vtk.vtkCellDerivatives()
    cell_derivatives.SetVectorModeToPassVectors()
    cell_derivatives.SetTensorModeToComputeGradient()
    cell_derivatives.SetInputData(ugrid_hex)
    cell_derivatives.Update()
    vector_gradient = numpy.reshape(cell_derivatives.GetOutput().GetCellData().GetArray("VectorGradient").GetTuple(0), (3,3))
    print("VectorGradient = " + str(vector_gradient))
    assert numpy.allclose(vector_gradient, GU)

    # Compute small strain tensor using vtkCellDerivatives
    cell_derivatives = vtk.vtkCellDerivatives()
    cell_derivatives.SetVectorModeToPassVectors()
    cell_derivatives.SetTensorModeToComputeStrain()
    cell_derivatives.SetInputData(ugrid_hex)
    cell_derivatives.Update()
    strain = numpy.reshape(cell_derivatives.GetOutput().GetCellData().GetArray("Strain").GetTuple(0), (3,3))
    print("Strain = " + str(strain))
    assert numpy.allclose(strain, e)

    # Compute Green-Lagrange strain tensor using vtkCellDerivatives
    cell_derivatives = vtk.vtkCellDerivatives()
    cell_derivatives.SetVectorModeToPassVectors()
    cell_derivatives.SetTensorModeToComputeGreenLagrangeStrain()
    cell_derivatives.SetInputData(ugrid_hex)
    cell_derivatives.Update()
    green_lagrange_strain = numpy.reshape(cell_derivatives.GetOutput().GetCellData().GetArray("GreenLagrangeStrain").GetTuple(0), (3,3))
    print("GreenLagrangeStrain = " + str(green_lagrange_strain))
    assert numpy.allclose(green_lagrange_strain, E)