Subject: [FFC-dev] Nonlinear elasticity

From: Marco Morandini <[email-masked]>
Date: Tue, 21 Feb 2006 10:09:23 +0100

With my classical structural FEM background,
my interpretation of "multilinear form"
is the linearized equation of the weks form, ready to be used
for a Newton iteration. Please correct me
if I'm wrong.

I'm interested in the equation of nonlinear elasticity.
I've come up with the attached forms, but ffc
requires more than 800 Mb for NonlinElasticity.form and
NonlinElasticity1.form, and 1200 Mb for NonlinElasticity2.form
and NonlinElasticity3.form, and I have to kill it
(I only have 500 Mb of RAM).

The question is: am I doing something wrong?
Is there a way to reduce the amount of memory required by ffc
for this problem?

Problem description:

weak form:
\int \delta \epsilon : \sigma dV = \int \delta x' f dV

where

x' is the current configuration of the body
\epsilon = 0.5(transp(F)*F -I) is the deformation tensor

F = grad(x')

(
alternatively, as in NonlinElasticity1.form,
the deformation is computed as
\epsilon = 0.5(transp(grad(u)) + grad(u) + transp(grad(u))*grad(u))
with u = x' - x the displacement from the undeformed position x
)

\sigma = E : \epsilon is the stress tensor (E is a fourth-order tensor)

In the .form files:

u: unknown increment of x'
    (or of displacement in NonlinElasticity1.form)
v: test functions for current position x'
w: current position (current x') from the previous Newton iteration

In NonlinElasticity2.form I'm assuming that the current stress \sigma
and elastic tensor E are computed elsewhere. In NonlinElasticiy3.form
I'm trying to exploit the symmetry of \sigma, \epsilon, and some of the 
symmetries of E.

Thanks in advance,

Marco
# Copyright (c) 2005 Johan Jansson ([email-masked])
# Licensed under the GNU GPL Version 2
#
# The bilinear form for classical linear elasticity (Navier)
# Compile this form with FFC: ffc Elasticity.form.

element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)

v = BasisFunction(element)
u = BasisFunction(element)
w = Function(element)

lmbda = Constant() # Lame coefficient
mu    = Constant() # Lame coefficient

f = Function(element) # Source

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

epsilon = 0.5 * (transp(F) * F - Identity(d) )

delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)

de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)

def E(e, lmbda, mu):
    Ee = 2.0 * mult(mu, e) + mult(lmbda, mult(trace(e), Identity(d)))
    
    return Ee

sigma = E(epsilon, lmbda, mu)

desigma = E(de_epsilon, lmbda, mu)

#def depsilon(v):
#    return 0.5 * (grad(v) + transp(grad(v)))

a = (0.25 * trace(deltaF * sigma * transp(deF) + deF * sigma * transp(deltaF) ) + \
    dot(desigma, delta_epsilon)) * dx
#a =  \
#    dot(desigma, delta_epsilon) * dx
L = (f[i] * v[i] - dot(sigma, de_epsilon) ) * dx 
# Copyright (c) 2005 Johan Jansson ([email-masked])
# Licensed under the GNU GPL Version 2
#
# The bilinear form for classical linear elasticity (Navier)
# Compile this form with FFC: ffc Elasticity.form.

element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)

v = BasisFunction(element)
u = BasisFunction(element)
w = Function(element)

lmbda = Constant() # Lame coefficient
mu    = Constant() # Lame coefficient

f = Function(element) # Source

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

#epsilon = 0.5 * (transp(F) * F - Identity(d) )
epsilon = 0.5 * (transp(F) + F + transp(F)*F )

#delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)
delta_epsilon = 0.5 * (transp(deltaF) + deltaF + transp(deltaF)*F + transp(F)*deltaF)

#de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)
de_epsilon = 0.5 * (transp(deF) + deF + transp(deF)*F + transp(F)*deF)

def E(e, lmbda, mu):
    Ee = 2.0 * mult(mu, e) + mult(lmbda, mult(trace(e), Identity(d)))
    
    return Ee

sigma = E(epsilon, lmbda, mu)

desigma = E(de_epsilon, lmbda, mu)

#def depsilon(v):
#    return 0.5 * (grad(v) + transp(grad(v)))

a = (0.25 * trace(deltaF * sigma * transp(deF) + deF * sigma * transp(deltaF) ) + \
    dot(desigma, delta_epsilon)) * dx
#a =  \
#    dot(desigma, delta_epsilon) * dx
L = (f[i] * v[i] - dot(sigma, de_epsilon) ) * dx 
# Copyright (c) 2005 Johan Jansson ([email-masked])
# Licensed under the GNU GPL Version 2
#
# The bilinear form for classical linear elasticity (Navier)
# Compile this form with FFC: ffc Elasticity.form.

from Numeric import *
from math import *

element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)
element2 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 81)

