Append triangulations

src/Arrays/AppendedArrays.jl

@@ -0,0 +1,118 @@
+
+"""
+"""
+function lazy_append(a::AbstractArray,b::AbstractArray)
+  AppendedArray(a,b)
+end
+
+"""
+"""
+function lazy_split(f::AbstractArray,n::Integer)
+  _lazy_split(f,n)
+end
+
+function lazy_split(f::CompressedArray,n::Integer)
+  _a , _b =_lazy_split(f,n)
+  a = _compact_values_ptrs(_a)
+  b = _compact_values_ptrs(_b)
+  a,b
+end
+
+function _lazy_split(f,n)
+  l = length(f)
+  @assert n <= l
+  ids_a = collect(1:n)
+  ids_b = collect((n+1):l)
+  a = reindex(f,ids_a)
+  b = reindex(f,ids_b)
+  a,b
+end
+
+function _compact_values_ptrs(a)
+  v = a.values
+  p = a.ptrs
+  touched = fill(false,length(v))
+  for i in p
+    touched[i] = true
+  end
+  if all(touched) || length(p) == 0
+    return a
+  else
+    j_to_i = findall(touched)
+    i_to_j = zeros(Int,length(v))
+    i_to_j[j_to_i] = 1:length(j_to_i)
+    _v = v[j_to_i]
+    _p = i_to_j[p]
+    return CompressedArray(_v, _p)
+  end
+end
+
+struct AppendedArray{T,A,B} <: AbstractVector{T}
+  a::A
+  b::B
+  function AppendedArray(a::AbstractArray,b::AbstractArray)
+    A = typeof(a)
+    B = typeof(b)
+    ai = testitem(a)
+    bi = testitem(b)
+    T = eltype(collect((ai,bi)))
+    new{T,A,B}(a,b)
+  end
+end
+
+Base.IndexStyle(::Type{<:AppendedArray}) = IndexLinear()
+
+function Base.getindex(v::AppendedArray,i::Integer)
+  l = length(v.a)
+  if i > l
+    v.b[i-l]
+  else
+    v.a[i]
+  end
+end
+
+Base.size(v::AppendedArray) = (length(v.a)+length(v.b),)
+
+function array_cache(v::AppendedArray)
+  ca = array_cache(v.a)
+  cb = array_cache(v.b)
+  (ca,cb)
+end
+
+function getindex!(cache,v::AppendedArray,i::Integer)
+  ca, cb = cache
+  l = length(v.a)
+  if i > l
+    getindex!(cb,v.b,i-l)
+  else
+    getindex!(ca,v.a,i)
+  end
+end
+
+function apply(f,a::AppendedArray...)
+  la = map(ai->length(ai.a),a)
+  lb = map(ai->length(ai.b),a)
+  if all(la .== first(la)) && all(lb .== first(lb))
+    c_a = apply(f,map(ai->ai.a,a)...)
+    c_b = apply(f,map(ai->ai.b,a)...)
+    lazy_append(c_a,c_b)
+  else
+    s = common_size(a...)
+    apply(Fill(f,s...),a...)
+  end
+end
+
+function apply(::Type{T},f,a::AppendedArray...) where T
+  la = map(ai->length(ai.a),a)
+  lb = map(ai->length(ai.b),a)
+  if all(la .== first(la)) && all(lb .== first(lb))
+    c_a = apply(T,f,map(ai->ai.a,a)...)
+    c_b = apply(T,f,map(ai->ai.b,a)...)
+    lazy_append(c_a,c_b)
+  else
+    s = common_size(a...)
+    apply(T,Fill(f,s...),a...)
+  end
+end
+
+

src/Arrays/Arrays.jl

@@ -86,6 +86,10 @@ export unpair_arrays
 
 export matvec_muladd!
 
+export lazy_append
+export lazy_split
+export AppendedArray
+
 import Base: size
 import Base: getindex, setindex!
 import Base: similar
@@ -118,4 +122,6 @@ include("SubVectors.jl")
 
 include("ArrayPairs.jl")
 
+include("AppendedArrays.jl")
+
 end # module

src/Arrays/Interface.jl

