### how to export Gradient of function to VTK file

130

views

0

I ran simple poisson equestion solver and it correctly calculated "u" and expored it to VTK file,

I was trying to export grad(u) to VTK file to see vectors, but I have confused with it,

you will find me code below:

I was trying to export grad(u) to VTK file to see vectors, but I have confused with it,

you will find me code below:

```
from __future__ import print_function
from fenics import *
import numpy as np
def solver(f, Nx, Ny, degree=1):
"""
Solve -Laplace(u) = f on [0,1] x [0,1] with 2*Nx*Ny Lagrange
elements of specified degree and u = u_D (Expresssion) on
the boundary.
"""
# Create mesh and define function space
mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'P', degree)
# Define boundary condition
def boundary(X, on_boundary):
x = X[0]
y = X[1]
if x > 0.7 and x < 0.8:
if y > 0.2 and y < 0.8:
return True
return False
def boundry2(X,jafar):
x=X[0]
y=X[1]
if x>0.2 and x<0.3:
if y>0.2 and y<0.8:
return True
return False
u_D2 = Expression('-10.0 * (1 - 10*(x[1] - 0.5) * (x[1] - 0.5))',degree=2)
u_D = Expression('10.0 * (1 - 10*(x[1] - 0.5) * (x[1] - 0.5))',degree=2)
bc = DirichletBC(V, u_D, boundary)
bc2 = DirichletBC(V,u_D2,boundry2)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, [bc,bc2])
vtkfile = File('poisson_solver/solution.pvd')
vtkfile << grad(u)
return u
def run_solver():
"Run solver to compute and post-process solution"
# Set up problem parameters and call solver
# u_D = Expression('1 - 4*x[0]*x[0] + 23*x[0]*x[1]', degree=2)
u_D = Constant(0.0)
f = Constant(-15.0)
u = solver(f, 15, 15, 2)
# Plot solution and mesh
# plot(u)
# plot(u.function_space().mesh())
# Save solution to file in VTK format
# vtkfile = File('poisson_solver/solution.pvd')
# vtkfile << u
def test_solver():
"Test solver by reproducing u = 1 + x^2 + 2y^2"
# Set up parameters for testing
tol = 1E-10
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
f = Constant(0.0)
# Iterate over mesh sizes and degrees
for Nx, Ny in [(3, 3), (3, 5), (5, 3), (20, 20)]:
for degree in 1, 2, 3:
print('Solving on a 2 x (%d x %d) mesh with P%d elements.'
% (Nx, Ny, degree))
# Compute solution
u = solver(f, u_D, Nx, Ny, degree)
# Extract the mesh
mesh = u.function_space().mesh()
# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_D - \
vertex_values_u))
# Check maximum error
msg = 'error_max = %g' % error_max
assert error_max < tol, msg
if __name__ == '__main__':
run_solver()
```

Community: FEniCS Project

1

Could you please properly format your code (using the "Insert code" functionality) so as to ease the process of copying and running your code?

written
11 weeks ago by
klunkean

### 3 Answers

5

The expression

You can project that expression onto a suitable finite element space which then can be exported:

This is the FEniCS internal equivalent to what Matthieu Diaz is doing manually in his answer.

`grad(u)`

is a UFL statement and essentially just a symbolic representation of what is to be done with the expression once it gets evaluated in a variational form. You can project that expression onto a suitable finite element space which then can be exported:

```
V_grad = VectorFunctionSpace(mesh, 'DG', degree-1)
grad_u = project(grad(u), V_grad)
vtkfile << grad_u
```

This is the FEniCS internal equivalent to what Matthieu Diaz is doing manually in his answer.

0

Hello,

You can only register solutions of variational formulations or meshes easily on VTK files. I had the same problem as you previously, I tried to save a non linear application of the solution u of my problem.

The tip is to solve a mini variational problem. Below I show you on my example how I succeeded. I solved a complex problem to get the temperature theta_ and I want to save on a VTK file a function of theta called phi I declared before in my code.

I think you can do the same, introducing a new function u2 for instance and write a problem like this : Fc = dot((u2 - grad(u)),v)* dx where v is a test function.

You can only register solutions of variational formulations or meshes easily on VTK files. I had the same problem as you previously, I tried to save a non linear application of the solution u of my problem.

The tip is to solve a mini variational problem. Below I show you on my example how I succeeded. I solved a complex problem to get the temperature theta_ and I want to save on a VTK file a function of theta called phi I declared before in my code.

```
# Post-treatment : phi(theta_T(x,y))
C = FunctionSpace(mesh, 'P', 1) # a new space for this mini problem
concentration = Function(C,name="Concentration") # the function I will save
c = TestFunction(C) # test function for the mini variational formulation
Fc = dot((concentration-phi(theta_)),c)*dx # The problem here
solve(Fc == 0, concentration)
vtkfile4 << (concentration, t)
```

I think you can do the same, introducing a new function u2 for instance and write a problem like this : Fc = dot((u2 - grad(u)),v)* dx where v is a test function.

Don't remember to tell if this has been able to help you ! :)

written
11 weeks ago by
Matthieu Diaz

This is essentially what

`project(...)`

does (see my answer).
written
11 weeks ago by
klunkean

-2

but these solutions, draw magnitude of grad(u), I want to have gard(u) as vector in VTK file,

as i know VTK file supports Vector too, I need to have vector not ||grad(u)||.

thanks for replies

as i know VTK file supports Vector too, I need to have vector not ||grad(u)||.

thanks for replies

No they don't. The file contains the whole vectorial data.

written
10 weeks ago by
klunkean

Please login to add an answer/comment or follow this question.