Hi everybody,
Does someone know how to print values of the tensor obtained as follow :
T = symgrad(u)
Where u is a Function structure
Many thanks
You can project UFL objects like T onto finite element FunctionSpaces for visualization, e.g.,
from dolfin import *
N = 16
mesh = UnitSquareMesh(N,N)
x = SpatialCoordinate(mesh)
V = VectorFunctionSpace(mesh,"CG",1)
# Obtain nonzero vector-valued Function for testing:
u = project(as_vector(2*[sin(pi*x[0])*sin(pi*x[1]),]),V)
# Define symmetric gradient of it in UFL:
T = sym(grad(u))
# Project to a tensor-valued FunctionSpace
# and visualize in ParaView:
VT = TensorFunctionSpace(mesh,"DG",0)
T_proj = project(T,VT)
File("T.pvd") << T_proj
# Alternatively, project an individual
# component to a scalar-valued FunctionSpace:
V0 = FunctionSpace(mesh,"DG",0)
T00_proj = project(T[0,0],V0)
File("Txx.pvd") << T00_proj
The Paraview visualization of tensor-valued functions appears to assume that tensors are 3x3, and gives a list of components indexed 0–8, with 4 of them nonzero, so it may take some reverse engineering to figure out which component is which.
You can also use dolfinx.Expression to interpolate into an appropriate function space, as described in: Implementation — FEniCSx tutorial
maybe I exposed badly my issue. I would like to see inside the object tensor Sym{} the number which correspond to the entries of the matrix. I can visualize it into paraview with all its components. What I can’t is to exploit matrix substantially.
Once you have the field as a Function object (following my post in old FEniCS or Jorgen’s example for FEniCSx), you can evaluate that Function at points. Continuing from my example, you could get the components at the point (0.3,0.3) with
print(T_proj([0.3,0.3]))
I don’t know why I cant print out the matrix. I did the things you recommended
# Define symmetric part of gradient di u_y
def symgrad(f):
return sym(grad(f))
T = TensorFunctionSpace(mesh, 'CG', 1)
t = project(symgrad(f1), T)
# Save for postprocessing in MATLAB
sc.io.savemat('ug_1.mat', {'t': t.vector().get_local()}) # Save coefficient of t
Then I print as you
print(t([0.3,0.3]))
but I get the following error:
TypeError: expected the geometry argument to be of length 3
It looks like you’re probably using a 3D mesh instead of a 2D one, so it’s expecting three coordinates for the point at which to evaluate t.
Yes, but if I want to show the matrix entries (which is a 3x3) what sense has the number 0.3 ? I have three possibilities: 0,1,2 which corresponds to x,y,z. Is there any way as in MATLAB for example to visualize the matrix ?
t is not a single 3x3 matrix, but rather a matrix-valued function. It will, in general, have different matrix entries at different points in the domain. To get a 3x3 matrix of actual floating-point numbers, you need to pick a specific point in space at which to evaluate t. I chose the point (0.3,0.3) arbitrarily in my example; there’s no special significance to those coordinates, aside from being inside of the mesh that I used in my example.
Ok. I understand, and if I chose the point (x0,y0,z0) how can I extract the tensor entries in this point?
Lastly, if my domain correspond to a sphere immersed in cylindrical channel, obtained from subtracting the sphere to the channel in the cad model, can I assume that the element dx of volume is correspondent to all the region except the sphere? or I need to specify this is some way? In GMSH i got only one volume label, 1.
if I chose the point (x0,y0,z0) how can I extract the tensor entries in this point?
If x0, y0, and z0 are Python floats, you can pass them as a list to t, e.g.,
print(t([x0,y0,z0]))
can I assume that the element dx of volume is correspondent to all the region except the sphere?
By default, dx, adds up integrals over all elements in the mesh. If the interior of the sphere is not meshed, then it will not be integrated over.
1 Like