Incremental implementation of hyperelasticity

variational formulation
1

Greetings,

A fenics newbie presenting himself. I would like to implement a hyperelastic problem using the linearized form I deduced myself using directional derivative, i.e., I do not want to use the fancy nonlinear problem function or other wrapper that does everything under the hood.

To this end, I implemented the fourth order elastic tensor and all its dependencies, but I’m really struggling with the solution of the variational problem, for I have the solution of the previous iteration, u, and the current update direction, Delta_u, but I do not know exactly how to update u with Delta_u inside the incremental loop. Do you have any ideas?

Besides, as I derived the fourth order elastic tensor, there is no need to compute the tangent operator using automatic differentiation, is there a way to guarantee that fenics does not use it?

2

Not sure to understand everything. What is Delta_u ? A correction in a home made Newton solver ? Please post a MWE so that we can help.

3

u_{k+1} = u_{k}+Delta_u,

where u is the displacement field and Delta_u is the update direction in the Newton-Raphson procedure. I send below a snippet of my code, where the home made Newton solver is defined. I can’t do the update rule expressed above (at least the residual norm does not shrink at all, which means either Delta_u is null or the update of u is not working).

’ ’ ’

Forma variacional

a = ((inner((F.T)*grad(delta_u),tensores_elasticos.dupla_contracao(
C_elastico,(F.T)*grad(Delta_u)))*dx))+(inner(grad(delta_u)*S,grad(
Delta_u))*dx)

L = (-(inner(grad(delta_u),F*S)*dx)+(dot(B,delta_u)dx)+(dot(T,delta_u)
ds))

Delta_u = Function(V)

Itera ao longo das iterações do método de Newton-Raphson

for i in range(n_maxIter):

print("\nInicia a iteração", i+1, "do método de Newton-Raphson\n")

# Soluciona o problema linearizado

solve(a==L, Delta_u, CCs_dirichlet)

# Atualiza o campo de deslocamento

u = project(u+Delta_u, V)
#u.vector[:] = u.vector[:]+Delta_u.vector[:]

# Calcula o resíduo

R = assemble(L)

for CC in CCs_dirichlet:

    CC.apply(R)

print("Norma do resíduo é ", norm(R))

’ ’ ’

4

Please format the code correctly and provide a full MWE with a simple mesh and loading so that people can help. Translating to english would also be easier

5

I will work on it. Can’t I simply send the python file? Which is the correct way for formatting code in this forum?

6

Correct formatting is using 3x`,
Ie

```python

#add code here

