Space time by extrusion

src/Arrays/Arrays.jl

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

src/Arrays/KeyToVal.jl

@@ -0,0 +1,21 @@
+struct KeyToVal{T} <: Map
+  key_to_val::T
+end
+
+
+function return_cache(m::KeyToVal,keys)
+  K = eltype(keys)
+  length(keys) > 0 ? _key = keys[1] : _key = testvalue(K)
+  V = typeof(m.key_to_val(_key))
+  d = Dict{K,V}()
+end
+
+function evaluate!(cache,m::KeyToVal,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/DiscreteModelWithFEMaps.jl

@@ -0,0 +1,156 @@
+"""
+Given a Discrete Model and a reffe, builds a new grid in which the geometrical
+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,F} <: Grid{Dc,Dp}
+  grid::Grid{Dc,Dp}
+  fe_sp::A
+  scal_fe_sp::B
+  fe_map::C
+  node_coords::D
+  reffes::E
+  cell_type::F
+end
+
+function GridWithFEMap(model,reffe::ReferenceFE; kwargs...)
+  cell_type = Fill(1,num_cells(model))
+  _reffes = [reffe]
+  reffes = lazy_map(Reindex(_reffes),cell_type)
+  GridWithFEMap(model,cell_type,reffes; kwargs...)
+end
+
+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
+
+function GridWithFEMap(model,cell_type,reffes::AbstractArray{<:ReferenceFE}; 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))
+
+  Vₕ = FESpace(model,reffes;conformity=:H1,kwargs...)
+
+  Ts = eltype(T)
+  os = get_order.(reffes)
+  ps = get_polytope.(reffes)
+  f(a,b) = LagrangianRefFE(Float64,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ₕ)
+  nc = get_node_coordinates(model)
+  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))
+  _cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
+
+  GridWithFEMap(grid, Vₕ, Vₕ_scal, xh, nodes_coords, reffes, cell_type)
+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
+  _cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
+
+  return nodes_coords
+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) = grid.cell_type
+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,reffe::ReferenceFE; kwargs...)
+  mapped_grid = GridWithFEMap(model,reffe;kwargs...)
+  MappedDiscreteModel(model,mapped_grid)
+end
+
+function DiscreteModelWithFEMap(model,cell_type,reffes::AbstractArray{<:ReferenceFE}; kwargs...)
+  mapped_grid = GridWithFEMap(model,cell_type,reffes;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,104 @@
+# module ModelWithFEMaps
+
+# using Gridap
+# using Gridap.Arrays
+# using Gridap.Fields
+# using Gridap.ReferenceFEs
+# using Gridap.Geometry
+# using Gridap.CellData
+# using Gridap.FESpaces
+
+# import Gridap.ReferenceFEs.get_node_coordinates
+# import Gridap.Geometry.get_cell_node_ids
+# import Gridap.Geometry.get_reffes
+# import Gridap.Geometry.get_cell_type
+# import Gridap.Geometry.Grid
+
+# import Gridap.Geometry.get_cell_map
+# import Gridap.Geometry.get_grid
+# import Gridap.Geometry.get_grid_topology
+# import Gridap.Geometry.get_face_labeling
+
+# """
+#     MappedGrid
+
+# Reprepresent 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 _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
+
+    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
+
+Reprepresent 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
+
+# end # module