Pitfall with change_domain and AdaptiveTriangulations

[amartinhuertas]

Hi @JordiManyer,

when I run the following MWE:

using Gridap
using Gridap.Adaptivity

modelH=CartesianDiscreteModel((0,1,0,1),(1,1))
modelh=Gridap.Adaptivity.refine(modelH)

f(x)=x[1]+x[2]
cf=CellField(f,Triangulation(modelh))

cf_ref1=change_domain(cf,get_triangulation(cf),ReferenceDomain())
print_op_tree(cf_ref1.cell_field)

cf_ref2=change_domain(cf,ReferenceDomain())
print_op_tree(cf_ref2.cell_field)

I obtain:

LazyArray
├─ Fill
│  └─ Broadcasting{typeof(∘)}
├─ LazyArray
│  ├─ Fill
│  │  └─ Broadcasting{typeof(∘)}
│  ├─ LazyArray
│  │  ├─ Fill
│  │  │  └─ Broadcasting{typeof(∘)}
│  │  ├─ Fill
│  │  │  └─ Gridap.Fields.GenericField{typeof(f)}
│  │  └─ Gridap.Geometry.CartesianMap{2, Float64, 4}
│  └─ LazyArray
│     ├─ Fill
│     │  └─ typeof(Gridap.Arrays.inverse_map)
│     └─ Gridap.Geometry.CartesianMap{2, Float64, 4}
└─ Gridap.Geometry.CartesianMap{2, Float64, 4}

in the first call to print_op_tree, while

LazyArray
├─ Fill
│  └─ Broadcasting{typeof(∘)}
├─ Fill
│  └─ Gridap.Fields.GenericField{typeof(f)}
└─ Gridap.Geometry.CartesianMap{2, Float64, 4}

in the second.

You will see the sequence composition of the direct, inverse, and direct geometrical mappings. This is inefficient, as only a single composition with the direct mapping is enough, as you can see in the second case.

The first call to change_domain is triggered when you, e.g., substract cf and a FE function uh. Namely, from the _to_common_domain(a::CellField...) Gridap.jl function ... so it is a quite frequent case.

In regards to why it happens, it is because the way change_domain is defined for AdaptedTriangulations:

function change_domain(a::CellField,target_trian::Triangulation,target_domain::DomainStyle)
  change_domain(a,get_triangulation(a),DomainStyle(a),target_trian,target_domain)
end

for sdomain in [:ReferenceDomain,:PhysicalDomain]
  for (stype,ttype) in [(:AdaptedTriangulation,:AdaptedTriangulation),(:AdaptedTriangulation,:Triangulation),(:Triangulation,:AdaptedTriangulation)]
    @eval begin
      function CellData.change_domain(a::CellField,strian::$stype,::$sdomain,ttrian::$ttype,::PhysicalDomain)
        a_ref  = change_domain(a,ReferenceDomain())
        atrian = change_domain(a_ref,strian,ReferenceDomain(),ttrian,ReferenceDomain())
        return change_domain(atrian,PhysicalDomain())
      end
    end
  end
end

An immediate way around this, it to modify the definition of the first change_domain right above by:

function change_domain(a::CellField,target_trian::Triangulation,target_domain::DomainStyle)

strian=get_triangulation(a) 
if (strian===target_trian)
change_domain(a,DomainStyle(a),target_domain)
else
change_domain(a,get_triangulation(a),DomainStyle(a),target_trian,target_domain)
end

BUT ... is there something else I might be missing??? First, is the current implementation of change_domain for AdaptiveTriangulations reasonable, if yes why? Second, can it be that we are somehow assuming that the target and source triangulations have to be different for this method to be the way to go?

[JordiManyer]

Hi @amartinhuertas , I agree it is a problem. I think I can explain the reasoning behind the current implementation.

All of this comes from the assumption that we can only deal with reference coordinates when changing between adapted meshes.

When the target domain is the reference domain, I can simply do the change. I believe the functions that do that also check if the triangulations are the same and optimizes accordingly.

When the target domain is the physical domain, we have to do the change in the reference space then convert to physical coordinates. If the input cellfield is in physical domain, that indeed involves all the steps you mention.
From what I see, the function that does this (second one you posted) does not check if the triangulations are the same.

When I come back Ill have a look at this and see if we can make the whole implementation more compact/coherent. But I do think you have a point.

[amartinhuertas]

When I come back Ill have a look at this and see if we can make the whole implementation more compact/coherent. But I do think you have a point.

No hurries. Enjoy your holidays. We talk when you are back.

[amartinhuertas]

Related to PR #886

[JordiManyer]

PR #886 solves the issue discussed here, but we believe there is still some room for improvement. More precisely, we believe some edge/weird cases can still be optimized (for instance when changing between TriangulationViews). We will look at it in the future.