```

See for instance Extended Syntax | Markdown Guide

The reason for not having file uploads (except pictures) is to avoid dumping of large files that are far from minimal reproducible examples.

1 Like
7

I will send this reply in three parts, with the code split in two messages, following the code structure. First, I will send the routine for evaluation of the fourth order elastic tensors (following Holzapfel 2000 notation) and, then, I will send the main code.

I’ve translated everything to English and everything is (methodically) commented. I would really like to implement this as a home made Newton.

8

file: elastic_tensors.py


from fenics import *

import ufl_legacy as ufl

i,j,k,l = ufl.indices(4)

# Defines a function to perform the operation of double contraction be-
# tween a fourth order tensor and a second order tensor

def double_contraction(C,A):

    return as_tensor(C[i,j,k,l]*A[k,l], (i,j))

# Defines a function to evaluate the fourth order elastic tensor using
# hyperelastic models that do not use isochoric invariants of the right
# Cauchy-Green tensor

def C_hyperNonIsochoric(C, I1, kvol, J, dkvol_dJ, dPsi_dI2,
ddPsi_ddI1, ddPsi_ddI2, ddPsi_dI1dI2, dimension):
    
    # Initiates the Kronecker's delta
    
    delta = Identity(dimension)
    
    # Calculates the inverse of the right Cauchy-Green tensor, C
    
    C_inv = inv(C)

    # Calculates the fourth order identity tensor
    
    II = as_tensor(delta[i,k]*delta[j,l], (i,j,k,l))

    # Calculates the ball product of the inverse of C (see Holzapfel(2000)
    # page 43 for further information)
    
    C1 = as_tensor((-0.5*((C_inv[i,k]*C_inv[l,j])+(C_inv[i,l]*C_inv[k,j]
    ))),(i,j,k,l))

    # Calculates the tensorial product of the inverse of C by itself
    
    C2 = as_tensor(C_inv[i,j]*C_inv[k,l], (i,j,k,l))

    # Calculates the tensorial product of C by itself
    
    C3 = as_tensor(C[i,j]*C[k,l], (i,j,k,l))

    # Calculates the symmetric part of the tensorial product of the se-
    # cond order identity tensor by C
    
    C4 = as_tensor(0.5*((delta[i,j]*C[k,l])+(delta[k,l]*C[i,j])), 
    (i,j,k,l))

    # Calculates the tensorial product of the second order identity ten-
    # sor by itself
    
    IxI = as_tensor(delta[i,j]*delta[k,l], (i,j,k,l))

    ### Calculas the coefficients that multiply the tensor evaluated a-
    ### bove

    # Calculates the coefficient of the tensorial product of the second
    # order identity tensor by itself
    
    coeff_IxI = (dPsi_dI2+ddPsi_ddI1+(2*I1*ddPsi_dI1dI2)+(I1*I1*
    ddPsi_ddI2))

    # Calculates the coefficient of the symmetric part of the tensorial
    # product of the second order identity tensor by C

    coeff_IxC = -2*(ddPsi_dI1dI2+(I1*ddPsi_ddI2))

    # Sums everything up and returns it

    return ((4*coeff_IxI*IxI)-(4*coeff_IxC*C4)+(4*ddPsi_ddI2*C3)-(4*
    dPsi_dI2*II)+(2*J*dkvol_dJ*C2)-(4*kvol*C1))
9

file: hyperelasticity.py

from dolfin import *

import ufl_legacy as ufl

i,j,k,l = ufl.indices(4)

import elastic_tensors

########################################################################
#                                 Mesh                                 #
########################################################################

# Defines the degree of the interpolation polynomial 

polynomial_degree = 1

# Defines the mesh

mesh = UnitCubeMesh(16, 16, 24)#BoxMesh(Point(0,0,0), Point(1,1,1), 1,1,1)

# Defines the solution space using polynomial interpolation functions

V = VectorFunctionSpace(mesh, "P", polynomial_degree)

########################################################################
#                        Newton-Raphson settings                       #
########################################################################

# Defines the maximum number of iterations

n_maxIter = 5

# Defines the norm of the residual to be used as stopping criterion to 
# the Newton-Raphson method

residual_tol = 1E-5

########################################################################
#                    Dirichlet boundary conditions                     #
########################################################################

# Defines domain subregions

below =  CompiledSubDomain("near(x[2], side) && on_boundary", side=0.0)

above = CompiledSubDomain("near(x[2], side) && on_boundary", side=1.0)

# Defines the expressions for the Dirichlet boundary conditions, i.e.,
# their values

expression_below = Expression(("0.0", "0.0", "0.0"), degree=polynomial_degree)

expression_above = Expression(("0.0", "0.0", "1.0"), degree=polynomial_degree)

# Defines two objects of Dirichlet boundary conditions

bc_DirichletBelow = DirichletBC(V, expression_below, below)

bc_DirichletAbove = DirichletBC(V, expression_above, above)

BCs_dirichlet = [bc_DirichletBelow, bc_DirichletAbove]

########################################################################
#                 Definition of solutions and variation                #
########################################################################

# Defines the function for the update direction in the Newton-Raphson 
# scheme

Delta_u = TrialFunction(V)

# Defines the variation

delta_u = TestFunction(V)

# Defines the displacement field

u = Function(V)

# Defines the body forces vector in the reference configuration

B = Constant((0.0, 0.0, 1.0))

# Defines the boundary traction vector in the reference configuration

T = Constant((0.0, 0.0, 0.0))

########################################################################
#                          Tensors definition                          #
########################################################################

# Recovers the problem's dimension

dimension = Delta_u.geometric_dimension()

# Initializes the second order identity tensor

I = Identity(dimension)

# Calculates the deformation gradient tensor in terms of the solution u

F = grad(u)+I

# Calculates the right Cauchy-Green tensor

C = (F.T)*F

# Calculates the first two invariants of the right Cauchy-Green tensor

I1 = tr(C)

I2 = 0.5*((tr(C)**2)-tr(C*C))

# Calculates the jacobian of the kinematic mapping from the reference to
# the deformed configurations

J = det(F)

########################################################################
#                        Constitutive modelling                        #
########################################################################
 
# Saint Venant-Kirchhoff model

E = Constant(200E6)

nu = Constant(0.3) 

mu = E/(2*(1+nu))

lmbda = (E*nu)/((1+nu)*(1-2*nu))

# Defines the Helmholtz energy function

Psi = ((0.5*mu*(1-I2))+(0.25*mu*((I1-1)**2))+(0.125*lmbda*((I1-3)**2)))

# Defines the derivatives of the energy function w.r.t. the invariants 
# of the right Cauchy-Green tensor up to the second order

dPsi_dI1 = (0.5*mu*(I1-1))+(0.25*lmbda*(I1-3))

dPsi_dI2 = -0.5*mu

ddPsi_ddI1 = (0.5*mu)+(0.25*lmbda)

ddPsi_ddI2 = 0.0

ddPsi_dI1dI2 = 0.0

# Defines the k factor for the volumetric part of the second Piola-
# Kirchhoff stress tensor and the derivative of the former w.r.t. the
# jacobian

k_vol = 0.0

dkvol_dJ = 0.0

# Calculates the second Piola-Kirchhoff stress tensor

S = (2*(((dPsi_dI1+(I1*dPsi_dI2))*I)-(dPsi_dI2*C)+(k_vol*inv(C))))

# Defines the fourth order elastic tensor

C_elastic = elastic_tensors.C_hyperNonIsochoric(C, I1, k_vol, J,
dkvol_dJ, dPsi_dI2, ddPsi_ddI1, ddPsi_ddI2, ddPsi_dI1dI2, dimension)

########################################################################
#               Variational form and solution procedure                #
########################################################################

# Defines the linearized variational form in terms of the update direc-
# tion, a==L. Note that a is constituted by the directional derivative 
# of the variation of the power/work; while L is formed by the residual
# itself. This scheme states the Newton-Raphson method into a variational
# form

a = ((inner((F.T)*grad(delta_u),elastic_tensors.double_contraction(
C_elastic,(F.T)*grad(Delta_u)))*dx))+(inner(grad(delta_u)*S,grad(
Delta_u))*dx)

L = (-(inner(grad(delta_u),F*S)*dx)+(dot(B,delta_u)*dx)+(dot(T,delta_u)*
ds))

Delta_u = Function(V)

# Iterates through the loop of the Newton-Raphson method

for i in range(n_maxIter):

    print("\nInitiates iteration", i+1, "of the Newton-Raphson method\n")

    # Solves the linearized variational form
    
    solve(a==L, Delta_u, BCs_dirichlet)

    # Updates the displacement field

    u = project(u+Delta_u, V)
    #u.vector[:] = u.vector[:]+Delta_u.vector[:]

    # Evaluates the residual
    
    R = assemble(L)

    for CC in BCs_dirichlet:

        CC.apply(R)

    print("Residual norm is", norm(R))

# Saves the solution into a VTK file

file = File("displacement.pvd");

file << u;
10

With the line u = project(u + Delta_u, V) you are creating a new Function each time which is no longer associated with the oeiginal form. Also when implementing a custom newton method with non zero Dirichlet BCs you should pay attention to these as you will apply u=1 for each iteration correction.

1 Like
11

Thank you! Three questions:

  1. Which is the work around u = project(…)?
  2. Totally, the boundary conditions must be applied to u; Delta_u must be null to those DOFs. I will come up with an idea to this problem.
  3. Is there a way to guarantee that fenics does not compute the tangent operator using automatic differentiation? As I’ve already evaluated the fourth order elastic tensor.
12
  1. Update the array of values directly since all live in the same function space.
  2. Custom Newton Solver problem with Dirichlet conditions should help
  3. It does not if you don’t ask for it with ufl.derivative. Here you are defining the tangent form yourself so it’s OK
1 Like
13

Thank you very much for your response. I’ve made a series of alterations following your invaluable tips. There is the code:


from dolfin import *

import ufl_legacy as ufl

from petsc4py import PETSc

i,j,k,l = ufl.indices(4)

import elastic_tensors

########################################################################
#                                 Mesh                                 #
########################################################################

# Defines the degree of the interpolation polynomial 

polynomial_degree = 1

element_type = 'P'

# Defines the mesh

mesh = UnitCubeMesh(16, 16, 24)#BoxMesh(Point(0,0,0), Point(1,1,1), 1,1,1)

# Defines the solution space using polynomial interpolation functions

V = VectorFunctionSpace(mesh, element_type, polynomial_degree)

########################################################################
#                        Newton-Raphson settings                       #
########################################################################

# Sets the final time

final_time = 1.0

# Sets the number of pseudotime points for loading

n_pseudoTimePoints = 1

# Sets the step in pseudo time

delta_t = final_time/n_pseudoTimePoints

# Defines the maximum number of iterations

n_maxIter = 40

# Defines a pseudotime variable to control the incremental process

t = 0.0

# Defines a function of the pseudo time that gives the displacement a-
# long the incremental process

def u_x(t):
    
    return t*0.0

def u_y(t):
    
    return t*0.0

def u_z(t):
    
    return t*0.0

# Defines the norm of the residual to be used as stopping criterion to 
# the Newton-Raphson method

residual_tol = 1E-4

########################################################################
#                    Dirichlet boundary conditions                     #
########################################################################

# Defines domain subregions

below =  CompiledSubDomain("near(x[2], side) && on_boundary", side=0.0)

above = CompiledSubDomain("near(x[2], side) && on_boundary", side=1.0)

# Defines the expressions for the Dirichlet boundary conditions, i.e.,
# their values for the displacement field

expression_below = Expression(("0.0", "0.0", "0.0"), degree=
polynomial_degree)

def expression_above(t):
    
    return Expression((str(u_x(t)), str(u_y(t)), str(u_z(t))), degree=
    polynomial_degree)

# Defines two objects of Dirichlet boundary conditions for the displace-
# ment field

bc_DirichletBelowDisplacement = DirichletBC(V, expression_below, below)

def bc_DirichletAboveDisplacement(t):
    
    return DirichletBC(V, expression_above(t), above)

def BCs_dirichletDisplacement(t):
    
    return [bc_DirichletBelowDisplacement, 
    bc_DirichletAboveDisplacement(t)]

# The update direction must be null in the regions where Dirichlet boun-
# dary conditions are applied to the displacement field

bc_DirichletAboveDirection = DirichletBC(V, Constant(("0.0", "0.0",
"0.0")), above)

bc_DirichletBelowDirection = DirichletBC(V, Constant(("0.0", "0.0",
"0.0")), below)

BCs_dirichletDirection = [bc_DirichletBelowDirection,
bc_DirichletAboveDirection]

########################################################################
#                 Definition of solutions and variation                #
########################################################################

# Defines the function for the update direction in the Newton-Raphson 
# scheme

Delta_u = TrialFunction(V)

# Defines the variation

delta_u = TestFunction(V)

# Defines the displacement field

u = Function(V)

# Defines the body forces vector in the reference configuration

B = Constant((0.0, 0.0, 0.0))

# Defines the boundary traction vector in the reference configuration

T = Constant((0.0, 0.0, 1.0))

########################################################################
#                          Tensors definition                          #
########################################################################

# Recovers the problem's dimension

dimension = Delta_u.geometric_dimension()

# Initializes the second order identity tensor

I = Identity(dimension)

# Calculates the deformation gradient tensor in terms of the solution u

F = grad(u)+I

# Calculates the right Cauchy-Green tensor

C = (F.T)*F

# Calculates the first two invariants of the right Cauchy-Green tensor

I1 = tr(C)

I2 = 0.5*((tr(C)**2)-tr(C*C))

# Calculates the jacobian of the kinematic mapping from the reference to
# the deformed configurations

J = det(F)

########################################################################
#                        Constitutive modelling                        #
########################################################################
 
# Saint Venant-Kirchhoff model

E = Constant(200E6)

nu = Constant(0.3) 

mu = E/(2*(1+nu))

lmbda = (E*nu)/((1+nu)*(1-2*nu))

# Defines the Helmholtz energy function

Psi = ((0.5*mu*(1-I2))+(0.25*mu*((I1-1)**2))+(0.125*lmbda*((I1-3)**2)))

# Defines the derivatives of the energy function w.r.t. the invariants 
# of the right Cauchy-Green tensor up to the second order

dPsi_dI1 = (0.5*mu*(I1-1))+(0.25*lmbda*(I1-3))

dPsi_dI2 = -0.5*mu

ddPsi_ddI1 = (0.5*mu)+(0.25*lmbda)

ddPsi_ddI2 = 0.0

ddPsi_dI1dI2 = 0.0

# Defines the k factor for the volumetric part of the second Piola-
# Kirchhoff stress tensor and the derivative of the former w.r.t. the
# jacobian

k_vol = 0.0

dkvol_dJ = 0.0

# Calculates the second Piola-Kirchhoff stress tensor

S = (2*(((dPsi_dI1+(I1*dPsi_dI2))*I)-(dPsi_dI2*C)+(k_vol*inv(C))))

C_variable = variable(C)

def S_variableFunction(C_variable):

    I1 = tr(C_variable)

    I2 = 0.5*((I1*I1)-tr(C_variable*C_variable))
    
    dPsi_dI1 = (0.5*mu*(I1-1))+(0.25*lmbda*(I1-3))
    
    dPsi_dI2 = -0.5*mu

    k_vol = 0.0

    return (2*(((dPsi_dI1+(I1*dPsi_dI2))*Identity(3))-(dPsi_dI2*C_variable
    )+(k_vol*inv(C_variable))))

S_variable = S_variableFunction(C_variable)

# Defines the fourth order elastic tensor

C_elastic = elastic_tensors.C_hyperNonIsochoric(C, I1, k_vol, J,
dkvol_dJ, dPsi_dI2, ddPsi_ddI1, ddPsi_ddI2, ddPsi_dI1dI2, dimension)

########################################################################
#               Variational form and solution procedure                #
########################################################################

# Defines the linearized variational form in terms of the update direc-
# tion, a==L. Note that a is constituted by the directional derivative 
# of the variation of the power/work; while L is formed by the residual
# itself. This scheme states the Newton-Raphson method into a variational
# form

L = (-(inner(grad(delta_u),F*S)*dx)+(dot(B,delta_u)*dx)+(dot(T,delta_u)*
ds))

# Analytical derivative of the fourth order elastic tensor

#a = ((inner((F.T)*grad(delta_u),elastic_tensors.double_contraction(
#C_elastic,(F.T)*grad(Delta_u)))*dx))+(inner(grad(delta_u)*S,grad(
#Delta_u))*dx)

# Numerical derivative of the fourth order elastic tensor

#a = ((inner((F.T)*grad(delta_u),elastic_tensors.double_contraction(
#2*diff(S_variable, C_variable),(F.T)*grad(Delta_u)))*dx))+(inner(grad(delta_u)*S,grad(
#Delta_u))*dx)

# Numerical derivative of the residue entire

a = -derivative(L, u, Delta_u)

Delta_u = Function(V)

# Iterates through the pseudotime points

for j in range(n_pseudoTimePoints):

    # Updates time

    t += delta_t

    print("\n#########################################################"+
    "\n#           Initiates the", j,"pseudotime point            #\n"+
    "#########################################################\n")

    print("t =", t)

    # Applies boundary conditions to the displacement field

    BCs_dirichletU = BCs_dirichletDisplacement(t)

    for CC in BCs_dirichletU:

        CC.apply(u.vector())

    # Iterates through the loop of the Newton-Raphson method

    for i in range(n_maxIter):

        print("\nInitiates iteration", i+1, "of the Newton-Raphson met"+
        "hod\n")

        # Solves the linearized variational form
        
        #solve(a==L, Delta_u, BCs_dirichletDirection)
        Kt, R = assemble_system(a, L, BCs_dirichletDirection)

        ### PETSc solver

        # Gets the tangent matrix

        Kt_matrix = as_backend_type(Kt).mat()

        vectorDelta_u, R_vector = Kt_matrix.getVecs()

        # Gets the residual vector

        R_vector.setValues(range(R.size()), R)

        # Sets the Krylov linear solver. Uses conjugate gradient ('cg')

        krylov_linearSolver = PETSc.KSP().create()

        krylov_linearSolver.setType('cg')

        krylov_linearSolver.setOperators(Kt_matrix)

        krylov_linearSolver.solve(R_vector, vectorDelta_u)

        # Evaluates the norm of the residual

        norm_residual = norm(R, "l2")

        print("Residual norm is", norm_residual)

        # Verifies stop criterion

        if norm_residual<residual_tol:

            print("Residual criterion met at iteration", i)

            break

        # Updates the displacement field
        
        u.vector()[:] = u.vector()[:]+vectorDelta_u

    print("\n\n")

# Saves the solution into a VTK file

file = File("displacement.pvd");

file << u;

I haven’t found, though, anyone who evaluates the derivative analytically as I intend to do. So, I tested their approach of calculating the directional derivative using fenics’ derivative function. The Newton method with the latter technique converged much faster than with the analytical derivative. Thus, I tested my derivative and asked fenics to evaluate the fourth order elastic tensor numerically; the result was the same as with my analytical approach. Consequently, I’m clueless about this behavior. Is fenics encountering any issue with the whole linearized bilinear form I prescribed or what?

Edit: you can test yourself the behavior of the solution by simply choosing one of the three options for the bilinear form a.

14

I have changed my code and found a most interesting behavior:

```from dolfin import *

