Get boundary values

Hello,

I would like to access the values of a function a boundary. Here is a MWE in case of a MixedElement representing for example the velocity and the pressure:

import numpy as np
from mpi4py import MPI
import ufl
import dolfinx
from dolfinx import fem, mesh

msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = fem.functionspace(msh, ufl.MixedElement([ufl.VectorElement('P', msh.ufl_cell(), 2), ufl.FiniteElement('P', msh.ufl_cell(), 1)]))
W1 = V.sub(1)
Q, _ = W1.collapse()

facets_inflow = mesh.locate_entities_boundary(msh, 1, lambda x : np.isclose(x[0], 0.0))
bdofs_inflow = fem.locate_dofs_topological(Q, msh.topology.dim - 1, facets_inflow)
print(bdofs_inflow)

function = fem.Function(V)
p = function.sub(1).collapse()
print(p.x.array[bdofs_inflow])

This works fine, the script outputs the value of the pressure at the left boundary of the unit square.

Let’s try to do the same thing for with VectorElement:

import numpy as np
import pylab as plt
from mpi4py import MPI
import ufl
import dolfinx
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType

msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.FunctionSpace(msh, ufl.VectorElement("DG", msh.ufl_cell(), 2, 3))

W1 = V.sub(0)
Q, _ = W1.collapse()

facets_inflow = mesh.locate_entities_boundary(msh, 1, lambda x : np.isclose(x[0], 0.0))
bdofs_inflow = fem.locate_dofs_topological(Q, msh.topology.dim - 1, facets_inflow)
print(bdofs_inflow)

function = fem.Function(V)
p = function.sub(0).collapse()
print(p.x.array[bdofs_inflow])

In this case bdofs_inflow is empty, so I’m doing something wrong. Could you, please, help me to correct this MWE to get the value of p at the left boundary?

The issue isn’t that it is a VectorElement, but a discontinuous element. Such elements have no degrees of freedom associated with their facets.

You would have to use locate_dofs_geometrical to get those.

Thanks for the help!

What I did is to project first the solution on a continuous space and then I could get the solution at the boundary.

That is usually not a good idea, if the solution is truely discontinuous.

So far the solution is continuous. If it becomes discontinuous, I was thinking about (how?) creating a mesh of dimension 1 corresponding to the boundary and interpolating on this mesh.