Arc-Length Implementation

Question

Dear all,

I am a new user of Fenics...

I am trying to implement the Arc-Length method inside of Fenics using python interface for nonlinear elasticity purpose.

I am finding some difficulties to assemble the updated tangent stiffness matrix to evaluate the displacement vector from the Nonlinear elasticity problem;

Here one example how I am defining my nonlinear elasticity problem:

def neoHookian(Mat_Par, u):
     [mu, lmbda] = Mat_Par;
     C = RightCauchyGreen(u);
     Ic = tr(C);     # Invariants of deformation tensors
     J = Jacobian(u);

     # Stored strain energy density (compressible neo-Hookean model of Simo-Ciarlet)
     psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

     return psi

# Elasticity parameters
E, nu = 3000.0, 0.4
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

phi = neoHookian([lmbda,mu],u);
Ui  = phi*dx;                                   #Internal Energy
Ues = - dot(g_B, u)*dx - dot(g_T, u)*ds(2);     #External Energy
Ut  = Ui + Ues;

dUt = derivative(Ut, u, v)
d2Ut = derivative(dUt, u, du)

And the I would like to get the assemble_system from the variable:

d2Ut

to get matrix and vector

A, b

to solve this Au = b

solve(A, u, b)

Has anyone some example to share?

Regards,
Ricardo

Comments

How do you want to extract b from d2Ut? I suppose it is a matrix. Have you tried to assemble it into a tensor (with assemble())?

Answer 1

Hello,

Ut == 0 or F == 0 that automatically solve it by using:

dUt = derivative(Ut, u, v)
d2Ut = derivative(dUt, u, du) #Jacobian
problem = NonlinearVariationalProblem(dUt, u, bcs, J=d2Ut);
solver  = NonlinearVariationalSolver(problem);

I would like to get a(u,v), L(v) from F directly (is it possible?) and get the following matrix and vectors:

K, f = assemble_system(a, L, bcs)      # Assembly the system

and then solve the linear system to obtain u vector, K u = f by using

solve(A, x, b)

Any idea?

Best Regards,

Comments

I think you don't have L(v) in your problem.
When you solve your defined problem as you described, you solve the equation: a(u,v) == 0. By assembling your form 'a' into a tensor: assemble(a,tensor=A) you can solve it with solve(A, x, b).

---

UFL has the functions lhs and rhs which will extract bilinear and linear parts resp. from a general form.

---

When you write:

problem = NonlinearVariationalProblem(dUt, u, bcs, J=d2Ut);

you solve dUt == 0 (meaning that your rhs = 0). You didn't define any rhs in your code above. What do you want to do exactly?

---

I would like to perform a perturbation analysis from F and J like this:

Ut = phi*dx - dot(g_B, u)*dx - dot(g_T, u)*ds(2);

F = derivative(Ut, u, v)

J = derivative(F, u, du)

K, f = suggested_assemble_system(J, bcs)

and then solving it (K u = f) :

solve(K, u, f)

I would like to implement a Newton method like this:

u_i+1 = u_i + K(u_i)\f(u_i)

each iteration I would like to update K, f from u_i to solve

---

I tried, but it not working get the lhs and rhs from F:

F = derivative(Ut, u, v)

I am going to prepare a example code to post it here

---

Here is my example code:

"""
Created on Fri Nov 27 13:47:19 2015

@author: Ricardo Lahuerta
"""

from dolfin import * 
from dolfin_adjoint import *
import numpy as np
#from fenicstools import *

# Optimization options for the form compiler
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
w = 1000 # 1000 mm comprimento, x
h = 250  # 250 mm largura, y

# Element Average Size
el_size = 10;   

nx = int(w/el_size);
ny = int(h/el_size);

# Mesh for cantilever beam
p0 = Point(0,0);
p1 = Point(w, h);

p_load = Point(w, h/2);

model = RectangleMesh(p0, p1, nx, ny)

# Compute all entities and connectivity.
model.init()

V = VectorFunctionSpace(model, "Lagrange", 2)

# Define Dirichlet boundary
# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class load_point(SubDomain):
    def inside(self, x, on_boundary):
        delta = 2*el_size
        return near(x[0],w) and (x[1] > h/2 - delta) and (x[1] < h/2 + delta) 

left = Left()
pt_load = load_point()

dim = model.topology().dim()-1;

sub_domains = MeshFunction("size_t", model, dim)

# Set all sub_domains in 0
sub_domains.set_all(0)

# Mark the sub_domain selection
left.mark(sub_domains, 1)
pt_load.mark(sub_domains, 2)

#plot(sub_domains_double)

ds = Measure("ds")[sub_domains]

# Elasticity parameters
E, nu = 3000.0, 0.4
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

Load_Area = 2*el_size;

force_x = 0;
force_y = -500;

press_x = force_x/Load_Area;
press_y = force_y/Load_Area;

g_B = Constant((0.0, 0.0))            # Body force per unit volume
g_T = Constant((press_x, press_y))    # Traction force on the boundary

# Define functions
du = TrialFunction(V)                   # Incremental displacement
v  = TestFunction(V)                    # Test function
u  = Function(V, name = "Displacement")   # Displacement from previous iteration 

#Kinematics
I = Identity(u.geometric_dimension());  # Identity tensor
F = variable(I + grad(u));              # Deformation Gradient
J = det(F);
C = variable(F.T*F);                # Right Cauchy-Green tensor
Ic = tr(C);     # Invariants of deformation tensors

# Stored strain energy density (compressible neo-Hookean)
phi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2;

Ui  = phi*dx;                                   #Internal Energy
Ues = - dot(g_B, u)*dx - dot(g_T, u)*ds(2);     #External Energy, only surface loads
Ut  = Ui + Ues;                                 #Equilibrium, Total potential energy

# Compute first variation of Ut (directional derivative about u in the direction of v)
#dUtdu = derivative(Ut, u, du)

# Compute Jacobian of dUtdu, "Tangent Tensor"
#d2Ut = derivative(dUtdu, u, du)

P = diff(phi,F);

a = inner(P,grad(v));
L = inner(g_B,v)*dx + inner(g_T,v)*ds(2);

# Apply the boundary condition
left_end_displacement = Constant((0.0, 0.0));

bcs = [ DirichletBC(V, left_end_displacement, left) ];

# Assembly the system
A, b = assemble_system(a, L, bcs)

solve(A, u, b);

Why fenics is not able to assembly a and L?

---

In this code u has to be TrialFunction instead of Function.

---

I tried, but it is not working, this is the error message:

Traceback (most recent call last):
File "", line 1, in
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 243, in assemble_system
A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
raise TypeError("Invalid form type %s" % (type(form),))
TypeError: Invalid form type