labeling of CartesianDiscreteModel

[kit-simmonds]

Could some one talk me through, or point me to the right resources for understanding, the default labelling of the CartesianDiscreteModel?

I have been looking at the hyperelasticity demo https://github.com/gridap/Tutorials/blob/master/src/hyperelasticity.jl, specifically the lines

# Define new boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri_0",[1,3,7])
add_tag_from_tags!(labels,"diri_1",[2,4,8])

and have played around with the numbers to find

# Labelling's: 1: corner [0,0], 2: corner [1,0], 3: corner [0,1], 4: corner [1,1]
#            5: bottom edge, 6: top edge, 7: left edge, 8: right edge - (I think!).

My problem is that I would like to extend this to a 3D cylindrical mesh and I'd like to be able to use gridap to set my BCs but I don't know the labelling regime.

Is there a simple way to set labels myself according to some condition?

Thanks,
Kit

[JordiManyer]

This is the place where all the face labeling happens: https://github.com/gridap/Gridap.jl/blob/master/src/Geometry/FaceLabelings.jl

There are two main ways of creating tags from existing tags:

• Gridap.jl/src/Geometry/FaceLabelings.jl

  Line 216 in 3b87738

  function add_tag!(lab::FaceLabeling,name::String,entities::Vector{<:Integer})

• Gridap.jl/src/Geometry/FaceLabelings.jl

  Line 228 in 3b87738

  function add_tag_from_tags!(lab::FaceLabeling, name::String, tags::Vector{Int})

If you want to create completely new tags, or start a FaceLabeling objects from scratch, you will have to to the face tagging yourself. This can be quite hands on, and will depend on your problem/mesh.

A taste of it: The following function update_labels! modifies the labels of a model by adding a new label name to all entities where f_Γ(x) return true for all its defining nodes. Note that this will destroy your initial tags. If you want to not do that, there is extra work to be done.

function update_labels!(e::Integer,model::DiscreteModel,f_Γ::Function,name::String)
  mask = mark_nodes(f_Γ,model)
  _update_labels_locally!(e,model,mask,name)
  nothing
end

function _update_labels_locally!(e,model::DiscreteModel{2},mask,name)
  topo   = get_grid_topology(model)
  labels = get_face_labeling(model)
  cell_to_entity = labels.d_to_dface_to_entity[end]
  entity = maximum(cell_to_entity) + e
  # Vertices
  vtxs_Γ = findall(mask)
  vtx_edge_connectivity = Array(get_faces(topo,0,1)[vtxs_Γ])
  # Edges
  edge_entries = [findall(x->any(x .∈  vtx_edge_connectivity[1:end.!=j]),
    vtx_edge_connectivity[j]) for j = 1:length(vtx_edge_connectivity)]
  edge_Γ = unique(reduce(vcat,getindex.(vtx_edge_connectivity,edge_entries),init=[]))
  labels.d_to_dface_to_entity[1][vtxs_Γ] .= entity
  labels.d_to_dface_to_entity[2][edge_Γ] .= entity
  add_tag!(labels,name,[entity])
  return cell_to_entity
end

function _update_labels_locally!(e,model::DiscreteModel{3},mask,name)
  topo   = get_grid_topology(model)
  labels = get_face_labeling(model)
  cell_to_entity = labels.d_to_dface_to_entity[end]
  entity = maximum(cell_to_entity) + e
  # Vertices
  vtxs_Γ = findall(mask)
  vtx_edge_connectivity = Array(Geometry.get_faces(topo,0,1)[vtxs_Γ])
  vtx_face_connectivity = Array(Geometry.get_faces(topo,0,2)[vtxs_Γ])
  # Edges
  edge_entries = [findall(x->any(x .∈  vtx_edge_connectivity[1:end.!=j]),
    vtx_edge_connectivity[j]) for j = 1:length(vtx_edge_connectivity)]
  edge_Γ = unique(reduce(vcat,getindex.(vtx_edge_connectivity,edge_entries),init=[]))
  # Faces
  face_entries = [findall(x->count(x .∈  vtx_face_connectivity[1:end.!=j])>2,
    vtx_face_connectivity[j]) for j = 1:length(vtx_face_connectivity)]
  face_Γ = unique(reduce(vcat,getindex.(vtx_face_connectivity,face_entries),init=[]))
  labels.d_to_dface_to_entity[1][vtxs_Γ] .= entity
  labels.d_to_dface_to_entity[2][edge_Γ] .= entity
  labels.d_to_dface_to_entity[3][face_Γ] .= entity
  add_tag!(labels,name,[entity])
  return cell_to_entity
end

function mark_nodes(f,model::DiscreteModel)
  topo   = get_grid_topology(model)
  coords = get_vertex_coordinates(topo)
  mask = map(f,coords)
  return mask
end

Usages of this can be found for instance here:

https://github.com/zjwegert/GridapTopOpt.jl/blob/ded39d8dd0f976d3c0830cde183b343236bcd7bb/test/seq/ThermalComplianceALMTests.jl#L39