Problem with creating vector with scalar function components

Hello!

In my problem I solved the Poisson equation for a unidirectional flow V=(0, 0, u (x, y)). I solved for u = u(x, y) successfully and implemented the following code to form the vector V and the stress tensor S in fluid dynamics:

``````Vx=0.0
Vy=0.0
Vz=u
V=[Vx, Vy, Vz]
``````

`S = (mu/2)*(grad(V)+(grad(V)).T)`

This gives the following error:

``````File "shear.py", line 76, in <module>
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 369, in grad
f = as_ufl(f)
File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: [0.0, 0.0, Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 20)] can not be converted to any UFL type.
``````

I am not sure what the error means. Any suggestions would be very helpful.

Use `V = as_vector((Vx, Vy, Vz))`

1 Like

Hi dokken! That worked. Thanks.

However, when I am printing the vector V, it gives me `[0, 0, f_20]`. What does f_20 mean?

Also, when I am executing the command

`S = (mu/2)*(grad(V) + ((grad(V)).T))`

itâs giving me the following error message:

``````Solving linear variational problem.
Can't add expressions with different shapes.
Traceback (most recent call last):
File "shear.py", line 79, in <module>
File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 201, in _add
return Sum(self, o)
File "/usr/lib/python3/dist-packages/ufl/algebra.py", line 41, in __new__
error("Can't add expressions with different shapes.")
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different shapes.
``````

It is the unique name dolfin has given to your function u. You can change this by:
`u.rename("replace_with_new_name", "")`

Shouldnât this be either:

• `S= (mu/2) * (grad(V) + grad(V).T)`
Or equivalently
• `S=mu*sym(grad(V))`

For future problems please make a minimal reproducable example, following the guidelines in:

dokken, thanks for your response. Your suggested code partially worked but it gave rise to another error. Below is the description of the problem follwed by the code and the error:

Problem: Calculating the velocity profile of blood flow through an artery that follows the Poisson equation. V is the blood velocity (unidirectional) that gives rise to shear stress on the vessel wall. I need to calculate the velocity profile (V) and the shear stress (S).

Python Script:

``````from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh

import meshio
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
#print(assemble(1*ds_custom(1)) ,assemble(1*ds_custom(2)))

# Normal vector to the intima
n = FacetNormal(mesh)

# Create mesh and define function space

V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
G , mu = 10, 0.000004
u_D=Constant(0.0)

bc = DirichletBC(V, u_D, mf, 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Ux, Uy, u))
#print('shear stress  =', S)
plot(S)
plot(mesh)
plt.show()

# Save solution to file in PVD format for Paraview
file = File("lumenshear.pvd");
file << u;
``````

Error message:

``````Solving linear variational problem.
Traceback (most recent call last):
File "shear.py", line 79, in <module>
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 369, in grad
f = as_ufl(f)
File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 437, in as_ufl
raise UFLValueError("Invalid type conversion: %s can not be converted"
ufl.log.UFLValueError: Invalid type conversion: FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)) can not be converted to any UFL type.``````

That is because your V here should be U.

Also note that the problem above is not reproducible, as you rely on a non built in mesh.

I tried that and changed V to U. And here is the error message:

``````Traceback (most recent call last):
File "shear.py", line 79, in <module>
File "/usr/lib/python3/dist-packages/ufl/exproperators.py", line 201, in _add
return Sum(self, o)
File "/usr/lib/python3/dist-packages/ufl/algebra.py", line 41, in __new__
error("Can't add expressions with different shapes.")
File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different shapes.``````

If you are using a 2D mesh, you would have to not prune the z-coordinates if you are trying to combine `grad(U)` and `grad(U).T`. If you are using a mesh with pruned z-coordinate, grad(U) is a 3x2 matrix, while grad(U).T is a 2x3 matrix. grad is only interpreted as a 3D quantity if you have z-coordinates in your mesh.

Also note that you need to project `S` into an appropriate function space to plot it

Hi dokken! Thanks for the feedback. I understand your suggestion but I am trying to do it a bit differently. Here is what I am trying to do:

Step 1: Suppose a liquid is flowing through a tube. The Z-axis runs through the middle of the tube along its length. At a fixed point on the Z-axis, I cut out a cross-section of the tube. Now I get a circle and I figure out the velocity profile in that circle (which in my code is the V(x, y)).

Step 2: Next, I take that output V(x, y) and create the velocity vector V = (0, 0, V(x, y)) (here I donât need the geometry anymore). Then the shear stress is given by `S = (mu/2) * (grad(V) + grad(V).T)`

According to me, Step 2 doesnât need the geometry anymore. Itâs simply computing a vector out of the output and doing matrix operations to give the shear stress.

Please correct me if I am wrong.

If I donât prune the Z-coordinate by the following code:

``````import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
#if prune_z:
#out_mesh.prune_z_0()
#return out_mesh
``````

itâs giving me the following error:

``````Traceback (most recent call last):
File "shear.py", line 33, in <module>
meshio.write("mesh.xdmf", triangle_mesh)
File "/home/WIN/avishmj/.local/lib/python3.8/site-packages/meshio/_helpers.py", line 132, in write
for key, value in mesh.cells:
AttributeError: 'NoneType' object has no attribute 'cells'``````

You still need to return the out_mesh (the third line you have commented out).

That worked! Thanks dokken!

In your previous post, you said " Also note that you need to project `S` into an appropriate function space to plot it". I didnât get that. Why is it necessary?
Do you mean I have to plug in the command `S=FunctionSpace(V)` before `S = (mu/2) * (grad(V) + grad(V).T)`

