Visualising shape functions

Linear Algebra
1

Elaborating on Visualising shape functions - FEniCS Q&A

Here is my code, it fails for some cases. Results in paraview look better using “warp by scalar”.

# Use this, for example, as:
#
# $ python3 visualize_basis.py CG 2

import dolfin
import sys

family = sys.argv[1]
degree = int(sys.argv[2]) if len(sys.argv) >= 3 else None

mesh = dolfin.UnitSquareMesh(5,5)
pvd_mesh = dolfin.File('results/mesh.pvd')
pvd_mesh << mesh

V = dolfin.FunctionSpace(mesh, family, degree)
v = dolfin.Function(V)

mesh2 = dolfin.UnitSquareMesh(50,50)
V2 = dolfin.FunctionSpace(mesh2, 'CG', 1)

pvd_file = dolfin.File(f'results/basis_{family}_{degree}_.pvd')

N = len(v.vector())
print(f"This basis has {N} elements.")

for i in range(N) :
    if N > 10 and i % (N//10) == 0 : print(f"Working on element: {i}")
    if i != 0 :
        v.vector()[i-1] = 0
    v.vector()[i] = 1
    v2 = dolfin.interpolate(v, V2)
    v2.rename ("phi", "phi")
    pvd_file << v2
2

In your code you are interpolating an arbitrary function into a CG 1 function space. In most cases you will then loose a lot of information about the basis function.

Note that using pvd results in a CG-1 interpolation of your field. For arbitrary order CG/DG see: Loading xdmf data back in - #4 by dokken

Also, you can find plots of shape functions and explicit formulas for shape functions at: https://defelement.com/

3

Thank you, I’ll try to improve the code.