Pitfall with `change_domain` and `AdaptiveTriangulation`s

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 `AdaptedTriangulation`s:

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

```julia
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:

```julia
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 `AdaptiveTriangulation`s 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?