I am facing another problem. `S = (mu/2) * (grad(V) + grad(V).T)` gives me the shear on the entire domain. In order to figure out the shear on the boundary I implemented the command
`S = (mu/2) * dot((grad(V) + grad(V).T), n)*ds_custom(1)` where `n=FacetNormal(mesh)` and 1 is the physical tag of the boundary.

This gives me the error:

``````Traceback (most recent call last):
``````

File âshear.pyâ, line 79, in
TypeError: dot() missing 1 required positional argument: âbâ

How to overcome this and what does it mean?

For the projection you need to do something like:

``````StressSpace = TensorFunctionSpace(mesh, "DG", 0)
S_func = project(S, StressSpace)
``````

Please note that to assemble the normal stress you need to reduce it to considering a single component at the time. As long as you are not supplying a minimal example that reproduces the code, rather than modified snippets from the original post.

Ok, so I am providing the .geo script and the full python script below.

gmsh code:

``````//+
//+
Point(3) = {-1.6, 0.3, 0, 1.0};
//+
Point(4) = {-1.4, -0.7, 0, 1.0};
//+
Point(5) = {-1.1, -1.3, 0, 1.0};
//+
Point(6) = {-0.3, -1.7, 0, 1.0};
//+
Point(7) = {0.8, -2.1, 0, 1.0};
//+
Point(8) = {1.7, -1.3, 0, 1.0};
//+
Point(9) = {1.7, -0.4, 0, 1.0};
//+
Point(10) = {1.9, 0.5, 0, 1.0};
//+
Point(11) = {1.6, 1.2, 0, 1.0};
//+
Point(12) = {0.9, 1.5, 0, 1.0};
//+
Point(13) = {0, 1.5, 0, 1.0};
//+
Point(14) = {-0.6, 0.9, 0, 1.0};
//+
Spline(3) = {3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 3};
//+
Physical Curve("C1", 1) = {3};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {1};
//+
Hide "*";
//+
Show {
Point{3}; Curve{3};
}
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3};
}
//+
Show "*";
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2};
}
//+
Curve Loop(5) = {3};
//+
Plane Surface(3) = {5};
//+
Physical Surface("lumen", 6) = {3};//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{3};
}
``````

python script:

``````from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}

import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
#if prune_z:
#out_mesh.prune_z_0()
return out_mesh

import meshio
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
#print(assemble(1*ds_custom(1)) ,assemble(1*ds_custom(2)))

# Normal vector to the intima
n = FacetNormal(mesh)

# Create mesh and define function space

V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
G , mu = 10, 0.000004
u_D=Constant(0.0)

bc = DirichletBC(V, u_D, mf, 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Ux, Uy, u))
StressSpace = TensorFunctionSpace(mesh, "DG", 0)
S_func = project(S, StressSpace)
#print('shear stress  =', S)
plot(S)
plot(mesh)
plt.show()

# Save solution to file in PVD format for Paraview
file = File("lumenshear.pvd");
file << u;
``````

The lumen.msh is the msh file of the gmsh script provided.

How do you suppose these lines should work? S should be defined before projecting it. Please remove all unnecessary code that is not needed to reproduce your error message.

Hi dokken! I made the suggested changes as follows:

``````StressSpace = TensorFunctionSpace(mesh, "DG", 0)
S_func = project(S, StressSpace)
``````

But the same error persists:

`Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3,) and free indices with labels ().`

In your previous post, you mentioned that the stress components must be evaluated element-wise. But on trying the split command, itâs giving me the error `âGradâ object has no attribute âsplitâ

I am confused as what to do now.
`

Also tried

``````StressSpace = TensorFunctionSpace(mesh, "DG", 0)
S_func = project(Sz, StressSpace)
``````

which gave the error:

`'ListTensor' object has no attribute 'sub'`

Consider the following:

``````U = as_vector((Ux, Uy, u))

print(assemble(Sn[0]*ds), assemble(Sn[1]*ds), assemble(Sn[2]*ds))
``````

This assembles the three components of your dot product (over the whole boundary, as I have not used your mesh with mesh tags).

1 Like

Thanks dokken! This works.

Just to clarify, Sn is a vector of the form Sn=(Sn[0], Sn[1], Sn[2]) on the boundary of my region. Is that correct? The code `assemble(Sn[2]*ds)` calculates Sn[2]*ds at each nodal point on the boundary and adds them up to give the total value of Sn[2] on the boundary. Am I correct?
Now if I want to calculate the third entry only, i.e., Sn[2], at each point of the boundary, and want to visualize it in paraview, I need to plug the commands

``````file = File("lumenshear.pvd");
file << Sn[2];
``````

Is my understanding correct?

This command assembles (integrates) the function Sn[2] over the boundary `ds`.
If you want to visualize `Sn[2]` on your mesh, you need to project `Sn[2]` into a suitable function space (usually a discontinuous space of one lower degree than the displacement). As you are using the normal vector, which is only defined on facets, you need to do a custom projection of `Sn[2]`, see for instance How to plot normal unit vector of faces in a 2D mesh? - #2 by dokken
where you should change `l` to be `l = inner(Sn[2], v)*ds`.
As the stress is a discontinuous function if you follow my suggestion of appropriate function space, you should use `XDMFFile.write_checkpoint` so that you visualize the actual function and not a CG-` interpolation. See for instance: Read mesh from XDMF file (write_checkpoint) - #2 by WeiOng

1 Like