Accessing dof by face tag

[kutsjuice]

Dear Gridap team,

First of all, I would like to thank you for your work and express my admiration. I recently started using Gridap and have already fallen in love with it completely. It's a really great tool!
I am trying to perform mechanical homogenization on a unit cell with some complex geometry, built using GMSH. The mesh is represented as an unstructured grid. For this task, I need to implement periodic boundary conditions. I plan to use a penalty-based method for this. However, I have run into a problem. How can I correctly retrieve the node indices with specific tags, and how can I get the DOFs associated with these nodes?
Currently, I am using the following function to retrieve node indices:

function get_nodes_by_tag(model::DiscreteModel, tag::String)::Vector{Int64}
    labels = get_face_labeling(model)

    dim = 0 # we will look for nodes which has dimension of 0
    
    dface_to_entity = get_face_entity(labels, dim)    
    dface_to_isontag = BitVector(undef, num_faces(labels,dim))
    tag_entities = get_tag_entities(labels, tag)

    for entity in tag_entities
        for i in eachindex(dface_to_entity)
            dface_to_isontag[i] =  dface_to_entity[i] == entity
        end
    end
    return findall(dface_to_isontag)
end

but with this method, I miss nodes that lie on the boundary of the face.

The second question is whether there is a correct way to obtain the DOFs associated with a node. Currently, I simply use the formulas

x_dof = (node*3 - 2)
y_dof = (node*3 - 1)
z_dof = (node*3)

and it seems to work fine for unconstrained fe spaces, but I'm not sure if it will work for other cases.

And an additional question: As I understand it, FESpacesWithLinearConstraints implements Master-Slave MPC. Does it support non-homogeneous MPC? And is there an example of using this FESpace available, because, unfortunately, I did not find any, and it is not clear from the source code how to use it.

Thank you in advance for helping.

Sincerely,
Micheal

[kutsjuice]

The first part of the question was solved:
After carefully reading GridapGmsh readme.md, I found that I need to create physical groups not only for the 2D entities in gmsh, but for the 1D and 0D entities.
the correct function to get nodes by tag sis looks like as follow now:

function get_nodes_by_tag(model::DiscreteModel, tag::String)::Vector{Int64}

    
    labels = get_face_labeling(model)
    tag_id = get_tag_from_name(labels, tag);
    
    dim = 0 # we will look for nodes which has dimension of 0
    
    dface_to_entity = get_face_entity(labels, dim) 
    
    dface_to_isontag = BitVector(undef, num_faces(labels,dim))

    tag_entities = get_tag_entities(labels, tag)

    for i in eachindex(dface_to_entity)
        buf = false
        for entity in tag_entities
            buf += dface_to_entity[i] == entity
        end
        dface_to_isontag[i] = buf
    end

    return findall(dface_to_isontag)

end

But then I found that this approach don't work for 2nd order meshes - the nodes, which lays on element edges are missed. At this moment, I came up with a temporary solution: selecting nodes by their coordinates directly, skipping the tag information from gmsh.
The second problem, that for the 2nd order meshes, the formulas for dofs:

x_dof = (node*3 - 2)
y_dof = (node*3 - 1)
z_dof = (node*3)

also don't work any more.

So currently I get this infromation from the cell_nodes_ids and cell_dofs_ids arrays.
The function to get nodes in boundin box and routines to select associated dofs are as follows:

function get_nodes_in_bounding_box(model, bounding_box)
    node_coordinates = get_node_coordinates(model)
    nodes_in_bbox = BitVector(undef, length(node_coordinates))

    for i in eachindex(nodes_in_bbox)
        nodes_in_bbox[i] = 
            node_coordinates[i][1] >= bounding_box[1] && 
            node_coordinates[i][1] <= bounding_box[4] &&
            
            node_coordinates[i][2] >= bounding_box[2] && 
            node_coordinates[i][2] <= bounding_box[5] &&
            
            node_coordinates[i][3] >= bounding_box[3] && 
            node_coordinates[i][3] <= bounding_box[6];
    end
    return findall(nodes_in_bbox)
end

node_to_dofs = Vector{Tuple{Int64, Int64, Int64}}(undef, num_nodes(fe_model))
node_is_set = BitVector(undef, num_nodes(fe_model)); fill!(node_is_set, 0);
cell_dofs = get_cell_dof_ids(trial_space_U)
cell_nodes = get_cell_node_ids(fe_model)

for (dofs, nodes) in zip(cell_dofs, cell_nodes)
    for (i, node) in enumerate(nodes)
        if !node_is_set[node]
            node_to_dofs[node] = (dofs[i], dofs[i + length(nodes)], dofs[i + 2*length(nodes)]);
            node_is_set[node] = 1;
        else
            # println(dofs')
            # println(nodes')
            if (dofs[i], dofs[i + length(nodes)], dofs[i + 2*length(nodes)]) != node_to_dofs[node]
                println(dofs')
                println(nodes')
            end
        end
    end
end