2D slice from 3D solution

1

Hello everyone,
I need your help. Indeed I solved a transient problem of soft tissue mechanics on a 3D slab. However, I would like to have the solution on an xy plane of the domain for z=0 and in the middle (z=0.5).
I have used some solutions that I have found on the internet but, I get no result.
I have mainly tested the solutions in the links:

https://fenicsproject.org/qa/3257/building-2d-function-g-x-y-from-a-3d-function-f-x-y-z-constant/

Thank you in advance.

# Create a 2D unit square mesh from the boundary of mesh3D
mesh2 = BoundaryMesh(mesh3D, "exterior")
side = MeshFunction("size_t", mesh2, 0)
xyplane = AutoSubDomain(lambda x: x[2] < 1e-8)
xyplane.mark(side, 1)
mesh3 = SubMesh(mesh2, side, 1)
x = mesh3.coordinates()
x[:, 2] += 0.5
VH2 = FunctionSpace(mesh3, "DG", 1)
# solution obtained from mechanical problem

u,p=up.split()

# Interpolation for 2D slice from 3D

u1 = interpolate(u, VH2)
2

You can try defining an expression that you can interpolate to get rid of your third axis.

Something like this:

#Dummy problem

from fenics import *
import matplotlib.pyplot as plt

# Create mesh and define function space
mesh = UnitCubeMesh(8, 8, 8)
V = FunctionSpace(mesh, 'CG', 1)

# Define boundary condition
u_D = Constant(0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(100)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

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

# Plot solution and mesh
plt.colorbar(plot(u))
plt.title('u')

full_sol

Define an expression to evaluate u at some value of z:

class u_slice(UserExpression):
    def eval(self, value, x):

        value[0] = u(x[0],x[1],0.5)

    def value_shape(self):
        return () 

mesh_2D = UnitSquareMesh(8, 8)
V_2D = FunctionSpace(mesh_2D, 'CG', 1)
u_sl = interpolate(u_slice(), V_2D)
plt.colorbar(plot(u_sl))
plt.title('u in the middle plane')

plane_sol

3

Hi
Thank you very much for your reply.
I managed to do it differently. And it and your code work perfectly.
Thank you.

1 Like
4

Hi!

Thanks for this!

A few questions:

  1. Plotting 3D functions with fenics doesn’t seem possible anymore
    Your code yields:
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[2], line 30
     27 solve(a == L, u, bc)
     29 # Plot solution and mesh
---> 30 plt.colorbar(plot(u))
     31 plt.title('u')

File /opt/conda/lib/python3.10/site-packages/dolfin/common/plotting.py:438, in plot(object, *args, **kwargs)
    436 # Plot
    437 if backend == "matplotlib":
--> 438     return _plot_matplotlib(object, mesh, kwargs)
    439 elif backend == "x3dom":
    440     return _plot_x3dom(object, kwargs)

File /opt/conda/lib/python3.10/site-packages/dolfin/common/plotting.py:279, in _plot_matplotlib(obj, mesh, kwargs)
    277     from mpl_toolkits.mplot3d import axes3d  # noqa
    278     # Enabling the 3d toolbox requires some additional arguments
--> 279     ax = plt.gca(projection='3d')
    280 else:
    281     ax = plt.gca()

TypeError: gca() got an unexpected keyword argument 'projection'
  1. Is there a way to do what you did using a SubMesh of the 3D mesh?
5

You can also save the solution as a .pvd or .xdmf file and read it back in paraview. Then you use the slice filter and export the 2D data as a .csv file. Then you can read and manipulate it with pandas.

6

Try downgrading the matplotlib version.

7

Yep should work. However this should be downgrading to a very old version of matplotlib. Any chance this is another case where we should PR to dolfin to fix the deprecated syntax?

8

Its fine to make a PR to update it, we can get that merged rather quickly. Just note that only minor patches to external packages will be made.

If you want 3D plotting I would recommend vedo or pyvista (not matplotlib). I’ve made a recent extension to interface with pyvista

http://jsdokken.com/dolfin_pyvista_adapter/docs/demo.html

1 Like
9

Thanks for this Jorgen!

I know matplotlib is not ideal in this situation.
I’m trying to have 3D visualisations in a jupyter notebook in a Github Codespace for a workshop.
I’ve given a go at vedo but couldn’t get it to work in Codespaces.

Does dolfin_pyvista_adapter work in jupyter notebooks?

10

Yes, as shown in dolfin_pyvista_adapter/docs/demo.ipynb at main · jorgensd/dolfin_pyvista_adapter · GitHub

1 Like