Expanding adaptivity to non-BodyFittedTriangulation triangulations

[JordiManyer]

In PR #838 we decided to limit the inter-grid change_domain capabilities to BodyFittedTriangulations. However, we would like to expand Adaptivity to other types of Triangulation, notably BoundaryTriangulation and SkeletonTriangulation (which will be needed for hybrid discretisations).

This is the roadmap towards implementing these goals:

• First, we need to be able to produce AdaptivityGlues which contain adaptivity information on all facet dimensions. We should modify our refine implementation for CartesianDiscreteModel to fill all facet information. This will allow us to test new features as we add them. EDIT: We have decided NOT to implement this for now (see comments in the PR).
• I believe the key to be able to expand adaptivity to all Triangulation subtypes in a general way is to somehow find a way to couple the AdaptivityGlue with Gridap's triangulation glue. Then, two triangulations T1 and T2 are compatible iff get_glue(T1) and get_glue(T2) can be matched using the AdaptivityGlue between their two background models.

[codecov-commenter]

Codecov Report

Merging #868 (f6928a2) into master (672e3d9) will increase coverage by 0.01%.
The diff coverage is 79.26%.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##           master     #868      +/-   ##
==========================================
+ Coverage   88.56%   88.57%   +0.01%     
==========================================
  Files         172      172              
  Lines       19994    20035      +41     
==========================================
+ Hits        17707    17746      +39     
- Misses       2287     2289       +2     

Impacted Files	Coverage Δ
src/Adaptivity/AdaptedTriangulations.jl	82.55% <79.26%> (-0.65%)	⬇️
src/Adaptivity/AdaptivityGlues.jl	89.06% <0.00%> (-3.13%)	⬇️
src/Geometry/SkeletonTriangulations.jl	87.33% <0.00%> (+0.90%)	⬆️
src/Fields/FieldsInterfaces.jl	84.48% <0.00%> (+1.47%)	⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[amartinhuertas]

notably BoundaryTriangulation and SkeletonTriangulation (which will be needed for hybrid discretisations).

