Getting element/cell ID

[MichaelLoibl]

Hi,

I want to compare some results with another code. In this context, it would be helpful to get the quadrature points in combination with the respective global element ID. Are there element IDs in Gridap?

I have the following code below and it would be helpful to get the element ID for each cell data Qₕ_cell_data[i] (resp. q). Thank you for your help!

module example_circle_1

using Gridap
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.CellData
using Gridap.Fields
using GridapEmbedded
using GridapEmbedded.CSG
using GridapEmbedded.LevelSetCutters
using Printf

#using a cellfield
domain = (0,1,0,1)
n = 5
bgmodel = CartesianDiscreteModel(domain,(n,n))

reffe = ReferenceFE(lagrangian,Float64,1)
Ω_bg = Triangulation(bgmodel)
V_bg = FESpace(Ω_bg,reffe)
φh = interpolate(x->(x[1] - 1/2)^2 + (x[2] - 1/2)^2 - 1/25,V_bg)
geo = DiscreteGeometry(φh,bgmodel)

# Cut the background model
cutgeo = cut(bgmodel,geo)

Ω_act = Triangulation(cutgeo,ACTIVE)

writevtk(Ω_bg,"bg_trian")
writevtk(Ω_act,"act_trian")

# subtriangulating cut elements
Ω = Triangulation(cutgeo,PHYSICAL)
writevtk(Ω,"phys_trian")

order = 1
degree = 2*order
dΩ = Measure(Ω,degree)

# compute area
area = sum(∫(1.0)*dΩ)

# extract quadrature points
Qₕ = CellQuadrature(Ω,degree)
Qₕ_cell_data = get_data(Qₕ)
num_cell = length(Qₕ_cell_data)
q = Qₕ_cell_data[rand(1:num_cells(Ω))]
p = get_coordinates(q)
w = get_weights(q)

end # module

[JordiManyer]

Can you describe what you mean by element id? Are you referring to cells/faces/dofs/points?

[MichaelLoibl]

I think that I mean cells in this context. I don't want to confuse it with the subdivision which is done for cut cells (see picture below). I wanted to refer to the cells of the initial mesh without any subdivision of cut cells. Thus, I would say for example that the mesh below has 5x5 elements.

[JordiManyer]

All information on the subcell-to-cell map is in the SubCellData structure. In particular, you are asking for the cell_to_bgcell array. From a SubCellTriangulation (i.e your Ω), you can access it through Ω.subcells.cell_to_bgcell. It will also be part of your FaceToFaceGlue when you call get_glue(Ω,Dc) with Dc your cell dimension. In case of doubt, just go through the code in GridapEmbedded/Interfaces/SubCellTriangulations.jl.