Poisson equationVTK visualization and verification

1

Hi,
I work on the example of Poisson equation from here, Fundamentals: Solving the Poisson equation

My case is easy, all the boundary conditions are dirichlet boundary conditions, and they all equal to 0,
the equation is
\div u =10

For the boundary, I am using u_D = Expression(‘0’, degree=0), does the degree equal to zero, if the boundary condition is a constant？
Since this is 3d so I am using unitcubemesh and cg for the functionspace

from fenics import *

# Create mesh and define function space
mesh = UnitCubeMesh(8, 8, 8)
V = FunctionSpace(mesh, 'CG', 1)

# Define boundary condition
u_D = Expression('10', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(10)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

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

# Plot solution and mesh
plot(u)
plot(mesh)

# Save solution to file in VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
print('error_L2  =', error_L2)
print('error_max =', error_max)

I got warnings and relative high norms like

*** Warning: Degree of exact solution may be inadequate for accurate result in errornorm.
error_L2 = 0.23352226896393796
error_max = 0.5491766911624044

Process finished with exit code 0

and for paraview, how to see the plot inside the cube?

image
image961×680 14 KB

2

Hi @Miraboreasu,

you can use the slice filter in paraview to visualise the interior of a solid volume.

fig_for
fig_for1072×755 70.1 KB

I obtain that for the solution of your problem.

About the warning, it dissapears if you raise the order of the analytical solution although it does not affect the resulting error as u_D is constant. You could also write the boundary condition as:

u_D = Constant(10)
1 Like
3

If I want to defince the length of the 3D domain, what function I should use?

4

To my knowledge there is no built in mesh for 3D rectangular prisms, if your problem has a symmetry you can reduce it to 2D and use RectangleMesh. If this is not the case I think you should build a mesh with some other software (gmsh for example) and import it.

5

You can use BoxMesh as shown in: Built-in meshes — DOLFIN documentation

1 Like