Plotting vector fields with a resolute mesh

mesh
1

Hello,

I am new to FEniCS and I am trying to solve some Poisson problem in curved geometry.
My code is working fine to solve the weak formulation for the pressure field, however it requires a ‘highly’ resolute mesh.

Therefore when it comes to plotting the associated velocity field, there are too much quivers and they can’t be distinguished.
I would like to know if there is a way to plot only 1 quiver every 5 quivers for instance (like in a ‘typical’ Python code)?

Here is the full working code :

from __future__ import print_function
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np

#create mesh
N = 100
L = 10
domain = Rectangle(Point(-L/2,-L/2),Point(L/2,L/2))
mesh = generate_mesh(domain, N)
d = mesh.topology().dim() # dimensionality of the problem
markers = MeshFunction("size_t", mesh, d , mesh.domains())

#plot(mesh,linewidth=0.3)
#plt.show()

#define function space
degreeElements = 3
FS = FunctionSpace(mesh, 'Lagrange', degreeElements)
VFS = VectorFunctionSpace(mesh, 'Lagrange', degreeElements)

#define parameters of the problem
h_0 = 1
sigma = 1
u_0 = 1
K = 1

h = Expression('h_0*exp(-(pow(x[0],2)+pow(x[1],2))/(2*pow(sigma,2)))',degree=degreeElements, h_0=h_0, sigma=sigma)

dh_x = Expression('-x[0]/pow(sigma,2)*h',degree=degreeElements, h=h, sigma=sigma)
dh_y = Expression('-x[1]/pow(sigma,2)*h',degree=degreeElements, h=h, sigma=sigma)

g = Expression('sqrt(1+pow(dh_x,2)+pow(dh_y,2))',degree=degreeElements, dh_x=dh_x, dh_y=dh_y) 

alpha = Expression('pow(dh_x,2)/g',degree=degreeElements, dh_x=dh_x, g=g) 
beta = Expression('dh_x*dh_y/g',degree=degreeElements, dh_x=dh_x, dh_y=dh_y, g=g)  
gamma = Expression('pow(dh_y,2)/g',degree=degreeElements, dh_y=dh_y, g=g)  


#impose Dirichlet boundary conditions

def boundary(x, on_boundary):
    return on_boundary

P_bc= Expression('-u_0/K*(x[0]-L)',degree=degreeElements, u_0=u_0, K=K, L=L)
bc = DirichletBC(FS, P_bc, boundary)

#define function u and test function v
P = Function(FS)
v = TestFunction(FS)

# weak formulation of the problem
Res = -dot(g*grad(P),grad(v))*dx + alpha*P.dx(0)*v.dx(0)*dx + beta*P.dx(1)*v.dx(0)*dx + beta*P.dx(0)*v.dx(1)*dx + gamma*P.dx(1)*v.dx(1)*dx


# solve the problem
solve(Res == 0, P, bc)

# plot solution
c = plot(P,mode='color',title='$u$')
plt.colorbar(c)
#plot(mesh,linewidth=0.3)
plt.show()

# plot v
P_x  = project(P.dx(0),FS)
P_y= project(P.dx(1),FS)

alpha = project(alpha,FS)
beta = project(beta,FS)
gamma = project(gamma,FS)

u = Expression(('-K*P_x + K*( alpha*P_x + beta*P_y)', 
                '-K*P_y + K*( gamma*P_y + beta*P_x)'), P_x=P_x, P_y=P_y, K=K, g=g, alpha=alpha, beta=beta, gamma=gamma, degree=degreeElements)

n_x = Expression(('1.0','0.0'), degree=degreeElements)

n_y = Expression(('0.0','1.0'), degree=degreeElements)
u_x = np.asarray(project(dot(u,n_x),FS))
u_y = np.asarray(project(dot(u,n_y),FS))

sigma=project(u,VFS)
c=plot(sigma, title='$v$',width=.008)
plt.colorbar(c)
#plot(mesh,linewidth=0.3)
plt.show()

Thanks in advance.

2

Please post a fully working minimal example. I will then fiddle with the plotting part.

3

In general, custom visualization is better left to graphical interfaces such as paraview. Save your velocity field to a xdmf or pvd file and open them in paraview.

4

I edited my first message with the fully working code, thanks for your answer

5

Ok I will dive into Paraview then, thanks for your answer