how to export Gradient of function to VTK file


130
views
0
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):
        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
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.
0
11 weeks ago by
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.

    # 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
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 »