Tensor representation

Please provide a minimal example of what is not working to improve the
chances of getting help :)

Johan

On Mon, Mar 5, 2012 at 11:20 AM, Dupront
<email address hidden> wrote:
>
> New question #189712 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/189712
>
> Hello,
>
> Is the tensor representation supposed to work for any
> kinds of model ? I tried with different examples. It works
> for poisson, it does not work for linear elasticity.
>
> Thank you
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

This is the official python elasticity example in which I added the line:
  parameters["form_compiler"]["representation"] = "quadrature"
If I use "tensor" instead of "quadrature it does not work.
I am using Fenics on a ubuntu lucid platform.

"""This demo program solves the equations of static linear elasticity
for a gear clamped at two of its ends and twisted 30 degrees."""

__author__ = "Kristian B. Oelgaard (<email address hidden>)"
__date__ = "2007-11-14 -- 2009-10-07"
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
__license__ = "GNU LGPL Version 2.1"

# Modified by Anders Logg, 2008.

from dolfin import *

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"

# Load mesh and define function space
mesh = Mesh("gear.xml.gz")
mesh.order()
V = VectorFunctionSpace(mesh, "CG", 1)

# Sub domain for clamp at left end
def left(x, on_boundary):
    return x[0] < 0.5 and on_boundary

# Dirichlet boundary condition for rotation at right end
class Rotation(Expression):

    def eval(self, values, x):

        # Center of rotation
        y0 = 0.5
        z0 = 0.219

        # Angle of rotation (30 degrees)
        theta = 0.5236

        # New coordinates
        y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta)
        z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta)

        # Clamp at right end
        values[0] = 0.0
        values[1] = y - x[1]
        values[2] = z - x[2]

    def dim(self):
        return 3

# Sub domain for rotation at right end
def right(x, on_boundary):
    return x[0] > 0.9 and on_boundary

# Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
f = Constant((0.0, 0.0, 0.0))

E = 10.0
nu = 0.3

mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)

a = inner(grad(v), sigma(u))*dx
L = inner(v, f)*dx

# Set up boundary condition at left end
c = Constant((0.0, 0.0, 0.0))
bcl = DirichletBC(V, c, left)

# Set up boundary condition at right end
r = Rotation()
bcr = DirichletBC(V, r, right)

# Set up boundary conditions
bcs = [bcl, bcr]

# Set up PDE and solve
problem = VariationalProblem(a, L, bcs)
problem.parameters["symmetric"] = True
u = problem.solve()

# Save solution to VTK format
vtk_file = File("elasticity.pvd")
vtk_file << u

# Plot solution
#plot(u, mode="displacement", interactive=True)

info(parameters,True)

print "\n\n Times:"
summary()

