How to assign resolved shear stress for every slipsystem in Crystal plasticity

1

Hello,

I am new to Fenics. I am trying to implement a simple crystal plasticity model, where I need to calculate the resolved shear stress on each slip system as
tau(i)= 1/ J * m(i) . S * n(i)
where J = det(F_e), m(i) and n(i) are slip direction and slip plane normal of slip system i and S is second Piola Kirchoff stress. Here is a minimal working example.

from __future__ import print_function
from fenics import *
from ufl import *
import numpy as np 
from dolfin import *

# constants
size = 1.0; 
mu = 41.0
lamb = 87.2
nSlip = 12 


# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(size, size, size), 10,10,10) # refinement, refinement, refinement)

# function spaces
V = VectorFunctionSpace(mesh, 'P', 1)
deg_u = 1
deg_stress = 1
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=9, quad_scheme='default')
W = FunctionSpace(mesh, We)
Ze = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=nSlip, quad_scheme='default')
Z = FunctionSpace(mesh, Ze)
Z0e = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
Z0 = FunctionSpace(mesh, Z0e)


# Define variational problem
u = Function(V)

v = TrialFunction(V)

F_p=Function(W) # plastic deformation
tau = Function(Z) # resolved shear stress on the slip systems  
d = u.geometric_dimension()  # space dimension
I=Identity(d)

# slip directions
sno_ = as_tensor(Constant( ((1, 1, 1), \
                 (1, 1, 1), \
                  (1, 1, 1), \
                 (-1, -1, 1), \
                 (-1, -1, 1), \
                 (-1, -1, 1), \
                 (-1, 1, 1), \
                 (-1, 1, 1), \
                 (-1, 1, 1), \
                 (1, -1, 1), \
                 (1,-1, 1), \
                 (1, -1, 1))) )

sno = sno_ * 1.0/np.sqrt(3.0)

# slip plane normals
sdo_ = as_tensor(Constant ( ((0, 1, -1), \
                 (-1, 0, 1), \
                 (1, -1, 0), \
                 (0, -1, -1), \
                 (1, 0, 1), \
                 (-1, 1, 0), \
                 (0, 1, -1), \
                 (1, 0, 1), \
                 (-1, -1, 0), \
                 (0, -1, -1), \
                 (-1, 0, 1), \
                 (1, 1, 0)) ) )

sdo= sdo_ * 1.0/np.sqrt(2.0)

# stiffness tensor 
def C_int(d): # C_intermediate
    tmp = np.zeros((d,d,d,d))
    
    for i in range(0,d):
        for j in range(0,d):
            for k in range(0,d):
                 for l in range(0,d):
                     if ((i==k) and (j==l)): 
                         tmp[i][j][k][l] = tmp[i][j][k][l] + mu 

                     if ((i==l) and (j==k)): 
                         tmp[i][j][k][l] = tmp[i][j][k][l] + mu

                     if ((i==j) and (k==l)): 
                         tmp[i][j][k][l] = tmp[i][j][k][l] + lamb

    return as_tensor(tmp)

# get elastic deformation
def get_Fe(u,F_p):
    F= I + grad(u)
    G_p = inv(as_3D_tensor(F_p))
    F_e = F*G_p
    return F_e

# get Second Piola Kirchoff stress
def get_Sint(F_e):
    E_e = 0.5*(F_e.T*F_e - I)                       # Right Cauchy-Green tensor
    i,j,k,l = indices(4)
    S_int= as_tensor(C_int(d)[i,j,k,l]*E_e[k,l],(i,j))
    return S_int

# get Resolved shear stress on all Slip systems 
def get_tau(J,S):
    for i in range(0,nSlip):
        tau.sub(i).assign( project(1/J * dot(sdo[0,:], S*sno[0,:]),Z0 ) )
    
    return tau 

def as_3D_tensor(X):
    return as_tensor([[X[0], X[1], X[2]],
                      [X[3], X[4], X[5]],
                      [X[6], X[7], X[8]]])

# initialize F_p
F_p.interpolate(Constant((1,0,0,0,1,0,0,0,1)))

F_e=get_Fe(u,F_p)

