In FEniCS there are three different kinds of function space which we can use.

FS = FunctionSpace(mesh, "CG", 1)
VFS = VectorFunctionSpace(mesh, "CG", 1)
TFS =  TensorFunctionSpace(mesh, "CG", 1)

As the name implies:

  • FunctionSpace is to define scalar functions over the mesh and have one value per vertex.
  • VectorFunctionSpace is used to define vectors over the mesh and have d values over the vertex where d is the dimension of the problem. Thus in one dimension it will return a scalar, but it is necessary to define the quantities such as gradients on VectorFunctionSpace which are vectors in 2D and 3D.
  • TensorFunctionSpace is used to define tensor functions over the mesh.

In this post I am going to talk about FunctionSpace and VectorFunctionSpace that I mostly use to define scalar or vector fields over the mesh. Suppose, I have the following mesh which is discretized with traingular elements haveing two degrees of freedom at each node.

image-20210111185517604

This mesh is generate with the help of following commands

from dolfin import *
%matplotlib inline
mesh = UnitSquareMesh(1,1)

Next, we define the displacement vector field over the mesh. This could be achieved by first defining a VectorFunctionSpace and then defining the Functionover the function space.

VFS = VectorFunctionSpace(mesh, "CG", 1)
unew, uold = Function(VFS), Function(VFS)

We can then manually define the vector field based on custom data. Under the hood - ignoring all the complex methods present in the function - a fucntion contains the values of the field at the degrees of freedom in the order

[x0,y0,x1,y1,x2,y2...,xn,yn]

where n is the degree of freedom number. In the above mesh we have 2 degrees of freedom per vertex. Thus, we can modify the vectors

unew.vector().vec().array = np.array([0,0,0,0,1,0,2,0])
uold.vector().vec().array = np.array([0,0,0,0,0,0,0,0])

We can then find the error field over the mesh. Now, since we are finding the error between two vector fields the error will also be a vector field.

error_u = project(unew-uold,VFS)

The vector error_u will contain two values at each vertex, one for error in x and other for error in y. We can check the values by using the method compute_vertex_values

error_u_x, error_u_y = error_u.split()
error_u_x.compute_vertex_values()
error_u_y.compute_vertex_values()

Now, I want the magnitude of the error at each vertex. For this I can write an expression in UFL for magnitude and then project that to a FunctionSpace since I want a single value per vertex

FS = FunctionSpace(mesh, 'CG', 1)
error_u_mag = project(sqrt(inner(error_u,error_u)),FS)

This will give me a field containing the magnitude of the error at each vertex, now to average the error over each element I can take help of Discontinuous Galerkin element

DG = FunctionSpace(mesh, "DG", 0)
error_u_mag_cell =interpolate(error_u_mag,DG)
error_u_mag_cell.vector()[:]

Out: array([ 1. , 0.33333333])

This way I will get the error magnitude over each cell in the mesh.

image-20210111192255567