project(grad(u), V)

Question

Hi folks,

I am trying to project the gradient of my solution to a function space, say V. project(grad(u), V) used to work well and the output was a vector. But after updating FEniCS to its latest version, I get the following error:

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Shape mismatch.
Traceback (most recent call last):
File "test_cj_simple.py", line 30, in
g_cj = project(grad(cj), V)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/projection.py", line 100, in project
L = ufl.inner(w, v)ufl.dx()
File "/usr/lib/python2.7/dist-packages/ufl/operators.py", line 127, in inner
return Inner(a, b)
File "/usr/lib/python2.7/dist-packages/ufl/tensoralgebra.py", line 184, in new
ufl_assert(ash == bsh, "Shape mismatch.")
File "/usr/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert
if not condition: error(message)
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 154, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shape mismatch.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

This error is produced running the following code:

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

from dolfin import *
import numpy

n = 100
mesh = IntervalMesh(n,0,1)
V = FunctionSpace(mesh, 'CG', 1)
mesh_coor = mesh.coordinates()

cj = TrialFunction(V)
v_cj = TestFunction(V)

x0_cj = mesh_coor[50][0]
f_cj = PointSource(V, Point(x0_cj))
a_cj = cjv_cjdx + inner(grad(cj), grad(v_cj))*dx
A_cj = assemble(a_cj)
B_cj = dolfin.Vector(A_cj.size(0))
B_cj.array()[:] = 0
f_cj.apply(B_cj)
cj = Function(V)
solve(A_cj, cj.vector(), B_cj)

g_cj = project(grad(cj), V)

plot(cj)
interactive()

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Any help is deeply appreciated.

Thanks,
Mo

Answer 1

Hi,

There is a problem with the grad operator in 1D. Try this instead:

g_cj = project(cj.dx(0), V)

Best regards

Mikael

Comments

This is not a problem but consistent behaviour. You always obtain rank 1 object (vector) by taking grad of scalar. Use rather div to get scalar from scalar.

---

This has been discusseed several times before. See for example here or here. Not sure whether it was recently decided it should not work anymore (your consistent bahaviour), or whether this is a new bug? It used to work.

---

Consider that code

g_cj = project(grad(cj), V)

cannot work in 2D, 3D either. On the other hand

g_cj = project(div(cj), V)

will work in arbitrary dimension. That's the consistency. You should ask Martin whether it has been decided. I would say that it is a consequence of some clean-up in UFL. I think there was an issue with inconsistent behaviour of V.cell().x - you had to use V.cell().x in 1D but V.cell().x[0] in 2D, ... to get $x$ coordinate. This was a consequence of treating 1D vector as a scalar. To keep code portable between 1D and higher-D you need to precisely distinguish scalars (shape ()) and vectors (shape (1,) in 1D).

---

Mikael and Jan, thank you both for your detailed answers.

---

Dear Jan,

But When I run the code with g_cj = project(div(cj), V), I get the following error:

Invalid number of indices (1) for tensor expression of rank 0:
Coefficient(FiniteElement('Lagrange', Domain(Cell('interval', 1), 'interval_multiverse', 1, 1), 1, None), 0)

Traceback (most recent call last):
File "test_cj_simple.py", line 24, in
g_cj = project(div(cj), V)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/projection.py", line 104, in project
form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 281, in assemble_system
subdomains, form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 56, in init
common_cell)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 66, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 154, in jit
return jit_compile(form, parameters=p, common_cell=common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 77, in jit
return jit_form(ufl_object, parameters, common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 159, in jit_form
element_mapping=element_mapping)
File "/usr/lib/python2.7/dist-packages/ufl/form.py", line 123, in compute_form_data
element_mapping=element_mapping)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/preprocess.py", line 115, in preprocess
form = expand_derivatives(form, common_cell.geometric_dimension()) # EXPR
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/ad.py", line 107, in expand_derivatives
return transform_integrands(form, _expand_derivatives)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 252, in transform_integrands
integrand = transform(integrand)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/ad.py", line 88, in _expand_derivatives
expression = expand_compounds(expression, dim)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/expand_compounds.py", line 277, in expand_compounds1
return apply_transformer(e, CompoundExpander(dim))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 277, in apply_transformer
return transform_integrands(e, _transform, domain_type)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 268, in transform_integrands
return transform(expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 276, in _transform
return transformer.visit(expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/expand_compounds.py", line 174, in div
return a[...,i].dx(i)
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 407, in _getitem
a = Indexed(self, indices)
File "/usr/lib/python2.7/dist-packages/ufl/indexed.py", line 47, in __init__
% (len(indices), expression.rank(), expression))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 154, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Invalid number of indices (1) for tensor expression of rank 0:
Coefficient(FiniteElement('Lagrange', Domain(Cell('interval', 1), 'interval_multiverse', 1, 1), 1, None), 0)

---

Ah, I was wrong, of course. You can get scalar by applying div on vector. The codes

V = FunctionSpace(mesh, 'CG', 1)
v = Function(V)
u = project(tr(grad(u)), V)

W = VectorFunctionSpace(mesh, 'CG', 1)
w = Function(W)
v = project(div(w), V)

should work in arbitrary dimension. From your 1D example it is not obvious how to generalize it to higher dimension. If you don't need doing it then nevermind.

---

Jan can i know what means tr represents in code above you said? Can you explain me. I am also using this code in my program.

---

tr should mean trace - so the code above is probably nonsense.