S= get_Sint(F_e)

tau=get_tau(det(F_e),S)

So I am unable to populate the function tau (Resolved shear stress) which is a QuadratureElement function of dimension 12 (number of slip systems). I want to populate it in a loop for all the slip systems. So I try to assign the i-th component using sub function and try to project the scalar value inside into the scalar QuadratureElement function space Z0. I am receiving this error

Shapes do not match: <Argument id=139866885635192> and <Product id=139866885324232>.
Traceback (most recent call last):
File "crystalPlasticity_min.py", line 132, in <module>
tau=get_tau(det(F_e),S)
File "crystalPlasticity_min.py", line 116, in get_tau
tau.sub(i).assign( project(1/J * dot(sdo[0,:], S*sno[0,:]),Z0 ) )
File "/usr/lib/python3/dist-packages/dolfin/fem/projection.py", line 129, in project
L = ufl.inner(w, v) * dx
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 169, in inner
return Inner(a, b)
File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 161, in _new_
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
File "/usr/lib/python3/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <Argument id=139866885635192> and <Product id=139866885324232>.

Could you please let me know what I am doing wrong here. Thanks.

1 Like
2

There are several issues that I can see by quickly skimming through it.

  • Edit tau looks good to me. You should just change the projection and it should work
  • Quadrature elements are finicky and should be used only when necessary. In any case, you will have a mismatch of quadrature points when trying to do a L2 projection. Your u is in lagrange space, wherereas tau isn’t.
  • Edit: tau is in Z which has nSlip.
    Try starting with a simpler example to get acquainted with FEniCS and then move your way up. This seems fairly involved at beginner level.

Some general code changes to make it more readable :slight_smile:


# Not advisable to set it globally unless you know it's accurate enough for your form
parameters["form_compiler"]["quadrature_degree"] = deg_stress 

def C_int(d):
    I = np.eye(d)
    II = np.einsum("ik,jl->ijkl", I, I) + np.einsum("il,jk->ijkl", I, I)
    JJ = np.einsum("ij,kl->ijkl", I, I)
    return as_tensor(mu*II + lamb*JJ)

# slip plane normals
sno_ = np.array(( 
        ((1, 1, 1),
        (1, 1, 1),
        (1, 1, 1),
        (-1, -1, 1),
        (-1, -1, 1),
        (-1, -1, 1),
        (-1, 1, 1),
        (-1, 1, 1),
        (-1, 1, 1),
        (1, -1, 1),
        (1,-1, 1),
        (1, -1, 1)) 
)) * 1./3**0.5

sdo_ = np.array((
             (0, 1, -1),
             (-1, 0, 1),
             (1, -1, 0),
             (0, -1, -1),
             (1, 0, 1),
             (-1, 1, 0),
             (0, 1, -1),
             (1, 0, 1),
             (-1, -1, 0),
             (0, -1, -1),
             (-1, 0, 1),
             (1, 1, 0)
    )) * 1.0/np.sqrt(2.0)

# get Resolved shear stress on all Slip systems 
def get_tau(J,S):
    for i in range(0,nSlip):
        mi, ni = sdo_[i], sno_[i]
        tau.sub(i).assign( project(1/J * dot(as_vector(mi), dot(S, as_vector(ni))) , Z0) ) # this should change
    return tau
3 Likes
3

Thanks a lot for your help. I followed the above and created these spaces and function for tau:

Mg= VectorFunctionSpace(mesh, 'DG', 1, dim=nSlip)
Zg= FunctionSpace(mesh,'DG',1)
tau = Function(Mg)

and then

tau.sub(i).assign( project(1/J * dot(as_vector(mi), dot(S, as_vector(ni))) , Zg) )

and this seems to work. I avoided quadrature element function altogether.

A compact form like this also seems to works.

tau= project(as_vector( [ 1.0/ J * dot(as_vector(sdo[i,:]), dot(S,as_vector(sno[i,:])) ) for i in range(0,nSlip) ] ), Mg)

Could you please tell me if local_project instead of project will make this more efficient.

4

I think it should be more efficient given that you are doing a L2 projection cell wise when using local_project. See this and this for more details.