### Convert mesh coordinates to UFL type

225

views

0

Hi all,

I want to calculate something like

\(\begin{align*}

\int_{S_p} \mathbf{r}_0 \times \mathbf{u}\, dS

\end{align*}\)

on a boundary \(S_p\) of the loaded mesh where \(\mathbf{r}_0 \) is the origin vector and \(\mathbf{u}\) the velocity which is the solution of the stokes equation solved before.

My first approach was to get the coordinates of the mesh with mesh.coordinates()

I want to calculate something like

\(\begin{align*}

\int_{S_p} \mathbf{r}_0 \times \mathbf{u}\, dS

\end{align*}\)

on a boundary \(S_p\) of the loaded mesh where \(\mathbf{r}_0 \) is the origin vector and \(\mathbf{u}\) the velocity which is the solution of the stokes equation solved before.

My first approach was to get the coordinates of the mesh with mesh.coordinates()

```
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
#Boundary conditions + variational problem
.
.
.
#Compute solution
U = Function(W)
solver.solve(U.vector(), bb)
#Get sub-functions
(u, p) = U.split()
r = mesh.coordinates()
c = cross(r, u)
c = assemble(c * dS)
```

but this results in an error:```
Traceback (most recent call last):
File "test.py", line 100, in <module>
c = cross(r, u)
File "/home/anaconda2/envs/fenics_test/lib/python3.6/site-packages/ufl/operators.py", line 200, in cross
a = as_ufl(a)
File "/home/anaconda2/envs/fenics_test/lib/python3.6/site-packages/ufl/constantvalue.py", line 416, in as_ufl
" to any UFL type." % str(expression))
ufl.log.UFLValueError: Invalid type conversion: [[-10. 10. -10. ]
[ 10. 10. -10. ]
[-10. -10. 10. ]
...,
[ 5.09017804 9.01670323 7.85035017]
[ 4.51327702 8.93986382 6.95476418]
[ -4.11376381 5.89215415 -2.70712164]] can not be converted to any UFL type.
```

So it seems that mesh.coordinates has a numpy array as an ouput, but to use cross with \(\mathbf{u}\) it needs to be in an UFL type.

Is there a way to load the coordinates of the mesh in such a format or to convert them?

Thanks in advance.

Community: FEniCS Project

### 3 Answers

3

FEniCS cannot evaluate integrals of rank > 0. Integrate each component separately.

```
from dolfin import *
mesh = UnitCubeMesh(10,10,10)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((1.0, 1.0, 0.0))
a = inner(grad(u), grad(v))*dx + inner(u,v)*dx
L = inner(f, v)*dx
u_sol = Function(V)
solve(a==L, u_sol)
r = SpatialCoordinate(mesh)
def my_cross(v,w):
return (v[1]*w[2]-v[2]*w[1],v[2]*w[0]-v[0]*w[2],v[0]*w[1]-v[1]*w[0])
cross_product_1 = my_cross(r,u_sol)[0]*ds
cross_product_2 = my_cross(r,u_sol)[1]*ds
cross_product_3 = my_cross(r,u_sol)[2]*ds
print(assemble(cross_product_1))
print(assemble(cross_product_2))
print(assemble(cross_product_3))
```

Thank you for your fast reply. That was exactly I was looking for

written
6 months ago by
J.H.59

1

I think this should do the trick:

```
r = Expression(("x[0]","x[1]"),degree=1,domain=mesh)
```

EDIT: The SpatialCoordinate() function is a cleaner solution; see Miguel's answer.

1

To evaluate vector-valued integrals, use the answer to this question. Is \(\mathbf{r}_0\) constant throughout \(S_p\)? If so, you can calculate

\[\bar{\mathbf{u}} = \int_{S_p} \mathbf{u}\;dS\]

separately, and then use numpy to compute \(\mathbf{r}_0\times\bar{\mathbf{u}}\). Otherwise you can define a Function object such as the position in the question referenced above, and use that to integrate the cross product throughout \(S_p\).

\[\bar{\mathbf{u}} = \int_{S_p} \mathbf{u}\;dS\]

separately, and then use numpy to compute \(\mathbf{r}_0\times\bar{\mathbf{u}}\). Otherwise you can define a Function object such as the position in the question referenced above, and use that to integrate the cross product throughout \(S_p\).

Please login to add an answer/comment or follow this question.