v = BasisFunction(element)
u = BasisFunction(element)
w = Function(element)
sigma = Function(element1)
desigma = Function(element2)

def tomatrix(q):
    return [ [q[3 * i + j] for i in range(3)] for j in range(3) ]

def to4tensor(q):
    return [ [ [ [q[ 3 * (3 * (3 * i + j) + k) + l] for i in range(3)] \
        for j in range(3) ] \
            for k in range(3) ] \
                for l in range(3) ]

sigmam = tomatrix(sigma)
desigmam = to4tensor(desigma)

lmbda = Constant() # Lame coefficient
mu    = Constant() # Lame coefficient

f = Function(element) # Source

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

epsilon = 0.5 * (transp(F) * F - Identity(d) )

delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)

de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)

#def E(e, lmbda, mu):
#    Ee = 2.0 * mult(mu, e) + mult(lmbda, mult(trace(e), Identity(d)))
#    
#    return Ee

#desigma = E(de_epsilon, lmbda, mu)

#def depsilon(v):
#    return 0.5 * (grad(v) + transp(grad(v)))

#deltaF = tomatrix(deltaF)
#deF = tomatrix(deF)
#de_epsilon = tomatrix(de_epsilon)

dummy = 0
for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        print i,j,k,l
        dummy += (0.25 * (deltaF[l][i] * sigmam[i][j] * deF[l][j] + deF[l][i] * sigmam[i][j] * deltaF[l][j] ) + \
            desigmam[i][j][k][l] * delta_epsilon[i][j] * de_epsilon[k][l] )

print 'urca'
a = dummy * dx;
#a =  \
#    dot(desigma, delta_epsilon) * dx
dummy = 0
for i in range(3):
  for j in range(3):
    print i,j
    dummy += (sigmam[i][j] * de_epsilon[i][j])

L = (f[i] * v[i] + dummy ) * dx 
print 'urcala'
# Copyright (c) 2005 Johan Jansson ([email-masked])
# Licensed under the GNU GPL Version 2
#
# The bilinear form for classical linear elasticity (Navier)
# Compile this form with FFC: ffc Elasticity.form.

from Numeric import *
from math import *

element = FiniteElement("Vector Lagrange", "tetrahedron", 1)
element1 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 9)
element2 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0, 81)

v = BasisFunction(element)
u = BasisFunction(element)
w = Function(element)
sigma = Function(element1)
desigma = Function(element2)

def tomatrix(q):
    return [ [q[3 * i + j] for i in range(3)] for j in range(3) ]

def to4tensor(q):
    return [ [ [ [q[ 3 * (3 * (3 * i + j) + i) + j] for i in range(3)] \
        for j in range(3) ] \
            for k in range(3) ] \
                for l in range(3) ]

sigmam = tomatrix(sigma)
for i in range(3):
  for j in range(i+1,3):
    sigmam[j][i] = sigmam [i][j]

desigmam = to4tensor(desigma)

for i in range(3):
  for j in range(i+1,3):
    for k in range(3):
      for l in range(k+1,3):
        desigmam[j][i][l][k] = desigmam[i][j][k][l]

f = Function(element) # Source

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

epsilon = 0.5 * (transp(F) * F - Identity(d) )

delta_epsilon = 0.5 * (transp(deltaF)*F + transp(F)*deltaF)

de_epsilon = 0.5 * (transp(deF)*F + transp(F)*deF)

#def E(e, lmbda, mu):
#    Ee = 2.0 * mult(mu, e) + mult(lmbda, mult(trace(e), Identity(d)))
#    
#    return Ee

#desigma = E(de_epsilon, lmbda, mu)

#def depsilon(v):
#    return 0.5 * (grad(v) + transp(grad(v)))

#deltaF = tomatrix(deltaF)
#deF = tomatrix(deF)
#de_epsilon = tomatrix(de_epsilon)

dummy = 0
for i in range(3):
  for j in range(i,3):
    if ( i != j ):
       mul1 = 2.
    else:
       mul1 = 1.
    for k in range(3):
      for l in range(k,3):
        if ( k != l ):
           mul2 = 2.
        else:
           mul2 = 1.
        print i,j,k,l
        dummy += (0.25 * (deltaF[l][i] * sigmam[i][j] * deF[l][j] + deF[l][i] * sigmam[i][j] * deltaF[l][j] ) + \
            desigmam[i][j][k][l] * delta_epsilon[i][j] * de_epsilon[k][l] ) * mul1 * mul2

print 'urca'
a = dummy * dx;
#a =  \
#    dot(desigma, delta_epsilon) * dx
dummy = 0
for i in range(3):
  for j in range(i,3):
    if ( i != j ):
       mul1 = 2.
    else:
       mul1 = 1.
    print i,j
    dummy += (sigmam[i][j] * de_epsilon[i][j]) * mul1

L = (f[i] * v[i] + dummy ) * dx 
print 'urcala'

