DOLFIN

Gradient operator in one dimension returns scalar

Asked by Charl

Hello,

I have made a toy 1D problem to get a handle on SUPG for electro-diffusion of particles. The scalar potential is calculated in the usual manner (cf. FEniCS book), and looks good on a plot. Now, curiously, when I try to derive the E-field from this, using:

mesh_funcspace_vector = VectorFunctionSpace(mesh, 'Lagrange', 1)
efield = project(-grad(potential), mesh_funcspace_vector)

...I receive a UFL shape mismatch error. The projection works if I project to a scalar function space--but this baffles me, as obviously the gradient of the scalar potential is a (1D) vector E-field.

(I also tried manual projection, where I get the same error; see below.)

Thanks in advance,
Charl

---------

In [21]: potential
Out[21]: Coefficient(FiniteElement('Lagrange', Cell('interval', Space(1)), 1, None), 2)

In [22]: inner(grad(potential), TestFunction(mesh_funcspace_vector))*dx
Shape mismatch.

In [23]: inner(grad(potential), TestFunction(mesh_funcspace_scalar))*dx
Out[23]: Form([Integral(Product(Argument(FiniteElement('Lagrange', Cell('interval', Space(1)), 1, None), -2), SpatialDerivative(Coefficient(FiniteElement('Lagrange', Cell('interval', Space(1)), 1, None), 2), MultiIndex((FixedIndex(0),), {}))), Measure('cell', 0, None))])

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Charl
Solved:
Last query:
Last reply:
Revision history for this message
Martin Sandve Alnæs (martinal) said : 		#1

Your observation is correct, UFL handles grad (and facet normal) in 1D as scalars.

I don't recall the pros and cons I had before me when
making this decision, but I remember I had to clean up
some inconsistencies around 1D vector vs scalar last year.

Is there a particular reason you want it to be a vector field or are you just curious?

Revision history for this message
Charl (78luphr0rnk2nuqimstywepozxn9kl19tqh0tx66b5dki1xxsh5mkz9gl21a5rlwfnr8jn6ln0m3jxne2k9x1ohg85w3-launchpad) said : 		#2

Dear Martin,

Thanks for clearing that up. The reason for wishing to use a vector field is to make this test code more portable to higher dimensions. I tried to hack around it by copying the values to a vector function:

efield_ = project(-grad(potential), mesh_funcspace_scalar)
efield = Function(mesh_funcspace_vector, efield_)

But when I try to use this function later:

assemble(trialfunc*inner(efield, nabla_grad(testfunc))*dx)

... it works using the scalar, but not the vector! So I suppose this design philosphy (scalars instead of vectors in the 1D case) is quite pervasive and I should just stick to the scalar functions for my dummy problem.

Revision history for this message
Martin Sandve Alnæs (martinal) said : 		#3

You could introduce a helper function to map scalar expressions to 1D vector expressions:

def s2v(e): # Not tested!
    return e if e.shape() else as_vector((e,))

and then do
efield_ = project(s2v(-grad(potential)), mesh_funcspace_vector)

As you've discovered, this choice has quite a few side effects so it's not so easy to change.