I want to understand how the mapping of DOFs of Hdiv elements, e.g. RT or BDM for identical subdomains work in FEniCS. Consider the following MWE:
from dolfin import *
mesh = UnitSquareMesh(2,1,'left')
# consider the left and right parts as 2 subdomains, which are identical.
# let's say the right subdomain is the translation of the left one by (0.5,0)
# plot to visualize the subdomains.
plt.figure(figsize=(6,6))
plot(mesh)
# define the function space
V_vel = FunctionSpace(mesh,'RT',1)
# access the dofmap of the function space
dofmap = V_vel.dofmap()
# get the DOFs coordinates
dof_coordinates = V_vel.tabulate_dof_coordinates() # returns a numpy array of DOF coordinates
print(dof_coordinates)
# this should be the ouput:
# [[0.75 1. ]
# [1. 0.5 ]
# [0.75 0.5 ]
# [0.5 0.5 ]
# [0.75 0. ]
# [0.25 1. ]
# [0.25 0.5 ]
# [0. 0.5 ]
# [0.25 0. ]]
# create functions for left and right subdomains
V_left = Function(V_vel)
V_right = Function(V_vel)
# these are the DOFs indices for left and right subdomains as is evident from the print command above
dof_index_right = [3,5,6,7,8]
dof_index_left = [0,1,2,3,4]
# this is the mapping of coefficients of left subdomain onto the right one according to DOFs mapping
mapping_left_to_right = [1,0,2,3,4]
# assigning some values to V_1 function s.t. it take non-zero values only for the DOFs of left subdomain
for i,f in enumerate(dof_index_left):
V_left.vector()[f] = (i+1)*0.1
# assigning the same coefficients of V_1 to V_2 according to DOFs mapping
for i,f in enumerate(mapping_left_to_right):
V_right.vector()[f] = V_left.vector()[dof_index_left[i]]
print(norm(V_left), norm(V_right))
# this is the output:
# 0.4987484335815002 0.5408326913195984
Why is the norm different? Moreover, upon plotting these two functions, the two fields do not look identical. Should the norms and visualization of fields not be identical?
Assigning values directly to dofs in spaces with vector-valued basis functions is usually not straight-forward.
This is because the degree of freedom for such spaces is defined as a scalar value, dependent on the normal vector of the facet the degree of freedom is related to DefElement
Assigning values directly to the dof is not adviced, as you need to pull them back to the reference element DefElement: Raviart–Thomas using the contravariant Piola map
Note that if you inspect the actual vectors, the content of them is identical, but it has different meaning for the different cells.
For plotting, I strongly advice interpolating into a higher order DG space and visualize with paraview (plotting with matplotlib interpolates into first order continuous lagrange, causing all sorts of issues):
Thanks for your reply. What you mentioned makes totally sense.
However, for what I am trying to accomplish, I have to use DOFs of one subdomain for the other. I want to get identical fields in these identical subdomains in the visual and function norms, where the subdomain I am looking at has non-zero DOFs values for the corresponding DOF coordinates in that subdomain and zero DOFs value for the other one. (as I tried to show in the MWE)
Regarding the pull back contravariant Piola mapping, could you please show how it can be added to the code to get the desired results? I expect Jacobian = [[1, 0], [0, 1]] and det(Jacobian) = 1 for every cell in the MWE. Should it change anything then?
I also saw some posts on using mesh.init_cell_orientations(global_vector) to allign the cell orientations according to the reference element using a global vector. But I am not sure how to helpful it can be for what I am trying to achieve.
Finally, regarding the visualization in Paraview, it did not make any difference for the MWE.
It is not going to fix the problem, as the problem relates to the pull-back of values. Note that the plots will be different if you use DG-2 (and XDMFFile.write_checkpoint). They will of course not be identical.
I’ve not used vector valued basis functions alot in legacy dolfin, so I cant help you much there. I would probably look at the code from dolfin.FunctionSpace.interpolate to try to understand the logic there.
What is the application in question where you need to work directly on the degrees of freedom like this?
It is about solving a problem on a subdomain of a full domain and using the solution of this subdomain for the other subdomains. Since the subdomains are identical and loading conditions/driving forces for the model for each subdomains are also similar, the solution should also be identical for every subdomain.
It works fine for nodal finite elements, but for these vector-valued finite elements, it doesn’t seem to work as expected, although the idea doesn’t change. Think of it as pressure-velocity pair for Darcy’s equation. I get the correct pressure fields when mapping between the subdomains, but velocity seems to act funny in some areas when mapping between the subdomains.
