Fourth order tensor cpp expression

Question

Hi,

Does anybody know how to define a 4th order tensor expression in with cpp code? For some reason, the following does not work:

import dolfin
import numpy as np

mesh = dolfin.UnitCubeMesh(1, 1, 1)

dx = dolfin.Measure("dx", domain=mesh)

code2 = '''

namespace dolfin
{

class MyFunc2 : public Expression
{

  mutable Array<double> UX;

public:

  std::shared_ptr<MeshFunction<std::size_t> > cell_data;
  std::shared_ptr<const Function> U;
  std::vector<std::size_t> shape = {3,3,3,3};

  MyFunc2() : Expression(shape), UX(2)
  {
  }

  void eval(Array<double>& values, const Array<double>& x) const
  {

    values[0] = 1;
    values[1] = 2;
    values[2] = 3;

    values[3] = 4;
    values[4] = 1;
    values[5] = 6;

    values[6] = 7;
    values[7] = 8;
    values[8] = 9;


  }
};
}'''


Tensor4 = dolfin.Expression(code2)

Answer 1

Hi, consider

import dolfin
import numpy as np

mesh = dolfin.UnitCubeMesh(1, 1, 1)

dx = dolfin.Measure("dx", domain=mesh)

code2 = '''

namespace dolfin
{

  std::vector<std::size_t> shape = {3,3,3,3};

class MyFunc2 : public Expression
{

  mutable Array<double> UX;

public:

  std::shared_ptr<MeshFunction<std::size_t> > cell_data;
  std::shared_ptr<const Function> U;

  MyFunc2() : Expression(shape), UX(2)
  {
  }

  void eval(Array<double>& values, const Array<double>& x) const
  {
    for(int i=0; i<81;i++){
        values[i] = i;
    }
  }
};
}'''


Tensor4 = dolfin.Expression(code2, degree=1)
assert Tensor4.ufl_shape == (3, 3, 3, 3)
assert max(np.abs(Tensor4(0, 0, 0) - np.arange(81))) < 1E-14

Comments

Thanks MiroK !!

I was able to run this. One more question - how do I a double contraction of the 4th order tensor expression with a 2nd order tensor?

I tried adding these two lines:

Idty = dolfin.Identity(3)
Tensor4*Idty

but it doesn't seems to work with the following error:

Invalid ranks 4 and 2 in product.
Traceback (most recent call last):
File "testfenics4thorder.py", line 43, in
Tensor4Idty
File "/home/likchuan/Research/fenics/hashstack/fenics/lib/python2.7/site-packages/ufl/exproperators.py", line 164, in _mul
return _mult(self, o)
File "/home/likchuan/Research/fenics/hashstack/fenics/lib/python2.7/site-packages/ufl/exproperators.py", line 142, in _mult
error("Invalid ranks {0} and {1} in product.".format(r1, r2))
File "/home/likchuan/Research/fenics/hashstack/fenics/lib/python2.7/site-packages/ufl/log.py", line 151, in error
raise self._exception_type(self._format_raw(message))
ufl.log.UFLException: Invalid ranks 4 and 2 in product.

---

The following works with the code I posted in the answer.

from dolfin import *

I = Identity(3)
foo = Tensor4[i, j, k, l]*I[k, l]

Please refer to the UFL manual to see how to use the 'asterisk' operator.

---

Many thanks, Mirok!!

One more question - I would like to pass an ufl form example grad(u) in which u is a vector function into the expression. Is there a way to do that? I tried the following and it doesn't work.

The reason for doing so is that we want to compute a fourth order elasticity tensor that depends on the strain from a displacement field to form the Jacobian for the Newton Solver.

import dolfin
import numpy as np

mesh = dolfin.UnitCubeMesh(1, 1, 1)
S = dolfin.FunctionSpace(mesh, "CG", 2)
u = dolfin.Function(S)
gradu = dolfin.grad(u)

dx = dolfin.Measure("dx", domain=mesh)

code2 = '''

namespace dolfin
{

  std::vector<std::size_t> shape = {3,3,3,3};

class MyFunc2 : public Expression
{

  mutable Array<double> UX;

public:

  std::shared_ptr<MeshFunction<std::size_t> > cell_data;
  std::shared_ptr<const Function> U;

  MyFunc2() : Expression(shape), UX(2)
  {
  }

  void eval(Array<double>& values, const Array<double>& x) const
  {
    for(int i=0; i<81;i++){
        values[i] = i;
    }
  }
};
}'''


Tensor4 = dolfin.Expression(code2, degree=1)
Tensor4.U = gradu

---

lee - please open a new question instead of asking questions in the comments.

---

Thanks, johannr.
I have opened up a new question in https://fenicsproject.org/qa/11118/how-to-pass-a-ufl-form-e-g-grad-u-into-a-cpp-expression