Assign vertex values to function with Lagrange element of degree higher than 1

1

Dear Fenics community,

I am trying to assign vertex values to function in function space with Lagrange element of degree higher than 1. I am aware of the method to do this when the degree of Lagrange element is 1 as the following.

import numpy as np
from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'P', 1)
u = Function(V)

a = np.random.rand(mesh.num_vertices())
d2v = dof_to_vertex_map(V)
a_new = [a[d2v[i]] for i in range(mesh.num_vertices())]
u.vector().set_local(a_new)
np.array_equal(u.compute_vertex_values(), a)

How to do it when the degree of Lagrange element is higher than 1?

Thanks in advance.

2

Why would you set values to only the vertices of a higher order element? Shouldn’t you set values to all dofs, including those on the faces,edges and interior (depending on the degree of your function space).

3

@dokken Thanks for your response. My original problem is to compute the Laplacian of mean curvature of curved surface \nabla^2 H, a term appearing in the equation governing the shape of cell membrane. The known condition is the coordinates of mesh vertices of curved surface. By using the discrete operator defined in this paper, the mean curvature (as well as Gauss curvature) can be calculated vertex-wise. So I have obtained the mean curvature on vertices.

To compute the Laplacian of mean curvature using Fenics, I think something like

V = FunctionSpace(mesh, 'CG', 3)
V1 = FunctionSpace(mesh, 'DG', 1)
H = Function(V)
Laplacian_H = project(div(grad(H)), V1)

will do the job. The values of function H on vertices should agree with those computed previously. That’s why I want to set values on vertices of higher order element. Other dofs are not constrained, my guess is that they can be specified arbitrarily. I am new to Fenics and not that familiar with finite element method. Correct me if I’m wrong. Thanks.

4

If you want the curvature, I have posted an alternative approach here: How to compute curvature of a boundary

5

Hi, @dokken. If I understand it correctly, the main idea in your post is to compute the normal vector of boundary mesh by constructing a volume mesh. It may well apply to the case where the surface is regular like sphere. But for general surface, I don’t know if the generation of volume mesh from surface mesh is computationally costly or not. If the computational cost is high, then it may be not suitable for the problem where the surface are evolving, which is the case in my problem.

I am wondering if it is possible to compute the normal vector directly from the boundary mesh. If not, what prevents us to do so?

6

To me it would seem very strange if the other dofs could be specified arbitrarily.
However, if you just want the vertex dofs:

V.dofmap().entity_dofs(mesh, 0)

will give you the dofs of the vertices.

7

Now sounds strange to me, too :smiley:. I guess the dofs on edges and faces (assuming degree 3) should interpolate those on vertices, right? If yes, how to interpolate? I know if an expression as a function of coordinate can be interpolated like

e_x0 = Expression('x[0]', degree=3)
V = FunctionSpace(mesh, 'P', 3)
x0 = interpolate(e_x0, V)

How to interpolate a discrete function with values only on vertices? Please walk me through this. Thanks.

8

You would have to write a custom interpolation routine for that.
It will look like something like this (pseudocode):

v_3 = interpolate(e_x0, V)
v_custom = np.array(num_dofs_in_v_3, dtype=np.float64)
for cell in cells(mesh):
    cell_basis = V.element.evaluate_basis(....)
    for dof in dofs(cells):
        if dof is vertex_dof:
             v_custom[dof] = cell_basis[dof]
        else:
             for basis in cell_basis:
                 v_custom[dof]+= basis

where you get the basis functions with:
https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/python/src/fem.cpp?at=master#lines-169 as described in https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/python/test/unit/fem/test_manifolds.py?at=master#lines-235,235:243

Note that this is a non-trivial interpolation, and that the code above only illustrates the idea, and will not work out of the box

2 Likes