Access nodes and Dirichlet dofs indices in FESpace consistently with matrix assembly order

[IvanBioli]

Hi Gridap team,

I’m working with an $H^1$-conforming Lagrangian FE space and need to compute the nodes (i.e. dofs coordinates) in the same order as used internally for matrix assembly in get_matrix and get_vector. Specifically, I’d like to achieve two things:

1. Retrieve a matrix dofs_nodes of size (2, ndofs) that contains the coordinates of the dofss.
2. Extract a vector of indices, dirichletdofs, corresponding to the dofs on the Dirichlet boundary.

The goal is to explicitly assemble the system matrices, taking into account Dirichlet constraints without relying on Gridap's built-in assembly adjustments. In the following code snippet, I need to ensure A_cond == Acond_gridap and b_cond == bcond_gridap:

# Mesh and space setup
domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain, partition)
order = 3
reffe = ReferenceFE(lagrangian, Float64, order)

# System assembly (without Dirichlet constraints)
V0 = TestFESpace(model, reffe; conformity=:H1)
Ug = TrialFESpace(V0)
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)
a_(u, v) = ∫(∇(v) ⋅ ∇(u)) * dΩ
b_(v) = ∫(v * f) * dΩ
op = AffineFEOperator(a_, b_, Ug, V0)
A = get_matrix(op)  # Stiffness matrix
b = get_vector(op)  # Load vector

# System assembly (with Dirichlet constraints)
V0 = TestFESpace(model, reffe; conformity=:H1, dirichlet_tags=["boundary"])
Ug = TrialFESpace(V0, g)
op = AffineFEOperator(a_, b_, Ug, V0)
Acond_gridap = get_matrix(op)
bcond_gridap = get_vector(op)

# Explicitly imposed Dirichlet constraints
dofs_nodes = ???        # Coordinates of DOFs
dirichletdofs = ???     # Indices of Dirichlet DOFs
freedofs = ???     # Indices of Free DOFs
uh = zeros(ndofs)
uh[dirichletdofs] = dropdims(mapslices(g, dofs_nodes[:, dirichletdofs]; dims=1), dims=1)
A_cond = A[freedofs, freedofs]
b_cond = b[freedofs] - A[freedofs, dirichletdofs] * uh[dirichletdofs]

I had a look at Tutorial 13: Low-level API Poisson equation, but I did not find much related to this. I understood that negative indices are used for constrained dofs and positive indices for free dofs, but then how do these compare to the ordering employed in assembling the matrices.

Thanks!