Convert mesh coordinates to UFL type


225
views
0
6 months ago by
J.H.59  
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()
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
6 months ago by
Miguel  


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
6 months ago by
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
6 months ago by
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\).
Please login to add an answer/comment or follow this question.

Similar posts:
Search »