how to export Gradient of function to VTK file

11 weeks ago by
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:

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):
        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 - \

            # Check maximum error
            msg = 'error_max = %g' % error_max
            assert error_max < tol, msg

if __name__ == '__main__':
Community: FEniCS Project
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

11 weeks ago by
The expression 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.
11 weeks ago by

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  
11 weeks ago by
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
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.

Similar posts:
Search »