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

2 Answers


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()
        ax = fig.add_subplot(111, aspect='equal')

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

        the_divider = make_axes_locatable(ax)
        color_axis = the_divider.append_axes("right", size="5%", pad=0.1)

        # 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
See https://fenicsproject.org/qa/1760/output-a-surface-of-a-function-defined-over-a-3d-mesh/
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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »