Modeling a conservative flux source term with div

Question

I'd like to model a part of the comsol coefficient form for PDEs. First equation in eqn. 16-2. Specifically, I am working on the $\nabla \cdot \gamma$ term. To understand this I have the following script:

from dolfin import *                                                               

mesh = UnitIntervalMesh(3)                                                         
w = FunctionSpace(mesh, 'P', 1)                                                    
#u = TrialFunction(w)                                                              
v = TestFunction(w)                                                                

def gamma(x): return x                                                             

x = MeshCoordinates(mesh)                                                          
f = Constant(1.0, cell = mesh.ufl_cell()) * gamma(x)                               
R = assemble(inner(f, grad(v))*dx(domain = mesh))                                  

print(R.array())                                                                   

Running this returns

[ 0.83333333 -0.33333333 -0.33333333 -0.16666667]

What I do not understand is that if I change the def of gamma to, say,

def gamma(x) return x**2

I get an error message like:

Shapes do not match: <Product id=139744267574944> and <Grad id=139744267598984>.
Traceback (most recent call last):
  File "test_1D_2_grad.py", line 14, in <module>
    R = assemble(inner(f, grad(v))*dx(domain = mesh))

which I can not make sense of since both x and x**2 are scalar. What am I missing?

Edit

I found that this also works:

R = assemble(inner(gamma(x), grad(v))*dx(domain = mesh))

which makes me believe that x gets the entire vector of mesh nodes. How would I operate on that vector to get, say, a x**2? This does not seem to work:

def gamma(x): return np.array(x) * np.array(x)

Comments

Hi, see by x.ufl_shape that x here is a vector of length one and so is grad(f). On the other hand x**2 is a scalar, cf. (x**2).ufl_shape. The get the form with scalars only consider

gamma(x[0])*v.dx(0)*dx(domain=mesh)

---

Thanks for looking at this and copying your expression for gamma in assemble does seem to work. I, however, have trouble understanding your comment. What does x contain? I was under the impression it has all coordinates from the mesh. Why is it only length 1 then?

If you make your comment into an answer I could accept it.

---

The following works with FEniCS 2016.2.0

from dolfin import *                                                               

mesh = UnitIntervalMesh(3)                                                         
w = FunctionSpace(mesh, 'P', 1)                                                    
u = TrialFunction(w)                                                              
v = TestFunction(w)                                                                

def gamma(x): return x**2

x = MeshCoordinates(mesh)                                                          
f = Constant(1.0, cell = mesh.ufl_cell()) * gamma(x[0])                               
R = assemble(f*v.dx(0)*dx(domain = mesh))                                  

print(R.array())

MeshCoordinates(mesh) is an UFL object representing the concept of position vector for the purpose of defining UFL forms. In your case of interval mesh, the length of vector is one (only one coordinate per point). It is not coordinates of mesh vertices.

---

OK, I start understand. x is a representation of a coordinate from mesh. And I can also

def gamma(x): return x[0]**2

Now, how could I change the def to return a length one vector? That way, I think, I could use it with grad(v) and make the assemble dimension independent.

Once more thanks a lot.

---

What should gamma be in 2d and 3d? Also, if v is scalar, the length of grad(v) in $R^d$ is d.

---

I'd say gamma is a vector of length d then - but I'd need to think a bit more about this. Thanks once more.