### How to visualize surface values of 3d function?

206
views
0
4 months ago by
Greetings.
I have a function and it's visualization:
from dolfin import *
from dolfin.cpp.mesh import UnitCubeMesh
import matplotlib.pyplot as plt
omega = UnitCubeMesh(5,5,5)
finite_element = FiniteElement("CG", omega.ufl_cell(), 1)
u = interpolate(Expression("x[0]+x[1]+x[2]", degree=3), FunctionSpace(omega, finite_element))
plot(u)
plt.show()
​
But how can I get a plot of u function values on bottom surface of omega?
For square I'm using this function:
def get_2d_boundaries(v):
left = AutoSubDomain(lambda x, on_bnd: near(x[0], 0) and on_bnd)
top = AutoSubDomain(lambda x, on_bnd: near(x[1], 1) and on_bnd)
right = AutoSubDomain(lambda x, b: near(x[0],1) and b)
bottom = AutoSubDomain(lambda x, b: near(x[1],0) and b)
V = v.ufl_function_space()
bc0 = DirichletBC(V, 1, left)
bc1 = DirichletBC(V, 2, top)
bc2 = DirichletBC(V, 3, right)
bc3 = DirichletBC(V, 4, bottom)
u = Function(V)
bc0.apply(u.vector())
bc1.apply(u.vector())
bc2.apply(u.vector())
bc3.apply(u.vector())
# v = interpolate(Expression("x[0]+x[1]", degree=2), V)
print("Values on left = ", *v.vector()[u.vector() == 1][::-1])
print("Values on top = ", *v.vector()[u.vector() == 2])
print("Values on right = ", *v.vector()[u.vector() == 3])
print("Values on bottom = ", *v.vector()[u.vector() == 4][::-1])
return
​

It's not perfect solution though, but it works. But in case of 3d function array of values is not enough to draw the plot.

1. How can I get six plots of function u on cube boundaries?
2. Is it possible to combine it in one image?

Thank you in advance for help.

Community: FEniCS Project

3
4 months ago by

Actual solution is much simpler than I thought.
We can use the fact that function(Point(x,y,z)) returns us a value in this point.

What I currently use:

def get_3d_boundaries(v, name='solution'):
top = lambda x, y: v(Point(x,y,1))
bottom = lambda x, y: v(Point(x,y,0))
def front_surfaces(x, y):
if x <= 1:
return v(Point(0,x,y))
if x <= 2:
return v(Point(x-1,1,y))
if x <= 3:
return v(Point(1,3-x,y))
if x <= 4:
return v(Point(4-x,0,y))

for action in [(top, 'top'), (bottom, 'bottom'), (front_surfaces, 'belt')]:
X, Y = meshgrid(linspace(0,1, num=100)*(4 if action[1] =='belt' else 1),
linspace(0,1, num=100))
top_vals = array(list(map(action[0],
X.reshape(100**2), Y.reshape(100**2)))).reshape(100,100)
fig = plt.figure()

SC = ax.pcolor(X, Y, top_vals, cmap=cm.RdBu, vmin=top_vals.min(), vmax=top_vals.max())

the_divider = make_axes_locatable(ax)

# Colorbar
cbar = plt.colorbar(SC, cax=color_axis)
plt.savefig('{}_{}.png'.format(name, action[1]))
plt.gcf().clear()
return


Code is not perfect but for sure demonstrates main idea :)

1
4 months ago by
I'm still confused with it.
As far as I understand I have to create a boundary mesh and define a function over it:
bmesh = BoundaryMesh(omega, "exterior")
Vb = FunctionSpace(omega, "CG", 1)
ub = Function(Vb)
ub.interpolate(u)
​

But then I have to create a submesh that correlates for some face of the cube. How to implement it in the code I have no idea.

written 4 months ago by Pavel Mesenev
For the UnitCubeMesh, this seems to work:
from fenics import *
import matplotlib.pyplot as plt

mesh = UnitCubeMesh(10, 10, 10)
V = FunctionSpace(mesh, 'P', 1)
u = interpolate(Expression('(x[0]+x[1])*(1-x[2])', degree=1), V)
plt.figure()
plot(u)

mesh_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
bottom = CompiledSubDomain('on_boundary && near(x[2], 0.0)')
bottom.mark(mesh_markers, 42)

bmesh   = BoundaryMesh(mesh, 'exterior')
cellmap = bmesh.entity_map(2)
bmesh_markers = MeshFunction("size_t", bmesh, bmesh.topology().dim(), 0)

for c in cells(bmesh):
if mesh_markers[cellmap[c.index()]] == 42:
bmesh_markers[c] = 1

smesh = SubMesh(bmesh, bmesh_markers, 1)
Vb = FunctionSpace(smesh, 'P', 1)
ub = interpolate(u, Vb)

plt.figure()
plot(ub)

plt.show()​
written 4 months ago by Adam Janecka
This solution works ok for bottom and top surfaces of the cube: https://ibb.co/focb9n
But, it's not working for the front face: https://ibb.co/h9aSh7
Looks like some additional mesh transformation required. In both cases I just want to get flat colored image with the color bar.
written 4 months ago by Pavel Mesenev