Bug fixed restrict(::AbstractArray,::TriangulationPortion) for portions of BoundaryTriangulation(s)

[amartinhuertas]

Hi @fverdugo,

when working with GridapDistributed.jl, I figured out that the (current) implementation of restrict(::AbstractArray,::TriangulationPortion), available here is wrong whenever it is fed up with a portion of a triangulation that extends BoundaryTriangulation (or a portion of a portion of BoundaryTriangulation, a portion of a portion of a portion of BoundaryTriangulation, etc.). The current patch fixes the implementation of restrict in such cases. It also immediately fixes restrict(::AbstractArray,::TriangulationPortion) for portions of SkeletonTriangulations, as these rely on portions of BoundaryTriangulation's. I tested that this fix works with GridapDistributed.jl, but I did not find a way to extend the current tests in Gridap.jl in order to test this feature. Any idea?

@amartinhuertas

[codecov-commenter]

Codecov Report

Merging #267 into master will not change coverage.
The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master     #267   +/-   ##
=======================================
  Coverage   89.00%   89.00%           
=======================================
  Files         147      147           
  Lines        9646     9646           
=======================================
  Hits         8585     8585           
  Misses       1061     1061           

Impacted Files	Coverage Δ
src/Geometry/BoundaryTriangulations.jl	74.07% <ø> (ø)
src/Geometry/TriangulationPortions.jl	83.33% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1fde696...eca1dda. Read the comment docs.

[fverdugo]

@amartinhuertas thanks for reporting! I am sure that this was difficult to spot...

In any case, I think that this simpler fix is enough:

function restrict(f::AbstractArray,trian::TriangulationPortion)
  reindex(restrict(f,trian.oldtrian),trian)
end

Can you check if it works?

[amartinhuertas]

Can you check if it works?

I tested it. It did not work (see output below). Hard to say from this output why it did not work. You can try to reproduce the current error in your machine if u like. I pushed your proposed fix in the branch underlying this PR.

What's the rationale behind your proposed fix? (I do not fully understand it)

include("test/DistributedPoissonDGTests.jl")
WARNING: replacing module DistributedPoissonDGTests.
Test Failed at 
/home/amartin/git-repos/GridapDistributed.jl/test/DistributedPoissonDGTests.jl:106
  Expression: e_l2 < tol
   Evaluated: 0.0002790769102222115 < 1.0e-9
ERROR: LoadError: There was an error during testing
in expression starting at 
/home/amartin/git-repos/GridapDistributed.jl/test/DistributedPoissonDGTests.jl:111

[fverdugo]

Wow this is very strange... the tests in travis are passing.

In any case, I have found a "little" error in the fix I proposed.

@amartinhuertas can you try if the following works ?

function restrict(f::AbstractArray,trian::TriangulationPortion)
  reindex(restrict(f,trian.oldtrian),trian.cell_to_oldcell)
end

It also would be nice to find a test that only works with this last fix (if it ends up to be correct) and not with the previous one.

[fverdugo]

It also would be nice to find a test that only works with this last fix (if it ends up to be correct) and not with the previous one.

Try with a case in which num_cells(trian.oldtrian) is smaller than some id in get_cell_id(trian)

[amartinhuertas]

@amartinhuertas can you try if the following works ?

Great! Now it works.

But I still do not understand the rationale behind your modified proposal.

Can u please elaborate a little bit on it?

Once I understand, then I can proceed with selecting a judicious test.

Thanks!

[fverdugo]

Can u please elaborate a little bit on it?

Sure!

It is important to note that we have 3 different sets of cells:

1. cells in the background interpolation mesh
2. cells in the oldtrian
3. cells in the trian

Let's call bgcell, oldcell, and cell to a generic cell id for each of these sets of cells. Now, let me rewrite the function with this notation:

function restrict(bgcell_to_f::AbstractArray,trian::TriangulationPortion)
  oldcell_to_f = restrict(bgcell_to_f,trian.oldtrian)
  cell_to_f = reindex(oldcell_to_f,trian.cell_to_oldcell)
  return cell_to_f
end

Written in this way, the code is self reveling.

Another important thing to understand:

oldcell_to_bgcell = get_cell_id(oldtrian) 
cell_to_bgcell = get_cell_id(trian) 

That is, indices in get_cell_id() always point to cells in the interpolation mesh

[fverdugo]

The tests are passing. However, I would wait to have the new test we were commenting so that we do not have a regression in this part.

[amartinhuertas]

Sure!
It is important to note that we have 3 different sets of cells:
...

Thanks for the explanation! Now I got the point. Your modification is clearly the way to go. Let me complement your explanation with two extra points that I think are true and helped me to understand the whole picture (can you please explicitly confirm/check that they are true?):