@@ -188,8 +188,8 @@ function test_array(
   t = true
   for i in eachindex(a)
     ai = getindex!(cache,a,i)
-    t = t && (typeof(ai) == eltype(a))
-    t = t && (typeof(ai) == T)
+    t = t && (typeof(ai) <: eltype(a))
+    t = t && (typeof(ai) <: T)
   end
   @test t
   @test IndexStyle(a) == IndexStyle(b)

src/FESpaces/ExtendedFESpaces.jl

@@ -77,11 +77,15 @@ function reindex(a::ExtendedVector,ptrs::IdentityVector)
   a
 end
 
-function reindex(a::ExtendedVector,trian::Triangulation)
-  ptrs = get_cell_id(trian)
-  _extended_reindex(a,ptrs)
+function reindex(a::ExtendedVector,ptrs::SkeletonPair)
+  _extended_reindex(a,ptrs::SkeletonPair)
 end
 
+#function reindex(a::ExtendedVector,trian::Triangulation)
+#  ptrs = get_cell_id(trian)
+#  _extended_reindex(a,ptrs)
+#end
+
 function _extended_reindex(a,ptrs::SkeletonPair)
   left= _extended_reindex(a,ptrs.left)
   right = _extended_reindex(a,ptrs.right)

src/Fields/FieldArrays.jl

@@ -284,3 +284,26 @@ for op in (:+,:-)
   end
 end
 
+# More optimizations
+
+function evaluate_field_array(a::AppendedArray,b::AppendedArray)
+  if (length(a.a) == length(b.a)) && (length(a.b) == length(b.b))
+    c_a = evaluate_field_array(a.a,b.a)
+    c_b = evaluate_field_array(a.b,b.b)
+    lazy_append(c_a,c_b)
+  else
+    _evaluate_field_array(a,b)
+  end
+end
+
+function evaluate_field_array(a::AppendedArray,b::AbstractArray)
+  n = length(a.a)
+  _b = lazy_append(lazy_split(b,n)...)
+  evaluate_field_array(a,_b)
+end
+
+function field_array_gradient(a::AppendedArray)
+  c_a = field_array_gradient(a.a)
+  c_b = field_array_gradient(a.b)
+  lazy_append(c_a,c_b)
+end

src/Geometry/AppendedTriangulations.jl

@@ -0,0 +1,79 @@
+
+function lazy_append(a::Triangulation,b::Triangulation)
+  AppendedTriangulation(a,b)
+end
+
+struct AppendedTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
+  a::Triangulation{Dc,Dp}
+  b::Triangulation{Dc,Dp}
+end
+
+function get_cell_coordinates(trian::AppendedTriangulation)
+  a = get_cell_coordinates(trian.a)
+  b = get_cell_coordinates(trian.b)
+  lazy_append(a,b)
+end
+
+function get_reffes(trian::AppendedTriangulation)
+  vcat(get_reffes(trian.a),get_reffes(trian.b))
+end
+
+function get_cell_type(trian::AppendedTriangulation)
+  a = get_cell_type(trian.a)
+  b = get_cell_type(trian.b) .+ Int8(length(get_reffes(trian.a)))
+  lazy_append(a,b)
+end
+
+function reindex(f::AbstractArray,trian::AppendedTriangulation)
+  a = reindex(f,trian.a)
+  b = reindex(f,trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_id(trian::AppendedTriangulation)
+  a = get_cell_id(trian.a)
+  b = get_cell_id(trian.b)
+  lazy_append(a,b)
+end
+
+function restrict(f::AbstractArray,trian::AppendedTriangulation)
+  a = restrict(f,trian.a)
+  b = restrict(f,trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_reffes(trian::AppendedTriangulation)
+  a = get_cell_reffes(trian.a)
+  b = get_cell_reffes(trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_shapefuns(trian::AppendedTriangulation)
+  a = get_cell_shapefuns(trian.a)
+  b = get_cell_shapefuns(trian.b)
+  lazy_append(a,b)
+end
+
+function get_cell_map(trian::AppendedTriangulation)
+  a = get_cell_map(trian.a)
+  b = get_cell_map(trian.b)
+  lazy_append(a,b)
+end
+
+function get_normal_vector(trian::AppendedTriangulation)
+  cm = get_cell_map(trian)
+  a = get_array(get_normal_vector(trian.a))
+  b = get_array(get_normal_vector(trian.b))
+  GenericCellField(lazy_append(a,b),cm)
+end
+
+function CellQuadrature(trian::AppendedTriangulation,degree1::Integer,degree2::Integer)
+  quad1 = CellQuadrature(trian.a,degree1)
+  quad2 = CellQuadrature(trian.b,degree2)
+  lazy_append(quad1,quad2)
+end
+
+function CellQuadrature(trian::AppendedTriangulation,degree::Integer)
+  CellQuadrature(trian,degree,degree)
+end
+

src/Geometry/CartesianGrids.jl

@@ -215,6 +215,11 @@ function field_array_gradient(a::CartesianMap)
   Fill(j,length(a))
 end
 
+function field_array_gradient(a::Reindexed{T,N,A}) where {T,N,A<:CartesianMap}
+  g = field_array_gradient(a.i_to_v)
+  reindex(g,a.j_to_i)
+end
+
 function get_cell_map(grid::CartesianGrid{D,T,typeof(identity)} where {D,T})
   CartesianMap(grid.node_coords.data)
 end

src/Geometry/CellQuadratures.jl

@@ -77,6 +77,12 @@ function _get_coordinates(q::Fill{<:Quadrature})
   Fill(coords,length(q))
 end
 
+function _get_coordinates(q::AppendedArray)
+  a = _get_coordinates(q.a)
+  b = _get_coordinates(q.b)
+  lazy_append(a,b)
+end
+
 function _get_weights(q::AbstractArray{<:Quadrature})
   @notimplemented "Not implemented, since we dont need it"
 end
@@ -91,6 +97,12 @@ function _get_weights(q::Fill{<:Quadrature})
   Fill(w,length(q))
 end
 
+function _get_weights(q::AppendedArray)
+  a = _get_weights(q.a)
+  b = _get_weights(q.b)
+  lazy_append(a,b)
+end
+
 """
     integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
 
@@ -108,3 +120,8 @@ function integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
   integrate(get_array(f),q,w,j)
 end
 
+function lazy_append(quad1::CellQuadrature,quad2::CellQuadrature)
+  array = lazy_append(quad1.array,quad2.array)
+  CellQuadrature(array)
+end
+

src/Geometry/Geometry.jl

@@ -27,10 +27,13 @@ using Gridap.ReferenceFEs: _get_offsets
 using Gridap.ReferenceFEs: _get_offset
 using Gridap.ReferenceFEs: _find_unique_with_indices
 
+using Gridap.Arrays: Reindexed
+
 import Gridap.Arrays: array_cache
 import Gridap.Arrays: getindex!
 import Gridap.Arrays: reindex
 import Gridap.Arrays: get_array
+import Gridap.Arrays: lazy_append
 
 import Gridap.Fields: field_cache
 import Gridap.Fields: evaluate_field!
@@ -202,6 +205,8 @@ export QPointCellField
 
 export RestrictedDiscreteModel
 
+export AppendedTriangulation
+
 include("GridTopologies.jl")
 
 include("GridTopologyMocks.jl")
@@ -252,4 +257,6 @@ include("CellQuadratures.jl")
 
 include("QPointCellFields.jl")
 
+include("AppendedTriangulations.jl")
+
 end # module

src/Geometry/RestrictedTriangulations.jl

@@ -54,3 +54,9 @@ end
 function get_cell_id(trian::RestrictedTriangulation)
   trian.cell_to_oldcell
 end
+
+function get_cell_map(trian::RestrictedTriangulation)
+  cell_map = get_cell_map(trian.oldtrian)
+  reindex(cell_map,trian.cell_to_oldcell)
+end
+

src/Geometry/TriangulationPortions.jl

@@ -23,3 +23,7 @@ function get_cell_coordinates(trian::TriangulationPortion)
   reindex(get_cell_coordinates(trian.oldtrian),trian.cell_to_oldcell)
 end
 
+function get_cell_map(trian::TriangulationPortion)
+  cell_map = get_cell_map(trian.oldtrian)
+  reindex(cell_map,trian.cell_to_oldcell)
+end