import ufl_legacy as ufl

from petsc4py import PETSc

i,j,k,l = ufl.indices(4)

import elastic_tensors

########################################################################
#                                 Mesh                                 #
########################################################################

# Defines the degree of the interpolation polynomial 

polynomial_degree = 1

element_type = 'P'

# Defines the mesh

mesh = UnitCubeMesh(16, 16, 24)#BoxMesh(Point(0,0,0), Point(1,1,1), 1,1,1)

# Defines the solution space using polynomial interpolation functions

V = VectorFunctionSpace(mesh, element_type, polynomial_degree)

########################################################################
#                        Newton-Raphson settings                       #
########################################################################

# Sets the final time

final_time = 1.0

# Sets the number of pseudotime points for loading

n_pseudoTimePoints = 1

# Sets the step in pseudo time

delta_t = final_time/n_pseudoTimePoints

# Defines the maximum number of iterations

n_maxIter = 40

# Defines a pseudotime variable to control the incremental process

t = 0.0

# Defines a function of the pseudo time that gives the displacement a-
# long the incremental process

def u_x(t):
    
    return t*0.0

def u_y(t):
    
    return t*0.0

def u_z(t):
    
    return t*0.0

# Defines the norm of the residual to be used as stopping criterion to 
# the Newton-Raphson method

