Space time by extrusion

NEWS.md

@@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [0.17.13] - 2022-05-31
 
 ### Added
+- `KeyToValMap` lazy map that dinamically creates a `Dict` with the outputs of a function over an array of inputs. Since PR [#801](https://github.com/gridap/Gridap.jl/pull/801)
+- `MappedDiscreteModel` and `MappedGrid`, which are geometrical models with one extra geo map in the physical space. Since PR [#801](https://github.com/gridap/Gridap.jl/pull/801)
+- `GridWithFEMap`, which has a geometrical map defined by a FE function. Since PR [#801](https://github.com/gridap/Gridap.jl/pull/801)
 - Vertex bisection algorithm for refinement of triangular meshes in 2D. Since PR [#733](https://github.com/gridap/Gridap.jl/pull/733)
 - Generalized-α method for 1st order ODEs. Since PR [#781](https://github.com/gridap/Gridap.jl/pull/781) 
 - Implemented (generalised) ModalC0 Polynomial bases and reference FEs. Since PR [#777](https://github.com/gridap/Gridap.jl/pull/777)

src/Arrays/Arrays.jl

@@ -78,6 +78,8 @@ export PosNegPartition
 
 export FilterMap
 
+export KeyToValMap
+
 export MulAddMap
 
 export Table
@@ -150,6 +152,8 @@ include("PosNegReindex.jl")
 
 include("Reindex.jl")
 
+include("KeyToValMaps.jl")
+
 include("FilteredArrays.jl")
 
 include("SubVectors.jl")

src/Arrays/KeyToValMaps.jl

@@ -0,0 +1,19 @@
+struct KeyToValMap{T<:Function} <: Map
+  key_to_val::T
+end
+
+function return_cache(m::KeyToValMap,key)
+  K = typeof(key)
+  V = typeof(m.key_to_val(key))
+  d = Dict{K,V}()
+end
+
+function evaluate!(cache,m::KeyToValMap,key)
+  if haskey(cache,key)
+    return get(cache,key,nothing)
+  else
+    val = m.key_to_val(key)
+    cache[key] = val
+    return val
+  end
+end

src/FESpaces/AffineFEOperators.jl

@@ -2,7 +2,7 @@
 """
     AffineFEOperator
 
-Reprepresent a fully assembled affine (linear) finite element problem.
+Represent a fully assembled affine (linear) finite element problem.
 See also [FEOperator](@ref)
 """
 struct AffineFEOperator <: FEOperator

src/FESpaces/DiscreteModelWithFEMaps.jl

@@ -0,0 +1,149 @@
+"""
+Given a Discrete Model and a reffe, builds a new grid in which the geometrical
+map is a `FEFunction`. This is useful when considering geometrical maps that are
+the result of a FE problem (mesh displacement).
+"""
+struct GridWithFEMap{Dc,Dp,A,B,C,D,E} <: Grid{Dc,Dp}
+  grid::Grid{Dc,Dp}
+  fe_sp::A
+  scal_fe_sp::B
+  fe_map::C
+  node_coords::D
+  reffes::E
+end
+
+function GridWithFEMap(model,order; kwargs...)
+  orders = Fill(order,num_cells(model))
+  # _reffes = [reffe]
+  # reffes = lazy_map(Reindex(_reffes),cell_type)
+  GridWithFEMap(model,orders; kwargs...)
+end
+
+function GridWithFEMap(model,orders::AbstractArray; kwargs...)
+
+  # Create a FESpace for the geometrical description
+  # The FEFunction that describes the coordinate field
+  # or displacement can be a interpolation or a solution
+  # of a mesh displacement problem
+  T = eltype(get_node_coordinates(model))
+  Ts = eltype(T)
+
+  os = orders
+  _ps = get_polytopes(model)
+  ct = get_cell_type(model)
+  ps = lazy_map(Reindex(_ps), ct)
+
+  f(a,b) = LagrangianRefFE(T,a,b)
+  reffes = lazy_map(f,ps,os)
+  Vₕ = FESpace(model,reffes;conformity=:H1,kwargs...)
+
+  fs(a,b) = LagrangianRefFE(Ts,a,b)
+  s_reffes = lazy_map(f,ps,os)
+  Vₕ_scal = FESpace(model,s_reffes;conformity=:H1)
+
+  grid = get_grid(model)
+  geo_map = get_cell_map(grid)
+
+  cell_ctype = get_cell_type(grid)
+  c_reffes = get_reffes(grid)
+
+  # Create a fe_map using the cell_map that can be evaluated at the
+  # vertices of the fe space (nodal type)
+  # This returns a FEFunction initialised with the coordinates
+  # But this is to build the FEFunction that will be inserted, it is
+  # an advanced constructor, not needed at this stage
+
+  c_dofs = get_fe_dof_basis(Vₕ)
+  dof_basis = get_data(c_dofs)
+  c_nodes = lazy_map(get_nodes,get_data(c_dofs))
+
+  xh = zero(Vₕ)
+  c_dofv = lazy_map(evaluate,dof_basis,geo_map)
+
+  Uₕ = TrialFESpace(Vₕ)
+
+  fv = get_free_dof_values(xh)
+  dv = get_dirichlet_dof_values(get_fe_space(xh))
+  gather_free_and_dirichlet_values!(fv,dv,Uₕ,c_dofv)
+
+  c_xh = lazy_map(evaluate,get_data(xh),c_nodes)
+  c_scal_ids = get_cell_dof_ids(Vₕ_scal)
+
+
+  nodes_coords = Vector{eltype(eltype(c_xh))}(undef,num_free_dofs(Vₕ_scal))
+  Geometry._cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
+
+  GridWithFEMap(grid, Vₕ, Vₕ_scal, xh, nodes_coords, reffes)
+end
+
+function _compute_node_coordinates(grid,xh)
+  c_dofs = get_fe_dof_basis(grid.fe_sp)
+  c_nodes = lazy_map(get_nodes,get_data(c_dofs))
+
+  c_xh = lazy_map(evaluate,get_data(xh),c_nodes)
+  c_scal_ids = get_cell_dof_ids(grid.scal_fe_sp)
+
+  nodes_coords = grid.node_coords
+  Geometry._cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
+
+  return nothing
+
+end
+
+function add_mesh_displacement!(grid::GridWithFEMap,dh::FEFunction)
+
+  Xh = grid.fe_map
+
+  Xh_fv = get_free_dof_values(Xh)
+  Xh_dv = get_dirichlet_dof_values(get_fe_space(Xh))
+
+  Xh_fv .= Xh_fv .+ get_free_dof_values(dh)
+  Xh_dv .= Xh_dv .+ get_dirichlet_dof_values(get_fe_space(dh))
+
+  _compute_node_coordinates(grid,Xh)
+
+end
+
+function update_coordinates!(grid::GridWithFEMap,dh::FEFunction)
+
+  Xh = grid.fe_map
+
+  Xh_fv = get_free_dof_values(Xh)
+  Xh_dv = get_dirichlet_dof_values(get_fe_space(Xh))
+
+  Xh_fv .= get_free_dof_values(dh)
+  Xh_dv .= get_dirichlet_dof_values(get_fe_space(dh))
+
+  _compute_node_coordinates(grid,Xh)
+
+end
+
+# Interface
+get_node_coordinates(grid::GridWithFEMap) = grid.node_coords
+get_cell_node_ids(grid::GridWithFEMap) = get_cell_dof_ids(grid.scal_fe_sp)
+get_reffes(grid::GridWithFEMap) = grid.reffes
+get_cell_type(grid::GridWithFEMap) = get_cell_type(grid.grid)
+OrientationStyle(grid::GridWithFEMap) = OrientationStyle(grid.grid)
+RegularityStyle(grid::GridWithFEMap) = RegularityStyle(grid.grid)
+# santiagobadia: To think (does it make sense here, I think it is for surface problems)
+get_facet_normal(grid::GridWithFEMap) = get_facet_normal(grid.grid)
+
+
+# struct DiscreteModelWithFEMap{Dc,Dp} <: DiscreteModel{Dc,Dp}
+#   model::DiscreteModel{Dc,Dp}
+#   mapped_grid::GridWithFEMap{Dc,Dp}
+#   function DiscreteModelWithFEMap(model::DiscreteModel{Dc,Dp},phys_map) where {Dc,Dp}
+#     mapped_grid = GridWithFEMap(get_grid(model),phys_map)
+#     new{Dc,Dp}(model,mapped_grid)
+#   end
+# end
+
+function DiscreteModelWithFEMap(model::DiscreteModel,order; kwargs...)
+  mapped_grid = GridWithFEMap(model,order;kwargs...)
+  MappedDiscreteModel(model,mapped_grid)
+end
+
+function DiscreteModelWithFEMap(model,orders::AbstractArray; kwargs...)
+  mapped_grid = GridWithFEMap(model,orders;kwargs...)
+  MappedDiscreteModel(model,mapped_grid)
+end

src/FESpaces/FESpaces.jl

@@ -31,6 +31,15 @@ import Gridap.Arrays: return_cache
 import Gridap.Arrays: evaluate!
 import Gridap.Geometry: get_triangulation
 import Gridap.Geometry: get_cell_shapefuns
+import Gridap.Geometry: get_cell_type
+import Gridap.Geometry: MappedGrid
+import Gridap.Geometry: MappedDiscreteModel
+import Gridap.Geometry: get_node_coordinates
+import Gridap.Geometry: get_reffes
+import Gridap.Geometry: get_cell_node_ids
+import Gridap.Geometry: OrientationStyle
+import Gridap.Geometry: RegularityStyle
+import Gridap.Geometry: get_facet_normal
 import Gridap.CellData: attach_constraints_rows
 import Gridap.CellData: attach_constraints_cols
 import Gridap.CellData: CellField
@@ -184,6 +193,11 @@ export FESpaceWithLinearConstraints
 
 export FiniteElements
 
+export DiscreteModelWithFEMap
+export GridWithFEMap
+export add_mesh_displacement!
+export update_coordinates!
+
 include("FESpaceInterface.jl")
 
 include("SingleFieldFESpaces.jl")
@@ -230,6 +244,8 @@ include("DirichletFESpaces.jl")
 
 include("FESpacesWithLinearConstraints.jl")
 
+include("DiscreteModelWithFEMaps.jl")
+
 export get_free_values
 function get_free_values(args...)
   @unreachable "get_free_values has been removed. Use get_free_dof_values instead."

src/Geometry/Geometry.jl

@@ -124,6 +124,7 @@ export compute_reference_grid
 
 export GridPortion
 export UnstructuredGrid
+export MappedGrid
 
 export CartesianGrid
 export CartesianDescriptor
@@ -164,6 +165,7 @@ export CartesianDiscreteModel
 export GenericTriangulation
 export BoundaryTriangulation
 export DiscreteModelPortion
+export MappedDiscreteModel
 
 export Interior
 export Boundary
@@ -214,6 +216,8 @@ include("UnstructuredDiscreteModels.jl")
 
 include("CartesianDiscreteModels.jl")
 
+include("MappedDiscreteModels.jl")
+
 include("GridPortions.jl")
 
 include("DiscreteModelPortions.jl")

src/Geometry/MappedDiscreteModels.jl

@@ -0,0 +1,81 @@
+function _cell_vector_to_dof_vector!(dof_vector,cell_node_ids, cell_vector)
+  cache_cell_node_ids = array_cache(cell_node_ids)
+  cache_cell_vector   = array_cache(cell_vector)
+  for k=1:length(cell_node_ids)
+    current_node_ids = getindex!(cache_cell_node_ids,cell_node_ids,k)
+    current_values   = getindex!(cache_cell_vector,cell_vector,k)
+    for (i,id) in enumerate(current_node_ids)
+      dof_vector[current_node_ids[i]]=current_values[i]
+    end
+  end
+end
+
+# """
+#     MappedGrid
+
+# Represent a grid with a geometrical map that is the composition of
+# a reference to a physical space map (standard)
+# and a (vector-valued) map from physical space to physical space. E.g.,
+# it can be used to include a high order map implemented by any map that is
+# a `CellField`.
+# """
+struct MappedGrid{Dc,Dp,T,M,L} <: Grid{Dc,Dp}
+  grid::Grid{Dc,Dp}
+  geo_map::T  # Composition of old map and new one
+  phys_map::M # New map in the physical space
+  node_coords::L
+  function MappedGrid(grid::Grid{Dc,Dp},phys_map) where {Dc,Dp}
+
+    @assert length(phys_map) == num_cells(grid)
+    @assert eltype(phys_map) <: Field
+
+    function _compute_node_coordinates(grid,phys_map)
+      cell_node_ids = get_cell_node_ids(grid)
+      old_nodes = get_node_coordinates(grid)
+      node_coordinates = Vector{eltype(old_nodes)}(undef,length(old_nodes))
+      c_coor = get_cell_coordinates(grid)
+      map_c_coor = lazy_map(evaluate,phys_map,c_coor)
+      _cell_vector_to_dof_vector!(node_coordinates,cell_node_ids,map_c_coor)
+      return node_coordinates
+    end
+
+    model_map=get_cell_map(grid)
+    geo_map=lazy_map(∘,phys_map,model_map)
+    node_coords = collect(_compute_node_coordinates(grid,phys_map))
+    new{Dc,Dp,typeof(geo_map),typeof(phys_map),typeof(node_coords)}(grid,geo_map,phys_map,node_coords)
+  end
+end
+
+function MappedGrid(grid::Grid{Dc,Dp},phys_map::Function) where {Dc,Dp}
+  c_map = Fill(GenericField(phys_map),num_cells(grid))
+  MappedGrid(grid,c_map)
+end
+
+get_node_coordinates(grid::MappedGrid) = grid.node_coords
+get_cell_node_ids(grid::MappedGrid) = get_cell_node_ids(grid.grid)
+get_reffes(grid::MappedGrid) = get_reffes(grid.grid)
+get_cell_type(grid::MappedGrid) = get_cell_type(grid.grid)
+
+"""
+MappedDiscreteModel
+
+Represent a model with a `MappedGrid` grid.
+See also [`MappedGrid`](@ref).
+"""
+struct MappedDiscreteModel{Dc,Dp} <: DiscreteModel{Dc,Dp}
+  model::DiscreteModel{Dc,Dp}
+  mapped_grid::MappedGrid{Dc,Dp}
+  function MappedDiscreteModel(model::DiscreteModel{Dc,Dp},phys_map) where {Dc,Dp}
+    mapped_grid = MappedGrid(get_grid(model),phys_map)
+    new{Dc,Dp}(model,mapped_grid)
+  end
+end
+
+get_grid(model::MappedDiscreteModel) = model.mapped_grid
+get_cell_map(model::MappedDiscreteModel) = get_cell_map(model.mapped_grid)
+get_grid_topology(model::MappedDiscreteModel) = get_grid_topology(model.model)
+get_face_labeling(model::MappedDiscreteModel) = get_face_labeling(model.model)
+
+function Grid(::Type{ReferenceFE{d}},model::MappedDiscreteModel) where {d}
+  get_grid(model)
+end

test/ArraysTests/KeyToValMapsTests.jl

@@ -0,0 +1,48 @@
+module KeyToValMapsTests
+
+using Test
+using Gridap.Arrays
+using Gridap.ReferenceFEs
+
+using Gridap
+using FillArrays
+
+keys = [1,2,3,4,5]
+key_to_val(i) = i*ones(10)
+
+m = KeyToValMap(key_to_val)
+r = lazy_map(m,keys)
+
+r[1]
+c = return_cache(r)
+
+_r = key_to_val.(keys)
+@test all(_r .== collect(r))
+
+
+key_to_val(i) = LagrangianRefFE(Float64,HEX,i)
+
+m = KeyToValMap(key_to_val)
+keys = [1,2,3,4,5,4,3,2,1]
+r = lazy_map(m,keys)
+_r = key_to_val.(keys)
+@test all(get_order.(_r)-get_order.(r) .== 0)
+
+end #module
+
+# @time key_to_val(1)
+# @time for i in r end
+# @time for i in keys key_to_val(i) end
+
+# @time key_to_val(5)
+# dict = Dict(keys .=> key_to_val.(keys))
+# # fetching from dict is fast
+# @time dict[5]
+# @time r[5]
+# a = getkey(dict,5,1)
+# a
+# return_cache(m,keys)
+
+# haskey(dict,5)
+# get(dict,5,nothing)
+# get!(dict,5,nothing)