Skeleton for faces of a particular type in CartesianModel

[santiagobadia]

@fverdugo In order to implement variational space-time FE formulations with the recently developed CDLagrangianRefFE, we need a method that given the extrusion type extracts all faces in the skeleton with that orientation in a CartesianModel. This way, we can integrate the DG terms in the time faces. We could consider other ways to define this set of faces, e.g., by the non-zero component in the normal vector, etc, but I think extrusion is better. We can discuss how to implement this.

[fverdugo]

This should be easy.

At low level, we only need to build a vector of masks signaling which facets we want to include in the skeleton triangulation and then use SkeletonTriangulation(model,masks) which i think it is already implemented.

We have different options to build this mask, and using the extrusion is one. Perhaps a more general approach would be to detect in which facets a given interpolation is discontinuous, using the face own dof info of the reffes. This should be as efficient as using the extrusion since the information is compressed at the reffe level in both cases.

That is, given a model + fe interpolation on top of this, get a skeleton triangulation only for the facets where the interpolation is discontinuous. This is actually what you need in dg methods and it is more general than relying on the extrusion for the particular case of cartesian models.

What do you think?

[fverdugo]

The user API could be something like:

SketetonTriangulation(model,reffe)

Or

SkeletonTriangulation(model,reffes)

depending if you have one or several reffes.

[santiagobadia]

Hi @fverdugo

I don't store info that can be easiy consumed for this. I could store that info in a CDLagrangianRefFE struct, and use it here. Something like (CONT,CONT,DISC). But we do not need this struct now, we have enough with a GenericRefFE, where I don't have such info. As a result, I was thinking about

T = VectorValue{2,Float64}
reffe = CDLagrangianRefFE(T,QUAD,order,(DISC,CONT))
V = FESpace(model=model,reffe=reffe)
SketetonTriangulation(model,face_mask=(false,true))

Otherwise, we could assume that an empty face (no owned dofs) in this method is a face in which we want to integrate.

[fverdugo]

Sorry, the continuity info can be taken from "face dofs" instead of "face own dofs". Faces, whose "face dofs" are empty, they can be treated as discontinuous.

That is, this would return the "discontinuous facets" of the reffe:

isempty.(get_face_dofs(reffe,num_dims(reffe)-1))

BTW, I just realized, that the CDLagrangianRefFE has the same "face dofs" that the corresponding continuous lagrangian RefFE. I would say this is not correct.

[santiagobadia]

We have agreed that we want to separate the RefFE from the notion of conformity. It will require some developments to be described in a separate issue.

[fverdugo]

SkeletonTriangulation for "discontinuous facets" done! PR on its way.

using Gridap.ReferenceFEs
using Gridap
using Gridap.FESpaces

reffe = CDLagrangianRefFE(Float64,QUAD,(2,2),(CONT,DISC))
# Manual fix
reffe.face_own_dofs[7] = [1,3,7]
reffe.face_own_dofs[8] = [2,4,8]
reffe.face_own_dofs[9] = [5,6,9]

# Skeleton triangulation
face_own_dofs = get_face_own_dofs(reffe)
strian = SkeletonTriangulation(model,reffe,face_own_dofs)
ns = get_normal_vector(strian)
writevtk(strian,"strian",cellfields=["normal"=>ns])

# Random function for visualization purposes
model = CartesianDiscreteModel((0,1,0,1),(10,5))
V = FESpace(model=model,reffe=reffe,conformity=true)
trian = Triangulation(model)
vh = FEFunction(V,rand(num_free_dofs(V)))
writevtk(trian,"trian",nsubcells=20,cellfields=["vh"=>vh])

[fverdugo]

I close the issue. Reopen it if you wish.