How to plot the divergence of a solution object

1

Hello,
I would like to check what is the divergence of the solution function. (for other F forms as well). Here I get the following error:

AttributeError: ‘Div’ object has no attribute '_cpp_object’

I feel like somehow the solution u_ is not a function anymore by the time we get to the end? How can I still somehow plot and check the maximum value of the divergence of u_? How to treat solution objects as if they were still functions?

v = TestFunction(V)
u = TrialFunction(V)
u_ = Function(V)
F = inner(grad(u) * u, v) * dx
F = action(F, u_)
J = derivative(F, u_, u)
problem = NonlinearVariationalProblem(F, u_, bcu, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

from ufl import div
file_div = File(“nonlinear_u_div.pvd”)
file_div << div(u_)

Thank you so much for any suggestions. (By the way
u_1, u_2, u_3 = split(u_)
div = u_1.dx(0) + u_2.dx(1) + u_3.dx(2)
also did not work.)

Best, Nora

2

Hi, i would suggest to project the div(u) to a suitable function space, and then save the resulting function.

3

Thank you for your response. I don’t know if I understand your idea well, I am not sure at which point you mean the projection. div(u) itself I guess does not really make sense because u is a TrialFunction, and the solution goes into u_, not u. Projecting div(u_) on the other hand seems problematic, because it would try to project something that does not exist in the first place. I mean, project(div(u_), V) wouldn’t work because div(u_) is not there. Because div(_cpp_object) does not exist, according to the error message.

4

Sorry that i wasnt clear in my original response. You can consider the following. Project the divegence of a function to a DG space. (In your case projecting u_)

S = VectorFunctionSpace(mesh, "CG",1)
s=Function(S)
V= FunctionSpace(mesh, "DG", 0)
ds = project(div(s),V)
1 Like
5

Thank you. I tried it and it does seem to work.

6

Hello,

I have a divergence related problem, I would like to plot the divergence of gradient of a scalar field. But, I’m getting zero for any scalar field. As an example, I assumed a simple polynomial y = x^2 + y^2, where its gradient is (2x,2y) and I expect a constant scalar field of 4 for the divergence of gradient, but I’m still getting zero answer. Below is my simple code for this example:
I would be appreciated if you could help me!

from fenics import *
from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
mesh = RectangleMesh(Point(0,0), Point(100, 100), 10, 10)

V1 = FunctionSpace(mesh, “DG”, 1)
V2 = FunctionSpace(mesh, “DG”, 0)

m = Expression ( "x[1]*x[1] + x[0]*x[0] ", degree = 2 )
m = project ( m, V1)

plot(m)
plt.colorbar(plot(m))
plt.show()

grad_m = (grad(m))
plot(grad_m)
plt.colorbar(plot(grad_m))
plt.show()

div_grad_m = div(grad(m))
div_grad_m = project(div_grad_m,V2)
plot(div_grad_m)
plt.colorbar(plot(div_grad_m))
plt.show()

7

Hi,
when you project m to V1 it gets approximated by a piecewise affine function, hence div(grad(m)) becomes zero.
You should consider using a space with quadratic interpolation instead.