ffc team mailing list archive

Thread
Date

Nonlinear elasticity

To: ffc-dev@xxxxxxxxxx
From: Marco Morandini <morandini@xxxxxxxxxxxxxx>
Date: Tue, 21 Feb 2006 10:09:23 +0100
User-agent: Thunderbird 1.5 (X11/20051201)

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 thesymmetries of E.



Thanks in advance,

Marco

# Copyright (c) 2005 Johan Jansson (johanjan@xxxxxxxxxxxxxxxx)
# 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 (johanjan@xxxxxxxxxxxxxxxx)
# 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 (johanjan@xxxxxxxxxxxxxxxx)
# 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 (johanjan@xxxxxxxxxxxxxxxx)
# 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'

Follow ups

Re: Nonlinear elasticity
From: Anders Logg, 2006-02-21