residual_tol = 1E-4

########################################################################
#                    Dirichlet boundary conditions                     #
########################################################################

# Defines domain subregions

below =  CompiledSubDomain("near(x[2], side) && on_boundary", side=0.0)

above = CompiledSubDomain("near(x[2], side) && on_boundary", side=1.0)

# Defines the expressions for the Dirichlet boundary conditions, i.e.,
# their values for the displacement field

expression_below = Expression(("0.0", "0.0", "0.0"), degree=
polynomial_degree)

def expression_above(t):
    
    return Expression((str(u_x(t)), str(u_y(t)), str(u_z(t))), degree=
    polynomial_degree)

# Defines two objects of Dirichlet boundary conditions for the displace-
# ment field

bc_DirichletBelowDisplacement = DirichletBC(V, expression_below, below)

def bc_DirichletAboveDisplacement(t):
    
    return DirichletBC(V, expression_above(t), above)

def BCs_dirichletDisplacement(t):
    
    return [bc_DirichletBelowDisplacement, 
    bc_DirichletAboveDisplacement(t)]

# The update direction must be null in the regions where Dirichlet boun-
# dary conditions are applied to the displacement field

#bc_DirichletAboveDirection = DirichletBC(V, Constant(("0.0", "0.0",
#"0.0")), above)