---

From: Anders Logg <[email-masked]>
Date: Tue, 21 Feb 2006 13:32:18 -0600

I tried to compile the forms and FFC does indeed consume a lot of
memory. This is a result of the forms having many coefficients (like
F^2) and derivatives and the approach used to evaluate the forms in
FFC (by a special tensor representation) may not be suitable in this
case. A straightforward quadrature-based evaluation may be more
suitable, but it will be a while before I get to implementing it.

In the meanwhile, try isolating the term that takes longest to compile
and consumes the most memory (by removing terms one by one) and maybe
we can find a suitable work-around for it. It's easier to spot the
problem if the form is a simple as possible.

/Anders

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/cgi-bin/mailman/listinfo/ffc-dev

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/

---

From: Marco Morandini <[email-masked]>
Date: Thu, 23 Feb 2006 16:16:07 +0100

> problem if the form is a simple as possible.

On an Athlon XP 2400+ :

NonlinElasticitys1.form:  5m43.068s compilation time, ~51 Mb
NonlinElasticitys2.form: 28m58.479s compilation time, ~244 Mb
NonlinElasticitys3.form: (no compilation time, sorry),~473 Mb

Let me know if I can do something else

Thanks,

Marco

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

def E(e, lmbda, mu):
    Ee = 2.0 * mult(mu, e)
    
    return Ee

a1 = 0.25 * trace(deltaF * sigma * transp(deF) )
#time a1: user    5m43.068s
#memory: ~51 Mb
#
a = a1 * dx
L1 = f[i] * v[i] 
L = (L1) * dx 

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

def E(e, lmbda, mu):
    Ee = 2.0 * mult(mu, e)
    
    return Ee

a3 = dot(desigma, delta_epsilon)
#time a3: user    28m58.479s
#memory: ~244 Mb
a = a3 * dx
L1 = f[i] * v[i] 
L = (L1) * dx 
# Copyright (c) 2005 Johan Jansson ([email-masked])
# Licensed under the GNU GPL Version 2
#
# The bilinear form for classical linear elasticity (Navier)
# Compile this form with FFC: ffc Elasticity.form.

# Dimension of domain
d = element.shapedim()

deltaF =  grad(v)

deF =  grad(u)

F = grad(w)

def E(e, lmbda, mu):
    Ee = mult(lmbda, mult(trace(e), Identity(d))) #Ee2
    
    return Ee

a3 = dot(desigma, delta_epsilon)
#time a3: 
#memory: ~473 Mb
a = a3 * dx
L1 = f[i] * v[i] 
L = (L1) * dx

---

From: Anders Logg <[email-masked]>
Date: Mon, 27 Feb 2006 14:19:19 -0600

Hi!

I checked out the last of the form files, NonlinElasticitys3.form,
and as you say it takes a long time to compile (I didn't wait).

I thought about it and couldn't come up with a quick short-cut, but
here's a work-around until I figure out how to do this in general.

The main problem seems to be to compile a form of the following type:

#------------------------------------------------------------
element = FiniteElement("Vector Lagrange", "tetrahedron", 1)

F = grad(w)

a = dot(transp(F)*grad(v), transp(F)*grad(u))*dx
#------------------------------------------------------------

This is a stripped-down version of the form you are working with, just
keeping the essential terms that are making trouble. The reason it
takes long is that you have a product of four factors and lots of
derivatives (four in total) so the tensors that get computed have high
rank.

To speed this up, we need to make sure that all function spaces have
as low dimensionality as possible, so we try to compute with piecewise
constants where possible.

Let's assume for a minute that you have already computed the gradients
of the components in w (which are piecewise constants). Then you can
define the following simplified form:

#------------------------------------------------------------
P1 = FiniteElement("Vector Lagrange", "tetrahedron", 1)
P0 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0)

v = BasisFunction(P1)
u = BasisFunction(P1)

F0 = Function(P0) # F0 = grad(w[0])
F1 = Function(P0) # F1 = grad(w[1])
F2 = Function(P0) # F2 = grad(w[2])

F = [vec(F0), vec(F1), vec(F2)] # put them in a matrix

a = dot(transp(F)*grad(v), transp(F)*grad(u))*dx
#------------------------------------------------------------

Now, this form compiles in just 1-2 seconds instead of forever.
All that remains is now to precompute the gradients of the components
of w. The following form does the trick for the first component (so
you need to have two more for the other components):

#------------------------------------------------------------
# Compute projection of gradient (first component)

P1 = FiniteElement("Vector Lagrange", "tetrahedron", 1)
P0 = FiniteElement("Discontinuous vector Lagrange", "tetrahedron", 0)

v = BasisFunction(P0)
u = BasisFunction(P0)
w = Function(P1)

a = dot(v, u)*dx
L = dot(v, grad(w[0]))*dx
#------------------------------------------------------------

Hope this helps.

/Anders