1. Underlying the restrict(f::AbstractArray,t::Triangulation) generic function is the precondition that f has to be of the same length as the number of cells of dimension Dc in the DiscreteModel t stems from (either directly or indirectly). If you agree, perhaps we could add this statement as documentation for restrict, and/or as a precondition to be checked via @assert.

2. In the current (final) implementation of restrict(f::AbstractArray,t::TriangulationPortion), the recursion which is triggered with restrict(f,t.oldtrian) ends up in two possible scenarios (are there more?):
   2.1. t.oldtrian is the triangulation that directly stemmed from DiscreteModel, and of the same Dc as it. In such as case, restrict(f,t.oldtrian)=f.
   2.2. t.oldtrian is a triangulation that extends BoundaryTriangulation. In such a case, restrict(f,t.oldtrian)=compose_field_arrays(reindex(f,t.oldtrian),get_face_to_cell_map(t.oldtrian))

@amartinhuertas

[amartinhuertas]

The tests are passing. However, I would wait to have the new test we were commenting so that we do not have a regression in this part.

Done in 829bdb0!

However, I have been trying to design a more clever test, without success. In particular, a test that spots the BUG that gives name to this PR: TriangulationPortion of BoundaryTriangulation. See the code I wrote so far below. The BUG was that the type of v_btrian_portion.array was wrong. In particular, it ressembled more v.array, without being totally equal (because of a ReindexedArray in the middle of the awkarly large data type name). In other words, the test @test typeof(v.array) != typeof(v_btrian_portion.array) also passed with the BUG. How could we properly design the test?

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
btrian = BoundaryTriangulation(model,[true for i=1:num_faces(model,1)])
btrian_portion = TriangulationPortion(btrian,[2,4,6])
btrian_portion_portion = TriangulationPortion(btrian_portion,[3,1])
order = 1
V0 = TestFESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="boundary")
v=Gridap.FESpaces.get_cell_basis(V0)
v_btrian_portion=restrict(v,btrian_portion)
v_btrian_portion_portion=restrict(v,btrian_portion_portion)
@test length(v_btrian_portion)==3
@test length(v_btrian_portion_portion)==2
@test typeof(v.array) != typeof(v_btrian_portion.array)

[fverdugo]

@amartinhuertas can you solve the conflict?

[fverdugo]

Underlying the restrict(f::AbstractArray,t::Triangulation) generic function is the precondition that f has to be of the same length as the number of cells of dimension Dc in the DiscreteModel t stems from (either directly or indirectly). If you agree, perhaps we could add this statement as documentation for restrict, and/or as a precondition to be checked via @Assert.

Correct. Be also aware that Dc in the discrete model can be smaller than Dp in general.

In the current (final) implementation of restrict(f::AbstractArray,t::TriangulationPortion), the recursion which is triggered with restrict(f,t.oldtrian) ends up in two possible scenarios (are there more?):
2.1. t.oldtrian is the triangulation that directly stemmed from DiscreteModel, and of the same Dc as it. In such as case, restrict(f,t.oldtrian)=f.
2.2. t.oldtrian is a triangulation that extends BoundaryTriangulation. In such a case, restrict(f,t.oldtrian)=compose_field_arrays(reindex(f,t.oldtrian),get_face_to_cell_map(t.oldtrian))

I think that in Gridap these are the two cases, however the number can increase since the code is extensible. E.g., in GridapEmbedded there are other cases.

[amartinhuertas]

@amartinhuertas can you solve the conflict?

Yes, done.

[fverdugo]

However, I have been trying to design a more clever test, without success. In particular, a test that spots the BUG that gives name to this PR: TriangulationPortion of BoundaryTriangulation. See the code I wrote so far below. The BUG was that the type of v_btrian_portion.array was wrong. In particular, it ressembled more v.array, without being totally equal (because of a ReindexedArray in the middle of the awkarly large data type name). In other words, the test @test typeof(v.array) != typeof(v_btrian_portion.array) also passed with the BUG. How could we properly design the test?

I think I will merge the PR even though this is not fixed since we need these devs in GridapDistributed ASAP. We can keep track of the pending test in an issue.

[amartinhuertas]

I think I will merge the PR even though this is not fixed since we need these devs in GridapDistributed ASAP. We can keep track of the pending test in an issue.

Ok. Agreed. Can u create the issue?

[amartinhuertas]

I think that in Gridap these are the two cases, however the number can increase since the code is extensible. E.g., in GridapEmbedded there are other cases.

May the BUG also affect to GridapEmbedded? Have you think about that?