plot a section of a 3d solution

Question

Hi,

I wrote the following code:

from dolfin import *

mesh = BoxMesh (-4,-4,-4,4,4,4,10,10,10)

V = FunctionSpace (mesh, 'CG',1)

dt = 1
rho = 0.012  #giorni^(-1)
G = 0.0013   #giorni^(-1)
D = 0.0065   #cm^2 \ giorni
tf = 60      #giorni

#Initial condition
u_t0 = Expression ('exp(-(x[0]*x[0])/(2*sig0*sig0)-(x[1]*x[1])/(2*sig0*sig0)-(x[2]*x[2])/(2*sig0*sig0))',sig0=3, cell = tetrahedron)

u = TrialFunction (V)
v = TestFunction (V)
f = u_t0 - 0.5*dt*(G-rho)*u_t0  #OK
a = (u*v + dt*D*0.5*inner(nabla_grad(u),nabla_grad(v)) + dt*0.5*(G-rho)*u*v)*dx  #OK
L = f*v*dx - dt*0.5*D*inner(nabla_grad(u_t0),nabla_grad(v))*dx  #OK

u = Function(V)

for i in range(1,tf):
    solve(a == L, u)
        if (i!=(tf-1)):
            u_t0 = u
            f = u_t0 - 0.5*dt*(G-rho)*u_t0
            L = f*v*dx - dt*0.5*D*inner(nabla_grad(u_t0),nabla_grad(v))*dx

file = File('gbm.pvd')
file << u

plot(mesh)
plot(u,interactive=True,title="GBM")

Now I'd like to see the section of the solution, for example at z=0 or z=k with k in [-4,4].
I tried with Para View, but when I open gbm.pvd in it, I have all the filter disabled.
How can I do?

Thanks
Valentina

Comments

Eventually, Do you know a way for render the solution continue also inside the box?

Kind regards

Valentina

Answer 1

I think it may be possible to do without Paraview, by interpolating onto a desired mesh. This is my attempt, albeit I am not very confident about this method - in particular that I used the vertex_to_dof_map correctly.

Code

from dolfin import *

mesh = BoxMesh (-4,-4,-4,4,4,4,10,10,10)

V = FunctionSpace (mesh, 'CG',1)

dt = 1
rho = 0.012  #giorni^(-1)
G = 0.0013   #giorni^(-1)
D = 0.0065   #cm^2 \ giorni
tf = 60      #giorni

#Initial condition
u_t0 = Expression ('exp(-(x[0]*x[0])/(2*sig0*sig0)-(x[1]*x[1])/(2*sig0*sig0)-(x[2]*x[2])/(2*sig0*sig0))',sig0=3, cell = tetrahedron)

u = TrialFunction (V)
v = TestFunction (V)
f = u_t0 - 0.5*dt*(G-rho)*u_t0  #OK
a = (u*v + dt*D*0.5*inner(nabla_grad(u),nabla_grad(v)) + dt*0.5*(G-rho)*u*v)*dx  #OK
L = f*v*dx - dt*0.5*D*inner(nabla_grad(u_t0),nabla_grad(v))*dx  #OK

u = Function(V)

for i in range(1,tf):
        solve(a == L, u)
        if (i!=(tf-1)):
                u_t0 = u
                f = u_t0 - 0.5*dt*(G-rho)*u_t0
                L = f*v*dx - dt*0.5*D*inner(nabla_grad(u_t0),nabla_grad(v))*dx

#Create mesh & corresponding function spaces
mesh2 = RectangleMesh(-4,-4,4,4,10,10)
V2 = FunctionSpace (mesh2, 'CG',1)
u2 = Function(V2)
u = interpolate(u, V)

coords = mesh2.coordinates()
dof_map = vertex_to_dof_map(V2)

# Pick your values of z
z = 1.0

for i in range(V2.dim()):
    u2.vector()[dof_map[i]] = u(coords[i][0],coords[i][1],z)

plot(u2)
interactive()

Comments

Thank you for your reply.

Your program is running and it solves my problem.

Do you know a way for render the solution continue also inside the box?

Kind regards

Valentina

---

I personally don't know how to do it, but this question might be of use: Update plot each time step

Answer 2

ParaView should be able to do what you want.
Maybe try to save a Mesh first, and open that in ParaView (it can also be sliced and sectioned).
Or try the ".xdmf" format instead of ".pvd".

Comments

I tried with .xdmf but I have the same problem.
I tried also to save and open a mesh.pvd, but the problem is the same.
Can the problem be the installation of para view?