bc_DirichletBelowDirection = DirichletBC(V, Constant(("0.0", "0.0",
"0.0")), below)

#BCs_dirichletDirection = [bc_DirichletBelowDirection,
#bc_DirichletAboveDirection]

BCs_dirichletDirection = [bc_DirichletBelowDirection]

########################################################################
#                     Neumann boundary conditions                      #
########################################################################

# Creates a mesh function for the boundaries to apply Neumann boundary
# conditions. The dimension of the integration over the boundary is one
# less than the domain integration. Besides, the last argument, 0, is to
# initialize the mesh with this value; so it substitutes the command 
# set_all(0)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

# Tags the labels for the subregions in the boundary

above.mark(boundaries, 1)

# Creates the surface area differential for the above facet

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Defines the boundary traction vector in the reference configuration

T = Constant((0.0, 0.0, 1000.0))

########################################################################
#                 Definition of solutions and variation                #
########################################################################

# Defines the function for the update direction in the Newton-Raphson 
# scheme

Delta_u = TrialFunction(V)

# Defines the variation

delta_u = TestFunction(V)

# Defines the displacement field

u = Function(V)

# Defines the body forces vector in the reference configuration

B = Constant((0.0, 0.0, 0.0))

########################################################################
#                          Tensors definition                          #
########################################################################

