2D finite strain tensor functions

gistfile1.py

import numpy as np

# fst_var -> the swarm variable for the finite strain tensor
# vField  -> the velocity field

gradV      = vField.fn_gradient
srFn       = fn.tensor.symmetric( gradV )
a_velocity = 0.5*(gradV[1] - gradV[2]) # angular velocity

def rotateTensor2D(t, theta):
    # vectorized version of 
    # t' = [Q^T] [t] [Q]
    
    '''
    # In the case tensor is a 2,2 numpy array. eg
    # t  = [[10, 3],[3,5]]
    
    Q = np.array([
            [ np.cos(theta) , -np.sin(theta) ],
            [ np.sin(theta) ,  np.cos(theta) ]])
    
    x = np.matmul(t,Q)
    return np.matmul(Q.T,x)
    '''
    
    # for symmetric tensor t. Using Q as above this simplifies to a 3x3
    t_ = t.copy()
    t_[:,0] = t[:,0]*np.cos(theta)**2 + t[:,1]*np.sin(theta)**2 + t[:,2]*np.sin(2*theta[:])
    t_[:,1] = t[:,0]*np.sin(theta)**2 + t[:,1]*np.cos(theta)**2 - t[:,2]*np.sin(2*theta[:])
    t_[:,2] = ( t[:,1] - t[:,0] ) *np.cos(theta[:])*np.sin(theta[:])  + t[:,2]*np.cos(2*theta[:])
    
    return t_

def update(dt):
    """
    timestepping pseudocode:

    rotate(fs)
    fs += dt * strainRate
    particle advect

    """

    dtheta = dt * a_velocity.evaluate(swarm).reshape(1,-1)
    
    fst_var.data[:,0:3] = rotateTensor2D(fst_var.data[:,0:3], dtheta)
    fst_var.data[:,0:3] += srFn.evaluate(swarm)
   
    advector.integrate(dt)