### 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)
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, [bc,bc2])
vtkfile = File('poisson_solver/solution.pvd')

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

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)

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