*********************************
* Error Message
*********************************
Got expression dimension = 3
Automatic selection of expression element: <Lagrange vector element of degree ? on a <? of degree 1>: 3 x <CG? on a <? of degree 1>>>
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "demo.py", line 85, in <module>
    problem = VariationalProblem(a, L, bcs)
  File "/usr/lib/python2.6/dist-packages/dolfin/fem/variationalproblem.py", line 39, in __init__
    self.a = Form(a, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python2.6/dist-packages/dolfin/fem/form.py", line 34, in __init__
    (self._compiled_form, module, self.form_data) = jit(form, form_compiler_parameters, common_cell)
  File "/usr/lib/python2.6/dist-packages/dolfin/compilemodules/jit.py", line 44, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/python2.6/dist-packages/dolfin/compilemodules/jit.py", line 103, in jit
    return jit_compile(form, parameters=p, common_cell=common_cell)
  File "/usr/lib/python2.6/dist-packages/ffc/jitcompiler.py", line 58, in jit
    return jit_form(object, parameters, common_cell)
  File "/usr/lib/python2.6/dist-packages/ffc/jitcompiler.py", line 95, in jit_form
    compile_form(preprocessed_form, prefix=jit_object.signature(), parameters=parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/compiler.py", line 140, in compile_form
    ir = compute_ir(analysis, parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/representation.py", line 65, in compute_ir
    irs = [_compute_integral_ir(f, i, parameters) for (i, f) in enumerate(forms)]
  File "/usr/lib/python2.6/dist-packages/ffc/representation.py", line 162, in _compute_integral_ir
    parameters)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/tensorrepresentation.py", line 39, in compute_integral_ir
    monomial_form = extract_monomial_form(integrals)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/monomialextraction.py", line 41, in extract_monomial_form
    integrand = extract_monomial_integrand(integrand)
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/monomialextraction.py", line 58, in extract_monomial_integrand
    monomial_integrand = apply_transformer(integrand, MonomialTransformer())
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 782, in apply_transformer
    return transform_integrands(e, _transform)
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 773, in transform_integrands
    return transform(expr)
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 781, in _transform
    return transformer.visit(expr)
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ufl/algorithms/transformations.py", line 129, in visit
    r = h(o, *map(self.visit, o.operands()))
  File "/usr/lib/python2.6/dist-packages/ffc/tensor/monomialextraction.py", line 271, in expr
    raise MonomialException, ("No handler defined for expression %s." % o._uflclass.__name__)
ffc.tensor.monomialextraction.MonomialException: No handler defined for expression Division.

The tensor representation doesn't fully support division, which is a problem for the symmetric gradient (a fix for this was pushed a few days ago to the development version of FFC) and it does support the identity matrix. This is why it doesn't work for the elasticity demo.

Thank you.

For the time being I have to use the lucid version of Fenics.
Is it a huge modification to make, that would require a lot of
effort or it is a very easy local one ?

lucid is a supported platform for the dorsal script installer. Have a look at:

  http://fenicsproject.org/download/installation_using_dorsal.html#installation-using-dorsal

Johan

On Tue, Mar 6, 2012 at 2:05 PM, Dupront
<email address hidden> wrote:
> Question #189712 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/189712
>
>    Status: Answered => Open
>
> Dupront is still having a problem:
> Thank you.
>
> For the time being I have to use the lucid version of Fenics.
> Is it a huge modification to make, that would require a lot of
> effort or it is a very easy local one ?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Thank a lot.

I guess I could not avoid it any longer.
So I gave it a try.

Dorsal was able to finish a built even if four libraries (prerequisites) were missing:
 libcgal-dev
 libcln-dev
 libginac-dev
 libgmp3-dev

I am still trying to install them.
I have dependency problems.
For example when I try to install libgmp3-dev with synaptic,
I get the following message:
 " libgmpxx4ldbl (=2:4.3.2+dfsg-1ubuntu1) but 2:5.0.2+dfsg-0ubuntu1 has to be installed ".
But I cannot find a repositary that contains 2:5.0.2+dfsg-0ubuntu1.

How important are these libraries ?

I tried few demo examples and they run fine.
Also few tests in tist from the last version of cbc.solve are working fine.

None of those libraries are important for DOLFIN. You probably want to have CGAL, but libcgal-dev in Lucid is too old for DOLFIN, so this must be built by Dorsal anyway. This is the default in the development version of Dorsal (bzr branch lp:dorsal).

BTW: There is also a PPA (ppa:fenics-packages/fenics-dev) with packages from the development version of DOLFIN.

Thank you !

I have a 1.0.0 version that works pretty well.

But I am confused now with Mr Garth Wells comment.
I understand that a modification has been made to make
the tensor representation for elastic problems works.
But in which release is it available ? The snapshot one ?

On 7 March 2012 22:15, Dupront <email address hidden> wrote:
> Question #189712 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/189712
>
>    Status: Answered => Open
>
> Dupront is still having a problem:
> Thank you !
>
> I have a 1.0.0 version that works pretty  well.
>
> But I am confused now with Mr Garth Wells comment.
> I understand that a modification has been made to make
> the tensor representation for elastic problems works.

That's not what I said.

If you get the latest ffc (via 'bzr branch ffc") you will get support
for fractions like 1.0/2.0 with the tensor representation, which deals
with the symmetric gradient. However, the elasticity demo also uses
the identity matrix, which is not supported for the tensor
representation.

Garth

> But in which release is it available ? The snapshot one ?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

It is clear now.
Thank you very for your kindness.