Doubt about Fenics solution

General
1

Hello to everyone.
I’m having some issues understanding some concepts in fenics. In the following code, I didn’t enforce any Dirichlet Boundary conditions and the function g doesn’t exist, so my space is still H1 and has a Neumann zero condition on the boundary ? How do the applications of Neumann and Dirichlet boundary conditions work in the sense of H1 and H0 spaces (Sobolev Spaces)? Finally, I’m questioning this because in the following code, without specifying any Dirichlet boundary conditions and setting f as 1, I obtained the following:

image
image1131×865 111 KB

from dolfin import *
from mshr import *

# Create mesh and define function space
domain = Rectangle(Point(-0.5, -0.5), Point(0.5, 0.5))
mesh = generate_mesh(domain, 64)
mesh = refine(mesh)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("1", degree=2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx 

# Compute solution
u = Function(V)
solve(a == L, u)

# Save solution in VTK format
file = File("poisson.pvd")
file << u

# Plot solution
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plot(u)
plt.show()

mesh_coordinates = mesh.coordinates()
x, y = mesh_coordinates[:, 0], mesh_coordinates[:, 1]

# Extract the values of u at mesh vertices
u_values = u.compute_vertex_values(mesh)

# Plot 3D solution of u
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(x, y, u_values, cmap='viridis')
ax.set_title("3D Solution of u")

plt.show()

And I expected the solution closer to 1.
Thanks.

2

You end up in this case https://fenicsproject.org/olddocs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html . The solution is defined up to a constant, and the constant you get is -7 rather than 1.

3

Hello, thank you for responding. But, since I’m considering the ‘g’ function as zero, is it correct to say that a pure Neumann zero condition is being applied at the boundary of the H1 space?

4

Yes, that’s correct, the code you have is equivalent to a homogeneous Neumann boundary condition on the whole boundary of the domain.

1 Like
5

Ok, thanks, it helps me a lot. My last question is, when you define subdomains, do we divide our ‘mesh’ into different disjoint domains? And if I don’t specify the function g in each subdomain, will FEniCS understand it as Neumann zero? To define a function g for each subdomain, do I need to mark each of them with a constant “i” and use ds(i) for each one, right? Sorry for asking this, I was reading the documentation and not understanding these things.

6

in the following I will assume that by subdomains you mean “boundary parts”

yes, you are dividing the faces on the boundary in disjoint parts

Yes, and you see this by writing the weak formulation of the problem. There would be the sum of boundary integrals of g * test * ds(i): if g = 0, the boundary integral vanishes

that’s correct.

1 Like