Possible error in computing cell derivatives/Green Lagrange Strain

Development
1

I believe there is an error in computing Green Lagrange strain. In the attached code, I am deforming a triangulated sphere with a uniform dilation. The Green Lagrange strain should be diagonal, but it is not. I have checked the expression for the GL strain in vtkCellDerivatives.cxx and it seems right. I think the issue may be in how derivatives are computed on triangles and rotated back to the global modeling frame.

The code below uses pyvista (a Python wrapper on vtk), which can be installed with pip install pyvista.

import pyvista as pv
import numpy as np
import vtk
from pyvista import examples

# Define a sphere
#surf = pv.Pyramid()
surf = pv.Sphere()

# Define deformed points and displacements
points = surf.points
X1 = points[:,0]
X2 = points[:,1]
X3 = points[:,2]

# Dilate sphere uniformly
x1 = 1.5 * X1
x2 = 1.5 * X2
x3 = 1.5 * X3

# Compute displacements and add to mesh
d1 = x1 - X1
d2 = x2 - X2
d3 = x3 - X3

displacements = np.zeros_like(points)
displacements[:,0] = d1
displacements[:,1] = d2
displacements[:,2] = d3
surf.point_data['Displacement'] = displacements


# Compute Green-Lagrange strain tensor using vtkCellDerivatives
cell_derivatives = vtk.vtkCellDerivatives()
cell_derivatives.SetVectorModeToPassVectors()
cell_derivatives.SetTensorModeToComputeGreenLagrangeStrain()
cell_derivatives.SetInputData(surf)
surf.GetPointData().SetActiveVectors('Displacement')
cell_derivatives.Update()


# Wrap output in PyVista mesh and make it polydata by extract_surface()
surf = pv.wrap(cell_derivatives.GetOutput())
surf = surf.extract_surface()


# Explicitly create E matrix
Evec = surf['GreenLagrangeStrain']
E = np.zeros((np.size(Evec, 0), 3, 3))

E[:,0,0] = Evec[:,0]
E[:,0,1] = Evec[:,1]
E[:,0,2] = Evec[:,2]

E[:,1,0] = Evec[:,3]
E[:,1,1] = Evec[:,4]
E[:,1,2] = Evec[:,5]

E[:,2,0] = Evec[:,6]
E[:,2,1] = Evec[:,7]
E[:,2,2] = Evec[:,8]

print('E: \n', E[0,:,:])

#surf.save('./sphere_test_vtk_issue.vtk')
2

Hi, I am attaching a function that computes surface strain correctly for triangulated surfaces. Iā€™m not sure whether the VTK version is wrong or just makes different assumptions, but my version, taken from this paper, considers each triangle to be a shell element, and determines the Green-Lagrange strain in the global coordinate system assuming 1) no shear strain of the shell element in the normal plane, and 2) no normal strain in the normal plane.