Let me just point here that in GridapHybrid.jl we also have a data type called SkeletonTriangulation (see https://github.com/gridap/GridapHybrid.jl/blob/048853981cef5d417c9930224a3d5116c9a2c158/src/Skeleton.jl#L113).

HOWEVER, the one in GridapHybrid.jl is conceptually different from the one in Gridap.jl. The former is a cell-wise triangulation, while the latter is a facet-wise triangulation. The former lets you implement integrals over the cell's boundaries, on a cell-wise basis (as per required by static condensation approaches in hybridizable methods), while the later lets you implement integrals over the facets with jumps and averages, as per required by DG schemes.

[JordiManyer]

@amartinhuertas Here are some thoughts:

Let's consider we have two Triangulations, the source strian and the target ttrian (with corresponding glues sglue and tglue). We also consider a CellField defined in the source, i.e sfield.

In vanilla Gridap, we can change domain between the two trians iff sglue.mface_to_tface exists. Then (roughly) what Gridap does is map sfield to the backgroud model using sglue.mface_to_tface and compose it with the inverse coordinate map of the target trian (tglue.tface_to_mface_map).

This means that if strian is a cell-wise triangulation (Dc==Dp) we already have all the tools we need to create the transformation doing something like this:

# Coarse to Fine
v_Ωc = interpolate(sol,Vc)
aux  = Triangulation(ReferenceFE{num_cell_dims(Ωc)},get_adapted_model(Λf))
v_Ωf = change_domain(v_Ωc.plus,aux,ReferenceDomain())
v_Λf = change_domain(v_Ωf,Λf,ReferenceDomain())

# Fine to Coarse
v_Ωf = interpolate(sol,Vf)
aux  = Triangulation(ReferenceFE{num_cell_dims(Ωf)},get_background_model(Λc))
v_Ωc = change_domain(v_Ωf.plus,aux,ReferenceDomain())
v_Λc = change_domain(v_Ωc,Λc,ReferenceDomain())

In this example, Ω is cell-wise and Λ is a skeleton (but would work for any other facet-wise triangulation). I would say this could work for any strian that is cell-wise and any ttrian such that is_change_possible(aux,ttrian) == true.

In the case of model portions, I think we could do something similar by taking the already implemented concept of get_active_model(trian).

So the only case we are not covering is when strian is not cell-wise. In this case, the AdaptivityGlue is not straightforward. The main issue is that, while every fine cell is contained within a coarse cell, that is not the case for lower dimensional facets. Does that imply we can only transfer triangulations which contain facets that are present in both models? Or should we project fine facets into the coarse cell containing it?

[amartinhuertas]

@amartinhuertas Here are some thoughts:

@JordiManyer ... thanks for sharing your thoughts.

I think that the first thing we (at least I) would need in order to move forward is to (try to) fully understand how the change_domain mechanism is currently implemented in Gridap and the underlying rationale. I think I currently don't.

Do you think it might be possible to come up with an internal comprehensive note that fully details, for all possible scenarios (e.g., all possible combinations of strian and ttrian) that are covered, how does the change_domain mechanism works?

I think it would help to write a code snippet where all the cases are covered, and do some debugging to understand the underlying processes and help writing the explanations for each case. Let us start with a very simple case. strian is a BodyFittedTriangulation{D,D} resulting from Triangulation(model) and ttrian is also a BodyFittedTriangulation{D,D} resulting from Base.view(t::BodyFittedTriangulation,ids::AbstractArray).

For example, in the following comment:

In vanilla Gridap, we can change domain between the two trians iff sglue.mface_to_tface exists. Then (roughly) what Gridap does is map sfield to the backgroud model using sglue.mface_to_tface and compose it with the inverse coordinate map of the target trian (tglue.tface_to_mface_map).

1. The reasoning applies whenever sglue and tglue are of type FaceToFaceGlue. But FacetoFaceGlue is not the only type which can be returned by any get_glue method, right?
2. what do you mean by "map sfield to the backgroud model using sglue.mface_to_tface" ?
3. why tglue.tface_to_mface_map is the inverse coordinate map of the target trian? what do you mean?

Does that imply we can only transfer triangulations which contain facets that are present in both models? Or should we project fine facets into the coarse cell containing it?

Before asking ourselves these kind of questions I think we should ask ourselves what do we want/need to do (e.g., by writing high level API snippets of scenarios we will have to deal with). This would lead to the appropriate questions.

The facet triangulations we are interested in, i.e., BoundaryTriangulation and SkeletonTriangulation are glued to the cells in the model (see Face2CellGlue data type; note that this is not a type which is returned by get_glue, despite having Glue in its name). This may solve the issue as we could express everything in terms of cells. But this is just a rough though, I do not have the complete picture in mind.

[JordiManyer]

@amartinhuertas

The reasoning applies whenever sglue and tglue are of type FaceToFaceGlue. But FacetoFaceGlue is not the only type which can be returned by any get_glue method, right?

As far as I can see, we have:

• When calling get_glue on a BodyFittedTriangulation we get the FaceToFaceGlue that the structure contains.
• A BoundaryTriangulation, despite having a FaceToCellGlue, builds and returns a FaceToFaceGlue.
• A SkeletonTriangulation calls get_glue by default on its plus triangulation, which is a BoundaryTriangulation. This also returns a FaceToFaceGlue.

Also, the only definitions of change_domain that I could find doing a search across Gridap, GridapEmbedded and GridapHydrid are defined for FaceToFaceGlue.

[amartinhuertas]

As far as I can see, we have:

Ok, thanks for the update. I see, e.g., the following definition:

function get_glue(trian::SkeletonTriangulation,::Val{D},::Val{D}) where D
  plus = get_glue(trian.plus,Val(D))
  minus = get_glue(trian.minus,Val(D))
  SkeletonPair(plus,minus)
end

which does not actually return a FaceToFaceGlue. In any case, I agree that FaceToFaceGlue seems to be the basic building block data type on which get_glue methods rely.

Given that a single data type covers quite a lot of scenarios, I think it may happen that its three member variables may have different semantics depending on context. Can we come up with an explanation of the semantics of each FaceToFaceGlue instance depending on each possible context? I think this may help a lot to inform our discussions, as I mentioned in #868 (comment)

On the other hand:

When calling get_glue on a BodyFittedTriangulation we get the FaceToFaceGlue that the structure contains.

BodyFittedTriangulation does not actually contain a FaceToFaceGlue instance. (not a big deal, I understand what you mean)

A BoundaryTriangulation, despite having a FaceToCellGlue, builds and returns a FaceToFaceGlue.

Is there any relationship among both glues?

A SkeletonTriangulation calls get_glue by default on its plus triangulation, which is a BoundaryTriangulation. This also returns a FaceToFaceGlue.

Ok ... you refer to the following definition:

get_glue(t::SkeletonTriangulation{D},::Val{D}) where D = get_glue(t.plus,Val(D))

Do you understand the underlying rationale? Which is the meaning of the two dimensions that are passed to the most general get_glue method?

[JordiManyer]

Do you understand the underlying rationale? Which is the meaning of the two dimensions that are passed to the most general get_glue method?

Yeah sorry, that was an over-simplification. In any case, the SkeletonPair returned contains two FaceToFaceGlues, which then get dispatched depending on the plus/minus CellFieldAt. It ends up defaulting to the same implementation of change_domain.

[JordiManyer]

Given that a single data type covers quite a lot of scenarios, I think it may happen that its three member variables may have different semantics depending on context. Can we come up with an explanation of the semantics of each FaceToFaceGlue instance depending on each possible context?

The idea is that when we call get_glue(trian::Triangulation{Dt},Val{Dm}) we return a FaceToFaceGlue that maps the faces in the Triangulation (of dimension Dt) with the faces of the model containing them of dimension Dm. Some examples:

A) BodyFittedTriangulation: It's a cell triangulation, so can only create a glue for Dm==Dp==Dcells. Then

