This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Convergence issue with quadratic Lagrange element

+1 vote

Hi,

I am trying to solve a finite elasticity problem using a Neo-hookean material in which a bar gets stiffer with time and has inhomogeneous stiffness in the domain (in the x-direction). I am applying constraints (1) ux = 0 on the left and right side of the bar (x = 0 and x = 10), (2) uy = 0 on the y = 1 plane and (3) uz = 0 on the z = 1 plane.

I am applying the inhomogenity by modifying the stiffness at the gauss points by defining a the stiffness as a "Quadrature" function and then modifying contents of the array (see the variable "Tact").

For some reason, the solution does not converged when I use quadratic element and when the inhomogeneity gets large (given by the parameter "stiffness_diff"). However, the solution converge without any issue when a linear element is used (even with very inhomogeneous stiffness)

Any idea of what is happening? Have attached the code - any help is greatly appreciated! Thanks!!

from dolfin import *
import numpy as np
import math as math

quad_degree = 2

# Optimization options for the form compiler
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
           "eliminate_zeros": True, \
           "precompute_basis_const": True, \
           "precompute_ip_const": True}

# Create mesh and define function space
L = 10.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(L, 1.0, 1.0), 10, 2, 2)
dx = dolfin.dx(mesh, metadata = {"integration_order":quad_degree})

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = L)
top = CompiledSubDomain("near(x[2], side) && on_boundary", side = 1.0)
bot = CompiledSubDomain("near(x[2], side) && on_boundary", side = 0.0)
front = CompiledSubDomain("near(x[1], side) && on_boundary", side = 1.0)
back = CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.0)

facetboundaries = FacetFunction("size_t", mesh)     
facetboundaries.set_all(0) 
right.mark(facetboundaries, 2)
left.mark(facetboundaries, 1)
top.mark(facetboundaries, 3)
front.mark(facetboundaries, 4)

W = VectorFunctionSpace(mesh, "Lagrange", 1)
Quads = FunctionSpace(mesh, "Quadrature", quad_degree)

bcright = DirichletBC(W.sub(0), Constant(("0")), facetboundaries, 2)
bctop = DirichletBC(W.sub(2), Constant(("0")), facetboundaries, 3)
bcfront = DirichletBC(W.sub(1), Constant(("0")), facetboundaries, 4)
bcleft = DirichletBC(W.sub(0), Constant(("0")), facetboundaries, 1)

bcs = [bcleft, bctop, bcfront, bcright]

# Define functions
du = TrialFunction(W)            # Incremental displacement
v  = TestFunction(W)             # Test function
u  = Function(W)                 # Displacement from previous iteration
Tact = Function(Quads)

d = u.geometric_dimension()
I = Identity(d)             
F = I + grad(u)             
J = det(F)
Cmat = F.T*F            
Fdev = J**(-1.0/3.0)*F
Cdev = Fdev.T*Fdev
Kappa = Constant(1.0e5);
f0 = Constant(("1.0", "0.0", "0.0"))

W1 = Constant(1000.0)*(tr(Cdev) - Constant(3.0)) + Kappa/Constant(2.0)*(J - 1.0)**2.0
W2 = Tact/2.0 * (Cmat[0,0] - 1.0)
Wtotal = (W1 + W2)

# Total potential energy
Pi = Wtotal*dx 

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
Jac = derivative(F, u, du)

t_niter = 40
r_threshold = 1e-9

# Solve variational problem
dt = 2.0
t = 0.0
Tmax = 80e3
Tend = 100.0;
ngpt = len(Tact.vector().array()[:])
t_niter = 40
r_threshold = 1e-9

# Varying the inhomogenity of the stiffness
stiffness_diff = 50000

disp_vtk_file = File("nonlinearelasticity_block_disp_diag.pvd")
while(t < Tend):


# Modifying stiffness at the Gauss points
Tact.vector()[:] = np.ones(ngpt)*sin(t/Tend*math.pi)*80e3 + np.random.rand(1,ngpt)[0]*stiffness_diff
print "Max stiffness = ", max(Tact.vector().array()),  " Min stiffness = ", min(Tact.vector().array())

solve(F == 0, u, bcs, J=Jac,form_compiler_parameters=ffc_options)

t = t + dt
disp_vtk_file << (u, t)

# Plot and hold solution
plot(u, interactive=True, mode="displacement")
asked Apr 29, 2016 by lee FEniCS User (1,170 points)
edited Apr 29, 2016 by lee
...