# Recovers the problem's dimension

dimension = Delta_u.geometric_dimension()

# Initializes the second order identity tensor

I = Identity(dimension)

# Calculates the deformation gradient tensor in terms of the solution u

F = grad(u)+I

# Calculates the right Cauchy-Green tensor

C = (F.T)*F

# Calculates the first two invariants of the right Cauchy-Green tensor

I1 = tr(C)

I2 = 0.5*((tr(C)**2)-tr(C*C))

# Calculates the jacobian of the kinematic mapping from the reference to
# the deformed configurations

J = det(F)

########################################################################
#                        Constitutive modelling                        #
########################################################################
 
# Saint Venant-Kirchhoff model

E = Constant(200E6)

nu = Constant(0.3) 

mu = E/(2*(1+nu))

lmbda = (E*nu)/((1+nu)*(1-2*nu))

# Defines the Helmholtz energy function

Psi = ((0.5*mu*(1-I2))+(0.25*mu*((I1-1)**2))+(0.125*lmbda*((I1-3)**2)))

# Defines the derivatives of the energy function w.r.t. the invariants 
# of the right Cauchy-Green tensor up to the second order

dPsi_dI1 = (0.5*mu*(I1-1))+(0.25*lmbda*(I1-3))

dPsi_dI2 = -0.5*mu