• tface_to_mface contains for each triangulation cell the corresponding cell id in the background model.
• mface_to_tface contains for each model cell the corresponding triangulation cell. In the case of views (where not every model cell belongs to the triangulation) this becomes a PosNegPartition (positive->owned, negative->not owned).
• tface_to_mface_map contains coordinate maps from the triangulation cells to the model cells. In this case, this is always a vector of identity maps.

B) BoundaryTriangulation:
B.1 - In the case of Dm == Dt, the glue contains exactly the same information as described above but for edges (i.e the mappings are between edges in the triangulation and edges in the model). This is then an edge-wise glue. The maps are also identity.
B.2 - In the case Dm == Dcells == Dt+1 we have the following (always plus orientation):

• tface_to_mface contains for each triangulation edge the cell id in the background model that contains it.
• mface_to_tface is nothing. This means we cannot go from a BoundaryTriangulation to a BodyFittedTriangulation (as explained in above discussions).
• tface_to_mface_map contains coordinate maps from the triangulation edges to the model cells. This is no longer identity maps, but rather maps that take the 1D coordinate from the edge reference system and maps it into the 2D coordinates in the corresponding reference cell.

C) The same can be said for SkeletonTriangulations, except that the case where Dm == Dcells == Dt+1 returns a SkeletonPair with the two glues mapping edges in the triangulation into plus and minus cells in the model (following what has been said for BoundaryTriangulations).

D) Views: As already mentioned above, what happens with views is that mface_to_tface is a PosNegPartition. change_domain leverages this by inserting zero fields in model cells with negative indexes (i.e outside of the domain of the triangulation).

[amartinhuertas]

The idea is that when we call get_glue(trian::Triangulation{Dt},Val{Dm}) ...

Be careful here.

There is a more general version of get_glue in which there are to Val-type arguments, not just one, right?

Does it mean that, in the most general scenario, we may end having three different values for D?

get_glue(trian::Triangulation{D1},Val{D2},Val{D3}) where {D1,D2,D3}

I am trying to understand ...

[amartinhuertas]

No, I am quite sure only get_glue(trian,Val{Dm}) gets called in practice, which fetches Dt = num_cell_dims(background_model) and calls get_glue(trian,Dt,Dm). So in practice we only have two dimension parameters (which is the only thing that makes sense).

Ok. So we have to deal with the topological dimension of triangulation faces (Dt) and the topological dimension of model faces (Dm) we want the triangulation faces to be glued to. In the most general case, Dt and Dm values do not match.

[JordiManyer]

I think the bulk of the work is done! I have also added some tests for the cases that work, which are views and BoundaryTriangulations.

SkeletonTriangulations work in theory, but I will need to add some dispatching to correctly treat SkeletonPairs (which are necessary now that I use the triangulation glues). I will do that tomorrow.

I also need to set new @notimplementedif flags for the cases that are not taken care of, which are the cases where both triangulations are of lower dimension (Boundary-Boundary, Skeleton-Skeleton, etc...). I have discussed with @santiagobadia today that these cases are not obvious and that we have no real reason to tackle them at the moment.

[JordiManyer]

@amartinhuertas I believe this is now ready to be reviewed.
We now support all cases where the source triangulation is cell-based.

The case where the source and target trians are face-based (Skeleton-Skeleton, Boundary-Boundary, etc...) is quite complex and I don't think that it can be treated in general. The reason is that some fine edges do not belong to the coarse grid, which poses the question of how to extrapolate the coarse cellfield into those edges when going from coarse to fine fields. Dealing with this could very well be FE method-dependent, i.e by using static condensation for HDG or simple extrapolation in other cases.

In any case, it is probably better to leave it for future work.

[JordiManyer]

The recent commit d6a53e3 is related to this issue.