### Matrix assembly takes too long...

Hello,

I am solving the 2d Stokes equations with a viscoelastic constitutive equation in a square domain with some Dirichlet conditions (I use P2/P1/P1/P1 elements for the velocities/pressure/velocity gradient projections/stresses). In order to achieve numerical stability I use the log conformation representation of the stress tensor. In this formulation, one needs to solve for the logarithm of the stress tensor instead solving for the stress tensor itself. To do so, we diagonalize the stress tensor, and we also calculate two other tensors based on dot products of the velocity gradient tensor and a tensor that contains the eigenvectors of the stress tensor. Note that all these processes are implicit and all these tensors are used in the nonlinear variational form. In addition, these operations have to be performed at the Gauss points (not at the nodes of the triangles). I have written 3 functions to perform the diagonalization process, but when I use them the matrix assembly process becomes extremely slow... More specifically the matrix assembly process becomes approximately 30 times slower when I use these functions (log representation) compared to the time that the matrix assembly process takes when I solve directly for the stress tensor (and thus, I do not use these functions..).

Any idea how to speed up my code??? (I use optimization flags and mpi parallel)

Here are the diagonalization functions:

from dolfin import *

#---------------------------------------------------------------------------

# Function eigenspectum

#---------------------------------------------------------------------------

def eigenspectum( Sxx, Syx, Syy ):

# Define trace of the tensor

TrS = Sxx + Syy

# Define discriminant of the characteristic equation

Dis = (Sxx-Syy)*(Sxx-Syy) + 4.0*Syx*Syx

# Calculate eigenvalues

l1 = (TrS+sqrt(Dis))/2.0

l2 = l1-sqrt(Dis)

# Calculate eigenvectors

v1x = conditional(gt(abs(Syx), 1.e-9), Sxx-Syy+sqrt(Dis), 1.0)

v1y = conditional(gt(abs(Syx), 1.e-9), 2.0*Syx, 0.0)

# Normalize eigenvectors

V1_m = sqrt(v1x*v1x+v1y*v1y)

v1x = v1x/V1_m

v1y = v1y/V1_m

return (l1, l2, v1x, v1y)

#---------------------------------------------------------------------------

# Function stretch_rotational_tensors

#---------------------------------------------------------------------------

def stretch_rotational_tensors( Sxx, Syx, Syy, Gxx, Gyx, Gxy, Gyy ):

# Calculate eigenvalues & eigenvectors of S tensor

(l1, l2, v1x, v1y) = eigenspectum( Sxx, Syx, Syy )

# Define R & R transpose tensors

R = as_tensor([\

[v1x, -v1y],\

[v1y, v1x]\

])

RT = transpose(R)

# Define velocity gradient transpose tensor

GUT = as_tensor([\

[Gxx, Gyx],\

[Gxy, Gyy]\

])

# Define M tensor and its diagonal

M = dot(dot(RT,GUT),R)

Md = as_tensor([\

[M[0,0], 0.0 ],\

[0.0, M[1,1]]\

])

# Define B tensor

B = dot(dot(R,Md),RT)

# Define B tensor

Bxx = conditional(gt(abs(l1-l2), 1.e-9), B[0,0], Gxx )

Byx = conditional(gt(abs(l1-l2), 1.e-9), B[0,1], 0.5*(Gyx+Gxy))

Byy = conditional(gt(abs(l1-l2), 1.e-9), B[1,1], Gyy )

# Define W tensor

wmega = (exp(l2)*M[0,1]+exp(l1)*M[1,0])/(exp(l2)-exp(l1))

Wa = as_tensor([\

[0.0, wmega],\

[-wmega, 0.0 ]\

])

W = dot(dot(R,Wa),RT)

Wxy = conditional(gt(abs(l1-l2), 1.e-9), W[0,1], 0.0)

Wyx = conditional(gt(abs(l1-l2), 1.e-9), W[1,0], 0.0)

return (Bxx, Byx, Byy, Wxy, Wyx)

#---------------------------------------------------------------------------

# Function log_conformation

#---------------------------------------------------------------------------

def log_conformation( Sxx, Syx, Syy ):

# Calculate eigenvalues & eigenvectors of S tensor

(l1, l2, v1x, v1y) = eigenspectum( Sxx, Syx, Syy )

# Define R & R transpose tensors

R = as_tensor([\

[v1x, -v1y],\

[v1y, v1x]\

])

RT = transpose(R)

# Define conformation tensor

D = as_tensor([\

[exp(l1), 0.0 ],\

[0.0, exp(l2)]\

])

C = dot(dot(R,D),RT)

Cxx = C[0,0]

Cyx = C[0,1]

Cyy = C[1,1]

# Define e^(-S) tensor

Di = as_tensor([\

[exp(-l1), 0.0 ],\

[0.0, exp(-l2)]\

])

exp_mS = dot(dot(R,Di),RT)

exp_mSxx = exp_mS[0,0]

exp_mSyx = exp_mS[0,1]

exp_mSyy = exp_mS[1,1]

return (Cxx, Cyx, Cyy, exp_mSxx, exp_mSyx, exp_mSyy)