interpolate_boundary_values() will not work because shape functions corresponding to enriched DoFs generally won't satisfy delta Kronecker.
Even if they are by chance, you may have a scalar FE with two shape functions (enriched / non-enriched) at the same support point.
What you could do is L2 projection on the boundary to determine the values of DoFs. That amounts to inverting a mass matrix on the boundary.
That should work.
Otherwise, unless you need to enrich elements at the boundary, i would suggest instead mixing enriched and not-enriched elements.
Depending on your enrichment, you may severely deteriorate the condition number of matrices. At least in some applications it is better to
enriched only locally (in some part of the domain).
Finally, in order to know whether a local DoF corresponds to enriched, you can query finite element for the base number of this DoF:
FiniteElement< dim, spacedim >::system_to_base_table()[local_dof].first.first // base element number
0 is un-enriched.
Regards,
Denis.