Projection of 4th order Tensor

164
views
1
4 months ago by
Hello everyone,

I'm trying to project 4th order tensor on finite element space. Is there any way to do that ? Is it possible to do it with Quadrature finite element space ? I tried to create TensorFunctionSpace but it did not work. I think it only works with second order tensors.

V = TensorFunctionSpace(mesh, 'CG', 1, shape=(3,3,3,3), symmetry=True)​
Thanks.
Community: FEniCS Project

5
4 months ago by
The problem here is symmetry. There are more notions of symmetry for tensors with rank 4 as opposed to rank 2 tensor where there is only one notorious symmetry T_{ij} = T_{ji} in flat euclidean space.

Try to do
elem = TensorElement("CG", mesh.ufl_cell(), shape=(3,3), symmetry=True)
elem.symmetry()​
and you'll see how is symmetry mapping constructed for (3, 3) tensor. It gives {(2, 0): (0, 2), (1, 0): (0, 1), (2, 1): (1, 2)}.
So prepare your own "symmetry" element as, e.g.
sym = {(1, 0, 0): (0, 1, 0)}
elem = TensorElement("CG", mesh.ufl_cell(), 1, shape=(3, 3, 3), symmetry=sym)​

Warning: This should be well supported in UFL codebase, but you'll encounter some (many) problems in DOLFIN (IO, plotting, ...)
Hello Michal,

Thank you for your advice. That's true. But lets say that, I do not give any information about the symmetry and keep that region empty. Then when I am trying to project by using the 4th order function space, I got error as following :
from dolfin import *
from mshr import *
import numpy as np
import warnings

mesh = BoxMesh(Point(0, 0, 0), Point(4, 2, 2), 4, 2, 2)

V = TensorFunctionSpace(mesh, 'DG', 0, shape=(3,3,3,3))

I = Identity(3)
def delta(i,j):
elta=I[i,j]
return elta

I4=delta(i,j)*delta(k,l)

C = project(I4,V)​
​

I simply used to project 4th order identity tensor over the 4th order TensorFunctionSpace. But it gives an error that the shapes do not match. However for both I4 and V the shape is (3,3,3,3). What can be reason of that error ? Thanks

Shapes do not match: <Argument id=139987557997272> and <Product id=139987565580864>.
---------------------------------------------------------------------------
UFLException Traceback (most recent call last)
<ipython-input-1-ed452f40a979> in <module>()
18
19
---> 20 C = project(I4,V)

/usr/lib/python3/dist-packages/dolfin/fem/projection.py in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
136 Pv = TrialFunction(V)
137 a = ufl.inner(w, Pv)*dx
--> 138 L = ufl.inner(w, v)*dx
139
140 # Assemble linear system

/usr/local/lib/python3.5/dist-packages/ufl/operators.py in inner(a, b)
143 if a.ufl_shape == () and b.ufl_shape == ():
144 return a*b
--> 145 return Inner(a, b)
146
147

/usr/local/lib/python3.5/dist-packages/ufl/tensoralgebra.py in __new__(cls, a, b)
163 ash, bsh = a.ufl_shape, b.ufl_shape
164 if ash != bsh:
--> 165 error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
166
167 # Simplification

/usr/local/lib/python3.5/dist-packages/ufl/log.py in error(self, *message)
170 "Write error message and raise an exception."
171 self._log.error(*message)
--> 172 raise self._exception_type(self._format_raw(*message))
173
174 def begin(self, *message):

UFLException: Shapes do not match: <Argument id=139987557997272> and <Product id=139987565580864>.​

​
written 4 months ago by Christian
2
Do I4 = outer(Identity(3), Identity(3)) or for index notation I4 = as_tensor(delta(i, j) * delta(k, l), (i, j, k, l))
written 4 months ago by Michal Habera