ddPsi_ddI1 = (0.5*mu)+(0.25*lmbda)

ddPsi_ddI2 = 0.0

ddPsi_dI1dI2 = 0.0

# Defines the k factor for the volumetric part of the second Piola-
# Kirchhoff stress tensor and the derivative of the former w.r.t. the
# jacobian

k_vol = 0.0

dkvol_dJ = 0.0

# Calculates the second Piola-Kirchhoff stress tensor

S = (2*(((dPsi_dI1+(I1*dPsi_dI2))*I)-(dPsi_dI2*C)+(k_vol*inv(C))))

# Defines the Cauchy stress tensor by the push-forward of the second Pi-
# ola-Kirchhoff stress tensor

sigma = (1/J)*F*S*F.T

# Evaluates the deviatoric parcel of the Cauchy stress tensor

dev_sigma = sigma-((1/3)*(tr(sigma))*I)

# Defines a function to evaluate the von Mises stress

von_mises = sqrt((3/2)*inner(dev_sigma,dev_sigma))

# Defines the fourth order elastic tensor

C_elastic = elastic_tensors.C_hyperNonIsochoric(C, I1, k_vol, J,
dkvol_dJ, dPsi_dI2, ddPsi_ddI1, ddPsi_ddI2, ddPsi_dI1dI2, dimension)

