ffc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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


On Tue, Feb 21, 2006 at 10:09:23AM +0100, Marco Morandini wrote:
> 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 (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'

> _______________________________________________
> FFC-dev mailing list
> FFC-dev@xxxxxxxxxx
> 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/