def compute_GL_strain(surf):
    '''
    Calculates the Green-Lagrange strain of a surface mesh, assuming it has an array 
    called "Displacements" at each node. The surface mesh should be in the reference
    configuration.

    Args:
        surf: A PyVista polydata surface in the reference configuration

    Returns: 
        surf: A PyVista polydata surface with an added cell vector array containing the 
        components of the Green-Lagrange strain tensor.

    These calculations are identical to Razeghi 2021 (https://www-nature-com.stanford.idm.oclc.org/articles/s41598-021-84935-x)
    Supplementary. The equations are a specific case on shell kinematic theory, 
    assuming 1) no out-of-plane shear strain, and 2) no normal strain in normal direction.
    See surf_def_grad_derivation.pdf for details.
    
    The components of the resulting VTK array correspond to
    GreenLagrangeStrain[0] = E_xx
    GreenLagrangeStrain[1] = E_xy
    GreenLagrangeStrain[2] = E_xz
    GreenLagrangeStrain[3] = E_yx
    GreenLagrangeStrain[4] = E_yy
    GreenLagrangeStrain[5] = E_yz
    GreenLagrangeStrain[6] = E_zx
    GreenLagrangeStrain[7] = E_zy
    GreenLagrangeStrain[8] = E_zz

    For interpretation of these components: 
    E_xx positive, for example, indicates stretching in the x direction
               
               ----------                          ^ y
               |        |                          |
          <--- |        |  --->   E_xx > 0         |
               |        |                          ------>  x
               ----------
    E_xy positive, for example, indicates shearing in the xy plane
           
                 ------>
               ----------                          ^ y
             | |        | ^                        |
             | |        | |   E_xy > 0             |
             v |        | |                        ------>  x
               ----------
                <-------
    '''

    ## --- Get cell points array in REFERENCE configuration --- ##
    # Get list of point indices for each cell
    cell_indices = surf.faces
    # Make sure all triangles
    assert not cell_indices.size % 4 and np.all(cell_indices.reshape(-1, 4)[:,0] == 3)
    # Reshape to get ncells x 3 array of nodal indices
    cell_node_ids = cell_indices.reshape(-1, 4)[:,1:4].ravel()
    cell_node_ids = cell_node_ids.reshape(-1, 3) # ncells x 3
    #print(cell_node_ids)

    # Get an array of points for each cell
    # cell_nodes[cell_id, node_id, dim]. cell_id = [0, ncells-1], node_id = [0,2], dim = [0,2]
    cell_nodes = surf.points[cell_node_ids] # ncells x 3 x 3

    # Reference configuration nodal positions for each cell
    X = cell_nodes

    ## --- Get cell points array in CURRENT configuration --- ##
    # Warp geometry by displacement
    surf_cur = surf.warp_by_vector('Displacement', factor = 1)
    # Get list of point indices for each cell
    cell_indices = surf_cur.faces
    # Make sure all triangles
    assert not cell_indices.size % 4 and np.all(cell_indices.reshape(-1, 4)[:,0] == 3)
    # Reshape to get ncells x 3 array of nodal indices
    cell_node_ids = cell_indices.reshape(-1, 4)[:,1:4].ravel()
    cell_node_ids = cell_node_ids.reshape(-1, 3) # ncells x 3
    #print(cell_node_ids)

    # Get an array of points for each cell
    # cell_nodes[cell_id, node_id, dim]. cell_id = [0, ncells-1], node_id = [0,2], dim = [0,2]
    cell_nodes = surf_cur.points[cell_node_ids] # ncells x 3 x 3

    # Reference configuration nodal positions for each cell
    x = cell_nodes

    ncells = surf.n_faces
    
    ## --- Compute deformation gradient --- ##

    # Construct A1, A2, N, and a1, a2, N
    A1 = X[:, 2-1] - X[:, 1-1] # -1 to switch to 0-indexing
    A2 = X[:, 3-1] - X[:, 1-1]
    N = np.cross(A1, A2)
    N = N / np.linalg.norm(N, axis=1)[:, np.newaxis]
    # Check
    #surf.compute_normals(point_normals=False, inplace=True)
    #normals = surf['Normals']
    #print(np.linalg.norm(N - normals))

    a1 = x[:, 2-1] - x[:, 1-1]
    a2 = x[:, 3-1] - x[:, 1-1]
    n = np.cross(a1, a2)
    n = n / np.linalg.norm(n, axis=1)[:, np.newaxis]

    # Form AAN and aan matrices
    AAN = np.empty((ncells, 3, 3))
    AAN[:, :, 0] = A1
    AAN[:, :, 1] = A2
    AAN[:, :, 2] = N

    aan = np.empty((ncells, 3, 3))
    aan[:, :, 0] = a1
    aan[:, :, 1] = a2
    aan[:, :, 2] = n

    # Compute deformation gradient from aan = F * AAN -> F = aan * inv(AAN)
    F = aan @ np.linalg.inv(AAN)
    #print(F.shape)

    # Compute Green-Lagrange strain tensor
    FT = np.transpose(F, (0, 2, 1))
    E = 1/2 * (FT @ F - np.identity(3))
    #print(E.shape)

    # Populate surf with GL strain components
    E_vec = np.empty((ncells, 9))
    E_vec[:,0] = E[:, 0, 0]
    E_vec[:,1] = E[:, 0, 1]
    E_vec[:,2] = E[:, 0, 2]
    E_vec[:,3] = E[:, 1, 0]
    E_vec[:,4] = E[:, 1, 1]
    E_vec[:,5] = E[:, 1, 2]
    E_vec[:,6] = E[:, 2, 0]
    E_vec[:,7] = E[:, 2, 1]
    E_vec[:,8] = E[:, 2, 2]
    surf.cell_data["GreenLagrangeStrain"] = E_vec
    
    return surf
1 Like