########################################################################
#               Variational form and solution procedure                #
########################################################################

# Defines the linearized variational form in terms of the update direc-
# tion, a==L. Note that a is constituted by the directional derivative 
# of the variation of the power/work; while L is formed by the residual
# itself. This scheme states the Newton-Raphson method into a variational
# form

L = (-inner((F.T)*grad(delta_u),S)*dx+(dot(B,delta_u)*dx)+(dot(T,delta_u
)*ds(1)))

# Analytical derivative of the fourth order elastic tensor

a = ((inner((F.T)*grad(delta_u),elastic_tensors.double_contraction(
C_elastic,sym((F.T)*grad(Delta_u))))*dx))+(inner(grad(delta_u)*S,grad(
Delta_u))*dx)

# Numerical derivative of the residue entire. Use it to validate the
# analytical derivative only

#a = -derivative(L, u, Delta_u)

Delta_u = Function(V)

# Iterates through the pseudotime points

for j in range(n_pseudoTimePoints):

    # Updates time

    t += delta_t

    print("\n#########################################################"+
    "\n#           Initiates the", j,"pseudotime point            #\n"+
    "#########################################################\n")

    print("t =", t)

    # Applies boundary conditions to the displacement field

    BCs_dirichletU = BCs_dirichletDisplacement(t)

    for CC in BCs_dirichletU:

        CC.apply(u.vector())

    # Iterates through the loop of the Newton-Raphson method

    for i in range(n_maxIter):

        print("\nInitiates iteration", i+1, "of the Newton-Raphson met"+
        "hod\n")

        # Solves the linearized variational form
        
        Kt, R = assemble_system(a, L, BCs_dirichletDirection)

        ### PETSc solver

        # Gets the tangent matrix

        Kt_matrix = as_backend_type(Kt).mat()

        vectorDelta_u, R_vector = Kt_matrix.getVecs()

        # Gets the residual vector

        R_vector.setValues(range(R.size()), R)

        # Sets the Krylov linear solver. Uses conjugate gradient ('cg')

        krylov_linearSolver = PETSc.KSP().create()

        krylov_linearSolver.setType('cg')

        krylov_linearSolver.setOperators(Kt_matrix)

        krylov_linearSolver.solve(R_vector, vectorDelta_u)

        # Evaluates the norm of the residual

        norm_residual = norm(R, "l2")

        print("Residual norm is", norm_residual)

        # Verifies stop criterion

        if norm_residual<residual_tol:

            print("Residual criterion met at iteration", i)

            break

        # Updates the displacement field
        
        u.vector()[:] = u.vector()[:]+vectorDelta_u

    print("\n\n")

    # Saves the solution into a VTK file

    file = File("displacement_t_"+str(j)+"_.pvd")

    file << u

# Calculates the von Mises stress

V = FunctionSpace(mesh, element_type, polynomial_degree)

von_mises = project(von_mises, V)

file_vonMises = File("von_misesStress.pvd")

file_vonMises << von_mises

After the comment “# Analytical derivative of the fourth order elastic tensor”, I put the linearized form. My Newton scheme converged very slowly and tried almost everything to fix it until, out of a stroke of luck, I decided to put the sym() for the F.T*grad(Delta_u) bit and it worked beautifully. Nevertheless, this behavior should not appear, for the fourth order tensor C has minor simmetry and, consequently, it operates with the symmetric part only of any second order tensor regardless. Therefore, I suspect it is a numerical error or some kind of a bug. Could you enlighten me with any idea for the origin of this behavior?