Basis Functions


333
views
2
6 months ago by
Dear all,
This question is more related to the mathematical aspects.
I have been playing around with some 1D CG elements.
I wanted to see the values of the quadratic shape functions for example.
So following Visualising shape functions I developed the following MWE:

from dolfin import *

mesh = IntervalMesh(4, 0, 1)
V = FunctionSpace(mesh, 'CG', 2)
v = Function(V)

v.vector()[5] = 1 #Phi 2
v.vector()[6] = 1 #Phi 3
v.vector()[3] = 1 # Phi 4

v.vector().array()
mesh2 = IntervalMesh(40, 0, 1)
V2 = FunctionSpace(mesh2, 'CG', 2)
v2 = interpolate(v, V2)
plot(v2)
interactive()
for i in np.linspace(0,1,5):
    plt.axvline(i,0,2, ls='--', c='k')
plt.ylim(-0.25,1.15)
plt.axhline(0, 0, 1, c='k')
plt.show()​
The output shows (here comes the problem because I am a little bit doubtful of what I am in fact plotting) the shape functions for nodes 2, 3 and 4.
The idea was to reproduce this from Langtangen:

My MWE outputs:

From what I see, the values on this output are always zero because I get the absolute sum of all shape functions.
Clearly I am missing a lot of points and anyone that can explain a little this doubts will make a dumb guy really happy.
All the best, Murilo

Community: FEniCS Project

2 Answers


5
6 months ago by
You are setting degrees of freedom for the same function, resulting obviously in the sum of corresponding basis functions. Just do the thing for three different functions, as here

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

mesh = IntervalMesh(4, 0, 1)
V = FunctionSpace(mesh, 'CG', 2)
v1 = Function(V)
v2 = Function(V)
v3 = Function(V)

v1.vector()[5] = 1 #Phi 2
v2.vector()[6] = 1 #Phi 3
v3.vector()[3] = 1 #Phi 4

mesh2 = IntervalMesh(40, 0, 1)
V2 = FunctionSpace(mesh2, 'CG', 2)
v12 = interpolate(v1, V2)
v22 = interpolate(v2, V2)
v32 = interpolate(v3, V2)

plot(v12)
plot(v22)
plot(v32)
#interactive()
for i in np.linspace(0,1,5):
    plt.axvline(i,0,2, ls='--', c='k')
plt.ylim(-0.25,1.15)
plt.axhline(0, 0, 1, c='k')
plt.show()​
Dear Michal,
Thank you very much for your help.
Just to go through and end my doubts and confirm how this piece of code is working (if you have the kind and patience to do so).
So initially I set up a coarse mesh and a function space with piecewise quadratic basis functions.
In this function space I initialize 3 different functions and after that I set a specifically degree of freedom with a value.
Following this I create a second finer mesh and a new function space (and I made some tests and regardless the basis functions I still have the same results regardless the order of the basis functions of this second function space, V2).
The step that I am trying to understand better is the interpolation of a function from the initial function space (V) to the new function space (V2).
From my perspective I expected that in the interpolation the values of the degrees of freedom on each position are set to be the same. My doubt is how the information of the original basis function (from the function space V) are "passed" to this new function space V2.
This is my last doubt regarding this!
Thank you very much.
All the best, Murilo

written 6 months ago by Murilo Moreira  
1
Interpolation is not just copying the values of degrees of freedom from one vector to another. FEniCS's interpolation is exactly the interpolation operator as you know from theory of finite elements. For more explanation see https://fenicsproject.org/qa/6832/what-is-the-difference-between-interpolate-and-project/
In short: CG2 space needs also value of function in midpoint of element (1D), so you are evaluating the coarse CG2 function at all points where degrees of freedom of finer CG2 function are situated.
written 6 months ago by Michal Habera  
Dear Michal,
Thank you for the explanation and the reference!
Just in case anyone needs, it would also work if I evaluate the value of v1 at each point of the interval!
Turns out that the interpolation was not necessary!
I will now just study a little bit more on the interpolation operator.
Thank you for the link.
All the best, Murilo
written 6 months ago by Murilo Moreira  
3
6 months ago by
Hi,

I would consider playing around with the methods evaluate_basis and/or evaluate_basis_all. Consider the following MWE for plotting the (individual) shape functions on a prescribed cell in your 1D-mesh using evaluate_basis_all

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np 

# Define mesh and function space
mesh = IntervalMesh(4, 0, 1)
V = FunctionSpace(mesh, 'CG', 2)

# Visualize shape functions on cell 1
cidx  = 1 

# Sample shape functions at points 
xcoords  = np.linspace(0,1,20001)

element         = V.element()
num_subspace    = element.num_sub_elements()
space_dimension = element.space_dimension()

shapefun_values = []
xhit            = []

for x in xcoords:
    cell_id = mesh.bounding_box_tree().compute_first_entity_collision(Point(x))
    
    # If point in cidx, then collect shape function info
    if cell_id == cidx:
        cell = Cell(mesh, cell_id)
        coordinate_dofs = cell.get_vertex_coordinates()

        basis_matrix = np.zeros(space_dimension, dtype=float)
        
        # Evaluate all the basis functions at point x
        element.evaluate_basis_all(basis_matrix,
                                    np.array(x),
                                    cell.get_vertex_coordinates(), 
                                    cell.orientation())
        shapefun_values.append(basis_matrix)
        xhit.append(x)

# Plot shape functions at cell cidx        
plt.plot(xhit, shapefun_values[:])
plt.show()
    ​

Reproducing the desired figure now only requires some additional 'knitting'. But that shouldn't be too hard.

Hope this helps!

Dear Jakob, thank you very much!
Your MWE was really helpful.
It is a really through example on how to work with the cell object from the mesh.
I just would like to ask you a little bit more on why do we have this difference.
I know that with your methodology we are getting the value of each shape function and then plotting them.
When considering the methodology that I posted I am getting the values of the shape functions all together. I am just not able to see this analytically.
Anyway thank you for your time,
All the best, Murilo
written 6 months ago by Murilo Moreira  
1
Michal's explanation below exactly addresses this issue, so see below ;)
written 6 months ago by Jakob Maljaars  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »