get_cell_map returns array of AffineMap for linear grids of simplices

src/CellData/CellFields.jl

@@ -79,7 +79,7 @@ end
 
 function CellField(f::Number,trian::Triangulation,domain_style::DomainStyle)
   s = size(get_cell_map(trian))
-  cell_field = Fill(ConstantField(f),s)
+  cell_field = Fill(constant_field(f),s)
   GenericCellField(cell_field,trian,domain_style)
 end
 
@@ -89,7 +89,7 @@ function CellField(f::AbstractArray{<:Number},trian::Triangulation,domain_style:
   on a Triangulation with $(num_cells(trian)) cells. The length of the given array
   and the number of cells should match.
   """
-  cell_field = lazy_map(ConstantField,f)
+  cell_field = lazy_map(constant_field,f)
   GenericCellField(cell_field,trian,domain_style)
 end
 

src/Fields/AffineMaps.jl

@@ -3,14 +3,19 @@
 A Field with this form
 y = x⋅G + y0
 """
-struct AffineMap{D,T,L} <:Field
-  gradient::TensorValue{D,D,T,L}
-  origin::Point{D,T}
-  function AffineMap(gradient::TensorValue{D,D,T,L}, origin::Point{D,T}) where {D,T,L}
-    new{D,T,L}(gradient,origin)
+struct AffineMap{D1,D2,T,L} <:Field
+  gradient::TensorValue{D1,D2,T,L}
+  origin::Point{D2,T}
+  function AffineMap(
+    gradient::TensorValue{D1,D2,T,L},
+    origin::Point{D2,T}) where {D1,D2,T,L}
+
+    new{D1,D2,T,L}(gradient,origin)
   end
 end
 
+@inline affine_map(gradient,origin) = AffineMap(gradient,origin)
+
 @inline function evaluate!(cache,f::AffineMap,x::Point)
   G = f.gradient
   y0 = f.origin
@@ -77,3 +82,9 @@ function inverse_map(f::AffineMap)
   x0 = -y0⋅invJt
   AffineMap(invJt,x0)
 end
+
+function lazy_map(::typeof(∇),a::LazyArray{<:Fill{typeof(affine_map)}})
+  gradients = a.args[1]
+  lazy_map(constant_field,gradients)
+end
+

src/Fields/Fields.jl

@@ -48,6 +48,7 @@ export return_cache
 export Field
 export GenericField
 export ConstantField
+export constant_field
 export FieldGradient
 export FieldGradientArray
 export ZeroField
@@ -57,6 +58,7 @@ export Point
 export inverse_map
 
 export AffineMap
+export affine_map
 
 export gradient
 export ∇

src/Fields/FieldsInterfaces.jl

@@ -228,6 +228,8 @@ struct ConstantField{T<:Number} <: Field
   object::T
 end
 
+@inline constant_field(a) = ConstantField(a)
+
 Base.zero(::Type{ConstantField{T}}) where T = ConstantField(zero(T))
 
 @inline function evaluate!(c,f::ConstantField,x::Point)
@@ -270,6 +272,12 @@ function evaluate!(c,f::FieldGradient{N,<:ConstantField},x::AbstractArray{<:Poin
   c.array
 end
 
+function lazy_map(::Operation{typeof(inv)},a::LazyArray{<:Fill{typeof(constant_field)}})
+  v = a.args[1]
+  vinv = lazy_map(inv,v)
+  lazy_map(constant_field,vinv)
+end
+
 ## Make Function behave like Field
 
 return_cache(f::FieldGradient{N,<:Function},x::Point) where N = gradient(f.object,Val(N))

src/Geometry/BoundaryTriangulations.jl

@@ -212,7 +212,7 @@ function get_facet_normal(trian::BoundaryTriangulation)
   end
   ctype_lface_pindex_to_nref = map(f, get_reffes(cell_trian))
   face_to_nref = FaceCompressedVector(ctype_lface_pindex_to_nref,glue)
-  face_s_nref = lazy_map(ConstantField,face_to_nref)
+  face_s_nref = lazy_map(constant_field,face_to_nref)
 
   # Inverse of the Jacobian transpose
   cell_q_x = get_cell_map(cell_trian)

src/Geometry/CartesianGrids.jl

@@ -237,7 +237,7 @@ end
 
 # Cell map
 
-struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,T,L},D}
+struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,D,T,L},D}
   data::CartesianDescriptor{D,T,typeof(identity)}
   function CartesianMap(des::CartesianDescriptor{D,T}) where {D,T}
     L = D*D

src/Geometry/UnstructuredGrids.jl

@@ -48,9 +48,22 @@ function _compute_cell_map(node_coords,cell_node_ids,ctype_reffe,cell_ctype)
   cell_coords = lazy_map(Broadcasting(Reindex(node_coords)),cell_node_ids)
   ctype_shapefuns = map(get_shapefuns,ctype_reffe)
   cell_shapefuns = expand_cell_data(ctype_shapefuns,cell_ctype)
-  cell_map = lazy_map(linear_combination,cell_coords,cell_shapefuns)
+  default_cell_map = lazy_map(linear_combination,cell_coords,cell_shapefuns)
+  ctype_poly = map(get_polytope,ctype_reffe)
+  has_affine_map =
+    all(map(is_first_order,ctype_reffe)) &&
+    ( all(map(is_simplex,ctype_poly)) || all(map(p->num_dims(p)==1,ctype_poly)) )
+  if has_affine_map
+    ctype_q0 = map(p->zero(first(get_vertex_coordinates(p))),ctype_poly)
+    cell_q0 = expand_cell_data(ctype_q0,cell_ctype)
+    default_cell_grad = lazy_map(∇,default_cell_map)
+    origins = lazy_map(evaluate,default_cell_map,cell_q0)
+    gradients = lazy_map(evaluate,default_cell_grad,cell_q0)
+    cell_map = lazy_map(Fields.affine_map,gradients,origins)
+  else
+    cell_map = default_cell_map
+  end
   Fields.MemoArray(cell_map)
-  #cell_map
 end
 
 """

src/TensorValues/Operations.jl

@@ -133,12 +133,17 @@ dot(a::MultiValue,b::MultiValue) = @notimplemented
         bk = data_index(B,i,j)
         s *= "a.data[$ak]*b.data[$bk]+"
       end
-        push!(ss,s[1:(end-1)]*", ")
+      push!(ss,s[1:(end-1)]*", ")
     end
     str = join(ss)
     Meta.parse("VectorValue{$D2}($str)")
 end
 
+function dot(a::A,b::B) where {A<:MultiValue{Tuple{0}},B<:MultiValue{Tuple{0,D2}}} where D2
+  T = eltype(zero(eltype(a))*zero(eltype(b)))
+  zero(VectorValue{D2,T})
+end
+
 @generated function dot(a::A,b::B) where {A<:MultiValue{Tuple{D1,D2}},B<:MultiValue{Tuple{D2}}} where {D1,D2}
     ss = String[]
     for i in 1:D1

test/TensorValuesTests/OperationsTests.jl

@@ -251,6 +251,10 @@ c = t3 ⋅ t1
 r = ThirdOrderTensorValue{2,2,1}(1,2,3,6)
 @test c == r
 
+x = VectorValue{0,Float64}()
+G = TensorValue{0,2,Float64,0}()
+@test x⋅G == VectorValue(0,0)
+
 # Inner product (full contraction)
 
 c = 2 ⊙ 3